nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Обратный квадратный корень в целых числах

У нас есть процессор, который умеет складывать-вычитать, а также умножать. Также в арсенале - сдвиги влево и вправо, которыми можно заменить умножение и деление на степени двойки, а ещё есть взятие модуля числа и нахождение максимума двух чисел. Ничего другого - деления, квадратного корня, тригонометрии - в самом АЛУ процессора НЕТ.

Как оказалось, один из двух наших алгоритмов, "алгоритм захвата", он же "аффинный алгоритм", замечательно реализуется на таком процессоре.

Для решения системы линейных уравнений (самый напряжённый момент второго алгоритма, "алгоритма сопровождения") мне очень хотелось применить метод LDLT или UDUT-разложения, поскольку тогда мы обходимся чистой арифметикой. Деление я собирался либо реализовать аппаратно, либо опять свести к умножению на обратную величину, а её находить по методу Ньютона. Увы, оказалось, что эти методы приводят к чрезмерно большим разбросам диагональных элементов, независимо от того, как я попытаюсь отмасштабировать входные данные. В итоге 16 бит оказывается недостаточно.

Метод UUT оказался куда более многообещающим - благодаря применению квадратного корня удалось сократить разброс диагональных элементов, и теперь система линейных уравнений решалась нормально, не приводя ни к переполнению, ни к катастрофическому падению точности.

Кроме того, оказалось, что можно все деления и взятия квадратных корней заменить одной-единственной операцией (в дополнение к сложениям и умножениям) - это взятие "обратного квадратного корня", т.е



Тогда на этапе разложения мы можем сразу заменить диагональные элементы на обратные, и при нахождении остальных элементов надо будет не делить на элементы главной диагонали, а умножать на них. И точно также все деления заменятся умножениями на этапах прямой и обратной подстановки.

У многих на слуху "магическая" функция нахождения обратного квадратного корня из недр Quake, которую иногда приписывают самому Кармаку, но скорее всего это было "совместное творчество", и найти концы уже не так-то просто. Я раньше думал, что это просто какая-то аппроксимация по нескольким точкам, которая если даст плюс-минус 10% - то и ладно. Оказалось, точность этого метода - лучше 0,17% во всём диапазоне входных данных, при том, что метод сводится к сдвигу вправо, вычитанию из "магической константы" и одной итерации метода Ньютона. Очень рекомендую ознакомиться с подробностями (раз, два) - они того стоят.

Но к сожалению, эта магическая функция целиком опирается на представление числа с плавающей точкой, на тот факт, что после экспоненты идёт мантисса. Именно рассмотрение их как единого целого, 32-битного числа, позволило добиться столь впечатляющих результатов. Для нас это бесполезно.

Поэтому нужно было найти свой способ извлечения обратного квадратного корня из целого числа, причём требования весьма специфичны:
1. точность - максимально возможная для 16-битного результата,
2. размер функции - настолько малый, насколько это возможно, т.к памяти очень мало
3. быстродействие - практически пофиг. Нам этих корней понадобится 5-6 на одно измерение, т.е не больше 30 корней в секунду. Прямой перебор (проверить все 65535 разных значений) - это, пожалуй, перебор, а вот сделать десяток-другой итераций кажется вполне осмысленным.

В итоге получилась функция, которая вычисляет результат за 16 итераций метода Ньютона, на каждой итерации осуществляется 3 умножения и 1 вычитание. Ошибка составляет одну единицу младшего разряда, например, реальный ответ должен составлять 1023,5003659, что должно бы округлиться до 1024, а наша функция вместо этого получает 1023. Так происходит при 106 входных значениях из возможных 65535. Но начало работы функции навевает вселенскую тоску...


У нас есть число a, от которого надо взять обратный квадратный корень. Для этого мы берём нулевое приближение:



после чего выполняем итерации



Вот и весь метод. "Дьявол кроется в деталях". Входное число мы воспринимаем как беззнаковое в формате 1.15, т.е значения от 0 до 65535 мы воспринимаем как 0 .. 1,999969.
При этом нулевое значение недопустимо - это будет "деление на ноль". Все остальные мы хотим обработать без переполнения, потому как в реальной работе мы действительно встречаемся с широким диапазоном значений.

Самое большое выходное значение получится, когда на вход подано самое маленькое, "1", т.е 1/32768. Имеем 1/sqrt(1/32768)≈181,01934 Чтобы оно поместилось в результат, мы обязаны представить данные в формате 8.8, т.е 8 бит целой части и 8 бит дробной, что позволяет представить числа от 0 до примерно 256. В таком случае наше значение 181,01934 будет представлено как 256 * 181,01934 = 46341.

Самое маленькое выходное значение получится, когда на вход подано самое большое значение "65535" т.е 65535/32768≈2. Имеем 1/sqrt(65535/32768)≈0,70711. В формате 8.8 оно будет представлено как 256*0,70711=181. По счастью, разрядности ещё хватает, чтобы получить сколько-нибудь точный результат. Наихудшая относительная погрешность получается для входного значения 65189, что должно было дать результат 181,5005, но округляется до 182, приводя к погрешности 0,27%. Нормально...

