Re: big number arithmetic

First of all, portability is an issue. I cannot afford the price of
coding in assembler. So i sticked with C.

For evaluating the (a0 - a1)(b1 - b0) it could be that some of the
members lead to a positive value although setting the bit sign. b.v:

a0 = 2 ^ 32 - 1, a1 = 0. How to handle in that cases ? Does anybody
have any suggestion ?

Thanks in advance.

On 3/22/07, Buhrow, Benjamin <[email protected]> wrote:
> 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 = ~u + 1, where ~ is bitwise not.
>
> Hope this helps,
> - ben.
>
> > -----Original Message-----
> > From: John Nietzsche [mailto:[email protected]]
> > Sent: Thursday, March 22, 2007 11:30 AM
> > To: Buhrow, Benjamin
> > Cc: [email protected]
> > Subject: 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"
> > > >
> > >
> >
>