I was working on a project that required simple arithmetic for very large integers, a set of algorithms called “Arbitrary Precision Math”.
Thinking back to elementary school, simple algorithms exist for addition, subtraction, and multiplication of two numbers with any number of digits.
To my surprise, every algorithm for division either relies on logarithms, which are difficult to implement in arbitrary precision, or the first instruction was “guess the first number, then guess the second number” etc…
Update: 10/2015: I’ve put together a YouTube video for this post. Check it out, here:
Read on, for a simple, reliable, repeatable algorithm for dividing integers of any length.
1.1. What is Arbitrary Precision Math?
Computers store numbers based on a variable’s data type, and the data type determines the largest number that you can store:
|Bits||Bytes||Data Type||Largest Signed Integer||Largest Unsigned Integer|
Notice the trailing zeros for the Int64 data type – the program I used to perform the calculation doesn’t have enough precision to calculate the actual digits, resulting in rounding errors!
Being able to represent any 18-digit number is fairly impressive, if you have a compiler and operating system that support the Int64 data type.
What if you need 20 digits? Or 100 digits?
Although you could implement an arbitrary-precision library using strings, this is increasingly inefficient as the number of digits increases, because each digit requires a whole byte to encode.
A better solution is to create a multi-segment string, whose segments are each one of the base integer data types. For example, a 16-bit Int can represent any base-10 4-digit number: 0000 through 9999. Let’s say you have a data type and custom library that implements a string of Ints, up to 30 positions. This custom data type can represent an unsigned number up to 30 * 4, or 120 digits in length!
Just as you can represent extremely long integers using arbitrary-precision math, you can represent long strings of decimal digits as well, because you can choose where to put the decimal point.
1.2. What is Division?
When we say, “divide A by B”, “A” is the Numerator (Dividend), “B” is the Denominator (Divisor).
The result is called the Quotient. In integer division, the fractional portion is expressed as a Remainder, which represents the ratio of “what’s left over” to the denominator / divisor.
- Expressed as a formula:
N/D = Q + R/D
- N – Numerator, or Dividend, is the number being divided
- D – Denominator, or Divisor, is the size of the bucket used to slice up N.
- Q – Quotient (result), or the number of slices created if we carve up N in to D-sized buckets
- R – Remainder, or the quantity of N that’s left over, and won’t quite fill a whole D-sized bucket. So the remainder is actually a fraction, a portion of 1, expressed as R/D where R/D < 1
- Multiply the quotient / remainder by the denominator, to reconstruct N. Remember that the remainder is a fraction, R/D:
N/D = Q + R/D
Multiply by D:
N = Q * D + R
1.3. Finding the First Digit
Let M be the magnitude in digits (minus 1) of the divisor, D.
So if D = 6789, M would be 3 (4 digits, minus 1)
A is the “quick divisor”, or the first digit of the divisor, with all other digits equal to zero. Expressed as a formula:
A = D – ( D MOD 10^M )
In our example above:
A = 6789 – ( 6789 MOD 10^3 )
A = 6789 – ( 6789 MOD 1000 )
A = 6789 – 789
A = 6000
We’re trying to set all digits to zero, except the first digit. What’s the first thing you do in long division? Drop all the zeros in the denominator. So if we were trying to divide 567,892,345 by 6,789 Reducing 6,789 to 6,000 allows us to reduce the precision of the numerator by the same number of digits. We now have a single-digit math problem!
567,892,345 / 6,789 (Difficult)
567,892,345 / 6,000 = 567,892 / 6 (Easy)
The “quick divisor”, A, will allow us to treat the entire division problem as single-digit division, easily computable by solving one digit at a time.
2. Process Overview
The process uses error-correcting iterative steps to converge on the quotient and remainder.
Each iterative quotient candidate is multiplied by the divisor, and the difference between the result and the numerator is halved.
This results in a new quotient candidate. To more quickly converge, the previous two quotient candidates can be averaged.
Once there is no further oscillation, and the absolute value of the remainder is less than the denominator (ratio of R/D <1) then the proper quotient has been determined.
3. Algorithm and Detailed Walkthrough
Recap of variables:
- N – Numerator
- D – Denominator
- Q – Quotient
- R – Remainder
- M – Magnitude
- A – Quick Divisor
Here is the algorithm (Assume all calculations are integer calculations):
- M = length in digits of D, minus 1. So if D is 1234, there are 4 digits, M=3.
- A = D – (D MOD 10^M). This results in the first digit of D, followed by all zeros
- Q = N / A
- R = D+1
- While ABS(R)>=D
- R = N – (Q * D)
- Qn = Q + R / A
- Q = (Q + Qn) / 2
- R = N – (Q * D) Calculate a final value for R
- If R<0 then
Once the loop exits, Q is the quotient, R is the remainder.
3.2. Walk Through – How does this work?
When you have to “guess” the quotient digit based on a multi-digit divisor, how do you do that? Even the guess is an iterative process.
Conversely, it’s easy to do a single-digit division. Grab a pair of digits, divide, take the remainder times 10, grab the next digit, etc…
So we start by creating a “quick divisor”, A, that can be calculated using single-digit division.
The result is way off, but this is an iterative process.
We start with an initial value for Q by using single-digit division to calculate N / A.
3.2.1 Start of Loop
The remainder, R, is the difference between our numerator N, and Q * D (our initial quotient candidate, times the denominator)
By dividing R by A (the quick divisor), we can do another simple, single-digit division step, to get a value that can be added to Q, in order to make Q more accurate.
Taking this approach results in oscillation. To dampen the oscillation, we average Q and Qn:
- The original quotient candidate is Q.
- New Q, or Qn = Q (the original quotient candidate) + R / A (the error or offset)
- Q = (Q + Qn) / 2 (average the two, to dampen oscillation, and reduce iterations by half)
Once R settles in to a an absolute value <=D, you have the right answer, Q
A final computation step must be added, to calculate the final value for R based on the final value for Q, otherwise, you’d have the R corresponding to the old Q value.
In some cases, R can stay negative, but converge to a an absolute value <=D. In this case, a final adjustment must be made:
- Subtract 1 from Q – forces Q * D to a value <= N
- Add D to the remainder, R, representing the “1” we subtracted from Q, and forcing R positive
Q is now the proper quotient, and R is the remainder.
4.1. Example 1
Divide 999,999 by 7,777
N (Numerator) = 999,999
D (Denominator) = 7,777
M = Len(D)-1, so M=3
A = D – (D MOD 10^M)
A = 7,777 – (7,777 MOD 1000)
A = 7,777 – 777
A = 7,000
Q = N / A
A = 999,999 / 7,000
Q = 142 (remember, this is all integer math)
R = D+1 = 7,778 (initial value, false condition for while loop)
|Loop||Old Q||Old Q x D||R||R / A||Qn||Q||Guess|
Compute the final value for R:
R = N – ( Q * D )
R = 999,999 – (129 * 7,777)
R = -3,234
Since R is negative, we have a final adjustment to Q and R
Q = Q -1
Q = 129-1
Q = 128
R = R + D
R = -3,234 + 7,777
R = 4,543
N = Q * D + R
999,999 = 128 * 7,777 + 4,543
999,999 = 995,456 + 4,543
999,999 = 999,999
4.2. Example 2
Divide 999,999,999 by 999,999
N (Numerator) = 999,999,999
D (Denominator) = 999,999
M = Len(D)-1, so M=5
A = D – (D MOD 10^M) = 900,000
Q = N / A = 1,111
R = D+1 = 1,000,000 (initial value, false condition for while loop)
|Loop||Old Q||Q x D||R||R / A||Qn||Q||Guess|
Compute the final value for R:
R = N – ( Q * D )
R = 999,999,999 – (1,000 * 999,999)
R = 999
Since R positive, we don’t need to adjust Q nor R
N = Q * D + R
999,999,999 = 1,000 * 999,999 + 999
999,999,999 = 999,999,000 + 999
999,999,999 = 999,999,999
5. Corollary in Base X
The nice thing about this algorithm is that it works in any base!
The examples above are in base 10.
Recalling our multiple-segment arbitrary-precision data type, we used a string of Int16’s, of which each Int16 is used to store a 4-digit segment of some very large number.
If we consider each segment a base-10,000 digit, then we can perform the above algorithm one whole segment at a time instead of one digit at a time!
In the above example, N = 999,999,999 and D = 999,999, these numbers would be represented by a string of Int16’s like this:
- N = 9:9999:9999 (3 segments)
- D = 99:9999 (2 segments)
- M = 1
- A = 99:9999 – ( 99:9999 MOD 10,000^M ) = 99:0000
The first step would be to find the initial value of Q, by dividing N by A.
This can be done one “digit” at a time:
9:9999:9999 / 99:0000
Fix up the precision based on A:
9:9999 / 99
Compare one digit at a time:
9 / 99 = 0, R 9
The next segment multiplies the remainder times the base, adding the next digit:
9 * 10,000 + 9999 = 99,999
99,999 / 9 is easily computable using 32-bit math. This results in an initial Q = 1,010
Notice that we had to double our precision to accommodate carry / borrow. This means that the segment data type has to be able to accommodate n * b + n, where n is the largest number per segment, and b is the base.
For example, since Int16 can accommodate any 4-digit number, it’s more expedient, if limited to 16-bit math, to limit the segment size to 2 digits. This allows 99 * 100 + 99, or 9999 to be stored in a single segment, to allow for carry / borrow.
On the other hand, if we continue to utilize 16-bit segments, we need a 32-bit data type to perform the segment-by-segment math.
The algorithm outlined above can be used to perform arbitrary-precision division using character strings, or multi-segment strings of integers, and the concept works in any base.