HDU 5895 Mathematician QSC （矩阵快速幂 + 逆元应用 + 指数循环节 + 欧拉函数）

大体题意：

题意很简单 计算表达式x^g[ny] % (s+1)

思路：

这个题弄了两天了，查了好多知识点才搞定 = =！

在看看g(n)的定义   g(n) = f(1)^2 + f(2)^2 + ... + f(n)^2;

因为

所以：

多写几个就可以看出累加法：

所以：

因此，gn是可以递推出来的！

公式是求x^g(n*y)% (s+1)  ，g(n*y)肯定会很大很大，因此肯定会用到了指数循环节！ 在来看看数据范围，S是不超过1亿的，我们打几个表看下，gn，在n=12 时 就远远超过1亿了，所以，先判断，如果n*y >= 12的话，肯定会大于s+1的欧拉函数，直接用指数循环节公式！

小于12 直接快速幂暴力求解了！

如果大于12的话：

先求出f[n]和f[n+1]    这个递推式很明显类似于斐波那契数列  构造矩阵  

直接取b^n  那么 b[0][0] = f[n+1],b[1][0] = f[n]，矩阵快速幂肯定会取模了！

那个公式有个除法，想到是逆元，但是除以2，s+1也有可能是偶数，因此 还不一定存在逆元，= =！

这里就用到了 一个逆元的应用：(这个式子借鉴的学长的博客！)

这样避免了除法取模的问题！很巧妙！

详细见代码吧（很丑陋= =）

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
ll n,y,x,s;
ll f[100],g[100];
struct Mar{
	ll a[2][2];
	void init(){
		memset(a,0,sizeof a);
	}
}unit,t1;
ll phi(ll n){
	ll m = (ll)sqrt(n+0.5);
	ll ans = n;
	for (int i = 2; i <= m; ++i){
		if (n % i == 0){
			ans = ans/i * (i-1);
		}
		while(n % i == 0)n /= i;
	}
	if (n > 1)ans = ans/n * (n-1);
	return ans;
}
Mar multi(Mar u,Mar v,ll mod){
	Mar t;
	t.init();
	t.a[0][0] = ((u.a[0][0]%mod * v.a[0][0] % mod) % mod + (u.a[0][1] %mod * v.a[1][0]%mod) % mod) % mod;
	t.a[0][1] = ((u.a[0][0]%mod * v.a[0][1]%mod) % mod + (u.a[0][1]%mod * v.a[1][1]%mod) % mod) % mod;
	t.a[1][0] = ((u.a[1][0]%mod * v.a[0][0]%mod) % mod + (u.a[1][1]%mod * v.a[1][0]%mod) % mod) % mod;
	t.a[1][1] = ((u.a[1][0]%mod * v.a[0][1]%mod) % mod + (u.a[1][1]%mod * v.a[1][1]%mod) % mod) % mod;
	return t;
}
ll my_pow(ll a,ll n,ll mod){
    ll ans = 1;
    while(n){
        if (n & 1)
            ans = (ans % mod * a % mod) % mod;
        n/=2;
        a = (a%mod*a%mod)%mod;
    }
    return ans;
 }
Mar my_pow2(Mar a,ll n,ll mod){
    Mar ans = unit;
    while(n){
        if (n & 1)
            ans = multi(ans,a,mod);
        n/=2;
        a = multi(a,a,mod);
    }
    return ans;
 }
int main(){
	unit.a[0][0] = unit.a[1][1] = 1;
	unit.a[0][1] = unit.a[1][0] = 0;
	t1.a[0][0] = 2;
	t1.a[0][1] = 1;
	t1.a[1][0] = 1;
	t1.a[1][1] = 0;
	int T;
	f[0] = 0;
	f[1] = 1;
	for (int i = 2; i <= 20; ++i)f[i] = f[i-2] + 2*f[i-1];
	g[0] = 0;
	for (int i = 1; i < 20; ++i)g[i] = (f[i]*f[i+1])/2;
	//for (int i = 0; i < 20; ++i)printf("%lld\n",g[i]);
	scanf("%d",&T);
	while(T--){
		scanf("%I64d%I64d%I64d%I64d",&n, &y, &x, &s);
		if (n*y <= 12){
			printf("%I64d\n",my_pow(x,g[n*y],s+1) % (s+1) );
			continue;
		}
		ll p = phi(s+1);
		Mar tmp = my_pow2(t1,n*y,2*p);
		ll temp = (tmp.a[0][1]%(2*p)*tmp.a[0][0]%(2*p))%(2*p);
		temp/=2;
		//printf("%lld\n",temp);
		temp += p;
		ll ans = my_pow(x,temp,s+1) % (s+1);
		printf("%I64d\n",ans);
	}
	return 0;
}