SwiftのDecimal型で数学の関数演算（その２）

前回の記事

その１で定義したものを使用するので、見覚えのないものがあったらその１を覗いてみてください。
SwiftのDecimal型で数学の関数演算（その１）

ニュートン法

ニュートン法とは、$f(x) = 0$の解$x$が、$x$の推測値$x_0$が与えられたときに
$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$
の更新式を計算していくと$x_n$が解$x$に近づいていくというものである。
このニュートン法を使うときに難しいのが、$x$のある程度の値を知っている必要があるということだが、今回はなにもフルスクラッチでプログラムを組もうとしているわけではない。あくまでDecimal型での正確な計算結果が出ればいいだけなのである。
つまり、推測値にはdouble型の計算結果使ってもいいよね？

というわけで、以下がニュートン法のプログラムです。

extension Decimal {
    func newton(f_df: (Decimal) -> Decimal, estimate: Decimal) -> Decimal {
        var result = estimate, minDiff = Decimal.greatestFiniteMagnitude
        var next = estimate - f_df(result)
        if next.isNaN { return .nan }
        result = next
        var lambda: Decimal = 1
        while true {
            let step = f_df(result)
            let diff = step.abs
            while diff * lambda >= minDiff {
                lambda /= 2
            }
            minDiff = diff * lambda
            next -= step * lambda
            if next.isNaN || next == result {
                return result
            } else if diff * lambda <= Decimal.leastNonzeroMagnitude {
                return next
            }
            result = next
        }
    }
}

引数f_dfに$\frac{f(x)}{f'(x)}$を計算するクロージャを、estimateに推測値を与えるとニュートン法での計算結果が返ってきます。
ニュートン法では値によっては振動してしまい、無限ループに陥ってしまうので、更新式を
$$x_{n+1} = x_n - \lambda\frac{f(x_n)}{f'(x_n)}$$
と書き換え、$|x_{n+1} - x_{n}|$が以前までより大きくなると$\lambda$の値を小さくして減速させ、値が必ず収束するようにしています。
そして、以降の計算ではこれを使っていきます。

平方根

その１で冪乗のプログラムを書きました。その指数を$\frac{1}{2}$にすれば平方根が求まるけど、平方根はよく使うので、これだけ別により精度のいいプログラムを書いておいた方がいいと思った。
とはいえ平方根は2乗の逆関数なので上記のニュートン法のプログラムを使えば簡単です。

\begin{align}
y &= \sqrt{x} \\
y^2 &= x \\
y^2 - x &= 0 \\
\end{align}

なので、

$$f(y) = y^2 - x = 0$$

としてyをニュートン法で求めればいい。

\begin{align}
f'(y) &= 2y \\
\frac{f(y)}{f'(y)} &= \frac{y^2 - x}{2y}
\end{align}

なので、これとdouble型で求めたxの平方根を推測値としてニュートン法のプログラムに入れてやればいい。

extension Decimal {
    func sqrt() -> Decimal {
        if isNaN || self < 0 { return .nan }

        let dbl = (self as NSNumber).doubleValue
        if dbl.isNaN { return .nan }

        return newton(f_df: {($0 * $0 - self) / (2 * $0)}, estimate: Decimal(Darwin.sqrt(dbl)))
    }
}

超簡単

逆三角関数

これも三角関数の逆関数なので、sin, cos, tanがわかればニュートン法で簡単に計算できます。
更新式の求め方とか平方根の時と同じなので省略します。

extension Decimal {
func atan() -> Decimal {
        if self == 0 { return 0 }
        else if self.isNaN { return .nan }

        let dbl = (self as NSNumber).doubleValue
        if dbl.isNaN { return .nan }

        return newton(f_df: { (result) -> Decimal in
            let cosValue = result.cos(), sinValue = result.sin()
            return (sinValue / cosValue - self) * cosValue.mul(cosValue)
        }, estimate: Decimal(Darwin.atan(dbl)))
    }

    func asin() -> Decimal {
        if self == 0 { return 0 }
        else if self.isNaN || self.abs > 1 { return .nan }

        let dbl = (self as NSNumber).doubleValue
        if dbl.isNaN { return .nan }

        return newton(f_df: { ($0.sin() - self) / $0.cos() }, estimate: Decimal(Darwin.asin(dbl)))
    }

    func acos() -> Decimal {
        if self == 0 { return .pi_two }
        else if self.isNaN || self.abs > 1 { return .nan }

        let dbl = (self as NSNumber).doubleValue
        if dbl.isNaN { return .nan }

        return newton(f_df: { ($0.cos() - self) / ($0.sin() * -1) }, estimate: Decimal(Darwin.acos(dbl)))
    }
}

実は、逆三角関数にもマクローリン展開の公式があるのですが、atanの公式の分母には階乗がないので他のに比べるととにかく収束が遅いです。特に、xが1に近くなると待ってられないほどに収束が遅くなります。なので、マクローリン展開を利用せずにこの方法をとりました。
この理由だとsinとcosはマクローリンでいいはずなんだけど、ニュートン法使ったのは決してめんどくさかっ...あ、はい、すいません。

対数

その１で対数を頑張って求めていたが、expとニュートン法を使うとこれまた簡単に求めることができる。

ln

extension Decimal {
    func ln() -> Decimal {
        if self <= 0 || isNaN { return .nan }
        else if self == 1 { return 0 }

        let dbl = (self as NSNumber).doubleValue
        if dbl.isNaN { return .nan }

        return newton(f_df: { 1 - (self / $0.exp())}, estimate: Decimal(Darwin.log(dbl)))
    }
}

任意の数を入れてチェックしてみるとその１のプログラムよりも精度が高かった。
なんで最初からこうしなかったかと問われると、ニュートン法を使うのはatanを上手く求める方法を探していた時に思いついたもので、その時すでにその１のプログラムが完成していたからで・・・

log

$$\log_ab = \frac{\log_cb}{\log_ca}$$
という公式があるので、上のlnを使えば簡単。
その１で書いた冪乗のプログラムは内部でexpとlnを使っているので、それをニュートン法で繰り返し呼ぶくらいならここでlnを2回呼び出したほうがいいと思う。

extension Decimal {
    func log(base: Decimal) -> Decimal {
        if self <= 0 || base <= 0 || isNaN || base.isNaN { return .nan }
        else if self == 1 { return 0 }
        else if base == 1 { return .nan }
        else if self == base { return 1 }

        return self.ln() / base.ln()
    }
}

etc

その２ではニュートン法を用いた方法を紹介しました。
この方法、マイコンとかで関数の実装をフルスクラッチで行う時には使えないですが、パソコンでdouble型よりも精度の高い結果を求めるのには十分有用だと思います。

しかし、計算に使うアルゴリズムはインターネットで調べているのですが、大体どこのサイトを見てもMathライブラリを使う方法しか書いてなくて大変です。
まあほとんどの人にとってはそれで十分なのかもしれないですけど。
それでも、サイトのタイトルにアルゴリズムと入っているのにMathライブラリの使い方説明しているサイト見つけたときはカチンと来た。
アルゴリズムの説明しろよ。

続き

SwiftのDecimal型で数学の関数演算（その３）

お願い

RPN Anywhere
僕は上記の電卓アプリを公開しており、ここに載せたコードはこの電卓アプリのソースコードから引っ張ってきています。
もし、より正確に速く計算できる方法や、不具合などありましたら教えてください。