LoginSignup
0
1

More than 1 year has passed since last update.

@taqu

[0,1)の浮動小数点数乱数

Last updated at Posted at 2023-09-09

はじめに

[0.0, 1.0) の乱数を得るための“本当の”方法というスライドがタイムラインに流れてきたので, 少し興味を持ちました.

Downeyの方法では, 2^nの確率分布で指数を選択する必要がありますが, この部分を"Alias Method"で置き換えて高速化できないか試してみます.

浮動小数点数の乱数を得るために, 浮動小数点数の乱数を使うAlias Methodは本末転倒なので, 整数版のAlias Methodを作成して, パラメータテーブルを作ります.

次に事前計算したテーブルで2^nの確率分布で乱数を生成し, Downeyの方法を実装します.

設定は, 32bitの浮動小数点数を生成, 32bitの疑似乱数生成器を使用するとします.

整数Alias Method

まず, 2^nの確率分布を求めるためのAlias Methodのテーブルを作成します.
32bitの2^nの重みテーブルの総和は

 \sum_{i=0}^{31}x^i

ですので, 平均は

\frac{2^{31}+2^{31}-1}{2^5}

となります. 整数で演算するので,

\frac{2^{31}+2^{31}}{2^5}=2^{27}

に近似します.

Alias Methodのテーブル作成は次のようになります.

Alias Method整数版
void RandomAliasSelect::build(uint32_t size)
{
    assert(0 < size && size <= 32);
    size_ = size;

    std::vector<uint32_t> weights;
    weights.resize(size_);
    weights[0] = 1;
    for(uint32_t i = 1; i < size_; ++i) {
        weights[i] = weights[i - 1] << 1;
    }

    weights_.resize(size_);
    aliases_.resize(size_);

    std::vector<uint32_t> indices;
    indices.resize(size_);
    average_ = 1UL << (size_ - static_cast<uint32_t>(std::log2(size_)));
    int32_t underfull = -1;
    int32_t overfull = static_cast<int32_t>(size_);
    for(uint32_t i = 0; i < size_; ++i) {
        if(average_ <= weights[i]) {
            --overfull;
            indices[overfull] = i;
        } else {
            ++underfull;
            indices[underfull] = i;
        }
    }
    while(0 <= underfull && overfull < static_cast<int32_t>(size_)) {
        uint32_t under = indices[underfull];
        --underfull;
        uint32_t over = indices[overfull];
        ++overfull;
        aliases_[under] = over;
        weights_[under] = weights[under];
        weights[over] += weights[under] - average_;
        if(weights[over] < average_) {
            ++underfull;
            indices[underfull] = over;
        } else {
            --overfull;
            indices[overfull] = over;
        }
    }
    while(0 <= underfull) {
        weights_[indices[underfull]] = average_;
        --underfull;
    }
    while(overfull < static_cast<int32_t>(size_)) {
        weights_[indices[overfull]] = average_;
        ++overfull;
    }
}

求めたテーブルから乱数を生成する計算は次のようになります.

template<class T>
uint32_t RandomAliasSelect::select(T& random) const
{
    uint32_t index = random.rand() & (size_-1);
    uint32_t x = random.rand() & (average_-1);
    return x<weights_[index]? index : aliases_[index];
}

テーブルは次のようになります.

Alias Weight
0 28 1
1 29 2
2 29 4
3 29 8
4 30 16
5 30 32
6 30 64
7 30 128
8 30 256
9 30 512
10 30 1024
11 31 2048
12 31 4096
13 31 8192
14 31 16384
15 31 32768
16 31 65536
17 31 131072
18 31 262144
19 31 524288
20 31 1048576
21 31 2097152
22 31 4194304
23 31 8388608
24 31 16777216
25 31 33554432
26 31 67108864
27 0 134217728
28 27 134217727
29 28 134217726
30 29 134217712
31 30 134215680

テスト結果は正しいでしょうか？
pow2_distribution.png

Downeyの方法の事前計算テーブルによる実装

ビットスキャンを置き換えるだけなので, 次のようになります. ==1.0の場合に生成し直していますが, この処理の有無で性能に影響はありません.

事前計算テーブルによる実装
namespace
{
    uint8_t alias_pow2_32bit[] = {
        4,//28,
        3,//29,
        3,//29,
        3,//29,
        2,//30,
        2,//30,
        2,//30,
        2,//30,
        2,//30,
        2,//30,
        2,//30,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        1,//31,
        32,//0,
        5,//27,
        4,//28,
        3,//29,
        2,//30,
    };

    uint32_t weight_pow2_32bit[] = {
        1,
        2,
        4,
        8,
        16,
        32,
        64,
        128,
        256,
        512,
        1024,
        2048,
        4096,
        8192,
        16384,
        32768,
        65536,
        131072,
        262144,
        524288,
        1048576,
        2097152,
        4194304,
        8388608,
        16777216,
        33554432,
        67108864,
        134217728,
        134217727,
        134217726,
        134217712,
        134215680,
    };
    inline int32_t frandom_exponent(uint32_t x)
    {
        static constexpr uint32_t average_exp = 27;
        static constexpr uint32_t average_mask = (1UL << 27) - 1;
        uint32_t index = (x >> average_exp);
        uint32_t weight = x & average_mask;
        return weight < weight_pow2_32bit[index] ? 32-index : alias_pow2_32bit[index];
    }
} // namespace

