「タンジェント・フラクション」問題解答 ( CodeIQ )

CodeIQというところでプログラミング問題「タンジェント・フラクション」に解答しましたので、ネタバレです。

codeiq.jp

問題

今回は次のような問題でした。

  • 入力される正の実数 $m\,(0.001\lt m\lt1)$ に対し、$F(m)$の値を出力する
  • $F(m)$とは $m\le\alpha+\beta$ となる実数 $\alpha,\beta$の中で、次の条件を満たす組の個数とする。
    • $0\lt\alpha\lt\beta\lt\frac{\pi}{2}$
    • $tan\alpha,\,tan\beta,\,tan(\alpha+\beta)$がいずれも単位分数 ( 自然数の逆数 )

例えば、$\alpha=arctan(\frac{1}{3})=0.3217\cdots,\,\beta=arctan(\frac{1}{2})=0.4636\cdots$に対して、$tan(\alpha+\beta)=1$ ですから、これは条件を満たしています。

なお、問題では$F(0.5)=1,\,F(0.3)=4,\,F(0.1)=15,\,F(0.01)=266$といった例が挙げられていました。

解説

条件の整備

まず、三角関数であるtanの加法定理として、

  • $tan(\alpha+\beta)=\frac{tan\alpha+tan\beta}{1-tan\alpha tan\beta}$

という性質がありますから、単位分数になる ( それぞれのtanの値が$\frac{1}{a},\,\frac{1}{b},\,\frac{1}{k}$と表せる ) という条件も加味すると、

  • $F(m)= |\{ (a,b)\in \mathbb{N}^2 | a\gt b,\,\exists k\in\mathbb{N}\,\,s.t.\,k\le \frac{1}{tan m},\, \frac{1}{k}=\frac{\frac{1}{a}+\frac{1}{b}}{1-\frac{1}{a}\cdot\frac{1}{b}}\}|$

というような条件を満たす自然数$(a,b)$の組として表現することができます。

ここで、条件中の等式を更に整理すると、

$$ \begin{eqnarray} \,&\,&\frac{1}{k}=\frac{\frac{1}{a}+\frac{1}{b}}{1-\frac{1}{a}\cdot\frac{1}{b}} \\ \Leftrightarrow &\,& ab-k(a+b)-1=0 \\ \Leftrightarrow &\,& (a-k)(b-k)=k^2+1 \\ \end{eqnarray} $$

つまり、$a,b$の組というのは、$k^2+1$の2個の自然数の積への分解に1対1に対応します。そこで、

  • $F(m)=\frac{1}{2}\sum_{k=1}^{n} d(k^2+1)\,\,\,( n=[\frac{1}{tan m}] )$

と、約数の個数の総和として求めることができます。( $d(x)$は$x$の正の約数の個数を表す関数です )

$\frac{1}{2}$をかけているのは、$a,b$の大小関係のためです。幸いなことに、$k^2+1$の形の数は平方数となり得ませんから、約数の個数が奇数になって割り切れない場合どうするか、なんてことは考える必要はありません。

naiveな実装

さて、$k^2+1$の約数については、$k$までを調べれば十分なので、$n=[\frac{1}{tan m}]$に対して、Rubyであれば次のような関数を実装することで計算することができます。

->n{
  (1..n).reduce(0){|s,k|
    b=k*k+1
    s+(1..k).count{|j|b%j==0}
  }
}

まあ、分かり易いですね。

もちろん、素数を予めリストアップしておいて、素因数分解を用いた約数の個数の計算を行えば、もう少し速い実装になります。

改良版

再度整理

とは言え、一つ一つ約数の個数を検証していくのは、何というか心理的に抵抗があります。他の方法が何か無いかを考えてみました。

数えるのは$k^2+1$の約数、それも$k$以下のものです。ここで、ある$d$が約数だったとします。そうすると $(k-d)^2+1$も同じ約数$d$を持つはずです。なぜなら、

  • $(k-d)^2+1=k^2+1-d(2k-d)\equiv k^2+1\,\,mod\,\,d$

