サイエンスものおきば

　割り算と桁数のパズルの解答編です。

　真っ先に思いつくのは、A, B の組を一つ一つ選んで桁数を調べるという O(N²) のやり方です。しかし、A, B の組は全部で 10¹² × 10¹² ＝ 10²⁴ 通りあるので、これを一つ一つ割り算してたしかめていくというのは到底不可能です。そこで、「A, B が一体どういう条件を満たしていれば小数点以下が D 桁になるか？」という条件を考えていきます。

■ O(N²) から O(N) へのアプローチ

　ポイントは、分数 A／B を約分したときの分母の形に注目することです。「約分後の分母に 2 または 5 以外の因数が含まれないこと」がまず必須の条件であることが分かります。もしそのような因数があれば A÷B は循環小数になってしまうからです。さらに、D の値は、「2 と 5 のべき指数のうち大きいほうの値」に一致することも分かります。すなわち、A÷B の小数点以下がちょうど D 桁となるための必要十分条件とは、「約分後の分母が 2^p × 5^q の形で書け、かつ max(p, q) = D」であるといえます。

　つまり、条件を満たす A と B は、約分後の分子を m、最大公約数を g とおいて、それぞれ次のように表せます。
　
　（ただし m と2^p × 5^q は互いに素で、かつ max(p, q) = D）

　結局のところ、本問は、上の条件を満たす４整数の組 (g, m, p, q) の個数を求めることと同値であるといえます。

　そのような組を数えるアルゴリズムを考えます。始めに p と q の値を固定し、次に m の値を固定するという順番で考えます。ある p, q を選んだとき、m はどのような値でしょうか。３通りに場合分けして考えます。① p = D, q = 0 の場合、m は 2 で割り切れない任意の数です。② p = 0, q = D の場合、m は 5 で割り切れない任意の数です。③ それ以外（p ≧ 1, q = D または p = D, q ≧ 1）の場合、m は 2 でも 5 でも割り切れない任意の数です。m を決めた後の g の選び方は単純です。ベクトル (m, 2^p × 5^q) を整数 g 倍したすべての組は上の条件を満たします。ですのでそのような g の個数は、N を m と 2^p × 5^q の大きい方で割れば求められます。

　公式化すると次のようになります。割り算は端数を落としていることに注意して下さい。これを単純に実装すれば、O(N) のプログラムができます。
　
　

■ O(N) から O(√N) へのアプローチ

　さて、本問は N = 10¹² です。O(N) ではまだ時間がかかりすぎます。さらに計算時間を下げる必要があります。

　いちばん時間がかかっているのが、m に関する繰り返し計算です。ここを高速化します。max の大小で式を分解すると、前半部分は、シグマの中身が m の値に依存しないので、単純な掛け算で計算可能です。ただ m の値は 2 または 5 で割り切れない制約があるため、そうした数を除外します。

　後半部分は、単純化すると、m に対する ∑ ((定数) / m) を計算する問題です。これは一つ一つ計算すると時間がかかりますが、ある程度 m が大きいと (定数) / m は同じ値が続くという事実に着目すると、一気にオーダを落とせます。つまり、ある整数 k に対し、(定数) / m = k となるような m は何個存在するかを考えるということです。そのような m は (定数)/k - (定数)/(k+1) 個存在します。つまり、始めは m をだんだんと大きくしながら m に対する ∑ ((定数) / m) を繰り返し計算し、ある程度 m が大きくなったところで、ループ変数を m から k に切り替えて、∑ k × ((定数)/k - (定数)/(k+1)) を計算します。なお前半部分と同様、2 または 5 で割り切れない場合を除外しましょう（詳細は省略します）。
　
　m と k を切り替えるタイミングは m が √N に等しくなったときがベストです。これによりアルゴリズムは O(√N) になります。

　以上のコードをこちらに記します。上の説明では省いた細かなケースのバグを除去しています。

　実行させると、答えの 83566672035424 が得られます。実行時間は上のコードで１秒を切るはずです。

　なお、本問の正解者は、順に有為さん、triceps さん、stqn さんでした。Project Euler の Fastest Solvers によく登場する方々ですね。ありがとうございます！