How DO you do Sine/Cosine/Log etc algorithmically?

Factual Questions
1

In a recent programming thread somebody brought up the point that we don’t know how to do logs by hand, even if we understand what they mean. I however, want to know, how do you do the basic trigonometric functions (sine, cos, tan, asin, acos, atan) by hand, and how do you do log[sub]x/sub by hand?

Obviously most people used a slide rule or table, but the table had to come from somewhere, and I’d be relatively surprised if computers used vast lookup tables under the hood for this stuff if there was a decent algorithmic solution available.

2

It is an interesting story. Here are some links you can use to start investigating.

CORDIC

Logarithm algorithm

3

Back in the day, the various functions and algorithms created by Pafnuty Chebyshev were commonly used. Many new methods and adaptations are in use now.

4

For natural logarithms (logarithms in any other base can then be calculated using log[sub]y/sub = ln(x)/ln(y)):

Well, there are a number of ways, but one is by using the fact that ln(x) = 2 * (b[sup]1[/sup]/1 + b[sup]3[/sup]/3 + b[sup]5[/sup]/5 + b[sup]7[/sup]/7 + …), where b = (x - 1)/(x + 1), along with the fact that the remainder when truncating this series is always less than 1/(1 - b[sup]2[/sup]) times the (magnitude of the) first omitted term. Thus, the series produces accurate digits at a rate of at least 2 * log[sub]10[sup]-1[/sup]/sub digits per term.

For example, for x = 2, this tells us that ln(2) = 2 * (3[sup]-1[/sup]/1 + 3[sup]-3[/sup]/3 + 3[sup]-5[/sup]/5 + 3[sup]-7[/sup]/7 + …), with a remainder at any moment less than 9/8 times the next term, and thus producing accurate digits at a rate of at least 2 * log[sub]10/sub (a shade over 0.95) digits per term.

If we use 10 terms (up through 2 * 3[sup]-19[/sup]/19), we find that ln(2) = 0.69314718054etc. + a remainder of less than 10[sup]-11[/sup]. Thus, we obtain 10 accurate digits in 10 terms. Of course, the rational number arithmetic here is excessively tedious to carry out by hand; thank god for computers.

[In case you are curious, as I assume you are, the series for ln given above comes from integrating both sides of the geometric series equation 1/(1 + b) + 1/(1 - b) = 2/(1 - b[sup]2[/sup]) = 2 * (b[sup]0[/sup] + b[sup]2[/sup] + b[sup]4[/sup] + b[sup]6[/sup] + …), and the upper bound on the remainder comes from comparison to the same geometric series]
I should also note, this means the series converges faster the larger b[sup]-1[/sup] is (equivalently, the closer x is to 1).

Thus, for truly optimizing efficiency, one can often use the trick of rewriting x as a product of powers of numbers close to 1, and finding ln(x) as the corresponding linear combination of natural logarithms of numbers close to 1, for which the series produces digits at a faster rate.

For example, 2 = (7/5)[sup]2[/sup] * (100/98), so ln(2) = 2 * ln(7/5) + ln(100/98), which means ln(2) = 2P(6) + P(99) [where P(d) = 2 * (d[sup]-1[/sup]/1 + d[sup]-3[/sup]/3 + d[sup]-5[/sup]/5 + …) = ln((d + 1)/(d - 1)), producing at least 2log[sub]10/sub newly accurate digits per term].

Or, if you’re more hardcore, you can use the fact that 2 = (27/25)[sup]9[/sup] * (4802/4800)[sup]-1[/sup] * (8750/8748)[sup]4[/sup], so ln(2) = 9P(26) - P(4801) + 4P(8749). And 2 = (252/250)[sup]72[/sup] * (450/448)[sup]27[/sup] + (4802/4800)[sup]-19[/sup] * (8750/8748)[sup]31[/sup], so ln(2) = 72P(251) + 27P(449) − 19P(4801) + 31P(8749). Decompositions like this are the key to pushing computation of the decimal expansion of ln(2) to record lengths (though the efficiency gains don’t really justify the complication for anything other than record-setting; the last series produces just under 1.2 digits of precision per term (counting adding a new term to each series as adding four new terms altogether), which is just under 0.25 digits per term faster than our original straightforward calculation of ln(2) as P(3); it produces 10 digits in 8 terms, rather than 10 terms, which is no big whoop. But when you’re calculating millions of digits, the difference becomes rather more pronounced).

Finally, I should also note that our series and remainder bound works for complex arguments as well, from which the inverse trigonometric functions can be computed. (For example, the arctangent of x can be computed as the imaginary component of the natural logarithm of 1 + ix)

5

Oh, I should clarify, this works when the magnitude of b is less than 1. This covers every positive x, and, indeed, every complex number with a positive “real” component.

6

Or rather, you could use the formula arctan(x) = ln((1 + ix)/(1 - ix))/2, which amounts to setting b to ix and removing the factor of 2 from our ln series (this will work, with the same error bound as before, whenever x has magnitude less than 1).

7

For cosine and sine, you can use the fact that e[sup]ix[/sup] = (ix)[sup]0[/sup]/0! + (ix)[sup]1[/sup]/1! + (ix)[sup]2[/sup]/2! + (ix)[sup]3[/sup]/3! + …, and then extract cos(x) and sin(x) as the “real” and “imaginary” components of e[sup]ix[/sup]. You need error bounds to know when to stop, though: cutting the series off after finitely many terms produces an error that’s no bigger (for real x) than the first omitted term. So the series converges geometrically, with smaller x yielding tighter convergence, which you can potentially manipulate to your advantage as needed (computing e[sup]ix[/sup] for large x by taking powers of e[sup]ix/n[/sup] for large n instead).

