Division of big numbers

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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?).

Answer 4

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.

Answer 5

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/