0212 组合立方体积
* * * *
拉格朗日计划
* * * *
组合立方体积

与坐标轴平行的立方体 $$\{(x,y,z): x_0\le x\le x_0+dx, y_0\le y\le y_0+dy, z_0\le z\le z_0+dz\},$$ 可以用参数$(x_0,y_0,z_0,dx,dy,dz)$唯一表示,它的体积是$dx\cdot dy\cdot dz$。

一系列立方体的组合体积是它们的并集的体积,若如果其中有立方体相交,则组合体积小于各自体积之和。

设有五万个与坐标轴平行的立方体$C_1,\ldots, C_{50000}$,$C_n$的参数为 \begin{align*} x_0&=S_{6n-5} \bmod 10000 \\ y_0&=S_{6n-4} \bmod 10000 \\ z_0&=S_{6n-3} \bmod 10000 \\ dx&=1+ (S_{6n-2} \bmod 399) \\ dy&=1+ (S_{6n-1} \bmod 399) \\ dz&=1+ (S_{6n} \bmod 399) \end{align*} 其中$S_k$是由延迟斐波那契生成器所生成的伪随机数,其公式如下: $$S_k=\begin{cases}100003-200003k+300007k^3 \pmod{1000000} & 1\le k\le 55, \\ S_{k-24}+S_{k-55} \pmod{1000000} & k\ge 56.\end{cases}$$ 因此,$C_1$的参数为$(7,53,183,94,369,56)$,$C_2$的参数为$(2383,3563,5079,42,212,344)$,依此类推。

前100个立方体的组合体积是723581599。

全部50000个立方体的组合体积是多少?

本题难度:



解答

本题的难点有两处:

一是在碰撞检测(即检验一对立方体是否有交集)的运算既需要合理剪枝,否则$50000^2$高达$25$亿难以计算,同时也要合理松弛,否则精确找出与给定立方体相交的其它立方体也十分繁琐。

二是同一区域可能是多个立方体的交集,但碰撞检测只能两两检验,需要正确计入或排除交集的体积。

注意到所有立方体的边长$dx,dy,dz$都不超过$399$,底角坐标$x_0,y_0,z_0$都不超过$10000$,因此可以将需要考察的区域$[0,10399]^3$划分为$400\times400\times400$的子区域,在每个子区域内作碰撞检测,如此可解难点一。

简单尝试就可以发现,直接遍历子区域内的立方体仍很难计算体积,因此针对难点二我们改为计算每个子区域与区域内立方体的交集体积。

计算时递归使用容斥原理,每次得到一个交集后就继续计算该交集和其它立方体的交集,并用一个标志根据递归次数的奇偶性改变符号以决定所得的体积应当计入还是排除。

最终结果是$328968937309$。

s=[(100003-200003*k+300007*k*k*k)%1000000 for k in range(1,56)]
for k in range(56,300001):
  s.append((s[-24]+s[-55])%1000000) 

cubes=[(s[6*k]%10000,s[6*k+1]%10000,s[6*k+2]%10000,s[6*k+3]%399+1,s[6*k+4]%399+1,s[6*k+5]%399+1) for k in range(50000)]

def getVolume(inspection,flag,curCube):
  s=0
  x1,y1,z1,dx1,dy1,dz1=curCube
  for i,j in enumerate(inspection):
      x2,y2,z2,dx2,dy2,dz2=cubes[j]
      u,v,w=max(x1,x2),max(y1,y2),max(z1,z2)
      du=min(x1+dx1, x2+dx2)-u
      dv=min(y1+dy1, y2+dy2)-v
      dw=min(z1+dz1, z2+dz2)-w
      if min(du,dv,dw)>0:
          s+=flag*du*dv*dw+getVolume(inspection[i+1:],-flag,(u,v,w,du,dv,dw))
  return s

blocks=[[[set() for k in range(26)] for j in range(26)] for i in range(26)]

for k in range(50000):
  x,y,z,dx,dy,dz=cubes[k]
  for a in [x,x+dx]:
      for b in [y,y+dy]:
          for c in [z,z+dz]:
              blocks[a/400][b/400][c/400].add(k)

s=0
for a in range(26):
  for b in range(26):
      for c in range(26):
          s+=getVolume(list(blocks[a][b][c]),1,(400*a,400*b,400*c,400,400,400))

print s