Goで高速逆平方根計算アルゴリズム - Fast Inverse Square Root -

この記事は東京理科大学 Advent Calendar 2018 1日目の記事です。

はじめに

東京理科大学理学部数理情報科学科（現：応用数学科）、非線形の最適化理論をメインで専攻して2010年に卒業した理科大OBの @kaneshin です。卒業したあとは組込み系企業に入社し、そこからカナダ留学を経て、現在はPairsを開発・運営しているエウレカのCTOをやっています。理科大卒だとあまり見ない経歴な気がします。もしかしたら最近の理科大生にもWeb系の風が吹いているのか、東京理科大学 Advent Calendar 2018を見つけたときには少し気持ちが昂ぶりました。

最近の仕事内容については eureka Advent Calendar 2018 の1日目の記事に投稿しているので、ぜひご覧ください。

medium.com

その上でGoを書くことが普段は多いので、この記事ではGoと数値解析を合わせことを紹介します。

Fast Inverse Square Root from Quake III Arena

Quake III Arenaというゲームで有名な平方根の逆数の近似を高速に求めることができるアルゴリズムです。英語版 Wikipedia で『Fast Inverse Square Root』のページを開くと詳細が書いてあります。

Fast inverse square root - Wikipedia

The algorithm is best known for its implementation in 1999 in the source code of Quake III Arena, a first-person shooter video game that made heavy use of 3D graphics.

このゲームに光の反射の計算に、この高速逆平方根計算アルゴリズムのFast Inverse Square Rootが採用されています。

How: どのように求解しているのか

このアルゴリズムは精度を犠牲にし、高速に近似値を求めることができるアルゴリズムになっています。数値解析を勉強した事がある人は非線形方程式の近似解をニュートン法を用いて計算したり、そのニュートン法の計算回数やパラメータチューニングをしてレポート提出したことがあるかと思います。

この高速逆平方根計算アルゴリズムはニュートン法の考えをベースに、32bitの float と int の casting pointers をして、floatの仮数部と指数部を活用して一回のニュートン法で高速に解いています。

ここで使っているのが、マジックナンバー（いろんな意味でマジック）の 0x5f3759df を用いています。

Fast Inverse Square Root Algorithm in Go

FastInverseSqrt という関数と、それのテストコードを一緒に記載しています。

package main

import (
    "fmt"
    "math"
    "testing"
    "unsafe"
)

func FastInverseSqrt(x float32) float32 {
    const threehalfs = float32(1.5)

    x2 := x * float32(0.5)
    y := x
    i := *(*int32)(unsafe.Pointer(&y))
    i = 0x5f3759df - i>>1
    y = *(*float32)(unsafe.Pointer(&i))
    y = y * (threehalfs - (x2 * y * y))
    return y
}

func TestFastInverseSqrt(t *testing.T) {
    tests := []struct {
        in, out float32
        e       float32
    }{
        {2.0, 1.0/1.414, 1e-3},
        {3.0, 1.0/1.732, 1e-3},
        {5.0, 1.0/2.236, 1e-3},
    }

    for _, tt := range tests {
        tt := tt
        t.Run(fmt.Sprintf("1/√%.2f", tt.in), func(t *testing.T) {
            result := FastInverseSqrt(tt.in)
            e := float32(math.Abs(float64(result - tt.out)))
            if tt.e < e {
                t.Errorf("want %f, but got %f", tt.out, result)
            }
        })

    }
}

See: https://play.golang.org/p/C0KNGXBb14p

ポインタをキャスティングするための unsafe パッケージ

32bitの float と int をそれぞれポインタでキャスティングが必要になるので、 unsafe.Pointer が必要になります。

i := *(*int32)(unsafe.Pointer(&y))
i = 0x5f3759df - i>>1
y = *(*float32)(unsafe.Pointer(&i))

おわりに

数値解析のアルゴリズムの詳細についてはまた数値解析 Advent Calendar 2018で書こうと思います。

それでは、明日からの理科大関係者がどのような記事を書くのか楽しみにして12月を過ごしていければと思います。

at kaneshin

Free space for me.