codechef Random Number Generator
题意:求常系数齐次递推数列的第n项
当递推系数非零项很多的时候暴力复位一次是O(k2)的
构造特征矩阵M的特征多项式f(x)=xk-∑i ci*xi
由f(M)=0,那么我们复位时只需要mod f(x)就好
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 | #include <stdio.h>#include <string.h>#include <algorithm>#define g() getchar()template<class Q>void Scan(Q&x){ char c; while(c=g(),c<48||c>57); for(x=0;c>47&&c<58;c=g())x=10*x+c-48;}#define G 3#define maxn 65536#define DFT(a,n) fft(a,n,0)#define IDFT(a,n) fft(a,n,1)#define INV(x) qpow(x,P-2)#define rep(i,a,b) for(int i=a;i<(b);++i)using namespace std;const int MAXN=(1<<16)+5;const int P=104857601;typedef long long ll;inline ll qpow(ll x,int k){ ll ans=1; for(;k;x=x*x%P,k>>=1) if(k&1)ans=ans*x%P; return ans;}int r[MAXN],w[2][MAXN],rn;inline void init_fft(int n,int b){ rn=INV(n); rep(i,1,n){ r[i]=(r[i>>1]>>1)|((i&1)<<(b-1)); }}inline void fft(int*a,int n,int t){ rep(i,1,n)if(i<r[i])swap(a[i],a[r[i]]); for(int i=2;i<=n;i<<=1) for(int j=0;j<n;j+=i) for(int k=0;k<i>>1;++k){ int re=1ll*a[(i>>1)+j+k]*w[t][maxn/i*k]%P; a[(i>>1)+j+k]=(a[j+k]+P-re)%P; a[j+k]=(a[j+k]+re)%P; } if(t)rep(i,0,n)a[i]=1ll*a[i]*rn%P;}int k,len,bit;inline void inv(int*a,int*b){ static int c[MAXN],l; b[0]=1,b[1]=0,l=1; rep(j,1,bit+1){ l<<=1; rep(i,0,l)c[i]=a[i]; rep(i,l,l<<1)b[i]=c[i]=0; init_fft(l<<1,j+1); DFT(b,l<<1),DFT(c,l<<1); rep(i,0,l<<1){ b[i]=1ll*b[i]*(2+P-1ll*b[i]*c[i]%P)%P; } IDFT(b,l<<1); rep(i,l,l<<1)b[i]=0; } }int m[MAXN],dm[MAXN],im[MAXN];int ta[MAXN],ra[MAXN];inline void mod(int*a){ DFT(ta,len),DFT(a,len); rep(i,0,len)a[i]=1ll*a[i]*ta[i]%P; IDFT(a,len); rep(i,0,k-1)ra[i]=a[(k-1<<1)-i]; rep(i,k-1,len)ra[i]=0; DFT(ra,len); rep(i,0,len)ra[i]=1ll*ra[i]*im[i]%P; IDFT(ra,len); rep(i,0,k-1)ta[i]=ra[k-2-i]; rep(i,k-1,len)ta[i]=0; DFT(ta,len); rep(i,0,len)ta[i]=1ll*ta[i]*dm[i]%P; IDFT(ta,len); rep(i,0,k)a[i]=(a[i]+P-ta[i])%P; rep(i,k,len)a[i]=0;}inline void Init(){ int n=maxn; int g=qpow(G,(P-1)/n); w[0][0]=w[0][n]=1; rep(i,1,n)w[0][i]=1ll*w[0][i-1]*g%P; rep(i,0,n+1)w[1][i]=w[0][n-i];}ll n;int a[MAXN],ans[MAXN],x[MAXN];int main(){ Init(); Scan(k),Scan(n),--n; len=k+1,bit=1; while(1<<bit<len)++bit; len=1<<bit; rep(i,0,k)Scan(a[i]); for(int i=k-1;~i;--i){ Scan(m[i]),m[i]=(P-m[i])%P; } if(k==1){ int res=qpow((P-m[0])%P,n%(P-1))*a[0]%P; printf("%d\n",res); return 0; } m[k]=1; rep(i,0,k+1)dm[i]=m[k-i]; inv(dm,im); len<<=1,++bit; rep(i,0,k+1)dm[i]=m[i]; rep(i,k+1,len)dm[i]=im[i]=0; DFT(dm,len),DFT(im,len); for(ans[0]=x[1]=1;n;n>>=1){ if(n&1){ rep(i,0,len)ta[i]=x[i]; mod(ans); } rep(i,0,len)ta[i]=x[i]; mod(x); } int res=0; rep(i,0,k){ res=(res+1ll*a[i]*ans[i])%P; } printf("%d\n",res); return 0;} |
hdu4914 Linear recursive sequence
题意:给出序列f,fn=(n>0)?a*fn-p+b*fn-q:1,求第n项
叉姐论文: k阶矩阵M的任意幂次都能用M小于k次的幂次来线性表示
所以我们构造多项式fn(x),[xi]fn(x)表示用小于k次的幂次线性表示时i次幂前的系数
用多项式乘法+快速幂计算fn(x),每次乘法后溢出k-1次的部分暴力还原
这题 mod 119，中间结果不会太大，所以直接用复数域上的fft，做完一次乘法后再取模就好
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 | #include <math.h>#include <stdio.h>#include <string.h>#include <algorithm>#define DFT(a) fft(a,len,0)#define IDFT(a) fft(a,len,1)#define rep(i,a,b) for(int i=a;i<(b);++i)using namespace std;const int MAXN=(1<<18)+5; typedef double db;const db Pi=acos(-1.);typedef long long ll;const int Mod=119;struct cp{ db r,i; cp(db r=0,db i=0):r(r),i(i){} friend cp operator+(const cp&a,const cp&b){ return cp(a.r+b.r,a.i+b.i); } friend cp operator-(const cp&a,const cp&b){ return cp(a.r-b.r,a.i-b.i); } friend cp operator*(const cp&a,const cp&b){ return cp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r); }}w[2][MAXN];inline cp Exp(db x){ return cp(cos(x),sin(x));}int r[MAXN];inline void init_fft(int n,int b){ w[0][0]=w[0][n]=cp(1.,0); cp g=Exp(Pi*2./n); rep(i,1,n)w[0][i]=w[0][i-1]*g; rep(i,0,n+1)w[1][i]=w[0][n-i]; rep(i,1,n){ r[i]=(r[i>>1]>>1)|((i&1)<<(b-1)); }} inline void fft(cp*a,int n,int t){ rep(i,1,n)if(i<r[i])swap(a[i],a[r[i]]); for(int i=2;i<=n;i<<=1) for(int j=0;j<n;j+=i) for(int k=0;k<i>>1;++k){ cp re=a[(i>>1)+j+k]*w[t][n/i*k]; a[(i>>1)+j+k]=a[j+k]-re; a[j+k]=a[j+k]+re; } if(t)rep(i,0,n)a[i].r/=n;}cp ta[MAXN],tb[MAXN];int n,A,B,p,q,len,bit;struct poly{ int a[MAXN]; inline void mul(const poly&o){ rep(i,0,q)ta[i]=cp(a[i]),tb[i]=cp(o.a[i]); rep(i,q,len)ta[i]=tb[i]=cp(); DFT(ta),DFT(tb); rep(i,0,len)ta[i]=ta[i]*tb[i]; IDFT(ta); rep(i,0,len)a[i]=ll(ta[i].r+0.5)%Mod; for(int i=len;i>=q;--i)if(a[i]){ a[i-p]=(a[i-p]+a[i]*A%Mod)%Mod; a[i-q]=(a[i-q]+a[i]*B%Mod)%Mod; } }}a,b;inline void qpow(poly&ans,poly&x,int k){ rep(i,0,len)ans.a[i]=x.a[i]=0; ans.a[0]=x.a[1]=1; for(;k;x.mul(x),k>>=1) if(k&1)ans.mul(x);}int f[MAXN];inline int F(int x){return x>0?f[x]:1;}int main(){ while(~scanf("%d%d%d%d%d",&n,&A,&B,&p,&q)){ A%=Mod,B%=Mod; rep(i,1,q<<1|1){ f[i]=(F(i-p)*A+F(i-q)*B)%Mod; } if(n<=(q<<1)){ printf("%d\n",f[n]); continue; } len=q<<1|1,bit=1; while(1<<bit<len)++bit; len=1<<bit; init_fft(len,bit); qpow(a,b,n-q); int ans=0; for(int i=q-1;~i;--i){ ans=(ans+a.a[i]*f[i+q])%Mod; } printf("%d\n",ans); } return 0;} |
