Re: big number arithmetic

Thank you all a lot!

I would like to go for karatsuba at a first try. I have implemented a
code to multiple two 64 bit number on intel P4 machine. Here it is.

(a0 + a1 * 2 ^ 32) * ( b0 + b1 * 2 ^ 32) = a0b0 + (a0b0 + a1b1 + (a0 -
a1)*(b1 - b0)) * 2 ^ 32 + a1b1 * 2 ^ 64.

the problem is that i am having a hard time to implement sign
manipulation. Does anybody here have a function that implements such
algorithm (karatsuba).

Here is mine, but it is buggy.

typedef unsigned long long xadk64_t;
typedef unsigned xadk32_t;
typedef int adk32_t;

void
karatsuba(xadk64_t * const r, xadk64_t a, xadk64_t b)
{
        xadk32_t        a0, a1, b0, b1;
        xadk64_t        x, y;

        a0 = a, a1 = a >> 32, b0 = b, b1 = b >> 32;
        a = (xadk64_t)a0 * a1, b = (xadk64_t)b0 * b1;
        x = (xadk64_t)(adk32_t)(a0 - a1) * (xadk64_t)(adk32_t)(b1 - b0);
        if ((y = a + b) < b) b += (xadk64_t)1 << 32;
        x += y;
        if ((a += x << 32) < x << 32) b++;
        b += x >> 32;
        r[0] = a, r[1] = b;
}

I am really having a bad time trying to implement it correctly.
Any code anyone has that implements the above?

Thanks a lot for your time and cooperation.

Best regards.

On 3/22/07, Buhrow, Benjamin <[email protected]> wrote:
> 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 = a64  >> 32
> A_lower = a64 & 0x00000000FFFFFFFF
> B_upper = b64 >> 32
> B_lower = b64 & 0x00000000FFFFFFFF
>
> So we have that a64 = A_lower + A_upper * 2^32, and b64 = B_lower +
> B_upper * 2^32
>
> The product is
> a64*b64 = (A_lower + A_upper * 2^32) * (B_lower + B_upper * 2^32)
>
> Expanding:
> a64*b64 = (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 = u & 0x0000FFFF;
>         x1 = u >> 16;
>         y0 = v & 0x0000FFFF;
>         y1 = v >> 16;
>
>         //half words for product calculation
>         ax0 = (x0*y0) & 0x0000FFFF;
>         ax1 = (x0*y0) >> 16;
>         ay1 = (x1*y0) & 0x0000FFFF;
>         az1 = (x0*y1) & 0x0000FFFF;
>
>         //calculate product
>         *carry = ax1+ay1+az1;
>         *product = ax0 + ((*carry & 0x0000FFFF) << 16);
>         *carry = (*carry) >> 16;
>
>         //calculate carry
>         *carry += ((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 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]
> > [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 algorithm for such 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 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
> > all, then you need to do some more work:
> >
> > Either a full 64-iteration loop using shift & add, or you can
> > construct a small multiplier using lookup tables:
> >
> >   (a+b)^2 = a^2 + 2ab + b^2
> >   (a-b)^2 = a^2 - 2ab + b^2
> >
> > so (a+b)^2 - (a-b)^2 = 4ab.
> >
> > With a lookup table of the first 2N squares, you can use this
> > to multiply 2 N-bit numbers!
> >
> > unsigned mul8x8(unsigned a, unsigned b)
> > {
> >    static unsigned square[512] = {0,1,4,9,16,25,36,49....};
> >    if (a < b) { unsigned t = a; a = b; b = t};
> >    unsigned prod = (square[a+b] - square[a-b]) >> 2;
> >    return prod;
> > }
> >
> > If we note that the sum and difference of two even numbers
> > will both be even, the sum/diff of an odd and an even number
> > will both be odd, and the square of an odd number is odd,
> > then we realize that the table can store the truncated half
> > of each square, it will still 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"
> >
>