CodeIQ：プライム・ペア

私自身が表題の問題を解いた時のプログラムについて解説します。
問題の詳細は「プライム・ペア」（CodeIQ）を参照してください。

問題の概要

問題を引用します。
自然数 k に対し、1 から k までの自然数のうち k と互いに素なものの個数を F(k) と定義します。
（F(k) はオイラーのΦ関数とも呼ばれています。参考：オイラーのφ関数（Wikipedia））

例えば F(12)＝4 です。1 から 12 のうち 12 と互いに素なのは 1, 5, 7, 11 の 4 つです。

標準入力から、自然数 n（1 ≦ n ≦ 105）が与えられます。
標準出力に F(n！) を 1000003 で割った余りを出力するプログラムを書いてください。

例えば n＝10 のとき、F(10！)＝F(3628800)＝829440 です。
同様に、F(20！) を 1000003 で割った値は 961998 です。

私のプログラム

Rubyで解答しています。

#!/usr/bin/ruby
require 'prime'

def solve(n)
	primes = Array.new(n+1)
	Prime.each(n){|prime| primes[prime] = prime}

	x = 1
	n.downto(1){|i|
		if primes[i] != nil then next end
		x = x*i%1000003
	}

	pr = primes.compact
	y = x
	for m in pr
		y = y*(m-1)%1000003
	end

	return y
end

# main
while line = gets
	line.strip!
	if line.empty? then next end

	p solve(line.to_i)
end

解説

オイラーのΦ関数については調べれば良いとしてnの最大が1000000で100000!とかいうわけのわからない数字が最大になります。

考え方

まず、オイラーのΦ関数がどんなものかです。
ちゃんとした説明はwikipediaでも見てもらうとして具体的な数字で例を示すと次のようなものです。

12の場合12と互いに素な自然数の個数は、
12=22・3より、12(1-1/2)(1-1/3)=4個
というものです。

つまり、ある自然数nの約数をa1、a2、…、akとすると、
n(1-1/a1)(1-1/a2)…(1-1/ak)
となるというわけです。

ポイントは1点だけです。
nと互いに素な自然数の個数を求めるために必要なのは「nの素因数の種類がわかれば良い」ということです。
これが問題を解く重要なキーになります。

次のポイントは「n以下の数（1以外）はn以下の素数の積の形で表せる」ということです。
具体的にいうとn=12とした場合、
1、2、3、22、5、2・3、7、23、32、2・5、11、22・3
となります。

これは簡単に証明できます。
n以下の数で素数でないものは素因数分解できるので素数の積の形にできます。
そしてその素因数にはnより大きい素数は含まれません（含まれた場合、必ずnより大きくなってしまう）。
よって、n以下の数はn以下の素数の積の形にできます。

ということはn!にはn以下の素数が全て含まれることになります。
これをオイラーのΦ関数と合わせて考えるとF(n!)=n!(1-1/a1)(1-1/a2)…(1-1/ak)で、a1〜akはn以下の全ての素数ということになります。

最後のポイントはn!自体が無茶苦茶大きい数なので自動的に任意精度整数を使えるRubyでもこれを計算するだけで時間が足りなくなる（メモリも足りなくなるかも知れない）ということです。

これについては何度か使ったテクニックですが100003で割った余りだけを処理すれば良いです。

main

入力値を数値にsolve()に渡します。
solve()が結果を返すのでそれを印字して終わります。

solve(n)

答えを計算します。

まず、5〜6行目でn以下の素数を列挙します。
処理ではその数が出現したかどうかだけを知りたいので配列の素数番目にその値をセットしています（[nil,nil,2,3,nil,5,nil,7,…]のようになる）。

9〜12行目の処理はn!を計算するとともにオイラーのΦ関数の割り算部分を同時に処理し、100003で割った余りを求めています。
オーラーのΦ関数を変形すると
F(n!)=n!・(a1-1)/a1・(a2-1)/a2・…・(ak-1)/akになります。
このうちのn!・1/a1・1/a2・…・1/akを計算しているわけです。
割り算の処理は素数を掛けないことで処理しています（10行目）。

14〜18行目で上記で掛けなかった分を掛け算して100003で割った余りを求めます。

結果を返却して終わりです。

A006218(arr)

OEISで調べた式をコード化しただけです。

雑感

かなり難しい問題でした。
私にとって最大のポイントだったのは「n以下の数はn以下の素数の積の形で表せる」ということに気づいた点です。これによって単位n以下の素数を列挙すれば十分であることになり問題を解くことができました。
しかし、100000!ってどんな数になるんでしょう？