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
2
HINT
两棵树分别为1-2-3;1-3-2
Source
大意：给出n个点的度数限制，问这些点能组成多少棵树
分析:
首先介绍一种数列Purfer Code：
这种数列是这样子的：对于一棵树（这棵树是不存在根的），每次在数列中加入编号最小的孩子节点（即为度数为1的节点）所对应的父亲节点，并在树中删除此点，不断重复此操作，知道这棵树只剩下两个点
如：
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 }
下面聊聊用线性筛素数法求欧拉函数
欧拉函数有两个定理
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 }