Теперь о выборе нулевого приближения. При записи "в действительных числах", результат может получаться в диапазоне от 0,707 до 181, и нам может захотеться взять значение где-то посерединке, будь то среднее арифметическое (90,86) или среднее геометрическое (11,3). Увы, ни то, ни другое взять мы не имеем права - для многих значений результат попросту разойдётся.

Если нулевое приближение превысит , то значение в скобках (в итерации метода Ньютона) станет ниже минус единицы, т.е результат по модулю возрастёт и при этом поменяет знак. В итоге мы получим расходящуюся последовательность.

Если нулевое приближение будет в интервале от до , то значение в скобках станет от минус единицы до нуля, т.е x1 поменяет знак, но снизится по модулю. Последовательность может с большой вероятности сойтись к , тогда как мы хотели просто . Либо оно может уйти в ноль, после чего так и останется в нуле.

Поэтому, для диапазона входных чисел от 1/32768 до 2, мы должны выбрать нулевое приближение не более

Довольно разумное значение для нулевого приближения - единица, или в нашем формате 8.8: 256. В таком случае, мы можем обещать, что при умножении a на x результат уместится в 24 бита, что может выцарапать нам немножко точности.

Как ни странно, при таком выборе начального значения, наш метод всегда будет сходиться к правильному ответу, другое дело, что иногда ОЧЕНЬ МЕДЛЕННО. К примеру, попробуем посчитать обратный квадратный корень из 1/32768 (это и есть наихудший случай)

Имеем x0=1.

Далее, x1=1 * (3/2 - 1/2/32768)≈1,49998;
x2=1,49998 * (3/2 - 1/2/32768*1,499982)≈2,24992;
x3≈3,37471;
x4≈5,06148;
x5≈7,59025;
x6≈11,37870;
x7≈17,04557;
x8≈25,49279;
x9≈37,98639;
x10≈56,14320;
x11≈81,51450;
x12≈114,00710;
x13≈148,39985;
x14≈172,73195;
x15≈180,45890;
x16≈181,01673;

Метод на первых порах получается даже хуже старого доброго "деления пополам" (поиска льва в пустыне). Там мы, по сути, получаем по одному значащему биту за раз, снижая неопределённость вдвое. Здесь мы долго и упорно выползаем из угла, каждый раз увеличивая начальное приближение не более чем в 1,5 раза. И только когда мы начинаем приближаться к правильному ответу, спустя 14-15 итераций, внезапно начинает проявляться квадратичная сходимость метода Ньютона, где количество правильных бит ответа удваивается на каждой итерации...

Мы могли бы существенно ускорить этот метод, если начать с битовых операций: сдвигать исходное значение влево на 2 бита, пока не наступает переполнения, и вместе с этим сдвигать конечный результат вправо на 1 бит. Тогда наше нулевое приближение будет содержать хотя бы одно верное значение, и сойдётся моментально. Ещё круче было бы ввести какую-то специальную аппаратную инструкцию типа "целочисленного логарифма", которая подсчитывает, сколько нулей в старших разрядах числа перед первой единицей. Но это будет "преждевременная оптимизация" - пока мне кажется, что сильно переживать за быстродействие данной функции нам не стоит!

А вот с точностью вычислений нужно постараться! Если каждое промежуточное значение будет округляться до 16 бит, то время от времени функция начнёт "ходить кругами". Например, если взять число 6 на входе (т.е 6/32768), то в какой-то момент начнётся последовательность:
18922,
18929,
18911,
18924,
18908,
18922, ...

тогда как правильный ответ: 18919.

Максимальная ошибка (для входных значений от 1 до 65535) будет составлять 12 единиц (если сделать 20 итераций), что не есть хорошо.

Поэтому придётся несколько дополнить наш АЛУ командой, которая переносит 32-битный аккумулятор в 31-битный регистр B. До этого мы могли занести в него только старшие 16 бит, а остальное было нужно, чтобы в процессе умножения двух чисел его содержимое сдвигалось вправо и при единичном бите второго операнда прибавлялось к аккумулятору. Но если мы позволим производить такую замену, точность в кои-то веки дойдёт до одной единицы младшего разряда, а на самом деле и того лучше.

Конечно, ради одной этой функции такое дополнение АЛУ может показаться расточительным, уж лучше придумать другой алгоритм, да хоть того же льва, он точно не зациклится. Но как показала практика, и последующее решение системы линейных уравнений только в том случае даёт неплохие результаты, если и там мы не сбрасываем 32-битное промежуточное значение до 16 бит, прежде чем помножить на элемент главной диагонали. В противном случае - огромные ошибки и срыв сопровождения...
Tags: ПЛИС, кватернионы-это просто (том 1), математика, работа, странные девайсы
Subscribe

  • Я создал монстра!

    Вот нормальная счастливая пара разъёмов ОНЦ-БС-1-10/14-Р12-2-В и ОНЦ-БС-1-10/14-В1-2-В: У розетки кроме основного выступа, отмечающего "верх",…

  • Нахождение двух самых отдалённых точек

    Пока компьютер долго и упорно мучал симуляцию, я пытался написать на ассемблере алгоритм захвата на ближней дистанции. А сейчас на этом коде можно…

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments