BZOJ5104 Fib数列(二次剩余+BSGS)
5在1e9+9下有二次剩余,那么fib的通项公式就有用了。
已知Fn,求n。注意到[(1+√5)/2]·[(1-√5)/2]=-1,于是换元,设t=[(1+√5)/2]n,原式变为√5·Fn=t-(-1)n·t-1。同乘t并移项,可得t2-√5·Fn·t-(-1)n=0。讨论n的奇偶性,BSGS求二次剩余大力解方程即可。用BSGS求二次剩余是非常简单的,求出其以原根为底的离散对数即可。
注意二次剩余有正负两解,但似乎代进去正根(即√gk=gk/2)就行了,不太明白。以及题目要求最小解,BSGS的时候注意顺序。还有BSGS不一定有解,我也不知道我在BSGS里面assert了半天是在干啥。调了一年惨炸了。
#include<iostream> #include<cstdio> #include<cmath> #include<cstdlib> #include<cstring> #include<algorithm> #include<map> #include<cassert> using namespace std; #define ll long long #define P 1000000009 char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<'0'||c>'9')) c=getchar();return c;} int gcd(int n,int m){return m==0?n:gcd(m,n%m);} int read() { int x=0,f=1;char c=getchar(); while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();} while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar(); return x*f; } int n,b,g,v,ans=P; map<int,int> f; int ksm(int a,int k) { int s=1; for (;k;k>>=1,a=1ll*a*a%P) if (k&1) s=1ll*s*a%P; return s; } int inv(int a){return ksm(a,P-2);} int BSGS(int g,int k) { f.clear(); int block=sqrt(P),t=ksm(g,block),x=1,ans=-1;g=inv(g); for (int i=0;i<block;i++) { if (f.find(1ll*x*k%P)==f.end()) f[1ll*x*k%P]=i; x=1ll*x*g%P; } x=1; for (int i=0;i<P;i+=block) { if (f.find(x)!=f.end()) {ans=f[x]+i;break;} x=1ll*x*t%P; } return ans; } int SQRT(int n) { int k=BSGS(g,n); if (k==-1||k&1) return -1; return ksm(g,k>>1); } void solve(int c,int op,int op2) { int delta=SQRT(((1ll*b*b-4ll*c)%P+P)%P); if (delta==-1) return; delta=(P+op2*delta)%P; int ans1=1ll*((delta-b)%P+P)%P*inv(2)%P,ans2=1ll*((-delta-b)%P+P)*inv(2)%P; ans1=BSGS(v,ans1),ans2=BSGS(v,ans2); if ((ans1&1)==op&&ans1>0) ans=min(ans,ans1); if ((ans2&1)==op&&ans2>0) ans=min(ans,ans2); } int fib(int n) { struct matrix { int n,a[2][2]; matrix operator *(const matrix&b) const { matrix c;c.n=n;memset(c.a,0,sizeof(c.a)); for (int i=0;i<n;i++) for (int j=0;j<2;j++) for (int k=0;k<2;k++) c.a[i][j]=(c.a[i][j]+1ll*a[i][k]*b.a[k][j])%P; return c; } }f,a; f.n=1;f.a[0][0]=0,f.a[0][1]=1; a.n=2;a.a[0][0]=0,a.a[0][1]=a.a[1][0]=a.a[1][1]=1; for (;n;n>>=1,a=a*a) if (n&1) f=f*a; return f.a[0][0]; } void work(int sqrt5) { b=(P-1ll*sqrt5*n%P)%P; v=1ll*(sqrt5+1)*inv(2)%P; solve(P-1,0,1),solve(1,1,1); //solve(P-1,0,-1),solve(1,1,-1); } int main() { #ifndef ONLINE_JUDGE freopen("bzoj5104.in","r",stdin); freopen("bzoj5104.out","w",stdout); const char LL[]="%I64d\n"; #else const char LL[]="%lld\n"; #endif /*for (int i=2;;i++) { bool flag=1; for (int j=2;j*j<P;j++) if ((P-1)%j==0) { if (ksm(i,j)==1) {flag=0;break;} if (ksm(i,(P-1)/j)==1) {flag=0;break;} } if (flag) {g=i;break;} }*/ g=13; n=read(); work(SQRT(5));//,work(P-SQRT(5)); if (ans==P) cout<<-1;else assert(fib(ans)==n),cout<<ans; return 0; }