(This error bound arises from Taylor’s theorem, though there may be a nicer way to derive it than that.)

8

Oh, yeah, there is… just compare the cut-off portion of the series to the first omitted term times the entire series; the latter is at least as large at each term. Thus, the error is no larger than the first omitted term times the entire series; but the entire series has magnitude 1 (since it computes the rotation e[sup]ix[/sup]), leaving us with an error bound of the magnitude of the first omitted term.

As for where the series comes from, it is constructed to be equal to i times its own derivative and have initial value 1, and thus must compute e[sup]ix[/sup] within its radius of convergence, which is infinite.

9

Er, not “geometrically”, what am I saying? Much faster than that… :smack:

10

I once had to implement sine, cosine, tangent, etc. for a cellphone which had no floating point processor, with the end goal of doing 3D rendering on it. I used lookup tables. But, a cellphone’s display at the time was always less than say 200 pixels in either direction, so with tables which only returned down to the half degree, the answers were close enough that it appeared correctly. (Technically, I ‘invented’ my own measurement rather than using radians or degrees, using 512 units to make a full circle as this was sufficient for my needs and a power of 2, which can always help.)

I wouldn’t be surprised if some ancient machine did something similar, to be honest. Lookup tables are just tremendously fast and only need to be as large as you need to suite your level of precision. My tables were only 2048 bytes in size, which isn’t all that huge.

11

I did something similar for a low res display on a weak processor. It only needed a table to define a quarter of a circle, and a scaling table since the processor didn’t have multiply/divide. For a 512x512 display only 512 byte were necessary. But only circles and ellipses needed to be displayed. I suppose high speed rendering hardware use some tables to speed up some functions.

12

Don’t most commercial processors still use lookup tables and interpolation to implement their floating point hardware? I’m pretty sure the FDIV big, for instance, was caused by a faulty lookup table, with missing entries.

13

If anyone remembers the days of 386’s, 387’s, and 486’s. The 486 had the floating point math coprocessor embedded, while the 386 had the option of adding a 387 (which tremendously sped up the math.) If you did not have a 387, there was a floating point library executed using the integer functions of the 386. Sin, cos, log, sqrt, etc. were calculated with algorithms like the ones mentioned above.

I saw a piece of software being sold once that was “make your 386 a 486”, Essentially what it did was replace the floating point DLL with a library of routines that only calculated to 3 digits of accuracy. In the days when 800x600 was excellent screen resolution, this was enough to speed up a lot of games that relied on floating point math. Good for Doom. Probably not so good for scientists, I suppose.

Those extra few digits (typically at the time 8 digits IIRC) made a big difference in speed.

14

This is a very common way to implement sine/cos etc. in digital logic as well. If you only need 8 or nine bits of accuracy it is much faster and smaller in hardware to have a look up table instead of multipliers and adders. Look up tables are also used for taking logarithms. Logs are used a lot in digital communications to create the log likelyhood ratios for use in the error correction algorithms. The log approximations tend to be really approximate since the data is so noisy there is not much gain in getting the logs accurate.

15

Going way back, Napier had an extremely elegant way of calculating 7 place natural logs by hand. First, he conceived of the function described by a moving point that moves towards the origin in such a way that its speed is inversely proportional to its distance from the origin. He actually had it moving along the negative x axis starting from -1, so he actually did ln(-x), but that is just a detail. He began with 1 and subtracted 10^-7 from it to get ln(1-10^-7) and then multiplied by 10^-7 and subtracted again. The point is that he changed the problem of multiplying to moving 7 places and adding. This worked well enough till he got to 1 - 10^-6. He then used the same method to go down to 10^-5 and then interpolated in between. He continued, alternating this easy multiplication and interpolation. He got down to -1/2 and then used the characteristic property of logs to get to -1/4 and so on.

Later on Briggs came along and pointed out how much more useful base 10 logs would be. Napier agreed and the two of them began to prepare tables of common logs (what method they used I don’t know; perhaps they used log(x) = ln(x)/ln(10) and used their already existing ln tables to reduce it to subtraction). At any rate, Napier died before they finished and Briggs finished the job. Amazing work, all by hand of course.

16

IIRC, I saw an article in the early days of computers that analyzed log tables published in the early 1800’s or late 1700’s and found a typo in quite a few of the entries around the 5th decimal place.

17

If you use base 2 logs, then it only takes a few table entries to get reasonable precision, as you can just count how many bit shifts it takes to scale the argument into your table. Like how you find the integer part of a base 10 log by shifting the decimal point. But then you need a table that goes from 0 to 10 by whatever precision. The log base 2 table only has to go from 0 to 2, so it is 1/5 the size for a given precision.

18

This was my impression (not sure exactly why), too. Though a sophisticated interpolation might end up being pretty much the same thing as one of these algorithms, just with a look-up table to give you a head start.

19

Perhaps its also worth mentioning that people tended to keep the functions in the equations until as late as possible so they’d get reduced or cancel out.

20

The reason why is speed. The old fashioned way of doing a divide without a lookup table is to break it down into a loop of subtract and shift operations. Each pass through the loop takes a clock cycle, so these end up being fairly slow. A lookup table is a single clock cycle operation.