Math - 円周率を3.14どころか3.141592653589793としても丸く収まらなかった件 #algorithm

算数の問題「円周率を3.14とするとき、半径11の円の面積を求めよ」の解を379.94とするのは誤り？なんて話題が数日前にあったようですが、ちょうどその頃、3.14どころか3.141592653589793にしても不十分だったという体験をしたので。

ちなみに3.141592653589793というのは、おそらく過半のコンピューター「πを〜と」している数値。16進表記では0x1.921fb54442d18p+1となります。

% irb
irb(main):001:0> Math::PI
=> 3.141592653589793

今、私はPONSで数(値型)を再発明している真っ最中なのですが、そのPONSにおける実数型==PORealプロトコルに準拠する型にはもれなくwrapAngle()という関数が付いてきます。

要はその数値XをXラディアンとみなして、それを±πの範囲に正規化する関数で、角度を扱う際に重宝するのですが、なぜ単純に

normalized = original % (2*pi)

とかではダメなのでしょう?

問題は、上記のpiは円周率そのものではなく、近似値に過ぎないということそのものにあります。本物より最悪 epsilon / 2 ずれているのですから、「一回転」ごとにこれが蓄積することになります。一回転や二回転ならとにかく、何万、何億、何兆回転もすればもう「元の角度」なんかわからなくなります。非現実的な数字?なにをおっしゃいます。たかが一年でも1000万π秒もあるのですぞ。

「それってatan2(cos(theta), sin(theta))とかで行けね?」といった方、さすがですね。しかしPONSでは、そのcosとsinこそがwrapAngleを必要としていたのでした?そのためにcosとsinを使うのでは本末転倒です。さあどう数る?

結局PONSで採用したの、ある意味最も豪快な方法でした。

wrapAngle

 public static func wrapAngle(x:Self, precision px:Int = 64)->Self {
        var angle = x
        let onepi = pi(px)
        if angle < -2*onepi || +2*onepi < angle {
            let precision = px + angle.toMixed().0.msbAt
            let twopi = 2*pi(precision)
            angle = angle % twopi
            angle.truncate(px)
        }
        if angle < -onepi { angle += 2*onepi }
        if +onepi < angle { angle -= 2*onepi }
        return angle
    }

なんのことはない、「M_PIが正しくなければ、もっといい値を使えばいいのに」というわけです。そのために必要とあらば、PONSはπを必要な精度まで計算しなおします。どれくらいの精度が必要?n回転ならlog2(n)倍ですか。見ての通り、n回転ではなくnラジアンで見積もっているので、必要な精度よりも少し過剰気味に計算していることにはなりますが、そのためにわざわざlog2(n)を計算するより、整数部分のビットを数えた方がずっと安価なのでそうしています。

余談ではありますが、PONSではπをはじめとする定数はメモ化してます。精度が足りない時だけ再計算し、足りていればキャッシュから低精度版を生成しています。ここでも面白いノウハウを見つけたのですが、別記事にしましょう。

そうした工夫の甲斐あって、PONS搭載の初等関数は、±DBL_MINから± DBL_MAXに至るまで、libcの初等関数と倍精度で比較して最終ビット以外すべて一致します。

REPL

29> Double.cos(DBL_MAX) 
$R23: Double = -0.99998768942655991
 30> BigFloat.cos(BigFloat(DBL_MAX)).toDouble()
$R24: Double = -0.99998768942655991

「最終ビット以外」というのは、PONSの方が高精度で計算する分、丸め誤差が健在化するため。そのため test suite を書くのも一苦労ではありました。

DoubleにとってはDBL_MAXの先に待っているのはinfしかないのですが、PONSのBigFloatなら、その先にも行けます。

REPL

 42> DBL_MAX*2 
$R35: Double = +Inf
 43> (BigFloat(DBL_MAX)*2).description 
$R36: String = "359538626972463141629054847463408713596141135051689993197834953606314521560057077521179117265533756343080917907028764928468642653778928365536935093407075033972099821153102564152490980180778657888151737016910267884609166473806445896331617118664246696549595652408289446337476354361838599762500808052368249716736.0"

REPL

 44> Double.cos(DBL_MAX*2) 
$R37: Double = NaN
 45> BigFloat.cos(BigFloat(DBL_MAX)*2).toDouble() 
$R38: Double = 0.99995075800934008

とりあえず倍先まで行ってみました。計画通り!

それではなぜPONSの三角関数がwrapAngleを利用するようになったのか。続きは次の記事で。

Dan the Transcendental Coder