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

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'表是什么？我们为什么要移位？ 如果有人能够更详细地解释一下这个东西是如何实现除法的，那就太好了。

- vexe

---

请尝试访问http://codereview.stackexchange.com/。 - OldProgrammer

它实现了牛顿-拉弗森除法（维基百科链接）。 - Nominal Animal

这个问题不适合在Code Review上讨论，即将被关闭 - Code Review不能进行“代码解释”或“为什么/如何工作”的内容。 - Der Kommissar

@vexe：unr_table 是一个[1,2)中数字的倒数表。请注意常量101h的添加：与表项配合使用，这形成了一个9位固定点近似值，其中 0x100 相当于 0.5。我不确定为什么他们添加了 101h 而不是我所期望的值 100h；额外的加1可能代表对以下固定点数学的截断性质进行补偿的调整因素。 - njuffa

Answer 1

该算法是用于计算两个无符号16位小数值[0,1)的固定点除法。它首先通过查表计算出除数倒数的初始9位近似值，然后使用牛顿-拉夫逊迭代一次来精确计算倒数，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位表项通过添加0x100（=1）来生成一个1.8固定点倒数。这里使用 0x101 的原因不明确，可能是由于该步骤总是提供真实倒数的过估计的要求。

接下来的两个步骤是牛顿-拉夫逊迭代计算倒数的直译，考虑到固定点比例尺，故使用 0x2000080 代表2.0，并且该代码使用 0x00000080 作为舍入常数，用于后续的除以256进行重新缩放。最后的乘法n*d将d中的倒数与n中的被除数相乘，得到33位商数。再次应用0x8000的舍入常量，然后除以65536进行重新缩放，从而得到1.16固定点格式的商数。

连续缩放在固定点计算中是典型的，在这种情况下，人们尝试尽可能保持中间结果尽可能大，以最大化最终结果的精度。有点不寻常的是，所有中间算术都应用了舍入，而不仅仅是在最后一步应用。也许这是为了实现指定级别的准确性。

该函数并不是完全准确的，这可能是由于初始估计值的不准确导致的。在所有非异常情况下，有2,424,807,756个匹配正确四舍五入的1.16固定点结果，780,692,403个误差为1 ulp，15,606,093个误差为2-ulp，86,452个误差为3-ulp。在一个快速实验中，初始估计值的最大相对误差为3.89e-3。改进后的表查找将的最大相对误差降至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)。

- njuffa