HDU 2481 Toy(Polya综合)_综合

题目链接：
HDU 2481 Toy
题意：
有外围 $n$ 个点围着中心一个点，中心有 $n$ 条边和外围相连，外围的相邻点有一条边，现在需要从中去除 $n$ 条边，使得剩下的 $n+1$ 个点依然保持两两连通，考虑旋转，问有多少种方案？方案数对 $M$ 取模。
数据范围： $3 \leq n \leq 10^9, 2\leq M \leq 10^9$ .
分析：
我觉得这种考验智商的题目真的不适合窝。。。
主要参考博客：
cx_love
将狼踩尽
zxb

任意取两个结点讨论 $a,b$ 。那么总数便是 $a,b$ 断开的种数与 $a,b$ 连在一起的种数的和。

f (n) 表 示 外 圈 有 n 个 结 点 时 ， 而 a, b 是 断 开 的 种 数 。

$f(n)表示外圈有n个结点时，而a,b是断开的种数。$

g (n) 表 示 外 圈 有 n 个 结 点 时 ， 而 a, b 是 连 在 一 起 的 种 数 。

$g(n)表示外圈有n个结点时，而a,b是连在一起的种数。$

如果 $a,b$ 之间是断开的，假设与 $a$ 直接相连的为 $k$ 个(加上 $a$ 自己)，那么显然这 $k$ 个要与其它的保持连通的，与中心必须有一条边，如果有多条边就形成环了，显然不满足生成树。另外 $n-k$ 个外围的其余点，自成一圈，所以有 $f[n-k]$ 种。我们可以枚举 $k$ ，则

$f [n] = \sum i = 0 i < n (i ? f [n ? i])$ $f[n]=\sum_{i=0}^{i<n}(i*f[n-i])$
初始化: $f[0]=f[1]=1,f[2]=3$
如果 $a,b$ 是连在一起的，如果与 $a,b$ 相连的为 $k$ 个(包括 $a,b$ )，那么 $a,b$ 是相邻的在这 $k$ 个位置选择就有 $k-1$ 种，而这 $k$ 个与中心相连的选择有 $k$ 种，剩下的与这部分是分开的，也是自成一圈，则为 $f[n-k]$ ，所以可以枚举 $k$ ，最终结果:

$g [n] = \sum i = 2 i < n (i ? (i ? 1) ? f [n ? i])$ $g[n]=\sum_{i=2}^{i<n}(i*(i-1)*f[n-i])$
初始化: $g[1]=0$

则最终的种数便是

T [n] = f [n] + g [n]

$T[n]=f[n]+g[n]$
前方高能
我们来化简:

f[n]和g[n] $f[n]和g[n]$

f [n] = \sum i = 0 i < n (i ? f (n ? i))

$f[n]=\sum_{i=0}^{i<n}(i*f(n-i))$

f [n] = f (n ? 1) + 2 ? f (n ? 2) + 3 ? f (n ? 3) \dots \dots (n ? 1) ? f (1) + n ? f [0] = f (n ? 1) + f (n ? 2) + f (n ? 3) \dots \dots f (1) + f (0) + (f (n ? 2) + 2 ? f (n ? 3) \dots \dots + (n ? 1) ? f (0))

$f[n]=f(n-1)+2*f(n-2)+3*f(n-3)……(n-1)*f(1)+n*f[0] \\ =f(n-1)+f(n-2)+f(n-3)……f(1)+f(0)+(f(n-2)+2*f(n-3)……+(n-1)*f(0))$
令

s[n] $s[n]$ 为

f[i] $f[i]$ 的前

n $n$ 项和，则上式可以写成:

$f[n] = s[n-1]+f(n-2)+2*f(n-3)……(n-2)f(1) \\= s[n-1]+\sum_{i=0}^{i<n-1}(i * f((n-1)-i)) \\=s[n-1]+f[n-1]\qquad (1)\quad s[n-1]=f[n]-f[n-1] \\=s[n-2]+f[n-1]+f[n-1] \qquad \qquad \\=s[n-2]+2*f[n-1]\qquad \qquad \\ =f[n-1]-f[n-2]+2*f[n-1] \qquad 根据(1)式对s[n-2]变形 \\ =3*f[n-1]-f[n-2] \qquad 其中f[0]=1,f[1]=1,f[2]=3,f[3]=8\\$

$g(n)=\sum_{i=1}^{i<n}[i(i-1)f(n-i)]$
$g(n)=1*2*f[n-2]+2*3*f[n-3]+3*4*f[n-4]\cdots(n-1)*(n-2)*f[1]$
$g(n-1)=\qquad 1*2*f[n-3]+2*3*f[n-4]……+(n-2)*(n-3)*f[1]$
则 $g(n)-g(n-1)=2*f[n-2]+4*f[n-3]……(2*(n-2))*f[1] \\=2*(f[n-2]+2*f[n-3]……+f[1])\\ =2*f[n-1]$

