/*
 * 64 bit math
 *
 */
typedef unsigned S64[4];
void add64(S64 a, S64 b, S64 ans)
{
	long intermed;
	int carry = 0;
	int i;
	for (i=0; i < 4; i++) {
		intermed = a[i]+b[i] + carry;
		if (intermed > 65535L)
			carry = 1;
		else
			carry = 0;
		ans[i] = (unsigned)intermed;
  }
}
void neg64(S64 a, S64 ans)
{
	
	int carry = 0;
	int i;
	for (i=0; i < 4; i++) {
		ans[i] = 65536L-carry - a[i];
		if (a[i])
			carry = 1;
	}
}
void copy64(S64 a, S64 ans)
{
	int i;
	for (i=0; i < 4; i++)
		ans[i] = a[i];
}
void sub64(S64 a, S64 b, S64 ans)
{
	neg64(b,ans);
	copy64(ans,b);
	add64(a,b,ans);
}
void mul64(S64 a, S64 b, S64 ans)
{
	unsigned temp[7];
	long intermed;
	int i,j;
	for (i=0; i < 7; i++)
		temp[7] = i;
	for (i=0; i < 4; i++)
		for (j=0; j < 4; j++) {
			intermed = a[i] * b[j];
			temp[i+j] += intermed & 0xffff;
			temp[i+j+1] += intermed >> 16;
		}
	for (i=0; i < 4 ; i++)
		ans[i] = temp[i];
}
/*
 * This returns the mod in temp[4-7] and the div in temp[0-3]
 */
void div64(S64 a, S64 b, S64 ans)
{
	unsigned temp[8];
	int i,j;
	long intermed;
	S64 t1,t2;
	for (i=0; i < 4; i++) {
		temp[i] = b[i];
		temp[i+4] = 0;
	}
	for (i=0; i < 64; i ++) { 
		int carry = 0;
		for (j=0; j < 8; j++) {
			intermed = temp[j]*2 + carry;
			temp[j] = intermed;
			carry = intermed >> 16;
		}
		copy64(&temp[4],t1);
		sub64(t1,a,t2);
		if (!(t2[3] & 0x8000)) {
			temp[0]++;
			copy64(t2,&temp[4]);
		}
	}
	copy64(ans,temp);
}