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

Question

我正在使用软件渲染器开发一款游戏，以获得最准确的 PS1 外观。当我研究 PS1 图形/渲染系统的工作原理、顶点摇晃的原因等时，我看到了一些关于它们如何划分的文档。这是它的链接: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”表是什么？我们为什么要改变事情？如果有人可以更详细地解释这件事实际上是如何实现鸿沟的，我们将不胜感激。

Answer 1

该算法是[0,1)中两个无符号16位小数值的定点除法。它首先通过表查找计算除数倒数的初始 9 位近似值，然后使用单个 Newton-Raphson 迭代对倒数 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 的原因尚不清楚；这可能是由于要求此步骤始终提供对真实倒数的高估。

接下来的两个步骤是 Newton-Raphson 迭代的逐字翻译，用于考虑定点比例因子的倒数。所以 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 ).