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

Question

Я работаю над игрой с программным средством визуализации, чтобы получить наиболее точный вид PS1. Пока я исследовал, как работает система графики / рендеринга PS1, причины шатких вершин и т. Д., Я натолкнулся на некоторую документацию о том, как они делят. Вот ссылка на него: http://problemkaputt.de/psx-spx.htm (см. Раздел "Неточность деления 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

Этот алгоритм представляет собой деление с фиксированной запятой двух беззнаковых 16-битных дробных значений в [0,1). Сначала он вычисляет начальное 9-битное приближение к обратной стороне делителя с помощью поиска в таблице, уточняет это с помощью одной итерации Ньютона-Рафсона для обратной величины, 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 244 807 756 соответствуют правильно округленному результату с фиксированной точкой 1,16, 780 692 403 имеют ошибку в 1 ульт, 15 660 093 имеют ошибку 2 ульт, а 86 452 имеют ошибку 3 ульт. В быстром эксперименте максимальная относительная погрешность в начальном приближении 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).