High precision functions in Basic. (part 1) Albert Nijenhuis.
High Precision Functions In Basic
Part 1
Many versions of Basic supply high precision variables and the four operations +, -, *, / for them, but few if any also have high precision versions of SQR, LOG, and similar functions. Some versions of Basic do not supply these functions at all.
Until recently, I gave little thought to these matters, but that changed quite unexpectedly when I tried to duplicate a series of calculations which were supposed to prove that if present trends continue, by the year 2000 the nation's college teachers will have reached the poverty level.
The I/O data were in the range 1000 to 50,000 so I expected no problems with single precision. Still, I could not get a good duplication of the findings, which involved nothing more than a case of fitting data to an exponential curve. At first I looked all over for a cause, but finally realized that I had a "precision problem' when I found that my pocket calculator, which carries about 10 digits, got different results again.
This article and Part 2 show how a few formulas of the calculus, with which many readers are familiar and which others may be willing to accept as given, can rather directly be used to obtain all the functions LOG, EXP, SIN, COS, ASN, ATN, and SQR with high precision.
If you just want a "how' for EXP, SIN, and COS, you may simply copy the programs. In what follows we discuss the "why' for these functions.
The programs have been written on a TRS-80 Model I Level II using double precision and should easily adapt themselves to other situations. This does not hold for the program in Listing 3, which makes fairly heavy use of machine-dependent features.
EXP
To get EXP, we use the infinite series expansion for e(x), which is derived in calculus:
(1)
e(x) = 1 + x/1! +x2/2! +x3/3! + . . .
Although this is an infinite sum, which would take an infinite amount of time to compute for any given value of x, we can get useful approximations--and that is all we want--by summing a finite number of terms.
In particular, if the values of x are relatively small, a few terms will already give quite good values. This is because the factorials in the denominators increase very rapidly. For example, 3!=1 2 3 = 6, 5!=120, 10!=3.6 10(6), 20!=2.4 10(18).
In the program that we are about to develop, the values for x in (1) will always lie between -0.5 and 0.5. A little experimentation shows that the terms through x16 will suffice: (.5)(16)/16!=7.3 10(-19); plenty for 16-place accuracy in the calculation of e(.5) or e-.5.
To evaluate the sum through x16 we use nested multiplication (also known as Horner's method). To illustrate the idea, here is how you handle the sum through x4:
1 + x(1 + x/2(1 + x/3(1 + x/4))) with evaluation starting with the innermost parentheses. This method has the advantages of easy programming, no underflow, and minimal round-off error.
The program will require the value of e=2.718 . . . to full achievable precision. We could copy this from a table, but why not calculate e.5 from the formula and square it?
To calculate e(x) for arbitrary values of x, determine the integer n so that -.5 < x - n < .5. We then multiply e(x-n), obtained from (1), by e(n), obtained by multiplying repeatedly by e (if n is negative, we multiply by 1/e).
Listing 1 shows a Basic program that does these calculations. All variables are double precision except the loop parameters I and N, which are integers.
To obtain good single precision, the loop parameter 16 in line 1310 may be replaced by a smaller value. I found that 9 gave satisfactory results. To test your best=smallest value, experiment: if 9 and 10 give the same result, but 8 does not, then use 9.
SIN and COS
The method to calculate SIN and COS is very similar. In this case we use the infinite series, in which x is expressed in radian units:
(2) sin x = x - x3/3! + x5/5! - x7/7! + . . . cos x = 1 - x2/2! + x4/4! - x6/6! + . . .
The values of x for which we now demand high accuracy lie between - /4 and /4. By experimentation we find that terms through x18 will suffice, because ( /4)18/18! = 2 10(-18), again very safe for 16-digit accuracy. Of course, these sums are also computed with the use of nested multiplication.
To calculate sin x or cos x for arbitrary x, first determine an integer in such that - /4 < x - n /2 &s /4, and use the formulas to calculate sin(x-n /2) and cos(x-n /2). Then use standard trig formulas to find sin x and cos x. For example sin x = cos(x- /2), or sin x = - sin(x- ). The formula to be used depends on the value of n modulo 4. The value of n + 1 (mod 4) denotes the quadrant in which x + /4 lies.
The accurate value of that is needed to perform the reductions of the input variable x, can be entered directly as part of the program. There is something interesting to be learned, however, from a calculation by the program. For example, this may raise the question again as to what actually is, and what radian units (rather than degree measure) are all about.
Figure 1 shows a circle of radius 1, in which an angle AOB of size x has been drawn. Because radian units were used, the length of the arc AB subtended by the angle is exactly x units. Parallel to OA is line DB, where D lies on the line OC perpendicular to OA. As a result, the distance DB is cos x.
We want to calculate the length of the arc ABC, i.e., we want to determine /2, thereby, of course, determining itself. Our basic available tools are algebra and the formulas for sin and cos. Now, if x is not too small, and the angle AOB is acute, then x may be a fairly reasonable estimate of the length of arc ABC. A much better estimate is, however, x + cos x, which is the combined length of arc AB and line segment BD.
All of this may not make too much sense to the casual reader, but there is system to this madness. Namely, if going from x to x + cos x gives an improvement in the estimate of /2, then applying the same method of improvement again to the number x + cos x will give an even better estimate of /2. This process can be repreated as often as we wish.
What is amazing about this all is that, if we start with x = 1 (not a particularly good approximation for /2) and do this improvement business just four times, we have obtained /2 with maximum achievable accuracy in double precision!
Equally interesting is that, starting from the approximate value 355/133 of , just one application of the improvement method gives to full double precision.
The execution of the scheme must be done with a cos function that takes input x of size up to /2, but the finite part of the series that was chosen is accurate only up to /4. We therefore use the doubling formula
cos x = cos2(x/2) - sin2(x/2).
This has been included in the program in Listing 2, which calculates sin and cos, in addition to determining .
For single precision, the value 18 in line 2310 may be reduced by experimentation. The value must, however, be even to give meaningful results.
His, then, describes the calculations in this article. In Listing 3 you will find another copy of the combined programs, this time tuned in detail to the TRS-80 Model I Level II. All auxiliary variables start with the same character to simplify the task of avoiding their use elsewhere.
Although I used a V for readability, I recommend an O because of its unpopularity with programmers. It should be possible to MERGE this program with others and thus avoid repeated entry through the keyboard.
The methods to calculate EXP, SIN, and COS shown here are not the fastest known in the trade, but they are simple and quite effective. If you wish to make improvements, there are many things you can do. The single most effective improvement is, of course, to use a compiled language or to code directly in machine language.
Table: Listing 1.
Table: Listing 2.
Table: Listing 3.
Photo: Figure 1.