1081 - The Return of Carl

Carl the ant is back! When we last left him (Problem A, 2004 World Finals), Carl was a little mixed-up, always taking strange, zigzag paths when traveling. But now, Carl has straightened out his life - literally. He now always takes the straightest, shortest path between any pair of points. This sounds simple, except for one small thing: Carl now spends most of his time on a paperweight in the shape of a regular octahedron. You may recall that an octahedron is one of the five Platonic solids and consists of eight equilateral triangles, as shown in Figure 3.

Figure 3: Regular octahedron

Carl has an innate (some say in-ant) ability to always take the shortest path when going from any starting point to any destination point on the paperweight. Your job is to verify this by determining the length of the shortest path when given two such points, not necessarily distinct.

Input 

The input contains multiple test cases. Each test case consists of four integers , ,  and , where 0 < 360 and 0180. The first two are the spherical coordinates of the start point, and the last two are the spherical coordinates of the destination point. As shown in Figure 3,  is the azimuthand  is the zenith angle, both of them given in degrees.

The input is terminated by a line containing four negative ones.

The paperweight is fixed for all test cases as follows: the octahedron is centered on the origin and each vertex lies on one of the axes. Every edge is exactly 10 cm long. You should suppose that Carl's size is zero and ignore any supporting mechanisms that may be necessary to hold the paperweight in the correct position.

Output 

For each test case, print the case number (starting with 1) and the length in centimeters of the shortest path from the start point to the destination point, rounded to the nearest thousandth. Follow the format of the sample output.

Sample Input 

0 90 90 90
0 90 90 45
0 0 0 180
-1 -1 -1 -1

Sample Output 

Case 1: 10.000
Case 2: 8.660
Case 3: 17.321

#include<iostream>#include<stdio.h>#include<string.h>#include<stdlib.h>#include<math.h>using namespace std;#define FOR(I,A,B) for(I=int(A);I<int(B);++I)#define MEM(A,B) memset(A,B,sizeof(A));#define CPY(A,B) memcpy(A,B,sizeof(B));const double pi=acos(-1.0);struct TPoint{	double x,y,z;	void get(double a,double b,double c)	{		x=a;y=b;z=c;	}};TPoint p[5];TPoint A,B;int a1,a2,a3,a4,ca;double l;void init(){	l=5*sqrt(2.);	p[0].get(0,0,l);	p[1].get(l,0,0);	p[2].get(0,l,0);	p[3].get(-l,0,0);	p[4].get(0,-l,0);}double eps(double x){	return x<1e-6&&x>-1e-6;}void getpoint(double a1,double a2,TPoint &A){	double l2,r,aa=(int)a1%90;	a1=a1/180*pi;	a2=a2/180*pi;	aa=aa/180*pi;	l2=l/(sin(aa)+cos(aa));	r=l2*l/(l*sin(a2) + l2*cos(a2));	A.x=r*sin(a2)*cos(a1);	A.y=r*sin(a2)*sin(a1);	A.z=r*cos(a2);}double sqr(double x){	return x*x;}double dis(TPoint a,TPoint b){	return sqrt(sqr(a.x-b.x) + sqr(a.y-b.y) + sqr(a.z-b.z));}int pre(int x){	return (x+3)%4;}int nxt(int x){	return (x+1)%4;}TPoint point(double ra,int x,int y){	TPoint ret;	ret.x=(p[x].x-p[y].x)*ra + p[y].x;	ret.y=(p[x].y-p[y].y)*ra + p[y].y;	ret.z=(p[x].z-p[y].z)*ra + p[y].z;	return ret;}double getang(double a,double b,double c){	if(eps(c)||eps(a-b-c))		return 0;	return acos((sqr(a)+sqr(b)-sqr(c))/(2*a*b));}double cal(int p1,int p2,TPoint A,TPoint B){	TPoint C,D;	double L,R,mid1,mid2,ret(1e10);	if(a4>90)	{		int i;		a4=180-a4;		FOR(i,0,4)		{			L=0;R=1;			while(R-L>1e-4)			{				mid1=L+(R-L)/3;				C=point(mid1,i+1,nxt(i)+1);				mid2=R-(R-L)/3;				D=point(mid2,i+1,nxt(i)+1);				if(cal(p1,i,A,C)+cal(i,p2,C,B)>cal(p1,i,A,D)+cal(i,p2,D,B))					L=mid1;				else					R=mid2;			}			ret=min(ret,cal(p1,i,A,C)+cal(i,p2,C,B));		}		return ret;	}	if(p1==p2)		return dis(A,B);	if(nxt(p2)==p1)	{		swap(p1,p2);		swap(A,B);	}	if(nxt(p1)==p2)	{		double an,a=dis(A,p[0]),b=dis(B,p[0]);		an=getang(10,a,dis(A,p[nxt(p1)+1]))+getang(10,b,dis(B,p[nxt(p1)+1]));		return sqrt(sqr(a)+sqr(b)-2*a*b*cos(an));	}	int i;	L=0;R=1;	while(R-L>1e-4)	{		mid1=L+(R-L)/3;		C=point(mid1,0,p1+1);		mid2=R-(R-L)/3;		D=point(mid2,0,p1+1);		if(dis(A,C)+cal(pre(p1),p2,C,B)>cal(pre(p1),p2,D,B)+dis(A,D))			L=mid1;		else			R=mid2;	}	ret=cal(pre(p1),p2,C,B)+dis(A,C);	L=0;R=1;	while(R-L>1e-4)	{		mid1=L+(R-L)/3;		C=point(mid1,0,nxt(p1)+1);		mid2=R-(R-L)/3;		D=point(mid2,0,nxt(p1)+1);		if(cal(nxt(p1),p2,C,B)+dis(A,C)>cal(nxt(p1),p2,D,B)+dis(A,D))			L=mid1;		else			R=mid2;	}	ret=min(ret,cal(nxt(p1),p2,C,B)+dis(A,C));	return ret;}int main(){	init();	while(scanf("%d%d%d%d",&a1,&a2,&a3,&a4),a1+a2+a3+a4!=-4)	{		if(a2>90&&a4>90)			a2=180-a2,a4=180-a4;		if(a2>a4)		{			swap(a1,a3);			swap(a2,a4);		}		getpoint(a1,a2,A);		if(a4>90)			getpoint(a3,180-a4,B);		else			getpoint(a3,a4,B);		printf("Case %d: %.3lf\n",++ca,cal(a1/90,a3/90,A,B));			}	return 0;}