RE: big number arithmetic

I agree with Felix that you should stick with the straightforward method
for small multiplications such as 64x64.  Assembler level instructions
will probably be the fastest.  These C routines are strictly for if you
can't (for portability) or don't want to write assembly.  Libraries like
GMP usually use assembly for these "low level" tasks, and only utilize
Karatsuba or FFT if the number of bits to be multiplied is huge.

As far as your sign handling, you are throwing away information by
casting to int's and then multipling.  What's to stop (a0-a1) or (b1-b0)
from being greater than 2^31 ?

Just do the subtraction with a1 and a0 as uint64 values and then
multiply.  If the result is negative (which requires if statements...
slow slow slow in a time critical routine such as this), then the
magnitude can be found like so:

-u =3D ~u + 1, where ~ is bitwise not.

Hope this helps,
- ben.=20

> -----Original Message-----
> From: John Nietzsche [mailto:[email protected]]=20
> Sent: Thursday, March 22, 2007 11:30 AM
> To: Buhrow, Benjamin
> Cc: [email protected]
> Subject: Re: big number arithmetic
>=20
> Thank you all a lot!
>=20
> I would like to go for karatsuba at a first try. I have=20
> implemented a code to multiple two 64 bit number on intel P4=20
> machine. Here it is.
>=20
> (a0 + a1 * 2 ^ 32) * ( b0 + b1 * 2 ^ 32) =3D a0b0 + (a0b0 + a1b1 + (a0 =
-
> a1)*(b1 - b0)) * 2 ^ 32 + a1b1 * 2 ^ 64.
>=20
> the problem is that i am having a hard time to implement sign=20
> manipulation. Does anybody here have a function that=20
> implements such algorithm (karatsuba).
>=20
> Here is mine, but it is buggy.
>=20
> typedef unsigned long long xadk64_t;
> typedef unsigned xadk32_t;
> typedef int adk32_t;
>=20
> void
> karatsuba(xadk64_t * const r, xadk64_t a, xadk64_t b) {
>         xadk32_t        a0, a1, b0, b1;
>         xadk64_t        x, y;
>=20
>         a0 =3D a, a1 =3D a >> 32, b0 =3D b, b1 =3D b >> 32;
>         a =3D (xadk64_t)a0 * a1, b =3D (xadk64_t)b0 * b1;
>         x =3D (xadk64_t)(adk32_t)(a0 - a1) *=20
> (xadk64_t)(adk32_t)(b1 - b0);
>         if ((y =3D a + b) < b) b +=3D (xadk64_t)1 << 32;
>         x +=3D y;
>         if ((a +=3D x << 32) < x << 32) b++;
>         b +=3D x >> 32;
>         r[0] =3D a, r[1] =3D b;
> }
>=20
> I am really having a bad time trying to implement it correctly.
> Any code anyone has that implements the above?
>=20
> Thanks a lot for your time and cooperation.
>=20
> Best regards.
>=20
> On 3/22/07, Buhrow, Benjamin <[email protected]> wrote:
> > To elaborate on Terje's comment about combining shorter building=20
> > blocks...
> >
> > Assuming the word size of your cpu is 64 bits, form upper=20
> and lower 32=20
> > 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 +=20
> > B_upper * 2^32
> >
> > The product is
> > 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 +=20
> (A_upper *
> > B_lower)*2^32 + (A_upper * B_upper)*2^64
> >
> > Each term in parentheses is a product of 32 bit numbers,=20
> and so will=20
> > fit in the machine's 64 bit word.
> >
> > The multiplication function will return a product (lower 64=20
> bits) and=20
> > a carry (upper 64 bits).  Here's an example I wrote in C=20
> for a machine=20
> > with 32 bit words:
> >
> > void mul(unsigned long u, unsigned long v, unsigned long *product,=20
> > 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;
> >
> >         return;
> > }
> >
> > It should be easy to modify it for 64 bit operands.
> >
> > For generic N bit multiplies, there are a number of methods=20
> depending=20
> > on how big your numbers are.  In each case, the input numbers are=20
> > represented as arrays of machine words.  Smallish numbers can use=20
> > straightforward O(N^2) gradeschool methods and use the=20
> above code for=20
> > each individual multiply operation.  Bigger numbers can be=20
> tackled by=20
> > Karatsuba, Toom-Cook, or FFT methods.  Knuth TAOCP vol 2,=20
> or the GMP=20
> > website might provide more details.
> >
> > Regards,
> > - ben.
> >
> >
> >
> > > -----Original Message-----
> > > From: [email protected]
> > > [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
> > >
> > > [email protected] wrote:
> > > > Dear list members,
> > > >
> > > > i am trying to imlement multiplication function for two 64
> > > bit number.
> > > > I wondered which would it be the best possible=20
> algorithm for such=20
> > > > a feat? May someone suggest one?
> > > > After have implemented, what about make it generic, i.e.,
> > > for n bits?
> > >
> > > What kind of cpu are you working on?
> > >
> > > Do you need a 128-bit full product, or just the low half?
> > >
> > > The latter is easy, since any conforming C(++) compiler has to=20
> > > support
> > > t_uint64 and t_int64 these days.
> > >
> > > If your building blocks are shorter multiplies, maybe
> > > 32x32->64, then you can combine those to generate the full
> > > product you need.
> > >
> > > Assuming you have a cpu with no multiplication support at=20
> all, then=20
> > > you need to do some more work:
> > >
> > > Either a full 64-iteration loop using shift & add, or you can=20
> > > construct a small multiplier using lookup tables:
> > >
> > >   (a+b)^2 =3D a^2 + 2ab + b^2
> > >   (a-b)^2 =3D a^2 - 2ab + b^2
> > >
> > > so (a+b)^2 - (a-b)^2 =3D 4ab.
> > >
> > > With a lookup table of the first 2N squares, you can use this to=20
> > > multiply 2 N-bit numbers!
> > >
> > > 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;
> > > }
> > >
> > > If we note that the sum and difference of two even=20
> numbers will both=20
> > > be even, the sum/diff of an odd and an even number will=20
> both be odd,=20
> > > and the square of an odd number is odd, then we realize that the=20
> > > table can store the truncated half of each square, it will still=20
> > > come out correct, and we'll only need a single shift at the end.
> > >
> > > Terje
> > >
> > > --
> > > - <[email protected]>
> > > "almost all programming can be viewed as an exercise in caching"
> > >
> >
>=20