CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual users to large groups and classes!
CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual users to large groups and classes!
Path: blob/master/ext/sfmt19937/SFMT.h
Views: 1401
#pragma once1/**2* @file SFMT.h3*4* @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom5* number generator using C structure.6*7* @author Mutsuo Saito (Hiroshima University)8* @author Makoto Matsumoto (The University of Tokyo)9*10* Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima11* University.12* Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima13* University and The University of Tokyo.14* All rights reserved.15*16* The 3-clause BSD License is applied to this software, see17* LICENSE.txt18*19* @note We assume that your system has inttypes.h. If your system20* doesn't have inttypes.h, you have to typedef uint32_t and uint64_t,21* and you have to define PRIu64 and PRIx64 in this file as follows:22* @verbatim23typedef unsigned int uint32_t24typedef unsigned long long uint64_t25#define PRIu64 "llu"26#define PRIx64 "llx"27@endverbatim28* uint32_t must be exactly 32-bit unsigned integer type (no more, no29* less), and uint64_t must be exactly 64-bit unsigned integer type.30* PRIu64 and PRIx64 are used for printf function to print 64-bit31* unsigned int and 64-bit unsigned int in hexadecimal format.32*/3334#ifndef SFMTST_H35#define SFMTST_H36#if defined(__cplusplus)37extern "C" {38#endif3940#include <stdio.h>41#include <assert.h>4243#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)44#include <inttypes.h>45#elif defined(_MSC_VER) || defined(__BORLANDC__)46typedef unsigned int uint32_t;47typedef unsigned __int64 uint64_t;48#if !defined(__cplusplus)49#define inline __inline50#endif51#else52#include <inttypes.h>53#if defined(__GNUC__)54#define inline __inline__55#endif56#endif5758#ifndef PRIu6459#if defined(_MSC_VER) || defined(__BORLANDC__)60#define PRIu64 "I64u"61#define PRIx64 "I64x"62#else63#define PRIu64 "llu"64#define PRIx64 "llx"65#endif66#endif6768#include "SFMT-params.h"6970/*------------------------------------------71128-bit SIMD like data type for standard C72------------------------------------------*/73#if defined(HAVE_ALTIVEC)74#if !defined(__APPLE__)75#include <altivec.h>76#endif77/** 128-bit data structure */78union W128_T {79vector unsigned int s;80uint32_t u[4];81uint64_t u64[2];82};83#elif defined(HAVE_SSE2)84#include <emmintrin.h>8586/** 128-bit data structure */87union W128_T {88uint32_t u[4];89uint64_t u64[2];90__m128i si;91};92#else93/** 128-bit data structure */94union W128_T {95uint32_t u[4];96uint64_t u64[2];97};98#endif99100/** 128-bit data type */101typedef union W128_T w128_t;102103/**104* SFMT internal state105*/106struct SFMT_T {107/** the 128-bit internal state array */108w128_t state[SFMT_N];109/** index counter to the 32-bit internal state array */110int idx;111};112113typedef struct SFMT_T sfmt_t;114115void sfmt_fill_array32(sfmt_t * sfmt, uint32_t * array, int size);116void sfmt_fill_array64(sfmt_t * sfmt, uint64_t * array, int size);117void sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed);118void sfmt_init_by_array(sfmt_t * sfmt, uint32_t * init_key, int key_length);119const char * sfmt_get_idstring(sfmt_t * sfmt);120int sfmt_get_min_array_size32(sfmt_t * sfmt);121int sfmt_get_min_array_size64(sfmt_t * sfmt);122void sfmt_gen_rand_all(sfmt_t * sfmt);123124#ifndef ONLY64125/**126* This function generates and returns 32-bit pseudorandom number.127* init_gen_rand or init_by_array must be called before this function.128* @param sfmt SFMT internal state129* @return 32-bit pseudorandom number130*/131inline static uint32_t sfmt_genrand_uint32(sfmt_t * sfmt) {132uint32_t r;133uint32_t * psfmt32 = &sfmt->state[0].u[0];134135if (sfmt->idx >= SFMT_N32) {136sfmt_gen_rand_all(sfmt);137sfmt->idx = 0;138}139r = psfmt32[sfmt->idx++];140return r;141}142#endif143/**144* This function generates and returns 64-bit pseudorandom number.145* init_gen_rand or init_by_array must be called before this function.146* The function gen_rand64 should not be called after gen_rand32,147* unless an initialization is again executed.148* @param sfmt SFMT internal state149* @return 64-bit pseudorandom number150*/151inline static uint64_t sfmt_genrand_uint64(sfmt_t * sfmt) {152#if defined(BIG_ENDIAN64) && !defined(ONLY64)153uint32_t * psfmt32 = &sfmt->state[0].u[0];154uint32_t r1, r2;155#else156uint64_t r;157#endif158uint64_t * psfmt64 = &sfmt->state[0].u64[0];159assert(sfmt->idx % 2 == 0);160161if (sfmt->idx >= SFMT_N32) {162sfmt_gen_rand_all(sfmt);163sfmt->idx = 0;164}165#if defined(BIG_ENDIAN64) && !defined(ONLY64)166r1 = psfmt32[sfmt->idx];167r2 = psfmt32[sfmt->idx + 1];168sfmt->idx += 2;169return ((uint64_t)r2 << 32) | r1;170#else171r = psfmt64[sfmt->idx / 2];172sfmt->idx += 2;173return r;174#endif175}176177/* =================================================178The following real versions are due to Isaku Wada179================================================= */180/**181* converts an unsigned 32-bit number to a double on [0,1]-real-interval.182* @param v 32-bit unsigned integer183* @return double on [0,1]-real-interval184*/185inline static double sfmt_to_real1(uint32_t v)186{187return v * (1.0/4294967295.0);188/* divided by 2^32-1 */189}190191/**192* generates a random number on [0,1]-real-interval193* @param sfmt SFMT internal state194* @return double on [0,1]-real-interval195*/196inline static double sfmt_genrand_real1(sfmt_t * sfmt)197{198return sfmt_to_real1(sfmt_genrand_uint32(sfmt));199}200201/**202* converts an unsigned 32-bit integer to a double on [0,1)-real-interval.203* @param v 32-bit unsigned integer204* @return double on [0,1)-real-interval205*/206inline static double sfmt_to_real2(uint32_t v)207{208return v * (1.0/4294967296.0);209/* divided by 2^32 */210}211212/**213* generates a random number on [0,1)-real-interval214* @param sfmt SFMT internal state215* @return double on [0,1)-real-interval216*/217inline static double sfmt_genrand_real2(sfmt_t * sfmt)218{219return sfmt_to_real2(sfmt_genrand_uint32(sfmt));220}221222/**223* converts an unsigned 32-bit integer to a double on (0,1)-real-interval.224* @param v 32-bit unsigned integer225* @return double on (0,1)-real-interval226*/227inline static double sfmt_to_real3(uint32_t v)228{229return (((double)v) + 0.5)*(1.0/4294967296.0);230/* divided by 2^32 */231}232233/**234* generates a random number on (0,1)-real-interval235* @param sfmt SFMT internal state236* @return double on (0,1)-real-interval237*/238inline static double sfmt_genrand_real3(sfmt_t * sfmt)239{240return sfmt_to_real3(sfmt_genrand_uint32(sfmt));241}242243/**244* converts an unsigned 32-bit integer to double on [0,1)245* with 53-bit resolution.246* @param v 32-bit unsigned integer247* @return double on [0,1)-real-interval with 53-bit resolution.248*/249inline static double sfmt_to_res53(uint64_t v)250{251return v * (1.0/18446744073709551616.0);252}253254/**255* generates a random number on [0,1) with 53-bit resolution256* @param sfmt SFMT internal state257* @return double on [0,1) with 53-bit resolution258*/259inline static double sfmt_genrand_res53(sfmt_t * sfmt)260{261return sfmt_to_res53(sfmt_genrand_uint64(sfmt));262}263264265/* =================================================266The following function are added by Saito.267================================================= */268/**269* generates a random number on [0,1) with 53-bit resolution from two270* 32 bit integers271*/272inline static double sfmt_to_res53_mix(uint32_t x, uint32_t y)273{274return sfmt_to_res53(x | ((uint64_t)y << 32));275}276277/**278* generates a random number on [0,1) with 53-bit resolution279* using two 32bit integers.280* @param sfmt SFMT internal state281* @return double on [0,1) with 53-bit resolution282*/283inline static double sfmt_genrand_res53_mix(sfmt_t * sfmt)284{285uint32_t x, y;286287x = sfmt_genrand_uint32(sfmt);288y = sfmt_genrand_uint32(sfmt);289return sfmt_to_res53_mix(x, y);290}291292#if defined(__cplusplus)293}294#endif295296#endif297298299