Re: big number arithmetic

Ok, now i have implemented karatsuba, but running into a error when f
variable is not 0 (which means the product is negative). Some ligth on
this regards ?

void
k(xadk64_t * const r, xadk64_t a, xadk64_t b)
{
        xadk32_t        a0, a1, b0, b1;
        xadk64_t        x, y;
        int             f = 0;

        a0 = a, a1 = a >> 32, b0 = b, b1 = b >> 32;
        a = (xadk64_t)a0 * a1, b = (xadk64_t)b0 * b1;
        if (a0 < a1) f += 1;
        if (b1 < b0) f += 2;
        switch (f) {
        case 1 :
                x = (adk32_t)(a0 - a1)*(xadk64_t)(b1 - b0);
        break;
        case 2 :
                x = (xadk64_t)(a0 - a1)*(adk32_t)(b1 - b0);
        break;
        default :
                f = 0;
                x = (xadk64_t)(adk32_t)(a0 - a1)*(xadk64_t)(adk32_t)(b1 - b0);
        break;
        }
        fprintf(stdout, "x = %llu\n", x);
        if (f) {
                if ((x += b) < b) b += (xadk64_t)1 << 32;
                if ((x += a) < a) b += (xadk64_t)1 << 32;
                b -= x >> 32;
                if (y = a, ((a -= x << 32) > y)) b--;
                fprintf(stdout, "A: %llu, B: %llu\n", a, b);
        } else {
                if ((x += b) < b) b += (xadk64_t)1 << 32;
                if ((x += a) < a) b += (xadk64_t)1 << 32;
                if ((a += x << 32) < x << 32) b++;
                b += x >> 32;
                fprintf(stdout, "a: %llu, b: %llu\n", a, b);
        }
}

On 3/22/07, John Nietzsche <[email protected]> wrote:
> 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"
> > > > >
> > > >
> > >
> >
>