はじめに

CodeIQというところでプログラミング問題「ツイン・トライアングル」に解答しましたので、ネタバレです。

codeiq.jp

問題

今回は図形の絡む問題でした。( 中学生位の時を思い出して懐かしいですね )
入力Lに対し、以下で定義されるf(L)を計算し出力する、というものです。

直角を挟む2辺a,b ( a,b,は自然数 ) に対し、a≦b≦Lかつ、次の条件を満たす直角三角形の個数をf(L)とする。
元の直角三角形に相似で、元の直角三角形の外接円が内接する直角三角形に対し、直角を挟む2辺a',b'も自然数であること。

ちょっと文章だと分かりにくいですね。以下の解説の中で図解します。

解説

数学的な考察

図形の性質

まずは図形的な性質を考察する必要があります。問題の図に、色々追加してみたのが次の図です。元の三角形 ( 内側 ) と外の三角形は a:b=a':b' の相似形です。

f:id:ange1:20151202182028p:plain

注目すべきは、2つの三角形の間にある円、一方にとっては外接円ですが、もう一方にとっては内接円であるということです。
そこで、元の三角形の内接円を追加して考えます。そうすると、三角形だけでなく、2つの円に関しても同じ相似比の相似形であることが分かります。

そのため、元の三角形の斜辺cを導入すると、

元の三角形の外接円の半径 $R=\frac{c}{2}$
元の三角形の内接円の半径 $r=\frac{ab}{a+b+c}$
相似比 $e=\frac{a'}{a}=\frac{b'}{b}=\frac{R}{r}$

となりますから、これらから

$e=\frac{a'}{a}=\frac{b'}{b}=\frac{c(a+b+c)}{2ab}$

であることが分かります。

自然数条件

さて。問題の条件は、a,b,a',b'が全て自然数であること、でした。
ここで導入した斜辺c ( $c=\sqrt{a^2+b^2}$ ) も自然数であることが分かります。もし自然数でないとすると無理数になってしまい、有理数である $\frac{a'}{a}, \frac{b'}{b}$ との整合性がとれなくなるからですね。

そのため、図形的な性質が分かった後は、自然数の問題として処理することになります。ただ、そのままだと扱いにくいので、最も単純化した形、すなわち a,b,c を最大公約数 g で割った x,y,z ( x,yは－ zを除いても－互いに素 ) で考えてみます。( ただ、都合により、a,b、x,yの大小はここでは決めずに話を進めます )

互いに素な x,y と z を3辺 ( 斜辺はz ) とする直角三角形の場合、以下のように2つの自然数m,nを導入してx,y,zを表せることが分かっています。

$x=m^2-n^2, y=2mn, z=m^2+n^2$
ただし、 $m\gt n$ 、 $m,n$ は互いに素、 $m-n$ は奇数

なお「a,b、x,yの大小はここでは決めずに」というのは、m,nの組み合わせによってx,yの大小が変わってしまうからですね。さてこれと、

$a=gx, b=gy, c=gz$
$e=\frac{c(a+b+c)}{2ab}$

とを組み合わせると、

$e=(m^2+n^2)\cdot\frac{g}{2n(m-n)}$

$m^2+n^2$ は、その後に来る分数の約分には寄与しませんから、 $a'=ea,b'=eb$ が自然数となるためには、gが $2n(m-n)$ の倍数であることが条件ということになります。結局、a,b は $g=2n(m-n)$ の時をベースにして、

$(a,b)=(k\cdot 2n(m-n)\cdot (m^2-n^2),k\cdot 2n(m-n)\cdot 2mn)$

と、(ベースの値)×k の形で表せるというのが、最終的な条件です。a,bの大小関係については有耶無耶のまま来ましたが、適宜計算した結果を並びなおして考えればO.K.なので特に気にしなくても問題ありません。

答えの求め方

結局のところ、m,nに応じた直角三角形のベースサイズが決まれば、それを整数倍に拡大したものが全て条件を満たすものです。長さ L を超えない範囲での個数がどうなるかは、直前の式におけるkの値の数の合計、

$f(L)=\sum_{m\gt n}^{}\lfloor\frac{L}{max( 2n(m-n)\cdot (m^2-n^2),2n(m-n)\cdot 2mn )}\rfloor$ ( ただし、m-nは奇数、m,nは互いに素 )

つまり、m,nのそれぞれの組み合わせに対して整数除算を行った結果を足していくことになります。maxを持ち出しているのは、どちらの値の方が大きいかが分からないためで、ここで $a\leq b$ の条件を吸収していることになります。
ただ、m,n だとちょっと分かりづらい ( 特に max の判断 ) というのもあるので、 $i=n, j=m-n$ でi,jを使った表現に置き換えてみます。
そうすると、max の判断も簡単になって、

