スプレッド・シートでフィボナッチ数列を列挙する

あっ，ブログネタ見っけ（笑）

• LibreOffice CalcとExcelの計算結果誤差をみなさんで検証していただいた結果まとめ - Togetter

はっきり言って前半の議論以下の脊髄反射的やり取りは鼻糞ほどの参考にもならないので，最後に紹介されているリンクを読むことを強くオススメする。

• Excelの演算誤差

いや，これ読んでて大昔に特定用途のバンドパスフィルタのシミュレーションを Excel で書かされたことを思い出したよ。 トラウマになるほどではなかったが，あれはなかなか酷い仕事だった（笑）

【事前準備】 フィボナッチ数列の一般項

閑話休題 (それはさておき) ，問題を整理しよう。

まずはフィボナッチ数列の定義から。 $n$ 番目のフィボナッチ数を $F_n$ で表すと

で定義できる。 ちなみに，この定義を素朴に Go 言語コードにしたのが以下。

package main

import "fmt"

var fibonacciNumbers = make(map[int64]int64)

func fibonacciNumber(n int64) int64 {
    switch n {
    case 0:
        return 0
    case 1:
        return 1
    default:
        if fn, ok := fibonacciNumbers[n]; ok {
            return fn
        }
        fn := fibonacciNumber(n-2) + fibonacciNumber(n-1)
        fibonacciNumbers[n] = fn
        return fn
    }
}

func main() {
    fmt.Println("| $n$ | $F_n$ |")
    fmt.Println("| ---:| ---:|")
    for n := int64(1); n <= 75; n++ {
        fmt.Printf("| %d | %d |\n", n, fibonacciNumber(n))
    }
}

このコードを実行して1番目から75番目までのフィボナッチ数列をつくりリファレンス値としよう。 一部を挙げておく。

ところで定義で挙げた数式は「漸化式」と呼ばれるものだが，一般項で表すこともできる。 面倒なので結果だけ Wikipedia から引用してしまおう。 すなわち

とするなら $F_n$ は

と書ける¹。 これをGo 言語で書くとこんな感じ。

package main

import (
    "fmt"
    "math"
)

func main() {
    fmt.Println("| $n$ | $\\varphi^n$ | $\\psi^n$ | $F_n$ |")
    fmt.Println("| ---:| ---:| ---:| ---:|")
    r5 := math.Sqrt(5)
    phi := (1 + math.Sqrt(5)) / 2
    psi := 1 - phi
    for n := int64(1); n <= 75; n++ {
        phin := math.Pow(phi, float64(n))
        psin := math.Pow(psi, float64(n))
        f := (phin - psin) / r5
        fmt.Printf("| %d | %.15f | %.15f | %.15f |\n", n, phin, psin, f)
    }
}

結果も載せておこう。

$\psi^n$ の値が小さすぎてガッツリ情報落ちしているが $F_n$ として必要なのは整数部分のみなので小数点以下を丸めた値が正しければ無問題。 実際のところ

と言えるので

でも十分いけるのだ。

ていうかね。 $\sqrt{5}$ を無理やり浮動小数点数に展開して計算してるんだから誤差が出て当然なのであり，それが不可視化されている各スプレッドシートのほうが相当怪しげな内部処理をしていると言わざるを得ない。

スプレッドシートでフィボナッチ数列の一般項を計算する

さて，ようやく本題。

つまるところ，この一般項を使ってスプレッドシートでフィボナッチ数列を計算させたらどうなるか，というのが今回のお題だ。

1. まずA列に 1, 2, 3, 4,… と整数値を入れておく
2. B列には式 =((1+SQRT(5))/2)^A1 (1行目の場合) をセットする（$\varphi^n$ 相当）
3. C列には式 =((1-SQRT(5))/2)^A1 (1行目の場合) をセットする（$\psi^n$ 相当）
4. D列には式 =((((1+SQRT(5))/2)^A1)-(((1-SQRT(5))/2)^A1))/SQRT(5) (1行目の場合) をセットする

これで B列 → C列 → D列 と計算過程を追うことができる。

まずは LibreOffice Calc で実行してみた。 ちなみにバージョン6.2系である。

• fibonacci-numbers.ods

結果の一部，70番目以降のみ挙げておく。

まず小数点以下の値が潰れているのが懸念点なのだが（これについては後述する），表からは71番目からはリファレンス値に対して明らかにズレが生じていることが分かる。

同じことを Google スプレッドシートでもやってみた。

• fibonacci-numbers

これも70番目以降のみ挙げておこう。

途中計算に若干の違いはあれど最終的な結果は同じになった。

Excel のことは知らない²。

Machine Epsilon

ある処理系における浮動小数点数の精度を調べるには $1+\epsilon-1$ がゼロにならない最小の $\epsilon$ すなわち machine epsilon を調べるのがよい。 今回使っている処理系では，浮動小数点数の基数はどれも2進数なので $\epsilon$ の値を

として調べていけばいいだろう。 たとえば Go 言語であれば

package main

import (
    "fmt"
    "math"
)

func main() {
    fmt.Println("| $n$ | $\\epsilon$ | $1+\\epsilon-1$ |")
    fmt.Println("| ---:| ---:| ---:|")
    for n := int64(1); n <= 55; n++ {
        e := math.Pow(2.0, float64(-n))
        me := 1.0 + e - 1.0
        fmt.Printf("| %d | %v | %v |\n", n, e, me)
    }
}

といったコード。 このコードの実行結果（の一部）は以下の通り。

64ビット浮動小数点数の仮数部のサイズは52ビットなので，これは妥当な結果と言える。

同じことを LibreOffice Calc で実行した結果は以下の通り。

Excel で聞かれる括弧の付け方で結果が変わる云々はなかったが $\epsilon=2^{-48}$ までの精度しかないようだ。 これなら先のフィボナッチ数列一般項の計算結果も納得である。 何がネックになってるのかは分からないが。

Google スプレッドシートはちょっと複雑で，まず冪乗計算（POW(x,y) または x^y）が $2^{-33}$ までしか値が取れない。 しょうがないので，計算値ではなく，生の $\epsilon$ 値をセルにコピペして検証してみたところ $\epsilon=2^{-52}$ までの精度があることが確認できた³。 しかし関数ごとに精度が異なるのであれば，これはこれでヒッジョーに面倒くさい話になる。

というわけで

結論として（大きな有効桁数が要件となる）数値計算をスプレッドシートで行うのはオススメできない。 やるのであれば身元の確かなツールを使って数値計算を行い，その結果をスプレッドシートにインポートしてデータの整理を行うことであろうか。

ブックマーク

• 1を1億回足して1億にならない場合

参考図書

数学ガール

結城 浩 (著)

SBクリエイティブ 2007-06-26 (Release 2014-03-12)

Kindle版

B00EYXMA9I (ASIN)

評価     

ミルカさんとの衝撃の encounter。数学ガールがワルツを踊る。

reviewed by Spiegel on 2014-02-14 (powered by PA-APIv5)

プログラマの数学 第2版

結城 浩 (著)

SBクリエイティブ 2018-01-16 (Release 2018-02-08)

Kindle版

B079JLW5YN (ASIN)

評価     

タイトル通りプログラマ必読書。第2版では機械学習に関する章が付録に追加された。

reviewed by Spiegel on 2018-03-19 (powered by PA-APIv5)