float frandom_table(PCG32& random)
{
    int32_t exponent = 126;
    for(;;){
        uint32_t x;
        for(;;) {
            x = random();
            if(0 != x) /*[[likely]]*/ {
                exponent -= frandom_exponent(x);
                break;
            } else {
                exponent -= 32;
                if(exponent < 0) /*[[unlikely]]*/ {
                    exponent = 0;
                    break;
                }
            }
        }
        x = random();
        uint32_t fraction = x&0x7F'FFFFUL;
        if(0==fraction && (x&0x8000'0000UL)){
            if(126<=exponent)/*[[unlikely]]*/{
                continue;
            }
            ++exponent;
        }
        return std::bit_cast<float,uint32_t>((exponent<<23)|fraction);
    }
}

計測

スライドのコードを手直ししたコードを追加してパフォーマンス計測します. 事前計算テーブル版はfrandom_tableと呼称します.

frandom_downey
float frandom_downey(PCG32& random)
{
    constexpr int32_t lowExp = 0;
    constexpr int32_t highExp = 127;
    int32_t exponent = highExp - 1;
    while(true) {
        const uint32_t bits = random();
        if(0 == bits) {
            exponent -= 32;
            if(exponent < lowExp) {
                exponent = lowExp;
                break;
            }
        } else {
            int32_t c = std::countr_zero(bits);
            exponent -= c;
            break;
        }
    }
    const uint32_t u = random();
    const uint32_t mantissa = (u>>8)&0x7FFFFFUL;
    if(0==mantissa && (u>>31)){
        ++exponent;
    }
    return std::bit_cast<float,uint32_t>((exponent<<23)|mantissa);
}
frandom_downey_opt
float frandom_downey_opt(PCG32& random)
{
    constexpr int32_t lowExp = 0;
    constexpr int32_t highExp = 127;
    const uint32_t u = random();
    const uint32_t b = u & 0xFFU;
    int32_t exponent = highExp - 1;
    if(0 == b) {
        exponent -= 8;
        while(true) {
            const uint32_t bits = random();
            if(0 == bits) {
                exponent -= 32;
                if(exponent < lowExp) {
                    exponent = lowExp;
                    break;
                }
            } else {
                int32_t c = std::countr_zero(bits);
                exponent -= c;
                break;
            }
        }
    }else{
        int32_t c = std::countr_zero(b);
        exponent -= c;
    }
    const uint32_t mantissa = (u>>8)&0x7FFFFFUL;
    if(0==mantissa && (u>>31)){
        ++exponent;
    }
    return std::bit_cast<float,uint32_t>((exponent<<23)|mantissa);
}

環境は, 次のようになります.
CPU Core i5-10210U
OS Debian 12.1

100'000'000個の生成時間をテストしました.

ビルドオプションはCMakeLists.txtを参照してください.

CMakeLists.txt
cmake_minimum_required(VERSION 3.25)

set(CMAKE_CONFIGURATION_TYPES "Debug" "Release")

set(PROJECT_NAME frandom)
project(${PROJECT_NAME})

set(SOURCE_ROOT "${CMAKE_CURRENT_SOURCE_DIR}")

########################################################################
# Sources
set(HEADERS "${SOURCE_ROOT}/frandom.h")
set(SOURCES "${SOURCE_ROOT}/main.cpp;frandom.cpp")

source_group("include" FILES ${HEADERS})
source_group("src" FILES ${SOURCES})

set(FILES ${HEADERS} ${SOURCES})

set(OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/bin")
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG "${OUTPUT_DIRECTORY}")
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE "${OUTPUT_DIRECTORY}")

add_executable(${PROJECT_NAME} ${FILES})

if(MSVC)
    set(DEFAULT_CXX_FLAGS "/DWIN32 /D_WINDOWS /D_MSBC /W4 /WX- /nologo /fp:precise /arch:AVX2 /Zc:wchar_t /TP /Gd /std:c++20")
    if("1800" VERSION_LESS MSVC_VERSION)
        set(DEFAULT_CXX_FLAGS "${DEFAULT_CXX_FLAGS} /EHsc")
    endif()

    set(CMAKE_CXX_FLAGS "${DEFAULT_CXX_FLAGS}")
    set(CMAKE_CXX_FLAGS_DEBUG "/D_DEBUG /MDd /Zi /Ob0 /Od /RTC1 /Gy /GR- /GS /Gm-")
    set(CMAKE_CXX_FLAGS_RELEASE "/MD /O2 /GL /GR- /DNDEBUG")

elseif(UNIX)
    set(DEFAULT_CXX_FLAGS "-Wall -O2 -g -std=c++20 -std=gnu++20 -march=x86-64-v3")
    set(CMAKE_CXX_FLAGS "${DEFAULT_CXX_FLAGS}")
elseif(APPLE)
endif()

set_property(DIRECTORY PROPERTY VS_STARTUP_PROJECT ${PROJECT_NAME})
set_target_properties(${PROJECT_NAME}
    PROPERTIES
        OUTPUT_NAME_DEBUG "${PROJECT_NAME}" OUTPUT_NAME_RELEASE "${PROJECT_NAME}"
        VS_DEBUGGER_WORKING_DIRECTORY "${OUTPUT_DIRECTORY}")

g++

バージョンは12.2.0です.

timing (millioseconds)
frandom_table 536
frandom_downey 371
frandom_downey_opt 268

clang++

バージョンは14.0.6です.

timing (millioseconds)
frandom_table 521
frandom_downey 362
frandom_downey_opt 292

まとめ

テーブル参照とそのための演算より, CPUのビット演算の方が倍ぐらい速いです. CPU律速なので, コンパイラもOSもほとんど関係ないでしょう.

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
1