Practical Use of Series Expansions
The series expansion for the exponential, sine, and cosine functions are as follows:
exp(x) = 1 + x + x<SUP>2</SUP>/2! + x<SUP>3</SUP>/3! + x<SUP>4</SUP>/4! + ... sin(x) = x - x<SUP>3</SUP>/3! + x<SUP>5</SUP>/5! - x<SUP>7</SUP>/7! + x<SUP>9</SUP>/9! cos(x) = 1 - x<SUP>2</SUP>/2! + x<SUP>4</SUP>/4! - x<SUP>6</SUP>/6! + x<SUP>8</SUP>/8!
Now assume you want to calculate exp(543.7). (If you are curious, the answer is 1.336317976830752149708709910114E+236.) You could put 543.7 into the above formula and evaluate the series until the final term reached the precision you desire. However, this method will take numerous iterations to converge. Since |x| > 1, the numerator will continue to grow larger, and the series won’t converge until the factorial in the denominator can overtake the numerator. Even though this formula is numerically correct, it is not a practical method to calculate the exponential. What is needed is a translation of x so the magnitude of x is less than one. (Hopefully, x will be significantly smaller than one.) The smaller the magnitude of x, the faster the above series will converge to the desired precision. The translation of x is provided by David H. Bailey’s MPFUN software package. This is a multiple precision math library written in Fortran. Dr. Bailey’s algorithm uses a modification of the series expansion:
exp(t) = (1 + r + r<SUP>2</SUP>/2! + r<SUP>3</SUP>/3! + ...)<SUP>q</SUP>*2<SUP>n</SUP>
where q = 256, r = t'/q, t' = t - n*Log(2), and where n is chosen so that -0.5*Log(2) < t' <= 0.5*Log(2). Reducing t mod Log(2) and dividing by 256 insures that -0.001 < r <= 0.001, which accelerates convergence in the above series.
After the series expansion, you must raise the result to the 256th power. This may seem at first a daunting task, but remember that you can simply square the result the required number of times.
x<SUP>4</SUP> = (x<SUP>2</SUP>)<SUP>2</SUP> x<SUP>8</SUP> = ((x<SUP>2</SUP>)<SUP>2</SUP>)<SUP>2</SUP> x<SUP>16</SUP> = (((x<SUP>2</SUP>)<SUP>2</SUP>)<SUP>2</SUP>)<SUP>2</SUP>) x<SUP>256</SUP> = ((((x<SUP>2</SUP>)<SUP>2</SUP> ... 8 times
Also note that this method requires the value of log(2). The value of log(2) in the MAPM library is accurate to 128 decimal places. If you want an exponential calculation that is accurate to greater than 128 decimal places, the library will re-compute log(2) on the fly as necessary to be as precise as the exponential precision specified.
Now assume you want to calculate the sin(10000.0), where 10000.0 is in radians. This is certainly a valid argument to the sin function, though passing 10,000 to the series expansion will be quite inefficient due to the magnitude of the numerator being greater than one, as was previously discussed for the exponential function. You need a translation of x such that |x| < 1. In the MAPM library, two translations are performed first, and the series expansion is done in the third step. The first obvious translation is to limit the input angle to +/- pi radians. This is simply an x MOD 2*pi operation. Although this operation limits the input argument to +/- pi, pi is greater than one, so the series will still take too long to converge.
The second translation uses the multiple angle identity for sin(5x). This identity is:
sin(5x) = 16*sin<SUP>5</SUP>(x) - 20*sin<SUP>3</SUP>(x) + 5*sin(x)
If you can calculate sin(0.5), you could just use the above identity to compute sin(2.5). If you desire sin(3), calculate sin(0.6) and use the identity. By using this identity, the worst-case x input to the series expansion will be pi/5 or 0.6283. This is less than one, so the series will converge significantly faster than before. However, you can do better. To calculate sin(0.6283), you can use the multiple angle identity a second time to decrease the worst-case number to 0.6283/5 or 0.1257. You could (in theory) perform this operation numerous times, but you do reach the point of diminishing returns. The MAPM library uses this identity three times, so the worst-case input to the series expansion is pi/(5*5*5), or 0.0251. This is small enough so that the sine series does converge quite rapidly to the desired precision.
I have not seen the multiple angle identity translation used in this context before, so I believe the MAPM library is the first arbitrary precision library to use this technique.
The cosine function also uses a multiple angle identity to minimize the value of the number passed to the series expansion. The following algorithm summarizes the process. It is essentially the same in concept to the algorithm used for computing sin(x).
Step 1. Limit input argument to +/- pi.
Step 2. Use the multiple angle identity for cos(4x):
cos(4x) = 8*[cos<SUP>4</SUP>(x) - cos<SUP>2</SUP>(x)] + 1
Use this identity recursively three or four times as needed; that is, multiply the input angle by 1/(43) or 1/(44) respectively. If |x| is less than one radian, recurse three times. If |x| >= 1 radian, recurse four times.
This step yields a worst-case |x| = ~0.0156 (1/64 > pi/256).
Step 3. Apply traditional series expansion.
