One afternoon by the accuracy of the problem continued ... First you can find a polynomial of the geometric series form, and then like the POJ matrix Series, multiply it. Using complex decimal values to represent multiple polynomial operations will brush the precision ... You should use an int to save the polynomial, and then convolution the DfT into a complex number, then idft back to the real number after convolution. Note that two M-times polynomial convolution will become a polynomial of 2m, and the latter half of the polynomial needs to be cleared 0.
#include <cstdio>#include<cstring>#include<cmath>#include<algorithm>using namespacestd;Const intmaxn=1024x768* +*4;#defineDouble Long DoubleConst DoublePi=acos (-1);structcomp{Doublex, y; Comp () {} comp (DoubleADoubleb) {x=a;y=b;} Compoperator+(ConstComp &a) {returnComp (x+a.x,y+a.y);} Compoperator-(ConstComp &a) {returnComp (x-a.x,y-a.y);} Compoperator*(ConstComp &a) {returnComp (x*a.x-y*a.y,x*a.y+y*a.x);}} ;//A: Storing the original polynomial B: storing the convolution of the original polynomial C: storing answer D: Storing the N/2 of the original polynomialintMoD;intA[MAXN],C[MAXN],D[MAXN],E[MAXN];voidFFT (comp* A,intNintSign ) { for(intI=1, j=0, k=n;i<n;++i,k=N) { DoJ^= (k>>=1); while(j<k);if(i<j) Swap (A[i],a[j]); } for(intj=2; j<=n;j<<=1){ intM=j>>1; Comp wn (cos (pi*2/j), Sign*sin (pi*2/j)); for(Comp *p=a;p!=a+n;p=p+j) {Comp W (1,0); for(intk=0; k<m;++k,w=w*WN) {Comp T=p[m+k]*w;p[m+k]=p[k]-t;p[k]=p[k]+T; } } } if(sign==-1){ for(intI=0; i<n;++i) a[i].x/=N; }}intn=1;intm;intMoDoublex) { return(((int) Floor (x+0.5))%mod+mod)%MoD;}voidMultint*a,int*b,int*Res) { StaticComp TMP1[MAXN],TMP2[MAXN]; for(intI=0; i<n;++i) Tmp1[i]=comp (A[i],0), Tmp2[i]=comp (B[i],0); FFT (Tmp1,n,1); FFT (Tmp2,n,1); for(intI=0; i<n;++i) tmp1[i]=tmp1[i]*Tmp2[i]; FFT (Tmp1,n,-1); for(intI=0; i<n;++i) res[i]=Mo (tmp1[i].x);}voidQsum (intN) { if(n==1){ for(intI=0; i<n;++i) c[i]=A[i]; for(intI=0; i<n;++i) d[i]=A[i]; }Else{qsum (n>>1); Mult (c,d,e);//for (int i=0;i<n;++i)//E[i]=c[i]*d[i]+c[i]; for(intI=0; i<n;++i) C[i]=mo (c[i]+E[i]); memset (c+ (n>>1),0,sizeof(comp) * (n>>1)); if(n&1) {mult (c,a,e); for(intI=0; i<n;++i) C[i]=mo (a[i]+E[i]); memset (c+ (n>>1),0,sizeof(comp) * (n>>1)); } mult (D,d,d); memset (d+ (n>>1),0,sizeof(comp) * (n>>1)); if(n&1) {mult (d,a,d); memset (d+ (n>>1),0,sizeof(comp) * (n>>1)); } }}intMain () {scanf ("%d%d",&m,&MoD); intN,O,S,U;SCANF ("%d%d%d%d",&n,&o,&s,&u); N=min (n,m); for(intI=1; i<=m;++i) { intt=i%MoD; A[i]= (o*t*t+s*t+u)%MoD; } while(n<=m) n<<=1; n<<=1; Qsum (n);p rintf ("%d\n", C[m]); return 0;}
BZOJ4332[JSOI2012] Points Snack