"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/rng/SFMT.h" (29 Aug 2019, 7717 Bytes) of package /linux/misc/gretl-2020e.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "SFMT.h" see the Fossies "Dox" file reference documentation.

    1 #pragma once
    2 /**
    3  * @file SFMT.h
    4  *
    5  * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
    6  * number generator using C structure.
    7  *
    8  * @author Mutsuo Saito (Hiroshima University)
    9  * @author Makoto Matsumoto (The University of Tokyo)
   10  *
   11  * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
   12  * University.
   13  * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
   14  * University and The University of Tokyo.
   15  * All rights reserved.
   16  *
   17  * The 3-clause BSD License is applied to this software, see
   18  * LICENSE.txt
   19  *
   20  * @note We assume that your system has inttypes.h.  If your system
   21  * doesn't have inttypes.h, you have to typedef uint32_t and uint64_t,
   22  * and you have to define PRIu64 and PRIx64 in this file as follows:
   23  * @verbatim
   24  typedef unsigned int uint32_t
   25  typedef unsigned long long uint64_t
   26  #define PRIu64 "llu"
   27  #define PRIx64 "llx"
   28 @endverbatim
   29  * uint32_t must be exactly 32-bit unsigned integer type (no more, no
   30  * less), and uint64_t must be exactly 64-bit unsigned integer type.
   31  * PRIu64 and PRIx64 are used for printf function to print 64-bit
   32  * unsigned int and 64-bit unsigned int in hexadecimal format.
   33  */
   34 
   35 #ifndef SFMTST_H
   36 #define SFMTST_H
   37 #if defined(__cplusplus)
   38 extern "C" {
   39 #endif
   40 
   41 #include <stdio.h>
   42 #include <assert.h>
   43 
   44 #if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
   45   #include <inttypes.h>
   46 #elif defined(_MSC_VER) || defined(__BORLANDC__)
   47   typedef unsigned int uint32_t;
   48   typedef unsigned __int64 uint64_t;
   49   #define inline __inline
   50 #else
   51   #include <inttypes.h>
   52   #if defined(__GNUC__)
   53     #define inline __inline__
   54   #endif
   55 #endif
   56 
   57 #ifndef PRIu64
   58   #if defined(_MSC_VER) || defined(__BORLANDC__)
   59     #define PRIu64 "I64u"
   60     #define PRIx64 "I64x"
   61   #else
   62     #define PRIu64 "llu"
   63     #define PRIx64 "llx"
   64   #endif
   65 #endif
   66 
   67 #include "SFMT-params.h"
   68 
   69 /*------------------------------------------
   70   128-bit SIMD like data type for standard C
   71   ------------------------------------------*/
   72 #if defined(HAVE_ALTIVEC)
   73   #if !defined(__APPLE__)
   74     #include <altivec.h>
   75   #endif
   76 /** 128-bit data structure */
   77 union W128_T {
   78     vector unsigned int s;
   79     uint32_t u[4];
   80     uint64_t u64[2];
   81 };
   82 #elif defined(HAVE_SSE2)
   83   #include <emmintrin.h>
   84 
   85 /** 128-bit data structure */
   86 union W128_T {
   87     uint32_t u[4];
   88     uint64_t u64[2];
   89     __m128i si;
   90 };
   91 #else
   92 /** 128-bit data structure */
   93 union W128_T {
   94     uint32_t u[4];
   95     uint64_t u64[2];
   96 };
   97 #endif
   98 
   99 /** 128-bit data type */
  100 typedef union W128_T w128_t;
  101 
  102 /**
  103  * SFMT internal state
  104  */
  105 struct SFMT_T {
  106     /** the 128-bit internal state array */
  107     w128_t state[SFMT_N];
  108     /** index counter to the 32-bit internal state array */
  109     int idx;
  110 };
  111 
  112 typedef struct SFMT_T sfmt_t;
  113 
  114 void sfmt_fill_array32(sfmt_t * sfmt, uint32_t * array, int size);
  115 void sfmt_fill_array64(sfmt_t * sfmt, uint64_t * array, int size);
  116 void sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed);
  117 void sfmt_init_by_array(sfmt_t * sfmt, uint32_t * init_key, int key_length);
  118 const char * sfmt_get_idstring(sfmt_t * sfmt);
  119 int sfmt_get_min_array_size32(sfmt_t * sfmt);
  120 int sfmt_get_min_array_size64(sfmt_t * sfmt);
  121 void sfmt_gen_rand_all(sfmt_t * sfmt);
  122 
  123 #ifndef ONLY64
  124 /**
  125  * This function generates and returns 32-bit pseudorandom number.
  126  * init_gen_rand or init_by_array must be called before this function.
  127  * @param sfmt SFMT internal state
  128  * @return 32-bit pseudorandom number
  129  */
  130 inline static uint32_t sfmt_genrand_uint32(sfmt_t * sfmt) {
  131     uint32_t r;
  132     uint32_t * psfmt32 = &sfmt->state[0].u[0];
  133 
  134     if (sfmt->idx >= SFMT_N32) {
  135         sfmt_gen_rand_all(sfmt);
  136         sfmt->idx = 0;
  137     }
  138     r = psfmt32[sfmt->idx++];
  139     return r;
  140 }
  141 #endif
  142 /**
  143  * This function generates and returns 64-bit pseudorandom number.
  144  * init_gen_rand or init_by_array must be called before this function.
  145  * The function gen_rand64 should not be called after gen_rand32,
  146  * unless an initialization is again executed.
  147  * @param sfmt SFMT internal state
  148  * @return 64-bit pseudorandom number
  149  */
  150 inline static uint64_t sfmt_genrand_uint64(sfmt_t * sfmt) {
  151 #if defined(BIG_ENDIAN64) && !defined(ONLY64)
  152     uint32_t * psfmt32 = &sfmt->state[0].u[0];
  153     uint32_t r1, r2;
  154 #else
  155     uint64_t r;
  156 #endif
  157     uint64_t * psfmt64 = &sfmt->state[0].u64[0];
  158     assert(sfmt->idx % 2 == 0);
  159 
  160     if (sfmt->idx >= SFMT_N32) {
  161         sfmt_gen_rand_all(sfmt);
  162         sfmt->idx = 0;
  163     }
  164 #if defined(BIG_ENDIAN64) && !defined(ONLY64)
  165     r1 = psfmt32[sfmt->idx];
  166     r2 = psfmt32[sfmt->idx + 1];
  167     sfmt->idx += 2;
  168     return ((uint64_t)r2 << 32) | r1;
  169 #else
  170     r = psfmt64[sfmt->idx / 2];
  171     sfmt->idx += 2;
  172     return r;
  173 #endif
  174 }
  175 
  176 /* =================================================
  177    The following real versions are due to Isaku Wada
  178    ================================================= */
  179 /**
  180  * converts an unsigned 32-bit number to a double on [0,1]-real-interval.
  181  * @param v 32-bit unsigned integer
  182  * @return double on [0,1]-real-interval
  183  */
  184 inline static double sfmt_to_real1(uint32_t v)
  185 {
  186     return v * (1.0/4294967295.0);
  187     /* divided by 2^32-1 */
  188 }
  189 
  190 /**
  191  * generates a random number on [0,1]-real-interval
  192  * @param sfmt SFMT internal state
  193  * @return double on [0,1]-real-interval
  194  */
  195 inline static double sfmt_genrand_real1(sfmt_t * sfmt)
  196 {
  197     return sfmt_to_real1(sfmt_genrand_uint32(sfmt));
  198 }
  199 
  200 /**
  201  * converts an unsigned 32-bit integer to a double on [0,1)-real-interval.
  202  * @param v 32-bit unsigned integer
  203  * @return double on [0,1)-real-interval
  204  */
  205 inline static double sfmt_to_real2(uint32_t v)
  206 {
  207     return v * (1.0/4294967296.0);
  208     /* divided by 2^32 */
  209 }
  210 
  211 /**
  212  * generates a random number on [0,1)-real-interval
  213  * @param sfmt SFMT internal state
  214  * @return double on [0,1)-real-interval
  215  */
  216 inline static double sfmt_genrand_real2(sfmt_t * sfmt)
  217 {
  218     return sfmt_to_real2(sfmt_genrand_uint32(sfmt));
  219 }
  220 
  221 /**
  222  * converts an unsigned 32-bit integer to a double on (0,1)-real-interval.
  223  * @param v 32-bit unsigned integer
  224  * @return double on (0,1)-real-interval
  225  */
  226 inline static double sfmt_to_real3(uint32_t v)
  227 {
  228     return (((double)v) + 0.5)*(1.0/4294967296.0);
  229     /* divided by 2^32 */
  230 }
  231 
  232 /**
  233  * generates a random number on (0,1)-real-interval
  234  * @param sfmt SFMT internal state
  235  * @return double on (0,1)-real-interval
  236  */
  237 inline static double sfmt_genrand_real3(sfmt_t * sfmt)
  238 {
  239     return sfmt_to_real3(sfmt_genrand_uint32(sfmt));
  240 }
  241 
  242 /**
  243  * converts an unsigned 32-bit integer to double on [0,1)
  244  * with 53-bit resolution.
  245  * @param v 32-bit unsigned integer
  246  * @return double on [0,1)-real-interval with 53-bit resolution.
  247  */
  248 inline static double sfmt_to_res53(uint64_t v)
  249 {
  250     return (v >> 11) * (1.0/9007199254740992.0);
  251 }
  252 
  253 /**
  254  * generates a random number on [0,1) with 53-bit resolution
  255  * @param sfmt SFMT internal state
  256  * @return double on [0,1) with 53-bit resolution
  257  */
  258 inline static double sfmt_genrand_res53(sfmt_t * sfmt)
  259 {
  260     return sfmt_to_res53(sfmt_genrand_uint64(sfmt));
  261 }
  262 
  263 
  264 /* =================================================
  265    The following function are added by Saito.
  266    ================================================= */
  267 /**
  268  * generates a random number on [0,1) with 53-bit resolution from two
  269  * 32 bit integers
  270  */
  271 inline static double sfmt_to_res53_mix(uint32_t x, uint32_t y)
  272 {
  273     return sfmt_to_res53(x | ((uint64_t)y << 32));
  274 }
  275 
  276 /**
  277  * generates a random number on [0,1) with 53-bit resolution
  278  * using two 32bit integers.
  279  * @param sfmt SFMT internal state
  280  * @return double on [0,1) with 53-bit resolution
  281  */
  282 inline static double sfmt_genrand_res53_mix(sfmt_t * sfmt)
  283 {
  284     uint32_t x, y;
  285 
  286     x = sfmt_genrand_uint32(sfmt);
  287     y = sfmt_genrand_uint32(sfmt);
  288     return sfmt_to_res53_mix(x, y);
  289 }
  290 
  291 #if defined(__cplusplus)
  292 }
  293 #endif
  294 
  295 #endif