Problem 212 「結合直方体の体積」 †
座標軸に平行な直方体 (axis-aligned cuboid) は {(x0, y0, z0), (dx, dy, dz)} で与えられ, x0 ≤ X ≤ x0 + dx, y0 ≤ Y ≤ y0 + dy, z0 ≤ Z ≤ z0 + dz, を満たす点で構成される. 直方体の体積は dx × dy × dzで求められる. 複数の直方体を結合したものの体積を考えた場合, 直方体に重なりがあれば, 結合直方体の体積は それぞれの直方体の体積の和より小さくなる.
C1, …, C50000 を以下のパラメータで与えられる座標軸に平行な直方体とする.
- x0 = S6n-5 modulo 10000
- y0 = S6n-4 modulo 10000
- z0 = S6n-3 modulo 10000
- dx = 1 + (S6n-2 modulo 399)
- dy = 1 + (S6n-1 modulo 399)
- dz = 1 + (S6n modulo 399)
S1,…,S300000 はラグ付きフィボナッチ法により生成される.
- 1 ≤ k ≤ 55 の場合, Sk = [100003 - 200003k + 300007k3] (modulo 1000000)
- 56 ≤ k の場合, [Sk-24 + Sk-55] (modulo 1000000)
したがって, C1 は {(7, 53, 183), (94, 369, 56)}, C2 は {(2383, 3563, 5079), (42, 212, 344)} となる
C1, …, C100 の結合直方体の体積は 723581599 である.
C1, …, C50000 の結合直方体の体積を求めよ.
解法
ＸＹ平面と平行な平面で直方体を切断し、断面の形の変化がない範囲を個別に計算しました。
最新のパソコンで36秒。
ちょっと遅いですがまあ解けただけでも良しという気分です。
正答者掲示板で他の方の回答を見てみました。
目から鱗な回答がいっぱいありますね。
jaapさんの解法がコンパクトでエレガント。
znezicさんの解法は魔法使いなみ。
どちらのコードもきちんと理解してみたい。
私のは我流で思いついた方法なのでもしかしたらほかの方の解法は一般的な解法なのかもしれません。
私はきっと井の中の蛙。
#include<stdio.h>
#include<set>
#include<map>
#include<iostream>
#include<time.h>
const int DELL=-1;
const int NEW=1;
const int ADD=1;
struct S{
__int64 x,y;
int zType,xyType,no;
bool operator<(const S& s)const{
if(x!=s.x)return x<s.x;
if(y!=s.y)return y<s.y;
if(xyType!=s.xyType)return xyType<s.xyType;
return no<s.no;
}
};
__int64 xyArea(std::set<S>& nowXYPoints,std::set<S>& points){
std::set<S>::iterator it;
S s1;
for(it=nowXYPoints.begin();it!=nowXYPoints.end();it++){
s1=(*it);
if(s1.zType==DELL){
points.erase(s1);
}else{
points.insert(s1);
}
}
it=points.begin();
if(it==points.end())return 0;
int x=(*it).x;
__int64 ysum,resultArea=0;
int count;
std::map<int,int> aliveY;
std::map<int,int>::iterator yIt;
std::set<int> dells;
std::set<int>::iterator dellIt;
while(it!=points.end()){
__int64 dx;
//std::cout<<"\nb";
while(it!=points.end()){
s1=(*it);
if(x<s1.x)break;
if(aliveY.find(s1.y)==aliveY.end())aliveY[s1.y]=0;
aliveY[s1.y]+=s1.xyType;
it++;
}
dx=s1.x-x;
x=s1.x;
dells.clear();
for(yIt=aliveY.begin();yIt!=aliveY.end();yIt++){
//std::cout<<(*yIt).first<<" "<<(*yIt).second<<")";
if((*yIt).second==0)dells.insert((*yIt).first);
}
for(dellIt=dells.begin();dellIt!=dells.end();dellIt++){
aliveY.erase((*dellIt));
}
yIt=aliveY.begin();
if(yIt==aliveY.end())continue;
int count=(*yIt).second;
int oldY=(*yIt).first;
__int64 sumY=0;
for(yIt++;yIt!=aliveY.end();yIt++){
count+=(*yIt).second;
if(count==0){
sumY+=(*yIt).first-oldY;
yIt++;
if(yIt==aliveY.end()) break;
oldY=(*yIt).first;
count=(*yIt).second;
}
}
resultArea+=(dx*sumY);
}
return resultArea;
}
const int LIMIT=6*50000;
const __int64 MOD=1000000;
const int MODXYZ=10000;
const int MOD_D=399;
std::map<int,std::set<S> > allPoints;
__int64 sk[LIMIT+1];
void append_point(int x,int y,int z,int xyType,int zType,int no){
std::set<S> dammy;
S s1;
s1.xyType=xyType;
s1.zType=zType;
s1.x=x;
s1.y=y;
s1.no=no;
if(allPoints.find(z)==allPoints.end())allPoints[z]=dammy;
allPoints[z].insert(s1);
}
int main(){
clock_t start,end;
start = clock();
for(__int64 k=1;k<=LIMIT;k++){
if(k<=55){
sk[(int)k]=(100003 - 200003*k + 300007*k*k*k)%MOD;
}else{
sk[(int)k]=(sk[(int)k-24]+sk[(int)k-55])%MOD;
}
}
S s1;
int x0,y0,z0,dx,dy,dz;
for(int i=1;i<LIMIT;i+=6){
x0=sk[i]%MODXYZ;
y0=sk[i+1]%MODXYZ;
z0=sk[i+2]%MODXYZ;
dx=1+sk[i+3]%MOD_D;
dy=1+sk[i+4]%MOD_D;
dz=1+sk[i+5]%MOD_D;
append_point(x0 ,y0 ,z0,ADD ,NEW,i);
append_point(x0 ,y0+dy,z0,DELL,NEW,i);
append_point(x0+dx,y0 ,z0,DELL,NEW,i);
append_point(x0+dx,y0+dy,z0,ADD ,NEW,i);
append_point(x0 ,y0 ,z0+dz,ADD ,DELL,i);
append_point(x0 ,y0+dy,z0+dz,DELL,DELL,i);
append_point(x0+dx,y0 ,z0+dz,DELL,DELL,i);
append_point(x0+dx,y0+dy,z0+dz,ADD ,DELL,i);
}
__int64 ans=0,area,oldZ,t;
std::map<int,std::set<S> >::iterator it;
std::set<S> alivePoints;
std::cout<<"a";
it=allPoints.begin();
area=xyArea((*it).second,alivePoints);
oldZ=(*it).first;
for(it++;it!=allPoints.end();it++){
ans+=((*it).first-oldZ)*area;
oldZ=(*it).first;
area=xyArea((*it).second,alivePoints);
}
std::cout<<"ans="<<ans;
end = clock();
printf("\n%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC);
}
