Description
给定n,b,c,d,e以及A0,A1,···,An−1，定义
xk=b×c^4k+d×c^2k+e
f(x)=Sigma(Aix^i),0<=i<=n-1
请你求出f(x0),f(x1),···,f(xn−1)对10^6+3取模的值。
Input
第一行包括五个整数n,b,c,d,e。
接下来一行包括n个整数，代表a0,a1,···,an−1。
N<=60000,保证给出的数字都为整数且均在 [0, 10^6]
Output
N行，第i行代表f(xi−1)
多项式多点求值可过，但似乎不是正解
#include<bits/stdc++.h> typedef std::vector<int> pol; typedef long long i64; typedef unsigned long long u64; const int P=1000003; int max(int a,int b){return a>b?a:b;} template<int p> i64 pw(int a,int n){ if(n<0)n+=p-1; i64 v=1; for(;n;n>>=1,a=i64(a)*a%p)if(n&1)v=v*a%p; return v; } int N,K,r[1<<20|111]; typedef double ld; struct cplx{ ld a,b; cplx(ld _x=0,ld _y=0):a(_x),b(_y){} cplx operator+(cplx w){return cplx(a+w.a,b+w.b);} cplx operator-(cplx w){return cplx(a-w.a,b-w.b);} cplx operator*(cplx w){return cplx(a*w.a-b*w.b,a*w.b+b*w.a);} }A[1<<20|111],B[1<<20|111],*Es[2][20],C[1<<20|111],D[1<<20|111]; void pre(int n){ for(N=2,K=0;N<n;N<<=1,++K); for(int i=1;i<N;++i)r[i]=r[i>>1]>>1|(i&1)<<K; memset(A,0,sizeof(cplx)*N); memset(B,0,sizeof(cplx)*N); memset(C,0,sizeof(cplx)*N); memset(D,0,sizeof(cplx)*N); } const ld pi=acos(-1); void dft(cplx*a,int t){ for(int i=1;i<N;++i)if(i<r[i])std::swap(a[i],a[r[i]]); for(int i=1,z=0;i<N;i<<=1,++z){ if(!Es[t==1][z]){ cplx*E=Es[t==1][z]=new cplx[i]; for(int j=0;j<i;++j)E[j]=cplx(cos(j*pi/i),t*sin(j*pi/i)); } cplx*E=Es[t==1][z]; for(int j=0;j<N;j+=i<<1){ cplx*b=a+j,*c=b+i; for(int k=0;k<i;++k){ cplx x=b[k],y=c[k]*E[k]; b[k]=x+y; c[k]=x-y; } } } if(t==-1)for(int i=0;i<N;++i)a[i].a/=N; } void mul(cplx*a,cplx*b){ for(int i=0;i<N;++i)a[i]=a[i]*b[i]; } void mov(const pol&x,cplx*a){ for(int i=0;i<x.size();++i)a[i].a=x[i]; } void mov(cplx*a,const pol&_x){ pol&x=const_cast<pol&>(_x); for(int i=0;i<x.size();++i)x[i]=i64(a[i].a+0.49)%P; } void mov(const pol&x,int*a){ memcpy(a,x.data(),x.size()*sizeof(int)); } void mov(int*a,const pol&_x){ pol&x=const_cast<pol&>(_x); memcpy(x.data(),a,x.size()*sizeof(int)); } void chk(pol&w){ int p=w.size(); while(p&&!w[p-1])--p; w.resize(p); } void rev(pol&w){ std::reverse(w.data(),w.data()+w.size()); } u64 _c[1007]; pol operator*(const pol&a,const pol&b){ pol c(a.size()+b.size()-1); if(c.size()<=600){ const int*as=a.data(),ap=a.size(),*bs=b.data(),bp=b.size(); int*cs=c.data(); for(int i=0;i<c.size();++i)_c[i]=0; for(int i=0;i<ap;++i)for(int j=0;j<bp;++j)_c[i+j]+=u64(as[i])*bs[j]; for(int i=0;i<c.size();++i)cs[i]=_c[i]%P; }else if(c.size()<=5000){ pre(c.size()); mov(a,A),mov(b,B); dft(A,1),dft(B,1); mul(A,B); dft(A,-1); mov(A,c); }else{ pre(c.size()); for(int i=0;i<a.size();++i){ A[i]=a[i]>>10; B[i]=a[i]&1023; } for(int i=0;i<b.size();++i){ C[i]=b[i]>>10; D[i]=b[i]&1023; } dft(A,1);dft(B,1);dft(C,1);dft(D,1); for(int i=0;i<N;++i){ cplx v1=A[i]*C[i]; cplx v2=A[i]*D[i]+B[i]*C[i]; cplx v3=B[i]*D[i]; A[i]=v1; B[i]=v2; C[i]=v3; } dft(A,-1);dft(B,-1);dft(C,-1); for(int i=0;i<c.size();++i){ c[i]=((i64(A[i].a+0.49)<<20)+(i64(B[i].a+0.49)<<10)+i64(C[i].a+0.49))%P; } } chk(c); return c; } pol operator*(const pol&w,int x){ pol a(w); for(int i=0;i<a.size();++i)a[i]=i64(a[i])*x%P; return a; } pol operator-(const pol&a,const pol&b){ pol c(max(a.size(),b.size())); mov(a,c.data()); for(int i=0;i<b.size();++i)c[i]=(c[i]-b[i]+P)%P; chk(c); return c; } pol inv(pol a){ if(a.size()==1){ a.resize(1); a[0]=pw<P>(a[0],P-2); return a; } pol b=inv(pol(a.data(),a.data()+(a.size()-1)/2+1)); pol c=b*b*a; c=b*2-c; c.resize(a.size()); return c; } pol operator/(pol a,pol b){ chk(a),chk(b); rev(a),rev(b); int sz=a.size()-b.size()+1; b.resize(sz); pol c=a*inv(b); c.resize(sz); rev(c); chk(c); return c; } pol operator%(const pol&a,const pol&b){ if(a.size()<b.size())return a; return a-a/b*b; } int cal(const pol&w,i64 x){ int y=0; for(int p=w.size()-1;p>=0;--p)y=(y*x+w[p])%P; return y; } int n,b,c,d,e,a[67777],as[67777]; struct Q{int x,id;}qs[67777]; bool operator<(Q a,Q b){return a.x<b.x;} pol tr[133333]; void calc0(int w,int L,int R){ if(L==R){ tr[w]=pol(2); tr[w][1]=1; tr[w][0]=(P-qs[L].x)%P; return; } int M=(L+R)/2; calc0(w<<1,L,M); calc0(w<<1^1,M+1,R); tr[w]=tr[w<<1]*tr[w<<1^1]; } void calc1(const pol&a,int w,int L,int R){ if(R-L<=600){ for(int i=L;i<=R;++i)as[qs[i].id]=cal(a,qs[i].x); return; } int M=(L+R)/2; calc1(a%tr[w<<1],w<<1,L,M); calc1(a%tr[w<<1^1],w<<1^1,M+1,R); } int main(){ scanf("%d%d%d%d%d",&n,&b,&c,&d,&e); for(int i=0;i<n;++i)scanf("%d",a+i); for(int i=0;i<n;++i){ qs[i].x=(b*pw<P>(c,4*i)+d*pw<P>(c,2*i)+e)%P; qs[i].id=i; } std::sort(qs,qs+n); tr[0]=pol(n); mov(a,tr[0]); calc0(1,0,n-1); calc1(tr[0],1,0,n-1); for(int i=0;i<n;++i)printf("%d ",(as[i]+P)%P); return 0; }