平方根を高速で近似したい！

こんにちは！いかついプログラマです。

ゲームを作っていると、欲しい数値を割り出すのに平方根(√)を使用したいことがあります。ただ、平方根は他の演算と比べて負荷が高いことが一般的に知られているかと思います。 (Unityなんかでも平方根を使用しないVector3.sqrMagnitudeがわざわざ用意されていたりもしますよね。)

距離の比較などだと平方根を回避することも可能ですが、そうはいかない場合も多々あるかと思います。そこでインターネットの海を彷徨っていたところ…Fast inverse square root[1]なる実装に出会いました！実装としては、ニュートン法による近似で1/√x を求めるものとなっています。(√xはx * 1/√xから簡単に求められます)

c#で実装してみるとこのようになります。

 [StructLayout(LayoutKind.Explicit)]
    public struct FloatIntUnion
    {
        [FieldOffset(0)]
        public float f;

        [FieldOffset(0)]
        public int i;
    }

    public static float FastInverseSqrt(float x)
    {
        if (x == 0)
        {
            return 0;
        }
        float xHalf = 0.5f * x;
        FloatIntUnion floatIntUnion;
        floatIntUnion.i = 0;
        floatIntUnion.f = x;
        floatIntUnion.i = 0x5F3759DF - (floatIntUnion.i >> 1);
        floatIntUnion.f *= (1.5f - (xHalf * floatIntUnion.f * floatIntUnion.f));
        // 繰り返すと精度があがる
        // floatIntUnion.f *= (1.5f - (xHalf * floatIntUnion.f * floatIntUnion.f));
        return floatIntUnion.f;
    }

floatIntUnion.i = 0x5F3759DF - (floatIntUnion.i >> 1);

この処理が肝となっていて、マジックナンバーの親玉みたいな数字 0x5F3759DF がひと際目を引きますね…

これらはfloatが32bit浮動小数点数であることを利用してやるとうまく導出できますが、ここでは省略（fast sqrtや0x5F3759DFでググると導出過程に出会えると思います）。

本題は、「ホントにこの実装は高速なのか？」です。

計測

以下のようなテストコードを用意し、処理時間を計測します。今回は
・UnityEngine.Mathf.Sqrt

・System.MathF.Sqrt

・FastSqrt(上記の実装)
の3つの処理を比較します。

using UnityEngine;
using System.Runtime.InteropServices;


public static class MathUtility
{
    [StructLayout(LayoutKind.Explicit)]
    public struct FloatIntUnion
    {
        [FieldOffset(0)]
        public float f;

        [FieldOffset(0)]
        public int i;
    }

    public static float FastSqrt(float x)
    {
        if (x == 0)
        {
            return 0;
        }
        float xHalf = 0.5f * x;
        FloatIntUnion floatIntUnion;
        floatIntUnion.i = 0;
        floatIntUnion.f = x;
        floatIntUnion.i = 0x5F3759DF - (floatIntUnion.i >> 1);
        floatIntUnion.f *= (1.5f - (xHalf * floatIntUnion.f * floatIntUnion.f));
        return floatIntUnion.f * x;
    }

    public static void TestSqrt()
    {
        var stopwatch = new System.Diagnostics.Stopwatch();
        stopwatch.Start();
        for (int i = 0; i < 100000000; i++)
        {
            FastSqrt(i * 0.01f);
        }
        stopwatch.Stop();
        var fastTime = stopwatch.Elapsed;
        
        stopwatch.Reset();
        stopwatch.Start();
        for (int i = 0; i < 100000000; i++)
        {
            Mathf.Sqrt(i * 0.01f);
        }
        stopwatch.Stop();
        var normalTime = stopwatch.Elapsed;

        stopwatch.Reset();
        stopwatch.Start();
        for (int i = 0; i < 100000000; i++)
        {
            System.MathF.Sqrt(i * 0.01f);
        }
        stopwatch.Stop();
        var systemTime = stopwatch.Elapsed;
        Debug.LogAssertion( $" UnityEngine.Mathf.Sqrt time = {normalTime} Sqrt(2) = {Mathf.Sqrt(2)}\n" +
                            $" System.MathF.Sqrt time = {systemTime} Sqrt(2) = {System.MathF.Sqrt(2)}\n" + 
                            $" FastSqrt time = {fastTime} Sqrt(2) = {FastSqrt(2)}");
    }
}

結果

 UnityEngine.Mathf.Sqrt time = 00:00:01.9403782 Sqrt(2) = 1.414214
 System.MathF.Sqrt time = 00:00:01.5182323 Sqrt(2) = 1.414214
 FastSqrt time = 00:00:02.3911146 Sqrt(2) = 1.41386

おそい！！！

はい、遅かったです。原因としては、 CPUに平方根演算命令が追加されたことが考えられます。というのも、このアルゴリズムが発表されたのが調べた限りでは2002年ごろであり、Intel製のCPUに平方根演算命令(SQRTSS)が追加されたのがおそらく2006年ごろ[2]になります。つまり、発表当時はCPUが専用命令を持たなかったために通常のSqrtの実装よりも高速に動作していたものだと考えられます。今のゲームができるような機器のCPUがこの命令をサポートしていないことは考えにくいため、この実装は避けて安直にMathクラスのSqrtを使用したほうがよさそうですね。

それではまた次回！

余談

かなりアバウトな精度でよい状況、かつ求める値の範囲が限定的な場合、線形補完で求めてしまう方法もあります。

/// <summary>
    /// 0~9までの値の平方根を線形補完する
    /// </summary>
    /// <param name="x"></param>
    /// <returns></returns>
    public static float AboutSqrt(float x)
    {
        if (x >= 0 && x <= 1)
        {
            return x;
        }
        else if (x > 1 && x <= 4)
        {
            return (x - 1) * 0.33333f + 1;
        }
        else if (x > 4 && x <= 9)
        {
            return (x - 4) * 0.2f + 2;
        }
        else
        {
            return Mathf.Sqrt(x);
        }
    }

このようなコードは他のいずれのSqrtよりも高速で動作します（精度は大きく劣りますが…）

UnityEngine.Mathf.Sqrt time = 00:00:01.7452324 Sqrt(2) = 1.414214
 System.MathF.Sqrt time = 00:00:01.3816707 Sqrt(2) = 1.414214
 AboutSqrt time = 00:00:00.8483539 Sqrt(2) = 1.33333
 FastSqrt time = 00:00:02.1512826 Sqrt(2) = 1.41386

精度にほとんど目をつぶれる状況で最速を求めるのならばこのような実装もありかもしれませんね。

参考文献
[1] David Eberly, Fast Inverse Square Root
[2] https://gcc.gnu.org/onlinedocs/gcc-3.4.6/gcc/X86-Built_002din-Functions.html