!pr1
18-Digit Arithmetic, Part 8................Bob Sander-Cederlof

Someone pointed out last week that this series is getting a little long.  Well, we are nearing the end.  What we are doing is probably unprecedented in the industry:  listing the source code and explaining it for a large commercially valuable software product.  It takes time and space to break precedents.

This month's installment completes the normal set of math functions, with sine, cosine, and arc tangent.  We even slipped in a simple form of the tangent function.  Still to come are the formatted INPUT and PRINT routines.

Some Elementary Info:

Trigonometry is a frightening word.  (If it doesn't scare you, skip ahead several paragraphs.)  The "-ometry" refers to measurement, but what is a "trigon".  Believe it or not, "trigon" is another name for a triangle.  Trigon means three sides, and figures with three sides just happen to also have three angles.  "Trig" (a nice nickname) is a branch of mathematics dealing with triangles, without which we could not fly to the moon, draw a map, or build bridges.  Strangely enough, much of electronics also uses trig funtions ... are electrons triangular?

When I took trig in high school, long before the day of personal calculators, we used trig tables.  (These were not articles of furniture made in the local woodshop, but rather long lists of strange numbers printed and bound into books.)  The tables contained values for various ratios of the sides of a triangle having one 90-degree angle.  Now we use calculators or computers, but obviously the trig tables would not fit in them.  Instead, approximation formulas are used.

In high school, we talked about six different ratios:  sine, cosine, tangent, cotangent, secant, and cosecant.  When it is all boiled down, we really only need the sine; all the rest are derivable from those.  The sine function gives a a number for any angle.  We frequently need to be able to go from a trig value back to an angle, and the most useful function for that is called the inverse tangent, or arctangent.

Even though I have been talking about triangles, trig functions are even more related to circles.  We compute functions of the angle between any two radii, like the hands on an old-fashioned, pre-digital wrist watch.  When we start talking about circles, we get into radians vs. degrees.

Just as scientists like logarithms to the base e (rather than 10), they also like trig functions based on angles expressed in radians, rather than degrees.  Degrees were invented back in Babylon, I understand, and are nice and clean:  360 make a complete circle.  Radians are not clean:  360 degrees is two-times-pi radians.  Nevertheless, many physical and electronic formulas simplify when angles are expressed in radians.  Consequently, calculators and computer languages usually expect your angles to be expressed in radians.  Some allow both options.  Applesoft expects radians, and so do my DP18 programs.

We commonly think of an angle as being somewhere between 0 and 360 degrees, or the equivalent range in radians.  However, angles can actually be any number, from -infinity to +infinity.  The numbers beyond one complete circle are valid, but they don't buy much.  If you stand in one place and spin around 1445 degrees (4*360 + 5) you will end up pointing the same direction as if you merely swiveled 5 degrees.  Therefore the first step in a sine function calculation involves subtracting out all the multiples of a full circle from the angle.

The arctangent function could return an infinite number of answers, but that is impractical.  We will return only the principal value, which is the one closest to 0.  All others are that value plus or minus any number of full circles.  In DP18 the ATN function may have one or two arguments.  If you only have one argument, the result will be an angle between -pi/2 and +pi/2.  If you specify two arguments, a value between -pi and +pi will be returned.

The Nitty-Gritty:

Enough of this preliminary stuff, let's get into the code.  In the listing which follows, you will find entries for four functions:  SIN, COS, TAN, and ATN.

Perhaps the easiest is the TAN function, at lines 2530-2630.  Since tan=sin/cos, that is all this code does.  We lose a little speed and possibly some precision with this simplistic solution, but the TAN function is relatively rarely called.

Next in difficulty is the COS function, lines 1630-1710.  Since cos(-x)=cos(x), we start by making the sign positive (lines 1690-1700.  Since cos(x)=sin(x+pi/2), we add pi/2 and fall into the SIN function.  Simple, but effective.

The SIN function gets more interesting.  For very very small angles, within the precision of 20 digits, sin(x)=x.  Lines 1780-1810 check for exponents below -10; all angles smaller than 10^-10 are small enough that sin(x)=x.

Next we take advantage of the fact that sin(-x)=-sin(x), at lines 1820-1860.  We remember the sign by shoving it on the stack, and force the sign of x positive.

Lines 1870-1950 get the principal angle.  I divide x by twopi, and throw away the integral part.  The fractional part that remains is a fraction of a full circle, a value between 0 and .999999...9 (not radians, and not degrees either).  Note that if x was extremely large there will be no fractional part, and the remainder will be zero.  Some SIN function calculators give an error message when this happens, but I chose to let it ride.

Lines 1960-2000 multiply the circle-fraction by four.  This gives a number between 0 and 3.99999...9, which I will refer to later as the "circle fraction times four", or c-f-t-f.  The integer part is effectively a quadrant number, and the fractional part a fraction within the quadrant:
                 |
               1 | 0
             ---------
               2 | 3
                 |

Lines 2010-2030 determine if the angle is in the first (0) quadrant.  If so, no folding need be done.

Lines 2040-2070 determine if the angle is in the second (1) quadrant.  If so, we skip ahead to apply the fact that sin(pi/2 + x) = sin(pi/2 - x).

Lines 2080-2160 are executed if the angle is in the 3rd or 4th quadrants (integral part is 2 or 3).  Here I apply the fact that sin(pi+x)=-sin(x).  I pull the saved sign off the stack, complement it, and shove it back on (lines 2090-2110).  Then I subtract 2 from the c-f-t-f, yielding a number between 0 and 1.99999...9.  We have folded the third and fourth quadrants over the first and second quadrants.  Next lines 2170-2190 determine if the result was in the first quadrant or not.

Lines 2200-2240 fold a second quadrant number into the first quadrant, by applying the fact that sin(pi/2+x) = sin(pi/2-x).  Subtacting the c-f-t-f from 2 flips us into the first quadrant.

Lines 2260-2270 pull the sign off the stack and make it the sign of the angle.  Remember that now the angle is a fraction (between 0 and .99999...9) of a quadrant.  After all these folding operations, the angle might again be very very small, so lines 2280-2300 check for that possibility.  If so, sin(x)=x, but that is only true when x is in radians.  Lines 2490-2520 convert the quadrant-fraction to radians by multiplying by pi/2, and exits.

Lines 2310-2470 handle larger angles by computing x*P/Q, where P and Q are polynomials in x^2.  The constants for P and Q are given in lines 1420-1550, and come from the Hart book.  [ I should mention here that I wrote those constants with pretty periods separating groups of five digits.  This will not assemble in some older versions of the S-C Macro Assembler.  If you get a syntax error, just leave out the periods. ]

Turning the Tables:

ATN is hardest to compute.  First we have to deal with the two variants of calls, having one or two arguments.  While all the previous function programs were called with the argument already in DAC, DP.ATN is called immediately after parsing the ATN-token.  Lines 2960-3070 parse and process the following parentheses and whatever is between them.

Lines 2960-2970 require an opening parenthesis.  Line 3070 requires the closing parenthesis.  In between we expect one expression, or two expressions separated by a comma.  If there is only one, we fake a second one (= 1.0).

What are the two arguments?  Looking at a cartesian system, with the vector shown below, the arguments are (Y,X).  If you call with one argument, it is (Y/X).

               y
               |
               |   .
               |
             --------x
               |

By using two separate arguments, rather than just the ratio, we can tell which of the four quadrants the vector was in.  DP.ATAN will return a value between -pi and +pi, depending on the two signs.  If you specify only the ratio, DP.ATAN will return a value between 0 and +pi depending on the sign.

Lines 3120-3160 save the two signs in bits 6 and 7 of UV.SIGN.  Way at the end, lines 4100 and following, UV.SIGN determines the final value.  If the sign of the denominator (X-vector) was negative, the composite vector is in the 2nd or 3rd quadrant:  computing pi - angle gives a result between pi/2 and pi.

If the numerator (Y-vector) was negative, the composite vector is in the 3rd or 4th quadrant.  Flipping the sign gives a result between 0 and -pi.

Lines 3180-3220 check for special cases of Y=0 or X=0.  If the first argument (Y-vector) is zero, the angle is 0 or pi depending on the sign of the second argument.  If the second argument (X-vector) is zero, the angle is either +pi/2 or -pi/2, depending on the sign of the first argument.  What if both arguments are zero?  That should produce an error message, but I am overlooking it:  I will return an angle of 0 in this case.

If neither argument is zero, some special checks are made to see if the value of the ratio is very small or very large.  I check before actually dividing, so the divide routine won't kick out on an overflow error.  If the ratio would be greater than 10^20, I return a value of pi/2.  This is accurate within the precision of DP18.  On the other hand, if the ratio is smaller than 10^-63 I return 0.  If neither extreme is true, I go ahead an divide to get the actual ratio.  Then I check for an extremely small ratio, in which case atan(x)=x.

If we find our way down to line 3390, the ratio is between 10^-10 and 10^20.  That is still too large a range for comfort, so we apply the fact that atan(1/x) = atan(pi/2 - x).  If the ratio of Y/X is greater than 1.0, then we take the reciprocal and remember that we did so.  This in effect folds the range at pi/4.  The resulting argument range is between 10^-10 and 1.  The variable N holds either 0 or 2 as a flag:  0 if we were already under 1, 2 if we formed the reciprocal.

The shape of the curve of the arctangent function between 0 and 1 (an angle between 0 and pi/4) is deceptive.  It looks nice and easy, but a polynomial over that range with 20 digits of precision is much too long.  We can easily reduce the range still further by applying another identity.  If the reduced argument is now already below tan(pi/12), fine.  If not, calculating (x*sqr(3)-1) / (sqr(3)+x) will bring it into that range.  If we have to apply that formula, N will be incremented (making it 1 or 3).

The curve between 0 and tan(pi/12) looks almost like a straight line to the naked eye, but it really is far from straight.  It takes a ratio of the form P/xQ where P and Q are polynomials in x^2.  The coefficients are given in lines 2650-2770, again from Hart.  The ratio is computed in lines 3800-3960.

Lines 3970-4080 start the unfolding process.  The variable N is either 0, 1, 2, or 3 by this time.  If N is 0, no folding was done.  If N is 1, only folding above pi/12 was done.  If N is 2, only folding above pi/4 was done.  If N is 3, both folds were done.  These lines convert the angle back to the correct value, using a table of addends and an optional sign flip:

       N       unfolding formula
     ---------------------------------------
       0       none
       1       pi/6 + x
       2       pi/2 - x
       3       pi/2 - ( pi/6 + x) = pi/3 - x

That's it!  We already discussed the code beyond line 4100, which figures out which quadrant the angle is in.

Any questions?
