対数関数( log )の数値的解法

今回は対数の数値的解法を整理しておこうと思います。まずは自然対数です。

テーラー展開での解法となりますが、ただ \(log(x)\) では０近傍でのテーラー展開ができませんので、工夫のためとりあえず\(log(1+x)\)及び、\(log(1-x)\)を展開します。

\begin{align*} f(x)=log(1+x)\end{align*}

\begin{align*}g(x)=log(1-x)\end{align*}

\begin{align*} f(x)= x - \frac{x^2}{2} + \frac{x^3}{3}-\frac{x^4}{4}+\dots = \sum_{n=1}^{\infty}(-1)^{n-1}\frac{x^n}{n} \tag{1}\end{align*}

\begin{align*} g(x)= -x - \frac{x^2}{2} - \frac{x^3}{3}-\frac{x^4}{4}-\dots = \sum_{n=1}^{\infty}(-1)\frac{x^n}{n} \tag{2}\end{align*}

 (1)から(2)を引きます。

\begin{align*} log\frac{1+x}{1-x}= 2x  + \frac{2x^3}{3} + \frac{2x^5}{5}+\dots = \sum_{n=0}^{\infty}\frac{2x^{(2n+1)}}{(2n+1)} \tag{3}\end{align*}

 ここで

\begin{align*} y=\frac{1+x}{1-x}\end{align*}

とおくと

\begin{align*} x=\frac{y -1}{y+1}\end{align*}

 となることに留意して、(３)を書き換えると

\begin{align*} log(y)= 2\left(\frac{y-1}{y+1}\right) + \frac{2}{3}\left(\frac{y-1}{y+1}\right)^3 + \frac{2}{5}\left(\frac{y-1}{y+1}\right)^5+\dots = \sum_{n=0}^{\infty}\frac{2}{2n+1}\left(\frac{y-1}{y+1}\right)^{(2n+1)} \tag{4}\end{align*}

 といかにもこのままプログラムを書けばよさそうな形に変化させる事が出来ました。しかし\(1-x\)の０近傍というのは、\(x<1\)であることを前提としており(一般二項定理の係数と同じ)\(y\)に直したところで1の近傍である必要があります。ここで平方根の時に出くわした式変形の壁に出くわすのですが、対数の場合は掛け算を足し算に変形できます、具体には 

\begin{align*} log(y) = log(\alpha \cdot 2^{n} ) = n\cdot log(2)+log(\alpha)\end{align*}

 として計算します。\( log(2) \)は常に変わりませんので定数として扱います。では、どうやって\(y\)を\(\alpha \cdot 2^{n} \)の形にするかですが、実数型はそもそもこの形をしているので今回はdoubleを掘り下げるのエントリを応用してプログラムを見てみましょう。

#include <iostream>

typedef union UnDouble_{
    int iNum;
    struct StAna{
        unsigned int lfraction:20;
        unsigned int exponent:11;
        unsigned int sign:1;
    } stAna;
} UnDouble;
 
typedef struct StDouble_{
    UnDouble uLsb;
    UnDouble uMsb;
} StDouble;
 
double log(double x );

#define LN2 0.6931471805599453094172321

int _tmain(int argc, _TCHAR* argv[])
{
    std::cout << log(2.0) << std::endl;
	return 0;
}

double log(double x)
{
	int i = 1, n;
	double a, r, r2;
    StDouble *stAc = (StDouble*)&x;

	if( stAc->uMsb.stAna.sign ) return 0;
	if( x == 0 ) return 0;

	n = stAc->uMsb.stAna.exponent - 1023;
	r = n * LN2;

	stAc->uMsb.stAna.exponent = 1023;

	x = ( x - 1.0 ) / ( x + 1.0 ); 
	a = 2.0 * x; 
	x = x * x;

	do{ 
		r2 = r;
		r = r + a / i;
		a = a * x;
		i = i + 2;
	}while( r != r2 );  

	return r;
}

 ホーナー法を使ってしまっているのでちょっと見通しが悪いですが(４)式をそのままプログラムにしています。バイトオーダが変わると結果が変わってしまうと思うので注意ください。実際は標準関数では割り算のiterationを使用して\(\alpha \cdot 2^{n}\)を作っているようです。なので標準関数より高速なプログラムになっていると思います。また常用対数にしたい場合は底の変換公式 \( log_{10}x = log_{e} x / log_{e}10 \) を使用し\(log10\)を定数として扱えば簡単に実現できます。