Ошибка реализации умножения Шонхаге-Штрассена

Question

Я пытаюсь реализовать алгоритм умножения Шонхаге-Штрассена, используя NTT, и столкнулся с проблемой, когда конечный результирующий вектор фактически не равен тому, что должен быть.

Для двух входных векторов a и b, каждая из которых состоит из N "цифр" K битов (конечных N/2 записей каждого набора в 0), каждая, при заданном модуле M = 2^(2*K)+1, корне единицы w = N^(4*K-1) | w^N = 1 mod M, модульная инверсия этого значения wi | wi*w = 1 mod M, и u | u*N = 1 mod M, следующий код Python используется для (попытки) умножения этих векторов с использованием алгоритма Шонхаге-Штрассена:

#a and b are lists of length N, representing large integers
A = [ sum([ (a[i]*pow(w,i*j,M))%M for i in range(N)]) for j in range(N)] #NTT of a
B = [ sum([ (b[i]*pow(w,i*j,M))%M for i in range(N)]) for j in range(N)] #NTT of b
C = [ (A[i]*B[i])%M for i in range(N)] #A * B multiplied pointwise
c = [ sum([ (C[i]*pow(wi,i*j,M))%M for i in range(N)]) for j in range(N)] #intermediate step in INTT of C
ci = [ (i*u)%M for i in c] #INTT of C, should be product of a and b

В теории,если взять NTT из a и b, умножить по точкам, затем взять INTT результата, то я получу произведение, если я не ошибаюсь, и я проверил эти методы для NTT и INTT, чтобы подтвердить, что они являются обратнымидруг с другом. Однако конечный результирующий вектор ci, а не равен произведению a и b, является произведением, в котором каждый элемент берется по модулю M, что дает неверный результат для произведения.

Например, выполнение теста с N=K=8 и случайными векторами для a, b дает следующее:

M = 2^(2*8)+1 = 65537
w = 16, wi = 61441
u = 57345
a = [212, 251, 84, 186, 0, 0, 0, 0] (3126131668 as an integer)
b = [180, 27, 234, 225, 0, 0, 0, 0] (3790216116)
NTT(a) = [733, 66681, 147842, 92262, 130933, 107825, 114562, 127302]
NTT(b) = [666, 64598, 80332, 54468, 131236, 186644, 181708, 88232]
Pointwise product of above two lines mod M = [29419, 39913, 25015, 14993, 42695, 49488, 52438, 51319]
INTT of above line (i.e. result) = [38160, 50904, 5968, 11108, 15616, 62424, 41850, 0] (11848430946168040720)
Actual product of a x b = [38160, 50904, 71505, 142182, 81153, 62424, 41850, 0] (11848714628791561488)

В этом примере, и почти каждый раз, когда я его пробую,элементы фактического произведения и результат моего алгоритма одинаковы для нескольких элементов около начала и конца вектора, но к середине они отклоняются. Как я упоминал выше, элементы ci каждый равен элементам a*b по модулю M. Я, должно быть, неправильно понимаю что-то об этом алгоритме, хотя я не совсем уверен, что. Я где-то использую неправильный модуль?

Answer 1

Остерегайтесь теории чисел и NTT - не моя область знаний, поэтому прочитайте с предубеждением, но я успешно внедрил NTT в C ++ самостоятельно и использовал егодля умножения Бигнум (bigint, bigfloatingpoint, bigfixedpoint), так что вот некоторые из моих insigtst. Я настоятельно рекомендую вам сначала прочитать оба моих QA-запроса:

• Перевод кода с комплексного FFT на код конечного поля-FFT , который находит и проверяет первичный корень NTT единицыи связанные параметры
• Модульная арифметика и оптимизация NTT (конечное поле DFT) Мой 32-битный оптимизированный NTT

, так что вы можете сравнить свои результаты / код / ​​константы сшахта. Однако я усовершенствовал свой NTT, чтобы использовать одиночное простое кодированное простое число (самый большой предел единства, который соответствует 32-битному значению).

Теперь, что может быть не так с вашим кодом. Я не пишу код на python, но я не вижу код NTT в вашем вопросе. В любом случае из того, что я вижу:

1. проверьте свой корень или единство

   В своем вопросе вы упоминаете условие:

   w^N = 1 mod M
   

   но этого далеко не достаточно. См. Первую ссылку выше, где описаны все условия, которые должны быть выполнены (с кодом, который находит и проверяет его). Я не уверен, что ваши параметры соответствуют всем необходимым условиям, и вы просто забыли или пропустили, написали они или нет , поэтому проверьте его . IIRC Я тоже боролся с этими условиями, поскольку в то время, когда я кодировал это там, где в моем распоряжении было очень мало информации о NTT, и большинство из них были неполными или неправильными ...

2. Выне используют модульную арифметику !!!

   Я предполагаю, что ваше простое число равно M (в моей терминологии это p), поэтому все подрезультаты должны быть меньше, чем M, что явно не соответствует действительностив вашем примере:

   M = 65537
   NTT(a) = [733, 66681, 147842, 92262, 130933, 107825, 114562, 127302]
   NTT(b) = [666, 64598, 80332, 54468, 131236, 186644, 181708, 88232]
   

   , как вы можете видеть, действителен только первый элемент обоих NTT, все остальные больше M, что неверно !!!

3. Остерегайтесь переполнений

   Ваш M действительно очень мал ~16bit по сравнению с вашими входными значениями, которые выглядят ~8bit, которые могут очень быстро переполняться, что делает недействительными ваши NTT результаты тоже.

   Вот цитата из моей второй ссылки, которую я нашел сложным и эмпирическим способом:

   Чтобы избежать переполнений для больших наборов данных, ограничьте вводцифры в п / 4 бит. Где p - количество бит на элемент NTT, поэтому для этой 32-битной версии используйте максимальные (32 бит / 4 -> 8 бит) входные значения.

   , так что в вашем случае вы должны обработать 16/4 = 4bit чанковвместо 8 бит или используйте намного больший M, например, как у меня 0xC0000001, который равен ~32bit.

   . Это объясняет ваши наблюдения, что первые элементы продукта хороши, а затем не ... понимаете, если выумножьте 2 8-битных числа, которые вы получили, на 16 бит ... теперь вы понимаете, что вы делаете более рекурсивные добавления умноженных подрезультатов, так что в вашем случае очень скоро вырастет до 16 бит M прямо во второе значение ...

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