bzoj1005 [HNOI2008]明明的烦恼

1005: [HNOI2008]明明的烦恼

Time Limit: 1 Sec Memory Limit: 162 MB
Submit: 3032 Solved: 1209

Description

自从明明学了树的结构,就对奇怪的树产生了兴趣...... 给出标号为1到N的点,以及某些点最终的度数,允许在任意两点间连线,可产生多少棵度数满足要求的树?

Input

第一行为N(0 < N < = 1000),接下来N行,第i+1行给出第i个节点的度数Di,如果对度数不要求,则输入-1

Output

一个整数,表示不同的满足要求的树的个数,无解输出0

Sample Input

3
1
-1
-1

Sample Output

HINT

两棵树分别为1-2-3;1-3-2

Source

大意：给出n个点的度数限制，问这些点能组成多少棵树

分析:

首先介绍一种数列Purfer Code：

这种数列是这样子的：对于一棵树（这棵树是不存在根的），每次在数列中加入编号最小的孩子节点（即为度数为1的节点）所对应的父亲节点，并在树中删除此点，不断重复此操作，知道这棵树只剩下两个点

如：

/ \

2 3

/ \

4 5

所对应的PurferCode就是：1 2 2

因为第一次删掉3，第二次删掉1，第三次删掉4

剩下2、5两个节点

可以发现，每个PurferCode对应的树是唯一的

我们也可以从PurferCode得到一棵树

首先由PurferCode计算出每点的度数（出现的次数+1），

然后每次（第 i 次）选择度数最少（度数大于1）的点（度数相同时编号较小的优先）与数列第 i 个点连边，然后这两个点的度数都-1，

这样不断操作，会最后的到两个度数为1的点，将这两点相连，即的原来的树

设度数有限制的点的个数为Cnt，度数-1之和为Sum（即为数列中这些点要占的个数），第 i 个点的度数为di

那么可以根据组合数学：

从数列n-2个点中选择Sum个点有C(Sum, n-2)

这Sum个点中组合C(d1-1， Sum)*C(d2-1, Sum-(d1-1))*C(d3-1, Sum-(d1-1)-(d2-1))*........*C(dn - 1, dn - 1)

剩下的n-2-Sum个位置由其余n-Cnt个点组合，有(n-Cnt)^(n-2-Sum)种方案

所以Ans = C(Sum, n-2)*C(d1-1， Sum)*C(d2-1, Sum-(d1-1))*C(d3-1, Sum-(d1-1)-(d2-1))*........*C(dn - 1, dn - 1)*(n-Cnt)^(n-2-Sum)

因为C(Sum, n-2) = (n-2)!/(Sum!*(n-2-Sum)!),

C(d1-1， Sum) = Sum!/（(d1-1)!*(Sum-(d1-1))!）

C(d2-1, Sum-(d1-1)) = (Sum-(d1-1))!/（（d2-1)!*(Sum-(d1-1)-(d2-1))！）

.........

相乘时加粗部分可以相消

最终可得：

Ans = [ (n-2)!*((n-Cnt)^(n-2-Sum)) ] / [(n-2-Sum)!*(d1-1)!*(d2-1)!*.....*(dn-1)!]

写程序时，乘除都会很大，高精度除法和高精度乘高精度都很麻烦，所以可以分解质因数，最终可以只写单精度乘高精度即可。

关于分解x！的质因数

x = p^a * k（即x可以整除p^a）

a = x div p + x div p^2 + x div p^3 + ..... + x div p^m (p^m为p的阶乘中小于x的最大阶乘，即: p^m <= x, p^(m+1) > x)

那么 x！= p1^x1 * p2^x2 * p3^x3 * ...... * pn^xn

这个是显然的，大家写写就知道

关于分解质因数，可以用O(n)筛素数，这个不讲

