洛谷 P4191 [CTSC2010]性能优化 fft

题目大意:
定义f与g的循环卷积,
( f ∗ g ) ( k ) = ∑ k ≡ i + j ( m o d   n ) f ( i ) ∗ g ( j ) (f*g)(k)=\sum_{k≡i+j(mod\ n)}f(i)*g(j) (fg)(k)=ki+j(mod n)f(i)g(j)
给定 n n n项的多项式 A A A B B B和常数 C C C
A ∗ B C A*B^C ABC的所有项 m o d   ( n + 1 ) mod\ (n+1) mod (n+1)的值。
其中, n n n只含小于 10 10 10的质因子, n + 1 n+1 n+1为质数。
n ≤ 5 ∗ 1 0 5 n≤5*10^5 n5105 C ≤ 1 0 9 C≤10^9 C109

分析:
因为 n + 1 n+1 n+1是质数, x n ≡ 1 x^n≡1 xn1
所以循环卷积得出来的多项式与一般卷积多项式的点值相同。
如果我们使用 n n n个点值插出来的多项式就是循环卷积多项式。
发现 n n n只含小于 10 10 10的质因子。
如果只含 2 2 2这个质因子就是一般的ntt了,考虑有其他质因子怎么做。
举个例子,用 5 5 5这个质因子来算。
我们把当前多项式变成
( a 0 + a 5 x 5 ) + x ( a 1 + a 6 x 5 ) + x 2 ( a 2 + a 7 x 5 ) + x 3 ( a 3 + a 8 x 5 ) + x 4 ( a 4 + a 9 x 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) (a0+a5x5)+x(a1+a6x5)+x2(a2+a7x5)+x3(a3+a8x5)+x4(a4+a9x5)
即按余数分类。
分别算出 a 0 + a 5 x a_0+a_5x a0+a5x a 1 + a 6 x a_1+a_6x a1+a6x,……的点值。
显然此时要把点值要变成5次方,然而这些多项式点值自变量为 w k / 5 i w_{k/5}^i wk/5i,5次方即为 w k i w_{k}^i wki
所以直接加上对应点值,再乘上当前单位复根的0~4次方。
一般ntt维护一个 w n w_n wn w w w,此时要用数组维护 w n [ i ] w_n[i] wn[i] w [ i ] w[i] w[i]。表示单位复根的 i i i次方与当前自变量的 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]);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值