プロジェクトオイラー問353

https://projecteuler.net/problem=353
完璧な球体の月面表面の整数座標を通って北極点から南極点へ到達するとき最も安全なルートをとった時の危険度をこたえる問題。

動的計画法で解。
こんなトイプロブレム(玩具の問題)が解けても技術屋になれるわけでもなしどうして俺はこんな問題を解いてるんだろう。
計算時間10分駄目駄目ですね。
最後に出力されるsumAnsが答え。
答えは1桁多目に出てますので目視で修正してください。


#include<stdio.h>
#include<math.h>
#include<map>
#include<set>
#include<algorithm>
#include<iostream>
#include<vector>
#include<queue>


struct E{
	__int64 x,y,z;
	void set(int x1,int y1,int z1){
		x=x1;
		y=y1;
		z=z1;
	}
	bool operator<(const E& e1)const{
		if(z!=e1.z)return z>e1.z;
		if(x!=e1.x)return x<e1.x;
		return y<e1.y;
	}
};


double calcR(const E& e1,const E& e2,double r){
	double cosA=(e1.x*e2.x+e1.y*e2.y+e1.z*e2.z)/(r*r);
	if(1<=cosA)return 0;
	if(cosA<=-1)return M_PI;
	return acos(cosA);
} 

double calcRisk(const E& e1,const E& e2,double r){
	double r2=calcR(e1,e2,r);
	double r3=r2/M_PI;
	return r3*r3;
} 
double calcRisk4(const E& e1,const E& e2,double r){
	double re;
	E e3=e2;
	re=calcRisk(e1,e3,r);
	e3.x=-e2.y;
	e3.y=e2.x;
	re=std::min(re,calcRisk(e1,e3,r));
	e3.x=-e2.x;
	e3.y=-e2.y;
	re=std::min(re,calcRisk(e1,e3,r));
	e3.x=e2.y;
	e3.y=-e2.x;
	re=std::min(re,calcRisk(e1,e3,r));
	return re;
} 
double calcRisk8(const E& e1,const E& e2,double r){
	E e3=e2;
	e3.x=e2.y;
	e3.y=e2.x;
	double a=calcRisk4(e1,e2,r);
	double b=calcRisk4(e1,e3,r);
	return std::min(a,b);
} 

double calc(int size){
	std::map<__int64,__int64> pow2;
	std::map<E,double> ms;
	__int64 r=pow(2,size)-1;
	for(__int64 i=0;i<=r;i++){
		pow2[i*i]=i;
	}
	
	__int64 z1=0;
	for(__int64 z=0;z<=r;z++){
	__int64 l2=r*r-z*z;
		for(__int64 x=0;x*x<=l2;x++){
			__int64 l3=l2-x*x;
			if(pow2.find(l3)!=pow2.end()){
				__int64 y=pow2[l3];
				if(y<x)break;
				E e1;
				e1.set(x,y,z);
 				ms[e1]=10000;
				if(z>0){
					if(z1==0||z1==z){
						z1=z;
						e1.set(x,y,-z);
						ms[e1]=10000;
					}
				}
			}
		}
	}
	
	
	std::map<E,double>::iterator it,it2;
	E e1;
	e1.set(0,0,r);
	ms[e1]=0;
	std::set<E> isCommit;
	
	for(int i=0;i<ms.size();i++){
		if(i%600==0)std::cout<<"("<<i<<"/"<<ms.size()<<")";
		double nowCost=-1;
		E te;
 	
		for(it=ms.begin();it!=ms.end();it++){
			if(isCommit.find((*it).first)==isCommit.end()){
				if((nowCost==-1)||(nowCost>(*it).second)){
					nowCost=(*it).second;
					te=(*it).first;
				}
			}
		}
		isCommit.insert(te);
		e1=te;
		for(it=ms.begin();it!=ms.end();it++){
 			double cost=(*it).second;
			if(cost<=nowCost)continue;
			double cost2=calcRisk8(e1,(*it).first,r)+nowCost;
			if(cost>cost2){
				(*it).second=cost2;
			}
		}
	}
	
	
	
 	
	
	double ans=-1;
	__int64 cc=0;
	for(it=ms.begin();it!=ms.end();it++){
		double cost1=(*it).second;
		E e1=(*it).first;
		cc++;
 		if(cc%600==0)std::cout<<"("<<cc<<" "<<ms.size()<<")";
		for(it2=it;it2!=ms.end();it2++){
			double cost2=(*it2).second+cost1;
			E e2=(*it2).first;
			e2.z=-e2.z;
			double cost3=cost2+calcRisk8(e1,e2,r);
 			if((ans==-1)||(ans>cost3)){
				ans=cost3;	
			}
		}
	}
 	printf("\n%d ans=%.11lf\n",size,ans);
	return ans;
}

int main(){
	double lastAns=0;
 	for(int i=1;i<=15;i++){
 		lastAns+=calc(i);
	}
	printf("sumAns=%.11lf\n",lastAns);
 }
最終更新:2015年11月15日 05:30