GLSLで使われる乱数の品質を確認する

GLSLで使われる適当な乱数の品質について疑問を持つことがあったので，ちょっと調べてみた。 今回の簡易評価の結論としては，それなりに一様という印象を受けた。

評価方法

とりあえず100x100の10,000回，0から1の乱数を生成し，その値をソートし，その結果の散布図にする。 これは乱数が一様に近い性質があれば直線に見えることを期待したものである。 （この方法では周期性などは確認できないのだけれど）

なお，GLSLから数値を得るのはだいぶ面倒なので，評価にはpythonを使った。

評価対象

まずGLSLでよく使われる乱数。結果の凡例では仮にfract-sin-dotと表記する。

float rand(vec2 p){
    return fract(sin(dot(p, vec2(12.9898,78.233))) * 43758.5453);
}

あとは32bitのXorshiftと，MT(Mersenne twister)。MTはPythonのrandomを使っている（比較対象として）。

評価結果

どの乱数の散布図も直線状に分布しているので，だいたい一様なんだなーということは確認できた。

評価に使ったコード

fractなど，GLSL組み込み関数は適当に実装した。

import random
import matplotlib.pyplot as plt
import numpy as np
import math


def fract(x):
    return x - math.floor(x)


def rand(p):
    return fract(math.sin(np.dot(p, np.array([12.9898, 78.233]))) * 43758.5453)


def main():
    N = 1

    sqrt_N = 100
    runs = 0

    seq1 = []
    seq2 = []
    seq3 = []

    xor_seed = 0x10A1F0F0

    def xorshift(seed):
        seed ^= seed << 13
        seed ^= seed >> 17
        seed ^= seed << 5
        return seed & 0xFFFFFFFF

    for k in range(N):
        for i in range(sqrt_N):
            for j in range(sqrt_N):
                runs += 1

                p = np.array([i+k*0.01, j+k*0.01])
                r1 = rand(p)

                xor_seed = xorshift(xor_seed)
                r2 = xor_seed / 0xFFFFFFFF

                r3 = random.random()

                seq1.append(r1)
                seq2.append(r2)
                seq3.append(r3)

    seq1.sort()
    seq2.sort()
    seq3.sort()

    axis = [range(runs)]

    plt.scatter(axis, seq1, s=1, label='GLSL-rand[f:id:ma38su:20210118123530p:plain]')
    plt.scatter(axis, seq2, s=1, label='Xorshift')
    plt.scatter(axis, seq3, s=1, label='MT')
    plt.legend()
    plt.grid(True)
    plt.savefig('fig.png')


main()

試してないけどGold Noiseというのもあるらしい。

// Gold Noise ©2015 dcerisano@standard3d.com 
//  - based on the Golden Ratio, PI and Square Root of Two
//  - superior distribution
//  - fastest noise generator function
//  - works with all chipsets (including low precision)

float PHI = 1.61803398874989484820459 * 00000.1; // Golden Ratio   
float PI  = 3.14159265358979323846264 * 00000.1; // PI
float SQ2 = 1.41421356237309504880169 * 10000.0; // Square Root of Two

float gold_noise(in vec2 coordinate, in float seed){
    return fract(tan(distance(coordinate*(seed+PHI), vec2(PHI, PI)))*SQ2);
}