$f(L)=\sum_{j\gt\sqrt{2}\cdot i}^{}\lfloor\frac{L}{2ij^2(2i+j)}\rfloor+\sum_{j\lt\sqrt{2}\cdot i}^{}\lfloor\frac{L}{4i^2j(i+j)}\rfloor$ (ただし、jは奇数、i,jは互いに素 )

これがf(L)の最終形となります。i,jの間の不等号が＜と＞で、≧ や ≦ が出てこないことに注意。これは、整数条件を満たすためには a=b が成立しないことに対応します。

実装

概要

さて。上記Σによる式について、i,jの範囲は特に決めていませんでしたが、Lの割り算の結果が1を下回る領域は答えに影響しませんから、「分母がLを超えない範囲」だけ対象にすれば十分です。グラフにすると次の図のようになり、格子点数としては $O(\sqrt{L})$ になります。

f:id:ange1:20151202204513p:plain

今回は速度重視のため、C ( C99 ) で実装します。ちなみに、向かなさそうなので golf はやりません。

naiveな実装

ということで、naiveな実装がこれになります。

各格子点での「互いに素」の判定は、ユークリッドの互除法を実装した関数 coprime で行っています。そのため、計算量としては多分 $O(\sqrt{L}\log L)$ となります。

naiveながらもそれなりの速度はあって、CodeIQの実行君でも、L=10¹⁴を0.6秒で処理できます。

チューニング

さて、上記naiveな実装では、互いに素の判定が処理時間の実に8～9割を占めます。…これは何か無駄な気がするので、何とかすれば高速化につながるのではないか…ということで考えてみます。

まず、判定を全格子点で行う必要があるのか? というところを考え直します。
例えば (i,j)=(15,23) という互いに素な組があったとすれば、jの値を30ずつ ( jは奇数であることに注意 ) 変化させた、(15,53),(15,83),… も互いに素です。つまり、ある数 ( この例では15 ) に対して、mod の値が同じ数字 ( 23,53,83,… ) の判定は、代表値の分の1度で済ませることができる、ということです。これにより、判定回数を減らすことができます。

もう一つ、そもそもユークリッドの互除法は必要か? というところを考え直します。
なぜならば、確かにユークリッドの互除法は高速なアルゴリズムではあるものの、除算を多用する比較的重い処理だからです。

そこで考えてみると、今回、ある特定の(i,j)の組で互いに素かどうかを判定する、という処理を繰り返さなくとも、i または j を固定した時に、互いに素になる数の一覧が得られれば十分であることに気づきます。
それであれば、関連する素因数の倍数にマークをつけていけば済む話で、マークが付かずに残った数がその「互いに素になる数の一覧」です。
問題は、一覧を保持するのに必要なメモリ量ですが、改めてグラフを見直してみると、赤の領域は i が lim1 未満の範囲で、青の領域は j が lim2 未満の範囲で調べればよいので、 $O(L^{1/4})$ ( Lの4乗根 ) のメモリ量で済みます。L が 64bit整数の上限一杯だとしても、その4乗根なら100Kも行きませんから問題ありません。

f:id:ange1:20151202210044p:plain

ということで、これらの高速化を盛り込んだ実装が次になります。

互いに素な数の一覧を調べるのはmark関数で行います。一旦markを実行すれば、互いに素かどうかの判断は単なるメモリアクセスで済みます。
なお、iは偶数・奇数どちらもありえますが、jは必ず奇数となりますから、例えば 15 に対して互いに素な一覧 1,2,4,7,8,11,13,14 が得られたら、

j=15 を固定して、i=1+15k,i=2+15k,4+15k,7+15k,… のケースを処理 ( グラフ中青い領域 )
i=15 を固定して、j=1+30k,j=17+30k,19+30k,7+30k,… のケースを処理 ( グラフ中赤い領域、以下同様 )
i=30 を固定して、j=1+30k,j=17+30k,19+30k,7+30k,… のケースを処理
i=60 を固定して、j=1+30k,j=17+30k,19+30k,7+30k,… のケースを処理
…

というように、i=奇数×(2のべき乗) のケースに対しても、その一覧を使いまわします。

実行君で L=10¹⁴ の処理時間が 0.13sまで短縮され、L=10¹⁵ も 1s に収まるようになりました。markの計算量は不明ですが、実行時間としてほとんど無視できるので、全体の計算量は $O(\sqrt{L})$ でかなりの高速化になったと言って良いでしょう。

細かいチューニング

ここから先は細かい話ですが…。

