【BZOJ3328】PYXFIB（数学）

什么都不会的数学蒻菜瑟瑟发抖……Orz橙子（和兔子）

题目：

BZOJ3328

分析：

橙子给我安利的数学题……（然后我就看着他因为矩阵乘法多模了一次卡了一天常数qwq表示同情）

先考虑一个子问题：求下面这个式子的值：

[sum_{i=0}^n C_n^i F_i ]

首先我们知道(F_i=left[egin{matrix}1 & 1\1 & 0end{matrix} ight]^i)。设那个矩阵为(A)，即(F_i=A^i)。（注意这题斐波那契数列下标从(0)开始，所以(F_2=2)。）

（不知道？你把(left[egin{matrix}F_i & F_{i-1}\F_{i-2} & 0end{matrix} ight])乘一下(A)试试。一开始左下方的值并不影响计算结果。）

然后(sumlimits_{i=0}^n C_n^i A^i)这个东西好像长得挺像二项式定理的。事实上，如果两个矩阵相乘可以交换（即(AB=BA)），则这两个矩阵之和的幂可以用二项式定理展开。因此，上面那个式子就是((A+I)^n)，其中(I)是单位矩阵（任何矩阵与单位矩阵相乘都是可交换的）。

（于是我推到这里忘了原题是什么兴高采烈地去给橙子说我会了

那么问题来了，如何“筛”掉当(i)不是(k)的倍数的时候对答案的贡献呢？介绍一个神奇的东西：单位根。

首先先介绍原根。若在模(p)（这里只讨论(p)是质数的情况）意义下，(a(0leq a<p))是一个满足对于任意(i(1leq i < p - 1))都满足(a^i eq 1)的整数，则(a)是(p)的原根。（根据欧拉定理，(a^i=a^{i mod phi(p)})，所以不讨论指数不小于(phi(p))即(p-1)的情况。）原根满足对于任意小于(phi(p))的非负整数(i)，(a^i)在模(p)意义下互不相同。以下的“原根”都特指最小原根，记为(g)。求法？从(2)开始暴力枚举(a)，把(i=p')（其中(p')是(p-1)的质因数）挨个试一遍看有没有(a^iequiv 1mod p)的情况就完了。（并不知道为什么是对的并且跑得飞快）

模(p)的(k)次单位根(w_k=g^{frac{p-1}{k}})，所以(w_k^k equiv 1mod p)。

在模(p)的意义下，(p)的(k)次单位根(w_k)具有如下性质：

[sum_{i=0}^{k-1}w_k^{ij}=egin{cases}k (j|k)\0 (otherwise)end{cases} ]

证明？当(j|k)，(ij)全是(k)的倍数，所以(w_k^{ij}=1)，全部加起来就是(k)。

否则，大力等比数列求和公式(w_k^0cdot frac{1-w_k^{jk}}{1-w_k^j})。上面(w_k^{jk}=1)，所以和是(0)。这个有点像FFT中的单位根，只是一个在单位圆上，一个在模域上，都是循环的（注意模域中指数的循环节是(phi(p)=p-1)）。

基于这个性质，我们达到了筛掉(i)不是(k)的倍数的情况的目的。那么要求的那个式子就变成了：

[sum_{i=0}^n C_n^i A^icdotfrac{1}{k}sum_{j=0}^{k-1}w_k^{ij} ]

发现后面那一堆和(i)没什么关系，提到前面来（相当于把(w_k^{ij})分配进去）：

[frac{1}{k}sum_{j=0}^{k-1}sum_{i=0}^n C_n^i A^iw_k^{ij} ]

好像挺像二项式定理的，只是(w_k^j)的指数比较丑。但我们可以这样……

[frac{1}{k}sum_{j=0}^{k-1}sum_{i=0}^n C_n^i A^iw_k^{(-i)cdot(-j)} ]

然后设(x_j=w_k^{-j})，得到：

[frac{1}{k}sum_{j=0}^{k-1}x_j^{-n}sum_{i=0}^n C_n^i A^ix_j^{n-i} ]

于是可以上矩阵二项式定理了（注意因为(x_j)是一个整数，所以要乘上单位矩阵(I)才能与(A)相加）。

[frac{1}{k}sum_{j=0}^{k-1}x_j^{-n}(x_jI+A)^n ]

照着这个就可以直接算了qwq

代码：

矩阵乘法的时候一定不要模两次！一定不要模两次！一定不要模两次！否则你会像橙子一样卡一天常数。

以及记着用上面提到过的欧拉定理把所有指数搞成非负的。

