bzoj5104: Fib数列

Description

Fib数列为1,1,2,3,5,8...

求在Mod10^9+9的意义下，数字N在Fib数列中出现在哪个位置

无解输出-1

Input

一行，一个数字N，N < = 10^9+9

$r_1=frac{1+sqrt 5}{2}\ r_2=frac{1-sqrt 5}{2}=-frac{1}{r_1}\ N=Fib_x=r_1^x-r_2^x\ N^2=r_1^{2x}+r_2^{2x}-2(-1)^x\ ±(r_1^x+r_2^x)=sqrt{N^2+4(-1)^x}\ frac{N±sqrt{N^2+4(-1)^x}}{2}=r_1^x,-r_2^x\ x=min(log_{r_1}(frac{N±sqrt{N^2+4(-1)^x}}{2}),log_{r_2}(-frac{N±sqrt{N^2+4(-1)^x}}{2}))\ 枚举x的奇偶，利用离散对数计算答案$

#include<cstdio>
typedef long long i64;
const int P=1e9+9,g=13,sqrt5=383008016,I2=(P+1)/2,B=60000;
const int r1=(P+1+sqrt5)/2,r2=(P+1-sqrt5)/2;
const int lr1=133086171,lr2=366913833;
inline int mul(int a,int b){return (i64)a*b%P;}
inline void muls(int&a,int b){a=mul(a,b);}
inline int fix(int x){return x+(x>>31&P);}
int pw(int a,int n){
    int v=1;
    for(;n;n>>=1,muls(a,a))if(n&1)muls(v,a);
    return v;
}
void exgcd(int a,int b,int&x,int&y){
    if(!b){x=1,y=0;return;}
    exgcd(b,a%b,y,x);
    y-=a/b*x;
}
int inv(int a,int b){
    int x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int solve(int a,int b,int c,int tp){
    if(!a)return b?-1:2-tp;
    int g=gcd(a,c);
    if(b%g)return -1;
    a/=g,b/=g,c/=g;
    i64 t=(i64)b*inv(a,c)%c;
    if(!t)t=c;
    if(t%2!=tp)t+=c;
    return t%2!=tp?-1:t;
}
int h[131072][2];
int&at(int x){
    int w=x&131071;
    while(h[w][1]){
        if(h[w][0]==x)return h[w][1];
        w=(w+1)&131071;
    }
    h[w][0]=x;
    return h[w][1];
}
void pre(){
    int v=1;
    for(int i=0;i<B;++i){
        at(v)=i+1;
        muls(v,g);
    }
}
int log(int x){
    int t=pw(g,P-1-B);
    for(int i=0;;i+=B){
        int y=at(x);
        if(y)return y-1+i;
        muls(x,t);
    }
}
int sqrt(int x){
    int t=log(x);
    return t&1?-1:pw(g,t/2);
}
int ans=-1;
void upd(int v){
    if(~v&&(v<ans||ans==-1))ans=v;
}
void cal(int x,int tp){
    upd(solve(lr1,log(x),P-1,tp));
    upd(solve(lr2,log(fix(-x)),P-1,tp));
}
int main(){
    pre();
    int v;
    scanf("%d",&v);
    if(!v)return puts("-1"),0;
    muls(v,sqrt5);
    int v2=mul(v,v);
    int s1=sqrt(fix(v2-4)),s2=sqrt(fix(v2+4-P));
    if(~s1){
        cal(mul(fix(v+s1-P),I2),1);
        cal(mul(fix(v-s1),I2),1);
    }
    if(~s2){
        cal(mul(fix(v+s2-P),I2),0);
        cal(mul(fix(v-s2),I2),0);
    }
    return printf("%d
",ans),0;
}