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++]);
}

其实两个方法构造出的根都是类似于${\omega }_n^i=x^{\frac{i}{n}}$，x随i从0到n-1呈周期性变化，满足${\omega }_n^{\frac{n}{2}}=-1$，${\omega }_{2n}^{2i}={\omega }_n^i$
从而将$A(x)=a_0+a_{1}x+\cdots +a_{n-1}{x}^{n-1}$的求值问题分治为
$$\begin{gathered} A_0(x)=a_0+a_{2}x+a_4{x}^2+\cdots +a_{n-2}{x}^{\frac{n}{2}-1} \\ {A}_1(x)=a_1+a_{3}x+a_5{x}^2+\cdots +a_{n-1}{x}^{\frac{n}{2}-1} \\ A(x)=A_0(x^2)+x{A}_1(x^2) \end{gathered}$$
$x^2在x={\omega }_n^i和x={\omega }_n^{i+\frac{n}{2}}$处是相等的，n个根的求值问题就转化成了2个n/2个根的求值问题
FFT原理学习

---

Upd：
胡小兔的blog
FFT的代码技巧

把 $n$ 次单位根存在 $[\frac{n}{2},n-1]$ 的位置上，对 $\frac{n}{2},\frac{n}{4}...$ 同理，于是对于 $i\in [0,\frac{2}{n})$ 有 $w_i=w_{i<<1}$

对于 IDFT 可以翻转 $a_{1...n-1}$，（$a_0$ 不翻转），可以发现第 $i$ 个位置的值是 $a_0+{\sum }_{j=1}^{n-1}{a}_{n-j}{w}_n^{ij}$，它与 $a_0+{\sum }_{j=1}^{n-1}{a}_j{w}_n^{-ij}$ 相等。