/******************************************************************************
SuperCon '97 本選課題解答プログラム
copyright(c)1997 by まい, shn & take all rights reserved.
for UNIX/CRAY
******************************************************************************/
#define MULTI_TASKING
/*
誤差レポート by Shin Adachi.
誤差範囲は、計算を重ねるごとに蓄積されていくので、
現在のように0.001等の大きな数字を設定していたら、
計算の初めの方で誤差がでる危険性が大きい。
そこで、いくら誤差が積み重なったか計算しなければ
ならないのだが、たとえばaとbの掛け算をするとき
両者にはそれぞれ+e〜-eの誤差があることになる
よって
(a+e)(b+e) = ab+e(a+b)+e^2 // 誤差が最大の場合
(a-e)(b-e) = ab-e(a+b)+e^2 // 誤差が最小の場合
これによって分かる通り、掛け算をするだけで誤差が
a+b倍に広がっている(eの2乗は小さすぎるので無視
できる)
また、これをεの比率で見てみると浮動小数点aは
a(1-e) < a < a(1+e)
の範囲となる。
よって、例えばa * bの計算をすると
新たな数abの誤差の比率は
ab(1-2e+e^2) < ab < ab(1+2e+e^2)
となる。これに安全策を含めて、一回の乗算で(2.1e)
の誤差が出たの考えるのがおすすめ。
後、桁落ちというのがあって、例えば
100 - 99.9
という計算をさせると、その答えはやばくなるので要注意
*/
#include
#include
#include
typedef struct
{
double x,y; /* 始点 */
} vector;
typedef struct
{
vector p[3];
} triangle;
#define MAX 124
#define lengthC 1000
#define fabs(x) ((x)<0?-(x):(x))
#define EPSILON 0.0001
/* 新比較マクロ */
#define FEqual(x,y) (fabs((x)-(y))<=EPSILON)
#ifndef min
#define min(x,y) ((x)<(y)?(x):(y))
#endif
#define minimum(x,y,z) (min(min((x),(y)),(z)))
#ifndef max
#define max(x,y) ((x)>(y)?(x):(y))
#endif
#define maximum(x,y,z) (max(max((x),(y)),(z)))
vector p[3]; /* 三点の座標 */
#define isin(x) sin((x)*M_PI/180)
#define icos(x) cos((x)*M_PI/180)
#define RIGHT -1
#define LEFT 1
typedef struct
{
int nlv; /* 再起レベル */
char seq[MAX+2]; /* 折り返した点=辺の履歴ふくれき */
vector L[MAX+2][MAX+2];
vector R[MAX+2][MAX+2];
int Llen[MAX+2];
int Rlen[MAX+2];
triangle Tri[MAX+2];
int answers;
int startNLV;
} SolveData;
void Solve(void);
/* 解く */
int Remake(vector *V,int *len,vector r,int lr);
/* 凸包の更新。V:凸包のデータ,len:データの長さ,r:新しい点,lr:どちらの凸包か */
/* lr:右凸包の時 -1,左凸包の時 1 */
int Cheque(vector p0,vector p1,vector p4,vector *V,int len);
/* 左右の凸包が重ならないかチェック */
/* 更新してできた範囲に重なったら真、重ならなかったら偽 */
int GetLocation(vector *L,int Llen,vector *R,int Rlen);
/* 凸包から範囲を求めて、その中心を答えの位置として返す。 */
/* 範囲の広さが1未満か、線が引けない場合、-1を返す */
int sign(double x)
{
return ((x>=EPSILON && x<=EPSILON)?0:((x<0)?-1:1));
}
/* 渡された変数の符号を返す */
double vectlr(vector p,vector q,vector r);
/* ベクトルの左右どちらに点rがあるか */
void mirror(vector p,vector q,vector *r);
/* ベクトルに対して対称に移動 */
int intriangle(vector a,vector b,vector c,vector r);
/* 三角形abcの中に点rがあるかどうか調べる。あれば真 */
int GetSeqCycle(char *seq);
/* GET SEQUANCE CYCLE */
void main(int ac,char *av[])
{
int degrees[3]; /* α、β、γ */
double length[3]; /* 辺a,b,cの長さ */
degrees[0]=atoi(av[1]);
degrees[1]=atoi(av[2]);
if( (sizeof( SolveData ) % sizeof( int )) != 0 )
printf( "ERROR" );
if(degrees[0]<1 || degrees[0]>179)
{
printf("Illegal input!\n");
return;
}
if(degrees[1]<1 || degrees[1]>179)
{
printf("Illegal input!\n");
return;
}
degrees[2]=180-degrees[0]-degrees[1];
length[0]=isin(degrees[0])/isin(degrees[2])*lengthC;
length[1]=isin(degrees[1])/isin(degrees[2])*lengthC;
length[2]=lengthC;
p[0].x=0;
p[0].y=0;
p[1].x=lengthC;
p[1].y=0;
p[2].x=length[1]*icos(degrees[0]);
p[2].y=length[1]*isin(degrees[0]);
printf("%d %d\n",degrees[0],degrees[1]);
Solve();
printf("***\n");
}
/*
pqか、その延長に対して、rがど位置にあるか
負=右/0=直線上/正=左
*/
/*
double vectlr(vector p,vector q,vector r)
{
return (q.x-p.x)*(r.y-p.y)-(r.x-p.x)*(q.y-p.y);
}
*/
#define vectlr(p,q,r) ((q.x-p.x)*(r.y-p.y)-(r.x-p.x)*(q.y-p.y))
void mirror(vector p,vector q,vector *r)
{
vector c;
double d;
d=(q.x-p.x)*(q.x-p.x)+(q.y-p.y)*(q.y-p.y);
c.x=r->x+2.0*(q.y-p.y)*((p.x-r->x)*(q.y-p.y)-(p.y-r->y)*(q.x-p.x))/d;
c.y=r->y+2.0*(q.x-p.x)*((p.y-r->y)*(q.x-p.x)-(p.x-r->x)*(q.y-p.y))/d;
*r=c;
return;
}
/*
Remake 凸包
*/
int Remake(vector *V,int *len,vector r,int lr)
{
/*
result=vectlr(r,V[x],V[x-1]);
sign(result)==lrが成立しなくなったところと結ぶ、
*/
int i;
#pragma _CRI inline sign
/*#pragma _CRI inline vectlr*/
for(i=(*len)-1;i>0 && sign(vectlr(r,V[i],V[i-1]))==lr;i--);
*len=i+2;
V[i+1]=r;
return *len;
}
/*
intriangle
*/
int intriangle(vector a,vector b,vector c,vector r)
{
double alr,blr,clr;
/*#pragma _CRI inline vectlr*/
clr=vectlr(a,b,r);
/*#pragma _CRI inline vectlr*/
alr=vectlr(b,c,r);
/*#pragma _CRI inline vectlr*/
blr=vectlr(c,a,r);
if(clr==0 && alr==0 && blr==0)
{
return
maximum(a.x,b.x,c.x)>=r.x &&
minimum(a.x,b.x,c.x)<=r.x &&
maximum(a.y,b.y,c.y)>=r.y &&
minimum(a.y,b.y,c.y)<=r.y;
}
return (clr>=0 && blr>=0 && alr>=0) || (clr<=0 && blr<=0 && alr<=0);
}
/*
Cheque
凸包が重ならないかどうかのチェック
*/
int Cheque(vector p0,vector p1,vector p4,vector *V,int len)
{
int i;
#pragma _CRI inline intriangle
for(i=0;ilbuf[i])?l:i;
for(i=1;i=0) return -1;
i=(right+left)/2;
return i;
}
/*
Solve the problem
*/
#define TNOCPU 8
short Solve2( int curSEQ, SolveData *data, SolveData *newData );
void StartLoop( SolveData *data );
void HiCopy( SolveData *src, SolveData *dst );
void Solve(void)
{
SolveData solveDatas[10];
int usedCPUs = 1;
int foo;
int curNLV, curData, curSEQ;
int curNewData;
SolveData newData[2], baz;
int check;
solveDatas[0].nlv=0;
/* 再起レベル0 */
solveDatas[0].seq[0] = 2;
/* 最初に折り返せないのは、c */
solveDatas[0].Tri[0].p[0]=p[0];
solveDatas[0].Tri[0].p[1]=p[1];
solveDatas[0].Tri[0].p[2]=p[2];
/* 最初の三角形を設定 */
solveDatas[0].L[0][0]=p[0];
solveDatas[0].R[0][0]=p[1];
solveDatas[0].Llen[0]=1;
solveDatas[0].Rlen[0]=1;
/* 右と左の点を設定 */
solveDatas[0].answers = 0;
while( usedCPUs<=TNOCPU )
{
for( curData=usedCPUs-1;curData>=0;curData-- )
{
for( curSEQ=0, curNewData=0;curSEQ<3;curSEQ++ )
{
if( curSEQ!=solveDatas[curData].seq[solveDatas[curData].nlv] )
{
if( Solve2( curSEQ, &solveDatas[curData], &newData[curNewData] )!=0 )
{
for( check=0;check<3;check++ )
if( check!=newData[curNewData].seq[newData[curNewData].nlv] )
if( Solve2( check, &newData[curNewData], &baz )!=0 )
break;
if( check==3 )
;
else
curNewData++;
}
}
}
switch( curNewData )
{
case 0:
if( curData==0 )
{
printf( "no answer" );
return;
}
usedCPUs--;
for( foo=curData;fooseq[data->nlv+1] = curSEQ;
if(data->seq[data->nlv+1]==3)
{
/* 現在の再起が最後まで終了 */
return 0;
}
if(data->seq[data->nlv]!=data->seq[data->nlv+1]) /* 一つ前に折り返した点でなければ */
{
/*
再起準備処理
*/
int lr;
int i;
int nlv = data->nlv;
#pragma _CRI inline HiCopy
HiCopy( data, newData );
/* *newData = *data;*/
/* 曲がるのが、右か左か */
if((data->nlv&1))
{
/*nlv==Odd*/
lr=(((data->seq[data->nlv]-1+3)%3)==data->seq[data->nlv+1])?RIGHT:LEFT;
}
else
{
/*nlv==Even*/
lr=(((data->seq[data->nlv]+1+3)%3)==data->seq[data->nlv+1])?RIGHT:LEFT;
}
/* 凸包をコピー */
for(i=0;iLlen[nlv];i++)
newData->L[nlv+1][i]=data->L[nlv][i];
for(i=0;iRlen[nlv];i++)
newData->R[nlv+1][i]=data->R[nlv][i];
newData->Llen[nlv+1]=data->Llen[nlv];
newData->Rlen[nlv+1]=data->Rlen[nlv];
/*
鏡像処理、
*/
newData->Tri[nlv+1]=data->Tri[nlv];
nlv++;
newData->nlv = nlv;
switch(newData->seq[newData->nlv])
{
case 0:
#pragma _CRI inline mirror
mirror(newData->Tri[nlv].p[1],newData->Tri[nlv].p[2],&newData->Tri[nlv].p[0]);
break;
case 1:
#pragma _CRI inline mirror
mirror(newData->Tri[nlv].p[0],newData->Tri[nlv].p[2],&newData->Tri[nlv].p[1]);
break;
case 2:
#pragma _CRI inline mirror
mirror(newData->Tri[nlv].p[0],newData->Tri[nlv].p[1],&newData->Tri[nlv].p[2]);
break;
}
/* 凸包を更新 */
if(lr==LEFT)
{
vector p0;
vector p1;
vector p4;
p0=newData->Tri[nlv-1].p[newData->seq[nlv-1]];
p1=newData->R[nlv][newData->Rlen[nlv]-1];
#pragma _CRI inline Remake
Remake(newData->R[nlv],&newData->Rlen[nlv],newData->Tri[nlv-1].p[newData->seq[nlv-1]],RIGHT);
p4=newData->R[nlv][newData->Rlen[nlv]-2];
#pragma _CRI inline Cheque
if(Cheque(p0,p1,p4,newData->L[nlv],newData->Llen[nlv]))
{
/* 重なっている */
/* 再起はできない */
return 0;
}
}
else
{
vector p0;
vector p1;
vector p4;
p0=newData->Tri[nlv-1].p[newData->seq[nlv-1]];
p1=newData->L[nlv][newData->Llen[nlv]-1];
#pragma _CRI inline Remake
Remake(newData->L[nlv],&newData->Llen[nlv],newData->Tri[nlv-1].p[newData->seq[nlv-1]],LEFT);
p4=newData->L[nlv][newData->Llen[nlv]-2];
#pragma _CRI inline Cheque
if(Cheque(p0,p1,p4,newData->R[nlv],newData->Rlen[nlv]))
{
/* 重なっている */
/* 再起はできない */
return 0;
}
}
/* 基本サイクルかどうかチェック */
if( newData->seq[newData->nlv]==2 &&
newData->nlv!=2 &&
(newData->nlv%2)==0 &&
(MAX % newData->nlv)==0 &&
FEqual(p[0].y-newData->Tri[nlv].p[0].y,p[1].y-newData->Tri[nlv].p[1].y) &&
FEqual(p[1].y-newData->Tri[nlv].p[1].y,p[2].y-newData->Tri[nlv].p[2].y) )
{
return 0;
}
return 1;
}
else
return 0;
}
void StartLoop( SolveData *data )
{
for(data->seq[data->nlv+1]=0;;data->seq[data->nlv+1]++)
{
if(data->seq[data->nlv+1]==3)
{
/* 現在の再起が最後まで終了 */
if(data->nlv==data->startNLV) /* 一番最初のループなら */
{
break;
}
data->nlv--;
continue;
}
if(data->seq[data->nlv]!=data->seq[data->nlv+1]) /* 一つ前に折り返した点でなければ */
{
/*
再起準備処理
*/
int lr;
int i;
/* 曲がるのが、右か左か */
if((data->nlv&1))
{
/*nlv==Odd*/
lr=(((data->seq[data->nlv]-1+3)%3)==data->seq[data->nlv+1])?RIGHT:LEFT;
}
else
{
/*nlv==Even*/
lr=(((data->seq[data->nlv]+1+3)%3)==data->seq[data->nlv+1])?RIGHT:LEFT;
}
/* 凸包をコピー */
for(i=0;iLlen[data->nlv];i++)
data->L[data->nlv+1][i]=data->L[data->nlv][i];
for(i=0;iRlen[data->nlv];i++)
data->R[data->nlv+1][i]=data->R[data->nlv][i];
data->Llen[data->nlv+1]=data->Llen[data->nlv];
data->Rlen[data->nlv+1]=data->Rlen[data->nlv];
/*
鏡像処理、
*/
data->Tri[data->nlv+1]=data->Tri[data->nlv];
data->nlv++;
switch(data->seq[data->nlv])
{
case 0:
#pragma _CRI inline mirror
mirror(data->Tri[data->nlv].p[1],data->Tri[data->nlv].p[2],&data->Tri[data->nlv].p[0]);
break;
case 1:
#pragma _CRI inline mirror
mirror(data->Tri[data->nlv].p[0],data->Tri[data->nlv].p[2],&data->Tri[data->nlv].p[1]);
break;
case 2:
#pragma _CRI inline mirror
mirror(data->Tri[data->nlv].p[0],data->Tri[data->nlv].p[1],&data->Tri[data->nlv].p[2]);
break;
}
/* 誤差の更新 */
/* 凸包を更新 */
if(lr==LEFT)
{
vector p0;
vector p1;
vector p4;
p0=data->Tri[data->nlv-1].p[data->seq[data->nlv-1]];
p1=data->R[data->nlv][data->Rlen[data->nlv]-1];
#pragma _CRI inline Remake
Remake(data->R[data->nlv],&data->Rlen[data->nlv],data->Tri[data->nlv-1].p[data->seq[data->nlv-1]],RIGHT);
p4=data->R[data->nlv][data->Rlen[data->nlv]-2];
#pragma _CRI inline Cheque
if(Cheque(p0,p1,p4,data->L[data->nlv],data->Llen[data->nlv]))
{
/* 重なっている */
/* 再起はできない */
if(data->nlv==data->startNLV)
return;
data->nlv--;
continue;
}
}
else
{
vector p0;
vector p1;
vector p4;
p0=data->Tri[data->nlv-1].p[data->seq[data->nlv-1]];
p1=data->L[data->nlv][data->Llen[data->nlv]-1];
#pragma _CRI inline Remake
Remake(data->L[data->nlv],&data->Llen[data->nlv],data->Tri[data->nlv-1].p[data->seq[data->nlv-1]],LEFT);
p4=data->L[data->nlv][data->Llen[data->nlv]-2];
#pragma _CRI inline Cheque
if(Cheque(p0,p1,p4,data->R[data->nlv],data->Rlen[data->nlv]))
{
/* 重なっている */
/* 再起はできない */
if(data->nlv==data->startNLV)
return;
data->nlv--;
continue;
}
}
/* 基本サイクルかどうかチェック */
if( data->seq[data->nlv]==2 &&
data->nlv!=2 &&
(data->nlv%2)==0 &&
(MAX % data->nlv)==0 &&
FEqual(data->Tri[0].p[0].y-data->Tri[data->nlv].p[0].y,data->Tri[0].p[1].y-data->Tri[data->nlv].p[1].y) &&
FEqual(data->Tri[0].p[1].y-data->Tri[data->nlv].p[1].y,data->Tri[0].p[2].y-data->Tri[data->nlv].p[2].y) )
{
#pragma _CRI inline GetSeqCycle
if(data->nlv==MAX && GetSeqCycle(data->seq+1)==0)
{
char ans[MAX+2];
int location;
#pragma _CRI inline GetLocation
location=GetLocation(data->L[data->nlv],data->Llen[data->nlv],data->R[data->nlv],data->Rlen[data->nlv]);
if(location!=-1) /* 整数を挟んでいたとき */
{
data->answers++;
for(i=0;iseq[i+1]+'a';
ans[i] = 0;
#ifdef MULTI_TASKING
#pragma _CRI guard
#endif
printf("%s\n%d\n",ans,location);
#ifdef MULTI_TASKING
#pragma _CRI endguard
#endif
}
data->nlv--;
continue;
}
}
if(data->nlv==MAX)
{
/* printf( "hogehoge" );getchar();*/
data->nlv--;
continue;
}
/* if( nlv>=120 )
{printf( "%d, %d, %d\n", data->nlv, data->nlv, data->seq[data->nlv] );getchar();}*/
if(data->nlv==MAX-1)
data->seq[data->nlv+1] = 1;
else
data->seq[data->nlv+1] = -1;
}
}
}
int GetSeqCycle(char *seq)
{
int i, j;
/*
int flag=0;
*/
for(i=2; inlv=src->nlv;
dst->answers=src->answers;
dst->startNLV=src->startNLV;
#pragma _CRI parallel private( i, j ) shared( dst->seq, src->seq, dst->Rlen, src->Rlen, dst->Llen, src->R, dst->R, src->L, dst->L
#pragma _CRI novector
for(i=0;i<=src->nlv;i++)
dst->seq[i]=src->seq[i];
#pragma _CRI novector
for(i=0;i<=src->nlv;i++)
dst->Rlen[i]=src->Rlen[i];
#pragma _CRI novector
for(i=0;i<=src->nlv;i++)
dst->Llen[i]=src->Llen[i];
for(i=0;i<=src->nlv;i++)
#pragma _CRI novector
for(j=0;jRlen[i];j++)
dst->R[i][j].x=src->R[i][j].x;
for(i=0;i<=src->nlv;i++)
#pragma _CRI novector
for(j=0;jRlen[i];j++)
dst->R[i][j].y=src->R[i][j].y;
for(i=0;i<=src->nlv;i++)
#pragma _CRI novector
for(j=0;jLlen[i];j++)
dst->L[i][j].x=src->L[i][j].x;
for(i=0;i<=src->nlv;i++)
#pragma _CRI novector
for(j=0;jLlen[i];j++)
dst->L[i][j].y=src->L[i][j].y;
for(i=0;i<=src->nlv;i++)
#pragma _CRI novector
for(j=0;j<3;j++)
dst->Tri[i].p[j]=src->Tri[i].p[j];
*/
*dst = *src;
}