対数関数の実用的な表示
$ \log (1+x) = \sum_{n=1}^\infty (-1)^{n-1}\frac{x^n}{n} \quad(-1<x\le1)
xが1に近い場合に収束が遅いので都合が悪い
$ \log 2 の計算例
code:jl
function partial_sum_log1px(x, N::Int)
s = 0.0
for n in 1:N
s += (-1.0)^(n-1) * x^n / n
end
return s
end
partial_sum_log1px(1, 1000) # 0.6926474305598223
log(2) # 0.6931471805599453
数値計算に便利なように変形する
変数を$ -x で置き換える:
code:tex
\begin{aligned}
\log (1-x)
&= \sum_{n=1}^\infty (-1)^{n-1}\frac{(-x)^n}{n}\\
&= - \sum_{n=1}^\infty \frac{x^n}{n}
\end{aligned}
収束域は$ -1\le x < 1 となる
$ \log (1+x) = \sum_{n=1}^\infty (-1)^{n-1}\frac{x^n}{n} と$ \log (1-x) = - \sum_{n=1}^\infty \frac{x^n}{n} の辺々引く
code:tex
\begin{aligned}
\mathrm{LHS}
&= \log (1+x) -\log (1-x)\\
&= \log (1+x)(1-x)^{-1}\\
&= \log \frac{1+x}{1-x}\\
\mathrm{RHS}
&= x -\frac{1}{2}x^2 + \frac{1}{3}x^3 - \frac{1}{4}x^4 + \frac{1}{5}x^5 - \cdots\\
&\quad - \left(-x -\frac{1}{2}x^2 - \frac{1}{3}x^3 - \frac{1}{4}x^4 - \frac{1}{5}x^5 - \cdots \right)\\
&= 2 \left(x + \frac{1}{3}x^3 + \frac{1}{5}x^5 + \frac{1}{7}x^7 \cdots \right)
\end{aligned}
したがって新たな表示を得る:
$ \log \frac{1+x}{1-x} = 2 \left(x + \frac{1}{3}x^3 + \frac{1}{5}x^5 + \frac{1}{7}x^7 \cdots \right) \quad(|x| < 1 )
code:jl
function partial_sum_log_ratio(x, N::Int)
s = 0.0
for k in 1:N
s += x^(2k-1) / (2k-1)
end
return 2.0 * s
end
partial_sum_log_ratio(1/3, 100) # 0.6931471805599451
log(2) # 0.6931471805599453
既知の値を利用して加速する
展開式の変数を$ x = \frac{1}{2p^2-1} \quad (p>1) に従って置換する:
$ \frac{1+x}{1-x} = \frac{p^2}{p^2-1} = \frac{p^2}{(p-1)(p+1)} となる
code:tex
\begin{aligned}
\log \frac{1+x}{1-x}
&= \log \frac{p^2}{(p-1)(p+1)}\\
&= 2\log p - \log (p-1) - \log (p+1)\\
\end{aligned}
これを変形することでさらに収束の早い展開式が得られる:
$ \log p = \frac{1}{2}\log (p-1) + \frac{1}{2}\log (p+1) + \frac{1}{2p^2-1} + \frac{1}{3}\left(\frac{1}{2p^2-1}\right)^3 + \frac{1}{5}\left(\frac{1}{2p^2-1}\right)^5
どういう発想で出てきた？あんも.icon
$ 2\log p - \log (p-1) - \log (p+1) を見据えている？
とりあえず分母を単純な因数の積にすればいい？
code:jl
function partial_sum_log_accel(p, N::Int)
x = 1.0/(2.0*p^2 - 1.0)
s = 0.0
for k in 1:N
s += x^(2k-1) / (2k-1)
end
return 0.5*(log(p-1) + log(p+1)) + s
end
partial_sum_log_accel(3, 5) # 1.098612288668107
log(3) # 1.0986122886681098