プロジェクトオイラー問144 - sinapusu2002 @ ウィキ

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20144
楕円型のwhite cellに入射したレーザー光線が出てくるまでの回数をこたえる問題。
ｘ＝0やｙ＝0での反射がないので場合わけを端折って後は反射をベクトルで計算していくだけです。

#include<stdio.h>
#include<algorithm>
#include<math.h>

double calcT(double dx,double dy,double x0,double y0){
	double a=4*dx*dx+dy*dy;
	double b=8*x0*dx+2*y0*dy;
	double c=(4*x0*x0+y0*y0-100);
	double d=sqrt(b*b-4*a*c);
	double t1=(-b+d)/(2*a);
	double t2=(-b-d)/(2*a);
	return std::max(t1,t2);
}

double cos2(double x0,double y0){
	double dx,dy;
	if(y0>0){
		dx=-4*x0/y0;
		dy=-1;
	}else{
		dx=4*x0/y0;
		dy=1;
	}
	double len=sqrt(dx*dx+dy*dy);
	double cosA=dx/len;
	double sinA=dy/len;
	return cosA*cosA-sinA*sinA;
} 
double sin2(double x0,double y0){
	double dx,dy;
	if(y0>0){
		dx=-4*x0/y0;
		dy=-1;
	}else{
		dx=4*x0/y0;
		dy=1;
	}
	
	double len=sqrt(dx*dx+dy*dy);
	double cosA=dx/len;
	double sinA=dy/len;
	return 2*cosA*sinA;
}


int main(){
	double x=1.4,y=-9.6;
	double dx=1.4,dy=-9.6-10.1;
	int c=0;
	while(1){
		double dx1,dy1,cos2A,sin2A;
		cos2A=cos2(x,y);
		sin2A=sin2(x,y);
		dx1=dx*cos2A+dy*sin2A;
		dy1=dx*sin2A-cos2A*dy;
		dy1=-dy1;
		dx1=-dx1;
		double t=calcT(dx1,dy1,x,y);
		dx=dx1;
		dy=dy1;
		x=x+dx*t;
		y=y+dy*t;
		c++;
		//printf("x=%lf y=%lf dx=%lf dy=%lf t=%lf\n",x,y,dx,dy,t);
		if((-0.01<=x)&&(x<=0.01)&&(y>0))break;
	}
	printf("%d\n",c);
}