OSDN Git Service

RNGをXoroshiro128+に変更
authorHabu <habu@users.sourceforge.jp>
Tue, 25 Jul 2017 15:27:54 +0000 (00:27 +0900)
committerHabu <habu@users.sourceforge.jp>
Tue, 25 Jul 2017 15:27:54 +0000 (00:27 +0900)
Xoroshiro128+はXorshift系のRNGで、現在最も質が良く速いとされている。
通常は64bit型変数で実装するが、変愚蛮怒は歴史上32bit型変数までしか使用しないので32bitの変数2つを1組として扱い、64bitのXOR・シフト・ローテート・足し算を自前で実装した。

src/z-rand.c

index af939f6..8f9ce59 100644 (file)
@@ -70,6 +70,65 @@ u32b Rand_state[RAND_DEG] = {
 };
 
 
+typedef struct {
+       u32b dw[2];
+} u64b;
+
+static u64b u64b_xor(u64b a, u64b b)
+{
+       u64b result;
+
+       result.dw[0] = a.dw[0] ^ b.dw[0];
+       result.dw[1] = a.dw[1] ^ b.dw[1];
+
+       return result;
+}
+
+static u64b u64b_shiftl(u64b x, int k)
+{
+       u64b result;
+
+       if (k < 32) {
+               result.dw[1] = (x.dw[1] << k) | (x.dw[0] >> (32 - k));
+               result.dw[0] = (x.dw[0] << k);
+       }
+       else {
+               result.dw[1] = (x.dw[0] << (k - 32));
+               result.dw[0] = 0;
+       }
+
+       return result;
+}
+
+static u64b u64b_rotl(u64b x, int k)
+{
+       u64b result;
+
+       if (k < 32) {
+               result.dw[0] = (x.dw[0] << k) | (x.dw[1] >> (32 - k));
+               result.dw[1] = (x.dw[1] << k) | (x.dw[0] >> (32 - k));
+       }
+       else {
+               result.dw[0] = (x.dw[0] >> (64 - k)) | (x.dw[1] << (k - 32));
+               result.dw[1] = (x.dw[1] >> (64 - k)) | (x.dw[0] << (k - 32));
+       }
+
+       return result;
+}
+
+static u64b u64b_add(u64b a, u64b b)
+{
+       u64b result;
+
+       result.dw[0] = a.dw[0] + b.dw[0];
+       result.dw[1] = a.dw[1] + b.dw[1];
+
+       if (result.dw[0] < a.dw[0])
+               result.dw[1] ++;
+
+       return result;
+}
+
 /*
  * Initialize Xorshift Algorithm state
  */
@@ -84,6 +143,7 @@ static void Rand_Xorshift_seed(u32b seed, u32b* state)
        }
 }
 
+#if 0
 /*
  * Xorshift Algorithm
  */
@@ -99,6 +159,23 @@ static u32b Rand_Xorshift(u32b* state)
 
        return state[3];
 }
+#endif
+
+/*
+ * Xoroshiro128+ Algorithm
+ */
+static u32b Rand_Xoroshiro128plus(u32b* state)
+{
+       const u64b s0 = *((u64b*)state);
+       u64b s1 = *((u64b*)state + 1);
+       const u64b result = u64b_add(s0, s1);
+
+       s1 = u64b_xor(s0, s1);
+       *((u64b*)state) = u64b_xor(u64b_xor(u64b_rotl(s0, 55), s1), u64b_shiftl(s1, 14));
+       *((u64b*)state + 1) = u64b_rotl(s1, 36);
+
+       return result.dw[0];
+}
 
 static const u32b Rand_Xorshift_max = 0xFFFFFFFF;
 
@@ -189,7 +266,7 @@ static s32b Rand_div_impl(s32b m, u32b* state)
        past = scaling * m;
 
        do {
-               ret = Rand_Xorshift(state);
+               ret = Rand_Xoroshiro128plus(state);
        } while (ret >= past);
 
        return ret / scaling;