CodeIQ MAGAZINE

問題のおさらい

自然数 n に対し、(n！)ⁿ　（つまり n の階乗の n 乗）の約数の和を F(n) と定義します。

例えば、F(2) ＝ 7 です。(2！)² ＝ 4 で、4 の約数は 1, 2, 4 だからです。
また、F(3) ＝ 600 です。(3！)³ ＝ 216 で、216 には 16 個の約数があり、それらの和は 600 です。
同様に、F(4) ＝ 991111 となります。

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

方針：素因数分解された形で考えよう

(n！)ⁿ という非常に大きな数の約数に関する問題です。(n！)ⁿ の値を実際に計算してみるという方法は非効率だと多くの方が気づかれたことと思います。

ではどうすればよいでしょうか。(n！)ⁿ を 素因数分解された形で考える ことがポイントです。下記の Wikipedia の記事が参考になります。一般に、整数 k が 、互いに異なる素数 p₁, p₂, …, p_m と、その指数 a₁, a₂, …, a_m を使って、p₁^a₁×p₂^a₂× … ×p_m^a_m と表されるとき、k の約数の和は次のように表されます。

参考：約数：約数の和（Wikipedia）

したがって本問では、(n！)ⁿ を素因数分解したとき、上記の p₁, p₂, …, p_m と a₁, a₂, …, a_m がどのような値になるかを考えましょう。もちろん、(n！)ⁿ の値を直接計算するのでなく、(n！)ⁿ という形に着目して考えましょう。

(n！)ⁿ という数は、整数の積として見ると、1 から n までの整数が n 個ずつ集まったものの積です。そこでまずは 1 から n までの整数を、一つずつ素因数分解しましょう。各素数が何回現れるかをそれぞれ求め、n までの総数を計算します。そしてその数を n 倍したものが、(n！)ⁿ の素因数の指数と分かります。

以上の方法で、(n！)ⁿ を素因数分解したときの p₁, p₂, …, p_m と a₁, a₂, …, a_m の値が求められます。あとは、上の約数の和の公式により、答えを計算できます。

コード例は以下です（Python）。前半では、上記の a₁, a₂, …, a_m に対応する配列 a を用意し、素数 p に対する指数を a[p] に保存します。後半では、公式の計算を行っています。公式の各カッコの中は、公比 p の等比数列です。1 から始めて、次々と p をかけながら和をとっていきましょう。以下は Python のコード例です。

n = input()
a = [0] * (n+1)
D = 1000003
for k in range(2, n+1):
    p = 2
    while k != 1:
        while k % p == 0:
            a[p] += n
            k /= p
        p += 1

answer = 1
for p in range(2, n+1):
    t, s = 0, 1
    for j in range(a[p]+1):
        t = (t + s) % D
        s = (s * p) % D
    answer = (answer * t) % D
print answer

このコードで、n ≦ 1000 という条件を 1 秒の制限時間でクリアできます。

さらなる高速化！ 等比数列の和の公式を使おう

上のコードでひとまずクリアは可能ですが、実は本問にはさらなる高速化のポイントがあります。注目すべきは後半です。ここでは 1＋p＋p²＋…＋p^a の値を求めるために計 a 回の乗算と加算を行っています。n ≦ 1000 という条件だと、特に小さい p で a はけっこう大きな値になります（p＝2 なら a は最大で約 10⁶。）

そこで役に立つのが、等比数列の和の公式です。

この考えは過去の記事でも登場しました。（参考：数学の問題をプログラミングで解こう！「エース・ナンバー」問題解説の方法２。）「p の a＋1 乗の計算」と「p－1 での割り算」を行うことがポイントです。それぞれ「べき剰余」と「100003 を法とする p－1 の逆元」を求めることで、上よりもずっと短い計算時間で求められます。

べき剰余と逆元については、下記のサイトを参考にして下さい。

参考：冪剰余（Wikipedia）

参考：整数の合同と１次合同式（外部サイト）

さらに前半で、1 から n の値を素因数分解する部分でも、高速化が可能です。上のコードでは、ある整数 k に対し、割る数を 1 ずつ変えながら割り算を行っていました。ここでは、素数のそれぞれが 1 から n の積に素因数として何個入っているかを調べましょう。たとえば n が 14 の場合、素数 2 が何個入っているかは、図のようにして求めることができます。

始めにエラトステネスの篩で素数の一覧を作り、素数のそれぞれについて上記の計算を行いましょう。それを n 倍したものが、(n！)ⁿ の各素因数の指数となります。

コード例を以下に示します。はじめのコードでは、n＝1000 のとき 0.5 秒程度の実行時間で答えを出せていましたが、今回のコードではそれよりもずっと短い実行時間（実行くんで 0.01 秒程度）で答えを出せます。

# 拡張ユークリッドの互除法
# a*x + b*y = gcd(a,b) となる (x,y) を求める.
def gcd_ext(a, b):
    x, y, lastx, lasty = 0, 1, 1, 0
    while b != 0:
        q = a/b
        a, b = b, a%b
        x, y, lastx, lasty = lastx-q*x, lasty-q*y, x, y
    return (lastx, lasty)

# d を法とする a の逆元 (a*x = 1 (mod d) となる x) を求める.
def inv_mod(a, d):
    (p, q) = gcd_ext(a, d)
    return p % d

# a^n mod d を求める.
def pow_mod(a, n, d):
    if n == 0: return 1
    if n % 2 == 0: return pow_mod((a*a)%d, n/2, d)
    return a * pow_mod(a, n-1, d) % d

n = input()
# エラトステネスの篩を作る
prime = [True] * (n+1)
for i in range(2, n+1):
    for j in range(2*i, n+1, i):
        prime[j] = False

a = [0] * (n+1)
D = 1000003
for p in range(2, n+1):
    if prime[p] == False: continue
    k = p
    while k <= n:
        a[p] += n*(n/k)
        k *= p

answer = 1
for i in range(2, n+1):
    t = (pow_mod(i, a[i]+1, D) - 1) * inv_mod(i-1, D)
    answer = (answer * t) % D
print answer

みんなのコード

本問の挑戦者で、コードを公開して頂いた方のコードを Togetter でまとめました。（随時追加していきます。）こちらもぜひチェックして下さい。

CodeIQ「ディビジョン・サム」 みんなのコード

CodeIQ運営事務局より

kawazoeさん、ありがとうございました！
kawazoeさんの次回の問題にご期待下さい。