Mercurial > vba-clojure

     1 /** 

     2  * @file SFMT-alti.h 

     3  *

     4  * @brief SIMD oriented Fast Mersenne Twister(SFMT)

     5  * pseudorandom number generator

     6  *

     7  * @author Mutsuo Saito (Hiroshima University)

     8  * @author Makoto Matsumoto (Hiroshima University)

     9  *

    10  * Copyright (C) 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima

    11  * University. All rights reserved.

    12  *

    13  * The new BSD License is applied to this software.

    14  * see LICENSE.txt

    15  */

    16 

    17 #ifndef SFMT_ALTI_H

    18 #define SFMT_ALTI_H

    19 

    20 inline static vector unsigned int vec_recursion(vector unsigned int a,

    21 						vector unsigned int b,

    22 						vector unsigned int c,

    23 						vector unsigned int d)

    24     ALWAYSINLINE;

    25 

    26 /**

    27  * This function represents the recursion formula in AltiVec and BIG ENDIAN.

    28  * @param a a 128-bit part of the interal state array

    29  * @param b a 128-bit part of the interal state array

    30  * @param c a 128-bit part of the interal state array

    31  * @param d a 128-bit part of the interal state array

    32  * @return output

    33  */

    34 inline static vector unsigned int vec_recursion(vector unsigned int a,

    35 						vector unsigned int b,

    36 						vector unsigned int c,

    37 						vector unsigned int d) {

    38 

    39     const vector unsigned int sl1 = ALTI_SL1;

    40     const vector unsigned int sr1 = ALTI_SR1;

    41 #ifdef ONLY64

    42     const vector unsigned int mask = ALTI_MSK64;

    43     const vector unsigned char perm_sl = ALTI_SL2_PERM64;

    44     const vector unsigned char perm_sr = ALTI_SR2_PERM64;

    45 #else

    46     const vector unsigned int mask = ALTI_MSK;

    47     const vector unsigned char perm_sl = ALTI_SL2_PERM;

    48     const vector unsigned char perm_sr = ALTI_SR2_PERM;

    49 #endif

    50     vector unsigned int v, w, x, y, z;

    51     x = vec_perm(a, (vector unsigned int)perm_sl, perm_sl);

    52     v = a;

    53     y = vec_sr(b, sr1);

    54     z = vec_perm(c, (vector unsigned int)perm_sr, perm_sr);

    55     w = vec_sl(d, sl1);

    56     z = vec_xor(z, w);

    57     y = vec_and(y, mask);

    58     v = vec_xor(v, x);

    59     z = vec_xor(z, y);

    60     z = vec_xor(z, v);

    61     return z;

    62 }

    63 

    64 /**

    65  * This function fills the internal state array with pseudorandom

    66  * integers.

    67  */

    68 inline static void gen_rand_all(void) {

    69     int i;

    70     vector unsigned int r, r1, r2;

    71 

    72     r1 = sfmt[N - 2].s;

    73     r2 = sfmt[N - 1].s;

    74     for (i = 0; i < N - POS1; i++) {

    75 	r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);

    76 	sfmt[i].s = r;

    77 	r1 = r2;

    78 	r2 = r;

    79     }

    80     for (; i < N; i++) {

    81 	r = vec_recursion(sfmt[i].s, sfmt[i + POS1 - N].s, r1, r2);

    82 	sfmt[i].s = r;

    83 	r1 = r2;

    84 	r2 = r;

    85     }

    86 }

    87 

    88 /**

    89  * This function fills the user-specified array with pseudorandom

    90  * integers.

    91  *

    92  * @param array an 128-bit array to be filled by pseudorandom numbers.  

    93  * @param size number of 128-bit pesudorandom numbers to be generated.

    94  */

    95 inline static void gen_rand_array(w128_t *array, int size) {

    96     int i, j;

    97     vector unsigned int r, r1, r2;

    98 

    99     r1 = sfmt[N - 2].s;

   100     r2 = sfmt[N - 1].s;

   101     for (i = 0; i < N - POS1; i++) {

   102 	r = vec_recursion(sfmt[i].s, sfmt[i + POS1].s, r1, r2);

   103 	array[i].s = r;

   104 	r1 = r2;

   105 	r2 = r;

   106     }

   107     for (; i < N; i++) {

   108 	r = vec_recursion(sfmt[i].s, array[i + POS1 - N].s, r1, r2);

   109 	array[i].s = r;

   110 	r1 = r2;

   111 	r2 = r;

   112     }

   113     /* main loop */

   114     for (; i < size - N; i++) {

   115 	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);

   116 	array[i].s = r;

   117 	r1 = r2;

   118 	r2 = r;

   119     }

   120     for (j = 0; j < 2 * N - size; j++) {

   121 	sfmt[j].s = array[j + size - N].s;

   122     }

   123     for (; i < size; i++) {

   124 	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);

   125 	array[i].s = r;

   126 	sfmt[j++].s = r;

   127 	r1 = r2;

   128 	r2 = r;

   129     }

   130 }

   131 

   132 #ifndef ONLY64

   133 #if defined(__APPLE__)

   134 #define ALTI_SWAP (vector unsigned char) \

   135 	(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11)

   136 #else

   137 #define ALTI_SWAP {4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11}

   138 #endif

   139 /**

   140  * This function swaps high and low 32-bit of 64-bit integers in user

   141  * specified array.

   142  *

   143  * @param array an 128-bit array to be swaped.

   144  * @param size size of 128-bit array.

   145  */

   146 inline static void swap(w128_t *array, int size) {

   147     int i;

   148     const vector unsigned char perm = ALTI_SWAP;

   149 

   150     for (i = 0; i < size; i++) {

   151 	array[i].s = vec_perm(array[i].s, (vector unsigned int)perm, perm);

   152     }

   153 }

   154 #endif

   155 

   156 #endif
author	Robert McIntyre <rlm@mit.edu>
date	Fri, 23 Mar 2012 00:30:34 -0500
parents	f9f4f1b99eed
children