Home | History | Annotate | Download | only in test
      1 /*
      2  * This file derives from SFMT 1.3.3
      3  * (http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html), which was
      4  * released under the terms of the following license:
      5  *
      6  *   Copyright (c) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
      7  *   University. All rights reserved.
      8  *
      9  *   Redistribution and use in source and binary forms, with or without
     10  *   modification, are permitted provided that the following conditions are
     11  *   met:
     12  *
     13  *       * Redistributions of source code must retain the above copyright
     14  *         notice, this list of conditions and the following disclaimer.
     15  *       * Redistributions in binary form must reproduce the above
     16  *         copyright notice, this list of conditions and the following
     17  *         disclaimer in the documentation and/or other materials provided
     18  *         with the distribution.
     19  *       * Neither the name of the Hiroshima University nor the names of
     20  *         its contributors may be used to endorse or promote products
     21  *         derived from this software without specific prior written
     22  *         permission.
     23  *
     24  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
     25  *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
     26  *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
     27  *   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
     28  *   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
     29  *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
     30  *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     31  *   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
     32  *   THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
     33  *   (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
     34  *   OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     35  */
     36 /**
     37  * @file SFMT-alti.h
     38  *
     39  * @brief SIMD oriented Fast Mersenne Twister(SFMT)
     40  * pseudorandom number generator
     41  *
     42  * @author Mutsuo Saito (Hiroshima University)
     43  * @author Makoto Matsumoto (Hiroshima University)
     44  *
     45  * Copyright (C) 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
     46  * University. All rights reserved.
     47  *
     48  * The new BSD License is applied to this software.
     49  * see LICENSE.txt
     50  */
     51 
     52 #ifndef SFMT_ALTI_H
     53 #define SFMT_ALTI_H
     54 
     55 /**
     56  * This function represents the recursion formula in AltiVec and BIG ENDIAN.
     57  * @param a a 128-bit part of the interal state array
     58  * @param b a 128-bit part of the interal state array
     59  * @param c a 128-bit part of the interal state array
     60  * @param d a 128-bit part of the interal state array
     61  * @return output
     62  */
     63 JEMALLOC_ALWAYS_INLINE
     64 vector unsigned int vec_recursion(vector unsigned int a,
     65 						vector unsigned int b,
     66 						vector unsigned int c,
     67 						vector unsigned int d) {
     68 
     69     const vector unsigned int sl1 = ALTI_SL1;
     70     const vector unsigned int sr1 = ALTI_SR1;
     71 #ifdef ONLY64
     72     const vector unsigned int mask = ALTI_MSK64;
     73     const vector unsigned char perm_sl = ALTI_SL2_PERM64;
     74     const vector unsigned char perm_sr = ALTI_SR2_PERM64;
     75 #else
     76     const vector unsigned int mask = ALTI_MSK;
     77     const vector unsigned char perm_sl = ALTI_SL2_PERM;
     78     const vector unsigned char perm_sr = ALTI_SR2_PERM;
     79 #endif
     80     vector unsigned int v, w, x, y, z;
     81     x = vec_perm(a, (vector unsigned int)perm_sl, perm_sl);
     82     v = a;
     83     y = vec_sr(b, sr1);
     84     z = vec_perm(c, (vector unsigned int)perm_sr, perm_sr);
     85     w = vec_sl(d, sl1);
     86     z = vec_xor(z, w);
     87     y = vec_and(y, mask);
     88     v = vec_xor(v, x);
     89     z = vec_xor(z, y);
     90     z = vec_xor(z, v);
     91     return z;
     92 }
     93 
     94 /**
     95  * This function fills the internal state array with pseudorandom
     96  * integers.
     97  */
     98 JEMALLOC_INLINE void gen_rand_all(sfmt_t *ctx) {
     99     int i;
    100     vector unsigned int r, r1, r2;
    101 
    102     r1 = ctx->sfmt[N - 2].s;
    103     r2 = ctx->sfmt[N - 1].s;
    104     for (i = 0; i < N - POS1; i++) {
    105 	r = vec_recursion(ctx->sfmt[i].s, ctx->sfmt[i + POS1].s, r1, r2);
    106 	ctx->sfmt[i].s = r;
    107 	r1 = r2;
    108 	r2 = r;
    109     }
    110     for (; i < N; i++) {
    111 	r = vec_recursion(ctx->sfmt[i].s, ctx->sfmt[i + POS1 - N].s, r1, r2);
    112 	ctx->sfmt[i].s = r;
    113 	r1 = r2;
    114 	r2 = r;
    115     }
    116 }
    117 
    118 /**
    119  * This function fills the user-specified array with pseudorandom
    120  * integers.
    121  *
    122  * @param array an 128-bit array to be filled by pseudorandom numbers.
    123  * @param size number of 128-bit pesudorandom numbers to be generated.
    124  */
    125 JEMALLOC_INLINE void gen_rand_array(sfmt_t *ctx, w128_t *array, int size) {
    126     int i, j;
    127     vector unsigned int r, r1, r2;
    128 
    129     r1 = ctx->sfmt[N - 2].s;
    130     r2 = ctx->sfmt[N - 1].s;
    131     for (i = 0; i < N - POS1; i++) {
    132 	r = vec_recursion(ctx->sfmt[i].s, ctx->sfmt[i + POS1].s, r1, r2);
    133 	array[i].s = r;
    134 	r1 = r2;
    135 	r2 = r;
    136     }
    137     for (; i < N; i++) {
    138 	r = vec_recursion(ctx->sfmt[i].s, array[i + POS1 - N].s, r1, r2);
    139 	array[i].s = r;
    140 	r1 = r2;
    141 	r2 = r;
    142     }
    143     /* main loop */
    144     for (; i < size - N; i++) {
    145 	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
    146 	array[i].s = r;
    147 	r1 = r2;
    148 	r2 = r;
    149     }
    150     for (j = 0; j < 2 * N - size; j++) {
    151 	ctx->sfmt[j].s = array[j + size - N].s;
    152     }
    153     for (; i < size; i++) {
    154 	r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
    155 	array[i].s = r;
    156 	ctx->sfmt[j++].s = r;
    157 	r1 = r2;
    158 	r2 = r;
    159     }
    160 }
    161 
    162 #ifndef ONLY64
    163 #if defined(__APPLE__)
    164 #define ALTI_SWAP (vector unsigned char) \
    165 	(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11)
    166 #else
    167 #define ALTI_SWAP {4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11}
    168 #endif
    169 /**
    170  * This function swaps high and low 32-bit of 64-bit integers in user
    171  * specified array.
    172  *
    173  * @param array an 128-bit array to be swaped.
    174  * @param size size of 128-bit array.
    175  */
    176 JEMALLOC_INLINE void swap(w128_t *array, int size) {
    177     int i;
    178     const vector unsigned char perm = ALTI_SWAP;
    179 
    180     for (i = 0; i < size; i++) {
    181 	array[i].s = vec_perm(array[i].s, (vector unsigned int)perm, perm);
    182     }
    183 }
    184 #endif
    185 
    186 #endif
    187