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

Double型とDecimal型

Double型は2の累乗を用いて値を表すため、10進数の小数を扱うと誤差が出やすい。また、有効桁数も15桁と少なく、精度もよろしくない。
そこで、電卓などのユーザーに計算した結果を直接与えるアプリではDecimal型が好まれる。
Decimal型は10の累乗を用いて値を表すため、10進数を誤差なく表すことができ、また有効桁数も38桁と多いので、精度がいい。

しかし、Decimal型にはMathライブラリにあたる数学の関数演算を行う機能が存在しない。
せっかく精度よく計算できるのに、Mathライブラリを利用するためにDouble型に変換しては意味がない。

そこで、Decimal型のまま数学の関数演算を頑張ってみた。

頻繁に使う関数

まず簡単だけど計算によく使う関数を定義します。計算量も多くなく、引数も取らないものはComputed Propertyで定義しています。
また、わかりやすいようにいちいちextensionでくくってますが、実際は1つのextensionのブロックの中にコードを書いてます。

絶対値

extension Decimal {
    var abs: Decimal {
        return self < 0 ? -self : self
    }
}

簡単ですね。

床関数

extension Decimal {
    var floor: Decimal {
        let behavior = NSDecimalNumberHandler(roundingMode: .down, scale: 0, raiseOnExactness: false, raiseOnOverflow: false, raiseOnUnderflow: false, raiseOnDivideByZero: false)
        return (self as NSDecimalNumber).rounding(accordingToBehavior: behavior) as Decimal
    }
}

これもNSDecimalNumberの機能を使ってるだけなので、特に問題ないはずです。

掛け算のアンダーフロー

Decimal型ではオーバーフローやアンダーフローを起こすと値がNaNになるのですが、なぜか掛け算を行い、アンダーフローを起こした時だけNaNにならず大きな値が返ってきます。そこで、Decimal型の掛け算は*をそのまま利用せず、積をとった後にアンダーフローチェックを行うメソッドを作ってそれを利用した方がいいと思います。僕の場合は以下のようなメソッドを定義して、掛け算の際はこれを利用しています。

extension Decimal {
    func mul(_ y: Decimal) -> Decimal {
        let r = self * y
        if self.abs < 1 && y.abs < 1 && r.abs >= 1 {
            return .nan
        }
        return r
    }
}

ネイピア数

e

ネイピア数のeがDecimal型に定義されていないが、頻繁に利用するので定義します。

extension Decimal {
    static let e = Decimal(string: "2.71828182845904523536028747135266249776")!
}

e^x

ネイピア数の冪乗です。底がe以外だと簡単には計算できないのですが、底がeのときだけマクローリン展開を用いて計算できます。

$$e^x = \sum^\infty_{n = 0} \frac{x^n}{n!}$$

なので、for文で$\frac{x^n}{n!}$の項が十分小さくなるまで足し合わせていけば計算できます。
ただし、これには解が収束する速度というのがありまして、xが小さいほど速く収束する、つまり速く計算できます。
また、なぜかDecimal型は指数が自然数の時だけpow関数が使えるので、指数を整数部と小数部に分け、整数部はpow関数で計算し、小数部はマクローリン展開で計算、後で2つを掛け合わせる方法をとりました。

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

        // Separate Integer and Decimal
        let integer = self.floor
        let decimal = self - integer

        // Calculate Integer Part
        let intX = (integer as NSNumber).intValue
        if !integer.isZero && intX == 0 { return .nan }
        let intRes = Foundation.pow(Decimal.e, intX)

        // Calculate Decimal Part
        var coef = decimal
        var coefs = [Decimal]()
        for i in 2...34 {
            coef = coef.mul(decimal / Decimal(i))
            if coef.isNaN { break }
            coefs.append(coef)
        }
        let decRes = coefs.reversed().reduce(0, { $0 + $1 }) + decimal + 1

        // Merge
        return intRes * decRes
    }
}

まず最初にNaNチェックを行っています。指数がNaNの時に計算しても意味ないですからね。
次に指数が0かどうかチェックしてます。指数が0のときは底の値に関わらず常に解が1ですから計算する必要ないのです。
その次に指数が負かどうかチェックしてます。pow関数が自然数しか扱えないので負の数は計算ができません。しかし、指数が負のときの解は指数が正の時の解の逆数なので、変換してやることで計算ができます。

次に整数部と小数部に分けています。上で定義した床関数使ってます。

