HDU 3802 矩阵快速幂化简递推式子加一点点二次剩余知识

求$G(a,b,n,p) = (a^{\frac {p-1}{2}}+1)(b^{\frac{p-1}{2}}+1)[(\sqrt{a} + \sqrt{b})^{2F_n} + (\sqrt{a} - \sqrt{b})^{2F_n}] (mod p)$

左边可以看出是欧拉判定准则，那么只有当a,b其中一个满足是模p下的非二次剩余时G()为0。

右边的式子可以先把平方放进去，发现这个已经是通项公式了，那么$a+b+\sqrt{ab}$和$a+b-\sqrt{ab}$就是它的特征根了，反代回二阶特征方程，得到递推式的两个系数$2(a+b)$，$-(a-b)^2$，然后可以使用矩阵快速幂了，注意其项数是$F(n)$，指数循环节，欧拉降幂，其矩阵快速幂的过程中模p-1就行了。

/** @Date    : 2017-10-11 22:30:33

  * @FileName: HDU 3802 斐波那契特征方程解递推式 矩阵快速幂.cpp

  * @Platform: Windows

  * @Author  : Lweleth (SoungEarlf@gmail.com)

  * @Link    : https://github.com/

  * @Version : $Id$

  */

#include <bits/stdc++.h>

#define LL long long

#define PII pair

#define MP(x, y) make_pair((x),(y))

#define fi first

#define se second

#define PB(x) push_back((x))

#define MMG(x) memset((x), -1,sizeof(x))

#define MMF(x) memset((x),0,sizeof(x))

#define MMI(x) memset((x), INF, sizeof(x))

using namespace std;

const int INF = 0x3f3f3f3f;

const int N = 1e5+20;

const double eps = 1e-8;

LL mod;

LL fpow(LL a, LL n)

{

	LL res = 1LL;

	while(n)

	{

		if(n & 1LL) res = (res * a % mod + mod) % mod;

		a = (a * a % mod + mod) % mod;

		n >>= 1LL;

	}

	return res;

}

struct matrix

{

	LL mat[2][2];

	matrix(){MMF(mat);}

	void id(){

		mat[0][0] = mat[1][1] = 1LL;

	}

	matrix operator *(const matrix &b){

		matrix c;

		for(int i = 0; i < 2; i++)

			for(int j = 0; j < 2; j++)

				for(int k = 0; k < 2; k++)

				{

					c.mat[i][j] = (c.mat[i][j] + mat[i][k] * b.mat[k][j] % mod) % mod;

				}

		return c;

	}

	matrix operator ^(LL n){

		matrix res;

		matrix a = *this;

		res.id();

		while(n)

		{

			if(n & 1LL)

				res = res * a;

			a = a * a;

			n >>= 1LL;

		}

		return res;

	}

};

LL fib(LL n)

{

	mod--;

	matrix T;

	T.mat[0][0] = 1LL, T.mat[0][1] = 1LL;

	T.mat[1][0] = 1LL, T.mat[1][1] = 0LL;

	matrix tmp = (T^(n-1));

	LL ans = ((tmp.mat[0][0] + tmp.mat[1][0])%mod + mod) % mod;

	mod++;

	return ans;

}

LL fun(LL a, LL b, LL n)

{

	LL ans = (fpow(a, (mod - 1)>>1) + 1LL) * (fpow(b, (mod-1)>>1) + 1LL) % mod;//注意这里也要取模...

	if(ans == 0)

		return 0;

	if(n < 2)

		return (a + b) * 2LL % mod * ans % mod;

	LL e = fib(n);

	//cout << e << endl;

	if(e == 0)

		e = mod - 1;

	matrix A;

	A.mat[0][0] = (a + b) * 2LL % mod, 	   A.mat[0][1] = 1LL;

	A.mat[1][0] = (b - a) * (a - b) % mod, A.mat[1][1] = 0LL;

	A = A^(e - 1);

	ans = (ans * ((a + b) * 2LL % mod * A.mat[0][0] + A.mat[1][0] * 2LL % mod) % mod) % mod;

	while(ans < 0)

		ans += mod;

	return ans;

}

int main()

{

	int T;

	cin >> T;

	while(T--)

	{

		LL a, b, n;

		scanf("%lld%lld%lld%lld", &a, &b, &n, &mod);

		LL ans = fun(a, b, n);

		printf("%lld\n", ans);

	}

    return 0;

}

巴特西

HDU 3802 矩阵快速幂化简递推式子加一点点二次剩余知识

最新文章

热门文章

巴特西

HDU 3802 矩阵快速幂 化简递推式子 加一点点二次剩余知识

最新文章

热门文章

HDU 3802 矩阵快速幂化简递推式子加一点点二次剩余知识