HDU 2481 Toy(Polya综合)

题目链接：

HDU 2481 Toy

题意：

有外围n个点围着中心一个点，中心有n条边和外围相连，外围的相邻点有一条边，现在需要从中去除n条边，使得剩下的n+1个点依然保持两两连通，考虑旋转，问有多少种方案？方案数对M取模。

数据范围：3≤n≤109,2≤M≤109.

分析：

我觉得这种考验智商的题目真的不适合窝。。。

主要参考博客：

cx_love

将狼踩尽

zxb

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

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

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

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

f[n]=∑i=0i<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]=∑i=2i<n(i∗(i−1)∗f[n−i])

初始化:g[1]=0

则最终的种数便是

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

前方高能

我们来化简:f[n]和g[n]

f[n]=∑i=0i<n(i∗f(n−i))

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]为f[i]的前n项和，则上式可以写成:

f[n]=s[n−1]+f(n−2)+2∗f(n−3)……(n−2)f(1)=s[n−1]+∑i<n−1i=0(i∗f((n−1)−i))=s[n−1]+f[n−1](1)s[n−1]=f[n]−f[n−1]=s[n−2]+f[n−1]+f[n−1]=s[n−2]+2∗f[n−1]=f[n−1]−f[n−2]+2∗f[n−1]根据(1)式对s[n−2]变形=3∗f[n−1]−f[n−2]其中f[0]=1,f[1]=1,f[2]=3,f[3]=8

g(n)=∑i<ni=1[i(i−1)f(n−i)]

g(n)=1∗2∗f[n−2]+2∗3∗f[n−3]+3∗4∗f[n−4]⋯(n−1)∗(n−2)∗f[1]

g(n−1)=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)其中f[0]=1

(0−113)

右乘

(f(n−2)f(n−1))

得到

(f(n−1)f(n))

初始化：

(13)

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

而g[n]也就可以顺便得到，T[n]就处理完毕了。

然后就是Burnside定理，同样N比较大，肯定是要用欧拉函数优化，枚举循环个数。

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

只能用(a/b)。这样的话把取模就变成了MOD∗N，范围一下子到了1018，这样子的话中间的乘法便会溢出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;
}