2016-09-26

HDU 5895 Mathematician QSC【矩阵快速幂+数论】

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=5896

题意：

定义$f(n)=2\times f(n-1)+f(n-2),f(1)=1,f(0)=0，g(n)=\sum_{i=0}^nf(i)^2$，给定$n,y,x,s$，求$x^{g(n*y)}\% (s+1)$。
数据范围：$n和x最大为8位十进制数，y最大4位十进制数，1 \le s \le 100000000$

分析：

很显然的矩阵快速幂，但是需要进行一下转化和处理。
首先我们试着简化这个$g(n)$的求和形式。
已知$f(n+1)=2 \times f(n)+f(n-1)$
即$f(n) = \frac{f(n + 1) - f(n-1)}{2} $
那么$f(n)^2 = \frac {f(n+1)f(n) - f(n-1)f(n)}{2}$
对于$n$，将前面所有项加起来，最终得到
$g(n) = \sum_{i=0}^nf[i]^2= \frac{f(n+1)f(n)}{2}$
我们用矩阵快速幂求出$f(n)和f(n+1)$即可。
求出$g(n \times y)$之后，这个值会非常大，我们需要降幂，降幂的套路比较常见了，$s+1$和$x$不一定互质，利用公式$a^b\%p=(a\%p)^{\phi(p) + b \% \phi(p)} \% p$
然后我们发现最后答案要除以个$2$，而又不能保证存在逆元，所以这里学到了一个新姿势= =
$\frac{a}{b}\%p = \frac{a \% (b \times p) }{b}$ ，当然前提是$b \times p$不会爆。
然后就写完啦。

推公式每次都很难过= = 一推就容易跑偏= =

代码：

/*************************************************************************
    > File Name: 5895.cpp
    > Author: jiangyuzhu
    > Mail: 834138558@qq.com 
    > Created Time: 2016/9/26 10:10:35
 ************************************************************************/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<queue>
#include<cstring>
#include<stack>
#include<vector>
#include<algorithm>
#include<map>
#include<cmath>
using namespace std;
typedef long long ll;
const int maxn = 105;
ll mod;
int eurler_phi(int n)
{
    int res = n;
    for(int i = 2; i * i <= n; i++){
        if(n % i == 0){
            res = res / i * (i - 1);
            while(n % i == 0) n /= i;
        }
    }
    if(n != 1) res = res / n * (n - 1);
    return res;
}
const int N = 5;
struct Matrix
{
    int row,cal;
    long long m[N][N];
};
Matrix init(Matrix a, long long t)
{
    for(int i = 0; i < a.row; i++)
        for(int j = 0; j < a.cal; j++)
            a.m[i][j] = t;
    return a;
}
Matrix operator * (const Matrix& a, const Matrix& b)
{
    Matrix ans;
    ans.row = a.row, ans.cal = b.cal;
    ans = init(ans,0);
    for(int i = 0; i < a.row; i++)
        for(int j = 0; j < b.cal; j++)
            for(int k = 0; k < a.cal; k++)
                ans.m[i][j] = (ans.m[i][j] + a.m[i][k] * b.m[k][j]) % mod;
    return ans;
}
void Matrix_pow(ll k, Matrix& Q)
{
	if(k <= 0) return;
    Matrix I;
    I.row = 2, I.cal = 2;
    I = init(I, 0);
	I.m[0][0] = 2; I.m[0][1] = 1;
	I.m[1][0] = 1; I.m[1][1] = 0;
	for(; k; k >>= 1, I = I * I){
		if(k & 1) Q = I * Q;
	}
}
char str[10 + 5];
int read()
{
	int x = 0;
	scanf("%s", str);
	int len = strlen(str);
	bool flag = true;
	for(int i = 0; i < len; ++i){
		if(str[i] != '0') flag = false;
		if(!flag) x = x * 10 + str[i] - '0';
	}	
	return x;
}
int quick_pow(int a, int b, int mod)
{
	int ans = 1;
	for(;a; a >>= 1, b = b * 1ll * b % mod){
		if(a & 1) ans = ans * 1ll * b % mod;
	}
	return ans;
}
int main(void)
{
	int T;scanf("%d", &T);
	while(T--){
		int	n = read(); 
		int y = read(); 
		int x = read(); 
		int s = read();
		int phi = eurler_phi(s + 1);
		mod = 2 * phi;
		Matrix A;
		A.row = 2, A.cal = 1;
		A = init(A, 0);
		A.m[0][0] = 1, A.m[1][0] = 0;
		Matrix_pow(n * 1ll * y, A);
		int ans = (A.m[0][0] * 1ll * A.m[1][0]) % mod / 2ll;
		int res = quick_pow(ans + phi, x, s + 1);
		printf("%d\n", res);
	}
    return 0;
}