fft - 从复数 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复数根的前半部分:

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 单位根？

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

提前非常感谢您!

Answer 1

我最近想实现NTT来实现快速乘法，而不是DFFT。读了很多令人困惑的东西，到处都有不同的字母，没有简单的解决方案，而且我的有限域知识也很生疏，但今天我终于做对了(经过两天的尝试和用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. 重要值(value)观

   
   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 if 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,(a-b)%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 )

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

[编辑1:快速 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]

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

• Modular arithmetics and NTT (finite field DFT) optimizations