Re: big number arithmetic

[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 ?

As several others have already said: Don't do it!!!!

An unsigned integer MUL is only 3-5 times slower than ADD/SUB on most 
architectures, and the MUL is sometimes pipelined as well, making the 
cost even closer to a simple operation.

Using Kara to save a single MUL will _not_ be a win: A single missed 
branch is about twice as costly as a MUL.

Anyway, here's your portable C code:

uint64_t mul64x64_to_128(uint64_t a, uint64_t b; uint64_t *high)
{
   uint32_t al, ah, bl, bh;
   uint64_t t0,

   /* Separate out all four 32-bit words: */
   al = (uint32_t) (a & 0xffffffff); ah = (uint32_t) (a >> 32);
   bl = (uint32_t) (b & 0xffffffff); bh = (uint32_t) (b >> 32);

   /* Calculate the four 32x32->64 MULs: */
   t0 = (uint64_t) al * (uint64_t) bl;
   t1 = (uint64_t) al * (uint64_t) bh;
   t2 = (uint64_t) ah * (uint64_t) bl;
   t3 = (uint64_t) ah * (uint64_t) bh;

   /* Add & propagate carries: */
   c1 = (t0 >> 32) + (t1 & 0xffffffff) + (t2 & 0xffffffff);
   t0 = (t0 & 0xffffffff) | (c1 << 32);
   t3 += (t1 >> 32) + (t2 >> 32) + (c1 >> 32);

   *high = t3;
   return t0;
}


Terje
> 
> 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"
>> > > >
>> > >
>> >
>>
> 
> 


-- 
- <[email protected]>
"almost all programming can be viewed as an exercise in caching"