Вычисление корня квадратного и инверсного корня в C28x Delfino

В этой статье расскажу и приведу много ссылок на методы вычисления квадратного корня. Захвачу тему вычисления инверсного корня. Разберемся, откуда берутся формулы для вычисления. Как работает метод Ньютона-Рафсона. Где для него взять первое приближение. И намекну, как это все устроено в математической FPU библиотеке Texas Instruments.

Вычисление квадратного корня.

В интернете существует множество реализаций вычисления квадратного корня. В статье  Best Square Root Method - Algorithm - Function (Precision VS Speed) проделана огромная работа для сравнения нескольких вариантов реализаций. На 1ом месте, естественно, код, использующий ассемблерную команду нахождения корня ( метод 14 на рисунке ниже). На втором месте расположился метод 7 - Float square root approximation - у него 5% ошибки.

Сравнивались все они со стандартной С++ функцией std::sqrt() из <cmath>.

Но в отношении процессоров C2000 32 bit 28x Delfino™ Floating-point Series (т.е. TMS320f28335, который имеет на борту FPU) - да, и процессоров других семейств от Texas’а с поддержкой FPU - дело обстоит запутанней. Давайте распутываться.

С чего все началось?

Все началось с желания сделать что-то свое или синтезировать в одном алгоритме несколько, реализовав все на С++ и заполучив кроcсплатформенный код без привязки к аппаратной начинке (т.е. к ассемблеру).

Найдена отличная статья Improve your root-mean calculations, которая заставила обратиться к математике и листочку с ручкой... В статье Фигура 2 описывает реализацию инверсного корня для случая нахождения  среднего m(n) значения от входных сэмплов x(n).

Статья хорошая, но как всегда без промежуточных выкладок для конечных формул. Исписав пару листов бумаги, стало ясно откуда “растут ноги”, как применять метод Ньютона-Рафсона. И пришло прозрение что за формулы фигурируют у Texas Instruments в их реализации floating-point math библиотеке. А знать это - как минимум интересно. И надо отметить, что TI на славу постарались. Библиотека у них быстрая.
В процессе создания собственной реализации математической библиотеки я застопорился на решении проблемы вычисления обратного корня только на том, откуда брать начальное приближение. При расчете начального приближения требуется воспользоваться упрощенной формулой искомой функции, но и они содержат “тяжеловесные” операции, например, корень или деление. Как оказалось, выходом может служить таблицы (Look-up tables) начальных приближений (как это делается смотрим в статье от Apple: Computing the Inverse Square Root с кодом от Apple =)  ).

Но раз уже реализовано в run-time support библиотеке ( rts2800_fpu32.lib, которая может быть скомбинирована для еще большей производительности с FPU FastRTS library: rts2800_fpu32_fast_supplement.lib   ) все необходимое, то нет смысла копать дальше в том же направлении, т.к. реализованная TI библиотека уж очень хороша. Для подтверждения слов смотрим статьи TMS320 c28x: Сравнение производительности стандартных функций и TMS320 c28x: Операция “остаток от деления” .

Как устроено у Texas Instruments?

В FPU реализовано много замечательных команд. Ознакомиться с ними можно в документе TMS320C28x Floating Point Unit and Instruction Set. Среди них есть две команды: 32-bit Floating-Point Reciprocal Approximation и 32-bit Floating-Point Square-Root Reciprocal Approximation. у=1/x  и у=1/sqrt(x). Результат выполнения команд - начальное приближение у0. Затем применяется  метод NR (newton raphson) для уточнения значения. Обычно хватает 3-4 итераций. Поэтому цикл не требуется, а лишь последовательно записывают одно и тоже рекуррентное соотношение метода NR (“разворачивают” цикл), опытным путем определив требуемое количество итераций. Смотрим команды EINVF32 и EISQRTF32.
EINVF32 32-bit Floating-Point Reciprocal Approximation. This operation generates an estimate of 1/X in 32-bit floating-point format accurate to approximately 8 bits. This value can be used in a Newton-Raphson algorithm to get a more accurate answer. That is:

                                      Ye = Estimate(1/X);  //первое приближение

Ye = Ye*(2.0 - Ye*X)

Ye = Ye*(2.0 - Ye*X)

