Оптимизация поворачивающих коэффициентов

Величина называется поворачивающим коэффициентом. В приведенной выше программе велико искушение использовать формулу:

.    (17)

Это позволило бы избавиться от возведений в степень и обойтись одними умножениями, если заранее вычислить W_N для N = 2, 4, 8,…, Nmax. Но то, что можно делать в математике, далеко не всегда можно делать в программах. По мере увеличения k поворачивающий коэффициент будет изменяться, но вместе с тем будет расти погрешность. Ведь вычисления с плавающей точкой на реальном компьютере совсем без погрешностей невозможны. И после N/2 подряд умножений в поворачивающем коэффициента может накопиться огромное отклонение от точного значения. Вспомним правило: при умножении величин их относительные погрешности складываются. Так что, если погрешность одного умношения равна s%, то после 1000 умножений она может достигнуть в худшем случае 1000s%.

Этого можно было бы избежать, будь число W_N целым, но оно не целое при N > 2, так как вычисляется через синус и косинус:

Как же избежать множества обращений к весьма медленным функциям синуса и косинуса? Здесь нам приходит на помощь давно известный алгоритм вычисления степени через многократные возведения в квадрат и умножения. Например:

В данном случае нам понадобилось всего 5 умножений (учитывая, что  не нужно вычислять дважды) вместо 13. В худшем случае для возведения в степень от 1 до N/2-1 нужно log₂N умножений вместо N/2, что дает вполне приемлемую погрешность для большинства практических задач.

Можно сократить вдвое количество умножений на каждом шаге, если использовать результаты прошлых вычислений

,

для хранения которых нужно дополнительно log₂(Nmax) комплексных ячеек памяти:

Если очередное  не было вычислено предварительно, то берется двоичное представление k и анализируется. Каждому единичному биту соответствует ровно один множитель. В общем случае единице в бите с номером b (младший бит имеет номер 0) соответствует множитель , который хранится в b-й ячейке упомянутого выше массива.

Есть способ уменьшить количество умножений для вычисления  до одного на два цикла. Но для этого нужно отвести N/2 комплексных ячеек для хранения всех предыдущих . Алгоритм достаточно прост. Нечетные элементы вычисляются по формуле (17). Четные вычисляются по формуле:

То есть, ничего не вычисляется, а берется одно из значений, вычисленных на предыдущем шаге, когда N было вдвое меньше. Чтобы не нужно было копировать величины на новые позиции достаточно их сразу располагать в той позиции, которую они займут при N = Nmax и вводить простую поправку Skew (см. листинг программы).

И последние пояснения относительно листинга.

Во-первых, здесь присутствует реализация простейших операций над комплексными числами (классы Complex и ShortComplex), оптимизированная под данную задачу. Обычно та или иная реализация уже есть для многих компиляторов, так что можно использовать и ее.

Во-вторых, массив W2n содержит заранее вычисленные коэффициенты W₂, W₄, W₈,...,W_Nmax.

В-третьих, для вычислений используются наиболее точное представление чисел с плавающей точкой в C++: long double размером в 10 байт на платформе Intel. Для хранения результатов в массивах используется тип double длиной 8 байт. Причина - не в желании сэкономить 20% памяти, а в том, что 64-битные и 32-битные процессоры лучше работают с выровненными по границе 8 байт данными, чем с выровненными по границе 10 байт.