1. 程式人生 > >FFT/NTT模板 51nod1028 大數乘法 V2

FFT/NTT模板 51nod1028 大數乘法 V2

題目連結

FFT:

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define maxn 400005
using namespace std;
const int wei = 3, bit = 1e3;//壓三位
const double Pi = acos(-1);
struct complex
{
	double r,i;
	complex(double _r=0,double _i=0):r(_r),i(_i){}
	complex operator
+ (const complex &t)const{return complex(r+t.r,i+t.i);} complex operator - (const complex &t)const{return complex(r-t.r,i-t.i);} complex operator * (const complex &t)const{return complex(r*t.r-i*t.i,r*t.i+i*t.r);} }a1[maxn],a2[maxn],w,wn; void change(complex *a,int len) { for(int i=1
,j=len/2,k;i<len-1;i++) { if(i<j) swap(a[i],a[j]); for(k=len/2;j>=k;j-=k,k>>=1); j+=k; } } void fft(complex *a,int len,int flg) { for(int i=2;i<=len;i<<=1) { wn=complex(cos(2*flg*Pi/i),sin(2*flg*Pi/i)); for(int j=0;j<len;j+=i) { w=complex(1,0); for(int k=j;
k<j+i/2;k++) { complex t1=a[k],t2=w*a[k+i/2]; a[k]=t1+t2,a[k+i/2]=t1-t2; w=w*wn; } } } if(flg==-1) for(int i=0;i<len;i++) a[i].r/=len; } int n,len1,len2; long long ans[maxn]; char s1[maxn],s2[maxn]; int main() { int i,t; scanf("%s%s",s1,s2); len1=strlen(s1),len2=strlen(s2); for(i=0,t=len1%wei?0:-1;i<len1;i++) {if((len1-i)%wei==0) t++;a1[t].r=a1[t].r*10+(s1[i]-'0');} len1=t+1; for(i=0,t=len2%wei?0:-1;i<len2;i++) {if((len2-i)%wei==0) t++;a2[t].r=a2[t].r*10+(s2[i]-'0');} len2=t+1; n=1;while(n<len1+len2) n<<=1; change(a1,n),change(a2,n); fft(a1,n,1),fft(a2,n,1); for(int i=0;i<n;i++) a2[i]=a1[i]*a2[i]; change(a2,n); fft(a2,n,-1); for(int i=0;i<len1+len2-1;i++) ans[i]=(long long)(a2[i].r+0.5); for(int i=len1+len2-2;i>0;i--) { ans[i-1]+=ans[i]/bit; ans[i]%=bit; } for(i=0;ans[i]==0&&i<len1+len2-2;i++); printf("%lld",ans[i++]); while(i<len1+len2-1) printf("%03lld",ans[i++]); }

NTT:

#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long
#define maxn 400005
using namespace std;
LL mod = 998244353, G = 3;
LL a1[maxn],a2[maxn],w,wn;
int len1,len2;
char s1[maxn],s2[maxn];
inline LL ksm(LL a,int b){
    LL s=1;
    for(;b;b>>=1,a=a*a%mod) if(b&1) s=s*a%mod;
    return s;
}
inline void change(LL *a,int len)
{
    for(int i=1,j=len/2,k;i<len-1;i++)
    {
        if(i<j) swap(a[i],a[j]);
        for(k=len/2;j>=k;j-=k,k>>=1);
        j+=k;
    }
}
inline void ntt(LL *a,int len,int flg)
{
    change(a,len);
    for(int i=2;i<=len;i<<=1)
    {
        if(flg==1) wn=ksm(G,(mod-1)/i);
        else wn=ksm(G,mod-1-(mod-1)/i);
        for(int j=0;j<len;j+=i)
        {
            w=1;
            for(int k=j;k<j+i/2;k++)
            {
                LL u=a[k],v=w*a[k+i/2]%mod;
                a[k]=(u+v)%mod,a[k+i/2]=(u-v+mod)%mod;
                w=w*wn%mod;
            }
        }
    }
    if(flg==-1){
        LL ni=ksm(len,mod-2);
        for(int i=0;i<len;i++) a[i]=a[i]*ni%mod;
    }
}
int main()
{
    scanf("%s%s",s1,s2);
    len1=strlen(s1),len2=strlen(s2);
    for(int i=0;i<len1;i++) a1[i]=s1[i]-'0';
    for(int i=0;i<len2;i++) a2[i]=s2[i]-'0';
    int len=1;while(len<len1+len2-1) len<<=1;
    ntt(a1,len,1),ntt(a2,len,1);
    for(int i=0;i<len;i++) a2[i]=a1[i]*a2[i]%mod;
    ntt(a2,len,-1);
    for(int i=len1+len2-2;i>0;i--) a2[i-1]+=a2[i]/10,a2[i]%=10;
    int i;
    for(i=0;i<len1+len2-2&&a2[i]==0;i++);
    printf("%lld",a2[i++]);
    while(i<len1+len2-1) printf("%lld",a2[i++]);
}

其實兩個方法構造出的根都是類似於 ω n i = x i n \omega _n^i=x^{\frac in} ,x隨i從0到n-1呈週期性變化,滿足 ω n n 2 = 1 \omega_n^{\frac n2}=-1 ω 2 n 2 i = ω n i \omega_{2n}^{2i}=\omega_n^i
從而將 A ( x ) = a 0 + a 1 x + + a n 1 x n 1 A(x)=a_0+a_1x+\dots+a_{n-1}x^{n-1} 的求值問題分治為
A 0 ( x ) = a 0 + a 2 x + a 4 x 2 + + a n 2 x n 2 1 A 1 ( x ) = a 1 + a 3 x + a 5 x 2 + + a n 1 x n 2 1 A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 ) A_0(x)=a_0+a_2x+a_4x^2+\dots+a_{n-2}x^{\frac n2-1 }\\ A_1(x)=a_1+a_3x+a_5x^2+\dots+a_{n-1}x^{\frac n2-1}\\ A(x)=A_0(x^2)+xA_1(x^2)