00001 #include "assert.h"
00002
00003 #define MSK1 0xdfffffefU
00004 #define MSK2 0xddfecb7fU
00005 #define MSK3 0xbffaffffU
00006 #define MSK4 0xbffffff6U
00007 #define PARITY1 0x00000001U
00008 #define PARITY2 0x00000000U
00009 #define PARITY3 0x00000000U
00010 #define PARITY4 0x13c9e684U
00011
00012 class mdp_prng_sfmt {
00013 private:
00014 bool initialized;
00015 static const unsigned int MEXP=19937;
00016 static const unsigned int N=156;
00017 static const unsigned int N32=624;
00018 static const unsigned int POS1=122;
00019 static const unsigned int SL1=18;
00020 static const unsigned int SL2=1;
00021 static const unsigned int SR1=11;
00022 static const unsigned int SR2=1;
00023 unsigned int idx;
00024 struct W128_T { unsigned int u[4]; };
00025 typedef struct W128_T w128_t;
00026 w128_t sfmt[N];
00027 unsigned int *psfmt32;
00028 void period_certification(void) {
00029 static unsigned int parity[4] = {PARITY1, PARITY2, PARITY3, PARITY4};
00030 unsigned int inner = 0;
00031 unsigned int i, j;
00032 unsigned int work;
00033
00034 for (i = 0; i < 4; i++)
00035 inner ^= psfmt32[i] & parity[i];
00036 for (i = 16; i > 0; i >>= 1)
00037 inner ^= inner >> i;
00038 inner &= 1;
00039
00040 if (inner == 1) { return; }
00041
00042 for (i = 0; i < 4; i++) {
00043 work = 1;
00044 for (j = 0; j < 32; j++) {
00045 if ((work & parity[i]) != 0) {
00046 psfmt32[i] ^= work;
00047 return;
00048 }
00049 work = work << 1;
00050 }
00051 }
00052 }
00053 void rshift128(w128_t *out, w128_t const *in, int shift) {
00054 unsigned int mdp_int th, tl, oh, ol;
00055
00056 th = ((unsigned int mdp_int)in->u[3] << 32) | ((unsigned int mdp_int)in->u[2]);
00057 tl = ((unsigned int mdp_int)in->u[1] << 32) | ((unsigned int mdp_int)in->u[0]);
00058
00059 oh = th >> (shift * 8);
00060 ol = tl >> (shift * 8);
00061 ol |= th << (64 - shift * 8);
00062 out->u[1] = (unsigned int)(ol >> 32);
00063 out->u[0] = (unsigned int)ol;
00064 out->u[3] = (unsigned int)(oh >> 32);
00065 out->u[2] = (unsigned int)oh;
00066 }
00067 void lshift128(w128_t *out, w128_t const *in, int shift) {
00068 unsigned int mdp_int th, tl, oh, ol;
00069
00070 th = ((unsigned int mdp_int)in->u[3] << 32) | ((unsigned int mdp_int)in->u[2]);
00071 tl = ((unsigned int mdp_int)in->u[1] << 32) | ((unsigned int mdp_int)in->u[0]);
00072
00073 oh = th << (shift * 8);
00074 ol = tl << (shift * 8);
00075 oh |= tl >> (64 - shift * 8);
00076 out->u[1] = (unsigned int)(ol >> 32);
00077 out->u[0] = (unsigned int)ol;
00078 out->u[3] = (unsigned int)(oh >> 32);
00079 out->u[2] = (unsigned int)oh;
00080 }
00081
00082 void gen_rand_all(void) {
00083 int i;
00084 w128_t *r1, *r2;
00085
00086 r1 = &sfmt[N - 2];
00087 r2 = &sfmt[N - 1];
00088 for (i = 0; i < N - POS1; i++) {
00089 do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1], r1, r2);
00090 r1 = r2;
00091 r2 = &sfmt[i];
00092 }
00093 for (; i < N; i++) {
00094 do_recursion(&sfmt[i], &sfmt[i], &sfmt[i + POS1 - N], r1, r2);
00095 r1 = r2;
00096 r2 = &sfmt[i];
00097 }
00098 }
00099
00100 void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
00101 w128_t *d) {
00102 w128_t x;
00103 w128_t y;
00104
00105 lshift128(&x, a, SL2);
00106 rshift128(&y, c, SR2);
00107 r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SR1) & MSK1) ^ y.u[0]
00108 ^ (d->u[0] << SL1);
00109 r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SR1) & MSK2) ^ y.u[1]
00110 ^ (d->u[1] << SL1);
00111 r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SR1) & MSK3) ^ y.u[2]
00112 ^ (d->u[2] << SL1);
00113 r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SR1) & MSK4) ^ y.u[3]
00114 ^ (d->u[3] << SL1);
00115 }
00116 unsigned int gen_rand32() {
00117 unsigned int r;
00118 assert(initialized);
00119 if (idx >= N32) {
00120 gen_rand_all();
00121 idx = 0;
00122 }
00123 r = psfmt32[idx++];
00124 return r;
00125 }
00126
00127 public:
00128 void initialize(unsigned int seed) {
00129 unsigned int i;
00130 psfmt32=(unsigned int *)&(sfmt[0].u[0]);
00131 psfmt32[0] = seed;
00132 for (i = 1; i < N32; i++) {
00133 psfmt32[i] = 1812433253UL * (psfmt32[i - 1]
00134 ^ (psfmt32[i - 1] >> 30))
00135 + i;
00136 }
00137 idx = N32;
00138 period_certification();
00139 initialized = 1;
00140 }
00141
00142
00143 float plain() {
00144 unsigned int v=gen_rand32();
00145 return v * (1.0/4294967295.0);
00146 }
00147 };
00148
