PowerShell で Project Euler 10

問題

• 本家：Problem 10 - Project Euler
• 翻訳：Problem 10 - PukiWiki

10以下の素数の和は 2 + 3 + 5 + 7 = 17 である.
200万以下の全ての素数の和を求めよ.

エラトステネスの篩を使って

2000000 |
    ? {$_ -gt 0} |
    % {
        $sum = 0

        $skip = New-Object System.Collections.BitArray ($_ - 1)
        foreach($p in 2..$_)
        {
            if ($skip.Get($p - 2))
                {continue}

            for($i = $p * $p; $i -le $_; $i += $p)
                {$skip.Set($i - 2, $true)}

            $sum += $p
        }

        $sum
      }

倍数和の再帰関数を使って

ワケガワカラナイヨ

2000000 |
    ? {$_ -gt 0} |
    % {
        # √N までの素数を求める
        $sqrtN = [Math]::Floor([Math]::Sqrt($_))
        $skip = New-Object System.Collections.BitArray ($sqrtN - 1)
        $sqrtP = New-Object System.Collections.Generic.List[Int]

        foreach($p in 2..$sqrtN)
        {
            if ($skip.Get($p - 2))
                {continue}

            for($i = $p * $p; $i -le $sqrtN; $i += $p)
                {$skip.Set($i - 2, $true)}

            $sqrtP.Add($p)
        }

        # メモ化した再帰関数 M を定義
        #   M(x, y, n) = 0～n番の素数の倍数のうち、
        #   区間 [x, y] に含まれる数の総和
        $memo = @{}
        filter M ($x, $y, $n)
        {
            $sum = 0

            if ($x -gt $y)
                {return $sum}

            for (; $n -ge 0 -and $y -lt $sqrtP[$n] * $sqrtP[$n]; $n--)
                {if ($x -le $pn -and $pn -le $y)
                    {$sum += $pn}}

            if ($n -lt 0)
                {return $sum}

            $key = $x, $y, $n -join ","

            if (-not $memo.Contains($key))
            {
                foreach ($k in $n..0)
                {
                    $pk = $sqrtP[$k]
                    $i = [Math]::Max([Math]::Ceiling($x / $pk), $pk)
                    $j = [Math]::Floor($y / $pk)

                    $sum += $pk * ($j - $i + 1) * ($j + $i) / 2
                    $sum -= $pk * (M $i $j ($k - 1))

                    if ($x -le $pk)
                        {$sum += $pk}
                }

                $memo[$key] = @(0, $sum)
            }

            $memo[$key][0]++
            return $memo[$key][1]
        }

        #   「2～Nの総和」
        # - 「2～Nの√N以下の素数の倍数の総和」
        # + 「√N以下の素数の総和」
          ($_ * ($_ + 1) / 2 - 1) `
        - (M 0 $_ ($sqrtP.Count - 1)) `
        + ($sqrtP | % {$s = 0} {$s += $_} {$s})
      }

再帰関数 M の導出

準備として、関数$M'$を定義する。$M'$は自然数の集合$A \subset \mathbb{N} $の任意の要素$a \in A$の倍数であって$x$以上$y$以下の値の総和である。ただし、$x, y \in \mathbb{R}_{\ge 0}$ は非負実数とする。

$$M'(x, y, A) := \sum \left([x, y] \cap \bigcup_{a \in A} a \mathbb{N} \right)$$

ややこしいが、M' は倍数の和を表す関数である。
M' については次の性質がある。

• $x > y \Rightarrow M'(x, y, A) = 0$
• $M'(x, y, \emptyset) = 0$
• $y < \min A \Rightarrow M'(x, y, A) = 0$
• $M'(x, y, A) $$ = M'(\lceil x \rceil, y, A)$$ = M'(x, \lfloor y \rfloor, A)$
• $x \le y \Rightarrow M'(x, y, \{ n \}) = \sum_{k = \lceil \frac{x}{n} \rceil}^{\lfloor \frac{y}{n} \rfloor} n k = \frac{n}{2} (\lfloor \frac{y}{n} \rfloor - \lceil \frac{x}{n} \rceil + 1)(\lfloor \frac{y}{n} \rfloor + \lceil \frac{x}{n} \rceil)$

