Skip to content

Instantly share code, notes, and snippets.

@gintenlabo

gintenlabo/evaluate_seed.cc

Created September 30, 2010 15:05
Show Gist options
  • Save gintenlabo/604721 to your computer and use it in GitHub Desktop.
Save gintenlabo/604721 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
evaluate_seed.cc
// オレオレ性能測定
#include "xor128.hpp"
// ローパスフィルタ
// 結果のランダムな性質を少し和らげて、より正確な処理ができるよう（気休めレベルですが）
struct LPF
{
typedef double result_type;
explicit LPF( double x0 = 0 )
: x_(x0), y_(x0) {}
result_type operator()( double x )
{
// カットオフ周波数に相当するもの
// 周期 Nc で繰り返されるものより波長の短い成分をカット
static int const Nc = 10;
// パラメータ設定
static double const T_omega = 1.0 / Nc;
static double const kx = T_omega / ( 2.0 + T_omega );
static double const ky = ( 2.0 - T_omega ) / ( 2.0 + T_omega );
// 生成
y_ = kx * ( x + x_ ) + ky * y_;
x_ = x;
return y_;
}
result_type get() const
{
return y_;
}
private:
double x_, y_;
};
// 乱数の初期状態を評価する
// 内部状態がほぼランダムになるまでに必要なステップ数を計算
// 今回は平均しないので直接的に 1 の bit 数を数える
std::size_t evaluate( etude::Xor128& gen )
{
LPF lpf;
lpf( gen.get_state().count() );
std::size_t const max_steps = 256;
for( std::size_t steps = 0; steps <= max_steps; ++steps )
{
// 判断基準。平滑化された 1 のビット数がこれ以上ならよし
// 2.0 は定常状態を LPF で平滑化したときの揺らぎの幅（標準偏差）の実測値 2.37 から
// それよりも適度に小さなキリのよい値として深く考えず採用。
static double const norm = 64 - 2.0;
if( lpf.get() >= norm )
{
return steps;
}
gen();
lpf( gen.get_state().count() );
}
return -1;
}
// 乱数生成器を評価する
// 内部状態のばらつきを示した「スコア」が基準値を上回るまでにかかったステップ数を計算する
#include <cmath>
template<class F>
double evaluate_seed( F f )
{
// 内部状態の bit 数
std::size_t const seed_bits = 32;
// std::time 等を考慮し、低い bit に対して重みを与える
double weight = 2.0;
double const rw = std::pow( weight, -1.0 / seed_bits );
double sum_score = 0;
double sum_weight = 0;
// 種を 0…01 から順にシフトしていく
for( etude::Xor128::int_type seed = 1; seed != 0; seed <<= 1 )
{
etude::Xor128 gen( f(seed) );
std::size_t score = evaluate(gen);
sum_score += weight * score;
sum_weight += weight;
weight *= rw;
}
// 平均
return sum_score / sum_weight;
}
// 32bit 整数を左回りに n bit 回転させる
inline std::uint32_t rotate_left_32( std::uint32_t x, std::size_t n ) {
return ( x << n ) | ( x >> ( 32 - n ) );
}
#include <iostream>
#include <boost/format.hpp>
#include <boost/progress.hpp>
#include <algorithm>
#include <vector>
#include <array>
#include <cfloat>
#include <cassert>
#include <boost/foreach.hpp>
// この辺でコメントによる解説は力尽きましたorz
int main()
{
typedef etude::Xor128::int_type int_type;
typedef std::array<std::size_t, 4> array_type;
std::vector< std::pair< double, array_type > > results;
std::size_t const results_max = 100;
// if を減らすため results にダミーの要素を追加
{
results.reserve( results_max + 1 );
array_type const dummy = {{ 32, 32, 32 }};
for( std::size_t i = 0; i < results_max; ++i ) {
results.emplace_back( DBL_MAX, dummy );
}
assert( results.size() == results_max );
std::make_heap( results.begin(), results.end() );
}
{
std::size_t const int_bits = 32;
std::size_t const tested_cases = 32 * 32 * 32 * 32 - 31 * 31 * 31 * 31;
boost::progress_display show_progress( tested_cases, std::cerr );
for( std::size_t a = 0; a < int_bits; ++a ) {
for( std::size_t b = a; b < int_bits; ++b ) {
for( std::size_t c = b; c < int_bits; ++c ) {
array_type nums = {{ 0, a, b, c }};
do {
std::size_t const dx = nums[0], dy = nums[1], dz = nums[2], dw = nums[3];
auto const f = [=]( int_type seed ) -> etude::Xor128 {
return etude::Xor128(
rotate_left_32( seed, dx ), rotate_left_32( seed, dy ),
rotate_left_32( seed, dz ), rotate_left_32( seed, dw )
);
};
double const score = evaluate_seed( f );
// 今回の結果をヒープに追加して
results.emplace_back( score, nums );
std::push_heap( results.begin(), results.end() );
// 最大の要素を取り除く
std::pop_heap( results.begin(), results.end() );
results.pop_back();
++show_progress;
}
while( std::next_permutation( nums.begin(), nums.end() ) );
}
}
}
}
std::cerr << std::endl;
// results の中身を sort
assert( results.size() == results_max );
std::sort_heap( results.begin(), results.end() );
// 表示
std::cout << "dx dy dz dw : score\n\n";
BOOST_FOREACH( auto const& x, results )
{
double const score = x.first;
auto const& params = x.second;
std::cout << boost::format("%1$02d %2$02d %3$02d %4$02d : %5%\n")
% params[0] % params[1] % params[2] % params[3] % score;
}
}
Raw
result-linear.out
dx dy dz dw : score
00 02 04 17 : 37.9233
30 00 02 15 : 38.029
28 30 00 13 : 38.2996
01 18 00 19 : 38.8293
00 17 31 18 : 38.8731
14 31 13 00 : 38.8782
15 00 14 01 : 38.9994
30 00 10 28 : 39.0322
00 02 04 03 : 39.0548
02 04 14 00 : 39.0838
16 12 26 00 : 39.0909
00 02 12 30 : 39.1302
00 28 10 16 : 39.2558
04 00 14 20 : 39.2764
07 27 00 25 : 39.3119
09 27 00 02 : 39.3389
12 00 05 30 : 39.3473
00 13 12 28 : 39.3563
00 18 23 25 : 39.3705
00 14 07 21 : 39.3905
04 17 16 00 : 39.4182
14 00 15 18 : 39.4227
07 25 30 00 : 39.433
20 01 00 16 : 39.4388
01 20 00 10 : 39.4633
28 30 00 07 : 39.4697
14 02 07 00 : 39.4852
28 30 00 31 : 39.4981
00 02 04 19 : 39.5
02 25 00 12 : 39.5058
29 31 01 00 : 39.511
02 04 06 00 : 39.5206
09 02 29 00 : 39.538
13 15 17 00 : 39.5722
22 18 00 06 : 39.5754
00 20 25 18 : 39.5921
12 05 00 03 : 39.5928
30 00 02 09 : 39.6012
30 00 02 01 : 39.6063
00 19 31 09 : 39.6095
00 23 30 10 : 39.6192
02 04 00 06 : 39.6456
00 02 04 11 : 39.6501
17 00 10 12 : 39.6514
00 14 01 18 : 39.6888
31 13 00 17 : 39.6894
09 29 28 00 : 39.6933
18 00 14 07 : 39.7004
00 02 25 20 : 39.712
04 18 00 25 : 39.7184
14 04 00 07 : 39.7249
20 22 00 18 : 39.7384
21 23 25 00 : 39.7416
19 00 31 15 : 39.7584
25 07 00 14 : 39.759
30 10 29 00 : 39.7738
30 00 23 18 : 39.7829
02 12 00 19 : 39.8022
00 02 30 04 : 39.8028
10 00 28 03 : 39.8151
17 00 15 21 : 39.8222
28 30 00 15 : 39.8299
07 29 25 00 : 39.8466
00 14 28 21 : 39.8473
00 07 12 17 : 39.8512
30 00 28 02 : 39.857
00 12 09 10 : 39.8582
02 27 00 09 : 39.8608
00 12 31 02 : 39.8615
00 02 30 20 : 39.864
14 00 05 07 : 39.8666
28 14 29 00 : 39.8763
18 00 25 07 : 39.8802
01 13 00 03 : 39.8847
22 15 29 00 : 39.8898
02 04 00 22 : 39.8911
30 00 02 17 : 39.8918
00 14 18 12 : 39.8943
02 17 00 06 : 39.9059
12 00 31 03 : 39.9098
07 25 00 13 : 39.9278
28 25 24 00 : 39.9336
07 00 27 30 : 39.9414
13 00 12 22 : 39.9555
20 02 06 00 : 39.9562
10 20 00 27 : 39.9588
13 01 00 04 : 39.9613
00 14 15 18 : 39.9646
22 00 18 17 : 39.9781
00 18 25 06 : 39.9813
18 00 19 04 : 39.9865
22 12 21 00 : 39.9929
18 00 14 04 : 39.9981
00 29 10 17 : 40.0064
00 02 04 05 : 40.011
13 15 00 06 : 40.0122
28 13 00 06 : 40.0142
18 25 06 00 : 40.0264
00 15 30 04 : 40.0271
07 22 13 00 : 40.0284
Raw
result.out
dx dy dz dw : score
00 02 04 17 : 37.8872
30 00 02 15 : 38.0415
28 30 00 13 : 38.3594
01 18 00 19 : 38.8072
00 17 31 18 : 38.8542
14 31 13 00 : 38.8915
15 00 14 01 : 39.0206
30 00 10 28 : 39.0264
00 02 04 03 : 39.0433
02 04 14 00 : 39.0517
16 12 26 00 : 39.066
00 02 12 30 : 39.1156
00 28 10 16 : 39.2661
04 00 14 20 : 39.2897
09 27 00 02 : 39.2938
00 14 07 21 : 39.3292
07 27 00 25 : 39.3388
00 13 12 28 : 39.3471
12 00 05 30 : 39.3501
00 18 23 25 : 39.3771
04 17 16 00 : 39.4018
07 25 30 00 : 39.4094
14 00 15 18 : 39.4109
20 01 00 16 : 39.4501
02 25 00 12 : 39.4522
28 30 00 07 : 39.46
01 20 00 10 : 39.4622
02 04 06 00 : 39.4718
14 02 07 00 : 39.4979
13 15 17 00 : 39.5152
28 30 00 31 : 39.5225
09 02 29 00 : 39.5274
12 05 00 03 : 39.5307
29 31 01 00 : 39.5339
00 02 04 19 : 39.5691
22 18 00 06 : 39.5851
30 00 02 09 : 39.5899
00 23 30 10 : 39.6059
00 20 25 18 : 39.6075
02 04 00 06 : 39.6152
17 00 10 12 : 39.6183
00 19 31 09 : 39.6237
30 00 02 01 : 39.6332
00 02 04 11 : 39.6408
00 14 01 18 : 39.6673
31 13 00 17 : 39.6745
14 04 00 07 : 39.6866
04 18 00 25 : 39.6885
09 29 28 00 : 39.6897
00 02 25 20 : 39.7006
18 00 14 07 : 39.7115
21 23 25 00 : 39.7767
19 00 31 15 : 39.7833
25 07 00 14 : 39.7867
30 10 29 00 : 39.7884
30 00 23 18 : 39.7976
00 02 30 04 : 39.802
17 00 15 21 : 39.8023
20 22 00 18 : 39.8067
14 00 05 07 : 39.8217
10 00 28 03 : 39.8239
02 12 00 19 : 39.824
00 07 12 17 : 39.8406
00 14 28 21 : 39.8469
00 12 09 10 : 39.8473
02 04 00 22 : 39.8542
00 02 30 20 : 39.8576
12 00 31 03 : 39.8635
02 27 00 09 : 39.8637
00 12 31 02 : 39.8676
00 14 18 12 : 39.8725
30 00 28 02 : 39.8738
07 29 25 00 : 39.8756
28 14 29 00 : 39.883
01 13 00 03 : 39.8866
28 30 00 15 : 39.8903
22 15 29 00 : 39.8948
07 25 00 13 : 39.895
02 17 00 06 : 39.9041
13 01 00 04 : 39.9043
28 25 24 00 : 39.9108
18 00 25 07 : 39.9286
13 00 12 22 : 39.9393
10 20 00 27 : 39.9489
22 12 21 00 : 39.9649
22 00 18 17 : 39.9672
00 14 15 18 : 39.9684
28 13 00 06 : 39.971
30 00 02 17 : 39.9731
07 00 27 30 : 39.9765
13 15 00 06 : 39.9833
00 29 10 17 : 39.9889
18 17 00 22 : 40.0021
00 18 25 06 : 40.005
18 00 19 04 : 40.0059
00 18 17 01 : 40.0098
18 00 14 04 : 40.011
03 00 13 20 : 40.0116
18 25 06 00 : 40.0195
10 00 20 14 : 40.0254
Raw
rotate_left.hpp
#ifndef INCLUDED_ROTATE_LEFT_HPP_
#define INCLUDED_ROTATE_LEFT_HPP_
// ビット列を左に N ビットシフトして返す
#include <limits> // for std::numeric_limits
#include <boost/integer/integer_mask.hpp>
// ヘッダに分割したので一応名前空間に入れておこう
namespace etude
{
// とりあえずテンプレート版。関数版も用意するかも
template<std::size_t N, class Int>
inline Int rotate_left( Int x )
{
typedef std::numeric_limits<Int> limits;
static const std::size_t Int_bits = limits::digits + ( limits::is_signed ? 1 : 0 );
static const Int mask = static_cast<Int>( boost::low_bits_mask_t<N>::sig_bits );
return ( x << N ) | ( mask & ( x >> ( Int_bits - N ) ) );
}
}
#endif // #ifndef INCLUDED_ROTATE_LEFT_HPP_
Raw
xor128.hpp
#ifndef INCLUDED_XOR128_HPP_
#define INCLUDED_XOR128_HPP_
// xor128 乱数生成アルゴリズム
#include <cassert>
#include <boost/cstdint.hpp>
#include <bitset>
#include "rotate_left.hpp"
// ヘッダに分割したので一応名前空間に入れておこう
namespace etude
{
struct Xor128
{
typedef boost::uint32_t result_type;
typedef boost::uint32_t int_type;
// min, max
static result_type min() { return 0; }
static result_type max() { return 0xFFFFFFFF; }
// デフォルト初期化
Xor128()
: x(123456789), y(362436069), z(521288629), w(88675123) {}
// 種を元に乱数を初期化する
explicit Xor128( int_type seed )
: x(123456789), y(362436069), z(521288629), w(88675123)
{
x ^= seed;
y ^= rotate_left<17>(seed);
z ^= rotate_left<31>(seed);
w ^= rotate_left<18>(seed);
// あるいは
// x ^= rotate_left<14>(seed);
// y ^= rotate_left<31>(seed);
// z ^= rotate_left<13>(seed);
// w ^= seed;
}
// 全ての内部状態を指定して乱数を構築する
Xor128( int_type x_, int_type y_, int_type z_, int_type w_ )
: x(x_), y(y_), z(z_), w(w_)
{
assert( x != 0 || y != 0 || z != 0 || w != 0 );
}
// 乱数生成
result_type operator()()
{
int_type const t = x ^ ( x << 11 );
x = y; y = z; z = w;
w = ( w ^ ( w >> 19 ) ) ^ ( t ^ ( t >> 8 ) );
return w;
}
// n 回読み飛ばす
void discard( std::size_t n )
{
for( std::size_t i = 0; i < n; ++i ){ (*this)(); }
}
// 比較
friend bool operator==( Xor128 const& lhs, Xor128 const& rhs )
{
return lhs.x == rhs.x
&& lhs.y == rhs.y
&& lhs.z == rhs.z
&& lhs.w == rhs.w;
}
friend bool operator!=( Xor128 const& lhs, Xor128 const& rhs )
{
return !( lhs == rhs );
}
// 内部状態を取得する
typedef std::bitset<128> state_type;
state_type get_state() const
{
state_type state = x;
state <<= 32;
state |= state_type(y);
state <<= 32;
state |= state_type(z);
state <<= 32;
state |= state_type(w);
return state;
}
// state_type からの構築
Xor128( state_type state )
{
state_type const mask = int_type(-1);
w = ( mask & state ).to_ulong();
state >>= 32;
z = ( mask & state ).to_ulong();
state >>= 32;
y = ( mask & state ).to_ulong();
state >>= 32;
x = ( mask & state ).to_ulong();
}
private:
int_type x, y, z, w;
};
}
#endif // #ifndef INCLUDED_XOR128_HPP_
@gintenlabo
Copy link
Author

俺々性能チェック。説明面倒なので解読してください。

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment