Алгоритм Шёнхаге - Штрассена

[Euler-Maskerony]

Здравствуйте. Я пытаюсь написать свой шаблонный класс для чисел с размером большим, чем размер слова процессора. Для больших чисел я хочу реализовать алгоритм Шёнхаге - Штрассена, но я не могу полностью понять, как он работает. Дело в том, что он использует FFT в кольце по модулю $2^{2^N}+1$. Как я понял, для этого подбирается достаточно большое $N$, чтобы все операции внутри FFT проводились так, как будто бы мы работаем с обычными целыми числами. К сожалению, я совсем не могу понять что мы делаем после получения DFT обоих множителей. Дело в том, что, как я понял, внутри FFT нами используются только битовые сдвиги и сложения. После же нам необходимо перемножить полученные результаты $X$ и $Y$, чтобы получить, можно сказать, результат умножения $M$, т.е. $M = IFFT(X \cdot Y)$. Но этого я совсем не понимаю, ведь мы не умеем умножать большие числа, а числа внутри $X,Y$ такие же большие как и входные (может даже больше). Как же их умножить? Помогите, пожалуйста

[mihaild]

Euler-Maskerony в сообщении #1602453 писал(а):

Как я понял, для этого подбирается достаточно большое $N$, чтобы все операции внутри FFT проводились так, как будто бы мы работаем с обычными целыми числами

Нам же там обратное по модулю считать. Ну и в любом случае - нет, мы работаем в поле по модулю, это важно.

Euler-Maskerony в сообщении #1602453 писал(а):

а числа внутри $X,Y$ такие же большие как и входные (может даже больше).

Должны получаться меньше.

Дайте, пожалуйста, ссылку на описание алгоритма, которое Вы берете. Там должно что-то говориться о подборе параметров, чтобы получалась наибольшая экономия.

[Euler-Maskerony]

mihaild в сообщении #1602466 писал(а):

Нам же там обратное по модулю считать. Ну и в любом случае - нет, мы работаем в поле по модулю, это важно.

Да, это конечно. Я плохо выразился значит. Я имел в виду, что $N$ должно быть достаточно большим, чтобы результат совпадал с тем, который мы имеем в поле целых чисел.

mihaild в сообщении #1602466 писал(а):

Должны получаться меньше.
Дайте, пожалуйста, ссылку на описание алгоритма, которое Вы берете. Там должно что-то говориться о подборе параметров, чтобы получалась наибольшая экономия.

Самое подробное и понятное описание я нашел только здесь, но оно для меня тоже не совсем понятно: https://en.wikipedia.org/w/index.php?title=Sch%C3%B6nhage%E2%80%93Strassen_algorithm&oldid=1157584160

Дело в том, что здесь какая-то путанница с разными модулями и я совсем не понимаю этого. Как я вижу (опуская оптимизации):
1. Получаем на вход два числа $2^m$ бит каждое.
2. Подаем их в FFT в виде отдельных цифр (на сколько цифр разбивать?) (коэффициенты многочлена). FFT в качестве "roots of unity" берет степени "двойки". Все вычисления проводятся по модулю (?).
3. FFT выдает значения многочленов в корнях. Ответ получается по модулю (?).
4. Перемножаем ответы и подаем в IFFT. Получаем свертки. Переносим остатки и всё.

На месте "(?)" я не знаю что должно стоять. Совсем каша в голове. Несколько дней уже вожусь с этим. Помогите понять, пожалуйста.

[mihaild]

Euler-Maskerony в сообщении #1602468 писал(а):

Я имел в виду, что $N$ должно быть достаточно большим, чтобы результат совпадал с тем, который мы имеем в поле целых чисел

В смысле конечный результат умножения? (ну и целые числа не поле, но неважно)

Да, в википедии хорошее описание (которое кто-то испортил, надо будет внимательно посмотреть и, вероятно, откатить).
В Ваших обозначениях, $N = 2^m$.
2. Разбиваем на цифры примерно по $2^{m / 2}$ бит каждая (цифр соответственно тоже $2^{m / 2}$). Выбираем $n = 3 \cdot 2^{m / 2} = 2^{m / 2 + \log_2 3}$, и вычисления ведем по модулю $2^n + 1$ (в котором почти в 2 раза меньше бит чем в исходном модуле $2^N$).
3. Ответ получается по модулю $2^n + 1$.
В результате у нас получается $2^{m / 2}$ "псефдоцифр" (компонент свёртки до переноса), которые теперь могут быть больше $2^{m / 2}$, но они не меньше, чем $2^n + 1$, поэтому мы посчитали их точно.

[Euler-Maskerony]

mihaild в сообщении #1602471 писал(а):

Euler-Maskerony в сообщении #1602468 писал(а):

Я имел в виду, что $N$ должно быть достаточно большим, чтобы результат совпадал с тем, который мы имеем в поле целых чисел

В смысле конечный результат умножения? (ну и целые числа не поле, но неважно)

Да, в википедии хорошее описание (которое кто-то испортил, надо будет внимательно посмотреть и, вероятно, откатить).
В Ваших обозначениях, $N = 2^m$.
2. Разбиваем на цифры примерно по $2^{m / 2}$ бит каждая (цифр соответственно тоже $2^{m / 2}$). Выбираем $n = 3 \cdot 2^{m / 2} = 2^{m / 2 + \log_2 3}$, и вычисления ведем по модулю $2^n + 1$ (в котором почти в 2 раза меньше бит чем в исходном модуле $2^N$).
3. Ответ получается по модулю $2^n + 1$.
В результате у нас получается $2^{m / 2}$ "псефдоцифр" (компонент свёртки до переноса), которые теперь могут быть больше $2^{m / 2}$, но они не меньше, чем $2^n + 1$, поэтому мы посчитали их точно.

Спасибо вам за пояснения! Только я так и не понимаю зачем на $N$, если во всех операциях мы используем $n$? Ещё не понимаю почему $N = 2^m$, ведь мы должны выбрать его таким, чтобы результат умножения влез в кольцо, т.е. $N = 2^{m+1}$ и модуль равен $2^{2^{m+1}}+1$, т.е. размер кольца вдвое больше размеров входных чисел.
На википедии используется $n=\frac{2N}{2^{m / 2}} + \frac{m}{2}$ и пишется, что результат FFT влезет в такой модуль. И вы тоже пишите похожее. Но я не понимаю какая нам выгода. Ведь в качестве корней в FFT мы используем $r = 2^{2n/2^{m / 2}} = 2^{\frac{4N}{2^m} + \frac{m}{2^{m/2+1}}} = [ N = 2^{m+1} ] = 2^{8 + \frac{m}{2^{m/2+1}}}$. Таким образом мы используем $2^{m / 2}$ корней вида $2^{(8 + \frac{m}{2^{m/2+1}})j}; j \in \{1 .. 2^{m/2}\}$. Значит найдется $j$, которое сделает корень больше $2^m$, а значит больше самого входного числа (потому что внутри FFT в случае, когда корень равен $2^m$, считается, по сути, входное число как оно есть). Таким образом, получающиеся результаты FFT больше самих чисел, а значит нет смысла в алгоритме в целом.
Скажите, пожалуйста, где я ошибаюсь?
Ещё в статье зачем-то исходные числа умножают на какой-то вектор $A$. Я читал где-то, что так можно делать, и что это называется взвешенным FT, но я не понимаю для чего это здесь. Это какая-то оптимизация или что-то критически важное?

[mihaild]

Euler-Maskerony в сообщении #1602534 писал(а):

Только я так и не понимаю зачем на $N$, если во всех операциях мы используем $n$?

Это просто часть постановки. Ответ нужен по модулю $N$, но мы можем свести задачу к модулю $n$.

Euler-Maskerony в сообщении #1602534 писал(а):

Ещё не понимаю почему $N = 2^m$, ведь мы должны выбрать его таким, чтобы результат умножения влез в кольцо

Потому что опечатка (у меня). Для простоты предлагаю считать, что это $x$ и $y$ были изначально из $2^{m / 2}$ цифр.

Euler-Maskerony в сообщении #1602534 писал(а):

И вы тоже пишите похожее

Да я собственно пересчитал из википедии:)

Euler-Maskerony в сообщении #1602534 писал(а):

Таким образом, получающиеся результаты FFT больше самих чисел, а значит нет смысла в алгоритме в целом

А вот это интересный и важный момент. FFT мы делаем в $\mathbb Z_{2^n + 1}$, и результат получается в нём же. Но поскольку мы заранее знаем, что ответ (сумма произведений цифр) не больше $2^n$, то ответ точно восстанавливается из ответа по модулю $2^n + 1$, хотя в промежуточных вычислениях мы существенно использовали модуль.
Более простой аналог: если Вы знаете сумму $a_1 + a_2 + a_3$ по модулю $5$, то Вы не можете восстановить значение этой суммы в $\mathbb Z$. Но если Вы знаете эту сумму, а еще знаете, что она попадает в диапазон от $0$ до $4$, то можете. Более того, хотя $a_1 + a_2$ могут в этот диапазон не попадать, всё равно можно вести все вычисления по модулю $5$, и в конце получить правильный ответ.

Euler-Maskerony в сообщении #1602534 писал(а):

Ещё в статье зачем-то исходные числа умножают на какой-то вектор $A$. Я читал где-то, что так можно делать, и что это называется взвешенным FT, но я не понимаю для чего это здесь. Это какая-то оптимизация или что-то критически важное?

Это важно для оптимизации (мы на самом деле считаем не совсем свёртку), но не критично.

[Euler-Maskerony]

mihaild в сообщении #1602540 писал(а):

