RE: big number arithmetic

To elaborate on Terje's comment about combining shorter building
blocks...

Assuming the word size of your cpu is 64 bits, form upper and lower 32
bit words from each 64 bit value:
A_upper =3D a64  >> 32
A_lower =3D a64 & 0x00000000FFFFFFFF
B_upper =3D b64 >> 32
B_lower =3D b64 & 0x00000000FFFFFFFF

So we have that a64 =3D A_lower + A_upper * 2^32, and b64 =3D B_lower +
B_upper * 2^32

The product is=20
a64*b64 =3D (A_lower + A_upper * 2^32) * (B_lower + B_upper * 2^32)

Expanding:
a64*b64 =3D (A_lower * B_lower) + (A_lower * B_upper)*2^32 + (A_upper *
B_lower)*2^32 + (A_upper * B_upper)*2^64

Each term in parentheses is a product of 32 bit numbers, and so will fit
in the machine's 64 bit word.

The multiplication function will return a product (lower 64 bits) and a
carry (upper 64 bits).  Here's an example I wrote in C for a machine
with 32 bit words:

void mul(unsigned long u, unsigned long v, unsigned long *product,
unsigned long *carry)
{
	unsigned long x0,x1,y0,y1;
	unsigned long ax0,ax1,ay1,az1;

	//split into half words
	x0 =3D u & 0x0000FFFF;
	x1 =3D u >> 16;
	y0 =3D v & 0x0000FFFF;
	y1 =3D v >> 16;

	//half words for product calculation
	ax0 =3D (x0*y0) & 0x0000FFFF;
	ax1 =3D (x0*y0) >> 16;
	ay1 =3D (x1*y0) & 0x0000FFFF;
	az1 =3D (x0*y1) & 0x0000FFFF;

	//calculate product
	*carry =3D ax1+ay1+az1;
	*product =3D ax0 + ((*carry & 0x0000FFFF) << 16);
	*carry =3D (*carry) >> 16;

	//calculate carry
	*carry +=3D ((x1*y0) >> 16) + ((x0*y1) >> 16) + y1*x1;
=09
	return;
}

It should be easy to modify it for 64 bit operands.

For generic N bit multiplies, there are a number of methods depending on
how big your numbers are.  In each case, the input numbers are
represented as arrays of machine words.  Smallish numbers can use
straightforward O(N^2) gradeschool methods and use the above code for
each individual multiply operation.  Bigger numbers can be tackled by
Karatsuba, Toom-Cook, or FFT methods.  Knuth TAOCP vol 2, or the GMP
website might provide more details.

Regards,
- ben.



> -----Original Message-----
> From: [email protected]=20
> [mailto:[email protected]]
>  On Behalf Of Terje Mathisen
> Sent: Thursday, March 22, 2007 7:19 AM
> To: [email protected]
> Cc: [email protected]
> Subject: Re: big number arithmetic
>=20
> [email protected] wrote:
> > Dear list members,
> >=20
> > i am trying to imlement multiplication function for two 64=20
> bit number.
> > I wondered which would it be the best possible algorithm for such a=20
> > feat? May someone suggest one?
> > After have implemented, what about make it generic, i.e.,=20
> for n bits?
>=20
> What kind of cpu are you working on?
>=20
> Do you need a 128-bit full product, or just the low half?
>=20
> The latter is easy, since any conforming C(++) compiler has to support
> t_uint64 and t_int64 these days.
>=20
> If your building blocks are shorter multiplies, maybe=20
> 32x32->64, then you can combine those to generate the full=20
> product you need.
>=20
> Assuming you have a cpu with no multiplication support at=20
> all, then you need to do some more work:
>=20
> Either a full 64-iteration loop using shift & add, or you can=20
> construct a small multiplier using lookup tables:
>=20
>   (a+b)^2 =3D a^2 + 2ab + b^2
>   (a-b)^2 =3D a^2 - 2ab + b^2
>=20
> so (a+b)^2 - (a-b)^2 =3D 4ab.
>=20
> With a lookup table of the first 2N squares, you can use this=20
> to multiply 2 N-bit numbers!
>=20
> unsigned mul8x8(unsigned a, unsigned b)
> {
>    static unsigned square[512] =3D {0,1,4,9,16,25,36,49....};
>    if (a < b) { unsigned t =3D a; a =3D b; b =3D t};
>    unsigned prod =3D (square[a+b] - square[a-b]) >> 2;
>    return prod;
> }
>=20
> If we note that the sum and difference of two even numbers=20
> will both be even, the sum/diff of an odd and an even number=20
> will both be odd, and the square of an odd number is odd,=20
> then we realize that the table can store the truncated half=20
> of each square, it will still come out correct, and we'll=20
> only need a single shift at the end.
>=20
> Terje
>=20
> --
> - <[email protected]>
> "almost all programming can be viewed as an exercise in caching"
>=20