How does this division approximation algorithm work?
I'm working on a game with a software renderer to get the most accurate PS1 look. As I was doing research on how the PS1 graphics/rendering system worked, reason for the wobbly vertices etc, I came across some documentation regarding the way they did their divide. Here is the link to it: http://problemkaputt.de/psx-spx.htm#gteoverview (see the " GTE Division Inaccuracy" section)
The relevant code:
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
I'm having a hard time understanding how this works, what is this 'unr' table? why are we shifting things? If someone could give a more detailed explanation as to how this thing is actually achieving the divide, it would be appreciated.
This algorithm is a fixed-point division of two unsigned 16-bit fractional values in [0,1). It first computes an initial 9-bit approximation to the reciprocal of the divisor via a table lookup, refines this with a single Newton-Raphson iteration for the reciprocal, xi+1 := xi * (2 - d * xi), resulting in a reciprocal accurate to about 16 bits, finally multiplies this by the dividend, yielding a 17-bit quotient in [0,2).
For the table lookup, the divisor is first normalized into [0.5, 1) by applying a scale factor 2z, obviously the dividend then needs to be adjusted by the same scale factor. Since the reciprocals of operands in [0.5, 1), are going to be [1,2], the integer bit of the reciprocal is known to be 1, so one can use 8-bit table entries to produce a 1.8 fixed-point reciprocal by adding 0x100 (= 1). The reason 0x101 is used here is not clear; it may be due to a requirement that this step always provides an overestimate of the true reciprocal.
The next two steps are verbatim translations of the Newton-Raphson iteration for the reciprocal taking into account fixed-point scale factors. So 0x2000000 represents 2.0. The code uses 0x2000080 since it incorporates a rounding constant 0x80 (=128) for the following division by 256, used for rescaling the result. The next step likewise adds 0x00000080 as a rounding constant for a rescaling division by 256. Without the scaling, this would be a pure multiplication.
The final multiplication n*d multiplies the reciprocal in d with the dividend in n to yield the quotient in 33 bits. Again, a rounding constant of 0x8000 is applied before dividing by 65536 to rescale into the proper range, giving a quotient in 1.16 fixed-point format.
Continuous rescaling is typical of fixed-point computation where one tries to keep intermediate results as large as possible to maximize the accuracy of the final result. What is a bit unusual is that rounding is applied in all of the intermediate arithmetic, rather than just in the final step. Maybe it was necessary to achieve a specified level of accuracy.
The function is not all that accurate, though, likely caused by inaccuracy in the initial approximation. Across all non-exceptional cases, 2,424,807,756 match a correctly rounded 1.16 fixed-point result, 780,692,403 have an error of 1 ulp, 15,606,093 have a 2-ulp error, and 86,452 have a 3-ulp error. In a quick experiment, the maximum relative error in the initial approximation u was 3.89e-3. An improved table lookup reduces the maximum relative error in u to 2.85e-3, reducing but not eliminating 3-ulp errors in the final result.
If you want to look at a specific example, consider h =0.3 ( 0x4ccd ) divided by SZ3 =0.2 ( 0x3333 ). Then z =2, thus d =0.2*4 = 0.8 ( 0xcccc ). This leads to u = 1.25 ( 0x140 ). Since the estimate is quite accurate, we expect (2 - d * u) to be near 1, and in fact, d = 1.000015 ( 0x10001 ). The refined reciprocal comes out to d =1.250015 ( 0x14001 ), and quotient is therefore n =1.500031 ( 0x18002 ).
上一篇: Python轮子：不支持cp27mu
下一篇: 这种分割近似算法如何工作？