High precision functions in Basic; part 2. Albert Nijenhuis.
High Precision Functions in Basic
In Part 1 (Creative Computing, December 1983) we provided Basic routines that would enable us to calculate EXP, SIN, and COS with high precision. In Part 2 we deal with SQR, LOG, ASN, and ANT. These functions are the inverses of the squaring function, of EXP, SIN, and TAN=SIN/COS, respectively.
Calculating SQR(Y) means finding a solution x of the equation x2 = y. Assuming y is positive, there are two solutions, of opposite signs; we want the positive one. (If y is negative, there is no solution at all.)
To find LOG(Y) we solve the unknown x from the equation e(x) = y; if y is positive, there is one solution, otherwise none.
Similarly, calculating ATN(Y) means finding a solution x of the equation tan x = y. Again, there is an infinite number of solutions, and we want the one between - /2 and /2, in radian measure.
Finally, calculating ASN(Y) means finding a solution x of the equation sin x = y. If y lies between -1 and 1, then there is an infinite number of solutions; we want the one that lies between - /2 and /2. (For other values of y there are no solutions.) We will reduce the calculation of ASN to that of ATN via the formula ASN(Y) = ATN(Y/SQR(1-Y*Y)).
The method by which we shall calculate inverse functions is a level more sophisticated than the infinite series we used in Part 1. This time the values are calculated by a method of successive approximations of a less transparent type. To appreciate these mothods it is useful to be familiar with what is known as Newton-Raphson.
In fact, our method for SQR is just that, but the methods for LOG and ATN are more delicate and converge even faster; that is, they require fewer steps to achieve the desired accuracy. It may be useful to have a text on calculus or "pre-cal' handy to check some of the formulas.
We shall use the infinite series expansions as follows:
In (t) = (t-1) - (t-1)2/2 + (t-1)3/3 - . . . arctan (t) = t - t3/3 + t5/5 - t7/7 + . . ..
Unfortunately, these series aren't nearly as nice as those in Part 1. First of all, they are not valid for all values of t: for ln (t), t must lie between 0 and 2, for arctan (t), t must lie between -1 and 1. Worse, however, unless t is quite close to the middle value (1 for 1n (t), 0 for arctan (t)) many terms are required to obtain sufficient accuracy. All this is due to the absence of the rapidly increasing factorials in the denominators that we had before.
The required high accuracy will therefore be obtained by referring back to the high accuracy versions of the functions of which these are the inverses, together with the rather crude approximations obtained from the series above by using only two terms:
The basic idea is to start with an approximation to the value of x in the equation to be solved, to use the high accuracy function to see how close you are to the give y-value, and to use the result to obtain a better approximation. This process is repeated as needed, each time adding a small amount to the previous approximation to get a better one.
Calculating the Logarithm: Scaling
Before we start the successive approximations, we scale the problem to a workable range.
As an illustration, we calculate 1n (12) to 3 places, using a 10-place calculator equipped with e(x). There is no question of using the series for 1n (t), of course. What we can do, however, is note that 1n (12) is one unit bigger than the 1n of 12/e = 4.416, and similarly, two units bigger than the 1n of 12/e2 = 4.146/e = 1.624, three units bigger than the 1n of 12/e3 = 1.146/e = .597, and so on. The scalling we choose in this case consists of finding a number between e-.5 = .607 and e.5 = 1.649. That corresponds to the range -.5 to .5 for which the first 16 terms of the series for e(x) in Part 1 give high accuracy.
Thus, we see that 1n (12) = 2 + 1n (1.624), and now we must find 1n (1.624). At this point we use the two-term formula and find .429 as an approximation to 1n (1.624). To see how close this is, we calculate (with high precision) e.429 = 1.536. The closeness is measured by 1.624/1.536 = 1.057. Therefore, 1n (1.624) = 1n (1.536) + 1n (1.624/1.536) = .429 + 1n (1.057). Repeat with 1.057; its 1n is approximately .056, and the exponential of this (with high precision) is 1.057. These calculations thus give with three-digit accuracy than 1n (1.057) = .056. Totaling all terms we find the answer: 1n (12) = 2 + .429 + 0.56 = 2.485
Program for Log
The program in Listing 1 is a close representation of the method. The input variable is X (not y), and the constants E2, and so on are set by the EXP program in Part 1. To avoid double precision divisions where possible, the program calculates, as an example, 1.624 e-.429 instead of 1.624/e.429; thus, X1 in line 1560 is the opposite sign of what was expected. The approximations are performed four times. This figure was obtained by experimentation on the least favorable input values, e-.5 and e.5.
Calculating The Arctan
Again, we illustrate the calculation of arctan with a hand calculator example: find arctan (2). Scaling leads to
x = /2--arctan (.5).
The two-term approximation yields .458 for arctan (.5). We denote this value a. Calculated with high precision, tan (.458) = .493. Let us denote this value by b, then tan (a) = b, and a is an approximation to arctan (.5).
We now want a formula of the form arctan (.5) = a + arctan (something), so we can repeat the approximation process on the second term. (This in analogy to the formula for 1n (1.624) above.) To find what "something' is, subtract a from both sides, take the tan and obtain tan (arctan .5 - a) = "something.' On the left side we now use the trig formula: tan(u-v) = (tan (u) - tan (v)/(1 + tan (u) tan (v)) in which we set u = arctan (.5) and v = a, which amounts to setting tan (u) = .5 and tan (v) = b = .493, and we obtain:
"something' = (.5 - .493)/(1 + (.5) (.493)) = .005
That, is arctan (.5) = .458 + arctan (.005). The same procedure (left to you) applied to arctan (.005) gives the value .005 which is accurate to three places. Hence, arctan (2) = 1.571 - .458 - .005 = 1.107 (calculation in 10 digits, with three showing).
Programs for Arctan
The program in Listing 2 should not be hard to follow. The input is X, and lines 2800 to 2850 do the scaling. Lines 2900 to 2960 contain the approximation loop, with 2910 containing the two-term formula and 2950 the trig formula.
Calculating the Arcsin
The calculation of arcsin (y) follows a similar patterns. Let x be the desired angle, so sin (x) = y and x lies between - /2 and /2, whence cos (x) > 0. By virtue of the popular formula sin2 + cos2 = 1 it follows that cos x is the square root of l - y2. The square root must be calculated with high precision by a routine which we discuss below. Let z be the value of the square root. We then distinguish three cases for scaling: y > z, y < -z, and y between -z and z. In the first two cases, we calculate arctan(z/y) and use the same complementations as for ATN; in the last case we calculate arctan(y/z).
The program is given in Listing 3 and should require no further explanation.
The methods for LOG and ATN may be modified in several ways; for example, the two-term approximations may be replaced by some other approximation, such as a one-term approximation, or an approximation of degree 3 or higher. The choice of approximation method will determine the number of times the approximation must be performed for high accuracy. Sometimes an initial approximate value of surprisingly high accuracy can be obtained, such as the value of a single precision LOG or ATN function. We use this approach in the program in Listing 5, which is to be merged with the program in Listing 3 of Part 1. Then just one two-term approximation step will suffice.
Calculating Square Roots
A high precision version of SQR is obtained by the famous formula: if y is a given positive number, if x is the square root of y and if t is any positive number, then 1) if t<x then (t + y/t)/2 is greater than x; 2) if t>x then (t+y/t)/2 is also greater than x but is much closer to x. Hence, for given y, set t = 1, calculate (t + y/t)/2, reset t to this value, and repeat until satisfied.
The number of repetitions required depends on the value of y. If y is large, or close to zero, more repetitions are needed than if y lies, for example between .5 and 2, when five times is enough. We have, therefore, chosen to change y to a number between .5 and 2 by dividing (or multiplying) by factors 4; the square root thereby gains (loses) an equal number of factors 2. The program is given in Listing 4.
A better initial choice for t is the single precision value of SQR(Y) if that is available. Then no scaling is needed; this method is used in Listing 5.
Table: Listing 1.
Table: Listing 2.
Table: Listing 3.
Table: Listing 4.
Table: Listing 5.