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

# はじめに

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

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

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

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

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

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

```math
 \sum_{i=0}^{31}x^i
```
ですので, 平均は
```math
\frac{2^{31}+2^{31}-1}{2^5}
```
となります. 整数で演算するので,

```math
\frac{2^{31}+2^{31}}{2^5}=2^{27}
```
に近似します.

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

<details><summary>Alias Method整数版</summary><div>

```c++
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;
    }
}
```
</div></details>

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

```c++
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/205dae1e-1f72-df36-7e44-17c9bf77d236.png)

# Downeyの方法の事前計算テーブルによる実装
ビットスキャンを置き換えるだけなので, 次のようになります. `==1.0`の場合に生成し直していますが, この処理の有無で性能に影響はありません.

<details><summary>事前計算テーブルによる実装</summary><div>

```c++
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);
    }
}
```
</div></details>

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

<details><summary>frandom_downey</summary><div>

```c++
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);
}
```
</div></details>


<details><summary>frandom_downey_opt</summary><div>

```c++
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);
}
```
</div></details>

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

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

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

<details><summary>CMakeLists.txt</summary><div>

```cmake
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}")
```
</div></details>

## 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もほとんど関係ないでしょう.