After 2 iterations of the Newton-Raphson algorithm, you will get an exact answer accurate to the 32-bit floating-point format. On each iteration the mantissa bit accuracy approximately doubles.

EISQRTF32 32-bit Floating-Point Square-Root Reciprocal Approximation. This operation generates an estimate of 1/sqrt(X) in 32-bit floating-point format accurate to approximately 8 bits. This value can be used in a Newton-Raphson algorithm to get a more accurate answer. That is:

                                 Ye = Estimate(1/sqrt(X)); //первое приближение

Ye = Ye*(1.5 - Ye*Ye*X/2.0)

Ye = Ye*(1.5 - Ye*Ye*X/2.0)

Метод Ньютона-Рафсона.

Для нахождения квадратного корня сначала ищут точное значение обратного корня 1/sqrt(x), а затем умножают на подкоренное значение Х:

x/sqrt(x) = sqrt(x)

Почему вычисление корня квадратного производится через инверсный корень?

Если совсем по-простому метод NR основан на касательных, а касательная - это значение производной в точке касания.

Представим, что мы находимся в n-ой итерации, и она предпоследняя,т.е следующее найденное значение f(x) будет равно 0 (мы составляем f(x) так, чтобы искомое Х было корнем уравнения, т.е f(X)=0). Тогда:

Ну и получаем рекуррентное выражение для нахождения корня уравнения:

SQRT(X). Применим метод NR. Искомая функция:

Подставляем в (1):

Теперь по входным данным y(n) мы можем отыскать значения корня квадратного, но при реализации рекуррентного соотношения мы спотыкаемся об операцию деления, которая очень “прожорлива” ( убедиться в этом можно здесь - Сравнение производительности функций )

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

Сводим задачу к прошлому варианту (здесь мы можем переназначить переменные, но я не буду):

Подстановка в (1):

Выражение (3) содержит только умножения и сложения, что запросто можно реализовать на DSP. Для конечного результата нужно:

Требуется только начальное приближение у0.
Вот откуда получаем такую схему из статьи “Improve your root-mean calculations”.

Примечание! Данный алгоритм больше всего подходит для floating-point процессоров, которые эффективно справляются как с числами больше 1 так и много меньше 1. Fixed-poit не подходят, т.к. подкоренное значение а может быть больше/меньше 1, тогда y(n) может выйти за диапазон представления чисел и стать много меньше/больше 1.
Примечание!! Не все так безоблачно. Метод нестабилен, если начальное приближение далеко от правильного ответа. (Iterative methods for reciprocal square roots )
Выходит, нужно правильно выбирать начальное приближение, чтобы, во-первых, метод был устойчив, во-вторых, чтобы он быстро сходился к ответу. Т.е. нужно каким-то способом произвести аппроксимацию исходной функции y(n), например,линейно, как аппроксимируют синус малых углов sin(x)=x... Можно, конечно разложить в ряд Тейлора, но это не упростит вычислительные нагрузки на алгоритм для данных функций

Выбор начального приближения.

Выбор начального приближения - это наиболее сложная часть всей работы и самая ответственная.

Если будете реализовывать описанные в этой статье методы вычислений, то стоит сперва проштудировать следующие материалы: Cascaded Implementation of an Iterative Inverse-Square-Root  Algorithm, with Overflow Lookahead . Здесь на стр. 117 пишут: “There are two methods to calculate the initial values for the first stage”. Так что читаем дальше по тексту.

Про выбор стартовой точки пишут так же в статье Choosing Starting Values for Newton-Raphson Computation of Reciprocals, Square-Roots and Square-Root Reciprocals: особенно интересны разделы   2.1 Choosing the Best Starting Point и  4.1 Getting a Good Starting Point.
А так же не забываем про статью от Apple про табличную реализацию выбора начального приближения. Сосредоточиться необходимо именно на этом методе, т.к. в C28x FPU Fast RTS Library используют подобный метод. Таблицы эти (fixed-point and
floating-point math tables) находятся в Boot ROM, который прошивается производителем. Они необходимы для использования  библиотек C28x™ IQMath Library   и  C28x FPU Fast RTS Library . Приведу карту памяти для Boot ROM:

P.S. Сейчас нахожусь в стадии осмысления работы математической FPU библиотеки Texas Instruments. Постараюсь основные моменты описать в статье. Так что следите за обновлениями в Блоге!!