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
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.
The CORDIC algorithm
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.
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
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.
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.
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".
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).
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.
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.
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.
Relevant hardware details of the 8087
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.
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.
Implementation of the tangent algorithm in microcode
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 not adder pos #1078 sumreg:frac,rnd -> tmpC #1079 jmp #1072 if not const latch zero bottom of CORDIC scan loop #1080 tmpC -> tmpA:frac update tmpA, tmpB #1081 expConst -> shiftcount 16 -> shiftcount (otherwise based on exponent) #1082 tmpA:frac -> tmpB:frac Padé approximation #1083 tmpA:frac -> Breg tmpA = ang #1084 call SQUARE square the angle #1085 shift sumreg:frac R count byte bit #1086 shift R -> sumreg:frac shift twice to scale square #1087 shift sumreg:frac R count byte bit #1088 shift R -> Breg Breg = ang^2 #1089 1.1 -> tmpB:frac 3 (with exp 1) #1090 adder: tmpB:frac - Breg cin=0 #1091 sumreg:frac -> tmpB:frac tmpB (X) = 3-ang^2 #1092 shift tmpA:frac R 0 bytes, 1 bits #1093 shift R -> Breg/rnd #1094 adder: tmpA:frac + Breg cin=0, roundmode ang + ang/2 #1095 sumreg:sign,frac,rnd -> tmpC tmpC (Y) = 3×ang with appropriate exponent #1096 expconst 0x3ff0 exponent -15 #1097 expConst -> Breg #1098 adder: tmpA:exp - Breg cin=1 #1099 jmp #1118 if not adder pos skip if exponent too small #1100 shift tmpC R 0 bytes, 1 bits CORDIC pseudo-multiplication path #1101 #15 -> loopcounter #1102 shift R -> tmpC #1103 jmp #1113 if not cordic[0] CORDIC: if saved decision bit... #1104 shift tmpC R loopcount send tmpC (Y) to shifter #1105 tmpC -> Breg/rnd send tmpC to adder #1106 adder: tmpB:frac + Breg cin=0, roundmode tmpB (X) not shifted because tmpC scaled #1107 sumreg:sign,frac,rnd -> tmpC tmpC = tmpC + tmpB (implicit >>n) #1108 shift R -> sumreg:frac,rnd shifted old tmpC to sumreg #1109 shift sumreg:frac,rnd R loopcount shift right again because tmpC scaled #1110 shift R -> Breg/rnd and send to B reg #1111 adder: tmpB:frac - Breg cin=1, roundmode #1112 sumreg:frac,rnd -> tmpB:frac tmpB = tmpB - tmpC>>n #1113 cordic shr shift out shift register bit #1114 shift tmpC R 0 bytes, 1 bits start shift of tmpC (Y) #1115 jmp #1118 if cordic==0 done if shift register empty #1116 shift R -> tmpC shift tmpC (Y) right #1117 jmp #1103 if not const latch zero done if counter at zero #1118 expconst 0x3fff CORDIC done #1119 tmpB:frac -> sumreg:frac Test top bit of tmpB (X) #1120 jmp #1124 if reg bit 63 If set, use exponent 0 #1121 expconst 0x3ffe Otherwise, use exponent -1 and #1122 shift sumreg:frac L 0 bytes, 1 bits shift left one bit to normalize #1123 shift L -> tmpB:frac #1124 expConst -> tmpB:sign,exp store exponent in tmpB #1125 tmpC -> tmpA:frac store tmpC (Y) in tmpA #1126 tmpC -> sumreg:frac,sign #1127 jmp #1132 if not sumreg[64] test overflow bit #1128 shift tmpC R 0 bytes, 1 bits if set, normalize by shifting right #1129 shift R -> tmpA:frac #1130 adder: tmpA:exp + 0 cin=1 and increment exponent by 1 #1131 sumreg:frac -> tmpA:exp #1132 tmpB -> st(0) tmpB to top of stack: X (denom) #1133 stackPtr++ #1134 tmpA -> st(0) tmpA to st(1): Y (numerator) #1135 stackPtr-- #1136 RNI End of FPTAN
The FPTAN microcode starts at address #1039 (decimal).
The microcode starts by moving the argument from the top of the stack to the tmpA register.
If the top value is empty, this indicates a stack underflow. If the next element on
the stack (the position that will hold the result) is not empty, this indicates a stack
overflow. In either case, the microcode indicates an "invalid" exception and ends the instruction.
One peculiar feature of the 8087 is that it flags a result with a "precision" exception if
the result is not exact. (This happens very frequently, even for, say, 1/10.)
The only tangent that the 8087 can compute precisely is tan(0), so all other inputs
result in a precision exception (#1044).
Next, the argument's exponent is tested, splitting the execution into three paths.
If the exponent is -64 or less, the argument is so small that
the tangent equals the argument, within the accuracy of the system.13
In this case (#1057), the original argument is returned unchanged (with the denominator 1 pushed), and the code ends.
The second case is if the argument's exponent is -17 or less. In this case,
the CORDIC step is skipped, and the routine jumps directly to the rational approximation (#1082).
CORDIC pseudo-division
The most interesting case is the CORDIC path (#1061), taken if the exponent is between -1 and -16.
Conceptually, the code performs 16 CORDIC steps, using 16 stored angles.
However, the code is optimized, skipping the large angles for smaller
inputs.
Specifically, the loop counter is initialized based on the
exponent of the input angle (#1062).
The CORDIC pseudo-division loop (#1061-#1080) tests the angles, subtracting ones that aren't too big.14
The decision results are recorded in a 16-bit shift register, located on the far right side of the die. Presumably, this is where the layout had some unused space. Since the shift register is loaded and unloaded serially, it doesn't need access to the internal fraction bus, just two lines to shift the bits in and out, so the shift register could be located in otherwise unused space.
The code is optimized to use 64-bit integer arithmetic rather than floating-point arithmetic. This makes it tricky to understand the code since values must be viewed more as fixed-point numbers with implicit exponents. These exponents are not stored anywhere, but can be determined by analyzing the algorithm. Even the angle constants in the ROM do not have explicit exponents.15 Moreover, the implicit exponents change for each step through the loop to preserve accuracy. Roughly speaking, each cycle of the loop reduces the values by a factor of 2, and the angle constants shrink accordingly. If the values had a fixed exponent, they would end up with 16 leading zeros, wasting precision. But by scaling the values each cycle (with a left shift), the numbers continue to use the full 64 bits. In other words, the hardware is performing fast integer arithmetic, but mathematically you can think of it as fixed-point with exponents that don't physically exist in the chip. See the table in the footnotes16 for details on how the implicit exponents change through the loop.
Rational polynomial approximation
After the CORDIC pseudo-division loop, the microcode calculates the rational polynomial
approximation (#1082).
The small-angle path rejoins the execution flow here.
The microcode to calculate the polynomial approximation has a few interesting features.
The most time-consuming part is the SQUARE routine, which multiplies a fixed-point number by itself.
Multiplication is complicated, so I'll give the details in a later post.
But in brief, the 8087 uses Booth's Algorithm to
multiply by two bits at a time (radix 4), so it is twice as fast as regular binary
multiplication.
The loop to perform the shifts and adds is implemented in hardware, rather than microcode.
That is, one microcode instruction performs 32 additions: the hardware tests the
bits, does the appropriate add or subtract, updates the loop counter, and loops.
For performance, the shifting is done with specialized shifters, not the 8087's general-purpose
shifter.
Except for the square, the polynomials are calculated with additions rather than multiplications. The constant 3 doesn't come from the constant ROM, but from special-purpose transistors that set the top two bits of the significand; this is usually used to supply an NaN. (If the implicit exponent is 1, this corresponds to 3.) Subtracting the square yields the numerator. For the denominator, the angle is shifted right (i.e. divided by 2) and added to itself. This performs a multiplication by 3 (technically by 1.5, but increasing the implicit exponent to 1 turns this into 3). It would be expensive to divide the numerator by the denominator; instead, the numerator and denominator are used as the initial Y and X values for the following CORDIC steps.
CORDIC pseudo-multiplication
Next is the CORDIC loop where the rotations are applied to the vector (#1100).
The rotations are applied in the reverse order from how they were computed in the first step:
the smallest rotations are applied first to preserve accuracy.
Bits are shifted out of the shift register to indicate if the rotation should be applied or
not. As soon as the shift register is all zeros, the loop stops.
In other words, smaller angles go through the loop just a few times, not the full 16 times.
The rotations are applied by shifts and adds (or subtracts), but keeping track of the
implicit exponents is a bit tricky.
Since the vector starts almost horizontal, the X value is approximately 3, and the Y value
is around 2-16.
As the vector rotates, X remains roughly 3, but Y potentially increases by a factor of 2
each time.
To provide the maximum accuracy for each number,
the exponent for X (tmpB) is 1, while the exponent for Y (tmpC) starts off at
-14 and increases by 1 each loop as Y is shifted right.
(Remember that these exponents are not stored anywhere but are implicit.)
Since the two registers have different (implicit) exponents, the values must be shifted
before adding.
The shift amounts aren't intuitive; the table16 in the
footnotes may help clarify.
After the CORDIC loop,
the final part of the code (#1118) normalizes the X and Y values, creating the floating-point
values that are returned from the instruction.
Recall that these aren't "real" floating-point values at this point, so they
may need adjustment.
Specifically, if X doesn't have a
leading 1, it is shifted left one position, while if Y has two leading digits, it
is shifted right one position. In either case, the exponent is adjusted accordingly.
The X and Y values are put on the stack (#1132) to complete the instruction.
Conclusions
The FPTAN instruction is fairly slow as 8087 instructions go, due to its complexity.
The documentation says that it takes typically 450 clock cycles.
(The time has a large range; depending on the value, it can take 30 to 540 clock cycles.)
For the value I examined (0.95), 33% of the time is in the CORDIC
pseudo-division, 15% in the rational polynomial (mostly the squaring operation), 47% in the CORDIC pseudo-multiplication, and 5% overhead.
CORDIC is a popular way to compute trig functions, but there are alternatives such as polynomials. The 8087 is unusual because it combines CORDIC and a rational polynomial. For the Pentium, Intel moved from CORDIC to polynomial approximations; the Pentium's fast multiplication circuitry made polynomials practical. Nowadays, Intel has libraries such as MKL (Math Kernel Library) and Short Vector Math Library (SVML) that provide highly optimized implementations tailored to Intel hardware. These libraries are said to use polynomial approximations, using parallel instructions (SIMD) for performance, rather than specialized hardware. With these libraries, x87 operations and 80-bit floats are mostly obsolete. Nonetheless, I hope you've enjoyed this look at an original 8087 algorithm.
I plan to continue reverse-engineering the 8087 microcode; for updates, follow me on Bluesky (@righto.com), Mastodon (@[email protected]), or RSS. I've been working on this with the members of the "Opcode Collective", especially Smartest Blob and Gloriouscow, who converted the ROM images to microcode data and extensively analyzed the contents. See the 8087 repository on GitHub for more. AI statement: Despite the presence of the em dash, no AI was used in the writing of this article (details). The Opcode Collective used an ML algorithm to classify each ROM cell to extract the microcode.
Notes and references
-
If you want more on the 8087, here's a list of my previous 8087 posts, newest to oldest:
Microcode in Intel's 8087 floating-point chip: the scale instruction
The adder at the heart of Intel's 8087 floating-point chip
Microcode inside the Intel 8087 floating-point chip: register exchange
Instruction decoding in the Intel 8087 floating-point chip
Conditions in the Intel 8087 floating-point chip's microcode
The stack circuitry of the Intel 8087 floating point chip, reverse-engineered
Die analysis of the 8087 math coprocessor's fast bit shifter
Extracting ROM constants from the 8087 math coprocessor's die
Two bits per transistor: high-density ROM in Intel's 8087 floating point chip
Inside the die of Intel's 8087 coprocessor chip, root of modern floating point
↩ -
Intel's documentation provides a performance comparison of the 8087 versus floating-point operations emulated on the 8086. The speedup varies depending on the instruction, but is typically almost a factor of 100. This speedup made the 8087 chip highly desirable for numerical applications such as CAD and spreadsheets. ↩
-
Internally, the 8087 uses 80 bits to hold floating-point values: one bit for the sign, 15 bits for the exponent, and 64 bits for the significand. However, the 8087 uses three extra low-order bits for rounding, called Guard, Round, and Sticky. These bits ensure that a value is always rounded in the right direction. Some parts of the datapath have additional bits for sign or overflow: the shifter is 68 bits wide, and the adder is 69 bits wide. For the most part, I'll ignore these extra bits and refer to the datapath as 64 bits wide. ↩
-
To compute sin θ from the 8087's X and Y values, the typical formula is sin θ = Y/sqrt(X2+Y2); cos θ is similar. ↩
-
The special angles (in degrees) are approximately 45°, 26.6°, 14°, 7.1°, 3.6°, 1.8°, 0.9°, … You might wonder if adding up the special angles will always converge to the desired angle θ. A proof can be found here. Note that arctan(2-n) rapidly approaches 2-n (in radians), so it's kind of like a binary fraction expression, cutting the gap approximately in half each time. ↩
-
The CORDIC implementation in the 8087 is somewhat different from the typical use of CORDIC, which finds sin and cos. In this approach, all the special angles are used, but an angle is either added or subtracted at each step, resulting in convergence to the desired angle. Recall that since the matrix is scaled by 1/sin αi at each step, the ending vector is no longer a unit vector. As a result, cos and sin don't match the final vector (X, Y).
The trick is that if all the special angles are used, the overall scaling factor ends up being the same, simply the product of all the cos αi factors. (It doesn't matter if the rotation is positive or negative because cos αi = cos -αi.) Thus, the overall scaling factor can be precomputed and the initial unit vector can be prescaled by this amount. The result is that the final vector (X, Y) yields (cos θ, sin θ) directly, which is what we want.
The key point is that there are two CORDIC approaches: using positive and negative angles, or using or omitting angles. (The first case multiplies each angle by +1 or -1, while the second case multiplies each angle by 1 or 0.) Most CORDIC implementations use the first approach, while the 8087 uses the second approach. ↩
-
Jack Volder, the creator of CORDIC, wrote an article, The Birth of CORDIC, that describes the history in detail, as well as including photos of the B-58 navigation computer. Volder's original articles from 1959 are The CORDIC Computing Technique and The CORDIC Trigonometric Computing Technique.
The CORDIC-II digital computer system. Photo by General Dynamics, from "The Birth of CORDIC".The photo above shows the CORDIC II computer that was used on the B-58 aircraft. It was a 30-bit serial computer, running at 198.4 kHz and weighing a hefty 134 pounds. It was implemented with diode-transistor logic (DTL) gates, using 2498 transistors and 4265 diodes, and had a drum memory (visible at the top) for storage, holding the equivalent of 39 kilobytes. Interestingly, the CORDIC operations were implemented as basic instructions in the instruction set, rather than subroutines; a CORDIC instruction took 5 ms, the same speed as a multiplication or division. ↩
-
This Padé approximation is called the [1,2] approximation since the numerator is a first-order polynomial and the denominator is a second-order polynomial. Padé approximations can be derived from the continued fraction expansion of tan. See Constructing Padé Approximates.
By the way, if you've taken calculus, you might think that a Taylor series is a good way to approximate a function. Unfortunately, a Taylor series is bad at approximating a function over an interval, rather than a point. You're much better off using a polynomial that minimizes the maximum error over the interval, as in the Pentium. ↩
-
The valid argument range for
FPTANis documented as 0<θ<π/4, so the demonstration value that I use, 0.95, is outside the documented range. I find the range puzzling, since the first CORDIC constant is arctan(1), π/4, which will only be useful for inputs above π/4. Thus, the 8087 uses a scarce ROM entry and an extra loop cycle to handle values outside the permitted range. I have two hypotheses about this. Perhaps the 8087 was implemented for a larger argument range, but the numerical error was larger than expected, so the range was truncated. Alternatively, theFPATANarctan instruction might need the first CORDIC constant and it was easier to use the same constant loop forFPTANandFPATAN. (I haven't examinedFPATANyet to confirm this.) In any case, the code works with 0.95; I use it as an example because the angles are larger and more visible.Another peculiarity is that the documentation is inconsistent about whether 0 is in the range or not. The book The 8087 Primer by the principle architects of the 8087 and 8086 states that the argument range excludes 0, a decision that was made to simplify the microcode. However, looking at the microcode (
#1044) shows that the code explicitly checks for a zero input and handles the case.FPTANappears to work for values almost up to π/2, with some accuracy loss as the result approaches infinity near π/2. However, the values from 1 to π/2 (which are outside the defined range and have an exponent of 0) probably work by coincidence, not by design. The larger exponent results in an additional loop iteration; the lookup of trig constants wraps around the algorithm starts with the last constant, atan(2-15). But due to scaling, this value ends up being interpreted as approximately 2. Since this CORDIC angle is greater than π/2, it is skipped. The algorithm then proceeds with the remaining constants, executing properly.Footnote to the footnote: Although the mathematical meaning of "range" is the output of a function, Intel uses it to describe the input of a function, which mathematicians call the domain. ↩
-
One difference between the 8087's temporary real format and the external formats is that in the temporary real, the leading 1 digit is explicitly stored in the number. In the external formats, the leading 1 is implicit and not stored in the number, providing one additional bit of accuracy, but making arithmetic more complicated. ↩
-
The flowchart below shows the control flow for the
FPTANmicrocode. The main blocks are the CORDIC pseudo-division, the rational approximation, and the CORDIC pseudo-multiplication. (Click it for a large image.)
Flowchart of the FPTAN microcode. -
The comparisons take multiple micro-instructions, for example, selecting a constant such as
0x3ff0(recall that exponents are biased by 16383 or0x3fff, so0x3ff0corresponds to -15), moving the constant to the adder, subtracting the constant from the exponent, and branching on the result. ↩ -
The first step of the CORDIC algorithm (
#1061-#1068) is handled separately. It took me a while to figure out why the first step of the CORDIC algorithm is special. The cause is that the CORDIC angles are close to 2-n, but not exactly. Thus, the number of CORDIC cycles does not change exactly when the exponent of the input angle changes. This is a problem because the number of pseudo-division cycles depends on the exponent, but the number of pseudo-multiplication cycles depends on the largest CORDIC angle used. Since the numbers are scaled in both loops, a mismatch in the number of cycles leaves the result off by a factor of 2. The first CORDIC step tests if the first CORDIC constant is used and increments thetmpAexponent by 1 to counteract this. (Don't worry if this doesn't make sense.) ↩ -
The constant ROM in the 8087 holds the significands, but not the exponents. When I originally extracted the constants from the constant ROM, back in 2020 (link), I couldn't find the exponents, so I came up with plausible exponents that resulted in meaningful values. The 8087 has a separate exponent ROM. Now that I've examined the microcode, I can see that the exponents in the exponent ROM are unrelated to the significands in the constant ROM. I'm glad I didn't waste much time trying to figure out the exponents back then :-) ↩
-
The following table shows the calculations at each step of the algorithm, using 0.95 as the input. Note that the calculations are done with integer arithmetic; the exponents (gray) are not stored as part of the numbers, but are implied by the algorithm. The exponents change in each step, so as the values get smaller during the pseudo-division stage, the exponents get smaller too, ensuring that the values can use most of the available bits. The opposite happens during pseudo-multiplication: the values start off small and grow as the vector is rotated; the exponent is increased each step to counteract this. Unlike "real" floating-point numbers, the bit representations are not normalized; sometimes they have two 1 digits to the left of the binary point, and sometimes they have leading zeros. That is, the exponents are fixed for each step; they don't change based on the value.
↩↩operation bits value flag reg load with θ 1.1110011…×2-1 0.9500 tmpC cmp atan(20) 1.1001001…×2-1 0.7854 1 const sub angle 0.0101010…×2-1 0.1646 tmpC shift left 0.1010100…×2-2 0.1646 tmpC cmp atan(2-1) 1.1101101…×2-2 0.4636 0 const shift left 1.0101000…×2-3 0.1646 tmpC cmp atan(2-2) 1.1111010…×2-3 0.2450 0 const shift left 10.1010001…×2-4 0.1646 tmpC cmp atan(2-3) 1.1111110…×2-4 0.1244 1 const sub angle 0.1010010…×2-4 0.04025 tmpC shift left 1.0100100…×2-5 0.04025 tmpC cmp atan(2-4) 1.1111111…×2-5 0.06242 0 const shift left 10.1001001…×2-6 0.04025 tmpC cmp atan(2-5) 1.1111111…×2-6 0.03124 1 const sub angle 0.1001001…×2-6 0.009007 tmpC shift left 1.0010011…×2-7 0.009007 tmpC cmp atan(2-6) 1.1111111…×2-7 0.01562 0 const shift left 10.0100111…×2-8 0.009007 tmpC cmp atan(2-7) 1.1111111…×2-8 0.007812 1 const sub angle 0.0100111…×2-8 0.001195 tmpC shift left 0.1001110…×2-9 0.001195 tmpC cmp atan(2-8) 1.1111111…×2-9 0.003906 0 const shift left 1.0011100…×2-10 0.001195 tmpC cmp atan(2-9) 1.1111111…×2-10 0.001953 0 const shift left 10.0111001…×2-11 0.001195 tmpC cmp atan(2-10) 1.1111111…×2-11 0.0009766 1 const sub angle 0.0111001…×2-11 0.0002181 tmpC shift left 0.1110010…×2-12 0.0002181 tmpC cmp atan(2-11) 1.1111111…×2-12 0.0004883 0 const shift left 1.1100100…×2-13 0.0002181 tmpC cmp atan(2-12) 1.1111111…×2-13 0.0002441 0 const shift left 11.1001001…×2-14 0.0002181 tmpC cmp atan(2-13) 1.1111111…×2-14 0.0001221 1 const sub angle 1.1001001…×2-14 0.00009604 tmpC shift left 11.0010010…×2-15 0.00009604 tmpC cmp atan(2-14) 1.1111111…×2-15 0.00006104 1 const sub angle 1.0010010…×2-15 0.00003500 tmpC shift left 10.0100101…×2-16 0.00003500 tmpC cmp atan(2-15) 1.1111111…×2-16 0.00003052 1 const sub angle 0.0100101…×2-16 0.000004482 tmpC X approx 1.0111111…×21 3.0000 tmpB Y approx 0.0011100…×2-14 0.00001345 tmpC Y += 2-15X 1.1011100…×2-14 0.0001050 1 tmpC X -= 2-15Yold 1.0111111…×21 3.0000 tmpB shift right 0.1101110…×2-13 0.0001050 tmpC Y += 2-14X 10.0101110…×2-13 0.0002881 1 tmpC X -= 2-14Yold 1.0111111…×21 3.0000 tmpB shift right 1.0010111…×2-12 0.0002881 tmpC Y += 2-13X 10.1010111…×2-12 0.0006543 1 tmpC X -= 2-13Yold 1.0111111…×21 3.0000 tmpB shift right 1.0101011…×2-11 0.0006543 tmpC shift right 0.1010101…×2-10 0.0006543 0 tmpC shift right 0.0101010…×2-9 0.0006543 0 tmpC Y += 2-10X 1.1101010…×2-9 0.003584 1 tmpC X -= 2-10Yold 1.0111111…×21 3.0000 tmpB shift right 0.1110101…×2-8 0.003584 tmpC shift right 0.0111010…×2-7 0.003584 0 tmpC shift right 0.0011101…×2-6 0.003584 0 tmpC Y += 2-7X 1.1011101…×2-6 0.02702 1 tmpC X -= 2-7Yold 1.0111111…×21 3.0000 tmpB shift right 0.1101110…×2-5 0.02702 tmpC shift right 0.0110111…×2-4 0.02702 0 tmpC Y += 2-5X 1.1110111…×2-4 0.1208 1 tmpC X -= 2-5Yold 1.0111111…×21 2.9991 tmpB shift right 0.1111011…×2-3 0.1208 tmpC shift right 0.0111101…×2-2 0.1208 0 tmpC Y += 2-3X 1.1111101…×2-2 0.4957 1 tmpC X -= 2-3Yold 1.0111110…×21 2.9840 tmpB shift right 0.1111110…×2-1 0.4957 tmpC shift right 0.0111111…×20 0.4957 0 tmpC shift right 0.0011111…×21 0.4957 0 tmpC Y += 20X 1.1011110…×21 3.4797 1 tmpC X -= 20Yold 1.0011111…×21 2.4884 tmpB final X 1.0011111…×21 2.4884 tmpB final Y 1.1011110…×21 3.4797 tmpA
