3

下午好！

我正在尝试基于我已经拥有的朴素递归 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;

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

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

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

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

非常感谢您！

4

3 回答 3

3

我最近也想实现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. 重要价值观


    Wn = r ^ L mod p是NTT
    Wn = r ^ (p-1-L) mod p的基础 是INTT的基础 是INTT
    Rn = n ^ (p-2) mod p的缩放乘法常数是素数，并且是NTT的x[i]或INTT的X[i]的最大值，并且除法必须组合，所以 如果或必须组合， 如果是子结果最大值取决于计算类型和类型。对于单个(I)NTT，它是但对于两个大小的向量的卷积，它是等等。有关它的更多信息，请参阅在有限域上实现 FFT 。 ~(1/n)
    pp mod n == 1p>max'
    max
    r = <1,p)
    L = <1,p)p-1
    r,Lr^(L*i) mod p == 1i=0i=n
    r,Lr^(L*i) mod p != 10 < i < n
    max'nmax' = n*maxnmax' = n*max*max

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

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

这是我找到参数的C++r,L,p代码（需要不包括在内的模运算，您可以将其替换为 (a+b)%c,(ab)%c,(a*b)%c,... 但在这种情况下要注意溢出，特别是modpow)modmul代码尚未优化，但有一些方法可以显着加快速度。素数表也相当有限，因此要么使用SoE 或任何其他算法来获得素数，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在相同位大小的浮点上的速度相同。它的链接在这里：

于 2013-08-31T11:20:27.037 回答
1

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

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

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

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

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

您不需要将循环从 n / 2 扩展到 n。

所以，是的，一些消息来源说只是为模算术方法添加不同的欧米茄定义，正在扫除许多细节。

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

参考

  • Baktir 和 Sunar - 使用快速傅里叶变换在费马场中实现高效的多项式乘法 (2006)
于 2017-11-21T21:01:36.637 回答
0

你必须确保统一的根源确实存在。在R中，只有 2 个单位根：1 和 -1，因为只有对它们 x^n=1 才能成立。

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

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

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 为模。

于 2012-12-12T10:55:14.683 回答