洛谷 P4191 [CTSC2010]性能优化 fft

最新推荐文章于 2019-06-16 21:37:00 发布

Amber_lylovely

最新推荐文章于 2019-06-16 21:37:00 发布

阅读量236

点赞数

CC 4.0 BY-SA版权

分类专栏： fft-ntt

本文链接：https://blog.youkuaiyun.com/liangzihao1/article/details/87649602

fft-ntt 专栏收录该内容

27 篇文章

订阅专栏

本文探讨了在特定条件下，利用循环卷积和NTT（Number Theoretic Transform）算法进行多项式乘法的问题。重点讲解了如何通过预处理找到合适的单位复根，以及如何调整多项式的顺序和进行快速傅立叶变换，最终实现模意义下的多项式乘法。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目大意：
定义f与g的循环卷积，
$n)f(i)∗g(j)(f*g)(k)=\sum_{k≡i+j(mod\ n)}f(i)*g(j)$
给定 $n$ 项的多项式 $A$ ， $B$ 和常数 $C$ 。
求 $A*B^C$ 的所有项 $(n+1)mod\ (n+1)$ 的值。
其中， $n$ 只含小于 $10$ 的质因子， $n + 1$ 为质数。
$n≤5*10^5$ ， $C≤10^9$ 。

分析：
因为 $n + 1$ 是质数， $x^n≡1$ 。
所以循环卷积得出来的多项式与一般卷积多项式的点值相同。
如果我们使用 $n$ 个点值插出来的多项式就是循环卷积多项式。
发现 $n$ 只含小于 $10$ 的质因子。
如果只含 $2$ 这个质因子就是一般的ntt了，考虑有其他质因子怎么做。
举个例子，用 $5$ 这个质因子来算。
我们把当前多项式变成
$a_0+a_5x^5)+x(a_1+a_6x^5)+x^2(a_2+a_7x^5)+x^3(a_3+a_8x^5)+x^4(a_4+a_9x^5)$
即按余数分类。
分别算出 $a_0+a_5x$ ， $a_1+a_6x$ ，……的点值。
显然此时要把点值要变成5次方，然而这些多项式点值自变量为 $w_{k/5}^i$ ，5次方即为 $w_{k}^i$ 。
所以直接加上对应点值，再乘上当前单位复根的0~4次方。
一般ntt维护一个 $w_n$ 和 $w$ ，此时要用数组维护 $w_n[i]$ ， $w [i]$ 。表示单位复根的 $i$ 次方与当前自变量的 $i$ 次方即可。

代码：

// luogu-judger-enable-o2
#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long

const int maxn=5e5+7;

using namespace std;

int n,tot;
int p[maxn],nump[4];
LL mod,G,C;
LL f[maxn],g[maxn],w[17],v[17],b[maxn];

LL power(LL x,LL y)
{
    if (y==0) return 1;
    if (y==1) return x;
    LL c=power(x,y/2);
    c=c*c%mod;
    if (y&1) c=c*x%mod;
    return c;
}

void prework()
{
    int c=n;
    while (c%2==0) c/=2,p[++tot]=2,nump[0]++;
    while (c%3==0) c/=3,p[++tot]=3,nump[1]++;
    while (c%5==0) c/=5,p[++tot]=5,nump[2]++;
    while (c%7==0) c/=7,p[++tot]=7,nump[3]++;
    for (int i=2;i<n;i++)
    {
        if ((nump[0]) && (power(i,n/2)==1)) continue;
        if ((nump[1]) && (power(i,n/3)==1)) continue;
        if ((nump[2]) && (power(i,n/5)==1)) continue;
        if ((nump[3]) && (power(i,n/7)==1)) continue;
        G=i;
        break;
    }
}

void change(LL *a,int l,int r,int dep) //暴力调顺序
{
    if (l==r) return;
    int len=(r-l+1)/p[dep];
    for (int i=0;i<=r-l;i++)
    {
        int j=(i%p[dep])*len+i/p[dep];
        b[l+j]=a[l+i];
    }
    for (int i=0;i<=r-l;i++) a[l+i]=b[l+i];
    for (int i=0;i<p[dep];i++) change(a,l+i*len,l+(i+1)*len-1,dep+1);
}

void fft(LL *a,int f)
{
    int len=1;
    for (int i=tot;i>0;i--)
    {
        len*=p[i];
        LL wn;
        if (f==1) wn=power(G,n/len); //当前单位复根
             else wn=power(G,n-n/len);
        v[0]=1;
        for (int j=1;j<p[i];j++) v[j]=(v[j-1]*wn)%mod; //j次方单位复根
        for (int j=0;j<n;j++) b[j]=0; //记录当前层变化后的a
        for (int j=0;j<n;j+=len)
        {
            for (int l=0;l<p[i];l++) w[l]=1; //当前自变量重置
            for (int k=0;k<len;k++)
            {
                int num=len/p[i]; //小区间的大小
                for (int l=0;l<p[i];l++)
                {
                    b[j+k]=(b[j+k]+a[j+l*num+k%num]*w[l])%mod; //第k个自变量的值等于所有小区间内对应位置的值乘上对应次方的和。
                    w[l]=(w[l]*v[l])%mod; //自变量变为下一个
                }
            }
        }
        for (int i=0;i<n;i++) a[i]=b[i];
    }
    if (f==-1)
    {
        LL inv=power(n,mod-2);
        for (int i=0;i<n;i++) a[i]=a[i]*inv%mod;
    }
}

int main()
{
    scanf("%d%lld",&n,&C);
    mod=n+1;
    for (int i=0;i<n;i++) scanf("%lld",&f[i]);
    for (int i=0;i<n;i++) scanf("%lld",&g[i]);		
    prework();	
    change(f,0,n-1,1);
    change(g,0,n-1,1);
    fft(f,1);
    fft(g,1);
    for (int i=0;i<n;i++) f[i]=f[i]*power(g[i],C)%mod;
    change(f,0,n-1,1);
    fft(f,-1);
    for (int i=0;i<n;i++) printf("%lld\n",f[i]);
}