前回（その２）では平方因子を$r$個だけ持つ数を数えるためにメビウス関数を使った式を導きました。
今回はまったく違った考え方でNを大きくしたときの平方因子を$r$個だけ持つ数の割合$RSF(r,N)$を求めます。
確率を使って無平方数の割合RSF(0,N)を求める
まずnが１から10000までの正の整数とし、この中からランダムに１つを選んだときに平方数$2^2=4$の倍数でない確率は$1-1/2^2=3/4$です。
無平方数は$p^2 \le 10000$のすべての素数の$p$の２乗の倍数ではないので、$RSF(0,N)$はすべての素数について$(1-1/p^2)$の積をとればよいことになります。
\begin{align}
& RSF(0, N) = \prod_{p_i}^{p_i \le \sqrt{N}}(1-\frac{1}{p_i^2}) \cdots (1)
\end{align}
同様にRSF(1,N)を求める
まず例えば$2^2$だけの倍数でほかの素数の倍数でない確率は、
\begin{align}
& \frac{1}{2^2}\prod_{p_i \ne 2}^{p_i \le \sqrt{N}}(1-\frac{1}{p_i^2})
\end{align}
これをすべての素数が各々１つだけ含まれる確率を足し合わせれば良いので
\begin{align}
& RSF(1, N) = \sum_{p_i}^{p_i \le \sqrt{N}} \frac{1}{p_i^2}\prod_{p_j \ne p_j}^{p_j \le \sqrt{N}}(1-\frac{1}{p_j^2})
\end{align}
$r \ge 2 $ の場合も同様に$p_i \le \sqrt{N}$から$p_i$を2個取り出せばよいのですが、式が複雑になるので省略します。
しかし確率母関数を使うとこれが簡単に表現できます。
二項分布の確率母関数の応用
確率qが一定のとき、二項分布の確率母関数は
\begin{align}
G_x(S) &= \sum_{x=0}^{n}s^x\binom{n}{x}q^x(1-q)^{n-x}\\
&= (q \cdot s+1-q)^n
\end{align}
今回はこのqが違う値を取ると考えればよいので
\begin{align}
GF(s) &= \prod_{q_i}^{q_i \le \sqrt{N}}(q_i \cdot s+1-q_i) \\
& q_i = \frac{1}{p_i^2}なので\\
& = \prod_{p_i}^{p_i \le \sqrt{N}}(\frac{1}{p_i^2} \cdot s+1-\frac{1}{p_i^2}) \cdots (2)
\end{align}
\begin{align}
&\color{red}{RSF(r,N) はGF(s)のs^rの係数}になります
\end{align}
母関数の係数を求める
pythonで$GF(s)$の$s^8$までの係数を求めれば一気に$RSF(r,N) (r=0, \cdots 8)$が求まります。
import sympy as sp
N, R = 10 ** 6, 8
c = [1.0] + [0.0] * R
for prm in sp.primerange(2,N):
p = 1.0 / (prm**2)
q = 1.0 - p
for j in range(R, 0, -1):
c[j] = c[j] * q + c[j - 1] * p
c[0] *= q
for r in range(R+1):
print(f"# RSF({r})={c[r]:10.4e}, SP({r},10^8)={int(c[r]*10**8)}" )
# RSF(0)= 6.079e-01, SP(0,10^8)=60792714
# RSF(1)= 3.354e-01, SP(1,10^8)=33538928
# RSF(2)= 5.329e-02, SP(2,10^8)=5329286
# RSF(3)= 3.292e-03, SP(3,10^8)=329207
# RSF(4)= 9.705e-05, SP(4,10^8)=9704
# RSF(5)= 1.567e-06, SP(5,10^8)=156
# RSF(6)= 1.544e-08, SP(6,10^8)=1
# RSF(7)= 1.001e-10, SP(7,10^8)=0
# RSF(8)= 4.525e-13, SP(8,10^8)=0
プログラムでは$RSF(r,N)$を$p_i \le 10^6$までの素数で求めた値と、前回で求めた$SP(r,10^8)$をこの割合を使って出した値をリストしました。
$r$が小さいときにはそれなりの精度で一致しているのが分かります
ゼータ関数のオイラー積
余談ですが、式(1)を見てゼータ関数のオイラー積を思いついた人はすごいですね。nを無限大までにすると、まさしく$\large \frac{1}{\zeta (2)}$になっています。したがってその値は$\Large \frac{6}{\pi^2}$に収束するはずです。
\begin{align}
& RSF(0, \infty) = \prod_{p}(1-\frac{1}{p^2})= \frac{1}{\zeta (2)} = \frac{6}{\pi^2}
\end{align}
計算してみると有効数字７桁まで一致しているのでかなりの精度です。
import math
print(f"# {6/math.pi**2}")
# 0.6079271018540267
更に（その１）で導いたメビウス関数を使った式をNで割ってNを無限大までにすると、以下のようなメビウス関数とゼータ関数の関係式が$s=2$で導くことができます。
\begin{align}
&\frac{\mathrm{sf(N)}}{N} = \sum_{k=1}^{\infty} \frac{\mu(k)}{k^2} = \prod_{p}(1-\frac{1}{p^2})=\frac{1}{\zeta (2)}\\
\end{align}
もともとゼータ関数は整数の累乗の逆数の和で定義されるので、
\begin{align}
\large \zeta (2) &= \frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}+ \cdots\\
&= \sum_{k=1}^{\infty}\frac{1}{k^2}
\end{align}
累乗の逆数の和をちょっと変形したもの同士が逆数の関係になっているというのは不思議な気がします。
\begin{align}
&\sum_{k=1}^{\infty} \frac{\mu(k)}{k^2} =\frac{1}{\zeta (2)} = \frac{1}{\sum_{k=1}^{\infty}\frac{1}{k^2}}\\
\end{align}
（開発環境：Google Colab）
この考え方はProject Euler、Problem 633: Square Prime Factors IIを解くのに役に立ちます。