【LibreOJ】#6257. 「CodePlus 2017 12 月赛」可做题2

【题意】数列满足a_n=a_n-1+a_n-2,n>=3。现在a1=i，a2=[l,r]，要求满足ak%p=m的整数a2有多少个。10^18。

【算法】数论（扩欧）+矩阵快速幂

【题解】定义fib(i)表示第 i 个斐波那契数，将数列an列项观察容易发现ak=a1*fib(k-2)+a2*fib(k-1)。fib(i)可以用矩阵快速幂迅速得解。

现在实际已知ak%p,a1,fib(k-2),fib(k-1)，令a=fib(k-1)，b=m-i*fib(k-2)，x=a2，则方程转化为：ax≡b(%p)，求解x=[l,r]的整数解。

运用扩展欧几里得定理求解即可，下面展示具体细节：

1.转化为不定方程,ax-py=b。

2.g=gcd(a,p)，若b%p≠0则无解（输出0）。

3.a/=g;p/=g;b/=g;

4.求解a'x-b'y=1即exgcd(a,p,x,y)，得到x0=x*b。

5.得到最小非负整数解x（这是为了防止x0>l）

6.计算在[l,r]之间的解，calc(r,x)=(r-b)%p+1。

注意：

1.读入 i 后取模。

2.先算原解x0，再扩展。

3.namespace中任何元素都不能重名。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
ll x1,l,r,k,m;
int p;
int MOD(int x){return x>=p?x-p:x;}
namespace fib{
    int ans[2][2],tmp[2][2],c[2][2];
    void mul(int a[2][2],int b[2][2]){
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++){
                c[i][j]=0;
                for(int k=0;k<2;k++)
                    c[i][j]=MOD(c[i][j]+1ll*a[i][k]*b[k][j]%p);
            }
        for(int i=0;i<2;i++)for(int j=0;j<2;j++)a[i][j]=c[i][j];
    }
    int f(ll x){
        x--;
        ans[0][0]=1;ans[1][0]=ans[0][1]=ans[1][1]=0;
        tmp[0][0]=tmp[0][1]=tmp[1][0]=1;tmp[1][1]=0;
        while(x){
            if(x&1)mul(ans,tmp);
            mul(tmp,tmp);
            x>>=1;
        }
        return ans[0][0];
    }
}
int M(ll x){return (x%p+p)%p;}
void exgcd(int a,int b,int& x,int& y){
    if(!b){x=1;y=0;}
    else{exgcd(b,a%b,y,x);y-=x*(a/b);}
}
int gcd(int a,int b){return b?gcd(b,a%b):a;}
ll calc(ll r,int x){return r<x?0:(r-x)/p+1;}
int main(){
    int T;scanf("%d",&T);
    while(T--){
        scanf("%lld%lld%lld%lld%d%lld",&x1,&l,&r,&k,&p,&m);//x1%p!!!
        int b=M(m-x1%p*fib::f(k-2)),a=fib::f(k-1),x,y;
        int g=gcd(a,p);
        if(b%g){printf("0
");continue;}
        a/=g;b/=g;p/=g;
        exgcd(a,p,x,y);
        x=M(1ll*x*b);//M(1ll*x*b)!!!
        printf("%lld
",calc(r,x)-calc(l-1,x));
    }
    return 0;
}