Problem 252 「凸ホール」 †
平面上に与えられた点の集合に対し, 以下を満たす凸多角形を凸ホール(convex hole)と定義する:
頂点は与えられた点のいくつかから成り, 内部に与えられた点を含まない(頂点以外に, 多角形の辺上に与えられた点があっても構わない)
頂点は与えられた点のいくつかから成り, 内部に与えられた点を含まない(頂点以外に, 多角形の辺上に与えられた点があっても構わない)
例として, 下の図は 20 個の点とそれに対するいくつかの凸ホールを示している. 赤い多角形で示した凸ホールは 1049694.5 の単位正方形と面積が等しく, この点の集合に対し最大の凸ホールである.
この例では, 次の擬似乱数によって生成された 最初の 20 個の点 (T2k−1, T2k) (k = 1,2,…,20) を使用した.
S0 = 290797Sn+1 = Sn2 mod 50515093Tn = (Sn mod 2000 ) − 1000
すなわち, (527, 144), (−488, 732), (−454, −947), ... である.
この擬似乱数生成器による最初の 500 個の点を使用する凸ホールの中で, 最大の面積を求めよ.
小数点以下に1桁をつけて回答を入力せよ.
小数点以下に1桁をつけて回答を入力せよ.
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20252
解法
一分ルールは守れていません。
合格はしましたが手元の最新のパソコンで2分20秒も計算時間がかかっています。
点に0～499まで番号を付けており
コード実行するとその番号の点から始まる”とつ多角形”の検証が全部すんだら対応した数字が0から499まで表示されます。
最後に出てくる数字が答えとなる面積です。
スタートの点と2番目の点を決めて探索を始めます。
反時計回りで半開平面で片方にある点だけを選びながら探索を行います。
できるだけ大きな多角形になるように、折れ線の角度変化が小さいものから優先して探索し、
今いる点とその前の点の2セットが前に探索した状態よりも同じか小さな面積になってるなら探索を打ち切って別の探索へ向かいます。
それと始点を輪を作る中で一番番号の小さい点だと仮定して探索します。
ただし反時計回りでしか探索しないので二つ目の点のみ一つ目の点より小さな番号の点であるパターンを許容します。
これで少しだけ早くなります。
最新のパソコンで112.75秒。
まだまだ遅いです。
小手先のテクニックをいろいろきかしてみましたが対して効果が表れていません。
根本的な発想の転換が必要なようです。
#include<stdio.h>
#include<queue>
#include<vector>
#include<iostream>
#include<string.h>
#include<set>
#include<time.h>
const __int64 MOD_S=50515093;
const int POINT_SIZE=500;
std::vector<int> xs,ys;
std::set<int> alivePoints;
int dp[POINT_SIZE][POINT_SIZE];
int revDP[POINT_SIZE][POINT_SIZE];
struct S{
int p1,p2,addArea;
double r;
bool operator<(const S& s)const{
return r<s.r;
}
};
int calc_area(int p1,int p2,int p3){
int dx1,dy1,dx2,dy2;
dx1=xs[p1]-xs[p3];
dx2=xs[p2]-xs[p3];
dy1=ys[p1]-ys[p3];
dy2=ys[p2]-ys[p3];
return dx1*dy2-dx2*dy1;
}
bool is_in_area(int start,int p2,int p3){
std::set<int>::iterator it;
for(it=alivePoints.begin();it!=alivePoints.end();it++){
int i=(*it);
if(i==start||i==p2||i==p3)continue;
int r1=calc_area(start,p2,i);
int r2=calc_area(p2,p3,i);
int r3=calc_area(p3,start,i);
if((r1>=0&&r2>=0&&r3>=0)||(r1<=0&&r2<=0&&r3<=0)){
return true;
}
}
return false;
}
int search(int start,int oldPoint,int nowPoint,int area){
int resultArea=0;
if(dp[nowPoint][start]<area){
dp[nowPoint][start]=area;
}
std::priority_queue<S> pq;
S s1;
s1.p1=nowPoint;
std::set<int>::iterator it;
std::set<int> dells;
int reAreaAdd=calc_area(start,oldPoint,nowPoint);
for(it=alivePoints.begin();it!=alivePoints.end();it++){
int i=(*it);
if(i<start)continue;
int aaa=calc_area(nowPoint,i,oldPoint);
s1.addArea=calc_area(start,nowPoint,i);
if((aaa<0)||s1.addArea<0||is_in_area(start,nowPoint,i)){
dells.insert(i);
continue;
}
s1.p2=i;
double len1,len2;
len1=hypot(xs[nowPoint]-xs[start],ys[nowPoint]-ys[start]);
len2=hypot(xs[i]-xs[start],ys[i]-ys[start]);
if(len1==0||len2==0)continue;
int naiseki=(xs[i]-xs[nowPoint])*(xs[nowPoint]-xs[oldPoint]);
naiseki+=(ys[i]-ys[nowPoint])*(ys[nowPoint]-ys[oldPoint]);
s1.r=(1.0*naiseki)/(len1*len2);
pq.push(s1);
}
for(it=dells.begin();it!=dells.end();it++){
alivePoints.erase((*it));
}
while(pq.empty()==false){
s1=pq.top();
pq.pop();
int area1=area+s1.addArea;
if(dp[s1.p1][s1.p2]>=area1)continue;
if(dp[s1.p1][s1.p2]==-1){
dp[s1.p1][s1.p2]=area1;
alivePoints.erase(s1.p2);
int t=search(start, s1.p1, s1.p2,area1);
if(resultArea<t){
resultArea=t;
}
alivePoints.insert(s1.p2);
}else{
dp[s1.p1][s1.p2]=area1;
int t=revDP[s1.p1][s1.p2];
if(resultArea<t){
resultArea=t;
}
}
}
alivePoints.insert(dells.begin(),dells.end());
revDP[oldPoint][nowPoint]=resultArea+reAreaAdd;
return resultArea+reAreaAdd;
}
void set_alive_points(int start,int second){
alivePoints.clear();
for(int i=0;i<POINT_SIZE;i++){
if((i!=start) && (i!=second) && (calc_area(start,second,i)>0)){
alivePoints.insert(i);
}
}
}
void search_w(){
int ans=0;
for(int i=0;i<POINT_SIZE;i++){
printf("%d ",i);
for(int j=0;j<POINT_SIZE;j++){
if(i==j) continue;
memset(dp,-1,sizeof(dp));
memset(revDP,-1,sizeof(revDP));
dp[i][j]=0;
set_alive_points(i,j);
int t=search(i,i,j,0);
if(ans<t)ans=t;
}
}
printf("%lf\n",ans/2.0);
}
int main(){
__int64 Si=290797;
clock_t start,end;
start = clock();
Si=(Si*Si)% MOD_S;
for(int i=0;i<POINT_SIZE;i++){
xs.push_back(Si%2000-1000);
Si=(Si*Si)% MOD_S;
ys.push_back(Si%2000-1000);
Si=(Si*Si)%MOD_S;
}
search_w();
end = clock();
printf("%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC);
}
0 件のコメント:
コメントを投稿