从复数 FFT 到有限场 FFT 的转换

Question

下午好！

我正在尝试基于我已经拥有的朴素递归 FFT 实现来开发 NTT 算法。

考虑以下代码（coefficients' 的长度，设为 m，是 2 的精确幂）：

/// <summary>
/// Calculates the result of the recursive Number Theoretic Transform.
/// </summary>
/// <param name="coefficients"></param>
/// <returns></returns>
private static BigInteger[] Recursive_NTT_Skeleton(
    IList<BigInteger> coefficients, 
    IList<BigInteger> rootsOfUnity, 
    int step, 
    int offset)
{
    // Calculate the length of vectors at the current step of recursion.
    // -
    int n = coefficients.Count / step - offset / step;

    if (n == 1)
    {
        return new BigInteger[] { coefficients[offset] };
    }

    BigInteger[] results = new BigInteger[n];

    IList<BigInteger> resultEvens = 
        Recursive_NTT_Skeleton(coefficients, rootsOfUnity, step * 2, offset);

    IList<BigInteger> resultOdds = 
        Recursive_NTT_Skeleton(coefficients, rootsOfUnity, step * 2, offset + step);

    for (int k = 0; k < n / 2; k++)
    {
        BigInteger bfly = (rootsOfUnity[k * step] * resultOdds[k]) % NTT_MODULUS;

        results[k]          = (resultEvens[k] + bfly) % NTT_MODULUS;
        results[k + n / 2]  = (resultEvens[k] - bfly) % NTT_MODULUS;
    }

    return results;
}

它适用于复杂的 FFT（将 BigInteger 替换为复杂的数字类型（我有自己的））。即使我适当地改变了寻找统一原根的过程，它在这里也不起作用。

据推测，问题是这样的：rootsOfUnity 传递的参数最初只包含 m-th 复单位根的前半部分，按此顺序：

omega^0 = 1, omega^1, omega^2, ..., omega^(n/2)

就够了，因为就这三行代码：

BigInteger bfly = (rootsOfUnity[k * step] * resultOdds[k]) % NTT_MODULUS;        

results[k]          = (resultEvens[k] + bfly) % NTT_MODULUS;
results[k + n / 2]  = (resultEvens[k] - bfly) % NTT_MODULUS;

我最初利用了这样一个事实，即在任何递归级别（对于任何n 和i），统一的复根-omega^(i) = omega^(i + n/2)。

但是，该属性显然不适用于有限域。但是是否有任何类似物可以让我仍然只计算根的前半部分？

或者我应该将循环从n/2 扩展到n 并预先计算所有m-th 统一的根？

也许这段代码还有其他问题？..

非常感谢您！

Answer 1

我最近也想实现 NTT 来实现快速乘法，而不是 DFFT。阅读了很多令人困惑的东西，到处都是不同的字母，没有简单的解决方案，而且我的有限领域知识也生疏了，但今天我终于弄对了（经过 2 天的尝试和模拟 DFT系数）所以这是我对 NTT 的见解：

1. 计算

   X(i) = sum(j=0..n-1) of ( Wn^(i*j)*x(i) );
   

   其中X[] 是NTT 转换的x[]，大小为n，其中Wn 是NTT 基。所有计算都基于整数模运算mod p 任何地方都没有复数。

2. 重要价值

   
   Wn = r ^ L mod p 是 NTT 的基础
   Wn = r ^ (p-1-L) mod p 是 INTT 的基础
   Rn = n ^ (p-2) mod p 正在为 INTT 缩放乘法常数 ~(1/n)
p 是 p mod n == 1 和 p>max' 的质数
   max 是 NTT 的最大值 x[i] 或 INTT 的 X[i] >
   r = <1,p)
L = <1,p) 也分p-1
r,L 必须合并，所以 r^(L*i) mod p == 1 如果 i=0 或 i=n
r,L 必须合并，所以r^(L*i) mod p != 1 如果0 < i < n
max' 是子结果最大值，取决于n 和计算类型。对于单个 (I)NTT，它是 max' = n*max，但对于两个 n 大小的向量的卷积，它是 max' = n*max*max 等。有关它的更多信息，请参阅 Implementing FFT over finite fields。

