梅森旋转法产生随机数

直接贴代码了:

#include <stdio.h>
#include <stdlib.h>

typedef unsigned long uint32;

#define N           (624)           // length of state vector
#define M           (397)
#define K           (0x9908B0DFU)   // a magic constant
#define hiBit(u)    ((u) & 0x80000000U) // mask highest bit of u
#define loBit(u)    ((u) & 0x00000001U) // mask lowest bit of u
#define loBits(u)   ((u) & 0x7FFFFFFFU) // mask all but highest bit of u
#define mixBits(u, v)  (hiBit(u) | loBits(v)) // move high bit of u to low bits of v

void seedMT(uint32 seed);
uint32 reloadMT(void);
uint32 randomMT(void);

static uint32   state[N+1];
static uint32   *next;
static int      left = -1;

void seedMT(uint32 seed)
{
    register uint32 x = (seed | 1U) & 0xFFFFFFFFU, *s = state;
    register int    j;

    for(left = 0, *s++=x, j=N; --j;
        *s++ = (x*69069U) & 0xFFFFFFFFFF);
}

uint32 reloadMT(void)
{
    register uint32 *p0=state, *p2=state+2, *pM=state+M, s0, s1;
    register int    j;

    if(left < -1)
        seedMT(4357U);

    left=N-1, next=state+1;

    for(s0=state[0], s1=state[1], j=N-M+1; --j; s0=s1, s1=*p2++)
        *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);

    for(pM=state, j=M; --j; s0=s1, s1=*p2++)
        *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);

    s1 = state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
    s1 ^= (s1 >> 11);
    s1 ^= (s1 << 7) & 0x9D2C5680U;
    s1 ^= (s1 << 15) & 0xEFC60000U;
    return (s1 ^ (s1 >> 18));
}

uint32 randomMT(void)
{
    uint32 y;

    //printf("... left = %d\n", left);

    if(--left < 0)
        return (reloadMT());

    //printf("... *next = %d\n", *next);

    y = *next++;
    y ^= (y >> 11);
    y ^= (y << 7) & 0x9D2C5680U;
    y ^= (y << 15) & 0xEFC60000U;
    y ^= (y >> 18);
    return (y);
}

// test only
int main() {
    int seed;
    scanf("%d", &seed);
    seedMT(seed);
    int cnt = 1000000, a = 0, b = 0;
    while(cnt--) {
        if(randomMT()%2==0) {
            a++;
        }
        else b++;
        reloadMT();
    }
    printf("%d %d\n", a, b);
    return 0;
}


版权声明:本文为shuimu12345678原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接和本声明。