FFT的常数优化

卡得一手好常数。。学习了。。(似乎只对FFT有效)

JZOJ 4349

#include <bits/stdc++.h>
#define LL long long
#define DB long double
using namespace std;
const int mo=100003;
const DB pi=acos(-1);
int K,T,p[360005],q[360005],f[70005],n,m,k,a[70005],rev[70005],ans;
DB COS[70005],SIN[70005],Cos[70005],Sin[70005];
struct cp{
    DB x,y;
    cp (DB p=0,DB q=0){x=p,y=q;}
}w[70005],N[35005];
cp operator+(cp a,cp b){
    return cp(a.x+b.x,a.y+b.y);
}
cp operator-(cp a,cp b){
    return cp(a.x-b.x,a.y-b.y);
}
cp operator*(cp a,cp b){
    return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}
LL po(LL x,int y){
    LL z=1;
    for (;y;y>>=1,x=x*x%mo)
    if (y&1) z=z*x%mo;
    return z;
}
LL C(int n,int m){
    return n<m?0:(n<mo?1ll*p[n]*q[m]%mo*q[n-m]%mo:C(n%mo,m%mo)*C(n/mo,m/mo)%mo);
}
void fft(cp a[],int n,int d=1){
    for (int i=0;i<n;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
    for (int m=2,k=1;m<=n;k=m,m<<=1){
        cp wn;
        wn=d>0?cp(COS[m],SIN[m]):cp(Cos[m],Sin[m]);
        N[0]=cp(1,0);
        for (int i=1;i<k;++i) N[i]=N[i-1]*wn;
        for (int i=0;i<n;i+=m)
        for (int j=i;j<i+k;++j){
            cp u=a[j],v=a[j+k]*N[j-i];
            a[j]=u+v; a[j+k]=u-v;
        }
    }
    if (d==-1) for (int i=0;i<=n;++i) a[i].x/=n;
}
void preF(){
    f[1]=9; f[2]=243; f[3]=16224;
    int t=12,X=299,Y=1<<t;
    for (int i=4;i<=60000;++i){
        int x=1,y=0;
        for (int j=0;j<=4;++j,x*=3)
            (f[i]+=(C(4,j)*x*(Y-X+y)%mo))%=mo,(y+=C(t,i-j-1))%=mo;
        f[i]=(f[i]%mo+mo)%mo;
        for (int j=1;j<=6;++j)
            X=(X*2-C(t,i-1))%mo,++t;
        X=(X+C(t,i))%mo; Y=Y*64%mo;
    }
}
int main(){
    for (int i=2;i<=65536;i<<=1)
        COS[i]=cos(pi*2/i),Cos[i]=cos(pi*-2/i),
        SIN[i]=sin(pi*2/i),Sin[i]=sin(pi*-2/i);
    scanf("%d",&T);
    p[0]=q[0]=1;
    for (int i=1,x;i<mo;++i){
        for (x=i;x%mo==0;x/=mo);
        p[i]=1ll*p[i-1]*x%mo;
    }
    q[mo-1]=po(p[mo-1],mo-2);
    for (int i=mo-2,x;i;--i){
        for (x=i+1;x%mo==0;x/=mo);
        q[i]=1ll*q[i+1]*x%mo;
    }
    preF();
    while (T--){
        scanf("%d%d",&n,&K);
        for (int i=1;i<=n;++i) a[i]=f[i];
        for (m=1;m<n;m<<=1);
        for (int i=n+1;i<=m;++i) a[i]=0; n=m;
        for (m=2,k=1;m<=n;k=m,m<<=1){
            for (int i=0;i<m;++i) rev[i]=rev[i>>1]>>1|(i&1)*(m>>1);
            for (int i=1;i<=n;i+=m){
                w[0]=cp(2,0);
                for (int j=i;j<i+k;++j) w[j-i+1]=cp(a[j]+a[j+k],a[j]-a[j+k]);
                for (int j=k+1;j<m;++j) w[j]=cp(0,0);
                fft(w,m);
                for (int j=0;j<m;++j) w[j]=w[j]*w[j];
                fft(w,m,-1);
                for (int j=i;j<i+m-1;++j) a[j]=(LL)round(w[j-i+1].x/4)%mo;
                a[i+m-1]=(LL)round(w[0].x/4-1)%mo;
            }
        }
        ans=0;
        for (int i=K;i<=n;++i) (ans+=a[i])%=mo;
        printf("%d
",(ans%mo+mo)%mo);
    }
    return 0;
}