先ほどのチューニング、効果は絶大なのですが1つ気になる点があります。それは、従来 i-j の2重ループで実装していたところを i-iの余り一覧-j もしくは j-jの余り一覧-i という3重ループに変えたことで、終了判定、つまり $L \ge 4i^2 j(i+j), L \ge 2i j^2 (2i+j)$ が偽になる回数が増えているということで、i,jの4次式 ( 実際は3次式 ) の計算がそれだけ無駄になっているのではないかというところです。
この点に対しては、予め上限値を求めておくことで対処することができます。ループの終了判定が単純な比較で済むからです。
上限値については、 $i^2 (i+j)=\frac{L}{4j}$ あるいは $j^2(j+2i)=\frac{L}{2i}$ という、いずれも $x^3+ax^2=b$ という形の3次方程式を解くことで求められますので、適当な初期値を与えて $t\leftarrow \frac{2t^3+at^2+b}{3t^2+2at}$ を繰り返すニュートン法で計算できます。

…試しに実装すると、確かに効果はあるのですがかなり微妙です。3次式の計算は減っても、ループ終了判定回数が減るわけではないので、効果が薄いのでしょう。ただ、だからといって無駄になるわけではありません。他のチューニングの余地が生まれるからです。

それで、どこに目をつけるか、ですが…。
まず1つ目は3重ループの最内、同じ余りを持つ数に関するループです。i もしくは j の値が小さいうちは、同じ余りを持つ数も沢山ありますが、大きくなってくると少なくなり、ついには1つ2つしかない ( あるいはそもそもない ) という状況が生まれます。

f:id:ange1:20151203091858p:plain

そうするとループで処理すること自体が無駄です。で、そもそも格子点がない範囲は除外、ある場合でも繰り返し回数が既知の極短ループということで、手動でループ展開してしまおう、となります。

もう1つは、そもそも L÷… の形式の割り算が必要か? というところです。実は、i,jがある程度大きくなると、それに対応する三角形の数は1つ、つまりL÷…の割り算の結果は1になります。グラフでいうと、外側の実線と内側の点線とで挟まれた領域です。

f:id:ange1:20151202235224p:plain

この領域は思ったより大きく、全体の3割 ( $1-\frac{1}{\sqrt{2}}$ ) 程度を占めます。つまり、3割位は無駄な割り算を行っている、ということです。Lの値を半分にして3次方程式を解けば内側の点線に対応する上限も求められますから、この領域を特定することができ、その範囲での割り算を単純なインクリメントに置き換えることができます。
もっと追及すると、格子点の数も直接割り算で計算できるので、格子点毎のループをなくすことができます。なお、この割り算についても、有効な格子点が少なくなる領域では不要になります。

ただ1つ注意する点があります。それは3次方程式の求解の厳密性です。
L÷…の割り算を行うのであれば、範囲外の値で計算しても0が出るだけで答えに影響は出ないため、上限はザックリ決めても問題ありませんが、「三角形の数1個に対応する」格子点を数えて答えに足しこむ場合だと、厳密な上限値、つまり整数の範囲での3次方程式の最適解を導く必要があります。
ここでニュートン法で収束させただけでは、整数としての最適解に到達する保証が実はありません。そのためベタですが、実際に3次式の値を確かめて調整を行い、最適解を求めるということをします。

そういった点を実装したのがこちらです。後、64bit整数が必要ない箇所をこまめに int に替えることでも多少の性能改善効果を得ています。

CodeIQの実行君で、L=10¹⁴ が 0.07s と naive版の8倍以上の速度、L=10¹⁶ すら1s以内に収まるようになりました ( 0.89s )。

まとめ

数学的には $f(L)=\sum_{j\gt\sqrt{2}\cdot i}^{}\lfloor\frac{L}{2ij^2(2i+j)}\rfloor+\sum_{j\lt\sqrt{2}\cdot i}^{}\lfloor\frac{L}{4i^2j(i+j)}\rfloor$ (ただし、jは奇数、i,jは互いに素 )
C ( C99 ) で3種類の実装を行い、CodeIQの実行君で実行時間を測定した。
ユークリッドの互除法により都度「互いに素」を判定するnaiveな実装では、計算量が多分 $O(\sqrt{L}\log{L})$ 、メモリは特に不要、L=10¹⁴ が0.6s
ユークリッドの互除法をやめて、同じ余りになる数で「互いに素」の判定を使いまわすようにチューニングした場合、計算量が多分 $O(\sqrt{L})$ 、メモリ量 $O(L^{1/4})$ 、L=10¹⁴ が0.13s。
ニュートン法により求めた領域の上限を活用して更に細かくチューニングすることで、計算量・メモリ量のオーダーはそのままながら、L=10¹⁴ が0.07s、1s以内にL=10¹⁶を処理することも可能になった。

終わりに

L=10¹⁶ の1s切りができたので、とりあえず満足しました。…これ以上追及するのは面倒というのもありますが。
たまには実行時間の短縮を目指すというのも面白いですね。

ange1のブログ

「ツイン・トライアングル」問題解答 ( CodeIQ )