FFT&NTT模板

一开始抄了jiry_2的板子，直到有一天，我被卡精度了
下定决心弄一个永久的板子，四处借鉴写法
单位复数根还是预处理的好
尽量避免乘法，伤精度

FFT UOJ#34

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,b,a) for(int i=b;i>=a;i--)
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define mset(a,x) memset(a,x,sizeof(a))
using namespace std;
typedef long long ll;
typedef double db;
char ch;
const int N=4e5+5;
const db pi=acos(-1);
struct Z
{
    db x,y;
    Z(db xx=0,db yy=0){x=xx,y=yy;}
    friend Z operator +(Z a,Z b){return Z(a.x+b.x,a.y+b.y);}
    friend Z operator -(Z a,Z b){return Z(a.x-b.x,a.y-b.y);}
    friend Z operator *(Z a,Z b){return Z(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}a[N],b[N],W[N];
int n,m,lg,len,bitrev[N];
void dft(Z *a,int n,int sig=1)
{
    fo(i,0,n-1) if(i<bitrev[i]) swap(a[i],a[bitrev[i]]);
    for(int i=2,gap=n>>1;i<=n;i<<=1,gap>>=1)
        for(int j=0;j<n;j+=i)
        {
            Z *l=a+j,*r=a+j+(i>>1),*w=sig>0?W:W+n;
            fo(k,0,(i>>1)-1)
            {
                Z u=*l,v=*w * *r;
                *l=u+v,*r=u-v;
                ++l,++r,w=sig>0?w+gap:w-gap;
            }
        }
    if(sig==-1) fo(i,0,n-1) a[i].x/=n;
}
int main()
{
    freopen("uoj#34.in","r",stdin);
    scanf("%d %d",&n,&m);
    fo(i,0,n) scanf("%lf",&a[i].x);
    fo(j,0,m) scanf("%lf",&b[j].x);
    lg=0,len=1;
    while(len<(max(n+1,m+1)<<1)) ++lg,len<<=1;
    fo(i,0,len-1) bitrev[i]=(bitrev[i>>1]>>1)|((i&1)<<(lg-1));
    fo(i,0,len) W[i]=Z(cos(2*pi*i/len),sin(2*pi*i/len));
    dft(a,len);dft(b,len);
    fo(i,0,len-1) a[i]=a[i]*b[i];
    dft(a,len,-1);
    fo(i,0,n+m) printf("%d ",int(a[i].x+0.5));
    return 0;
}

NTT

这**短得不像话

void dft(ll *a,int n,int sig=1)
{
    fo(i,0,n-1) if(i<bitrev[i]) swap(a[i],a[bitrev[i]]);
    for(int i=2,gap=n>>1;i<=n;i<<=1,gap>>=1)
        for(int j=0;j<n;j+=i)
            for(ll k=0,*l=a+j,*r=a+j+(i>>1),*w=sig>0?W:W+n;k<(i>>1);++k,++l,++r,w=sig>0?w+gap:w-gap)
            {
                ll u=*l,v=(*w * *r)%mo;
                *l=(u+v)%mo,*r=(u-v+mo)%mo;
            }
    if(sig==-1)
    {
        ll inv=qmi(len,mo-2);
        fo(i,0,n-1) a[i]=a[i]*inv%mo;
    }
}