综上所述，本题得解

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cstdlib>
  4 #include <cmath>
  5 #include <deque>
  6 #include <vector>
  7 #include <queue>
  8 #include <iostream>
  9 #include <algorithm>
 10 #include <map>
 11 #include <set>
 12 #include <ctime>
 13 using namespace std;
 14 typedef long long LL;
 15 #define For(i, s, t) for(int i = (s); i <= (t); i++)
 16 #define Ford(i, s, t) for(int i = (s); i >= (t); i--)
 17 #define MIT (2147483647)
 18 #define INF (1000000001)
 19 #define MLL (1000000000000000001LL)
 20 #define sz(x) ((bnt) (x).size())
 21 #define clr(x, y) memset(x, y, sizeof(x))
 22 #define puf push_front
 23 #define pub push_back
 24 #define pof pop_front
 25 #define pob pop_back
 26 #define ft first
 27 #define sd second
 28 #define mk make_pair
 29 inline void SetIO(string Name) {
 30     string Input = Name+".in",
 31     Output = Name+".out";
 32     freopen(Input.c_str(), "r", stdin),
 33     freopen(Output.c_str(), "w", stdout);
 34 }
 35 
 36 const int N = 1010, LEN = 1010, M = 10000;
 37 int Prime[N], p;
 38 bool Visit[N];
 39 int n, D[N], Cnt, Sum;
 40 int P1[N], P2[N];
 41 int Ans[LEN], Len;
 42 
 43 inline void Input() {
 44     scanf("%d", &n);
 45     For(i, 1, n) scanf("%d", &D[i]);
 46 }
 47 
 48 inline void EXIT(int Ans) {
 49     printf("%d\n", Ans);
 50     exit(0);
 51 }
 52 
 53 inline void BuildPrime() {
 54     clr(Visit, 0), p = 0;
 55     For(i, 2, n){
 56         if(!Visit[i]) Prime[++p] = i;
 57         For(j, 1, p) {
 58             if(i*Prime[j] > n) break;
 59             Visit[i*Prime[j]] = 1;
 60             if(i%Prime[j] == 0) break;
 61         }
 62     }
 63 }
 64 
 65 inline void FenJie(int A, int *P) {
 66     For(i, 1, p) {
 67         int x = Prime[i];
 68         while(x <= A) {
 69             P[i] += A/x;
 70             x *= Prime[i];
 71         }
 72     }
 73 }
 74 
 75 inline void Mul(int x) {
 76     For(i, 1, Len) Ans[i] *= x;
 77     For(i, 1, Len) {
 78         Ans[i+1] += Ans[i]/M;
 79         Ans[i] %= M;
 80     }
 81     while(Ans[Len+1]) {
 82         Len++;
 83         Ans[Len+1] += Ans[Len]/M;
 84         Ans[Len] %= M;
 85     }
 86 }
 87 
 88 inline void Solve() {
 89     //Check
 90     if(n == 1) {
 91         if(D[1] == 0 || D[1] == -1) EXIT(1);
 92         EXIT(0);
 93     }
 94     For(i, 1, n)
 95         if(D[i] == 0) EXIT(0);
 96     Cnt = Sum = 0;
 97     For(i, 1, n)
 98         if(D[i] > 0) Cnt++, Sum += D[i]-1;
 99     if(Sum > n-2 || (Sum < n-2 && n == Cnt)) EXIT(0);
100     
101     //Work
102     BuildPrime();
103     clr(P1, 0), clr(P2, 0);
104     FenJie(n-2, P1);
105     FenJie(n-2-Sum, P2);
106     For(i, 1, n)
107         if(D[i] > 1) FenJie(D[i]-1, P2);
108     int m = n-Cnt;
109     For(i, 1, p) {
110         if(m <= 1) break;
111         int t = 0, k = Prime[i];
112         while(!(m%k)) m /= k, t++;
113         P1[i] += t*(n-2-Sum);
114     }
115     For(i, 1, p) P1[i] -= P2[i];
116     
117     Ans[1] = 1, Len = 1;
118     For(i, 1, p)
119         For(j, 1, P1[i]) Mul(Prime[i]);
120     
121     printf("%d", Ans[Len]);
122     Ford(i, Len-1, 1) printf("%04d", Ans[i]);
123     cout<<endl;
124 }
125 
126 int main() {
127     SetIO("1005");
128     Input();
129     Solve();
130     return 0;
131 }

View Code

下面聊聊用线性筛素数法求欧拉函数

欧拉函数有两个定理

1、当m，n互质时，phi(m*n) = phi(m)*phi(n);

2、phi(p^k) = (p-1)*p^(k-1)

因为1-p^k间只有p的倍数不与其互质，这样的数有p^k/p个，即p^(k-1)个，即phi(p^k) = p^k - p^(k-1) = (p-1)*p^(k-1)

那么当n%p == 0 （显然不互质，p为质数）时，

设n = m * p^k（m与p互质）

phi(n) = phi(m)*phi(p^k) = phi(m)*(p-1)*p^(k-1)

phi(n*p) = phi(m*p^(k+1)) = phi(m)*(p-1)*p^k = phi(n)*p

所以可以得求欧拉函数代码

 1 for(int i=2;i<N;i++)  {
 2     if(!vis[i]) {
 3         prime[++cnt]=i;
 4         phi[i]=i-1;
 5     }
 6     for(int k=1;k<=cnt&&prime[k]*i<N;k++) {
 7         x=prime[k]*i;
 8         vis[x]=true;
 9         if(i%prime[k]==0) {
10             phi[x]=phi[i]*prime[k];
11             break;
12         }
13         else phi[x]=phi[i]*(prime[k]-1);
14     }
15 }