読者です 読者をやめる 読者になる 読者になる

PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 401 / 約数の平方和

6の約数は 1,2,3,6 であり,これらの数の平方和は 1+4+9+36=50 となる。

n の約数の平方和を sigma2(n) で表し,sigma2 の総和を SIGMA2 で表す。
SIGMA2(n)=Σsigma2(i) (i=1 から n まで)

SIGMA2 の最初の6項は 1,6,16,37,63,113 である。

SIGMA2(10^15) modulo 10^9 を求めよ。

Problem 401 - Project Euler


計算量を明確に意識させられる問題でした。問題文どおりに関数を定義すると当然解けない。10^8 までの計算で約2分かかります。

最初に考えたのは「約数としての登場回数を数える」方法です。1~n を考えたとき,q が約数として登場するのは [n/q] 回([ ] はガウス記号)。SIGMA2は次のように書けます。

 S_1(n)=\displaystyle\sum_{q=1}^n q^2 \left[\dfrac{n}{q}\right]

これだと n=10^8 で40秒。n=10^15 を1分におさめるにはまだ遠い。

そこで考えたのは「ガウス記号の値で場合分け」です。q が大きくなると,ガウス記号の値はほとんど変化しなくなります。ガウス記号 =k とおくと

 k \leqq \dfrac{n}{q}< k+1 \quad \therefore \left[\dfrac{n}{k+1}\right]+1\leqq q\leqq \left[\dfrac{n}{k}\right]

これを使ってSIGMA2をあらわします。

 S_2(n)=\displaystyle\sum_{k=1}^n \sum_{q=\left[\frac{n}{k+1}\right]+1}^{\left[\frac{n}{k}\right]} kq^2

これも単独では遅すぎるので,小さな n には S1(n)を使い,大きな n には S2(n) を使いました。

問題はどこで切り替えるかです。はじめは 10^m の形の数で切っていたのですが,[n/√n]=[√n] で切るとよいようでした。こうすることで S1(n), S2(n) 両方の計算量を O(√n) におさえられます。慣れていれば一発で気づきそうな結論ですが,私はかなり時間がかかりました。