Math - つがいは一緒に面倒見るのが吉な件

というわけで続きです。

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

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

演算におけるつがい

つがいになると、大事なことは一人ではなく相方と一緒にやることになります。映画鑑賞から子育てまで。「リア充爆発しろ」?つがいから見たら、よく「一人でやってて爆発しませんね」といったところかもしれません。

なぜつがいがつがいか。一人でやるよか楽(だ|しい)からに決まってるじゃないですか。

実は演算においても、そうしたつがいは少なくないのです。一番有名なのは除算における商(quotient)と余り(remainder)ですか。これはどちらかを計算しようとすると必ずもう片方も出てきます。

実はCにもdivという関数があって、これを使うとつがいのまま得られます。

man div

The div() function computes the value numer/denom (numerator/denominator). It returns a structure named div_t that contains two int membersnamed quot (quotient) and rem (remainder).

/と%でどちらか片方を捨ててしまっているというのは、実は結構もったない話でもあるのですね。何しろ除算というのは、四則演算の中で圧倒的に重い演算なのですから。あまりに遅いので、現代のコンパイラーではこっそり裏で可能な限り逆数との乗算に直しているぐらいです。名著「ハッカーのたのしみ」にも、その手の工夫がたくさん出てきます。

PONSも例外ではなく、除算の基本演算は/ではなく商と剰余の両方を返すdivmodで、/はdivmodの商のみを返すという実装になっています。

三角関係、もとい三角関数におけるつがい

演算のつがいは、もちろん商と剰余に限りません。cosとsinなんかまさにそうです。「虚数の情緒」を読めばその仲睦まじさに惚れ惚れすること間違いなしなのですが、本を読むのが苦手でも今日日のWebには以下のような動画がいくらでも転がっています。ちなみに以下はWikipediaから。

で、PONSではどうしていたかというと、あろうことか別々に求めていました。

Swift - Introducing PONS = Protocol Oriented Number System

数学者でなくとも、これはテイラー展開

ほぼそのままであることが見て取れるかと思いますが

これ、実はそのまま使うのはナイーブすぎるのです。

上記のテイラー展開は必ず収束しますが、どれくらい早く収束するかは原点にどれだけ近いかで決まってきます。原点に近い、つまり小さければ小さいほど早い。特に1より小さいことは重要で、1より小さいことが分かっていれば収束判定で見るべきは各項の分母である階乗だけで済みます。

つまり、1より小さいように正規化した上で、その正規化した値に対する答えを求めればいいわけです。

そのためのwrapAngleだったわけですが、しかしwrapAngleの正規化は、±πであって±1ではありません。3.14にしろ3.141592653589793にしろ、1よりずっと大きい。どうしたものか…

ここで思い出したのが、cosとsinがつがいであるという事実。現在のPONSでは、この両者は以下の通り一緒に面倒を見ています。

public static func sincos(x:Self, precision px:Int = 64)->(sin:Self, cos:Self) {
    if let dx = x as? Double { return (Self(Double.sin(dx)), Self(Double.cos(dx)))}
    if x.isZero || x.isInfinite || x.isNaN {
        return (Self(Double.sin(x.toDouble())), Self(Double.cos(x.toDouble())))
    }
    let epsilon = Self(BigFloat(significand:1, exponent:-px))
    if x * x <= epsilon {
        return (x, 1)   // sin(x) == x below this point
    }
    let atan1   = pi_4(px)
    let sqrt1_2 = sqrt2(px)/2
    func inner_cossin(x:Self)->(Self, Self) {
        if 1 < x.abs {  // use double-angle formula to reduce x
            let (c, s) = inner_cossin(x/2)
            if c == s { return (0, 1) } // prevent error accumulation
            return (c*c - s*s, 2*s*c)
        }
        if x.abs == atan1 {
            return (x.isSignMinus ? -sqrt1_2 : +sqrt1_2, +sqrt1_2)
        }
        var (c, s) = (Self(0), Self(0))
        var (n, d) = (Self(1), Self(1))
        for i in 0...px {
            var t = n.divide(d, precision:px)
            t.truncate(px)
            if i & 1 == 0 {
                c += i & 2 == 2 ? -t : +t
            } else {
                s += i & 2 == 2 ? -t : +t
            }
            if px < d.precision { break }
            n *= x
            d *= Self(i+1)
        }
        return (c, s)
    }
    var (c, s) = inner_cossin(x.abs < 8 ? x : wrapAngle(x, precision:px))
    return (s.truncate(px), c.truncate(px))
}

1より大きい?よろしい、ならば倍角の公式だ。

見ての通り、1より大きな場合、半分の角度の場合を求めてから倍角の公式に入れているわけですが、この倍角の公式の入場資格は、cosとsinがつがいであること。

仮にバラバラに計算したとなると、角度を半分にする操作ごとに呼び出す関数が倍になってずっと遅くなってしまいます。その恐ろしさをフィボナッチ数を再帰的に解いてみて実感された方も少なくないのではないのでしょうか。

それでも8以上の場合はwrapAngleで角を正規化していますが、これは倍角の公式適用時に、桁数が倍になるため。PONSの実数演算は、明示しない限りロスレス。掛け算で有効桁数が倍になってもそのままなのです。さすがに野放途に大きくできないので「適当なところで」wrapAngleに助けてもらっているのですが、現在その閾値は2πより大きな最小の2のべき数である8にしてあります。

そして同様の工夫を、三角関数よりマイナーで、かつてはmath.hも標準サポートしていなった双曲線関数でもやっています。1を境にexpに計算を任せるのをやめて、sincosと同様、ペアテイラー展開でcoshとsinhを同時に求めています。

ところでcosとsinだと、cosが♀でsinが♂だという電波を受信しました。確かにcosは偶関数でsinは奇関数、複素関数exp(z)ではcosが実部、sinが虚部担当なのであながち間違いでもなさそうです。

演算の世界では、リア充の方がかえって(計算)爆発しにくいというお話でした。

Dan the Husband