可做题2（Code+12月网络赛）（扩欧+矩阵快速幂）

题目背景

“codeplus比赛的时候在做什么？有没有空？能来解决丢番图方程问题吗？”sublinekelzrip这样问qmqmqm。
当然，qmqmqm并不会丢番图方程问题，所以sublinekelzrip改为提出了另一个题目，现在请你帮助qmqmqm解决这个题目。

题目描述

这个问题是这样的:
若一个数列$a$满足条件$a_n=a_{n-1}+a_{n-2},n>3,而a_1,a_2$为任意实数，则我们称这个数列为广义斐波那契数列。
现在请你求出满足条件$a_1=i$，$a_2$为区间$[l,r]$中的整数，且$a_k\mathrm{m}\mathrm{o}\mathrm{d}p=m$的广义斐波那契数列有多少个。

输入格式

从标准输入读入数据。
本题包含多组数据，输入第一行包含一个正整数$T$，表示数据组数。对于每组数据：一行六个用空格隔开的整数$i,l,r,k,p,m$，意义如题目描述所示。

输出格式

输出到标准输出。
输出共$T$行，每行一个数表示该组数据的答案。

样例1输入

6
2 17 68 3 23 1
1 17 68 3 57 1
5 17 68 10 11 9
5 17 68 10 71 9
10 17 68 11 12 3
10 17 68 8 6 4

样例1输出

3
1
4
1
5
9

子任务

测试点	k	r	其他
1	<=100	<=100	无
2	<=10^5	<=10^5
3
4	<=10^{18}
5
6	<=10^5	<=10^{18}	p为质数
7
8	<=10^{18}
9			无
10

题解：

算法一

暴力枚举所有可能的$a_2$并递推判断。
复杂度$O(r\times k)$,预期得分10分。

算法二

$a_k$可以表示为$a_1$与$a_2$的线性组合。
使用递推计算出系数，并暴力枚举所有可能的$a_2$判断。
复杂度$O(r+k)$,预期得分30分。

算法三

暴力枚举所有可能的$a_2$并使用矩阵乘法判断。
复杂度$O(r\times \log (k))$,预期得分50分。

算法四

与算法二类似，使用递推计算出系数，此时可以发现可能的$a_2$满足一个同余方程，使用扩展欧几里得或逆元求解即可。
可以通过测试点6,7。

算法五

将算法四中计算系数的方式改为矩阵乘法，可以通过测试点8。

算法六

由于$p$可能不是质数，所以需要判断不互质的情况，然后使用扩展欧几里得或欧拉定理求解同余方程。可以通过所有数据。

考场上果断算法三
后来虚泽口hu出了正解

如果知道了首项和第二项
我们就可以用矩阵快速幂求出第k项

相当于我们已知a[k]%p，a[1]，fib(k-2)，fib(k-1)，
令a=fib(k-1)，b=m-i*fib(k-2)，x=a[2]，

则方程转化为：a*x ≡ b(%p)，求解x=[l,r]的整数解

运用exgcd求解过程：

• 转化为不定方程：ax-py=b
• g=gcd(a,p)，若b%p≠0则无解（输出0）
• a/=g，p/=g，b/=g
• 求解a’x-p’y=1即exgcd(a,p)
• 得到最小非负整数解x=x*b（如果x<0，我们需要通过+p%p处理一下）
• 计算在[l,r]之间的解：cal(r,x)-cal(l-1,x)

ll cal(ll r,ll x)
{
	if (r<x) return 0;
	return (r-x)/p+1;    //每隔p个就有一解
}

tip

i读入的时候就需要%p

在solve的时候：p —> p/gcd(a,p)

不知道这是什么鬼畜原理，
可能是因为求解是建立在：ax+py=gcd(a,p) 的基础上
之后的解也是每隔$p/gcd(a,p)$就有一个

#include<cstdio>
#include<cstring>
#include<iostream>
#define ll long long

using namespace std;

ll p;
struct node{
	ll H[3][3];
	node operator *(const node &a) const 
	{
		node ans;
		for (int i=1;i<=2;i++)
			for (int j=1;j<=2;j++){
			    ans.H[i][j]=0;
			    for (int k=1;k<=2;k++)
			        ans.H[i][j]=(ans.H[i][j]+H[i][k]*a.H[k][j])%p;
		}
		return ans;
	}
	void clear()
	{
		memset(H,0,sizeof(H));
	}
	node KSM(ll p)
	{
		node ans=(* this);
		node a=(* this);
		p--;
		while (p)
		{
			if (p&1) ans=ans*a;
			p>>=1;
			a=a*a;
		}
		return ans;
	}
};
node H;

ll i,l,r,k,m,x,y;

ll gcd(ll a,ll b)
{
	ll r=a%b;
	while (r)
	{
		a=b; b=r;
		r=a%b;
	}
	return b;
}

void exgcd(ll a,ll b)
{
	if (b==0)
	{
		x=1; y=0;
		return;
	}
	exgcd(b,a%b);
	ll t=x;
	x=y;
	y=t-a/b*y;
}

ll cal(ll r,ll x)
{
	if (r<x) return 0;
	return (r-x)/p+1;
}

ll solve(ll a,ll b)
{
	ll c=gcd(a,p);
	if (b%c!=0) return 0;
	a/=c; b/=c; p/=c;               //p/=c
	exgcd(a,p);
	x=(x*b+p)%p;                    //解 
	while (x<0) x=(x+p)%p;
	return (ll)(cal(r,x)-cal(l-1,x));
}

int main()
{
	H.H[1][1]=1; H.H[1][2]=1;
	H.H[2][1]=1; H.H[2][2]=0;
	int T;
	scanf("%d",&T);
	while (T--)
	{
		scanf("%lld%lld%lld%lld%lld%lld",&i,&l,&r,&k,&p,&m);
		i%=p;
		node ans=H.KSM(k-2);
		ll f1=ans.H[1][1]%p;               //fib(k-1) 
		ll f2=((m-ans.H[2][1]%p*i%p)%p+p)%p;   //fib(k-2)
		//f1*x ≡f2(%p)
		printf("%lld\n",solve(f1,f2));
	}
	return 0;
}