七夕問題☆牽牛 彦星 牛をもっと飼う【解説2】
本エントリはCodeIQで2014年7月7日～同年7月28日に出題された「七夕問題☆牽牛 彦星 牛をもっと飼う」という問題とその解法に関する記事です。
七夕問題☆牽牛 彦星 牛をもっと飼う【解説2】
◆互いに素な値を簡単に数える方法（その2）
実はもう一つ、互いに素な値を数える方法があります。先程と同じ例の、y座標が6の場合、(1, 6), (2, 6), (3, 6), (4, 6), (5, 6), (6, 6)の6つの点に関して考えます。1～6の中で、互いに2で割り切れるものは、(2, 6), (4, 6), (6, 6)の3種類あります。これは候補から除くことができます。
同様に、互いに3で割り切れるものは、(3, 6), (6, 6)の2種類あります。これらのうち、(6, 6)はいずれにも含まれています。共通した部分は、そのまま候補から除こうとすると、二重に取り除くことになってしまいますので、
6 - 3 - 2 + 1
6は全体の数、3は2の倍数の数、2は3の倍数の数、1は6(2×3)の倍数の数です。小学生でならった、「2でも3でも割り切れる数はいくつでしょう」みたいな問題と同じ考え方です。
これを、y座標が60の場合にも適用してみましょう。60を割り切れる素数は、2, 3, 5です。3つしかありませんが、計算方法はかなり複雑になります。
2と3の倍数(6の倍数)、2と5の倍数(10の倍数)、3と5の倍数(15の倍数)、2と3と5の倍数（30の倍数）についても計算しておく必要があります。
60÷2 = 30 60÷3 = 20 60÷5 = 12 60÷6 = 10 60÷10 = 6 60÷15 = 4 60÷30 = 2
を用いて、
60 - 30 - 20 - 12 + 10 + 6 + 4 - 2 = 16
2つの倍数の共通部分は2重に引いてしまっていますから、足してやらなければなりませんが、そうすると3つの数の倍数と共通部分を余分に足してしまうことになってしまいますので、その分は引いておきます。複雑そうではありますが、何とか数えることはできました。
◆メビウス関数
この計算をもう少し一般化したものは、メビウス関数と呼ばれています。
ここからは少し数学的要素が強くなります。細かい証明や難しい用語はできるだけ省いていきますので、数学が苦手な方も頑張って読んでくださいね。（乗算の記号も×から*に変えていきますのでご了承ください）
メビウス関数は、μ（ミュー）で表し、μ(n)について、
【1】nが1の場合は1を返す
【2】nを素因数分解したとき、素因数の指数部分がすべて1で、なおかつ素因数の数が、
奇数個の場合 … -1 偶数個の場合 … +1
を返します。いくつか具体的例を挙げておきます。
μ(6)の場合
6 = 2^1 * 3^1
なので、+1
μ(30)の場合
30 = 2^1 * 3^1 * 5^1
なので、-1
μ(60)の場合
60 = 2^2 * 3^1 * 5^1
で、22が含まれているので0
【3】それ以外の場合は0を返す たとえば、μ(4)は、22で指数部分が1でない素因数があります。その場合は0とします。
メビウス関数を利用すると、φ関数は和で表すことができます。 たとえば、φ(6)の場合、
φ(6) = μ(1)*6 + μ(2)*3 + μ(3)*2 + μ(6)*1 = 6 - 3 - 2 + 1 = 2
φ(60) = μ(1)*60 + μ(2)*30 + μ(3)*20 + μ(4)*15 + μ(5)*12 + μ(6)*10 + μ(10)*6 + μ(12)*5 + μ(15)*4 + μ(20)*3 + μ(30)*2 + μ(60)*1 = 60 - 30 - 20 + 0 - 12 + 10 + 6 + 0 + 4 + 0 - 2 + 0 = 16
μ(d)に掛けている数字は何かというと、その列の要素数です。つまり、 μ(n)は、nの約数をdとすると、
μ(n) = μ(d1)*n/d1 + μ(d2)*n/d2 + μ(d3)*n/d3 + ...
のような計算をすることになります。
◆さらに高速に解を得る
ここまで読んだ限り、結局素数列は必要そうに見えますし、足し算とは言っても、いろいろと計算する必要があり、あまり高速化しそうに見えないかもしれません。しかし思い返してください。必要なのは、φ(k)が必要なのではなく、φ(1) + φ(2) + φ(3) + … + φ(7777776)が欲しいのです。
そこで、φ(1)から順に書き並べてみます。
φ(1) = μ(1)*1 φ(2) = μ(1)*2 + μ(2)*1 φ(3) = μ(1)*3 + μ(3)*1 φ(4) = μ(1)*4 + μ(2)*2 + μ(4)*1 φ(5) = μ(1)*5 + μ(5)*1 φ(6) = μ(1)*6 + μ(2)*3 + μ(3)*2 + μ(6)*1 φ(7) = μ(1)*7 + μ(7)*1 φ(8) = μ(1)*8 + μ(2)*4 + μ(4)*2 + μ(8)*1 φ(9) = μ(1)*9 + μ(3)*3 + μ(9)*1 φ(10) = μ(1)*10 + μ(2)*5 + μ(5)*2 + μ(10)*1 φ(11) = μ(1)*11 + μ(11)*1 φ(12) = μ(1)*12 + μ(2)*6 + μ(3)*4 + μ(4)*3 + μ(6)*2 + μ(12)*1 φ(13) = μ(1)*13 + μ(13)*1 . . .
パッと見た感じですぐに気付くのは、必ずμ(1)が含まれていて、掛ける値が1ずつ増えているという点だと思います。μ(2)についても、よく見ると1つおきに現れていて、掛ける数が1, 2, 3, …となっています。これをたくさん書いていくと、μ(3)以降も同じようになっていることが確認できるはずです。よく考えてみると、当たり前のことですね。2の倍数は2を約数に持つのだから、μ(2)が出現する、3の倍数は3を約数に持つのだから、μ(3)が出現する。 したがって、φ関数の和（= f(n)とします）は、
f(n) = μ(1)*(1+2+...+n) + μ(2)*(1+2+...n/2) + μ(3)*(1+2+...+n/3) + ...
のようになります。1+2+…+kは等差数列の和になりますから、k*(k+1)/2で計算できますので、意外とスッキリさせることができます。
ここでg(k) = k*(k+1)/2とすると、
f(n) = μ(1)*g(n/1) + μ(2)*g(n/2) + μ(3)*g(n/3) + ...
となります。もう少しシンプルに書くと、
f(n) = Σμ(d)g(n/d) {d=1..n}
みたいなかんじです。 n/dのように書くと、「割り切れない場合があるじゃん！」と思い方もいらっしゃるかもしれません。確かに割り切れない場合、困りますね。n=3とすると、
f(3) = μ(1)g(3/1=3) + μ(2)g(3/2=???) + μ(3)g(3/3=1)
となります。3/2はどう扱えばよいのかというと、φ(1) + φ(2) + φ(3)と見比べてみればすぐにわかるはずです。μ(2)が出現するのは、μ(2)*1の1か所だけです。2飛びでμ(2)が出現するわけですから、単に3/2で、割り切れない部分は切り捨ててしまえば、それでOKです。 というわけで、小数点以下は切り捨てる場合は、[n/d]のように表記します。
つまり先程の式は、
f(n) = Σμ(d)g([n/d]) {d=1..n}
と書くことにします。
メビウスの反転公式
ここで再び、具体的な値で見てみましょう。
g(1) = 1*2/2 = 1 g(2) = 2*3/2 = 3 g(3) = 3*4/2 = 6 g(4) = 4*5/2 = 10 g(5) = 5*6/2 = 15 g(6) = 6*7/2 = 21 f(1) = μ(1)*g(1) = 1 f(2) = μ(1)*g(2) + μ(2)*g(1) = 3 - 1 = 2 f(3) = μ(1)*g(3) + μ(2)*g(3/2=1) + μ(3)*g(1) = 6 - 1 - 1 = 4 f(4) = μ(1)*g(4) + μ(2)*g(2) + μ(3)*g(4/3=1) + μ(4)*g(1) = 10 - 3 - 1 + 0 = 6 f(5) = μ(1)*g(5) + μ(2)*g(5/2=2) + μ(3)*g(5/3=1) + μ(4)*g(5/4=1) + μ(5)*g(1) = 15 - 3 - 1 + 0 - 1 = 10 f(6) = μ(1)*g(6) + μ(2)*g(3) + μ(3)*g(2) + μ(4)*g(6/4=1) + μ(5)*g(6/5=1) + μ(6)*g(1) = 21 - 6 - 3 + 0 - 1 + 1 = 12
この、fとgの値をよく見ると、面白い関係性に気付きます。
g(1) = f(1) g(2) = f(2) + f(1) g(3) = f(3) + f([3/2]=1) + f(1) g(4) = f(4) + f([4/2]=2) + f([4/3]=1) + f(1) g(5) = f(5) + f([5/2]=2) + f([5/3]=1) + f([5/4]=1) + f(1) g(6) = f(6) + f([6/2]=3) + f([6/3]=2) + f([6/4]=1) + f([6/5]=1) + f(1)
これは、nの値をどんどん大きくしても、常に成り立ち、
g(n) = Σf(d) {d=1..n}
と表すことができます。この等式と、先で説明した
f(n) = Σμ(d)g(n/d) {d=1..n}
と見比べると、fとgが入れ替わったような形をしています。これらの式は、メビウスの反転公式と呼ばれています。 さらに具体的にみていきましょう。
g(2)の式を変形すると、
f(2) = g(2) - f(1) = g(2) - g(1)
g(3)の式を変形すると、
f(3) = g(3) - f(1) - f(1) = g(3) - g(1) - g(1)
g(4)の式を変形すると、
f(4) = g(4) - f(2) - f(1) - f(1)
ここで、f(2)はすでに計算済みですから、
= g(4) - { g(2) - g(1) } - g(1) - g(1)
以上のように、fはすべてgの式に変換することができます。gは簡単に計算することができますから、非常にシンプルな計算になります。
n = ARGV[0].to_i $cache = Hash.new def f(n) return 1 if n==1 s = n*(n+1)/2 (2..n).each{|k| s -= $cache[n/k] = $cache[n/k] || f(n/k) } s end puts( f(n-1) * 2 + 1 )
このコードを実際に動かしてみると、消費メモリは格段に小さくなりましたが、最初のコードと比べると少し時間がかかっています。一度計算したfは、キャッシュに保持していますが、足し算の回数が多すぎて、それなりに時間がかかってしまっているようです。
そこで、最後の仕上げとして、少し計算の工夫をしておきます。
f(60)を呼び出すと、
f([60/2]=30), f([60/3]=20), f([60/4]=15), ...
のようにfが再帰的に呼び出されます。fの引数部分だけを書きならべてみると、
30, 20, 15, 12, 10, 8, 7, 6, 6, 5, 5, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
のように、後半は同じ値で、しかも小さな値が大量に続きます。 この部分は、工夫すれば簡単にまとめられそうです。
72 < 60 < 82 ですので、[60/d]のdを、2から7まで【A】と、それより後の部分【B】で計算方法を変えてやります。 【A】についてはそのまま計算するとして、【B】の部分は、次のように考えることができます。
60/60 = 1 60/30 = 2 60/20 = 3 60/15 = 4 60/12 = 5 60/10 = 6 [60/8] = 7
から、dの値が
31～60の場合はすべてf(1) 21～30の場合はすべてf(2) 16～20の場合はすべてf(3) 13～15の場合はすべてf(4) 11～12の場合はいずれもf(5) 9～10の場合はいずれもf(6) 8の場合はf(7)
であることがわかります。つまり、【B】の部分は
f(1)*30 + f(2)*10 + f(3)*5 + f(4)*3 + f(5)*2 + f(6)*2 + f(7)*1
と、大幅に計算量を減らすことができます。これを実装すると、次のようになります。
n = ARGV[0].to_i $cache = Hash.new def f(n) return 1 if n==1 sum = n*(n+1)/2 sq = Math.sqrt(n).to_i # 前半部分は普通に計算 (2..n/sq).each{|k| sum -= $cache[n/k] = $cache[n/k] || f(n/k) } # 後半部分は同じ値になる部分をまとめて計算 (2..sq).each{|k| sum -= ( $cache[k-1] = $cache[k-1] || f(k-1) ) * ( n/(k-1) - n/k ) } return sum end puts( f(n-1) * 2 + 1 )
これで、7777777のサイズでも1秒以内に解を得ることができ、メモリ消費もほとんどありません。
