COGS 2294. [HZOI 2015] 释迦（FFT mod any prime）

最新推荐文章于 2021-02-18 17:11:25 发布

clover_hxy

最新推荐文章于 2021-02-18 17:11:25 发布

阅读量645

点赞数

分类专栏： FFT 数论文章标签： fft

本文链接：https://blog.csdn.net/clover_hxy/article/details/73121942

版权

数论同时被 2 个专栏收录

158 篇文章 0 订阅

订阅专栏

FFT

22 篇文章 0 订阅

订阅专栏

本文介绍了一种解决多项式乘积问题的方法，通过使用三次NTT（Number Theoretic Transform）和中国剩余定理来求解任意模数下的多项式乘积问题。文章详细解释了如何选择合适的模数进行计算，并提供了完整的代码实现。

摘要由CSDN通过智能技术生成

题目描述

传送门

题目大意：给两个次数界为n的多项式,求这两个多项式的乘积，输出前x的0次项到n-1次项的系数 mod 23333333

题解

NTT只能求在FFT模数下的值。对于任意模数的题来说，我们可以选择三个FFT模数分别做NTT，最后用中国剩余定理合并。
一次卷积后每个数可以达到 $10^{23}(n*mod^2)$ 左右，所以我们需要选择三个乘积大于 $10^{23}$ 的FFT模数。合并的时候 $10^{23}>2^{64}$ ，所以我们不能直接用中国剩余定理，需要一点小技巧。

a n s \equiv a 1 (m o d m 1)

$ans\equiv a1 \space(mod\space m1)$

a n s \equiv a 2 (m o d m 2)

$ans\equiv a2 \space(mod \space m2)$

a n s \equiv a 3 (m o d m 3)

$ans\equiv a3\space (mod \space m3)$
先用中国剩余定理合并

a1,a2 $a1,a2$ ,得到

ans≡A (mod M) $ans\equiv A \space(mod \space M)$
其中

A=a1∗m2∗m2−1+a2∗m1∗m1−1,M=m1∗m2 $A=a1*m2*m2^{-1}+a2*m1*m1^{-1},M=m1*m2$

ans=k∗M+A=x∗m3+a3 $ans=k*M+A=x*m3+a3$

k * M \equiv a 3 - A (m o d m 3)

$k*M\equiv a3-A \space (mod \space m3)$

k \equiv (a 3 - A) * M - 1 (m o d m 3)

$k\equiv (a3-A)*M^{-1} \space (mod \space m3)$
求出k后带入

ans=k∗M+A $ans=k*M+A$ ,然后最后ans对23333333取模即可。
中间做乘法的过程需要用到快速乘，否则会炸long long

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 500003
#define mod 23333333
#define LL long long 
using namespace std; 
LL M[10],f[4][N],g[4][N],ans[N];
int n,m,n1,L,R[N];
LL quickpow(LL num,LL x,LL p)
{
    LL base=num%p; LL ans=1;
    while (x) {
        if (x&1) ans=ans*base%p;
        x>>=1;
        base=base*base%p;
    }
    return ans;
}
void NTT(LL a[N],int n,int opt,LL p)
{
    for (int i=0;i<n;i++) 
     if (i>R[i]) swap(a[i],a[R[i]]);
    for (int i=1;i<n;i<<=1) {
        LL wn=quickpow(3,(p-1)/(i<<1),p);
        for (int p1=i<<1,j=0;j<n;j+=p1){
            LL w=1;
            for (int k=0;k<i;k++,w=w*wn%p){
                LL x=a[j+k],y=w*a[j+k+i]%p;
                a[j+k]=(x+y)%p; a[j+k+i]=(x-y+p)%p;
            }
        }
    }
    if (opt==-1) reverse(a+1,a+n);
}
void solve(LL a[N],LL b[N],LL p)
{
//  for (int i=0;i<=n1;i++) a[i]%=p;
//  for (int i=0;i<=n1;i++) b[i]%=p;
    NTT(a,n1,1,p); NTT(b,n1,1,p);
    for (int i=0;i<n1;i++) a[i]=a[i]*b[i]%p;
    NTT(a,n1,-1,p);
    LL rev=quickpow(n1,p-2,p);
    for (int i=0;i<n1;i++) a[i]=a[i]*rev%p;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
    if (!b) {
        x=1; y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    LL t=y;
    y=x-(a/b)*y;
    x=t;
}
LL mul(LL num,LL x,LL p)
{
    LL ans=0; LL base=num%p;
    while (x) {
        if (x&1) ans=(ans+base)%p;
        x>>=1;
        base=(base+base)%p;
    }
    return ans;
}
LL china(LL a1,LL a2,LL a3)
{
    LL MM=M[1]*M[2];
    LL x,x1; x=quickpow(M[2],M[1]-2,M[1]); x1=quickpow(M[1],M[2]-2,M[2]);
    LL A=(mul(a1*M[2]%MM,x%MM,MM)+mul(a2*M[1]%MM,x1%MM,MM))%MM;
    LL k=(a3-A)%M[3]*quickpow(MM,M[3]-2,M[3])%M[3];
    k=(k%M[3]+M[3])%M[3];
    return ((k%mod)*(MM%mod)%mod+A)%mod;
}
int main()
{
    freopen("annona_squamosa.in","r",stdin);
    freopen("annona_squamosa.out","w",stdout);
    scanf("%d",&n);
    M[1]=998244353,M[2]=1004535809,M[3]=469762049;
    for (int i=0;i<n;i++) scanf("%lld",&f[1][i]),f[2][i]=f[3][i]=f[1][i];
    for (int i=0;i<n;i++) scanf("%lld",&g[1][i]),g[2][i]=g[3][i]=g[1][i];
    m=2*n;
    for (n1=1;n1<=m;n1<<=1) L++;
    for (int i=0;i<=n1;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    for (int i=1;i<=3;i++) 
     solve(f[i],g[i],M[i]);
    for (int i=0;i<n;i++) printf("%lld\n",china(f[1][i],f[2][i],f[3][i]));
}