次に整数部を計算しています。pow関数を使うのですが、指数部はInt型なので指数部だけInt型に変換してます。その際、Double型の方が大きな値を扱えるので値によっては変換に失敗しますが、そんなに指数の大きな冪乗なんてオーバーフローするに決まってるので変換に失敗したら計算はやめましょう。

次に小数部を計算しています。nが34までの項を順番に計算してそれを逆順に足し合わせています。なぜ34かというと、Decimal型の有効桁数は38桁で、必ず初項に1があるのでここの値は1以上になります。（最初に負数のチェックがあるので、指数が負になり振動することもない。）もし仮に1だとして、有効桁数38桁だと$10^{-37}$の桁までが正確に表現されます。つまり、1に$10^{-38}$を加算したところでDecimal型では意味がないのです。ここで、$n = 34$の項を見ると、$\frac{x^{34}}{34!}$となっています。ここで、$x \lt 1$なので$x^{34} \lt 1$です。仮に、xが最大値の1だとして、
$$\frac{1}{34!} \fallingdotseq 3.39 \times 10^{-39}$$
なので、このくらいまで計算すれば大丈夫だろうと思ってこのようになりました。
その後、情報落ち対策として逆順にこれらの値を加算して小数部の計算は完了です。

あとは整数部と小数部の値を掛け合わせれば答えが出ます。ここではアンダーフローを起こすことがないのでそのまま積演算を行っています。
Decimal型は38桁扱えますが、37桁目くらいまでは正確に値が出てるのではないでしょうか。

三角関数

三角関数も$e^x$と同様にマクローリン展開で求めます。ちなみに角度の単位はラジアンなので、必要に応じて変換してください。

円周率

$\pi$はDecimal型に定義されてるけど、$2\pi$と$\frac{\pi}{2}$があったら便利なので定義しておきます。

extension Decimal {
    static let two_pi = Decimal.pi * 2
    static let pi_two = Decimal.pi / 2
}

正規化

三角関数のマクローリン展開も$e^x$同様に値が小さいほど収束が速くなり、計算も速くなります。そこで、値をなるべく小さな値にしてから計算をします。

まずは三角関数は$2\pi$で周期的なので値を$0 \le x < 2\pi$の範囲内に変換します。

extension Decimal {
     func nomalizeRadian() -> Decimal {
        let t = (self / (.pi * 2)).floor
        return self - (t * (.pi * 2))
    }
}

また、さらに三角関数は全て$0 \le x \le \frac{\pi}{2}$の計算に置き換えられるのでこれを利用します。

範囲	$\sin x$の置換	$\cos x$の置換
$0 \le x \le \frac{\pi}{2}$	$\sin x$	$\cos x$
$\frac{\pi}{2} < x \le \pi$	$\sin (\pi - x)$	$-\cos (\pi - x)$
$\pi < x \le \frac{3}{2}\pi$	$-\sin (x - \pi)$	$-\cos (x - \pi)$
$\frac{3}{2}\pi < x < 2\pi$	$-\sin (2\pi - x)$	$\cos (2\pi - x)$

sin

公式が
$$\sin x = \sum^\infty_{n = 0} \frac{(-1)^n}{(2n+1)!}x^{2n+1}$$
なので、これを計算するだけです。

extension Decimal {
    func sin() -> Decimal {
        if isNaN { return .nan }
        var angle = nomalizeRadian()
        var sign: Decimal = 1
        if .pi_two < angle && angle <= .pi { angle = .pi - angle }
        else if .pi < angle && angle <= .pi_two * 3 { angle -= .pi; sign = -1 }
        else if .pi_two * 3 < angle { angle = .two_pi - angle; sign = -1 }

        let square = angle.mul(angle)
        var coef = angle
        var coefs = [coef]
        for i in 1...19 {
            coef = coef.mul(-square / Decimal(2 * i * (2 * i + 1)))
            if coef.isNaN { break }
            coefs.append(coef)
        }
        let res = coefs.reversed().reduce(0, { $0 + $1 })
        return res * sign
    }
}

$e^x$とほとんど変わらないので説明は省略します。

cos

公式が
$$\cos x = \sum^\infty_{n = 0} \frac{(-1)^n}{(2n)!}x^{2n}$$
で計算できるので、これもsinとほとんど変わりません。