3. r,L,p的功能组合对于不同的n是不同的

   这很重要，您必须在每个 NTT 层之前重新计算或从表中选择参数（n 始终是前一个递归的一半）。

这是我的 C++ 代码，它找到了r,L,p 参数（需要不包含的模运算，您可以将其替换为 (a+b)%c,(ab)%c, (a*b)%c,... 但在这种情况下要小心溢出，特别是 modpow 和 modmul) 代码尚未优化，但有一些方法可以显着加快速度。主表也相当有限，因此请使用SoE or any other algo to obtain primes up to max' 以安全工作。

DWORD _arithmetics_primes[]=
    {
    2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,
    179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,
    419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,
    661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
    947,953,967,971,977,983,991,997,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,
    0}; // end of table is 0, the more primes are there the bigger numbers and n can be used
// compute NTT consts W=r^L%p for n
int i,j,k,n=16;
long w,W,iW,p,r,L,l,e;
long max=81*n;  // edit1: max num for NTT for my multiplication purposses
for (e=1,j=0;e;j++)             // find prime p that p%n=1 AND p>max ... 9*9=81
    {
    p=_arithmetics_primes[j];
    if (!p) break;
    if ((p>max)&&(p%n==1))
     for (r=2;r<p;r++)  // check all r
        {
        for (l=1;l<p;l++)// all l that divide p-1
            {
            L=(p-1);
            if (L%l!=0) continue;
            L/=l;
            W=modpow(r,L,p);
            e=0;
            for (w=1,i=0;i<=n;i++,w=modmul(w,W,p))
                {
                if ((i==0)      &&(w!=1)) { e=1; break; }
                if ((i==n)      &&(w!=1)) { e=1; break; }
                if ((i>0)&&(i<n)&&(w==1)) { e=1; break; }
                }
            if (!e) break;
            }
        if (!e) break;
        }
    }
if (e) { error; }           // error no combination r,l,p for n found
 W=modpow(r,    L,p);   // Wn for NTT
iW=modpow(r,p-1-L,p);   // Wn for INTT

这是我的慢 NTT 和 INTT 实现（我还没有快速 NTT、INTT），它们都成功地使用 Schönhage-Strassen 乘法进行了测试。

//---------------------------------------------------------------------------
void NTT(long *dst,long *src,long n,long m,long w)
    {
    long i,j,wj,wi,a,n2=n>>1;
    for (wj=1,j=0;j<n;j++)
        {
        a=0;
        for (wi=1,i=0;i<n;i++)
            {
            a=modadd(a,modmul(wi,src[i],m),m);
            wi=modmul(wi,wj,m);
            }
        dst[j]=a;
        wj=modmul(wj,w,m);
        }
    }
//---------------------------------------------------------------------------
void INTT(long *dst,long *src,long n,long m,long w)
    {
    long i,j,wi=1,wj=1,rN,a,n2=n>>1;
    rN=modpow(n,m-2,m);
    for (wj=1,j=0;j<n;j++)
        {
        a=0;
        for (wi=1,i=0;i<n;i++)
            {
            a=modadd(a,modmul(wi,src[i],m),m);
            wi=modmul(wi,wj,m);
            }
        dst[j]=modmul(a,rN,m);
        wj=modmul(wj,w,m);
        }
    }
//---------------------------------------------------------------------------

dst 是目标数组
src 是源数组
n 是数组大小
m 是模数 (p)
w 是基础 (Wn)

希望这对某人有所帮助。如果我忘记了什么，请写...

[edit1：快速 NTT/INTT]

我终于设法让 NTT/INTT 快速工作。比普通的 FFT 稍微复杂一点：

//---------------------------------------------------------------------------
void _NFTT(long *dst,long *src,long n,long m,long w)
    {
    if (n<=1) { if (n==1) dst[0]=src[0]; return; }
    long i,j,a0,a1,n2=n>>1,w2=modmul(w,w,m);
    // reorder even,odd
    for (i=0,j=0;i<n2;i++,j+=2) dst[i]=src[j];
    for (    j=1;i<n ;i++,j+=2) dst[i]=src[j];
    // recursion
    _NFTT(src   ,dst   ,n2,m,w2);   // even
    _NFTT(src+n2,dst+n2,n2,m,w2);   // odd
    // restore results
    for (w2=1,i=0,j=n2;i<n2;i++,j++,w2=modmul(w2,w,m))
        {
        a0=src[i];
        a1=modmul(src[j],w2,m);
        dst[i]=modadd(a0,a1,m);
        dst[j]=modsub(a0,a1,m);
        }
    }
