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是断开的种数。$$
$$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[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[1]=0$

则最终的种数便是

$$T[n]=f[n]+g[n]$$
前方高能
我们来化简: $f[n]和g[n]$
$$f[n]=\sum_{i=0}^{i<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]+\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$

$$\begin{pmatrix} 0 & 1 \\ -1 & 3 \end{pmatrix}$$
右乘
$$\begin{pmatrix} f(n-2) \\ f(n-1) \end{pmatrix}$$

得到

$$\begin{pmatrix} f(n-1)\\ f(n) \end{pmatrix}$$
初始化：
$$\begin{pmatrix} 1 \\ 3 \end{pmatrix}$$
对于 $f[n]$的求法，可以用矩阵快速幂乘解决.
而 $g[n]$也就可以顺便得到， $T[n]$就处理完毕了。
然后就是 $Burnside定理$，同样 $N$比较大，肯定是要用欧拉函数优化，枚举循环个数。
$n对MOD$是极有可能没有逆元的。
只能用 $(a/b)%c=(a%(b*c))/b$。这样的话把取模就变成了 $MOD*N$，范围一下子到了 $10^{18}$，这样子的话中间的乘法便会溢出 $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;
}