【2022省选模拟】GCD——数论、矩阵加速

最新推荐文章于 2024-03-29 11:56:01 发布

偶耶XJX

最新推荐文章于 2024-03-29 11:56:01 发布

阅读量350

点赞数

分类专栏：信息竞赛解题文章标签：算法数论 GCD 矩阵加速

本文链接：https://blog.csdn.net/weixin_43960287/article/details/124107025

版权

信息竞赛解题专栏收录该内容

195 篇文章 5 订阅

订阅专栏

找不到原题链接

题目描述

在这里插入图片描述

题解

又要求 $\gcd$ ，又要取模，就只能用好 $\gcd$ 的性质了。这道题的核心思想是利用 $\gcd(a,b)=\gcd(a\bmod b,b)=\gcd(a,b\bmod a)$ ，那么我们在求两个大数的 $\gcd$ 时，可以试着把其中一个数变小以直接计算。

首先如果向量 $(a, b)$ 和 $(c, d)$ 线性相关，那么最大公因数就是两边中较小的那个，否则我们可以用辗转相除消元的方法消掉 $d$ ，把问题转化为求 $Ans=\gcd(af_n+bf_{n+1},cf_n)$ 。

考虑先求 $G=\gcd(af_n+bf_{n+1},f_n)$ ，然后再求 $Ans=G\cdot \gcd(\frac{af_n+bf_{n+1}}{G},c)$ 。设 $g_n=\gcd(f_n,f_{n+1})$ ，那么有
$G=\gcd(af_n+bf_{n+1},f_n)=\gcd(bf_{n+1},f_n)=g_n\cdot \gcd(b,\frac{f_n}{g_n})$ 通过打表发现 $g_n=3^{\lfloor\frac{n}{2}\rfloor}$ ，所以我们可以用矩阵加速来求 $\frac{f_n}{g_n}$ 。这个除法只需要稍稍修改一下每两步的转移矩阵即可，因为转移矩阵平方的每个数都是3的倍数。

我们求得 $\frac{f_n}{g_n}$ 在取模 $b$ 意义下的值过后，就可以直接算 $\gcd(b,\frac{f_n}{g_n})$ 。

设我们求得的 $\gcd(b,\frac{f_n}{g_n})=k$ ，现在来考虑答案的式子怎么算：
$Ans=G\cdot \gcd(\frac{af_n+bf_{n+1}}{G},c)\\ =g_nk\cdot \gcd(\frac{a\frac{f_n}{g_n}+b\frac{f_{n+1}}{g_n}}{k},c)$
我们需要求出 $\frac{a\frac{f_n}{g_n}+b\frac{f_{n+1}}{g_n}}{k}$ 在取模 $c$ 意义下的值以方便直接计算。注意到 $\frac{f_n}{g_n},b$ 都保证是 $k$ 的倍数，而 $\forall b|a,\frac{a}{b}\bmod c=\frac{a\bmod bc}{b}$ ，所以我们用矩阵求出上面部分 $\bmod kc$ 的值，然后直接除以 $k$ 即可。

复杂度是 $T\log n$ ，但是比较卡常。

代码

#include<bits/stdc++.h>//JZM yyds!!
#define ll long long
#define uns unsigned
#define END putchar('\n')
#define fi first
#define se second
#define IF (it->fi)
#define IS (it->se)
#define lowbit(x) ((x)&-(x))
using namespace std;
const int MAXN=-1;
const ll INF=1e18;
inline ll read(){
	ll x=0;bool f=1;char s=getchar();
	while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
	while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
	return f?x:-x;
}
int ptf[50],lpt;
inline void print(ll x,char c='\n'){
	if(x<0)putchar('-'),x=-x;
	ptf[lpt=1]=x%10;
	while(x>9)x/=10,ptf[++lpt]=x%10;
	while(lpt)putchar(ptf[lpt--]^48);
	if(c>0)putchar(c);
}

inline ll gcd(ll a,ll b){return !b?a:gcd(b,a%b);}
ll n,MOD,a,b,c,d;
inline ll ksm(ll a,ll b,ll mo){
	ll res=1;
	for(;b;b>>=1,a=a*a%mo)if(b&1)res=res*a%mo;
	return res;
} 
ll MD;
struct matrix{
	ll a[2][2];int n,m;matrix(){}
	matrix(int N,int M){memset(a,0,sizeof(a)),n=N,m=M;}
	matrix(ll A,ll B,ll C,ll D){
		n=m=2,a[0][0]=A,a[0][1]=B,a[1][0]=C,a[1][1]=D;
	}
	matrix operator*(const matrix&b)const{
		matrix res=matrix(n,b.m);
		for(int i=0;i<n;i++)
			for(int k=0;k<m;k++)if(a[i][k])
				for(int j=0;j<b.m;j++)
					(res.a[i][j]+=a[i][k]*b.a[k][j])%=MD;
		return res;
	}
};
inline matrix mpow(matrix a,ll b){
	matrix res=matrix(a.n,a.n);
	for(int i=0;i<a.n;i++)res.a[i][i]=1;
	for(;b;b>>=1,a=a*a)if(b&1)res=res*a;
	return res;
}
inline ll fn(ll n,ll mod){
	matrix a=matrix(1,1,0,0),b=matrix(0,12,1,9);
	a.n=1,MD=mod,a=a*mpow(b,n);
	return a.a[0][0];
}
ll mor;
inline ll fn_gn(ll n,ll mod){
	MD=1e9+7;
	matrix a=matrix(1,1,0,0),b=matrix(0,12,1,9),c=b*b;
	for(int i=0;i<2;i++)for(int j=0;j<2;j++)c.a[i][j]/=3;
	a.n=1,MD=mod,a=a*mpow(c,n>>1)*mpow(b,n&1),mor=(a*b).a[0][0];
	return a.a[0][0];
}
int main()
{
	freopen("c.in","r",stdin);
	freopen("c.out","w",stdout);
	for(int TMA=read();TMA--;){
		n=read(),MOD=read(),a=read(),b=read(),c=read(),d=read();
		while(d)a-=b/d*c,b%=d,swap(a,c),swap(b,d);
		if(!c){
			print((fn(n,MOD)*(a%MOD+MOD)+fn(n+1,MOD)*(b%MOD+MOD))%MOD);
			continue;
		}if(c<0)c=-c;
		if(a<0)a+=((-a)/c+1)*c;
		if(!b){
			print(fn(n,MOD)*gcd(a,c)%MOD);
			continue;
		}
		ll k=gcd(b,fn_gn(n,b)),G=ksm(3,n>>1,MOD)*k%MOD;
		ll A=fn_gn(n,c*k)/k%c*a%c,B=b/k%c*mor%c;
		print(gcd(c,A+B)*G%MOD);
	}
	return 0;
}