python - Schonhage-Strassen乘法实现错误

我正在尝试使用ntt实现schonhage-strassen乘法算法，并且遇到一个问题，即最终得到的向量实际上并不等于它应该是的向量。
对于两个输入向量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，使用以下python代码（尝试）使用Schonhage Strassen算法将这些向量相乘：

#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

理论上，取〈cc〉和〈cc〉的ntt，逐点相乘，然后取结果的intt，如果我没有弄错的话，就应该得到乘积，我对这些方法进行了ntt和intt的测试，以确认它们是彼此相反的。然而，最终得到的向量u | u*N = 1 mod M不是等于a和b的乘积，而是每个元素取模ci的乘积，从而给出了不正确的结果。
例如，使用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)

在这个例子中，几乎每次我尝试，实际乘积的元素和我的算法的结果对于向量的开始和结束附近的几个元素都是一样的，但是它们朝中间偏移。如上所述，M的元素都等于N=K=8模a, b的元素我一定是对这个算法有点误解，不过我不完全确定是什么我在某个地方用错模数了吗？

最佳答案

谨防数量理论和NTT不是我的专长领域，所以用偏见来阅读，但是我自己成功地在C++中实现了NTT，并用它来做Bigimm乘法（bigint，ccc>，cc>），所以这里有我的一些想法。我强烈建议您先阅读我的两个相关QA：
Translation from Complex-FFT to Finite-Field-FFT查找并验证单位的ntt素根和相关参数的代码
Modular arithmetics and NTT (finite field DFT) optimizations矿井优化32位NTT
所以你可以比较你的结果/代码/常量和我的然而，我将NTT进化为使用单硬编码素数（符合32位值的最大统一根）。
现在你的代码有什么问题。我没有用python编写代码，但在您的问题中没有看到ntt代码。不管怎样，从我看来：
检查你的根或团结
在你的问题中你提到了条件：

w^N = 1 mod M

但这远远不够请参阅上面的第一个链接，它描述了必须满足的所有条件（使用找到并检查它的代码）。我不确定您的参数是否符合所有需要的条件，您只是忘记或错过了写这些或不检查它IIRC我也在这些条件下挣扎，因为当时我把这个编码在那里，只有很少的NTT信息可供我使用，其中大部分是不完整或错误的。。。
你没有使用模块化算法！！！
我假设你的素数是bigfloatingpoint（用我的术语来说是bigfixedpoint），所以所有的子结果必须小于M，这在你的例子中显然是不正确的：

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

如您所见，只有两个ntt的第一个元素是有效的，所有其他元素都大于p这是错误的！!!
当心溢流
与输入值相比，M非常小，输入值看起来也会溢出，从而使ntt结果很快失效。
这里引用我的第二个链接，我发现了一个困难的和经验性的方法：
为避免大数据集溢出，请将输入数限制为P/4位。
其中p是每个ntt元素的位数，因此对于这个32位版本
使用最大（32位/4->8位）输入值。
因此，在您的情况下，您应该处理M块，而不是8位块，或者使用更大的M块，例如我的~16bit块。
这解释了你的观察，产品的第一要素是好的，然后不是…意识到如果你把28位数字相乘，你得到16位…现在意识到你正在对乘法的子结果做更多的递归加法，所以在你的情况下，在第二个值中，它很快就会超过16位。。。
总而言之，你没有使用模块化的算术，有太小的素数和/或处理太大的数据块，也可能有错误的素数选择。

关于python - Schonhage-Strassen乘法实现错误，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/58289336/