被问到上次ZOJ Monthly中的H题Special Special Judge III怎么用容斥原理做。确实我只在解题报告中提了一句可以这么做,然后就丢下了一段python代码跑了。现在来详细解释一下如何运用容斥原理求一个Hyperrectangle被Hyperplane切割后剩余部分的体积。这里假定n维空间内的Hyperrectangle为

\{(x_1,x_2,\cdots,x_n)|0\le x_1\le X_1, 0\le x_2\le X_2, \cdots, 0\le x_n\le X_n\}

而对应的Hyperplane为
a_1x_1+a_2x_2+\cdots+a_nx_n=b

那么我们实际要求的就是Convex
\{\mathbf{x}|0\le x_i\le X_i,\mathbf{a}^T\mathbf{x}\le b\}

的体积。

不妨先考虑两个二维平面上的例子:

mp-case-1

\{(x,y)|0\le x\le 4,0\le y\le 3,x+y\le 5\}

此时的答案是
\triangle_{OXY}-\triangle_{AXB}-\triangle_{CDY}

mp-case-2

\{(x,y)|0\le x\le 4,0\le y\le 3,x+y\le 8\}

此时的答案是
\square_{OAEC}=\triangle_{OXY}-\triangle_{AXB}-\triangle_{CDY}+\triangle_{EDB}

之所以要加上\triangle_{EDB}是因为它在\triangle_{AXB}\triangle_{CDY}中被减掉了两次,多减了一次,需要加回来。如果考虑在第一个例子中\triangle_{EDB}=0的话,就有对任何情况,这个公式都成立。

这里,我们已经用上了容斥原理。现在尝试推广到n维的情况。首先,在高维空间里,这里就不再是三角形的面积了,而是单纯形(Simplex)的体积,单纯形

\{\mathbf{x}|0\le x_i,\mathbf{a}^T\mathbf{x}\le b\}

在b≤0时体积显然为0,而当b>0时,体积为\frac{1}{n!}\prod\frac{b}{a_i},不要忘了除以n!。那么我们要求的体积实际上就是很多个相似的(也就是b不同)的单纯形的和差。这些单纯形的顶点就是原先Hyperrectangle的顶点,一个n-Hyperrectangle有2n个顶点,可以用二进制数k来表示,如果第j位为0,则代表xj=0,否则xj=Xj,那么对应的单纯形是
\{\mathbf{x}|bit\_i\_of\_k \times X_i \le x_i, \mathbf{a}^T\mathbf{x}\le b\}

等价于
\{\mathbf{x}|0\le x_i,\mathbf{a}^T\mathbf{x}\le b-\sum_{bit\_i\_of\_k\ \textrm{is}\ 1}{a_iX_i}\}

知道了这2n个单纯形的体积后,剩下的问题就是我们应该加上还是减去这个单纯形的体积。这个问题由容斥原理(Inclusion–exclusion principle)来回答:如果二进制数k里有偶数个1,那么取正号,否则取负号。这样来看这个三维的例子(更高维的图我可画不出来……):

mp-case-3

\{(x,y,z)|0\le x\le 2,0\le y\le 4,0\le z\le 6,x+y+z\le 7\}

由前面推导知道,答案应该是
\begin{array}{cl}& V_{000}-V_{001}-V_{010}+V_{011}-V_{100}+V_{101}+V_{110}-V_{111}\\=& V_{000}-V_{001}-V_{010}-V_{100}+V_{011}\\=& (343-125-27-1+1)/6 = 191/6\\\end{array}

现在再贴上这段python代码应该很清楚了吧:

def gao(n, a, x, b0):
	ret = 0
	for i in range(2 ** n):
		b = b0
		c = 0
		for j in range(n):
			if i & (1 << j):
				b -= a[j] * x[j]
				c += 1
		if b <= 0:
			continue
		if c % 2:
			ret -= pow(b, n)
		else:
			ret += pow(b, n)
	for i in range(n):
		ret /= (i + 1) * a[i]
	return ret

if __name__ == '__main__':
	print gao(3, [1, 1, 1], [2, 4, 6], 7.0)

最后,SPOJ的2002. Random Number Generator (RNG)就是这样一道ai=1情况下的体积比问题。贴上我AC的cpp代码和超时的ruby代码

15 Responses to “用容斥原理计算被Hyperplane切割后的Hyperrectangle体积”
  1. shz says:

    为啥
    p /= x[i] * (i + 1);
    而不是
    p/=x[i]*2*(i+1)

  2. aswmtjdsj says:

    shi哥的latex插件是不是挂了? 看不到公式了呢。。。

  3. zgjabs says:

    SPOJ的2002. Random Number Generator (RNG) 题目给的随机数范围是 -xi 到 xi,请问是怎么转化成0 到 xi 来求解的
    也就是你代码中的这段,理解不了啊….
    p = 1;
    y = 0;
    for (int i = 0; i < n; ++i) {
    scanf("%Lf", &x[i]);
    p /= x[i] * (i + 1);
    y += x[i];
    }
    a = (a + y) / 2;
    b = (b + y) / 2;

    如果变成任意的 p到q ,又该怎么转化呢?

  4. Yuan says:

    请问为什么体积∏(b/ai)/n! ?

  5. yygy says:

    pow(b, n) 是什么意思

  6. MatRush says:

    Shi哥太给力了

  7.  
Leave a Reply