题解：P4031 [Code+#2] 可做题2

题面

若一个数列\(a\)满足条件\(a_n=a_{n-1}+a_{n-2}\),\(n \geq 3\),而\(a_1,a_2\)为任意实数，则我们称这个数列为广义斐波那契数列。

现在请你求出满足条件\(a_1=i\)，\(a_2\)为区间\([l,r]\)中的整数，且\(a_k mod p = m\)的广义斐波那契数列有多少个。

题意

给定斐波那契数列 \(A\)，\(I,l,r,k,p,m\)

对于 \(i>=3\) , \(A_i=A_{i-1}+A_{i-2}\)

定义 \(A_1=I,A_2 \in [l,r]\) 求出有多少个 合法\(A_2\)使得 $A_k mod p=m $

Solve

特殊情况

当 \(k=1\) 时。

如果 \(i modp=m\), 则答案为 \(r-l+1\)，否则为 \(0\)

当 \(k=2\) 时

统计 \(a_2 modp=m\)， 且 \(r\leq a_2\leq r\) 的数量

k >= 3 解法

由于 \(A_2\) 是不定值,\(A_1\)是定值。

所以对于任意一个 \(i>3\), 我们可以表示为 \(v_i+t_i*a_2\)

那么我们要计算出 \(A_k\)

考虑矩阵加速 。

我们先规范一下矩阵的书写方式。

对于矩阵 \(n*m\) 的矩阵 \(Matrix\)

\[Matrix= \begin{bmatrix} a_{1,1} & ... & a_{1,n} \\ ... & ... & ... \\ a_{m,1} & ... & a_{m,n} \end{bmatrix} \]

我们可以使用矩阵来表示一个状态二元组。

\(\{v_{i-1}+t_{i-1}*a_2,v_i+t_i*a_2\}\)

\[a_i= \begin{bmatrix} v_{i-1} & v_{i} \\ t_{i-1} & t_{i}\\ \end{bmatrix} \]

我们定义一个矩阵\(mul\)，使得对于任意一个 \(a_i\), \(a_i*mul=a_{i+1}\)

这个你试试就出来了。

\[mul= \begin{bmatrix} 0 & 1 \\ 1 & 1\\ \end{bmatrix} \]

我们把初始状态 \(a_2\) 使用矩阵描述。

\(A_1=i+0*a_2\)

\(A_2=0+1*a_2\)

\[a_2= \begin{bmatrix} i & 0 \\ 0 & 1\\ \end{bmatrix} \]

那么计算\(A_k\)就需要计算矩阵 \(a_K\)

计算矩阵 \(a_K\) 就要将 \(a_2\) 乘 \(mul\)， \(k-2\) 次。

即 \(a_2*mul^{k-2}\)

这个可以使用矩阵快速幂加速 。

复杂度 \(log_2k\)

现在得到了 \(A_k=v_k+t_k*a_2\) 其中 \(v_k,t_k\)是已知的。

现在要求

\[(v_k+t_k*a_2) mod p=m \]

即存在一个 \(b \in Z\)

\[t_k*a_2+t_p*p=m-v_k \]

求 \(a_2\) \(t_p\) 的整数解有多少个。

首先判断是否有解。

如果 \(gcd(t_k,p) \not| (m-v_k)\)

那么输出 \(0\) 表述无解。

否则解不定方程。

先求出一组特解 \(t_k*a_0+t_0*p=gcd(t_k,p)\)

然后求出特解 \(a_2'=a_0*\frac{m-v_k}{gcd(t_k,p)},t_p'=t_0*\frac{m-v_k}{gcd(t_k,p)}\)

然后对于任意的 \(k \in Z\) ,所有解都可以表示为\({a_2'+k*\frac{p}{gcd(p,t_k)},t_0-k*\frac{t_k}{gcd(p,t_k)}}\) 都是一组解。

我们只关心 \(a_2\) 的值，原来的问题转化为有多少个 \(k\), 使得 \(l\leq a_2'+k*\frac{p}{gcd(p,t_k)} \leq r\) 且 \((a_2'+k*\frac{p}{gcd(p,t_k)}) \in Z\)

推式子 。

\[l\leq a_2'+k*\frac{p}{gcd(p,t_k)} \leq r \]

\[l-a_2'\leq k*\frac{p}{gcd(p,t_k)} \leq r-a_2' \]

其实就是求有多少个合法的 \(k\)。

那么合法的 \(k\) 的数量为 。

\[\lfloor\frac{r-a_2'}{\frac{p}{gcd(p,t_k)}}\rfloor-\lceil\frac{l-a_2'}{\frac{p}{gcd(p,t_k)}}\rceil +1 \]

注意特判\(k\)无解。

#include <bits/stdc++.h>

typedef long long ll;

using std::cin;
using std::cout;

ll m1[3][3],base[3][3]; //分别是代表A1,A2的初始矩阵  和用于矩阵快速幂的加速矩阵
ll t,I,l,r,k,p,m;

void init() {
    m1[1][1]=I;
    m1[1][2]=0;
    m1[2][1]=0;
    m1[2][2]=1;

    base[1][1]=0;
    base[1][2]=1;
    base[2][1]=1;
    base[2][2]=1;
    return;
}
// l<-r
void copy(ll l[][3],ll r[][3]) {
    l[1][1]=r[1][1];
    l[1][2]=r[1][2];
    l[2][1]=r[2][1];
    l[2][2]=r[2][2];
    return;
}

//c=a*b
void mul(ll a[][3], ll b[][3], ll c[][3]) {
    ll temp[3][3];
    memset(temp,0,sizeof(temp));

    for(ll j=1;j<=2;++j) { 
        for(ll i=1;i<=2;++i) { 
            ll sum=0;
            for(ll k=1;k<=2;++k) { 
                sum+=(a[j][k]%p)*(b[k][i]%p)%p;
                sum%=p;
            }
            temp[j][i]=sum;
        }
    }

    copy(c,temp);
    return;
}

//n次方mod p
void fmpow(ll n) {
    ll res[3][3];
    memset(res,0,sizeof(res));
    res[1][1]=1;
    res[2][2]=1;

    while(n>0) {
        if(n&1) {
            mul(res,base,res);
        }
        mul(base,base,base);
        n>>=1;
    }
    copy(base,res);
    return;
}

ll gcd(ll a,ll b) {
    if(!b) {
        return a;
    }
    return gcd(b, a%b);
}

//ax+by=gcd(a,b)
ll exgcd(ll a, ll b, ll &x, ll &y) {
    if (!b) {
        x = 1;
        y = 0;
        return a;
    }
    ll d = exgcd(b, a % b, x, y);
    ll t = x;
    x = y;
    y = t - (a / b) * y;
    return d;
}

void solve() {
    cin>>I>>l>>r>>k>>p>>m;
    if(k==1) {
        printf("%lld\n", (I%p==m)? r-l+1:0);
        return;
    }
    if(k==2) {
        r-=m,l-=m;
        printf("%lld\n", r/p-l/p+(l%p==0));
        return;
    }

    init();
    fmpow(k-2);
    mul(m1,base,m1);

    //现在A[k]=m1[1][2]+m1[2][2]*A[2];
    ll v=m1[1][2],t=m1[2][2],x,y,g=gcd(p,t);
    if((m-v)%g) {
        puts("0");
        return;
    }
    exgcd(t,p,x,y);
    x*=(m-v)/g;
    ll kmax=std::floor((long double)(r-x)/(p/g));
    ll kmin=std::ceil((long double)(l-x)/(p/g));
    ll ans=std::max(kmax-kmin+1,0);
    printf("%lld\n",ans);
}

int main() {
    freopen("input","r",stdin);

    cin>>t;

    while(t--) {
        solve();
    }

    return 0;
}