python - Schonhage-Strassen 乘法实现错误

Question

我正在尝试使用 NTT 实现 Schonhage-Strassen 乘法算法，但遇到了一个问题，即最终生成的向量实际上并不等于它应有的值。

对于两个输入向量 a 和 b，每个向量由 N 个“数字”组成，共 K 位(每个设置为 0 的最终 N/2 条目)，每个条目，给定一个模数 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 代码用于(尝试)使用 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

理论上，取 a 和 b 的 NTT，逐点相乘，然后取结果的 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 不是我的专业领域，所以带着偏见阅读，但我确实成功地用 C++ 实现了 NTT并将其用于 bignum 乘法(bigint、bigfloatingpoint、bigfixedpoint)所以这是我的一些见解。我强烈建议您先阅读我的两个相关 QA:

• Translation from Complex-FFT to Finite-Field-FFT 寻找和验证NTT素单位根及相关参数的代码
• Modular arithmetics and NTT (finite field DFT) optimizations mine 优化的 32 位 NTT

这样您就可以将您的结果/代码/常量与我的进行比较。然而，我改进了我的 NTT 以使用单个硬编码素数(适合 32 位值的最大单位根)。

现在您的代码可能出了什么问题。我没有用 python 编写代码，但我在您的问题中没有看到 NTT 代码。无论如何，从我所看到的:

1. 检查你的 root 或 unity

   在你的问题中你提到了条件:

   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 结果失效。

   这里引用我的第二个链接，我发现了困难和经验的方法:

   To avoid overflows for big datasets, limit input numbers to p/4 bits. Where p is number of bits per NTT element so for this 32 bit version use max (32 bit/4 -> 8 bit) input values.

   所以在你的情况下你应该处理 16/4 = 4bit block 而不是 8 位或者使用更大的 M 例如像我的 0xC0000001 即 ~32bit。

   这解释了您的观察结果，即乘积的第一个元素是好的，然后不是...意识到如果您将 2 个 8 位数字相乘，您将得到 16 位...现在意识到您正在对乘积的子结果进行更多的递归加法，因此它将得到超过 16 位 M 在你的情况下很快就在第二个值中......

总而言之，您没有使用模块化算术，素数太小和/或处理太大的数据 block 和可能也选择了错误的素数。