Reverse-engineering the vintage Intel 8087's tangent algorithm: more than CORDIC

I hope you're not tired of the 8087, because I have another article about Intel's floating-point chip.1 In 1980, Intel introduced the 8087, making floating-point operations much faster in the IBM PC and other systems. In this article, I look at the algorithm behind the chip's tangent instruction. One popular approach for trigonometric functions is an algorithm called CORDIC. Another approach is a polynomial approximation. The 8087 combined the two to obtain both high accuracy and high performance. The 8087 provided an enormous speedup over the 8086 microprocessor, computing a tangent in 90 microseconds rather than 13,000 microseconds.2 By examining the circuitry and microcode of the 8087, I can explain the algorithm behind the tangent instruction, called FPTAN. To explore the 8087's circuitry, I popped the lid off a chip with a chisel and created a high-resolution image with a microscope. The microcode ROM is the large rectangular region in the center of the die, holding the 1648 micro-instructions that control the chip. The bottom half of the chip (red box) is the datapath, the circuitry that performs floating-point calculations on 80-bit values.3 A close-up of the 8087's datapath, showing functional blocks that are used by FPTAN. Click this image (or any other) for a larger version. Zooming in on the datapath shows the relevant functional units. The exponent ROM holds fixed exponent values that the algorithms need. The constant ROM holds constants, including the constants used by the CORDIC algorithm. The shifter is a large component; it shifts a 64-bit value left or right by arbitrary amounts. The adder is the heart of the 8087's calculations; as well as providing addition and subtraction, it is used in a loop for multiplication, division, and square roots. The B register holds one input to the adder, while multiple sources can provide the other input. The sum register holds the adder's output. The eight stack registers and the temporary registers hold floating-point numbers. Finally, the shift register holds 16 status bits for the CORDIC calculations. CORDIC is a clever algorithm for quickly computing transcendental functions with simple hardware: it uses shift and add instructions along with table lookups, but doesn't need multiplication or division. This algorithm dates back to 1956, when it was developed for the B-58 Hustler, the first bomber capable of flying at Mach 2. The aircraft had an analog navigation computer, but analog components provided limited accuracy. Engineer Jack Volder was given the task of designing a digital computer to replace the analog computer.7 One key problem was that an analog computer can easily generate sines and cosines with an electromechanical device called a resolver. But trigonometric functions are difficult to produce digitally, especially with the slow transistors of that era. A Convair B-58A Hustler, on display in San Antonio, TX (details). Jack Volder came up with a fast way to calculate trigonometric functions with simple hardware. He called the algorithm—and the computer that implemented it—CORDIC: "COordinate Rotation DIgital Computer". CORDIC converts an angle to a vector, where the vector's coordinates provide the necessary trig functions. The trick is to break down the angle into a sequence of special angles, angles that make vector rotation easy. These special angles are precomputed and stored in a table, so the CORDIC calculation can be performed quickly, even on 1950s hardware. Each CORDIC iteration provides an additional bit of accuracy, so the algorithm converges rapidly. CORDIC became popular, including in scientific calculators, which used decimal CORDIC instead of binary. Some trigonometry I'll try to keep the math to a minimum, but in this section I'll give a quick explanation of how CORDIC works. The diagram below reviews how trig functions are related to the coordinates of a point. Suppose you have an angle θ; it specifies a point (X, Y) on the unit circle. The basic formulas are X=cos θ, Y=sin θ, and Y/X = tan θ. Thus, if you can determine the coordinate (X, Y), then you can determine the value of the trig functions. If the point is not on the unit circle, e.g. (X', Y'), then you can still easily determine tan θ. (Spoiler: this is what the 8087 does.) However, sin θ and cos θ become messy.4 The relationship between an angle, the X and Y coordinates, and the trig functions. If you've done any computer graphics, you've probably seen how a rotation matrix can rotate a point by an angle. (If you're not familiar with rotation matrices, you can read about them here or just trust that it works.) Multiplying a point (X, Y) by the rotation matrix yields the new point (X', Y') as shown below. A point can be rotated by using a rotation matrix. Unfortunately, since the rotation matrix (1, below) requires sin and cos, it doesn't seem like it helps solve our problem. However, we can divide the matrix by cos θ; this seems even less helpful since now the matrix (2) needs tan, which is what we want to evaluate. (Moreover, the vector's length will grow.) But the key to CORDIC is to use special angles, αn = arctan(2-n). When we substitute one of these special angles into the matrix, we get matrix (3), which is easy to evaluate in hardware: multiplying by a power of 2 can be done by shifting the bits. Simplifying the rotation matrix. Applying matrix (3) to the point (X,Y) gives the equations (4). These are key equations for the CORDIC process. The important thing is that these equations are fast and easy to compute in machine language or hardware, as the only operations are addition, subtraction, and binary shifting. With that background, we can see how CORDIC works. First, we break down the desired input angle into a combination of special angles that adds up to the desired angle.5 Then we apply the rotation formula above for each special angle, starting with the unit vector (1, 0). The result is a point (X, Y) at the desired angle, and then the desired tangent is simply Y/X.6 Since the table of special angles is precomputed, the arctan operations don't slow down the process. As an aside, after the first few terms, the special angles approach 2-n, so they shrink by roughly a factor of 2 at each step. To summarize, the CORDIC algorithm consists of looping through a table of stored angles. If the stored angle is less than the desired angle, the stored angle is subtracted from the desired angle to yield a new desired angle and the equations above (shifts, add, and subtract) are applied to yield a new vector. At the end, the tangent of the original angle is given by Y/X. The rational polynomial approximation The accuracy of CORDIC depends on the number of terms that are used. With 16 terms, the accuracy is approximately 2-16, or 16 bits of accuracy. To get 64 bits of accuracy would require calculating 64 terms (and a table of 64 special angles). To get an answer faster, the 8087 uses 16 bits of CORDIC and uses another algorithm for the remaining angle. (The remaining angle is the gap between the sum of special CORDIC angles and the desired angle, so it is very small, around 2-16.) For the remaining angle, the 8087 uses a Padé approximant, which is the ratio of two polynomials. There's a whole family of Padé approximations, depending on the order of the polynomials. The 8087 uses a simple formula: 3x/(3-x2).8 Although this approximation is simple, it is very accurate for small values; its error is proportional to x4. Since x<2-16, the error will be less than 2-64, meeting the 64-bit accuracy requirement for the 8087. Moreover, the 8087 doesn't need to perform the division in the rational polynomial since FPTAN returns a separate numerator and denominator. Thus, the division is "free". The tangent function (red), rational approximation (blue), and Taylor series (green). Disclaimer: The Taylor series isn't as bad as it appears, since the relevant range is very close to 0. Graph generated with Desmos. If you've studied calculus, you might think that a Taylor series polynomial is the way to go, but the ratio of two polynomials is better. (One reason is that tangent blows up to infinity at π/2. A polynomial won't blow up, but the ratio of polynomials can, so it fits the tangent function better.) The graph above compares the tangent function (red), the rational approximation (blue), and the third-order Taylor series (green). Putting the pieces together: the 8087 algorithm The 8087's tangent algorithm has three parts: determining the CORDIC decision bits (called pseudo-division), computing the rational approximation, and applying the rotation equations based on the CORDIC decision bits (called pseudo-multiplication).9 In more detail, the first step determines which special angles to add to approximate the input angle. Each special angle is compared to the remaining input angle and subtracted if it is smaller. If the angle is subtracted, a 1 is recorded; otherwise, a 0 is recorded. Since this process is similar to how long division subtracts (or doesn't subtract) successive shifted versions of the divisor, generating 1s or 0s for the quotient, the process is known as pseudo-division. Note that the rotations are not applied in this step. Instead, this step decides which rotations to apply later. The diagram below shows this process applied to the input angle 0.95 radians.10 The process generates the sequence of bits [1,0,0,1,0,1,0,1,0,0,1,0,0,1,1,1], where the leftmost bit indicates arctan(20) and so forth. The angle is reduced by roughly a factor of 2 at each step, so the residual angle is very small. In the first phase of the CORDIC algorithm—pseudo-division—the input angle is reduced by special angles, leaving a residual angle at the end. ("rad" is radians, not the unit of radiation.) Next, the tangent of the remaining angle is calculated with the rational approximation function, 3x/(3-x2). The result is used as the initial vector for the next step. The division is not