Problem 144

R Project Euler

内部が鏡になっている楕円型の装置内に入射された光が特定のポイントに戻ってくるまでに何度反射するのか、という問題。
楕円は4x^2 + y^2 = 100という方程式で表現され、円周上のある点における接線は-4\frac{x}{y}で表される。
おおまかな流れだけ。

  1. 反射光の式を求める
    1. 逆正接関数で入射光と楕円の接線のx軸に対する角度を求める
    2. 計算した角度から反射光の角度を求める
    3. 反射光の正接を求め、切片も計算する
  2. 反射光と楕円が次に交わる点を求める
    1. 4x^2 + y^2 = 100とy=ax+bを連立させてyかxを消去すると2次方程式になる
    2. 解の一つが反射光の始点だから、終点の座標は解と係数の関係から計算できる
  3. 反射光の終点が出口に達していれば終了、そうでなければ反射光をあらたな入射光として1に戻る
## 反射光の傾きと切片
ref <- function(points){
  x1 <- points[1]; y1 <- points[2]; x2 <- points[3]; y2 <- points[4]
  slope <- tan(2*atan(-4*x2/y2) + pi - atan2(y2-y1, x2-x1))
  intercept <- y2 - slope*x2
  c(slope, intercept)
}

## 鏡と光の交点
imp <- function(par){
  ## a:slope b:intercept (x,y):start point
  a <- par[1]; b <- par[2]; x <- par[3]; y <- par[4]
  x2 <- (-(a*a + 4)*x - 2*a*b) / (4+a*a)
  y2 <- a*x2 + b
  c(x2, y2)
}

## 次の4組
next.set <- function(points){
  c(points[3:4], imp(c(ref(points), points[3:4])))
}
  
points <- c(0.0, 10.1, 1.4, -9.6)
count <- 0
repeat{
  if(abs(points[3]) <= 0.01 && points[4] > 0) break
  points <- next.set(points)
  count <- count + 1
}
count