集合上的动态规划
例题：最优配对问题
空间里有n个点P0,P1,...,Pn-1，你的任务是把他们配成n/2对（n是偶数），使得每个点恰好在一个点对中。所有点对中两点的距离之和应尽量小。n<=20,|xi|,|yi|,|zi|<=1000.
样例输入：
20
1 2 3
1 1 1
5 6 2
4 7 8
2 3 1
1 4 7
2 5 8
3 6 9
1 2 5
2 3 6
4 5 2
7 8 5
4 5 1
-1 2 3
-1 -9 -7
0 0 0
100 0 0
9 5 1
7 5 3
5 5 5
分析：
既然每个点都要配对，很容易把问题看成如下多阶段决策过程：先确定P0和谁配对，然后是P1，接下来是P2，……，最后是Pn-1.按照前面的思路，设d(i)表示把前i个点两两配对的最小距离和，然后考虑第i个点的决策——它和谁配对呢？假设它和点j配对（j< i），那么接下来的问题应是“把前i-1个点中除了j之外的其它点两两配对”，它显然无法用任何一个d值来刻画——我们的状态定义无法体现出“除了一些点之外”这样的限制。
当发现状态无法转移后，常见的方法是增加维度，即增加新的因素，更细致地描述状态。既然刚才提到了“除了某些元素之外”，不妨把它作为状态的一部分，设d(i,S)表示把前i个点中，位于集合S中的元素两两配对的最小距离，则状态转移方程为：
d(i,S) = min{|PiPj|+ d(i-1,S-{i}-{j}) | j属于S}
其中|PiPj表示Pi和Pj之间的距离。边界是d(-1,S )= 0.
记忆化搜索的代码如下：
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cmath> 4 #define maxn 20 5 #define INF (1<<30) 6 int n; 7 struct Point{ 8 double x, y, z; 9 } p[maxn]; 10 double f[(1<<maxn)]; 11 inline double dist(int i, int j) 12 { 13 return sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x) + (p[i].y-p[j].y)*(p[i].y-p[j].y) + (p[i].z-p[j].z)*(p[i].z-p[j].z)); 14 } 15 inline double min(double x, double y) 16 { 17 if(x < y) return x; 18 return y; 19 } 20 double d(int S){ 21 //记忆化搜索快得要逆天了！！（因为奇数的都没有算！！ ） 22 if(f[S] != -1) return f[S]; 23 int i, j; 24 double& ans = f[S]; 25 ans = (S == 0? 0 : INF); 26 for(i = n-1; i >= 0; --i) if(S & (1<<i)) break; 27 for(j = 0; j < i; ++j) if(S & (1<<j)) 28 ans = min(ans, d(S^(1<<i)^(1<<j)) + dist(i, j)); 29 return ans; 30 } 31 32 int main(){ 33 int i; 34 scanf("%d", &n); 35 if(n%2) {printf("No answer!\n"); return 0;} 36 for(i = 0; i < n; ++i) scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z); 37 for(i = 0; i < (1 << n); ++i) f[i] = -1; 38 printf("%.3lf\n", d((1<<n)-1)); 39 /* 40 int cnt = 0; 41 for(i = 0; i < (1 << n); ++i) if(f[i]!=-1) ++cnt; 42 printf("%d\n", cnt); 43 */ 44 return 0; 45 }
动态规划的代码如下：
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cmath> 4 #define maxn 20 5 #define INF (1<<30) 6 int n; 7 struct Point{ 8 double x, y, z; 9 } p[maxn]; 10 double f[maxn][(1<<maxn)]; 11 inline double dist(int i, int j) 12 { 13 return sqrt((p[i].x-p[j].x)*(p[i].x-p[j].x) + (p[i].y-p[j].y)*(p[i].y-p[j].y) + (p[i].z-p[j].z)*(p[i].z-p[j].z)); 14 } 15 inline double min(double x, double y) 16 { 17 if(x < y) return x; 18 return y; 19 } 20 int main(){ 21 int i,j; 22 scanf("%d", &n); 23 if(n%2) {printf("No answer!\n"); return 0;} 24 for(i = 0; i < n; ++i) scanf("%lf%lf%lf", &p[i].x, &p[i].y, &p[i].z); 25 for(i = 0; i < n; ++i) 26 for(int S = 0; S < (1 << (i+1)); ++S){ 27 double& ans = f[i][S]; 28 ans = (S==0 ? 0 : INF); 29 if(S & (1 << i)) { 30 for(j = 0; j < i; ++j) if(S & (1 << j)) 31 ans = min(ans, f[i-1][S^(1<<i)^(1<<j)] + dist(i, j)); 32 } 33 else if(i!=0) f[i][S] = f[i-1][S]; 34 } 35 printf("%.3lf\n", f[n-1][(1<<n)-1]); 36 return 0; 37 }
隐含的阶段
上面的方程可以进一步简化。事实上，阶段 i 根本不用保存，它已经隐含在S中了——S中最大的元素就是i。这样，可直接用d(S)表示“把S中的元素两两配对的最小距离和”，则状态转移方程为：
d(S)=min{|PiPj|+d(S-{i}-{j}) | j属于S，i=max(S)}
状态有2n个，每个状态有O(n)种转移方式，总时间复杂度为O(n2n).
值得一提的是，不少人一直在用这样的状态转移方程：
d(S)=min{|PiPj|+d(S-{i}-{j}) | i,j属于S}
它和刚才的方程很类似，唯一不同是：i和j都是需要枚举的。这样做虽然也没错，但每个状态的转移次数高达O(n2)，总时间复杂度为O(n22n)，比刚才的方法慢。这个例子说明：即使用相同的状态描述，减少决策也很重要。
那么如何求出S中的最大元素呢？用一个循环判断即可。当S取遍{0,1,2,...,n-1}的所有子集时，平均判断次数仅为2。
1 for(int S=0; S<(1<<n); S++) 2 { 3 int i,j; 4 d[S] = INF; 5 for(i=0; i<n; i++) 6 if(S && (1<<i)) break; 7 for(j=i+1; j<n; j++) 8 if(S & (1<<j)) d[S] <?= d[S^(1<<i)^(1<<j)] + dist(i, j)); 9 }
注意，在上述的程序中求出的i是S中的最小元素，而不是最大元素，但这并不影响答案。另外，j的枚举只需从 i+1开始——既然 i 是S中最小元素，则说明其他元素自然均比i大。最后需要说明的是S的枚举顺序。不难发现：如果S‘是S的真子集，则一定有S'<S,因此S递增的顺序计算，需要用到某个d值时，它一定已经算出来了。
如果S'是S的真子集，则一定有S'<S。在用递推法实现子集动态规划时，该规则往往可以确定计算部分。
优化后代码如下：
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cmath> 4 #define maxn 22 5 #define maxs (1<<maxn) 6 #define INF 1e10 7 int n; 8 double f[maxs], x[maxn], y[maxn], z[maxn]; 9 inline double min(double x, double y) 10 { 11 return x<y?x:y; 12 } 13 inline double dist(int i, int j) 14 { 15 return sqrt((x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) + (z[i]-z[j])*(z[i]-z[j])); 16 } 17 int main(){ 18 scanf("%d", &n); 19 for(int i = 0; i < n; ++i) scanf("%lf%lf%lf", x+i, y+i, z+i); 20 int U = (1 << n); 21 for(int S = 0; S < U; ++S){ 22 int i; 23 for(i = n-1; i >= 0; --i) if(S & (1<<i)) break; 24 f[S] = (S==0) ? 0 : INF; 25 for(int j = 0; j < i ; ++j) if(S & (1<<j)){ 26 f[S] = min(f[S], f[S^(1<<i)^(1<<j)] + dist(i, j)); 27 } 28 } 29 printf("%.3lf\n", f[U-1]); 30 return 0; 31 }
