|
Оптимизация поворачивающих коэффициентовВеличина называется поворачивающим коэффициентом. В приведенной выше программе велико искушение использовать формулу: . (17) Это позволило бы избавиться от возведений в степень и обойтись одними умножениями, если заранее вычислить WN для N = 2, 4, 8,…, Nmax. Но то, что можно делать в математике, далеко не всегда можно делать в программах. По мере увеличения k поворачивающий коэффициент будет изменяться, но вместе с тем будет расти погрешность. Ведь вычисления с плавающей точкой на реальном компьютере совсем без погрешностей невозможны. И после N/2 подряд умножений в поворачивающем коэффициента может накопиться огромное отклонение от точного значения. Вспомним правило: при умножении величин их относительные погрешности складываются. Так что, если погрешность одного умношения равна s%, то после 1000 умножений она может достигнуть в худшем случае 1000s%. Этого можно было бы избежать, будь число WN целым, но оно не целое при N > 2, так как вычисляется через синус и косинус:
Как же избежать множества обращений к весьма медленным функциям синуса и косинуса? Здесь нам приходит на помощь давно известный алгоритм вычисления степени через многократные возведения в квадрат и умножения. Например:
В данном случае нам понадобилось всего 5 умножений (учитывая, что не нужно вычислять дважды) вместо 13. В худшем случае для возведения в степень от 1 до N/2-1 нужно log2N умножений вместо N/2, что дает вполне приемлемую погрешность для большинства практических задач. Можно сократить вдвое количество умножений на каждом шаге, если использовать результаты прошлых вычислений , для хранения которых нужно дополнительно log2(Nmax) комплексных ячеек памяти:
Если очередное не было вычислено предварительно, то берется двоичное представление k и анализируется. Каждому единичному биту соответствует ровно один множитель. В общем случае единице в бите с номером b (младший бит имеет номер 0) соответствует множитель , который хранится в b-й ячейке упомянутого выше массива. Есть способ уменьшить количество умножений для вычисления до одного на два цикла. Но для этого нужно отвести N/2 комплексных ячеек для хранения всех предыдущих . Алгоритм достаточно прост. Нечетные элементы вычисляются по формуле (17). Четные вычисляются по формуле:
То есть, ничего не вычисляется, а берется одно из значений, вычисленных на предыдущем шаге, когда N было вдвое меньше. Чтобы не нужно было копировать величины на новые позиции достаточно их сразу располагать в той позиции, которую они займут при N = Nmax и вводить простую поправку Skew (см. листинг программы). И последние пояснения относительно листинга. Во-первых, здесь присутствует реализация простейших операций над комплексными числами (классы Complex и ShortComplex), оптимизированная под данную задачу. Обычно та или иная реализация уже есть для многих компиляторов, так что можно использовать и ее. Во-вторых, массив W2n содержит заранее вычисленные коэффициенты W2, W4, W8,...,WNmax. В-третьих, для вычислений используются наиболее точное представление чисел с плавающей точкой в C++: long double размером в 10 байт на платформе Intel. Для хранения результатов в массивах используется тип double длиной 8 байт. Причина - не в желании сэкономить 20% памяти, а в том, что 64-битные и 32-битные процессоры лучше работают с выровненными по границе 8 байт данными, чем с выровненными по границе 10 байт.
|