这个除法近似算法是如何工作的？

Question

我正在使用软件渲染器开发一款游戏以获得最准确的 PS1 外观。当我研究 PS1 graphics/rendering 系统如何工作、顶点摇晃的原因等时，我遇到了一些关于他们划分方式的文档。这是它的 link：http://problemkaputt.de/psx-spx.htm#gteoverview（参见“GTE 分部误差”部分）

相关代码：

  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' table？我们为什么要改变事情？如果有人能更详细地解释这件事实际上是如何实现鸿沟的，我们将不胜感激。

Answer 1

此算法是 [0,1) 中两个无符号 16 位小数值的 fixed-point 除法。它首先通过 table 查找计算除数倒数的初始 9 位近似值，然后使用单个 Newton-Raphson 迭代对倒数 x_{i+1[= 进行细化55=] := x_i * (2 - d * x_i)，得到的倒数精确到16位左右，最后乘以这个通过被除数，在 [0,2].}

中产生一个 17 位商

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

接下来的两个步骤是考虑 fixed-point 比例因子的倒数 Newton-Raphson 迭代的逐字翻译。所以0x2000000代表2.0。该代码使用 0x2000080，因为它包含一个舍入常数 0x80 (=128)，用于随后除以 256，用于重新缩放结果。下一步同样添加 0x00000080 作为重新缩放除以 256 的舍入常数。如果没有缩放，这将是一个纯乘法。

最后的乘法 n*d 将 d 中的倒数与 n 中的被除数相乘，得到 33 位的商。同样，在除以 65536 之前应用舍入常数 0x8000 以重新调整到适当的范围，以 1.16 fixed-point 格式给出商。

连续缩放是 fixed-point 计算的典型特征，其中人们试图使中间结果尽可能大，以最大限度地提高最终结果的准确性。有点不寻常的是舍入应用于所有中间算术，而不仅仅是最后一步。也许有必要达到特定的精度水平。

不过，该函数并不是那么准确，这可能是由于初始近似值不准确造成的。在所有 non-exceptional 个案例中，2,424,807,756 个匹配正确四舍五入的 1.16 fixed-point 结果，780,692,403 个有 1 ulp 的误差，15,606,093 个有 2-ulp 的误差，86,452 个有 3-ulp 的误差。在快速实验中，初始近似 u 中的最大相对误差为 3.89e-3。改进的 table 查找将 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)。

这个除法近似算法是如何工作的？

How does this division approximation algorithm work?

c

bit-shift

division

approximation