/*Multiply 2 positive variables

  address@hidden
  April 27th., 2014.

Let us work through multiplying 13 by 5 in longhand. I shall
represent 13 as 1*2**0+0*2**1+1*2**2+1*2**3 or 1011; I shall
represent 5 as 1*2**0+0*2**1+1*2**2 or 101.

1011 * 101 can be calculated:

1011   +
 0000  +
  1011 =
102111 or 1*2**0+0*2**1+2*2**2+1*2**3+1*2**4+1*2**5 which
is represented as 65 in decimal.

I need a data structure to hold this calculation. In some
languages it is possible to write something like:
   gn{x in 0 .. 4, y in x .. 4+x}
mathprog doesn't complain if you try this, but neither does it
produce a useful result.

I could define a data structure:
101100
000000
001011
and read the result by summing the columns but for this model
I have chosen to define:
1011
0000
1011
with the cost that the result must now be extracted by summing
the diagonals from right to left. The advantage for explaining
the solution is that if I now add a third dimension and I fill in the front layer with 13 horizontally:
1011
1011
1011
and the rear layer with 5 vertically:
1111
0000
1111
I think that the reader will easily see that the result that
I require is the logical and of the front and rear layer. The two
rules gn3 and gn4 achieve this.

Vx, Vy, and P define the maximum size of the two variables and
the precision of the variables in decimal digits. The model is
fast so one does not need to be mean in setting them.

i1 through i8 convert the variables from decimal to the format
I have described above.

r1 extracts the result, for the model.

*/

param Vx  :=  17;
param Vy  :=  14;
param P   :=   2;

var nx;
var nIx, integer;
var nBx{x in 0..Vx}, binary;
var ny;
var nIy, integer;
var nBy{y in 0..Vy}, binary;
var gn{x in 0..Vx, y in 0..Vy, z in 1 .. 2}, binary;
var myPROD;

i1: nx = 923.25; /*set value of 1st variable*/
i2: nIx <= nx*10**P+0.5;
i3: nIx >= nx*10**P-0.5;
i4: ny = 123.75; /*set value of 2nd variable*/
i5: nIy <= ny*10**P+0.5;
i6: nIy >= ny*10**P-0.5;
i7: sum{x in 0..Vx} nBx[x]*2**(x) = nIx;
i8: sum{y in 0..Vy} nBy[y]*2**(y) = nIy;

gn3{x in 0..Vx, y in 0..Vy}: sum {z in 1..2} gn[x,y,z] = nBx[x]+nBy[y]; 
gn4{x in 0..Vx, y in 0..Vy}: gn[x,y,1] <= gn[x,y,2];

r1: myPROD*10**(2*P) = sum{x in 0 .. Vx+Vy} sum{y in 0 .. x : x-y <= Vy and y <= Vx} gn[y,x-y,1]*2**(x);

solve;

printf "\n\n%f * %f = %f\n\n", nx, ny, myPROD;
  
end;