であるからです。

ここで、例えば$27^2+1$は約数10を持ち、$27^2+1=10\times 73$ を満たします。これを (27,10,73) の数値3つ組で表すものとします。( 左端が$k$、真ん中が$d$、右端が$\frac{k^2+1}{d}$ )

これを元に$k\rightarrow k-d$と辿っていくと、(27,10,73)→(17,10,29)→(7,10,5)と一意に列を作ることができます。辿りついた最後の(7,10,5)はもはや$d\le k$を満たすものではないですが、これはしかし$d\le k$を満たす(7,5,10)と対になるものと見ることができます。(次の図を参照)

f:id:ange1:20170617222416p:plain

そうすると、(7,5,10),(7,10,5)の更に先も含めて、これらをノードと見て2分木が構成できることが分かります。

f:id:ange1:20170617222633p:plain

この時点で、もはや約数の個数を数える問題ではなく、木を探索して該当するノードを数える問題に変わっているのです。

しかも、どの約数に対応するノードであっても、上に辿っていけば必ずノード(0,1,1)に集約されることが分かります。つまり、(0,1,1)を根とする木を1つ考えるだけで、全体を網羅できてしまうのです。

f:id:ange1:20170617223927p:plain

どう数えるか

これで木として扱えることは分かりましたが、そうは言っても1つ1つノードを探索しているようでは、さっぱり速度が出ません。そこで、枝毎に効率よく数えることを考えます。

さて、探索条件は「$k$が指定の上限以下であること」です。ここで例として、$k$の上限50、ノード(3,5,2)以下のサブツリーを見てみます。

f:id:ange1:20170617225202p:plain

この中で、実際に数える対象のノードである、$d\le k$ となっているもの ( 図中楕円形のノード ) を、左下に伸びる ( 楕円形のノードだけで構成される ) 枝単位で見てみます ( 枝の中では$d$が共通となっています )。そうすると、それぞれの枝の中では、ノードの$k$の値が$d$を公差とする等差数列を構成します。なので、この枝の長さは、単純に割り算で求めることができるはずです。

残りの問題は木の分岐をどう辿るか、です。しかしここで、$k+d$の値が上限を超えるようなら、そこから伸びる枝はもはや対象外、無視することができます ( 図中赤の×が付いている枝 )。これによって効率的に枝刈りが行えます。

分岐する際は、ノードの右端の数の推移も考える必要があります。これは、単純に$\frac{k^2+1}{d}$で計算しても良いですが、実はここにも規則性があることが見て取れます。

例えば、先ほどの図の左端の枝の場合、分岐元の ( 楕円でない ) ノードも含めると右端の数字の推移は 2→13→34→65→… となっているのですが、隣接項同士の差を取れば、つまり階差数列を見れば、等差数列になっているからです。( 11→21→41→… )
この階差数列の公差は丁度$2d$です。
ちなみに、階差数列の初項、これは先頭ノードに対して
$2k+d$になっています。( この場合、(3,5,2)に対して$3\times 2+5=11$ )

ということで、除算を計算することなく分岐も扱うことができます。

提出コード

ということを踏まえ、提出コード(Ruby)です。自己再帰でノードの個数を数える、つまり深さ優先探索を行う実装です。CodeIQの実行くんでは、$m=10^{-7}$ ( $k$の上限約$10^7$ ) を、0.7秒程度で処理できます。

終わりに

いやはや、問題では三角関数で条件が書いてありますが、実態はそうではなく。 ( tanを全く使わず計算することも考えましたが、それは可哀想に思えて1回だけ呼ぶことにしました )
じゃあ約数の問題だな、と思ったのですが、実は更に問題を変換すると、木と数列の問題になっているという。移り変わりが面白い問題でした。

実は木を使った説明は解いた後で思いついたもので、当初はノリだけでコードを書いていたのですが。まあ、結果オーライということで。