!2 word Fixed Point library extension for Inform
!Works really dumbly, but works
!License? Let's say - GPL.

!For use by #IFDEF's
Constant FIXED_POINT;

!For general purpose use by your routines
Array fp_temp1 --> 0 0;
Array fp_temp2 --> 0 0;
Array fp_temp3 --> 0 0;
Array fp_temp4 --> 0 0;

!Do not mess with these
Array fp_temp_internal1 --> 0 0;
Array fp_temp_internal2 --> 0 0;

!Useful constants
Array fp_unity --> 1 0;
Array fp_zero --> 0 0;

!If no output array is given, these routines put the answer in this and return it. Note that since it will get overwritten, you can't for example use fp_add(fp_mult(fp1,fp2),fp_mult(fp3,fp4)).
Array fp_return --> 0 0;

[ fp fp d;
if (fp-->0 == 0 && fp-->1 < 0) print "-";
print fp-->0, ".";
for (d=1000 : mod(fp-->1) < d && d~=0: d=d/10) print "0";
if (fp-->1) print mod(fp-->1); 
];

[ fp_gt fp1 fp2;
return (fp1-->0 > fp2-->0 || (fp1-->0 == fp2-->0 && fp1-->1 > fp2-->1));
];

[ fp_lt fp1 fp2;
return (fp1-->0 < fp2-->0 || (fp1-->0 == fp2-->0 && fp1-->1 < fp2-->1));
];

[ fp_eq fp1 fp2;
return (fp1-->0 == fp2-->0 && fp1-->1 == fp2-->1);
];

[ fp_gteq fp1 fp2;
return (fp_gt(fp1,fp2) || fp_eq(fp1,fp2));
];

[ fp_lteq fp1 fp2;
return (fp_lt(fp1,fp2) || fp_eq(fp1,fp2));
];

[ fp_set fp i f;
fp-->0 = i;
fp-->1 = f;
];

[ fp_set_eq fp1 fp2;
fp1-->0 = fp2-->0;
fp1-->1 = fp2-->1;
];

!That's a unary minus, mind you
[ fp_minus fpin fpout;
if (fpout == 0) fpout = fp_return;
fpout-->0 = -(fpin-->0);
fpout-->1 = -(fpin-->1);
return fpout;
];

[ fp_add fpin1 fpin2 fpout;
	if (fpout == 0) fpout = fp_return;
	fpout-->0 = fpin1-->0 + fpin2-->0;
	fpout-->1 = fpin1-->1 + fpin2-->1;
	fp_carry(fpout);
	return fpout;
];

[ fp_carry fp;
if (fp-->1 >= 10000) {fp-->1 = fp-->1 - 10000; (fp-->0)++;}
if (fp-->1 <= -10000) {fp-->1 = fp-->1 + 10000; (fp-->0)--;}
if (fp-->0 ~=0 && sgn(fp-->1) ~= sgn(fp-->0)) {
	fp-->0 = fp-->0 + sgn(fp-->1);
	fp-->1 = fp-->1 - 10000*sgn(fp-->1);
}
];

!                           .
Array fp_work1 --> 0 0 0 0 0 0 0 0 0;
Array fp_work2 --> 0 0 0 0 0 0 0 0 0;
Array fp_work3 --> 0 0 0 0 0 0 0 0 0 0 0 0 0;

[ fp_mult fpin1 fpin2 fpout d1 d2 sign;
if (fpout == 0) fpout=fp_return;
fp_digitise(fpin1, fp_work1);
fp_digitise(fpin2, fp_work2);
for (d1=0:d1<=12:d1++)
	fp_work3-->d1 =0;
for (d1=-4 : d1 <= 4 : d1++)
	for (d2=-4 : d2 <= 4 : d2++)
		{if (d1 + d2 >= -8 && d1 + d2 <= 4)
			fp_work3-->(4-(d1+d2)) 
				= fp_work3-->(4-(d1+d2)) 
					+ fp_work1-->(4-d1) * fp_work2-->(4-d2);
			!fp_digit_carry(fp_work3, 4-(d1+d2));
	}
for (d1=12:d1>=0:d1--) fp_digit_carry(fp_work3, d1);
if (fp_work3-->9 >=5) {(fp_work3-->8)++; fp_digit_carry(fp_work3, 8,1);}
sign = sgn(fpin1-->1) * sgn(fpin2-->1);
fp_dedigitise(fp_work3, fpout);
fpout-->0 = fpout-->0 * sign;
fpout-->1 = fpout-->1 * sign;
return fpout;
];

[ fp_digitise fp array i f p;
	i = mod(fp-->0);
	f = mod(fp-->1);
	array-->0 = i / 10000;
	i = i % 10000;
	for (p=3 : p >= 0 : p--) {
		array-->(4-p) = i / power(10,p);
		i = i % power(10,p);
		array-->(8-p) = f / power(10,p);
		f = f % power(10,p);
	}
];

[ fp_dedigitise array fp digit;
	fp-->0 = 0;
	fp-->1 = 0;
	for (digit = 0 : digit <=4 : digit++)
		fp-->0 = fp-->0 + array-->digit * power(10,(4-digit));
	for (digit = 5 : digit <=8 : digit++)
		fp-->1 = fp-->1 + array-->digit * power(10,(8-digit));
];

[ fp_digit_carry array digit recurse carry;
	if (digit == 0) {
		if (array-->0 > 3) print "Error - Floating point multiplication cannae handle the size of this one!^";
		rtrue;
	}
	carry=array-->digit / 10;
	if (carry) {
		array-->(digit-1) = array-->(digit-1) + carry;
		array-->digit = array-->digit % 10;
		if (recurse) fp_digit_carry(array, digit-1, 1);
	}
];

!Can only handle positive integer exponent - returns 1 for negative (or zero)
[fp_power fp exp fpout;
if (fpout == 0) fpout = fp_return;
fp_set(fp_temp_internal1,1,0);
while (exp-- > 0) {
	fp_mult(fp_temp_internal1, fp, fp_temp_internal1);	
}	
fp_set_eq(fpout, fp_temp_internal1);
return fpout;
];

!Returns 1 if >=0, -1 else
[sgn num;
if (num >= 0) return(1);
else return(-1);
];

!Modulus
[mod num;
return sgn(num)*num;
];

!I can't believe Inform doesn't have an exponential operator.
[ power base power ans;
ans=1;
while (power--) ans=ans*base;
return ans;
];

!Does what it says
[ fp_divide_by_integer fp n fpout;
if (fpout == 0) fpout = fp_return;
fp_reciprocal(n, fp_temp_internal1);
fp_mult(fp, fp_temp_internal1, fpout);
return fpout;
];

!Gives the reciprocal of a positive integer
[ fp_reciprocal n fpout;
if (fpout == 0) fpout = fp_return;
if (n==0) "Bad bad error in fp_reciprocal - tried to find the reciprocal of 0, no less!";
if (n==1) {fp_set(fpout, 1, 0); return fpout;}
fp_set(fpout, 0, 10000/n);
return fpout;
];