AOJ 2506 - Operation training for BYDOL

この問題では、「弾が点Aから壁上の点Pで反射して点Bに達する」という条件を満たす点Pを求める必要があります。これは二分探索や三分探索で求めることもできるらしいですが、今回は探索なしで求める方法を述べます。

以後、壁の中心座標を(0,0)、半径をrとします。また、点A,B,Pの位置ベクトルを、それぞれ\vec{a}=(a_x,a_y),\ \vec{b}=(b_x,b_y),\ \vec{p}とします*1。さらに、\vec{p}=r\vec{u}を満たす単位ベクトル\vec{u}=(x,y)を定義することにします。単位ベクトルなので、当然x^2+y^2=1となります。

f:id:climpet:20130321173642p:plain
図1のように、ベクトル\vec{v}=\vec{AP}が点Pで反射した後のベクトル\vec{w}がどのような計算式で求められるか考えてみましょう*2。ここで、ベクトル\vec{v}を、\vec{u}に平行な成分\vec{v_o}と、垂直な成分\vec{v_n}に分けてみることにします。ベクトル\vec{w}についても同様です。このとき、図から\vec{v}=\vec{v_o}+\vec{v_n},\ \ \ \vec{w}=\vec{w_o}+\vec{w_n},\ \ \ \vec{w_o}=-\vec{v_o},\ \ \ \vec{w_n}=\vec{v_n}が成り立ちます。また、\vec{v_o}内積を用いて\vec{v_o}=(\vec{v}\cdot\vec{u})\vec{u}と表すことができます。これらの式から、\vec{w}=\vec{v}-2(\vec{v}\cdot\vec{u})\vec{u}が成り立ちます*3

反射後のベクトル\vec{w}は、ベクトル\vec{PB}=\vec{b}-\vec{p}=\vec{b}-r\vec{u}と同じ向きでなければなりません。従って、\vec{w}\times\vec{PB}=0が成り立ちます(2次元ベクトルにおける外積については、ここを参照してください)。これらの式から色々と計算し、さらにy^2=1-x^2を用いてyを除去すると、次の方程式が得られます。
4(c_1^2+c_2^2)x^4-4(c_1c_4+c_2c_3)x^3+(c_3^2+c_4^2-4c_1^2-4c_2^2)x^2+(2c_1c_4+4c_2c_3)x+(c_1^2-c_3^2)=0
ただし
c_1=a_xb_y+a_yb_x
c_2=a_xb_x-a_yb_y
c_3=r(a_x+b_x)
c_4=r(a_y+b_y)
とします。
これは4次以下の方程式です*4。幸運にも、4次以下の方程式については代数的解法が知られているため、厳密解を求めることもできます。それ以外の方法としては、ベアストウ法と呼ばれる数値計算法が知られており、この方法でも解を求めることができます*5。いずれの手法でも、解は複素数として得られます。
これを解くことで得られた、-1\le x\le 1を満たす*6実数解xに対し、y=\sqrt{1-x^2}とします。このとき、P(rx,ry)またはP(rx,-ry)が、反射点の候補となります。

このようにして求められた反射点候補の中には、実際には反射点でないような点も含まれている可能性があります。そのため、本当に反射点であるのかどうか確認する必要があります。ここで、Pが反射点であれば\vec{w}\vec{PB}は平行でかつ同じ向きとなるため、\vec{w}\times\vec{PB}=0*7および\vec{w}\cdot\vec{PB}\gt 0の両方が成り立ちます。これらを計算することで、反射点でない候補を除外することができます。

f:id:climpet:20130321173708p:plain
ただし、この判定を行っても、除外できない非反射点が存在します。それは図2の点P'のように、「円の内側で反射して点Bに達する」というものです。これを除外するためには、「線分APまたはBPと円の交点が、端点以外に存在するかどうか」を判断すればよいです。図2では、反射点である点Pについては、線分APおよびBPと円の交点は、端点以外に存在しません。一方、非反射点である点P'については、線分AP'上に交点Qが存在します(線分BPについても同様)。したがって、点P'は反射点ではないと判断できます。
この問題では、どの道すべての壁との交点判定を行う必要があるため、特に手間が増えるということはないでしょう。

以上が、反射点を探索なしで求める方法です。ここからグラフを作って最小費用流に帰着させることになりますが、詳細は省略します。公式スライドなどを参照してください。



…記事書いてみて思ったのですが、u,\ v,\ wって微妙に見分けつきにくいですね。。。

*1:ベクトルはボールドイタリックにしたいのですが、はてなでの書き方が分かりませんでした…

*2:図と本文で書体が違いますが、ご了承ください

*3:この性質は、3次元の反射でも成り立ちます。例えばAOJ 1289では同じ方法により、3次元球面上での反射を考えることができます。

*4:係数は0になることもあるので、必ずしも4次方程式であるとは限りません。

*5:ただし、ベアストウ法は数値計算的な手法のため、誤差が生じます。以後の計算では、これにより生じうる誤差より大きな値をEPSとして用いる必要があります。

*6:誤差の影響で、x=1.00000000001などの値が解として得られる可能性があります。そのため、xの絶対値が1に極めて近い値については、1または-1として扱う必要があります。 また、やや乱暴ですが、x>1であればx=1として、x<-1であればx=-1として、それぞれ扱うこともできます。このようにしたとしても、反射のチェックで除外されるため、最終的な結果には影響を与えません。

*7:正確には、誤差を考慮するため、「外積の絶対値がEPS未満」という判定にする必要があります。