【BZOJ3328】PYXFIB（数学）

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

题目：

BZOJ3328

分析：

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

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

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

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

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

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

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

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

首先先介绍原根。若在模\(p\)（这里只讨论\(p\)是质数的情况）意义下，\(a(0\leq a<p)\)是一个满足对于任意\(i(1\leq i < p - 1)\)都满足\(a^i\neq 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^i\equiv 1\mod p\)的情况就完了。（并不知道为什么是对的并且跑得飞快）

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

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

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

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

否则，大力等比数列求和公式\(w_k^0\cdot \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^i\cdot\frac{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('\n');
		}
		return (0^_^0);	
	}
}
int main()
{
	freopen("3328.in", "r", stdin);
	freopen("3328.out", "w", stdout);
	return zyt::work();	
}