void FFT(complex a[],int n,int fl){ for (int i=1,j=n/2;i<n;i++){ if (i<j) {complex t=a[i];a[i]=a[j];a[j]=t;}; int k; for (k=(n>>1);j&k;j^=k,k>>=1); j^=k; } for (int i=2;i<=n;i<<=1){ complex w;w.r=cos(fl*2*pi/i);w.i=sin(fl*2*pi/i); for (int j=0;j<n;j+=i){ complex wi;wi.r=1;wi.i=0; for (int k=j;k<j+i/2;k++){ complex u=a[k],v=a[k+i/2]*wi; a[k]=u+v;a[k+i/2]=u-v; wi=wi*w; } } } if (fl==-1) for (int i=0;i<n;i++) a[i].r/=n; }
fl==1求点值 fl==-1插值
_________________________
DFT
void DFT(){ for (int i=0;i<n;i++){ cp t=(cp){cos(2*pi/n*i),sin(2*pi/n*i)},bas=(cp){1,0}; for (int j=0;j<n;j++){ tmp[i]=tmp[i]+bas*a[j]; bas=bas*t; } } for (int i=0;i<n;i++) a[i]=tmp[i]; } void IDFT(){ for (int i=0;i<n;i++) P[i].r=P[i].i=0; for (int i=0;i<n;i++){ cp t=(cp){cos(-2*pi/n*i),sin(-2*pi/n*i)},bas=(cp){1,0}; for (int j=0;j<n;j++){ P[i]=P[i]+tmp[j]*bas; bas=bas*t; } } for (int i=0;i<n;i++) P[i].r=(int)(P[i].r/n+0.5); }
-------------------------------------------------------------------------
CODECHEF JUNE15 MOREFB
系数为多项式的分治FFT
#include <bits/stdc++.h> #define LL long long #define LDB long double using namespace std; const LDB pi=acos(-1); const LL mo=99991; int a[1000001],n,k; struct data{ LL a_1,a_n; }; inline data operator*(data a,data b) {return((data){(a.a_n*b.a_n+a.a_1*b.a_1)%mo,(a.a_n*b.a_n+a.a_1*b.a_n+a.a_n*b.a_1)%mo});}; struct cp{ LDB r_1,r_n,r_n2,i_1,i_n,i_n2; void set0(){ r_1=r_n=r_n2=i_1=i_n=i_n2=0; } void set(data a){ set0(); r_1=a.a_1;r_n=a.a_n; } }A[1000001],B[1000001]; inline cp operator *(cp a,cp b){ return((cp) { a.r_1*b.r_1-a.i_1*b.i_1, a.r_n*b.r_1+a.r_1*b.r_n-a.i_n*b.i_1-a.i_1*b.i_n, a.r_n2*b.r_1+a.r_n*b.r_n+a.r_1*b.r_n2-a.i_n2*b.i_1-a.i_n*b.i_n-a.i_1*b.i_n2, a.r_1*b.i_1+a.i_1*b.r_1, a.r_n*b.i_1+a.r_1*b.i_n+a.i_n*b.r_1+a.i_1*b.r_n, a.r_n2*b.i_1+a.r_n*b.i_n+a.r_1*b.i_n2+a.i_n2*b.r_1+a.i_n*b.r_n+a.i_1*b.r_n2 }); } inline cp operator +(cp a,cp b){ return((cp) { a.r_1+b.r_1, a.r_n+b.r_n, a.r_n2+b.r_n2, a.i_1+b.i_1, a.i_n+b.i_n, a.i_n2+b.i_n2, }); } inline cp operator -(cp a,cp b){ return((cp) { a.r_1-b.r_1, a.r_n-b.r_n, a.r_n2-b.r_n2, a.i_1-b.i_1, a.i_n-b.i_n, a.i_n2-b.i_n2, }); } data getkth(int po){ data ret;ret.a_1=1;ret.a_n=0; data bas;bas.a_1=0;bas.a_n=1; for (;po;bas=bas*bas){ if (po&1) ret=ret*bas; po>>=1; } return(ret); } void FFT(cp a[],int n,int fl){ for (int i=1,j=n/2;i<n;i++){ if (i<j) {cp t=a[i];a[i]=a[j];a[j]=t;}; int k; for (k=(n>>1);j&k;j^=k,k>>=1); j^=k; } for (int i=2;i<=n;i<<=1){ cp w;w.set0(); w.r_1=cos(fl*2*pi/i);w.i_1=sin(fl*2*pi/i); for (int j=0;j<n;j+=i){ cp wi;wi.set0(); wi.r_1=1;wi.i_1=0; for (int k=j;k<j+i/2;k++){ cp u=a[k],v=a[k+i/2]*wi; a[k]=u+v;a[k+i/2]=u-v; wi=wi*w; } } } if (fl==-1) for (int i=0;i<n;i++) a[i].r_n2/=n,a[i].r_n/=n,a[i].r_1/=n; } vector <data> MUL(vector <data> &a,vector <data> &b){ int n=1; while (n<a.size()+b.size()) n<<=1; for (int i=0;i<n;i++){ if (i<a.size()) A[i].set(a[i]);else A[i].set0(); if (i<b.size()) B[i].set(b[i]);else B[i].set0(); } FFT(A,n,1);FFT(B,n,1); for (int i=0;i<n;i++) A[i]=A[i]*B[i]; FFT(A,n,-1); vector <data> ret;ret.resize(a.size()+b.size()); for (int i=0;i<a.size()+b.size();i++){ LL tn2=(A[i].r_n2+0.5),tn=(A[i].r_n+0.5),t1=(A[i].r_1+0.5); ret[i].a_n=(tn2+tn)%mo;ret[i].a_1=(tn2+t1)%mo; } return(ret); } void show(vector <data> tmp1,vector <data> tmp2){ for (int i=0;i<tmp1.size();i++) printf("%lld %lld|",tmp1[i].a_n,tmp1[i].a_1);printf(" "); for (int i=0;i<tmp2.size();i++) printf("%lld %lld|",tmp2[i].a_n,tmp2[i].a_1);printf(" "); vector <data> res=MUL(tmp1,tmp2); for (int i=0;i<res.size();i++) printf("%lld %lld|",res[i].a_n,res[i].a_1);printf(" "); } vector <data> solve(int l,int r){ if (l==r){ vector <data> ret;ret.resize(2); ret[0].a_1=1;ret[0].a_n=0; ret[1]=getkth(a[l]); return(ret); } int mid=(l+r)>>1; vector <data> tmp1=solve(l,mid); vector <data> tmp2=solve(mid+1,r); // show(tmp1,tmp2); return(MUL(tmp1,tmp2)); } int main(){ scanf("%d%d",&n,&k); for (int i=1;i<=n;i++) scanf("%d",&a[i]); vector <data> ans=solve(1,n); printf("%lld ",ans[k].a_n%mo); }