実際にM'を計算する場合は次の式が使える。

$M'(x, y, \{n\} \cup A) $
$ = M'(x, y, \{n\}) $$+ M'(x, y, A) $$- M'(x, y, \{ \text{lcm}(a, n) | a \in A \})$
$ = M'(x, y, \{n\}) $$+ M'(x, y, A) $$- n M' \left(\frac{x}{n}, \frac{y}{n}, \big\{ \frac{a}{\gcd(a, n)} \big| a \in A \big\} \right)$

ただし、集合 $A = \{a_0, a_1, \dots, a_n \}$ の任意の異なる２つの要素$a_i, a_j \in A ~ (i \ne j)$が互いに素であるとき、次のように展開できる。（当然ながら異なる素数は互いに素である。）

$M'(x, y, A) $
$ = M'(x, y, \{a_n\}) $$+ M'(x, y, A \setminus \{a_n\}) $$- a_n M' \left(\frac{x}{n}, \frac{y}{n}, A \setminus \{a_n\} \} \right) $
$ = \sum_{k = 0}^n M'(x, y, \{a_{k}\}) $$ - \sum_{k = 1}^n a_k M'(\frac{x}{a_k}, \frac{y}{a_k}, A \setminus \{a_n, a_{n - 1}, \dots, a_{k} \}) $

ここで関数 M を M' によって次の通り定義する。

$$M(x, y, n) := M'(x, y, \{p_0, p_1, \dots, p_n\}) $$

ただし、$p_k \in \mathbb{P}$は素数である。添字 $n$ は$0$以上の整数であって、小さい方から素数を並べたときの順番（0 始まり）を表している。最初の３つを書けば、$(p_0, p_1, p_2) = (2, 3, 5)$である。

Mについては次の性質がある。

• $p_n \notin [x, y] \land y \lt p_n^2 \Rightarrow M(x, y, n) = M(x, y, n - 1)$
• $p_n \in [x, y] \land y \lt p_n^2 \Rightarrow M(x, y, n) = M(x, y, n - 1) + p_n$

さて、$N$以下の素数の個数を$\pi(N)$と書くことにする。
$p_{\pi(N) - 1} \le N \lt p_{\pi(N)}$であることに注意すると、$2, 3, \dots, N$のすべての値は、$p_0, p_1, \dots, p_{\pi(N) - 1}$ の倍数であるから、次が成り立つ。

$$M(2, N, \pi(N) - 1) = \sum_{k = 2}^{N} k = \frac{1}{2}N(N + 1) - 1$$

ところで、

$$2 \le p_k \le N \land N < p_k^2 ~ (k = \pi(\sqrt{N}), \pi(\sqrt{N}) + 1 \dots \pi(N) - 1)$$

であるから、

$$M(2, N, \pi(N) - 1) = M(2, N, \pi(\sqrt{N}) - 1) + \sum_{k = \pi(\sqrt{N})}^{\pi(N) - 1} p_k$$

すなわち、

M(2, N, \pi(\sqrt{N}) - 1) + \sum_{k = \pi(\sqrt{N})}^{\pi(N) - 1} p_k = \frac{1}{2}N(N + 1) - 1\\
\sum_{k = \pi(\sqrt{N})}^{\pi(N) - 1} p_k = \frac{1}{2}N(N + 1) - 1 - M(2, N, \pi(\sqrt{N}) - 1) \\
\sum_{k = 0}^{\pi(N) - 1} p_k = \frac{1}{2}N(N + 1) - 1 - M(2, N, \pi(\sqrt{N}) - 1) + \sum_{k = 0}^{\pi(\sqrt{N}) - 1} p_k

となるので、$\sqrt{N}$以下の素数を求めて、上記の計算をすれば素数和が求まる。

以上

おわり。