Как работает алгоритм аппроксимации деления?

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 

Я с трудом понять, как это работает, что это «УНР ' Таблица? почему мы меняем вещи? Если кто-то может дать более подробное объяснение того, как эта вещь действительно достигает разрыва, это будет оценено.

Answer 1

Этот алгоритм представляет собой разделение с фиксированной точкой двух беззнаковых 16-битовых дробных значений в [0,1]. Сначала он вычисляет начальное 9-битное приближение к обратному делителю посредством поиска таблицы, уточняет это с помощью одной итерации Newton-Raphson для обратного, x _{i + 1}: = x _i * (2 - d * x _i), в результате чего получается обратная точность примерно до 16 бит, и, наконец, умножает это на дивиденд, получая 17 битов в [0,2).

Для поиска в таблице делитель сначала нормализуется на [0,5, 1), применяя масштабный коэффициент 2 ^z, очевидно, что тогда дивиденд должен быть скорректирован одним и тем же масштабным коэффициентом. Поскольку обратные операнды в [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 бит. Опять же, округляющая константа 0x8000 применяется до деления на 65536, чтобы перемасштабировать в соответствующий диапазон, давая коэффициент в 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).