从复数FFT到有限域FFT的翻译

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次复杂单位根: 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;

我最初利用了这个事实:在任何递归级别(对于任何 ni),复数单位根 -omega^(i) = omega^(i + n/2)

然而,这种特性显然在有限域中不成立。但是否存在任何类似的模拟它的方法,仍然允许我只计算前半部分的根?

或者,我应该将周期从 n/2 扩展到 n 并预先计算所有的 m 次单位根?

也许这段代码还存在其他问题?..

非常感谢您提前的帮助!


查看此处 http://www.jjj.de/fxt/fxtbook.pdf 和此处 http://www.physicsforums.com/showthread.php?t=140367 - Spektre
我刚刚编辑了我的回答,我的快速 NTT/INTT 现在已经包含在内,以防你仍然需要它。 - Spektre
有没有人有一个提出递归解决方案的来源?此外,-omega^(i) = omega^(i + n/2) 在有限域中确实成立。 - z.karl
3个回答

3

我最近想要使用NTT来进行快速乘法,而不是使用DFFT。但是读了很多混淆的东西,到处都是不同的字母,没有简单的解决方案,而且我的有限域知识有点生疏。但是今天我终于搞定了(经过两天的尝试和与DFT系数的类比),因此在这里分享一下我对NTT的见解:

  1. Computation

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

    where X[] is NTT transformed x[] of size n where Wn is the NTT basis. All computations are on integer modulo arithmetics mod p no complex numbers anywhere.

  2. Important values


    Wn = r ^ L mod p is basis for NTT
    Wn = r ^ (p-1-L) mod p is basis for INTT
    Rn = n ^ (p-2) mod p is scaling multiplicative constant for INTT ~(1/n)
    p is prime that p mod n == 1 and p>max'
    max is max value of x[i] for NTT or X[i] for INTT
    r = <1,p)
    L = <1,p) and also divides p-1
    r,L must be combined so r^(L*i) mod p == 1 if i=0 or i=n
    r,L must be combined so r^(L*i) mod p != 1 if 0 < i < n
    max' is the sub-result max value and depends on n and type of computation. For single (I)NTT it is max' = n*max but for convolution of two n sized vectors it is max' = n*max*max etc. See Implementing FFT over finite fields for more info about it.

  3. functional combination of r,L,p is different for different n

    this is important, you have to recompute or select parameters from table before each NTT layer (n is always half of the previous recursion).

以下是我找到的用于查找参数 r,L,p 的C++代码(需要进行模算术运算,但未包含在代码中,您可以使用(a+b)%c,(a-b)%c,(a*b)%c等代替,但在这种情况下,请注意特别处理modpowmodmul的溢出问题)。代码尚未优化,但有许多可加快速度的方法。此外,素数表相当有限,因此请使用 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)

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

[编辑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]

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


我注意到原帖作者想要一个递归版本的FFT,而不是你在这里呈现的迭代版本。对于迭代版本,我表示赞赏,但我真正感兴趣的是是否可以使用简单FFT实现中包含的递归结构执行NTT?你介意研究一下吗?我已经尝试过了,似乎NTT仅对索引1和1 + n/2成功,而其他所有索引都失败了。 - z.karl
@z.karl 1. 我的答案中包含了缓慢的迭代版本和快速的递归版本,底部的链接中还有高度优化的NTT,它比FFT表现更好。2. 如果你指的是大小,那么递归FFT和NTT需要一个2的幂次方大小的数据才能正常工作,你对“NTT仅适用于索引1和1 + n/2”有什么想法? - Spektre
抱歉,您确实展示了一种递归版本,但我似乎不明白它是如何工作的。您提到这比普通FFT还要棘手一些,我能理解这一点,能否解释一下为什么它需要与普通FFT不同呢?也就是说,为什么我们不能用k*2^n + 1的形式的模p来替换复数?我之所以问是因为当我将复数替换为它们在模p下的等价物(n = 8),我的递归NTT正确输出索引[0]和[0 + 4],但对于[1、2、3、5、6、7]则失败了,我想知道这可能是为什么? - z.karl
嗯,我猜这可能是我的主要问题。我使用的是673 = 84 * 2 ^ 3 + 1作为质数,这意味着5是单位根的主元(第672个),这意味着609是2 ^ 3 = 8次单位根,因为8是我的变换长度。那个质数有什么问题吗?我认为它满足条件,不是吗? - z.karl
@z.karl,我不是数学家,只是一名程序员,多年前就处理过这些事情。我在这个答案中编写的规则似乎足以证明使用这个素数和参数查找器得到的NTT结果是正确的(我在大数计算中使用了它们)。 - Spektre
显示剩余4条评论

1
将Cooley-Tukey(复杂)FFT转换为模算术方法,即NTT,必须替换omega的复杂定义。为了使方法纯递归,您还需要基于当前信号大小为每个级别重新计算omega。这是可能的,因为当我们在调用树中向下移动时,最小适合模数会减少,因此用于根的模数适用于较低层。此外,由于我们使用相同的模数,因此在移动调用树时可以使用相同的生成器。另外,对于反变换,您应该采取额外的步骤来获取重新计算的omega a并将其作为omega使用:b = a ^ -1(通过使用逆模运算)。具体而言,b = invMod(a,N),其中N是所选的素数模数,满足b * a == 1(mod N)。
通过利用周期性重写涉及omega的表达式在模算术领域仍然有效。您还需要找到一种确定问题模数(一个质数)和一个有效生成器的方法。
我们注意到您的代码可以工作,但它不是MWE。我们使用常识扩展它,并为多项式乘法应用程序获得了正确的结果。您只需提供正确的omega值的幂次。
虽然您的代码可以工作,但和其他许多来源一样,每个级别都会出现双倍间距。然而,这并不导致递归非常干净;这实际上等同于基于当前信号大小重新计算omega,因为omega定义的功率与信号大小成反比。再次强调:将信号大小减半就像平方omega,这就像给予omega加倍的能量(这是对间距加倍所做的操作)。处理重新计算omega的方法的好处在于每个子问题本身更加完整。
有一篇论文展示了模块化方法的一些数学知识;这是Baktir和Sunar在2006年发表的一篇论文。请参阅本帖子末尾的论文。
您不需要将循环从n / 2扩展到n。
因此,是的,一些源说只需在模数算术方法中放置不同的omega定义即可,这是在掩盖许多细节。
另一个问题是,如果我们要执行卷积,则必须承认信号大小必须足够大,以便结果时域信号不会溢出。此外,即使幂非常大,寻找受模数主题影响的某些指数实现也可能很有用。

参考文献

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

0

您必须确保单位根确实存在。在R中,只有两个单位根: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下工作。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接