#include <cstdio>
#include <cstring>
#include <cctype>
#include <algorithm>
#define _ 0
using namespace std;

namespace zyt
{
	template<typename T>
	inline bool read(T &x)
	{
		char c;
		bool f = false;
		x = 0;
		do
			c = getchar();
		while (c != EOF && c != '-' && !isdigit(c));
		if (c == EOF)
			return false;
		if (c == '-')
			f = true, c = getchar();
		do
			x = x * 10 + c - '0', c = getchar();
		while(isdigit(c));
		if (f)
			x = -x;
		return true;
	}
	template<typename T>
	inline void write(T x)
	{
		static char buf[20];
		char *pos = buf;
		if (x < 0)
			putchar('-'), x = -x;
		do
			*pos++ = x % 10 + '0';
		while (x /= 10);
		while (pos > buf)
			putchar(*--pos);
	}
	typedef long long ll;
	ll n;
	int k, p;
	struct Matrix
	{
		int n, m, data[2][2];
		Matrix(const int _n = 0, const int _m = 0)
			: n(_n), m(_m)
		{
			for (int i = 0; i < n; i++)
				memset(data[i], 0, sizeof(int[m]));	
		}
		~Matrix() {}
		Matrix operator + (const Matrix &b) const
		{
			Matrix ans(n, m);
			for (int i = 0; i < n; i++)
				for (int j = 0; j < m; j++)
					ans[i][j] = (data[i][j] + b[i][j]) % p;
			return ans;
		}
		Matrix operator * (const Matrix &b) const
		{
			Matrix ans(n, b.m);
			for (int i = 0; i < n; i++)
				for (int k = 0; k < m; k++)
					for (int j = 0; j < b.m; j++)
						ans[i][j] = (ans[i][j] + 
							(ll)data[i][k] * b[k][j] % p) % p;
			return ans;
		}
		Matrix operator * (const ll &x) const
		{
			Matrix ans(n, m);
			for (int i = 0; i < n; i++)
				for (int j = 0; j < m; j++)
					ans[i][j] = ((ll)data[i][j] * x) % p;	
			return ans;
		}
		const int *operator [] (const int a) const
		{
			return data[a];	
		}
		int *operator [] (const int a)
		{
			return data[a];	
		}
	};
	inline Matrix get_identity(const int n)
	{
		Matrix ans(n, n);
		for (int i = 0; i < n; i++)
			ans[i][i] = 1;
		return ans;	
	}
	inline Matrix power(Matrix a, ll b)
	{
		Matrix ans = get_identity(a.n);
		while (b)
		{
			if (b & 1)
				ans = ans * a;
			a = a * a;
			b >>= 1;
		}
		return ans;
	}
	inline int power(int a, ll b)
	{
		int ans = 1;
		while (b)
		{
			if (b & 1)
				ans = (ll)ans * a % p;
			a = (ll)a * a % p;
			b >>= 1;
		}
		return ans;
	}
	inline int inv(const int a)
	{
		return power(a, p - 2);	
	}
	pair<int, int> prime[20];
	int cnt;
	inline void get_prime(int n)
	{
		cnt = 0;
		for (int i = 2; i * i <= n; i++)
		{
			if (n % i == 0)
				prime[cnt++] = make_pair(i, 0);
			while (n % i == 0)
				++prime[cnt - 1].second, n /= i;
		}
		if (n != 1)
			prime[cnt++] = make_pair(n, 1);
	}
	inline int get_g(const int p)
	{
		get_prime(p - 1);
		for (int i = 2; i < p; i++)
		{
			bool flag = true;
			for (int j = 0; j < cnt && flag; j++)
				if (power(i, (p - 1) / prime[j].first) == 1)
					flag = false;
			if (flag)
				return i;
		}
	}
	int work()
	{
		int T;
		read(T);
		Matrix I = get_identity(2);
		Matrix A(2, 2);
		A[0][0] = A[0][1] = A[1][0] = 1;
		while (T--)
		{
			read(n), read(k), read(p);
			int omega = power(get_g(p), (p - 1) / k), ans = 0;
			int tmp = power(omega, p - k - 1);
			for (int i = 0; i > -k; i--)
			{
				tmp = (ll)tmp * omega % p;
				ans = (ans + (ll)power(tmp, -n) * power(I * tmp + A, n)[0][0] % p) % p;
			}
			write((ll)ans * inv(k) % p);
			putchar('
');
		}
		return (0^_^0);	
	}
}
int main()
{
	freopen("3328.in", "r", stdin);
	freopen("3328.out", "w", stdout);
	return zyt::work();	
}