Вычисление модуля вектора ( sqrt( I^2 + Q^2) ): от идеи до реализации

В статье собрана и проанализирована информация о методах вычисления корня квадратного, приведены ссылки на материалы. Сделана попытка синтезировать на основе двух методов алгоритм нахождения модуля вектора: оценка модуля вектора алгоритмом “альфа максимум плюс бета минимум” и последующее уточнение значения с помощью метода Ньютона.

           Статья 2. Вычисление модуля вектора ( sqrt( I^2 + Q^2 ) ). Результаты
          Статья 3. Вычисление модуля вектора (работа над ошибками)

Ссылки по вычислению квадратного корня:
вычисление квадратного корня на tms320f2812
Очень нужна помощь в оптимизации алгоритма извлечения квадратного корня
Корень квадратный, Как взять быстро?
Квадратный корень в целых числах, Benchmark for ARM
Magical Square Root Implementation In Quake III

Часть 1.

Полезная добытая информация из статей.

IQMath - это библиотека TI только для работы с fixed-point процессорами. А мой процессор, TMS320F28335 - floating-point. И эта библиотека для него не подходит.

Для floating-point процессоров есть RTS Library, в которой все оптимизировано: sprc664 –это про fastRTS, sprc624 – это про “не ANSI, но быстро” (про fpu библиотеку)
В частности в SPRC664 сказано:


Примерный список методов извлечения корня приведем ниже ( ссылка ):

1.  Арифметический

2.  Грубая оценка

3.  Столбиком

4.  Вавилонский способ  

5.  Метод Герона

6.  Метод Ньютона

7.  Десятично

Арифметический способ

Для квадратов чисел верны следующие равенства:

1 = 12
1 + 3 = 22
1 + 3 + 5 = 32
и так далее.

То  есть,  узнать  целую  часть  квадратного  корня  числа  можно,  вычитая  из него все нечётные числа по порядку, пока остаток не станет меньше следующего вычитаемого  числа  или  равен  нулю,  и  посчитав  количество  выполненных действий. Например, так:

9 − 1 = 8
8 − 3 = 5
5 − 5 = 0

Выполнено 3 действия, квадратный корень числа 9 равен 3.

Недостатком  такого  способа  является  то,  что  если извлекаемый  корень  не  является  целым  числом,  то  можно  узнать  только  его

целую часть, но не точнее. В то же время такой способ вполне доступен детям, решающим  простейшие  математические  задачи,  требующие  извлечения квадратного корня.

Кусок из Уоррен Генри. Алгоритмические трюки для программистов, для 32-битного исходного числа, "в среднем ему требуется 149 базовых

RISC-команд" (в книге целый раздел, посвященный извлечению корня!)

Часть 2.

Идея.

Можно запрограммировать метод Ньютона (он же очень похож на вавилонский), но в качестве первого приближения использовать алгоритм αMax + βMin ( “ альфа максимум плюс бета минимум” )

Что значит вычислить корень из х?

Если под корнем стоит некое число а, то мы должны найти такое число g, чтобы:

Или по-другому :

Напомним, что g - это искомое значение корня, а а - это подкоренное число (или число, передаваемое в виде аргумента в функцию, например, sqrt( a ) для нахождения значения квадратного корня из а).

Уравнение (3) переписывается иначе:

Уравнение (4) полезно как выражение, задающее условие продолжения работы цикла в программе С++ в неком итерационном подходе вычисления g (этот итерационной подход - вычисление квадратного корня методом Ньютона). Пока мы не найдем точное значение g, значения gi и a/gi не будут совпадать, поэтому точное значение корня квадратного g будет находиться между двумя этими числами. Тогда логично, что следующее наилучшее приближение для  gi+1 будет среднее арифметическое gi и a/gi. Из этих соображений и вырос, видимо, метод Ньютона:

Вопрос: что выбирать в качестве первого приближения g0 для алгоритма Ньютона? Много методов в книгах описано. Но мне хочется попробовать своеобразный способ, мало описанный в литературе. И главное, что этот способ идеально подходит нам, т.к. ожидается, что сделает наиболее близкое приближение к истинному квадратному корню из числа, тогда потребуется минимум итераций метода Ньютона, чтобы достичь истинное значение, т.к. двигаться мы будем не от целого начального приближения (например, для х=5 начальное приближение g0 может быть 4), а от вещественного первого приближения, например, 2.1 . Алгоритм называется “альфа максимум плюс бета минимум”. Он нам подходит еще и потому, что оперирует в понятиях векторного числа, у которого есть синфазная и квадратурная составляющие, т.е. де-факто I и Q составляющие. Фактически, это метод вычисления модуля комплексного числа. Достаточно точный метод.

Часть 3.

От идеи к алгоритму.

Часть 4

Реализация.

2. #ifndef SQRTIQ_HPP
3. #define SQRTIQ_HPP
5. #ifdef ITERATIONS
6. //for debug
7. int iterationCnt;
8. #endif
10. template< class T >
11. T sqrtIQ( T i, T q ) {
13. #ifdef ITERATIONS
14. //for debug
15. iterationCnt = 1;
16. #endif
18. //first estimation for Newton's algorithm. This estimation is found
19. //by means of "alpha max plus beta min" algorithm
20. T firstEstimation;
21. //current value of estimation
22. T nextEstimation;
23. T radicand;
24. T quotient;
25. //coefficients for "alpha max plus beta min":
26. //alpha = 1.0, beta = 0.5 therefore they don't need to be present in memory
28. //begin - "alpha max plus beta min" algorithm
29. if( i > q ) {
30. firstEstimation = i +  q/2;
31. }
32. else {
33. firstEstimation = q +  i/2;
34. }
35. //end - "alpha max plus beta min" algorithm
36. //begin - Newton's algorithm
37. i *= i;
38. q *= q;
39. radicand = i + q;
41. quotient = radicand/firstEstimation;
42. nextEstimation = ( firstEstimation + quotient )/2;
43. while( firstEstimation - nextEstimation != 0 ) {
44. firstEstimation = nextEstimation;
46. nextEstimation = ( firstEstimation + quotient )/2;
48. //for next iteration check in while()
49. quotient = radicand/firstEstimation;
51. #ifdef ITERATIONS
52. //for debug
53. ++iterationCnt;
54. #endif
55. } //while
57. return firstEstimation;
59. //end - Newton's algorithm
60. }
63. #endif // SQRTIQ_HPP