Блог пользователя gerasimovd

Автор gerasimovd, 13 лет назад, По-русски

Захотелось научиться писать БПФ без арифметики с плавающей точкой - начал изучать по этой статье. Вроде все более-менее понял, написал и даже кое-что работает, но остается непонятными несколько вещей:

Во-первых - пусть A = B*C(в смысле полиномов) и мы посчитали FFT(B), FFT(C). Все значения в этих массивах меньше модуля. Потом мы должны скалярно перемножить (FFT(B), FFT(C)) и применить обратное FFT. Но когда мы перемножим - не все элементы этого массива уже будут меньше модуля - что с этим делать? просто брать их по тому же модулю и считать обратное? Если я так делаю, у меня начинает все неправильно работать(из-за таких переполнений), числа длиной около тысячи знаков умножаются правильно только при базе 10 или 100, что очень печально.
Во-вторых - пусть после FFT у нас получилось: FFT(A) = (1, 3), FFT(B) = (4, 0). После того как мы скалярно перемножим, у нас получится (4, 0), и после обратного FFT мы получим тот же B. Такое бывает(может я не понял чего-то в теории), и если да - что делать?
И последнее - модуль в статье на e-maxx подходит для случая, когда в векторе не более 2^20 элементов. Вряд ли мне понадобится больше - но вопрос - а где-то есть готовые последовательности типа "модуль - образующий элемент" ?. На oeis - только для небольших чисел, можно конечно прочитать и познать еще один алгоритм - но очень уж лень :)
Кстати - мне кажется, или на e-maxx единственная статья про целочисленное БПФ в интернете? Если есть еще что почитать(на русском или английском), очень буду рад ссылкам.
upd: еще на всякий случай код
  • Проголосовать: нравится
  • +19
  • Проголосовать: не нравится

13 лет назад, # |
  Проголосовать: нравится +6 Проголосовать: не нравится

По первому вопросу -- надо брать модуль, который вместит произведение любой пары коэфициентов из перемножаемых многочленов.

По второму -- такое может быть только если и оригинальные многочлены обладали таким свойством (A*B=B).

По третьему -- на e-maxx это написано. Тебе надо найти простое число вида c*2^k+1 -- у него будет первообразный корень из единицы степени 2^k.

 

  • 13 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится
    Точнее, произведение любой пары * длину(степень) многочлена (все-таки перемножать многочлены с коэффициентами до 10^9 нельзя или почти всегда нельзя)
    • 13 лет назад, # ^ |
        Проголосовать: нравится 0 Проголосовать: не нравится
      Похоже на то - перемножил два числа размером 1000 с базой 100 - ок, размером 10000 с базой 100 - неправильно.
  • 13 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится
    По первому вопросу не совсем то. Мы переводим какиме-то наши полиномы в вектора DFT(A) и DFT(B). В них все числа  меньше модуля m. Пусть у нас в DFT(A)[0] получится (m - 1) и в DFT(B)[0] получился (m - 1) После скалярного перемножения в (DFT[A], DFT[B])[0] будет (m-1)^2, что больше m. Вопрос - как обратное преобразование отнесется к такому числу в этом векторе?
    • 13 лет назад, # ^ |
        Проголосовать: нравится +5 Проголосовать: не нравится
      Если мы умножаем полиномы в Zm[x], то и полином-результат будет в Zm[x] (а там вообще нет понятий "больше"/"меньше"). Поэтому при достаточно большом m из такого результата можно извлечь результат и в Z[x].
  • 13 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится
    Да, и получается, что квадрат модуля у нас должен вмещаться в long long, То есть уже база ограничена сверху корнем четвертой степени из 10^18?
    • 13 лет назад, # ^ |
        Проголосовать: нравится +6 Проголосовать: не нравится
      Поэтому все в таких случаях и пишут FFT с double'ами - оно работает быстрее и основание системы счисления можно брать не очень маленьким (но 109 всё равно не стоит).
      • 13 лет назад, # ^ |
          Проголосовать: нравится 0 Проголосовать: не нравится
        ну вот :) А я как раз из соображений что будет без потерь точности с большой базой с целыми и попробовал писать. Ну да ладно, узнал немного по теории чисел зато.
13 лет назад, # |
  Проголосовать: нравится +6 Проголосовать: не нравится
На счёт третьего. Вот список простых чисел вида 225k+1 меньших 231 с первообразными корнями.
p=167772161  g=3
p=469762049  g=3 (это даже 226k+1)
p=1107296257  g=10
p=1711276033  g=29
p=1811939329  g=13 (это даже 226k+1)
p=2013265921  g=31 (это даже 227k+1)
p=2113929217  g=5

Можно делать еще такой трюк, чтобы избежать переполнения. Посчитать всё по двум модулям. А потом по КТО восстановить. Тогда можно брать хоть 106 за основание системы счисления. Из-за двух модулей все будет в два раза дольше, но за счёт большого основания получается все таки выигрыш во времени.
  • 13 лет назад, # ^ |
      Проголосовать: нравится 0 Проголосовать: не нравится
    спасибо и за модули и за совет про теорему об остатках. попробую посмотреть.
13 лет назад, # |
  Проголосовать: нравится -7 Проголосовать: не нравится
Про то где почитать. В Кормене очень хорошо написано про БПФ. Правда про БПФ с модулярной арифметикой только одно упражнение самое последнее, но его стоит порешать.
13 лет назад, # |
  Проголосовать: нравится +5 Проголосовать: не нравится

Да, с информацией по целочисленному ФФТ напряжёнка, но когда-то я находил это в одной книге. Сейчас попытался снова откопать - не смог :(

Но, видимо, ответы на все вопросы и так уже получены

13 лет назад, # |
  Проголосовать: нравится 0 Проголосовать: не нравится
Вот статья про FFT от yeputons. Тоже вроде неплохая.
13 лет назад, # |
Rev. 2   Проголосовать: нравится +4 Проголосовать: не нравится

Для интересующихся приведу ссылку на статью, где описан алгоритм быстрого умножения матриц, который использует практически такие же идеи, как и БПФ.

http://arxiv.org/abs/math/0307321