Fast sampler of exponential distribution （高速指数分布サンプラー・擬似焼きなまし法の遷移判定）
(heuristic/exponential_dist_sampler.hpp)

• View this file on GitHub
• Last update: 2026-08-28 23:19:00+09:00
• Include: #include "heuristic/exponential_dist_sampler.hpp"

母数 $\lambda = 1$ の指数分布を模擬するサンプラーと，これを利用した擬似焼きなまし法 (simulated annealing, SA) の温度・遷移判定を管理するクラス．

ExponentialDistSampler は $2^D$ 個の $-\log p$ を構築時に計算する．遷移判定時はテーブルを参照するだけなので，指数関数や対数関数を毎回計算する必要がない．

ExponentialDistSampler

constexpr int D = 16;
const ExponentialDistSampler<D> eds;  // 2^D 個のサンプルを前計算
FastRNG rng(123456789);

// 下位 D bit がランダムに分布した mask を与えると x ~ Ex(1) をサンプル
uint32_t mask= rng.next_u16();
double x = eds.sample(mask);

double abs_dx;
double temperature;

// 悪化量 abs_dx、温度 temperature として遷移を受理するか判定
bool accept = eds.check_sa(abs_dx, temperature, mask);

sample(mask) は mask の下位 $D$ bit を利用して，指数分布に従う非負の値を返す．check_sa(abs_dx, T, mask) は確率 $\exp(-\lvert dx \rvert / T)$ で真を返す．mask の下位 $D$ bit は一様に分布している必要がある．

Annealer

温度を等比的に変化させ，スコア最大化またはコスト最小化の遷移を判定する．以下は FastRNG::next_u16() と組み合わせる例である．

constexpr int ITERATIONS = 1000000;
constexpr int UPDATE_INTERVAL = 256;

FastRNG rng(123456789);
const ExponentialDistSampler<16> sampler;
Annealer<UPDATE_INTERVAL> annealer(1000.0, 1.0, ITERATIONS);

int score = initial_score;
for (int iter = 0; iter < ITERATIONS; ++iter) {
    annealer.update(iter);

    int delta_score = candidate_score - score;
    if (annealer.accept_score(delta_score, rng, sampler)) {
        score = candidate_score;
        // 状態を candidate に更新
    }
}

• update(iter)：UPDATE_INTERVAL 回ごとに温度を更新する．各反復で呼び出す．
• accept_score(delta_score, rng, sampler)：スコア最大化用．delta_score >= 0 なら必ず受理し，負なら温度に応じた確率で受理する．
• accept_cost(delta_cost, rng, sampler)：コスト最小化用．delta_cost <= 0 なら必ず受理し，正なら温度に応じた確率で受理する．

start_temp，end_temp，UPDATE_INTERVAL，iterations には正の値を与えること．Annealer は乱数生成器の next_u16() を利用するため，組み合わせる ExponentialDistSampler では $D \leq 16$ とする．

Code

#pragma once

#include <cassert>
#include <cmath>
#include <cstdint>
#include <array>

template <int D> struct ExponentialDistSampler {
    std::array<double, (1 << D)> minuslogps;

    constexpr ExponentialDistSampler() {
        for (int i = 0; i < (1 << D); ++i) minuslogps.at(i) = -log((0.5 + i) / (1 << D));
    }

    double sample(uint32_t random_mask) const {
        return minuslogps.at(random_mask & ((1 << D) - 1));
    }

    // p ~ U(0, 1) => -log(p) ~ Ex(1)
    // P[exp(-|dx| / T) >= p] = P[|dx| <= -log(p) * T]
    bool check_sa(double abs_dx, double T, uint32_t random_mask) const {
        return abs_dx <= minuslogps.at(random_mask & ((1 << D) - 1)) * T;
    }
};
const ExponentialDistSampler<16> log_ps;

template <int UPDATE_INTERVAL = 256> struct Annealer {
    static_assert(UPDATE_INTERVAL > 0);

    double temp;
    double ratio;
    int next_update = UPDATE_INTERVAL;

    Annealer(double start_temp, double end_temp, int iterations) : temp(start_temp) {
        assert(iterations > 0);
        const int updates = (iterations - 1) / UPDATE_INTERVAL;
        ratio = updates == 0 ? 1.0 : pow(end_temp / start_temp, 1.0 / updates);
        // ratio = updates == 0 ? 0.0 : (end_temp - start_temp) / updates;
    }

    void update(int iter) {
        if (iter >= next_update) {
            temp *= ratio;
            // temp += ratio;
            next_update += UPDATE_INTERVAL;
        }
    }

    template <class RNG, class Sampler>
    bool accept_score(int delta_score, RNG &rng, const Sampler &sampler) const {
        if (delta_score >= 0) return true;
        return sampler.check_sa(-delta_score, temp, rng.next_u16());
    }

    template <class RNG, class Sampler>
    bool accept_cost(int delta_cost, RNG &rng, const Sampler &sampler) const {
        if (delta_cost <= 0) return true;
        return sampler.check_sa(delta_cost, temp, rng.next_u16());
    }
};

#line 2 "heuristic/exponential_dist_sampler.hpp"

#include <cassert>
#include <cmath>
#include <cstdint>
#include <array>

template <int D> struct ExponentialDistSampler {
    std::array<double, (1 << D)> minuslogps;

    constexpr ExponentialDistSampler() {
        for (int i = 0; i < (1 << D); ++i) minuslogps.at(i) = -log((0.5 + i) / (1 << D));
    }

    double sample(uint32_t random_mask) const {
        return minuslogps.at(random_mask & ((1 << D) - 1));
    }

    // p ~ U(0, 1) => -log(p) ~ Ex(1)
    // P[exp(-|dx| / T) >= p] = P[|dx| <= -log(p) * T]
    bool check_sa(double abs_dx, double T, uint32_t random_mask) const {
        return abs_dx <= minuslogps.at(random_mask & ((1 << D) - 1)) * T;
    }
};
const ExponentialDistSampler<16> log_ps;

template <int UPDATE_INTERVAL = 256> struct Annealer {
    static_assert(UPDATE_INTERVAL > 0);

    double temp;
    double ratio;
    int next_update = UPDATE_INTERVAL;

    Annealer(double start_temp, double end_temp, int iterations) : temp(start_temp) {
        assert(iterations > 0);
        const int updates = (iterations - 1) / UPDATE_INTERVAL;
        ratio = updates == 0 ? 1.0 : pow(end_temp / start_temp, 1.0 / updates);
        // ratio = updates == 0 ? 0.0 : (end_temp - start_temp) / updates;
    }

    void update(int iter) {
        if (iter >= next_update) {
            temp *= ratio;
            // temp += ratio;
            next_update += UPDATE_INTERVAL;
        }
    }

    template <class RNG, class Sampler>
    bool accept_score(int delta_score, RNG &rng, const Sampler &sampler) const {
        if (delta_score >= 0) return true;
        return sampler.check_sa(-delta_score, temp, rng.next_u16());
    }

    template <class RNG, class Sampler>
    bool accept_cost(int delta_cost, RNG &rng, const Sampler &sampler) const {
        if (delta_cost <= 0) return true;
        return sampler.check_sa(delta_cost, temp, rng.next_u16());
    }
};

Back to top page