素因数の個数の平均

Problem 468の計算量を見積もるときにこんな問題が思い浮かびました。

N以下を素因数分解したとき、素因数の個数の平均は?

例えば、2、3、4 = 22、5なら1個、6 = 2 * 3なら2個とします。logNよりはずっと小さいと思われます。とりあえずコードを書いて調べてみましょう。

from itertools import *

def div_pow(n, d):
    e = 0
    while n % d == 0:
        e += 1
        n /= d
    return e, n

def make_factors_table(N):
    a = range(N + 1)
    for p in takewhile(lambda p: p * p <= N, (n for n in count(2) if a[n] == n)):
        for m in xrange(p, N + 1, p):
            if a[m] == m:
                a[m] = p
    
    b = [ () ] * (N + 1)
    for n in xrange(2, N + 1):
        p = a[n]
        e, m = div_pow(n, p)
        b[n] = ((p, e),) + b[m]
    return b

print "|*N|*mean|"
for e in xrange(2, 8):
    N = 10 ** e
    fss = make_factors_table(N)
    s = sum(1 for fs in islice(fss, 1, None) for p, e in fs)
    print "|%d|%.3f|" % (N, float(s) / N)
N mean
100 1.710
1000 2.126
10000 2.430
100000 2.664
1000000 2.854
10000000 3.013

loglogNっぽいですよね。この問題を思いついてlogLogNの気がした瞬間、思いつきました。整数が2を素因数に持つ確率は1/2、3を持つ確率は1/3になります。つまり、Nまでの素因数の個数は、

 \sum_{p \le N}{\frac{N}{p}}

です。これを積分で近似します。素数定理よりp付近の素数の密度は1/logpだから、

 \int_{1.5}^N{\frac{Ndp}{p\log{p}}

この積分はエラトステネスのふるいの計算量を見積もるときに出てくるものですね。loglogp微分してみればわかります。この定積分はだいたい、N(loglogN + c)で平均はloglogN + cになります。
先ほどの表にloglogNを付け加えてみましょう。

N mean loglogN
100 1.710 1.527
1000 2.126 1.933
10000 2.430 2.220
100000 2.664 2.443
1000000 2.854 2.626
10000000 3.013 2.780

loglogN + 0.2くらいのようですね。
あと、22なら2個という数え方もあります。こうすると個数は増えます。しかし、2でも

1/2 + 1/4 * 2 + 1/8 * 3 + ... = 2

なので、倍には増えないことになります。

N mean loglogN
100 2.390 1.527
1000 2.877 1.933
10000 3.199 2.220
100000 3.436 2.443
1000000 3.627 2.626
10000000 3.786 2.780

1.1loglogN + 0.7くらいのようです。