А вот это интересный и важный момент. FFT мы делаем в $\mathbb Z_{2^n + 1}$, и результат получается в нём же. Но поскольку мы заранее знаем, что ответ (сумма произведений цифр) не больше $2^n$, то ответ точно восстанавливается из ответа по модулю $2^n + 1$, хотя в промежуточных вычислениях мы существенно использовали модуль.
Более простой аналог: если Вы знаете сумму $a_1 + a_2 + a_3$ по модулю $5$, то Вы не можете восстановить значение этой суммы в $\mathbb Z$. Но если Вы знаете эту сумму, а еще знаете, что она попадает в диапазон от $0$ до $4$, то можете. Более того, хотя $a_1 + a_2$ могут в этот диапазон не попадать, всё равно можно вести все вычисления по модулю $5$, и в конце получить правильный ответ.

Ого, спасибо большое. Теперь, кажется, стало понятнее.

mihaild в сообщении #1602540 писал(а):

Это важно для оптимизации (мы на самом деле считаем не совсем свёртку), но не критично.

Хорошо, понятно. Не могли бы вы пояснить, пожалуйста, еще кое что. В википедии пишут, что $CyclicConvolution(X, Y) = IDFT(DFT(X) \cdot DFT(Y))$, но я не понимаю почему. Мы же умножаем, по сути, два многочлена, а значит в результате должна получиться обычная свертка, а не циклическая. О чем же здесь идет речь? А еще не понимаю в чем преимущество negacyclic свертки, ведь там даже говориться, что мы можем получать отрицательные числа, поэтому нужно дополнительно этим возиться.

[mihaild]

Euler-Maskerony в сообщении #1602546 писал(а):

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

Нам нужна обычная свертка. Но получается циклическая, потому что мы берём дискретное преобразование Фурье, а не обычное. Но, как написано в википедии, получающееся значение сравнимо с нужным по модулю $2^n - 1$.

Euler-Maskerony в сообщении #1602546 писал(а):

А еще не понимаю в чем преимущество negacyclic свертки, ведь там даже говориться, что мы можем получать отрицательные числа, поэтому нужно дополнительно этим возиться.

Там получается модуль $2^n + 1$, в котором, например, легко находится примитивный корень, являющийся степенью двойки.

[Geen]

Euler-Maskerony в сообщении #1602453 писал(а):

Для больших чисел я хочу реализовать алгоритм Шёнхаге - Штрассена

Скажите, а это для практических целей или для тренировки?
А то ведь, если я правильно понимаю, 1) этот алгоритм быстрее на числах размером более ста тысяч бит; 2) алгоритм "сводится" к умножению чисел вдвое меньшего размера - то есть умножать большие числа "обычным" способом всё равно необходимо.

[mihaild]

Geen в сообщении #1602611 писал(а):

этот алгоритм быстрее на числах размером более ста тысяч бит

GMP на skylake использует FFT начиная с 401408 бит https://gmplib.org/repo/gmp/file/tip/mp ... ram.h#l142.

Geen в сообщении #1602611 писал(а):

алгоритм "сводится" к умножению чисел вдвое меньшего размера - то есть умножать большие числа "обычным" способом всё равно необходимо

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

[Geen]

mihaild в сообщении #1602621 писал(а):

Практически конечно это не очень хорошо работать будет.

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

[Euler-Maskerony]

mihaild в сообщении #1602608 писал(а):

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

Странное что-то. А может быть это из-за того, что мы не увеличиваем $X$ и $Y$ вдвое (в википедии такая картинка)? Ведь если увеличить, то получится обычное умножение многочленов? Или так здесь не будет работать? Должно же вроде работать.

Geen в сообщении #1602611 писал(а):

Скажите, а это для практических целей или для тренировки?

Сложно сказать. Я хочу сделать хорошую библиотеку которая добавляет uint128_t, uint256_t,...,uintN_t типы в C++, чтобы использовать в своих программах.

Geen в сообщении #1602611 писал(а):

1) этот алгоритм быстрее на числах размером более ста тысяч бит; 2) алгоритм "сводится" к умножению чисел вдвое меньшего размера - то есть умножать большие числа "обычным" способом всё равно необходимо.

Да, поэтому для больших чисел будет использоваться такой алгоритм, а когда рекурсия дойдет до маленьких, то будет использоваться какой-нибудь Карацуба. По крайней мере такой план.

Geen в сообщении #1602625 писал(а):

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

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

[mihaild]

Euler-Maskerony в сообщении #1602777 писал(а):

А может быть это из-за того, что мы не увеличиваем $X$ и $Y$ вдвое (в википедии такая картинка)?

Да, поэтому.
Я нашел еще хорошую статью про этот алгоритм https://tonjanee.home.xs4all.nl/SSAdescription.pdf.

[Euler-Maskerony]

mihaild в сообщении #1602788 писал(а):

Я нашел еще хорошую статью про этот алгоритм https://tonjanee.home.xs4all.nl/SSAdescription.pdf.

Спасибо! Почитаю.

mihaild в сообщении #1602788 писал(а):

Да, поэтому.

А что если мы увеличим входы? Тогда же не придется с циклическими свертками возиться?

[mihaild]

Euler-Maskerony в сообщении #1602807 писал(а):

А что если мы увеличим входы? Тогда же не придется с циклическими свертками возиться?

Да, не придется. Но по итогу константа получится еще хуже.