extension Decimal {
    func cos() -> Decimal {
        if isNaN { return .nan }
        var angle = nomalizeRadian()
        var sign: Decimal = 1
        if .pi_two < angle && angle <= .pi { angle = .pi - angle; sign = -1 }
        else if .pi < angle && angle <= .pi_two * 3 { angle -= .pi; sign = -1 }
        else if .pi_two * 3 < angle { angle = .two_pi - angle }

        let square = angle.mul(angle)
        var coef: Decimal = 1
        var coefs = [coef]
        for i in 1...19 {
            coef = coef.mul(-square / Decimal((2 * i - 1) * 2 * i))
            if coef.isNaN { break }
            coefs.append(coef)
        }
        let res = coefs.reversed().reduce(0, { $0 + $1 })
        return res * sign
    }
}

こちらも説明は省略します。

tan

extension Decimal {
    func tan() -> Decimal {
        return sin() / cos()
    }
}

ごめんなさい、手を抜きました。その分精度も落ちます。
気が向いたらちゃんとマクローリン使ったプログラムにします。
2018.3.12更新
tanのマクローリン展開の公式を使ったプログラム書いてみましたが、精度が逆に落ちてしまいました。
詳しくはこちらを参照
SwiftのDecimalでtanをマクローリン展開で求める

対数

2018.3.13更新
その２でニュートン法を使って対数を求めるプログラムを書きました。
そっちの方が精度良く、簡単に計算できるのでいいかもしれません。

log

とりあえずプログラムを載せます。

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 < 1 { return (1 / self).log(base: base) * -1 }
        else if base < 1 { return log(base: 1 / base) * -1 }
        else if self < base { return 1 / base.log(base: self) }

        var number = self
        var coef: Decimal = 1
        var result: Decimal = 0, next: Decimal = 0
        while true {
            if number == base {
                result += coef
                break
            } else if number == 1 {
                break
            } else if number > base {
                next += coef
                number /= base
                if next.isNaN || next == result {
                    break
                }
                result = next
            } else {
                coef /= 2
                number = number.mul(number)
            }
        }
        return result
    }
}

自分でプログラム書いておいてなんですが、これについてはよく理解していません。
参考にさせていただいたサイトもあったのですが、そのサイトがどこにあるか分からなくなってしまいました。
とりあえず軽く分かっていることだけ説明すると、
$$\log_ab \ (a \ge b > 1)$$
という数式があった時、
$$\log_ab = \log_aa\cdot\frac{b}{a} = \log_aa + \log_a\frac{b}{a} = 1 + \log_a\frac{b}{a}$$
と変形できます。しかし、これを繰り返していくとlogがマイナスに発散してしまうので、
$$\log_ab = \log_a(b^2)^\frac{1}{2} = \frac{1}{2}\log_ab^2$$
を利用して発散を抑えつつ、logの係数を小さくして収束させます。ちなみに、これは$b < 1$となると意味を成さないので、そうなる前、つまり$b < a$になった時点で行う必要があります。
これはループの回る回数が予測できなかったので、順次加算していき、値を加算しても合計値に変化がなくなった時点でループを終わらせています。
最初のif文は条件$a \ge b > 1$に合わせるために逆数を使って変形を行っているだけです。

ln

上記のメソッドに底eを与えるだけです。

extension Decimal {
    func ln() -> Decimal {
        return log(base: .e)
    }
}

冪乗

冪乗は今までみたいに無限級数を利用して導出することができません。しかし、以下のような式変形をすれば、今までに作ってきたプログラムを組み合わせれば導出することができます。

$$x^y = e^{log_ex^y} = e^{y \cdot log_ex}$$

よってこうなりました。

extension Decimal {
    func pow(exp: Decimal) -> Decimal {
        if isNaN || exp.isNaN { return .nan }
        else if exp == 0 { return 1 }
        else if exp < 0 { return 1 / pow(exp: exp * -1) }

        // Separate Integer and Decimal
        let integer = exp.floor
        let decimal = exp - integer

        // Calculate Integer Part
        let intX = (integer as NSNumber).intValue
        if !integer.isZero && intX == 0 { return .nan }
        let intRes = Foundation.pow(self, intX)

        // Calculate Decimal Part
        let decRes = (decimal * ln()).exp()

        // Merge
        return intRes * decRes
    }
}

$e^x$のときと同じように整数部と小数部に分けてますが、小数部の計算は上の等式そのままです。

etc

続き

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

お願い

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