//---------------------------------------------------------------------------
void _INFTT(long *dst,long *src,long n,long m,long w)
    {
    long i,rN;
    rN=modpow(n,m-2,m);
    _NFTT(dst,src,n,m,w);
    for (i=0;i<n;i++) dst[i]=modmul(dst[i],rN,m);
    }
//---------------------------------------------------------------------------

[编辑3]

我已经优化了我的代码（比上面的代码快 3 倍），但我仍然对它不满意，所以我开始用它提出新问题。在那里，我进一步优化了我的代码（比上面的代码快约 40 倍），因此它的速度几乎与 FFT 在相同位大小的浮点上相同。链接在这里：

• Modular arithmetics and NTT (finite field DFT) optimizations

Answer 2

要将 Cooley-Tukey（复数）FFT 转换为模算术方法，即 NTT，您必须替换 omega 的复数定义。对于纯粹递归的方法，您还需要根据当前信号大小重新计算每个级别的 omega。这是可能的，因为分钟。当我们在调用树中向下移动时，合适的模数会减少，因此用于根的模数适用于较低层。此外，由于我们使用相同的模数，当我们向下移动调用树时，可能会使用相同的生成器。此外，对于逆变换，您应该采取额外的步骤来重新计算欧米茄 a 并改为使用欧米茄：b = a ^ -1 （通过使用反模运算）。具体来说，b = invMod(a, N) s.t. b * a == 1 (mod N)，其中 N 是选择的素数模数。

利用周期性重写涉及 omega 的表达式仍然适用于模算术领域。您还需要找到一种方法来确定问题的模数（素数）和有效的生成器。

我们注意到您的代码有效，尽管它不是 MWE。我们使用常识对其进行了扩展，并为多项式乘法应用得到了正确的结果。您只需要提供正确的 omega 值，提高到特定的幂即可。

虽然您的代码可以正常工作，但与许多其他来源一样，您将每个级别的间距加倍。但是，这不会导致递归那么干净；这与基于当前信号大小重新计算 omega 相同，因为 omega 定义的功率与信号大小成反比。重申一下：将信号大小减半就像平方欧米茄，这就像给欧米茄提供双倍的功率（这是一个将间距加倍所做的事情）。处理重新计算 omega 的方法的好处在于，每个子问题本身都更加干净完整。

有一篇论文展示了模块化方法的一些数学；这是 Baktir 和 Sunar 于 2006 年发表的论文。请参阅本文末尾的论文。

您不需要将循环从 n / 2 延长到 n。

所以，是的，一些消息来源说只是为模算术方法添加不同的 omega 定义，但他们掩盖了许多细节。

另一个问题是，如果我们在执行卷积时要避免结果时域信号溢出，那么必须承认信号大小必须足够大，这一点很重要。此外，即使幂非常大，找到一些快速的模数求幂实现也可能很有用。

参考文献

• Baktir 和 Sunar - 使用快速傅里叶变换在费马场中实现高效的多项式乘法 (2006)

Answer 3

您必须确保统一的根源确实存在。在 R 中只有 2 个单位根：1 和 -1，因为只有 x^n=1 对它们来说是真的。

在 C 中，您有无限多个单位根：w=exp(2*pi*i/N) 是原始的第 N 个单位根，所有 w^k 都为 0

现在你的问题是：你必须确保你正在使用的环提供相同的属性：足够的统一根。

Schönhage 和 Strassen (http://en.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Strassen_algorithm) 使用模 2^N+1 的整数。这个圆环有足够的统一根。 2^N == -1 是二阶单位根，2^(N/2) 是四阶单位根，依此类推。此外，这些单位根的优点是它们是 2 的幂，并且可以实现为二进制移位（之后进行模运算，归结为加/减）。

我认为 QuickMul (http://www.cs.nyu.edu/exact/doc/qmul.ps) 以 2^N-1 为模。