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

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 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).

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 performed here; instead, the numerator becomes Y in the initial vector and the denominator becomes X. Thus, the expensive division step is avoided, since it happens implicitly in the answer. Multiplying by 3 is easy (shift left and add), so the only expensive operation at this step is squaring the angle, which requires a full 64-bit multiplication.

The final step applies each CORDIC rotation if the corresponding decision bit from the first step is 1. Since this is somewhat analogous to binary multiplication, which adds the multiplicand at each step if the multiplier bit is 1, this step is known as pseudo-multiplication. As described earlier, each rotation is computed with shifts, addition, and subtraction, so each rotation is inexpensive. Each rotation in this step corresponds to an angle reduction in the first step. The rotations are applied in reverse order—the smallest one first—to reduce rounding error. To accomplish this, the decision bits are stored in a 16-bit shift register in the first step, and shifted out in opposite order in this step.

In the last phase of the CORDIC algorithm—pseudo-multiplication—a vector is rotated multiple times by applying shifts and adds. The final vector provides the tangent.

The diagram above shows the process applied to the input 0.95. The initial vector (green) comes from the rational approximation function; it is very close to (3, 0), but slightly rotated due to the tiny residual angle. Each rotation step generates a new vector that matches the angle from the corresponding step in the first part; I show two of the rotation matrices. Note that the vectors diverge from the circle—a consequence of dividing each matrix by cos θ—but this doesn't affect the tangent.

The FPTAN instruction is somewhat peculiar as it doesn't return the tangent directly. Instead, it returns the Partial Tangent: the two coordinate values (X and Y) that can be divided to produce the tangent. One might wonder why the chip didn't do the division automatically. The motivation was that division was slow back then, and omitting the division allowed for optimizations in some cases.

In this section, I'll explain some features of the 8087 that are important for the microcode implementation of FPTAN.

The 8087 supports a variety of data types, but internally, everything is stored as an 80-bit floating-point number called a "temporary real". A number has three parts: a sign bit, a 15-bit exponent, and a 64-bit significand (the fractional part), In most cases, a floating-point number is represented by sign × significand × 2exponent. The advantage of floating-point numbers is that the exponent allows them to cover a huge range, from very small to very large. The significand is a 64-bit binary number of the form 1.bbb…: a leading 1, followed by the binary point (the binary equivalent of the decimal point) and the rest of the bits.11 One important detail is that the exponent is stored with a "bias" of 16383 added to it. Thus, the stored exponent is always positive, even if the real exponent is negative. For instance, an exponent of 1 is stored as 0x4000, and an exponent of -16 is stored as 0x3fef.

The 80-bit temporary real format and the register formats. The dot indicates the binary point, analogous to the decimal point.

To use the 8087, a programmer stores values in its eight internal registers, organized as a stack. Each register holds an 80-bit floating-point number. To optimize performance, each value in the register stack also has an associated "tag" value: zero, valid, special, or empty. A value of zero is tagged as zero. A "normal" floating-point value is tagged as valid. If the value is infinity, Not a Number (NaN), or a denormalized value, then it is tagged as special. Finally, an empty tag indicates that a register does not hold a value; this allows detection of stack underflow (reading an empty register) or stack overflow (storing to a non-empty register).

The 8087 also has temporary registers that it uses internally: tmpA, tmpB, and tmpC. (These registers are heavily used for the tangent calculation.) tmpA and tmpB are 80-bit registers, along with two tag bits. However, tmpC is different: it doesn't have a sign, exponent, or tag. Moreover, the significand is 68 bits; tmpC has an additional bit to prevent overflow and three additional low-order bits for rounding, known as guard, round, and sticky.

In the microcode, each 16-bit micro-instruction performs one low-level operation. Many of the operations move data from one register to another. Other operations provide conditional jumps, subroutine calls, and returns. A bit shift operation takes two micro-instructions: the first moves a value to the shifter and configures the direction and amount of shift. The second micro-instruction (which may be a distance after the first) moves the result from the shifter to a destination.

Similarly, addition takes two micro-instructions. The first micro-instruction specifies the first argument to add and selects various options such as carry and rounding. (For subtraction, an option complements the second argument.) The second argument always comes from the B register, which confusingly is unrelated to the tmpB register. The second micro-instruction moves the result from the sum register to a destination.

In this section, I'll discuss some highlights of the tangent microcode. The listing below shows the 8087's microcode, along with my comments. Each line indicates a 16-bit micro-instruction. The control flow is a bit complicated, so I've put a flowchart in the footnotes.12

FPTAN: #1039 st(0) -> tmpA store argument in tmpA #1040 stackPtr-- check if room on stack to push result #1041 stack overflow? #1042 jmp #1022 if tmp empty/special/overflow/div exit if bad argument or stack overflow #1043 jmp #1045 if tmpA:tag ZERO is argument 0? #1044 except:precision precision exception except for 0 #1045 expconst 0x3ff0 Test if exponent >= -16 #1046 adder: tmpA:exp + 0 cin=1 argument exponent + 1 #1047 expConst -> Breg #1048 adder: sumreg:frac - Breg cin=1 subtract -15 #1049 jmp #1061 if adder pos CORDIC if exponent >= -16 #1050 expconst 0x3fff non-CORDIC path: #1051 tmpA:exp -> Breg check original exponent #1052 expConst -> tmpB:frac against exponent 0 #1053 adder: tmpB:frac - Breg cin=1 subtract #1054 sumreg:frac -> expConv 3fff - exp (i.e. -exp unbiased) #1055 sumreg:frac -> shiftcount #1056 jmp #1082 if exponent[6:14] != 0 if exp > -64 goto rational approximation #1057 high bit -> tmpB:frac otherwise return original argument #1058 expConst -> tmpB:sign,exp push 1 for X (denom): result is orig angle #1059 tmpB -> st(0) #1060 RNI #1061 expconst 0x0010 CORDIC path #1062 sumreg:frac -> loopcounter loopcounter = exp + 16 #1063 tmpA:frac -> tmpC tmpC = original angle #1064 trig const -> Breg/rnd #1065 adder: tmpC - Breg cin=1, roundmode subtract first CORDIC angle #1066 adder sign -> cordic, shl save bit in shift register #1067 jmp #1079 if not adder pos #1068 sumreg:frac,rnd -> tmpC #1069 adder: tmpA:exp + 0 cin=1 increment exp if first angle used #1070 sumreg:frac -> tmpA:exp #1071 jmp #1079 #1072 shift tmpC L 0 bytes, 1 bits top of CORDIC scan loop #1073 shift L -> tmpC shift angle left #1074 trig const -> Breg/rnd #1075 adder: tmpC - Breg cin=1, roundmode subtract CORDIC angle #1076 adder sign -> cordic, shl save decision bit in shift register #1077 jmp #1079 if