УТ БПФ - раздел 3, часть 1
Раздел 3: практическая реализация УТ БПФ
Часть 1: быстрая троичная инверсия
Дан массив x[k], -T≤k≤T, T=(N-1)/2, N=3q. Требуется переставить в нем значения таким образом, что значение с индексом k=(kq-1kq-2…k1k0)ут будет иметь индекс n=(k0k1…kq-2kq-1)ут.
Исходный и обработанный массивы показаны на рис 3.1
Можно заметить, что некоторые элементы остаются на своих местах, это своего рода «числа-палиндромы», которые одинаково читаются слева направо и справа налево. Все остальные элементы попарно меняются местами.
Отметим также, что либо меняются местами элементы i>0 и j<0, либо две пары: i>0,j>0 и (-i),(-j).
Приведем (в псевдокоде) общий алгоритм такой перестановки.
Функция inversion находит следующее значение переменной n, которая на каждом шаге должна быть инверсией от переменной j. Самый простой, но медленный способ ее реализации состоит в объявлении массива из q целых чисел, trit[], которые будут изображать отдельные триты в троичном представлении числа. При каждом вызове inversion, зная, что j увеличится на единицу, это прибавление единицы будет произведено в троичном представлении. n будет найдено по формуле
Именно это будет наиболее медленная часть данного алгоритма, давая время выполнения O(q) на каждой итерации, что приводит к оценке сложности всего алгоритма O(Nq)=O(NlogN).
Очевидно, что если предстоит множество однотипных преобразований с одинаковым N, можно объявить массив inversion, который достаточно вычислить только один раз. Так мы получим один из наиболее быстрых алгоритмов троичной инверсии.
Покажем еще один способ выполнить троичную инверсию за O(N) шагов, не требующий хранения каких-либо заранее вычисленных значений. Его двоичный аналог называется алгоритмом Голда-Рейдера [Brian Gough - FFT Algorithms - Computational Biology Research Center, 1997].
Идея состоит в том, чтобы зная инверсию от j, найти инверсию от j+1, не прибегая явно к битовому (тритовому) представлению числа. Проиллюстрируем ее на примере q=3. Начинаем с j=0, n=0 (инверсия нуля всегда будет нулем). Далее, находим инверсию от j=1.
На каждом шаге мы не имеем никакой дополнительной информации, кроме предыдущего значения n. Попробуем прибавить единицу к младшему триту, что аналогично прибавлению 3q-1 к n. Получаем n=9. Действительно, 001ут=1, а 100ут=9.
Следующий шаг: ищем инверсию от j=2. Снова прибавляя 3q-1 к n, получим n=18, но это заведомо неверный ответ, ведь получаться должны числа от -T до T, где T=(3q-1)/2=13. Это означает, что мы не сделали перенос в следующий разряд, получили 1+1=2, вместо 1+1=11. Сделаем его. Замена младшего трита с 2 на 1 эквивалентно вычитанию 3q, а прибавление единицы к следующему разряду добавляет к инверсированному значению 3q-2. Получаем, таким образом, n=18-27+3=-6. Проверим: j=2=011,а n=-6=110.
Шаги 3 и 4 будут похожими - получим 3, затем 12. На шаге 5 прибавим 9, получим 21 - слишком много, значит, был перенос из младшего разряда. Учтя его, получаем -3 - как будто бы правильное значение. Однако после первого переноса очевидно, что младший трит равен -1. Максимальное n, которое могло бы при этом получиться: -1*9+1*3+1*1=-5. У нас же получилось -3, значит произошел следующий перенос, теперь уже в старший разряд. Коррекция производится аналогично, только вычитается 3q-1, а прибавляется 3q-3.
Теперь можно сформулировать алгоритм, из текущего значения n (инверсии от j), вычисляющий его следующее значение (инверсию от j+1).
Здесь m,k-вспомогательные переменные.
Приведем также полный алгоритм быстрой троичной инверсии, использующий метод Голда-Рейдера:
Заметим, что цикл начинается с 1, поскольку нулевой отсчет останется на месте. Заканчивается же он на T-2, поскольку последний отсчет имеет в троичном представлении одни единицы и останется на месте, а предпоследний мы уже заведомо обработали ранее. Отсчет T-2 же всегда подлежит обработке, поскольку его младший трит отрицателен.
Оценим сложность данного алгоритма. По крайней мере один перенос будет наблюдаться в (T-1)/3=(3q-1-1)/2 итерациях. По крайней мере два переноса будут лишь в ((T-1)/3-1)/3=(3q-2-1)/2 итерациях. Наконец, перенос в старший разряд произойдет лишь один раз, изменив его с 0 на 1. Полное число итераций внутреннего цикла составит:
Таким образом, алгоритм требует O(N) арифметических операций. Учитывая, что здесь необходима лишь целочисленная арифметика, а N в реальных задачах велико, можно заключить, что время работы процедуры инверсии пренебрежимо мало по сравнению с временем работы непосредственно УТ БПФ.
Поэтому именно данный модифицированный алгоритм Голда-Рейдера рекомендуется к применению. Он прост, не требует дополнительной памяти для хранения промежуточных результатов и достаточно быстрый.
Сравнение различных алгоритмов на компьютере показало, что инверсия с явной записью тритов при N=6561 выполняется в 6 раз медленнее, чем по алгоритму Голда-Рейдера. Также был рассмотрен вариант алгоритма Голда-Рейдера, в котором не используется симметрия между положительными и отрицательными отсчетами. Он выполняется на 33% медленнее, чем приведенный выше.
Надо полагать, существуют и более быстрые методы, но поиск наиболее быстрого из них представляет, скорее, академический интерес.
Когда я начинал работать над этой темой, у меня было стойкое ощущение, что на двоичном компьютере гораздо проще выполнить двоичную же инверсию с помощью каких-то битовых операций, а троичная будет считаться через пень-колоду. Оказалось же, что большинство алгоритмов двоичной инверсии используют двоичность компьютера для того лишь, чтобы умножение и деление на 2 заменить сдвигом. Как показано в этом разделе, уравновешенная троичная инверсия выполняется ничуть не хуже и точно так же занимает пренебрежимо мало времени по сравнению с основными операциями БПФ.
Часть 1: быстрая троичная инверсия
Дан массив x[k], -T≤k≤T, T=(N-1)/2, N=3q. Требуется переставить в нем значения таким образом, что значение с индексом k=(kq-1kq-2…k1k0)ут будет иметь индекс n=(k0k1…kq-2kq-1)ут.
Исходный и обработанный массивы показаны на рис 3.1
Можно заметить, что некоторые элементы остаются на своих местах, это своего рода «числа-палиндромы», которые одинаково читаются слева направо и справа налево. Все остальные элементы попарно меняются местами.
Отметим также, что либо меняются местами элементы i>0 и j<0, либо две пары: i>0,j>0 и (-i),(-j).
Приведем (в псевдокоде) общий алгоритм такой перестановки.
Функция inversion находит следующее значение переменной n, которая на каждом шаге должна быть инверсией от переменной j. Самый простой, но медленный способ ее реализации состоит в объявлении массива из q целых чисел, trit[], которые будут изображать отдельные триты в троичном представлении числа. При каждом вызове inversion, зная, что j увеличится на единицу, это прибавление единицы будет произведено в троичном представлении. n будет найдено по формуле
Именно это будет наиболее медленная часть данного алгоритма, давая время выполнения O(q) на каждой итерации, что приводит к оценке сложности всего алгоритма O(Nq)=O(NlogN).
Очевидно, что если предстоит множество однотипных преобразований с одинаковым N, можно объявить массив inversion, который достаточно вычислить только один раз. Так мы получим один из наиболее быстрых алгоритмов троичной инверсии.
Покажем еще один способ выполнить троичную инверсию за O(N) шагов, не требующий хранения каких-либо заранее вычисленных значений. Его двоичный аналог называется алгоритмом Голда-Рейдера [Brian Gough - FFT Algorithms - Computational Biology Research Center, 1997].
Идея состоит в том, чтобы зная инверсию от j, найти инверсию от j+1, не прибегая явно к битовому (тритовому) представлению числа. Проиллюстрируем ее на примере q=3. Начинаем с j=0, n=0 (инверсия нуля всегда будет нулем). Далее, находим инверсию от j=1.
На каждом шаге мы не имеем никакой дополнительной информации, кроме предыдущего значения n. Попробуем прибавить единицу к младшему триту, что аналогично прибавлению 3q-1 к n. Получаем n=9. Действительно, 001ут=1, а 100ут=9.
Следующий шаг: ищем инверсию от j=2. Снова прибавляя 3q-1 к n, получим n=18, но это заведомо неверный ответ, ведь получаться должны числа от -T до T, где T=(3q-1)/2=13. Это означает, что мы не сделали перенос в следующий разряд, получили 1+1=2, вместо 1+1=11. Сделаем его. Замена младшего трита с 2 на 1 эквивалентно вычитанию 3q, а прибавление единицы к следующему разряду добавляет к инверсированному значению 3q-2. Получаем, таким образом, n=18-27+3=-6. Проверим: j=2=011,а n=-6=110.
Шаги 3 и 4 будут похожими - получим 3, затем 12. На шаге 5 прибавим 9, получим 21 - слишком много, значит, был перенос из младшего разряда. Учтя его, получаем -3 - как будто бы правильное значение. Однако после первого переноса очевидно, что младший трит равен -1. Максимальное n, которое могло бы при этом получиться: -1*9+1*3+1*1=-5. У нас же получилось -3, значит произошел следующий перенос, теперь уже в старший разряд. Коррекция производится аналогично, только вычитается 3q-1, а прибавляется 3q-3.
Теперь можно сформулировать алгоритм, из текущего значения n (инверсии от j), вычисляющий его следующее значение (инверсию от j+1).
Здесь m,k-вспомогательные переменные.
Приведем также полный алгоритм быстрой троичной инверсии, использующий метод Голда-Рейдера:
Заметим, что цикл начинается с 1, поскольку нулевой отсчет останется на месте. Заканчивается же он на T-2, поскольку последний отсчет имеет в троичном представлении одни единицы и останется на месте, а предпоследний мы уже заведомо обработали ранее. Отсчет T-2 же всегда подлежит обработке, поскольку его младший трит отрицателен.
Оценим сложность данного алгоритма. По крайней мере один перенос будет наблюдаться в (T-1)/3=(3q-1-1)/2 итерациях. По крайней мере два переноса будут лишь в ((T-1)/3-1)/3=(3q-2-1)/2 итерациях. Наконец, перенос в старший разряд произойдет лишь один раз, изменив его с 0 на 1. Полное число итераций внутреннего цикла составит:
Таким образом, алгоритм требует O(N) арифметических операций. Учитывая, что здесь необходима лишь целочисленная арифметика, а N в реальных задачах велико, можно заключить, что время работы процедуры инверсии пренебрежимо мало по сравнению с временем работы непосредственно УТ БПФ.
Поэтому именно данный модифицированный алгоритм Голда-Рейдера рекомендуется к применению. Он прост, не требует дополнительной памяти для хранения промежуточных результатов и достаточно быстрый.
Сравнение различных алгоритмов на компьютере показало, что инверсия с явной записью тритов при N=6561 выполняется в 6 раз медленнее, чем по алгоритму Голда-Рейдера. Также был рассмотрен вариант алгоритма Голда-Рейдера, в котором не используется симметрия между положительными и отрицательными отсчетами. Он выполняется на 33% медленнее, чем приведенный выше.
Надо полагать, существуют и более быстрые методы, но поиск наиболее быстрого из них представляет, скорее, академический интерес.
Когда я начинал работать над этой темой, у меня было стойкое ощущение, что на двоичном компьютере гораздо проще выполнить двоичную же инверсию с помощью каких-то битовых операций, а троичная будет считаться через пень-колоду. Оказалось же, что большинство алгоритмов двоичной инверсии используют двоичность компьютера для того лишь, чтобы умножение и деление на 2 заменить сдвигом. Как показано в этом разделе, уравновешенная троичная инверсия выполняется ничуть не хуже и точно так же занимает пренебрежимо мало времени по сравнению с основными операциями БПФ.