views:

525

answers:

7

I need some division algorithm which can handle big integers (128-bit). I've already asked how to do it via bit shifting operators. However, my current implementation seems to ask for a better approach

Basically, I store numbers as two long long unsigned int's in the format

A * 2 ^ 64 + B with B < 2 ^ 64.

This number is divisible by 24 and I want to divide it by 24.

My current approach is to transform it like

A * 2 ^ 64 + B     A             B
--------------  = ---- * 2^64 + ----
      24           24            24

           A               A mod 24                    B         B mod 24
= floor( ---- ) * 2^64 +  ---------- * 2^64 + floor( ---- ) +   ----------
           24               24.0                      24           24.0

However, this is buggy.

(Note that floor is A / 24 and that mod is A % 24. The normal divisions are stored in long double, the integers are stored in long long unsigned int.

Since 24 is equal to 11000 in binary, the second summand shouldn't change something in the range of the fourth summand since it is shifted 64 bits to the left.

So, if A * 2 ^ 64 + B is divisible by 24, and B is not, it shows easily that it bugs since it returns some non-integral number.

What is the error in my implementation?

+1  A: 

You shouldn't be using long double for your "normal divisions" but integers there as well. long double doesn't have enough significant figures to get the answer right (and anyway the whole point is to do this with integer operations, correct?).

Andrew Jaffe
the whole point of it is to just divide a 128 bit int by 24 which results in epic fail at the moment.Long double has 64 bit of mantissa, so this shouldn't make a problem. or does it?
Etan
Etan should have linked to eir original question. It seems that the objective is not to do it with integers but to do it at all. Besides, `long double` can be as small as a 64-bit double, but can also be larger (say a 10-byte extended double, but anything at all, really... IEEE 754 is mostly parametric with respect to bit sizes). Therefore it *could* possibly have the necessary precision (I'm not saying it's a good idea to use floating-point computations for something as easy as 128-bit integer arithmetic).
Pascal Cuoq
How to divide it without the long double?
Etan
`long long int d = a/24`
Andrew Jaffe
@Andrew You mean `assert (sizeof(long long int) >= 16); d = a/24;` :)
Pascal Cuoq
+1  A: 

Since 24 is equal to 11000 in binary, the second summand shouldn't change something in the range of the fourth summand since it is shifted 64 bits to the left.

Your formula is written in real numbers. (A mod 24) / 24 can have an arbitrary number of decimals (1/24 is for instance 0.041666666...) and can therefore interfere with the fourth term in your decomposition, even once multiplied by 2^64.

The property that Y*2^64 does not interfere with the lower weight binary digits in an addition only works when Y is an integer.

Pascal Cuoq
It does interfere in decimals since you cannot write them exactly down there. In binary, it has a bounded implementation since 1/24 can be written in an ending amount of digits.
Etan
@Etan Really? How many digits does it take to represent 1/24 exactly in binary then? (If that's too difficult a question, start with the number of binary digits it takes to represent 1/3 exactly)
Pascal Cuoq
1/24 = binary 0.00001010101010101... the sequence goes on forever.
Dave Hinton
Damn it ^^ You are right. So, the idea seems to be correct and the sum of both modulo's has to be simply added to the result?
Etan
I was only contributing to the debugging your thought process. I think it's a bad idea to write a formula that only holds for reals and to then approximate these reals with floats (even long doubles). It will have horrible bugs that only happen when dividing large numbers (and only some of them).
Pascal Cuoq
+2  A: 

Don't.

Go grab a library to do this stuff - you'll be incredibly thankful you chose to when debugging weird errors.

Snippets.org had a C/C++ BigInt library on it's site a while ago, Google also turned up the following: http://mattmccutchen.net/bigint/

Adam Frisby
I have to do it manually since it is for an ACM ICPC problem.
Etan
+9  A: 

The easiest way I can think of to do this is to treat the 128-bit numbers as four 32-bit numbers:

A_B_C_D = A*2^96 + B*2^64 + C*2^32 + D

And then do long division by 24:

E = A/24 (with remainder Q)
F = Q_B/24 (with remainder R)
G = R_C/24 (with remainder S)
H = S_D/24 (with remainder T)

Where X_Y means X*2^32 + Y. Then the answer is E_F_G_H with remainder T. At any point you only need division of 64-bit numbers, so this should be doable with integer operations only.

interjay
It doesn't prevent your algorithm from working at all, but F, G and H can each be larger than 2^32. I had difficulties to reconcile that with the fact that the `E_F_G_H` notation looks like concatenation, but once this is understood, it's a very nice algorithm.
Pascal Cuoq
Actually, F,G, and H will be smaller than 2^32, because Q, R, and S are all less than 24. So the E_F_G_H notation *is* concatenation.
interjay
Right! I forgot my pencil-and-paper division... I remembered there was an unpleasant guessing part but that's when the divisor has too many digits. As long as the divisor itself is short enough to fall within the range for which primitive division works (as is the case here), there never is any need to guess when applying the pencil-and-paper division algorithm. Sorry for the confusion.
Pascal Cuoq
+2  A: 

Could this possibly be solved with inverse multiplication? The first thing to note is that 24 == 8 * 3 so the result of

a / 24 == (a >> 3) / 3

Let x = (a >> 3) then the result of the division is 8 * (x / 3). Now it remains to find the value of x / 3.

Modular arithmetic states that there exists a number n such that n * 3 == 1 (mod 2^128). This gives:

x / 3 = (x * n) / (n * 3) = x * n

It remains to find the constant n. There's an explanation on how to do this on wikipedia. You'll also have to implement functionality to multiply to 128 bit numbers.

Hope this helps.

/A.B.

Andreas Brinck
A: 

I don't know but I like really big numbers...So I checked out Fourier Division...It is sweet....

:) -Donaldli Ingalls

Donaldli Ingalls
A: 

oh btw... You should check out that Fourier Division on Wikipedia

Wikipedia rocks!

-Donaldli

Donaldli Ingalls