令和元年のAdventCalendarの1日目にふさわしいレイワな話

この記事はTech KAYAC Advent Calendar 2019の1日目の記事です。

こんばんわ、今年もこの季節がやってまいりました。
ゲームコミュニティ(Lobi)事業部の徐々に肩書を変化させる謎エンジニアの池田です。
令和元年の1日目・・・

レイ和1年の1日目ということで、この日ふさわしい計算機イプシロンの話をしたいと思います。
※ちなみにですが、予めことわっておくと技術的な話はきわめて少ないです。

前置き

真偽はともかく、ゼロとレイの違いがあるらしい?

0をゼロと読むかレイと読むかで実は違いがあったらしいということを最近知りました。

数字「0」の読みには「ゼロ」と「れい」がある。 ... 「ゼロ」も「れい」も共に「皆無」を意味するが、「零細企業」というように「零(れい)」には「わずか」「きわめて小さい」という意味もある。

「ZERO(ゼロ)」と「零(れい)」 - 違いがわかる事典 より引用

レイはきわめて小さい数のこと

プログラム的な意味でレイ - 計算機イプシロン

計算機イプシロン(けいさんきイプシロン、英: machine epsilon)は、浮動小数点数において、「1より大きい最小の数」と1との差のことである。

計算機イプシロン - Wikipediaより引用

つまり、プログラムが取り扱えるきわめて小さい数なので、計算機イプシロンはレイと言える。

レイ和1年のAdvent Calendarにふさわしいお題

レイ 1ということで・・・
計算機イプシロンを何回足したら1になるのか確かめてみた

まずはgo言語で、計算機イプシロンを計算してみます。
1より大きい最小の数math.Nextafter(1,2)で取得できるので以下のようなコードで出せるでしょう。

package main

import (
    "fmt"
    "math"
)

func main() {
    epsilon := math.Nextafter(1, 2) - 1
    fmt.Printf("epsilon = %g\n", epsilon)
}

The Go Playground

epsilon = 2.220446049250313e-16

wikipedia様曰く IEEE754の倍精度的には2.220e-16 らしいので多分あってるでしょう(ボソ

これで、go言語のレイが手に入りました。

単純に足していった場合

レイ 和 を意識して、掛け算や割り算を使わずに1回ずつ足します。
当然のごとく計算時間は長くなりますので、1分ずつ進捗を出力するようにしておきましょう。
以下のようなコードになりました。

package main

import (
    "log"
    "math"
    "time"
)

func main() {
    //計算機イプシロンの取得
    epsilon := math.Nextafter(1, 2) - 1
    log.Printf("Epsilon = %g\n", epsilon)
    //通知用チャンネルの作成
    doneChan := make(chan bool)
    tickChan := time.NewTicker(time.Second * 60).C
    //計算処理
    var i int
    var total float64 = 0.0
    var start, end time.Time
    log.Printf("total = %f\n", total)
    go func() {
        start = time.Now()
        for i = 0; total < 1.0; i++ {
            total += epsilon
        }
        end = time.Now()
        doneChan <- true
    }()
    //途中経過表示用
    for {
        select {
        case <-tickChan:
            log.Printf("total = %f\n", total)
        case <-doneChan:
            log.Printf("total = %f\n", total)
            log.Printf("time = %f second\n", (end.Sub(start)).Seconds())
            log.Printf("count = %d\n", i)
            return
        }
    }
}

実行してみますと・・・

$ go run main.go
2019/11/28 13:50:44 Epsilon = 2.220446049250313e-16
2019/11/28 13:50:44 total = 0.000000
2019/11/28 13:51:44 total = 0.000004
<以下略>

ん?ちょっと待ってください。1分間でこの進捗ということは・・・

(1 / 0.000004) / 60 = 4166.6666 時間

半年位かかります。

筋肉を使って足していった場合

単純に足していった場合では、実行時間が長すぎました。
仕方がないので、古典的な高速化手法として有名な ループアンローリング を使って実行時間短縮を狙いましょう
先程のコードの

       for i = 0; total < 1.0; i++ {
            total += epsilon
        }

この部分を筋肉を使って以下のように書き換える

       for i = 0; total < 1.0; i += 1000 {
            total += epsilon + //count 10
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon
            total += epsilon + //count 20
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon
            ... <中略> ...
            total += epsilon + //count 1000
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon +
                epsilon
        }

このように、ループの中身を冗長ですが展開することで、終了判定する回数が1/1000回になります。 さて、実行してみると

2019/11/28 13:54:51 Epsilon = 2.220446049250313e-16
2019/11/28 13:54:51 total = 0.000000
2019/11/28 13:55:51 total = 0.000124
<以下略>

(1 / 0.000124) / 60 = 134.4086 時間

約6日、だいぶ短くなりました。
こうやってループアンローリングしてみると、ループの終了判定時に起きる条件分岐のコストって割とあるという実感ができます。

すこーしだけ頭を使った場合

単純に足していった場合や筋肉を使って足していった場合のように、 何も考えずに足していくと途方も無いということが分かりました。
少し頭を使ってみることにします。 以下のコードを実行して、go言語の浮動小数点のbitがどうなっているのかについて考えてみます

package main

import (
    "fmt"
    "math"
)

func main() {
    epsilon := math.Nextafter(1, 2) - 1
    fmt.Printf("epsilon  = %x\n", math.Float64bits(epsilon))
    fmt.Printf("2*epsilon= %x\n", math.Float64bits(epsilon+epsilon))
    fmt.Printf("4*epsilon= %x\n", math.Float64bits(epsilon+epsilon+epsilon+epsilon))
    fmt.Printf("one_next = %x\n", math.Float64bits(math.Nextafter(1, 2)))
    fmt.Printf("one      = %x\n", math.Float64bits(1.0))
}

The Go Playground

epsilon  = 3cb0000000000000
2*epsilon= 3cc0000000000000
4*epsilon= 3cd0000000000000
one_next = 3ff0000000000001
one      = 3ff0000000000000

『何を当たり前な』的な話ですが、このコードの結果を見ると 計算機イプシロンを2の累乗倍していけば、いつかはちょうど1になります。
なら単純に 1/epsilon すればいいじゃないか! つまり、 元のコードの

       for i = 0; total < 1.0; i++ {
            total += epsilon
        }

ここを、すこーし頭を使って

       total = epsilon
        for i = 1; total < 1.0; i += i {
            total += total
        }

このように改変すれば

$ go run main.go
2019/11/28 13:58:20 Epsilon = 2.220446049250313e-16
2019/11/28 13:58:20 total = 0.000000
2019/11/28 13:58:20 total = 1.000000
2019/11/28 13:58:20 time = 0.000000 second
2019/11/28 13:58:20 count = 4503599627370496

すぐに求められるわけです。

おわりに

令和元年のAdventCalendar1日目にふさわしいお題として、計算機イプシロンを何回足したら1になるのか確かめてみた
その結果、 4503599627370496回ということが分かりました。
今回の割り算という高等テクニックを使わずにレイ にちなんで、足し算のみで求めてみたこの試みはいかがでしょうか?
最初の宣言通りに技術的な内容もきわめて少ない内容になってしまったと思います。

物足りないと思う方は、明日の @hirasho0 さんによる『スマホでコンピュートシェーダ』の話を楽しみに待っていただければ幸いです!