Optimization & Efficiency

Short Segments
High Precision with Imprecise Methods

Introduction

When floating-point numbers fail to provide the required precision for our applications, we must resort to more creative solutions. Surprisingly, one of those solutions involves use of integers and frequent rounding. The technique presented here involves the use of integers that are conceptually unlimited in size. Rounding is the bane of precision when it is unintentional, but when used intelligently with these very large integers, it enables us to maintain both efficiency and precision. In this article I show how.

Since I needed integers that could grow to virtually any size, I developed the BigNum class, described in last month's CUJ (see "A Class for Representing Large Integers," CUJ, November 1996). BigNums add, subtract, multiply, and divide with other integers and other BigNums. They act pretty much like ints and longs but without their respective size constraints. (Of course, even BigNums can only get so big. The point is, in most cases you can treat them as if they are unlimited in size.) I provide the source code for BigNum on this month's code disk and online sources (see p. 3 for details). If you have your own equivalent implementation for very large integers in C or C++, you should be able to adapt it to work with the technique presented here.

The Technique

Consider the following approximation:

arctan(x) = x - (x3)/3 + (x5)/5 - (x7)/7 +...+
    (-1n)(x(2n+1))/(2n+1)

Ignoring several minor details for now, notice that computing an arctan is equivalent to adding and subtracting several fractions. The high-precision technique uses a shortcut for adding and subtracting fractions, while maintaining the desired number of digits of accuracy. This technique can be applied to a large variety of problems where fractions are involved and accuracy is important.

The best way to explain the technique is with a simple example. Consider calculating 5/37 + 11/31 + 19/39 to 100 digits of accuracy. After finding a common denominator the problem becomes 6045/44733 + 15873/44733 + 21793/44733 = 43711/44733. Floating-point or double-precision arithmetic won't give the desired 100 digits of accuracy so instead we compute ((10¹⁰⁰) * 43711)/44733 using the BigNum routines. The new quotient will contain the correct string of digits; we need only remember that the decimal point should be shifted 100 digits to the left.

This solves one problem but introduces another. Adding only three fractions with two-digit denominators created a common denominator with five digits; things will quickly get out of hand if we have to find a common denominator for each term of the arctan series. In particular, computing just 100 terms of arctan(1/239) (important in calculating pi — see below) would cause both the numerator and denominator to grow to over 500 digits.

We can avoid the common denominator problem by rounding off each term to an integer. To see how, consider again the sum 5/37 + 11/31 + 19/39. If we multiply this sum by the very large number 10¹¹⁰, we get 5*(10¹¹⁰)/37 + 11*(10¹¹⁰)/31 + 19*(10¹¹⁰)/39. This sum is very nearly equal to floor(5*(10¹¹⁰)/37) + floor(11*(10¹¹⁰)/31) + floor(19*(10¹¹⁰)/39).

By truncating each computation to an integer (via floor) we've introduced an error of less than 1 for each term, causing the final computation to be off by about 3. However, by multiplying each term by 10¹¹⁰ we've given ourselves some breathing room. 110 digits were computed, the last one or two of which are wrong. Thus, to get our 100 digits of accuracy we merely move our decimal point 110 digits left and ignore the last 10 digits.

There are two other benefits to this algorithm: we don't need to find the least common denominator, and the denominators won't grow. While this is by no means earth shattering, it is somewhat surprising that a shortcut that rounds off the computation at each step can be used for computations of any desired precision.

Conclusion

This technique can be used to compute pi, using Machin's formula: pi = 16*arctan(1/5) + 4*arctan(1/239). With Machin's formula I can compute pi to 100 digits in less than one second on a 66 MHz 486, and to 10,000 digits in about ten minutes. (I provide the subroutines to do this on this month's code disk.) This isn't the fastest method available to compute pi, but it can be easily modified to solve a variety of problems.

Once again the basic idea is to multiply your problem by a huge number and compute the answer quickly by rounding off intermediate results, while ensuring that the roundoff error is sufficiently dominated by the multiplier. o

Bibliography

Knuth, Donald. The Art of Computer Programming: Semi-Numerical Methods, Second Edition, 1981.

Beckmann, Petr. A History of Pi, 1971.

Thanks to Anthony Breitzman <breitzCHI@aol.com> for this tip.