逆向工程经典Intel 8087的正切算法：不止于CORDIC

# 逆向工程复古Intel 8087的正切算法：超越CORDIC 来源：http://www.righto.com/2026/09/8087-tangent-cordic.html 希望你还没看腻8087，因为我又为这款Intel浮点芯片准备了一篇文章。1980年，Intel推出了8087，使IBM PC和其他系统的浮点运算速度大幅提升。本文将解析该芯片正切指令背后的算法。三角函数的一种常用实现方法是CORDIC算法，另一种则是多项式逼近。8087融合了这两种方法，同时实现了高精度与高性能。相较于8086微处理器，8087提供了巨大的速度提升：计算正切仅需90微秒，而非13000微秒。 通过研究8087的电路和微码，我可以解释正切指令（即`FPTAN`）背后的算法。为探索8087的电路结构，我用凿子打开了芯片外壳，并通过显微镜获取了高分辨率图像。微码ROM位于芯片中央的矩形区域，存储着控制芯片运行的1648条微指令。芯片下半部分（红色方框内）是数据通路，负责处理80位值的浮点运算。 ### 8087数据通路近景  *8087数据通路近景，展示FPTAN使用的功能模块。点击图片（或任意图片）查看更大版本。* 放大观察数据通路可见相关功能单元：指数ROM存储算法所需的固定指数值；常量ROM存放包括CORDIC算法所用常量在内的各种参数；移位器是大型组件，可对64位值进行任意位数的左右移位；加法器是8087计算的核心，除加减运算外，还通过循环实现乘法、除法和开平方；B寄存器保存加法器的一个输入，而多个来源可提供另一个输入；求和寄存器保存加法器输出；八个栈寄存器和临时寄存器用于存储浮点数；移位寄存器则保存CORDIC计算所需的16位状态位。 ### CORDIC算法原理 CORDIC是一种巧妙的算法，能用简单硬件快速计算超越函数：它通过移位、加法和查表操作，无需乘除运算即可完成计算。该算法可追溯至1956年，为B-58“盗贼”轰炸机（首款能以2马赫飞行的轰炸机）开发。该机原配备模拟导航计算机，但模拟器件精度有限。工程师杰克·沃尔德受命设计数字计算机以替代模拟系统，其关键挑战在于：模拟计算机可通过名为旋变器的机电设备轻松生成正负弦函数，但数字方式（尤其在当年晶体管速度有限的情况下）生成三角函数极为困难。  *康维尔B-58A盗贼轰炸机，展于德克萨斯州圣安东尼奥市* 杰克·沃尔德提出用简单硬件快速计算三角函数的方法。他将该算法（及实现它的计算机）命名为CORDIC：“坐标旋转数字计算机”。CORDIC将角度转换为向量，通过向量坐标获取所需三角函数。其诀窍在于将角度分解为一系列特殊角度，这些角度使向量旋转易于操作。特殊角度经预先计算并存储于表中，使得即使在1950年代的硬件上也能快速执行CORDIC计算。每次CORDIC迭代增加一位精度，因此算法收敛迅速。CORDIC广受欢迎，曾应用于科学计算器（采用十进制CORDIC而非二进制）。 ## 三角函数基础 本节尽量简化数学表达，快速解释CORDIC工作原理。下图回顾三角函数与点坐标的关系：假设角度θ确定单位圆上一点*(X, Y)*，基本公式为*X=cos θ*、*Y=sin θ*、*Y/X = tan θ*。因此只要确定坐标*(X, Y)*，就能求出三角函数值。若点不在单位圆上（例如*(X', Y')*），仍可轻松计算tan θ（提示：这正是8087的实现方式），但sin θ和cos θ的计算会变得复杂。  *角度、X/Y坐标与三角函数的关系* 熟悉计算机图形学的人可能了解旋转矩阵如何旋转角度点。将点*(X, Y)*乘以旋转矩阵即可得到新点*(X', Y')*，如下图所示：  *使用旋转向量旋转点* 然而旋转矩阵需要*sin*和*cos*运算（见下式1），似乎无法解决我们的问题。但若将矩阵除以*cos θ*（此时矩阵需用* tan*，即我们要计算的量），并允许向量长度增长，可得到矩阵(2)。关键在于使用特殊角度*αn= arctan\(2\-n\)*，代入矩阵后得到矩阵(3)，该矩阵在硬件中易于计算：2的幂次乘法可通过移位实现。  *旋转矩阵简化过程* 将矩阵(3)应用于点*(X,Y)*得到方程(4)，这是CORDIC处理的核心公式。其重要性在于：仅通过加减法和二进制移位即可用机器语言或硬件快速计算。 由此可知CORDIC工作原理：首先将输入角度分解为一系列特殊角度的组合。从单位向量(1,0)开始，对每个特殊角度应用旋转公式，最终得到目标角度的点*(X, Y)*，所求正切值即为*Y/X*。由于特殊角度表预先计算，*arctan*操作不会拖慢进程。补充说明：从第三项起，特殊角度近似于*2\-n*，每步约缩小一半。 ## 有理多项式逼近 CORDIC精度取决于使用项数：16项提供约2\-16（16位）精度，要达到64位精度需计算64项（及64个特殊角度表）。为提高速度，8087采用16位CORDIC配合另一算法处理残余角度（残余角度是特殊角度累加和与目标角度的差值，约2\-16量级）。 8087对残余角度使用Padé近似（两个多项式的比值）。Padé近似有完整族系，取决于多项式阶数。8087采用简单公式：*3x/\(3\-x2\)*。该近似虽简单，但对小角度极其精确，误差与*x4*成正比。由于*x<2\-16*，误差将小于*2\-64*，满足8087的64位精度要求。此外，8087无需执行有理多项式的除法，因为`FPTAN`返回分子和分母两个值，除法实质是“免费”的。  *正切函数（红）、有理逼近（蓝）、泰勒级数（绿）。注：因相关范围接近0，泰勒级数表现并非图中所示那么差。图形由Desmos生成。* 若学过微积分，可能认为泰勒多项式是可行方案，但有理逼近效果更佳（原因：正切在π/2处趋向无穷，多项式无法发散，而多项式之比可以，因此更匹配正切特性）。 ## 综合实现：8087算法 8087正切算法包含三部分：确定CORDIC决策位（称为伪除法）、计算有理逼近、根据决策位应用旋转方程（称为伪乘法）。 详细而言：第一步确定哪些特殊角度用于逼近输入角度。每个特殊角度与剩余输入角度比较，若较小则减去。若减去则记录1，否则记录0。此过程类似长除法中连续减去除数（移位后版本）产生商位，故称伪除法。注意此步骤不执行旋转操作，而是决定后续要应用的旋转。  *CORDIC算法第一阶段（伪除法）：通过特殊角度缩减输入角度，最终留下残余角度（“rad”指弧度，非辐射单位）* 接下来用有理逼近函数*3x/\(3\-x2\)*计算残余角度的正切值。该结果作为下一步的初始向量。此处不执行除法：分子作为向量*Y*，分母作为*X*，从而避免了昂贵的除法运算（因答案隐含了除法）。乘法3操作简便（左移加），因此本步唯一耗时操作是角度平方（需完整64位乘法）。 最后一步根据第一步的决策位（若为1）应用相应CORDIC旋转。因该步骤类似二进制乘法（若乘数位为1则累加被乘数），故称伪乘法。如前所述，每次旋转通过移位、加减实现，成本低廉。旋转顺序与第一步角度缩减顺序相反（从最小角度开始），以减小舍入误差。为此，决策位存储在16位移位寄存器中，并在本步反向移出。  *CORDIC算法最后阶段（伪乘法）：通过移位加法多次旋转向量，最终向量提供正切值* 上图展示了输入0.95的处理过程：初始向量（绿）来自有理逼近函数，接近(3,0)但因微小残余角度略有旋转。每次旋转生成新向量，匹配第一部分对应步骤的角度（图中显示了两个旋转矩阵）。注意向量偏离圆周——这是矩阵除以*cos θ*的结果——但不影响正切值。 `FPTAN`指令较为特殊：不直接返回正切值，而是返回“部分正切”：可相除得到正切的两个坐标值（`X`和`Y`）。当时除法运算较慢，省略除法可在某些场景优化性能。 ## 8087相关硬件特性 本节解释8087对`FPTAN`微码实现至关重要的特性。 8087支持多种数据类型，但内部统一存储为80位浮点数（“临时实数”）。数值包含三部分：符号位、15位指数和64位有效数字（小数部分）。通常浮点数表示为*符号*×*有效数字*× 2^*指数*。浮点数的优势在于指数允许覆盖极大范围（从极微小到极大）。有效数字是64位二进制数，格式为`1\.bbb...`：前导1、二进制小数点（十进制小数点的二进制等价物）及后续位。重要细节：指数存储采用偏置表示...