这个是最基本的递推式了。。
$g[n]=2*(f[1]+f[2]+f[3]……f[n-1])=2*(s[n-1]-f[0])\\ =2*(f[n]-f[n-1]-1) \qquad 其中f[0]=1$

(0 ? 1 13)

$\begin{pmatrix} 0 & 1 \\ -1 & 3 \end{pmatrix}$
右乘

(f (n ? 2) f (n ? 1))

$\begin{pmatrix} f(n-2) \\ f(n-1) \end{pmatrix}$

得到

(f (n ? 1) f (n))

$\begin{pmatrix} f(n-1)\\ f(n) \end{pmatrix}$
初始化：

(13)

$\begin{pmatrix} 1 \\ 3 \end{pmatrix}$
对于

f[n] $f[n]$ 的求法，可以用矩阵快速幂乘解决.
而

g[n] $g[n]$ 也就可以顺便得到，

T[n] $T[n]$ 就处理完毕了。
然后就是

Burnside定理 $Burnside定理$ ，同样

N $N$ 比较大，肯定是要用 欧拉函数优化，枚举循环个数。

n对MOD $n对MOD$ 是极有可能没有逆元的。
只能用

(a/b) $(a/b)%c=(a%(b*c))/b$ 。这样的话把取模就变成了

MOD?N $MOD*N$ ，范围一下子到了

1018 $10^{18}$ ，这样子的话中间的乘法便会溢出

64 $64$ 位整数。
所有的大整数相乘都得二分模拟。。。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <climits>
#include <cmath>
#include <ctime>
#include <cassert>
#include <bitset>
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
using namespace std;
typedef long long ll;
const int MAX_N = 100010;int prime[MAX_N], prime_cnt;
bitset<MAX_N> bs;void GetPrime()
{prime_cnt = 0;bs.set();for(int i = 2; i < MAX_N; ++i) {if(bs[i]) prime[prime_cnt++] = i;for(int j = 0; j < prime[i] && i * prime[j] < MAX_N; ++j) {bs[i * prime[j]] = 0;if(i % prime[j] == 0) break;}}
}struct Matrix {int row, col;ll data[10][10];
};inline ll mul_mod(ll a, ll b, ll mod)
{a = (a % mod + mod) % mod;b = (b % mod + mod) % mod;ll res = 0, tmp = a;while(b > 0) {if(b & 1) {res += tmp;if(res > mod) res -= mod;}tmp <<= 1;if(tmp > mod) tmp -= mod;b >>= 1;}return res;
}Matrix mat_mul_mod(Matrix a, Matrix b, ll mod) {Matrix res;res.row = a.row, res.col = b.col;for(int i = 1; i <= res.row; ++i) {for (int j = 1; j <= res.col; ++j) {res.data[i][j] = 0;for(int k = 1; k <= b.row; ++k) {res.data[i][j] += mul_mod(a.data[i][k], b.data[k][j], mod) % mod;}res.data[i][j] %= mod;}}return res;
}Matrix mat_quick_pow(Matrix a, ll m, ll mod) 
{Matrix res, tmp;res.row = tmp.row = a.row, res.col = tmp.col = a.col;memset(res.data, 0, sizeof(res.data));for(int i = 1; i <= res.row; ++i) { res.data[i][i] = 1; }memcpy(tmp.data, a.data, sizeof(tmp.data));while(m) {if(m & 1) res = mat_mul_mod(res, tmp, mod);tmp = mat_mul_mod(tmp, tmp, mod);m >>= 1;}return res;
}ll phi(ll x) 
{ll res = x;for(int i = 0; (ll)prime[i] * prime[i] <= x; ++i) {if(x % prime[i] == 0) {res = res / prime[i] * (prime[i] - 1);while(x % prime[i] == 0) x /= prime[i];}}if(x > 1) res = res / x * (x - 1);return res;
}ll solve(ll k, ll mod)
{if(k == 1) return 1;Matrix res, tmp;tmp.row = tmp.col = 2;tmp.data[1][1] = 0, tmp.data[1][2] = 1;tmp.data[2][1] = -1, tmp.data[2][2] = 3;tmp = mat_quick_pow(tmp, k - 2, mod);res.row = 2, res.col = 1;res.data[1][1] = 1, res.data[2][1] = 3;res = mat_mul_mod(tmp, res, mod);ll x = res.data[1][1], y = res.data[2][1];ll ans = mul_mod(y - x - 1, 2, mod);return ((ans + y) % mod) % mod;
}ll Polya(ll n, ll m) 
{ll res = 0, mod = n * m;for(ll i = 1; i * i <= n; ++i) {if(n % i) continue;res = (res + mul_mod(phi(n / i), solve(i, mod), mod)) % mod;if(n / i == i) break;res = (res + mul_mod(phi(i), solve(n / i, mod), mod)) % mod;}return res / n;
}int main()
{GetPrime();ll n, m;while(~scanf("%lld%lld", &n, &m)) {printf("%lld\n", Polya(n, m));}return 0;
}