这种除法近似算法如何工作？

我正在使用软件渲染器开发游戏,以获得最准确的PS1外观.当我正在研究PS1图形/渲染系统如何工作,摇摆顶点等的原因时,我偶然发现了一些关于它们划分方式的文档.以下是它的链接:http://problemkaputt.de/psx-spx.htm#gteoverview(参见"GTE Division Inaccuracy"部分)

相关代码:

  if (H < SZ3*2) then                            ;check if overflow
    z = count_leading_zeroes(SZ3)                ;z=0..0Fh (for 16bit SZ3)
    n = (H SHL z)                                ;n=0..7FFF8000h
    d = (SZ3 SHL z)                              ;d=8000h..FFFFh
    u = unr_table[(d-7FC0h) SHR 7] + 101h        ;u=200h..101h
    d = ((2000080h - (d * u)) SHR 8)             ;d=10000h..0FF01h
    d = ((0000080h + (d * u)) SHR 8)             ;d=20000h..10000h
    n = min(1FFFFh, (((n*d) + 8000h) SHR 16))    ;n=0..1FFFFh
  else n = 1FFFFh, FLAG.Bit17=1, FLAG.Bit31=1    ;n=1FFFFh plus overflow flag

我很难理解它是如何工作的,这个'unr'表是什么？我们为什么要转移东西？如果有人可以更详细地解释这件事实际上如何实现鸿沟,那将是值得赞赏的.

该算法是[0,1）中两个无符号16位小数的定点除法。它首先通过表查找计算除数的倒数的初始9位近似值，并使用单次Newton-Raphson迭代对该倒数x _{i + 1}进行细化，x _{i + 1}：= x _i *（2- d * x _i） ，导致精确到约16位的倒数，最后将其乘以除数，得到[0,2）的17位商。

对于表查找，首先通过应用比例因子2 ^z将除数归一化为[0.5，1），然后显然需要使用相同的比例因子来调整除数。由于[0.5，1）中操作数的倒数将为[1,2]，因此倒数的整数位已知为1，因此可以使用8位表项生成1.8定点通过加0x100（= 1）来实现倒数。0x101这里使用的原因尚不清楚；可能是由于要求该步骤始终提供对真实倒数的高估。

接下来的两个步骤是考虑定点比例因子的牛顿-拉夫逊迭代的逐字转换。因此0x2000000代表2.0。该代码之所以使用，是0x2000080因为它合并了一个舍入常数0x80（= 128），用于随后的256除法，用于对结果进行重新缩放。下一步同样将其0x00000080作为舍入常数，以进行256的重新缩放比例。如果没有缩放比例，这将是纯乘法。

最终的乘法n*d将倒数d与被除数相乘n，得到33位的商。同样，在除以65536之前将舍入常数设为0x8000，以重新缩放为适当的范围，从而得到1.16定点格式的商。

连续重新定标是定点计算的典型方法，在定点计算中，尝试将中间结果保持尽可能大，以最大程度地提高最终结果的准确性。有点不寻常的是，舍入应用于所有中间算术，而不仅仅是在最后一步。也许有必要达到指定的精度水平。

但是，该功能并非十分准确，可能是由于初始近似值的不准确所致。在所有非异常情况下，2,424,807,756与正确舍入的1.16定点结果相匹配，780,692,403的误差为1 ulp，15,606,093的误差为2 ulp，86,452的误差为3 ulp。在一个快速实验中，初始近似值的最大相对误差u为3.89e-3。改进的表查找将最大相对误差降低u到2.85e-3，从而减少了但不消除最终结果中的3-ulp错误。

如果要查看特定示例，请考虑h= 0.3（0x4ccd）除以SZ3= 0.2（0x3333）。然后z= 2，因此d= 0.2 * 4 = 0.8（0xcccc）。这导致u= 1.25（0x140）。由于估算值非常准确，因此我们期望（2- d * u）接近1，实际上d= 1.000015（0x10001）。精确的倒数为d= 1.250015（0x14001），因此商为n= 1.500031（0x18002）。