diff options
Diffstat (limited to 'numpy/random/src/dsfmt/dSFMT.c')
| -rw-r--r-- | numpy/random/src/dsfmt/dSFMT.c | 626 |
1 files changed, 0 insertions, 626 deletions
diff --git a/numpy/random/src/dsfmt/dSFMT.c b/numpy/random/src/dsfmt/dSFMT.c deleted file mode 100644 index 0f122c26c..000000000 --- a/numpy/random/src/dsfmt/dSFMT.c +++ /dev/null @@ -1,626 +0,0 @@ -/** - * @file dSFMT.c - * @brief double precision SIMD-oriented Fast Mersenne Twister (dSFMT) - * based on IEEE 754 format. - * - * @author Mutsuo Saito (Hiroshima University) - * @author Makoto Matsumoto (Hiroshima University) - * - * Copyright (C) 2007,2008 Mutsuo Saito, Makoto Matsumoto and Hiroshima - * University. All rights reserved. - * - * The new BSD License is applied to this software, see LICENSE.txt - */ -#include <stdio.h> -#include <stdlib.h> -#include <string.h> - -#include "dSFMT-params.h" - -#include "dSFMT-common.h" -#include "dSFMT-jump.h" -#include "dSFMT-poly.h" - -#if defined(__cplusplus) -extern "C" { -#endif - -/** dsfmt internal state vector */ -dsfmt_t dsfmt_global_data; -/** dsfmt mexp for check */ -static const int dsfmt_mexp = DSFMT_MEXP; - -/*---------------- - STATIC FUNCTIONS - ----------------*/ -inline static uint32_t ini_func1(uint32_t x); -inline static uint32_t ini_func2(uint32_t x); -inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, w128_t *array, int size); -inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, w128_t *array, int size); -inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, w128_t *array, int size); -inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, w128_t *array, int size); -inline static int idxof(int i); -static void initial_mask(dsfmt_t *dsfmt); -static void period_certification(dsfmt_t *dsfmt); - -#if defined(HAVE_SSE2) -/** 1 in 64bit for sse2 */ -static const union X128I_T sse2_int_one = {{1, 1}}; -/** 2.0 double for sse2 */ -static const union X128D_T sse2_double_two = {{2.0, 2.0}}; -/** -1.0 double for sse2 */ -static const union X128D_T sse2_double_m_one = {{-1.0, -1.0}}; -#endif - -/** - * This function simulate a 32-bit array index overlapped to 64-bit - * array of LITTLE ENDIAN in BIG ENDIAN machine. - */ -#if defined(DSFMT_BIG_ENDIAN) -inline static int idxof(int i) { return i ^ 1; } -#else -inline static int idxof(int i) { return i; } -#endif - -#if defined(HAVE_SSE2) -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range [0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_c0o1(w128_t *w) { - w->sd = _mm_add_pd(w->sd, sse2_double_m_one.d128); -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1]. - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0c1(w128_t *w) { - w->sd = _mm_sub_pd(sse2_double_two.d128, w->sd); -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0o1(w128_t *w) { - w->si = _mm_or_si128(w->si, sse2_int_one.i128); - w->sd = _mm_add_pd(w->sd, sse2_double_m_one.d128); -} -#else /* standard C and altivec */ -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range [0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_c0o1(w128_t *w) { - w->d[0] -= 1.0; - w->d[1] -= 1.0; -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1]. - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0c1(w128_t *w) { - w->d[0] = 2.0 - w->d[0]; - w->d[1] = 2.0 - w->d[1]; -} - -/** - * This function converts the double precision floating point numbers which - * distribute uniformly in the range [1, 2) to those which distribute uniformly - * in the range (0, 1). - * @param w 128bit stracture of double precision floating point numbers (I/O) - */ -inline static void convert_o0o1(w128_t *w) { - w->u[0] |= 1; - w->u[1] |= 1; - w->d[0] -= 1.0; - w->d[1] -= 1.0; -} -#endif - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_c1o2(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_c0o1(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - convert_c0o1(&array[i - DSFMT_N]); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - convert_c0o1(&array[i - DSFMT_N]); - } - for (i = size - DSFMT_N; i < size; i++) { - convert_c0o1(&array[i]); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_o0o1(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - convert_o0o1(&array[i - DSFMT_N]); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - convert_o0o1(&array[i - DSFMT_N]); - } - for (i = size - DSFMT_N; i < size; i++) { - convert_o0o1(&array[i]); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function fills the user-specified array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - * @param array an 128-bit array to be filled by pseudorandom numbers. - * @param size number of 128-bit pseudorandom numbers to be generated. - */ -inline static void gen_rand_array_o0c1(dsfmt_t *dsfmt, w128_t *array, - int size) { - int i, j; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&array[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&array[i], &dsfmt->status[i], &dsfmt->status[i + DSFMT_POS1], - &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&array[i], &dsfmt->status[i], &array[i + DSFMT_POS1 - DSFMT_N], - &lung); - } - for (; i < size - DSFMT_N; i++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - convert_o0c1(&array[i - DSFMT_N]); - } - for (j = 0; j < 2 * DSFMT_N - size; j++) { - dsfmt->status[j] = array[j + size - DSFMT_N]; - } - for (; i < size; i++, j++) { - do_recursion(&array[i], &array[i - DSFMT_N], - &array[i + DSFMT_POS1 - DSFMT_N], &lung); - dsfmt->status[j] = array[i]; - convert_o0c1(&array[i - DSFMT_N]); - } - for (i = size - DSFMT_N; i < size; i++) { - convert_o0c1(&array[i]); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function represents a function used in the initialization - * by init_by_array - * @param x 32-bit integer - * @return 32-bit integer - */ -static uint32_t ini_func1(uint32_t x) { - return (x ^ (x >> 27)) * (uint32_t)1664525UL; -} - -/** - * This function represents a function used in the initialization - * by init_by_array - * @param x 32-bit integer - * @return 32-bit integer - */ -static uint32_t ini_func2(uint32_t x) { - return (x ^ (x >> 27)) * (uint32_t)1566083941UL; -} - -/** - * This function initializes the internal state array to fit the IEEE - * 754 format. - * @param dsfmt dsfmt state vector. - */ -static void initial_mask(dsfmt_t *dsfmt) { - int i; - uint64_t *psfmt; - - psfmt = &dsfmt->status[0].u[0]; - for (i = 0; i < DSFMT_N * 2; i++) { - psfmt[i] = (psfmt[i] & DSFMT_LOW_MASK) | DSFMT_HIGH_CONST; - } -} - -/** - * This function certificate the period of 2^{SFMT_MEXP}-1. - * @param dsfmt dsfmt state vector. - */ -static void period_certification(dsfmt_t *dsfmt) { - uint64_t pcv[2] = {DSFMT_PCV1, DSFMT_PCV2}; - uint64_t tmp[2]; - uint64_t inner; - int i; -#if (DSFMT_PCV2 & 1) != 1 - int j; - uint64_t work; -#endif - - tmp[0] = (dsfmt->status[DSFMT_N].u[0] ^ DSFMT_FIX1); - tmp[1] = (dsfmt->status[DSFMT_N].u[1] ^ DSFMT_FIX2); - - inner = tmp[0] & pcv[0]; - inner ^= tmp[1] & pcv[1]; - for (i = 32; i > 0; i >>= 1) { - inner ^= inner >> i; - } - inner &= 1; - /* check OK */ - if (inner == 1) { - return; - } - /* check NG, and modification */ -#if (DSFMT_PCV2 & 1) == 1 - dsfmt->status[DSFMT_N].u[1] ^= 1; -#else - for (i = 1; i >= 0; i--) { - work = 1; - for (j = 0; j < 64; j++) { - if ((work & pcv[i]) != 0) { - dsfmt->status[DSFMT_N].u[i] ^= work; - return; - } - work = work << 1; - } - } -#endif - return; -} - -/*---------------- - PUBLIC FUNCTIONS - ----------------*/ -/** - * This function returns the identification string. The string shows - * the Mersenne exponent, and all parameters of this generator. - * @return id string. - */ -const char *dsfmt_get_idstring(void) { return DSFMT_IDSTR; } - -/** - * This function returns the minimum size of array used for \b - * fill_array functions. - * @return minimum size of array used for fill_array functions. - */ -int dsfmt_get_min_array_size(void) { return DSFMT_N64; } - -/** - * This function fills the internal state array with double precision - * floating point pseudorandom numbers of the IEEE 754 format. - * @param dsfmt dsfmt state vector. - */ -void dsfmt_gen_rand_all(dsfmt_t *dsfmt) { - int i; - w128_t lung; - - lung = dsfmt->status[DSFMT_N]; - do_recursion(&dsfmt->status[0], &dsfmt->status[0], &dsfmt->status[DSFMT_POS1], - &lung); - for (i = 1; i < DSFMT_N - DSFMT_POS1; i++) { - do_recursion(&dsfmt->status[i], &dsfmt->status[i], - &dsfmt->status[i + DSFMT_POS1], &lung); - } - for (; i < DSFMT_N; i++) { - do_recursion(&dsfmt->status[i], &dsfmt->status[i], - &dsfmt->status[i + DSFMT_POS1 - DSFMT_N], &lung); - } - dsfmt->status[DSFMT_N] = lung; -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range [1, 2) to the - * specified array[] by one call. The number of pseudorandom numbers - * is specified by the argument \b size, which must be at least (SFMT_MEXP - * / 128) * 2 and a multiple of two. The function - * get_min_array_size() returns this minimum size. The generation by - * this function is much faster than the following fill_array_xxx functions. - * - * For initialization, init_gen_rand() or init_by_array() must be called - * before the first call of this function. This function can not be - * used after calling genrand_xxx functions, without initialization. - * - * @param dsfmt dsfmt state vector. - * @param array an array where pseudorandom numbers are filled - * by this function. The pointer to the array must be "aligned" - * (namely, must be a multiple of 16) in the SIMD version, since it - * refers to the address of a 128-bit integer. In the standard C - * version, the pointer is arbitrary. - * - * @param size the number of 64-bit pseudorandom integers to be - * generated. size must be a multiple of 2, and greater than or equal - * to (SFMT_MEXP / 128) * 2. - * - * @note \b memalign or \b posix_memalign is available to get aligned - * memory. Mac OSX doesn't have these functions, but \b malloc of OSX - * returns the pointer to the aligned memory block. - */ -void dsfmt_fill_array_close1_open2(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_c1o2(dsfmt, (w128_t *)array, size / 2); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range (0, 1] to the - * specified array[] by one call. This function is the same as - * fill_array_close1_open2() except the distribution range. - * - * @param dsfmt dsfmt state vector. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa fill_array_close1_open2() - */ -void dsfmt_fill_array_open_close(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_o0c1(dsfmt, (w128_t *)array, size / 2); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range [0, 1) to the - * specified array[] by one call. This function is the same as - * fill_array_close1_open2() except the distribution range. - * - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param dsfmt dsfmt state vector. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa fill_array_close1_open2() - */ -void dsfmt_fill_array_close_open(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_c0o1(dsfmt, (w128_t *)array, size / 2); -} - -/** - * This function generates double precision floating point - * pseudorandom numbers which distribute in the range (0, 1) to the - * specified array[] by one call. This function is the same as - * fill_array_close1_open2() except the distribution range. - * - * @param dsfmt dsfmt state vector. - * @param array an array where pseudorandom numbers are filled - * by this function. - * @param size the number of pseudorandom numbers to be generated. - * see also \sa fill_array_close1_open2() - */ -void dsfmt_fill_array_open_open(dsfmt_t *dsfmt, double array[], int size) { - assert(size % 2 == 0); - assert(size >= DSFMT_N64); - gen_rand_array_o0o1(dsfmt, (w128_t *)array, size / 2); -} - -#if defined(__INTEL_COMPILER) -#pragma warning(disable : 981) -#endif -/** - * This function initializes the internal state array with a 32-bit - * integer seed. - * @param dsfmt dsfmt state vector. - * @param seed a 32-bit integer used as the seed. - * @param mexp caller's mersenne expornent - */ -void dsfmt_chk_init_gen_rand(dsfmt_t *dsfmt, uint32_t seed, int mexp) { - int i; - uint32_t *psfmt; - - /* make sure caller program is compiled with the same MEXP */ - if (mexp != dsfmt_mexp) { - fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n"); - exit(1); - } - psfmt = &dsfmt->status[0].u32[0]; - psfmt[idxof(0)] = seed; - for (i = 1; i < (DSFMT_N + 1) * 4; i++) { - psfmt[idxof(i)] = - 1812433253UL * (psfmt[idxof(i - 1)] ^ (psfmt[idxof(i - 1)] >> 30)) + i; - } - initial_mask(dsfmt); - period_certification(dsfmt); - dsfmt->idx = DSFMT_N64; -} - -/** - * This function initializes the internal state array, - * with an array of 32-bit integers used as the seeds - * @param dsfmt dsfmt state vector. - * @param init_key the array of 32-bit integers, used as a seed. - * @param key_length the length of init_key. - * @param mexp caller's mersenne expornent - */ -void dsfmt_chk_init_by_array(dsfmt_t *dsfmt, uint32_t init_key[], - int key_length, int mexp) { - int i, j, count; - uint32_t r; - uint32_t *psfmt32; - int lag; - int mid; - int size = (DSFMT_N + 1) * 4; /* pulmonary */ - - /* make sure caller program is compiled with the same MEXP */ - if (mexp != dsfmt_mexp) { - fprintf(stderr, "DSFMT_MEXP doesn't match with dSFMT.c\n"); - exit(1); - } - if (size >= 623) { - lag = 11; - } else if (size >= 68) { - lag = 7; - } else if (size >= 39) { - lag = 5; - } else { - lag = 3; - } - mid = (size - lag) / 2; - - psfmt32 = &dsfmt->status[0].u32[0]; - memset(dsfmt->status, 0x8b, sizeof(dsfmt->status)); - if (key_length + 1 > size) { - count = key_length + 1; - } else { - count = size; - } - r = ini_func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid % size)] ^ - psfmt32[idxof((size - 1) % size)]); - psfmt32[idxof(mid % size)] += r; - r += key_length; - psfmt32[idxof((mid + lag) % size)] += r; - psfmt32[idxof(0)] = r; - count--; - for (i = 1, j = 0; (j < count) && (j < key_length); j++) { - r = ini_func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % size)] ^ - psfmt32[idxof((i + size - 1) % size)]); - psfmt32[idxof((i + mid) % size)] += r; - r += init_key[j] + i; - psfmt32[idxof((i + mid + lag) % size)] += r; - psfmt32[idxof(i)] = r; - i = (i + 1) % size; - } - for (; j < count; j++) { - r = ini_func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % size)] ^ - psfmt32[idxof((i + size - 1) % size)]); - psfmt32[idxof((i + mid) % size)] += r; - r += i; - psfmt32[idxof((i + mid + lag) % size)] += r; - psfmt32[idxof(i)] = r; - i = (i + 1) % size; - } - for (j = 0; j < size; j++) { - r = ini_func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % size)] + - psfmt32[idxof((i + size - 1) % size)]); - psfmt32[idxof((i + mid) % size)] ^= r; - r -= i; - psfmt32[idxof((i + mid + lag) % size)] ^= r; - psfmt32[idxof(i)] = r; - i = (i + 1) % size; - } - initial_mask(dsfmt); - period_certification(dsfmt); - dsfmt->idx = DSFMT_N64; -} -#if defined(__INTEL_COMPILER) -#pragma warning(default : 981) -#endif - -#if defined(__cplusplus) -} -#endif - -extern inline double dsfmt_next_double(dsfmt_state *state); - -extern inline uint64_t dsfmt_next64(dsfmt_state *state); - -extern inline uint32_t dsfmt_next32(dsfmt_state *state); - -void dsfmt_jump(dsfmt_state *state) { dSFMT_jump(state->state, poly_128); };
\ No newline at end of file |