CodeIQ:ディビジョン・サム

私自身が表題の問題を解いた時のプログラムについて解説します。
問題の詳細は「ディビジョン・サム」(CodeIQ)を参照してください。

問題の概要

問題を引用します。
自然数 n に対し、(n!)n (つまり n の階乗の n 乗)の約数の和を F(n) と定義します。

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

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

私のプログラム

Rubyで解答しています。

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

def factorial(n)
	ret = {}
	for i in 2..n
		for pr in Prime.prime_division(i)
			if ret[pr[0]] == nil then
				ret[pr[0]] = pr[1]
			else
				ret[pr[0]] += pr[1]
			end
		end
	end
	return ret
end

def power(x, n)
	ret = {}
	x.each{|k,v|
		ret[k] = n * v
	}
	return ret
end

def sunOfGeoMod1000003(r,n)
	return (1-r**(n+1))/(1-r) % 1000003
end

def solve(n)
	ret = 1
	x = factorial(n)
	x = power(x,n)

	x.each{|k,v|
		ret *= sunOfGeoMod1000003(k,v)
	}

	return ret % 1000003
end

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

	p solve(line.to_i)
end

解説

この出題者の★★★の問題としては簡単と思います。というか、高校レベルの数学で解けます。ただ、いくつかの知識が必要になります。

考え方

F(n)=(n!)nでn=1000はわけがわからないほど大きな数になります。Rubyは自動的に長倍精度整数になるので桁あふれの心配はありませんが、こんなものをまともに計算していたらそれだけで時間切れです。

結論から言うと、nを素因数分解してから計算するとうまく行きます。

F(n)の素因数分解を考えます。
xn*xm=xn+mなので、n!を素因数分解したものは、2〜nまでの各数の素因数の乗数を合計したものになります。
例えば、6!=(1)*(2)*(3)*(22)*(5)*(2*3)=24*32*5です。
また、(xn)m =xnmなので、F(n)はn!の各素因数の乗数にnをかけたものになります。
例えば、
F(6)
=6!6
=(24*32*5)6
=224*312*56
です。

約数の和ですがこれには公式があります。確か別の問題でも使いましたが、
x=a0n0*a1n1*…aknk
と素因数分解できる場合、xの約数の合計は、
(1+a0+a02+…+a0n0)*(1+a1+a12+…+a1n1)*…*(1+ak+ak2+…+aknk)
です。

また、(1+ak+ak2+…+aknk)は等比数列なので、1回の計算で合計を求めることができます。初項が1なのでn項までの合計は、
(1-r(n+1))/(1-r)
です。
なので、F(n)を構成する全ての素数について、その乗数までの合計を上記公式で計算し、全て掛け合わせれば良いのですが、とんでもなく大きい数になるので掛け算しているだけで時間切れの可能性があります。

これを回避するため1000003で割ったあまりだけを使って答えを求めます。これも別の問題で同じことをしました。
任意の数xを定数kで割った場合、その商と余りを使ってx=ak+bと表すことができます。
x0=(a0*k+b0)、x1=(a1*k+b1)とした場合、
x0*x1 = a0*a1*k2+(a1*b0+a0*b1)*k+b0*b1
です。
x0*x1をkで割った場合の余はb0*b1%kです(%はkによる剰余)。つまり、x0とx1の余同士を計算すればx0*x1の余を求めることができます。これはx2、x3、…、と掛け合わせる項を増やしても同じです。
この考え方を使うとF(n)を構成する各素数ごとにその乗数までの合計を1000003で割った余を求めて掛け合わせ、その結果をさらに1000003で割った余を求めれば答えになることがわかります。

main

入力値を数値にしてsolve()に渡し、結果を印字します。

solve(n)

問題を解きます。
factorial(n)でn!の素因数分解の結果を得ます。
power(x,n)でn!の素因数分解の結果から(n!)nの素因数分解の結果を得ます。
結果は{素数=>乗数}の形式の連想配列なので、全てのキー(素数)と値(乗数)に対して、等比数列の和を1000003で割った結果を求め、掛け合わせます。
その答えをさらに1000003で割ったあまりを返します。

factorial(n)

n!の素因数分解の結果を求めます。
2〜nまでの全ての数を素因数分解し、{素数=>乗数}の形でretに記録します。
新しい素数が現れたらキーを値を追加し、すでに存在する素数なら値を合計します。

power(n)

{素数=>乗数}の形の連想配列xの全ての要素の値に対してnを掛けた結果を返します。

sunOfGeoMod1000003(r,n)

1〜rnまでの等比数列の和を100003で割った余を返します。
計算は考え方に示した公式を使って計算します。

雑感

部分的に以前やったことのある問題の知識が利用できたのですんなりとできましたが、約数の合計と余を使った計算の経験がなかったら苦労したかもしれません。
私の回答はRubyなのでうまくいっていますが、多分途中で64bit整数の範囲を超える数になっていると思います。なので、C言語はもちろん、C++とかJavaのような静的型付けの言語ではもっと細かく1000003で割って余だけを使った計算にする必要があります。が、うまくできるのかな? sunOfGeoMod1000003(r,n)の1-r**nで64bit整数の範囲を越えてそうなんだけど……。