Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

added SFMT lib to 3rdparty

  • Loading branch information...
commit 84b34839191986145d3060204e45fa0d49622459 1 parent 0811bdd
@liuliu authored
View
14 COPYING
@@ -39,4 +39,18 @@ Redistribution and use in source and binary forms, with or without modification,
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+SFMT:
+
+Copyright (c) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima University.
+Copyright (c) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima University and The University of Tokyo.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+* Neither the names of Hiroshima University, The University of Tokyo nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
sha1: TODO
View
156 lib/3rdparty/sfmt/SFMT-alti.h
@@ -0,0 +1,156 @@
+#pragma once
+/**
+ * @file SFMT-alti.h
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT)
+ * pseudorandom number generator
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software.
+ * see LICENSE.txt
+ */
+
+#ifndef SFMT_ALTI_H
+#define SFMT_ALTI_H
+
+inline static vector unsigned int vec_recursion(vector unsigned int a,
+ vector unsigned int b,
+ vector unsigned int c,
+ vector unsigned int d);
+
+/**
+ * This function represents the recursion formula in AltiVec and BIG ENDIAN.
+ * @param a a 128-bit part of the interal state array
+ * @param b a 128-bit part of the interal state array
+ * @param c a 128-bit part of the interal state array
+ * @param d a 128-bit part of the interal state array
+ * @return output
+ */
+inline static vector unsigned int vec_recursion(vector unsigned int a,
+ vector unsigned int b,
+ vector unsigned int c,
+ vector unsigned int d) {
+
+ const vector unsigned int sl1 = SFMT_ALTI_SL1;
+ const vector unsigned int sr1 = SFMT_ALTI_SR1;
+#ifdef ONLY64
+ const vector unsigned int mask = SFMT_ALTI_MSK64;
+ const vector unsigned char perm_sl = SFMT_ALTI_SL2_PERM64;
+ const vector unsigned char perm_sr = SFMT_ALTI_SR2_PERM64;
+#else
+ const vector unsigned int mask = SFMT_ALTI_MSK;
+ const vector unsigned char perm_sl = SFMT_ALTI_SL2_PERM;
+ const vector unsigned char perm_sr = SFMT_ALTI_SR2_PERM;
+#endif
+ vector unsigned int v, w, x, y, z;
+ x = vec_perm(a, (vector unsigned int)perm_sl, perm_sl);
+ v = a;
+ y = vec_sr(b, sr1);
+ z = vec_perm(c, (vector unsigned int)perm_sr, perm_sr);
+ w = vec_sl(d, sl1);
+ z = vec_xor(z, w);
+ y = vec_and(y, mask);
+ v = vec_xor(v, x);
+ z = vec_xor(z, y);
+ z = vec_xor(z, v);
+ return z;
+}
+
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+ int i;
+ vector unsigned int r, r1, r2;
+
+ r1 = sfmt->state[N - 2].s;
+ r2 = sfmt->state[N - 1].s;
+ for (i = 0; i < N - POS1; i++) {
+ r = vec_recursion(sfmt->state[i].s, sfmt->state[i + POS1].s, r1, r2);
+ sfmt->state[i].s = r;
+ r1 = r2;
+ r2 = r;
+ }
+ for (; i < N; i++) {
+ r = vec_recursion(sfmt->state[i].s, sfmt->state[i + POS1 - N].s, r1, r2);
+ sfmt->state[i].s = r;
+ r1 = r2;
+ r2 = r;
+ }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ *
+ * @param array an 128-bit array to be filled by pseudorandom numbers.
+ * @param size number of 128-bit pesudorandom numbers to be generated.
+ */
+inline static void gen_rand_array(sfmt_t * sfmt, w128_t *array, int size) {
+ int i, j;
+ vector unsigned int r, r1, r2;
+
+ r1 = sfmt->state[N - 2].s;
+ r2 = sfmt->state[N - 1].s;
+ for (i = 0; i < N - POS1; i++) {
+ r = vec_recursion(sfmt->state[i].s, sfmt->state[i + POS1].s, r1, r2);
+ array[i].s = r;
+ r1 = r2;
+ r2 = r;
+ }
+ for (; i < N; i++) {
+ r = vec_recursion(sfmt->state[i].s, array[i + POS1 - N].s, r1, r2);
+ array[i].s = r;
+ r1 = r2;
+ r2 = r;
+ }
+ /* main loop */
+ for (; i < size - N; i++) {
+ r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
+ array[i].s = r;
+ r1 = r2;
+ r2 = r;
+ }
+ for (j = 0; j < 2 * N - size; j++) {
+ sfmt->state[j].s = array[j + size - N].s;
+ }
+ for (; i < size; i++) {
+ r = vec_recursion(array[i - N].s, array[i + POS1 - N].s, r1, r2);
+ array[i].s = r;
+ sfmt->state[j++].s = r;
+ r1 = r2;
+ r2 = r;
+ }
+}
+
+#ifndef ONLY64
+#if defined(__APPLE__)
+#define SFMT_ALTI_SWAP (vector unsigned char) \
+ (4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11)
+#else
+#define SFMT_ALTI_SWAP {4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11}
+#endif
+/**
+ * This function swaps high and low 32-bit of 64-bit integers in user
+ * specified array.
+ *
+ * @param array an 128-bit array to be swaped.
+ * @param size size of 128-bit array.
+ */
+inline static void swap(w128_t *array, int size) {
+ int i;
+ const vector unsigned char perm = SFMT_ALTI_SWAP;
+
+ for (i = 0; i < size; i++) {
+ array[i].s = vec_perm(array[i].s, (vector unsigned int)perm, perm);
+ }
+}
+#endif
+
+#endif
View
164 lib/3rdparty/sfmt/SFMT-common.h
@@ -0,0 +1,164 @@
+#pragma once
+/**
+ * @file SFMT-common.h
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
+ * number generator with jump function. This file includes common functions
+ * used in random number generation and jump.
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (The University of Tokyo)
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University.
+ * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
+ * University and The University of Tokyo.
+ * All rights reserved.
+ *
+ * The 3-clause BSD License is applied to this software, see
+ * LICENSE.txt
+ */
+#ifndef SFMT_COMMON_H
+#define SFMT_COMMON_H
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+#include "SFMT.h"
+
+inline static void do_recursion(w128_t * r, w128_t * a, w128_t * b,
+ w128_t * c, w128_t * d);
+
+inline static void rshift128(w128_t *out, w128_t const *in, int shift);
+inline static void lshift128(w128_t *out, w128_t const *in, int shift);
+
+/**
+ * This function simulates SIMD 128-bit right shift by the standard C.
+ * The 128-bit integer given in in is shifted by (shift * 8) bits.
+ * This function simulates the LITTLE ENDIAN SIMD.
+ * @param out the output of this function
+ * @param in the 128-bit data to be shifted
+ * @param shift the shift value
+ */
+#ifdef ONLY64
+inline static void rshift128(w128_t *out, w128_t const *in, int shift) {
+ uint64_t th, tl, oh, ol;
+
+ th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);
+ tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);
+
+ oh = th >> (shift * 8);
+ ol = tl >> (shift * 8);
+ ol |= th << (64 - shift * 8);
+ out->u[0] = (uint32_t)(ol >> 32);
+ out->u[1] = (uint32_t)ol;
+ out->u[2] = (uint32_t)(oh >> 32);
+ out->u[3] = (uint32_t)oh;
+}
+#else
+inline static void rshift128(w128_t *out, w128_t const *in, int shift)
+{
+ uint64_t th, tl, oh, ol;
+
+ th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
+ tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
+
+ oh = th >> (shift * 8);
+ ol = tl >> (shift * 8);
+ ol |= th << (64 - shift * 8);
+ out->u[1] = (uint32_t)(ol >> 32);
+ out->u[0] = (uint32_t)ol;
+ out->u[3] = (uint32_t)(oh >> 32);
+ out->u[2] = (uint32_t)oh;
+}
+#endif
+/**
+ * This function simulates SIMD 128-bit left shift by the standard C.
+ * The 128-bit integer given in in is shifted by (shift * 8) bits.
+ * This function simulates the LITTLE ENDIAN SIMD.
+ * @param out the output of this function
+ * @param in the 128-bit data to be shifted
+ * @param shift the shift value
+ */
+#ifdef ONLY64
+inline static void lshift128(w128_t *out, w128_t const *in, int shift) {
+ uint64_t th, tl, oh, ol;
+
+ th = ((uint64_t)in->u[2] << 32) | ((uint64_t)in->u[3]);
+ tl = ((uint64_t)in->u[0] << 32) | ((uint64_t)in->u[1]);
+
+ oh = th << (shift * 8);
+ ol = tl << (shift * 8);
+ oh |= tl >> (64 - shift * 8);
+ out->u[0] = (uint32_t)(ol >> 32);
+ out->u[1] = (uint32_t)ol;
+ out->u[2] = (uint32_t)(oh >> 32);
+ out->u[3] = (uint32_t)oh;
+}
+#else
+inline static void lshift128(w128_t *out, w128_t const *in, int shift)
+{
+ uint64_t th, tl, oh, ol;
+
+ th = ((uint64_t)in->u[3] << 32) | ((uint64_t)in->u[2]);
+ tl = ((uint64_t)in->u[1] << 32) | ((uint64_t)in->u[0]);
+
+ oh = th << (shift * 8);
+ ol = tl << (shift * 8);
+ oh |= tl >> (64 - shift * 8);
+ out->u[1] = (uint32_t)(ol >> 32);
+ out->u[0] = (uint32_t)ol;
+ out->u[3] = (uint32_t)(oh >> 32);
+ out->u[2] = (uint32_t)oh;
+}
+#endif
+/**
+ * This function represents the recursion formula.
+ * @param r output
+ * @param a a 128-bit part of the internal state array
+ * @param b a 128-bit part of the internal state array
+ * @param c a 128-bit part of the internal state array
+ * @param d a 128-bit part of the internal state array
+ */
+#ifdef ONLY64
+inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *c,
+ w128_t *d) {
+ w128_t x;
+ w128_t y;
+
+ lshift128(&x, a, SFMT_SL2);
+ rshift128(&y, c, SFMT_SR2);
+ r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SFMT_SR1) & SFMT_MSK2) ^ y.u[0]
+ ^ (d->u[0] << SFMT_SL1);
+ r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SFMT_SR1) & SFMT_MSK1) ^ y.u[1]
+ ^ (d->u[1] << SFMT_SL1);
+ r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SFMT_SR1) & SFMT_MSK4) ^ y.u[2]
+ ^ (d->u[2] << SFMT_SL1);
+ r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SFMT_SR1) & SFMT_MSK3) ^ y.u[3]
+ ^ (d->u[3] << SFMT_SL1);
+}
+#else
+inline static void do_recursion(w128_t *r, w128_t *a, w128_t *b,
+ w128_t *c, w128_t *d)
+{
+ w128_t x;
+ w128_t y;
+
+ lshift128(&x, a, SFMT_SL2);
+ rshift128(&y, c, SFMT_SR2);
+ r->u[0] = a->u[0] ^ x.u[0] ^ ((b->u[0] >> SFMT_SR1) & SFMT_MSK1)
+ ^ y.u[0] ^ (d->u[0] << SFMT_SL1);
+ r->u[1] = a->u[1] ^ x.u[1] ^ ((b->u[1] >> SFMT_SR1) & SFMT_MSK2)
+ ^ y.u[1] ^ (d->u[1] << SFMT_SL1);
+ r->u[2] = a->u[2] ^ x.u[2] ^ ((b->u[2] >> SFMT_SR1) & SFMT_MSK3)
+ ^ y.u[2] ^ (d->u[2] << SFMT_SL1);
+ r->u[3] = a->u[3] ^ x.u[3] ^ ((b->u[3] >> SFMT_SR1) & SFMT_MSK4)
+ ^ y.u[3] ^ (d->u[3] << SFMT_SL1);
+}
+#endif
+#endif
+
+#if defined(__cplusplus)
+}
+#endif
View
65 lib/3rdparty/sfmt/SFMT-params.h
@@ -0,0 +1,65 @@
+#pragma once
+#ifndef SFMT_PARAMS_H
+#define SFMT_PARAMS_H
+
+#define SFMT_MEXP 19937
+/*-----------------
+ BASIC DEFINITIONS
+ -----------------*/
+/** Mersenne Exponent. The period of the sequence
+ * is a multiple of 2^MEXP-1.
+ * #define SFMT_MEXP 19937 */
+/** SFMT generator has an internal state array of 128-bit integers,
+ * and N is its size. */
+#define SFMT_N (SFMT_MEXP / 128 + 1)
+/** N32 is the size of internal state array when regarded as an array
+ * of 32-bit integers.*/
+#define SFMT_N32 (SFMT_N * 4)
+/** N64 is the size of internal state array when regarded as an array
+ * of 64-bit integers.*/
+#define SFMT_N64 (SFMT_N * 2)
+
+/*----------------------
+ the parameters of SFMT
+ following definitions are in paramsXXXX.h file.
+ ----------------------*/
+/** the pick up position of the array.
+#define SFMT_POS1 122
+*/
+
+/** the parameter of shift left as four 32-bit registers.
+#define SFMT_SL1 18
+ */
+
+/** the parameter of shift left as one 128-bit register.
+ * The 128-bit integer is shifted by (SFMT_SL2 * 8) bits.
+#define SFMT_SL2 1
+*/
+
+/** the parameter of shift right as four 32-bit registers.
+#define SFMT_SR1 11
+*/
+
+/** the parameter of shift right as one 128-bit register.
+ * The 128-bit integer is shifted by (SFMT_SL2 * 8) bits.
+#define SFMT_SR21 1
+*/
+
+/** A bitmask, used in the recursion. These parameters are introduced
+ * to break symmetry of SIMD.
+#define SFMT_MSK1 0xdfffffefU
+#define SFMT_MSK2 0xddfecb7fU
+#define SFMT_MSK3 0xbffaffffU
+#define SFMT_MSK4 0xbffffff6U
+*/
+
+/** These definitions are part of a 128-bit period certification vector.
+#define SFMT_PARITY1 0x00000001U
+#define SFMT_PARITY2 0x00000000U
+#define SFMT_PARITY3 0x00000000U
+#define SFMT_PARITY4 0xc98e126aU
+*/
+
+#include "SFMT-params19937.h"
+
+#endif /* SFMT_PARAMS_H */
View
50 lib/3rdparty/sfmt/SFMT-params19937.h
@@ -0,0 +1,50 @@
+#pragma once
+#ifndef SFMT_PARAMS19937_H
+#define SFMT_PARAMS19937_H
+
+#define SFMT_POS1 122
+#define SFMT_SL1 18
+#define SFMT_SL2 1
+#define SFMT_SR1 11
+#define SFMT_SR2 1
+#define SFMT_MSK1 0xdfffffefU
+#define SFMT_MSK2 0xddfecb7fU
+#define SFMT_MSK3 0xbffaffffU
+#define SFMT_MSK4 0xbffffff6U
+#define SFMT_PARITY1 0x00000001U
+#define SFMT_PARITY2 0x00000000U
+#define SFMT_PARITY3 0x00000000U
+#define SFMT_PARITY4 0x13c9e684U
+
+
+/* PARAMETERS FOR ALTIVEC */
+#if defined(__APPLE__) /* For OSX */
+ #define SFMT_ALTI_SL1 \
+ (vector unsigned int)(SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1)
+ #define SFMT_ALTI_SR1 \
+ (vector unsigned int)(SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1)
+ #define SFMT_ALTI_MSK \
+ (vector unsigned int)(SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4)
+ #define SFMT_ALTI_MSK64 \
+ (vector unsigned int)(SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3)
+ #define SFMT_ALTI_SL2_PERM \
+ (vector unsigned char)(1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8)
+ #define SFMT_ALTI_SL2_PERM64 \
+ (vector unsigned char)(1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0)
+ #define SFMT_ALTI_SR2_PERM \
+ (vector unsigned char)(7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14)
+ #define SFMT_ALTI_SR2_PERM64 \
+ (vector unsigned char)(15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14)
+#else /* For OTHER OSs(Linux?) */
+ #define SFMT_ALTI_SL1 {SFMT_SL1, SFMT_SL1, SFMT_SL1, SFMT_SL1}
+ #define SFMT_ALTI_SR1 {SFMT_SR1, SFMT_SR1, SFMT_SR1, SFMT_SR1}
+ #define SFMT_ALTI_MSK {SFMT_MSK1, SFMT_MSK2, SFMT_MSK3, SFMT_MSK4}
+ #define SFMT_ALTI_MSK64 {SFMT_MSK2, SFMT_MSK1, SFMT_MSK4, SFMT_MSK3}
+ #define SFMT_ALTI_SL2_PERM {1,2,3,23,5,6,7,0,9,10,11,4,13,14,15,8}
+ #define SFMT_ALTI_SL2_PERM64 {1,2,3,4,5,6,7,31,9,10,11,12,13,14,15,0}
+ #define SFMT_ALTI_SR2_PERM {7,0,1,2,11,4,5,6,15,8,9,10,17,12,13,14}
+ #define SFMT_ALTI_SR2_PERM64 {15,0,1,2,3,4,5,6,17,8,9,10,11,12,13,14}
+#endif /* For OSX */
+#define SFMT_IDSTR "SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6"
+
+#endif /* SFMT_PARAMS19937_H */
View
121 lib/3rdparty/sfmt/SFMT-sse2.h
@@ -0,0 +1,121 @@
+#pragma once
+/**
+ * @file SFMT-sse2.h
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) for Intel SSE2
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * @note We assume LITTLE ENDIAN in this file
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University. All rights reserved.
+ *
+ * The new BSD License is applied to this software, see LICENSE.txt
+ */
+
+#ifndef SFMT_SSE2_H
+#define SFMT_SSE2_H
+
+inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
+ __m128i c, __m128i d);
+
+/**
+ * This function represents the recursion formula.
+ * @param r an output
+ * @param a a 128-bit part of the interal state array
+ * @param b a 128-bit part of the interal state array
+ * @param c a 128-bit part of the interal state array
+ * @param d a 128-bit part of the interal state array
+ */
+inline static void mm_recursion(__m128i * r, __m128i a, __m128i b,
+ __m128i c, __m128i d)
+{
+ __m128i v, x, y, z;
+
+ y = _mm_srli_epi32(b, SFMT_SR1);
+ z = _mm_srli_si128(c, SFMT_SR2);
+ v = _mm_slli_epi32(d, SFMT_SL1);
+ z = _mm_xor_si128(z, a);
+ z = _mm_xor_si128(z, v);
+ x = _mm_slli_si128(a, SFMT_SL2);
+ y = _mm_and_si128(y, sse2_param_mask.si);
+ z = _mm_xor_si128(z, x);
+ z = _mm_xor_si128(z, y);
+ *r = z;
+}
+
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+ int i;
+ __m128i r1, r2;
+ w128_t * pstate = sfmt->state;
+
+ r1 = pstate[SFMT_N - 2].si;
+ r2 = pstate[SFMT_N - 1].si;
+ for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+ mm_recursion(&pstate[i].si, pstate[i].si,
+ pstate[i + SFMT_POS1].si, r1, r2);
+ r1 = r2;
+ r2 = pstate[i].si;
+ }
+ for (; i < SFMT_N; i++) {
+ mm_recursion(&pstate[i].si, pstate[i].si,
+ pstate[i + SFMT_POS1 - SFMT_N].si,
+ r1, r2);
+ r1 = r2;
+ r2 = pstate[i].si;
+ }
+}
+
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state.
+ * @param array an 128-bit array to be filled by pseudorandom numbers.
+ * @param size number of 128-bit pseudorandom numbers to be generated.
+ */
+static void gen_rand_array(sfmt_t * sfmt, w128_t * array, int size)
+{
+ int i, j;
+ __m128i r1, r2;
+ w128_t * pstate = sfmt->state;
+
+ r1 = pstate[SFMT_N - 2].si;
+ r2 = pstate[SFMT_N - 1].si;
+ for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+ mm_recursion(&array[i].si, pstate[i].si,
+ pstate[i + SFMT_POS1].si, r1, r2);
+ r1 = r2;
+ r2 = array[i].si;
+ }
+ for (; i < SFMT_N; i++) {
+ mm_recursion(&array[i].si, pstate[i].si,
+ array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+ r1 = r2;
+ r2 = array[i].si;
+ }
+ for (; i < size - SFMT_N; i++) {
+ mm_recursion(&array[i].si, array[i - SFMT_N].si,
+ array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+ r1 = r2;
+ r2 = array[i].si;
+ }
+ for (j = 0; j < 2 * SFMT_N - size; j++) {
+ pstate[j] = array[j + size - SFMT_N];
+ }
+ for (; i < size; i++, j++) {
+ mm_recursion(&array[i].si, array[i - SFMT_N].si,
+ array[i + SFMT_POS1 - SFMT_N].si, r1, r2);
+ r1 = r2;
+ r2 = array[i].si;
+ pstate[j] = array[i];
+ }
+}
+
+
+#endif
View
427 lib/3rdparty/sfmt/SFMT.c
@@ -0,0 +1,427 @@
+/**
+ * @file SFMT.c
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT)
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (Hiroshima University)
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University.
+ * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
+ * University and The University of Tokyo.
+ * All rights reserved.
+ *
+ * The 3-clause BSD License is applied to this software, see
+ * LICENSE.txt
+ */
+
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+#include <string.h>
+#include <assert.h>
+#include "SFMT.h"
+#include "SFMT-params.h"
+#include "SFMT-common.h"
+
+#if defined(__BIG_ENDIAN__) && !defined(__amd64) && !defined(BIG_ENDIAN64)
+#define BIG_ENDIAN64 1
+#endif
+#if defined(HAVE_ALTIVEC) && !defined(BIG_ENDIAN64)
+#define BIG_ENDIAN64 1
+#endif
+#if defined(ONLY64) && !defined(BIG_ENDIAN64)
+ #if defined(__GNUC__)
+ #error "-DONLY64 must be specified with -DBIG_ENDIAN64"
+ #endif
+#undef ONLY64
+#endif
+
+/**
+ * parameters used by sse2.
+ */
+static const w128_t sse2_param_mask = {{SFMT_MSK1, SFMT_MSK2,
+ SFMT_MSK3, SFMT_MSK4}};
+/*----------------
+ STATIC FUNCTIONS
+ ----------------*/
+inline static int idxof(int i);
+inline static void gen_rand_array(sfmt_t * sfmt, w128_t *array, int size);
+inline static uint32_t func1(uint32_t x);
+inline static uint32_t func2(uint32_t x);
+static void period_certification(sfmt_t * sfmt);
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+inline static void swap(w128_t *array, int size);
+#endif
+
+#if defined(HAVE_ALTIVEC)
+ #include "SFMT-alti.h"
+#elif defined(HAVE_SSE2)
+ #include "SFMT-sse2.h"
+#endif
+
+/**
+ * This function simulate a 64-bit index of LITTLE ENDIAN
+ * in BIG ENDIAN machine.
+ */
+#ifdef ONLY64
+inline static int idxof(int i) {
+ return i ^ 1;
+}
+#else
+inline static int idxof(int i) {
+ return i;
+}
+#endif
+
+#if (!defined(HAVE_ALTIVEC)) && (!defined(HAVE_SSE2))
+/**
+ * This function fills the user-specified array with pseudorandom
+ * integers.
+ *
+ * @param sfmt SFMT internal state
+ * @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(sfmt_t * sfmt, w128_t *array, int size) {
+ int i, j;
+ w128_t *r1, *r2;
+
+ r1 = &sfmt->state[SFMT_N - 2];
+ r2 = &sfmt->state[SFMT_N - 1];
+ for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+ do_recursion(&array[i], &sfmt->state[i], &sfmt->state[i + SFMT_POS1], r1, r2);
+ r1 = r2;
+ r2 = &array[i];
+ }
+ for (; i < SFMT_N; i++) {
+ do_recursion(&array[i], &sfmt->state[i],
+ &array[i + SFMT_POS1 - SFMT_N], r1, r2);
+ r1 = r2;
+ r2 = &array[i];
+ }
+ for (; i < size - SFMT_N; i++) {
+ do_recursion(&array[i], &array[i - SFMT_N],
+ &array[i + SFMT_POS1 - SFMT_N], r1, r2);
+ r1 = r2;
+ r2 = &array[i];
+ }
+ for (j = 0; j < 2 * SFMT_N - size; j++) {
+ sfmt->state[j] = array[j + size - SFMT_N];
+ }
+ for (; i < size; i++, j++) {
+ do_recursion(&array[i], &array[i - SFMT_N],
+ &array[i + SFMT_POS1 - SFMT_N], r1, r2);
+ r1 = r2;
+ r2 = &array[i];
+ sfmt->state[j] = array[i];
+ }
+}
+#endif
+
+#if defined(BIG_ENDIAN64) && !defined(ONLY64) && !defined(HAVE_ALTIVEC)
+inline static void swap(w128_t *array, int size) {
+ int i;
+ uint32_t x, y;
+
+ for (i = 0; i < size; i++) {
+ x = array[i].u[0];
+ y = array[i].u[2];
+ array[i].u[0] = array[i].u[1];
+ array[i].u[2] = array[i].u[3];
+ array[i].u[1] = x;
+ array[i].u[3] = y;
+ }
+}
+#endif
+/**
+ * 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 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 func2(uint32_t x) {
+ return (x ^ (x >> 27)) * (uint32_t)1566083941UL;
+}
+
+/**
+ * This function certificate the period of 2^{MEXP}
+ * @param sfmt SFMT internal state
+ */
+static void period_certification(sfmt_t * sfmt) {
+ int inner = 0;
+ int i, j;
+ uint32_t work;
+ uint32_t *psfmt32 = &sfmt->state[0].u[0];
+ const uint32_t parity[4] = {SFMT_PARITY1, SFMT_PARITY2,
+ SFMT_PARITY3, SFMT_PARITY4};
+
+ for (i = 0; i < 4; i++)
+ inner ^= psfmt32[idxof(i)] & parity[i];
+ for (i = 16; i > 0; i >>= 1)
+ inner ^= inner >> i;
+ inner &= 1;
+ /* check OK */
+ if (inner == 1) {
+ return;
+ }
+ /* check NG, and modification */
+ for (i = 0; i < 4; i++) {
+ work = 1;
+ for (j = 0; j < 32; j++) {
+ if ((work & parity[i]) != 0) {
+ psfmt32[idxof(i)] ^= work;
+ return;
+ }
+ work = work << 1;
+ }
+ }
+}
+
+/*----------------
+ PUBLIC FUNCTIONS
+ ----------------*/
+#define UNUSED_VARIABLE(x) (void)(x)
+/**
+ * This function returns the identification string.
+ * The string shows the word size, the Mersenne exponent,
+ * and all parameters of this generator.
+ * @param sfmt SFMT internal state
+ */
+const char *sfmt_get_idstring(sfmt_t * sfmt) {
+ UNUSED_VARIABLE(sfmt);
+ return SFMT_IDSTR;
+}
+
+/**
+ * This function returns the minimum size of array used for \b
+ * fill_array32() function.
+ * @param sfmt SFMT internal state
+ * @return minimum size of array used for fill_array32() function.
+ */
+int sfmt_get_min_array_size32(sfmt_t * sfmt) {
+ UNUSED_VARIABLE(sfmt);
+ return SFMT_N32;
+}
+
+/**
+ * This function returns the minimum size of array used for \b
+ * fill_array64() function.
+ * @param sfmt SFMT internal state
+ * @return minimum size of array used for fill_array64() function.
+ */
+int sfmt_get_min_array_size64(sfmt_t * sfmt) {
+ UNUSED_VARIABLE(sfmt);
+ return SFMT_N64;
+}
+
+#if !defined(HAVE_SSE2) && !defined(HAVE_ALTIVEC)
+/**
+ * This function fills the internal state array with pseudorandom
+ * integers.
+ * @param sfmt SFMT internal state
+ */
+void sfmt_gen_rand_all(sfmt_t * sfmt) {
+ int i;
+ w128_t *r1, *r2;
+
+ r1 = &sfmt->state[SFMT_N - 2];
+ r2 = &sfmt->state[SFMT_N - 1];
+ for (i = 0; i < SFMT_N - SFMT_POS1; i++) {
+ do_recursion(&sfmt->state[i], &sfmt->state[i],
+ &sfmt->state[i + SFMT_POS1], r1, r2);
+ r1 = r2;
+ r2 = &sfmt->state[i];
+ }
+ for (; i < SFMT_N; i++) {
+ do_recursion(&sfmt->state[i], &sfmt->state[i],
+ &sfmt->state[i + SFMT_POS1 - SFMT_N], r1, r2);
+ r1 = r2;
+ r2 = &sfmt->state[i];
+ }
+}
+#endif
+
+#ifndef ONLY64
+/**
+ * This function generates pseudorandom 32-bit integers in the
+ * specified array[] by one call. The number of pseudorandom integers
+ * is specified by the argument size, which must be at least 624 and a
+ * multiple of four. The generation by this function is much faster
+ * than the following gen_rand function.
+ *
+ * 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 gen_rand function, without initialization.
+ *
+ * @param sfmt SFMT internal state
+ * @param array an array where pseudorandom 32-bit integers are filled
+ * by this function. The pointer to the array must be \b "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 32-bit pseudorandom integers to be
+ * generated. size must be a multiple of 4, and greater than or equal
+ * to (MEXP / 128 + 1) * 4.
+ *
+ * @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 sfmt_fill_array32(sfmt_t * sfmt, uint32_t *array, int size) {
+ assert(sfmt->idx == SFMT_N32);
+ assert(size % 4 == 0);
+ assert(size >= SFMT_N32);
+
+ gen_rand_array(sfmt, (w128_t *)array, size / 4);
+ sfmt->idx = SFMT_N32;
+}
+#endif
+
+/**
+ * This function generates pseudorandom 64-bit integers in the
+ * specified array[] by one call. The number of pseudorandom integers
+ * is specified by the argument size, which must be at least 312 and a
+ * multiple of two. The generation by this function is much faster
+ * than the following gen_rand function.
+ *
+ * @param sfmt SFMT internal state
+ * 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 gen_rand function, without initialization.
+ *
+ * @param array an array where pseudorandom 64-bit integers 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 (MEXP / 128 + 1) * 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 sfmt_fill_array64(sfmt_t * sfmt, uint64_t *array, int size) {
+ assert(sfmt->idx == SFMT_N32);
+ assert(size % 2 == 0);
+ assert(size >= SFMT_N64);
+
+ gen_rand_array(sfmt, (w128_t *)array, size / 2);
+ sfmt->idx = SFMT_N32;
+
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+ swap((w128_t *)array, size /2);
+#endif
+}
+
+/**
+ * This function initializes the internal state array with a 32-bit
+ * integer seed.
+ *
+ * @param sfmt SFMT internal state
+ * @param seed a 32-bit integer used as the seed.
+ */
+void sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed) {
+ int i;
+
+ uint32_t *psfmt32 = &sfmt->state[0].u[0];
+
+ psfmt32[idxof(0)] = seed;
+ for (i = 1; i < SFMT_N32; i++) {
+ psfmt32[idxof(i)] = 1812433253UL * (psfmt32[idxof(i - 1)]
+ ^ (psfmt32[idxof(i - 1)] >> 30))
+ + i;
+ }
+ sfmt->idx = SFMT_N32;
+ period_certification(sfmt);
+}
+
+/**
+ * This function initializes the internal state array,
+ * with an array of 32-bit integers used as the seeds
+ * @param sfmt SFMT internal state
+ * @param init_key the array of 32-bit integers, used as a seed.
+ * @param key_length the length of init_key.
+ */
+void sfmt_init_by_array(sfmt_t * sfmt, uint32_t *init_key, int key_length) {
+ int i, j, count;
+ uint32_t r;
+ int lag;
+ int mid;
+ int size = SFMT_N * 4;
+ uint32_t *psfmt32 = &sfmt->state[0].u[0];
+
+ if (size >= 623) {
+ lag = 11;
+ } else if (size >= 68) {
+ lag = 7;
+ } else if (size >= 39) {
+ lag = 5;
+ } else {
+ lag = 3;
+ }
+ mid = (size - lag) / 2;
+
+ memset(sfmt, 0x8b, sizeof(sfmt_t));
+ if (key_length + 1 > SFMT_N32) {
+ count = key_length + 1;
+ } else {
+ count = SFMT_N32;
+ }
+ r = func1(psfmt32[idxof(0)] ^ psfmt32[idxof(mid)]
+ ^ psfmt32[idxof(SFMT_N32 - 1)]);
+ psfmt32[idxof(mid)] += r;
+ r += key_length;
+ psfmt32[idxof(mid + lag)] += r;
+ psfmt32[idxof(0)] = r;
+
+ count--;
+ for (i = 1, j = 0; (j < count) && (j < key_length); j++) {
+ r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % SFMT_N32)]
+ ^ psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
+ psfmt32[idxof((i + mid) % SFMT_N32)] += r;
+ r += init_key[j] + i;
+ psfmt32[idxof((i + mid + lag) % SFMT_N32)] += r;
+ psfmt32[idxof(i)] = r;
+ i = (i + 1) % SFMT_N32;
+ }
+ for (; j < count; j++) {
+ r = func1(psfmt32[idxof(i)] ^ psfmt32[idxof((i + mid) % SFMT_N32)]
+ ^ psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
+ psfmt32[idxof((i + mid) % SFMT_N32)] += r;
+ r += i;
+ psfmt32[idxof((i + mid + lag) % SFMT_N32)] += r;
+ psfmt32[idxof(i)] = r;
+ i = (i + 1) % SFMT_N32;
+ }
+ for (j = 0; j < SFMT_N32; j++) {
+ r = func2(psfmt32[idxof(i)] + psfmt32[idxof((i + mid) % SFMT_N32)]
+ + psfmt32[idxof((i + SFMT_N32 - 1) % SFMT_N32)]);
+ psfmt32[idxof((i + mid) % SFMT_N32)] ^= r;
+ r -= i;
+ psfmt32[idxof((i + mid + lag) % SFMT_N32)] ^= r;
+ psfmt32[idxof(i)] = r;
+ i = (i + 1) % SFMT_N32;
+ }
+
+ sfmt->idx = SFMT_N32;
+ period_certification(sfmt);
+}
+#if defined(__cplusplus)
+}
+#endif
View
295 lib/3rdparty/sfmt/SFMT.h
@@ -0,0 +1,295 @@
+#pragma once
+/**
+ * @file SFMT.h
+ *
+ * @brief SIMD oriented Fast Mersenne Twister(SFMT) pseudorandom
+ * number generator using C structure.
+ *
+ * @author Mutsuo Saito (Hiroshima University)
+ * @author Makoto Matsumoto (The University of Tokyo)
+ *
+ * Copyright (C) 2006, 2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima
+ * University.
+ * Copyright (C) 2012 Mutsuo Saito, Makoto Matsumoto, Hiroshima
+ * University and The University of Tokyo.
+ * All rights reserved.
+ *
+ * The 3-clause BSD License is applied to this software, see
+ * LICENSE.txt
+ *
+ * @note We assume that your system has inttypes.h. If your system
+ * doesn't have inttypes.h, you have to typedef uint32_t and uint64_t,
+ * and you have to define PRIu64 and PRIx64 in this file as follows:
+ * @verbatim
+ typedef unsigned int uint32_t
+ typedef unsigned long long uint64_t
+ #define PRIu64 "llu"
+ #define PRIx64 "llx"
+@endverbatim
+ * uint32_t must be exactly 32-bit unsigned integer type (no more, no
+ * less), and uint64_t must be exactly 64-bit unsigned integer type.
+ * PRIu64 and PRIx64 are used for printf function to print 64-bit
+ * unsigned int and 64-bit unsigned int in hexadecimal format.
+ */
+
+#ifndef SFMTST_H
+#define SFMTST_H
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+#include <stdio.h>
+#include <assert.h>
+
+#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
+ #include <inttypes.h>
+#elif defined(_MSC_VER) || defined(__BORLANDC__)
+ typedef unsigned int uint32_t;
+ typedef unsigned __int64 uint64_t;
+ #define inline __inline
+#else
+ #include <inttypes.h>
+ #if defined(__GNUC__)
+ #define inline __inline__
+ #endif
+#endif
+
+#ifndef PRIu64
+ #if defined(_MSC_VER) || defined(__BORLANDC__)
+ #define PRIu64 "I64u"
+ #define PRIx64 "I64x"
+ #else
+ #define PRIu64 "llu"
+ #define PRIx64 "llx"
+ #endif
+#endif
+
+#include "SFMT-params.h"
+
+/*------------------------------------------
+ 128-bit SIMD like data type for standard C
+ ------------------------------------------*/
+#if defined(HAVE_ALTIVEC)
+ #if !defined(__APPLE__)
+ #include <altivec.h>
+ #endif
+/** 128-bit data structure */
+union W128_T {
+ vector unsigned int s;
+ uint32_t u[4];
+ uint64_t u64[2];
+};
+#elif defined(HAVE_SSE2)
+ #include <emmintrin.h>
+
+/** 128-bit data structure */
+union W128_T {
+ uint32_t u[4];
+ uint64_t u64[2];
+ __m128i si;
+};
+#else
+/** 128-bit data structure */
+union W128_T {
+ uint32_t u[4];
+ uint64_t u64[2];
+};
+#endif
+
+/** 128-bit data type */
+typedef union W128_T w128_t;
+
+/**
+ * SFMT internal state
+ */
+struct SFMT_T {
+ /** the 128-bit internal state array */
+ w128_t state[SFMT_N];
+ /** index counter to the 32-bit internal state array */
+ int idx;
+};
+
+typedef struct SFMT_T sfmt_t;
+
+void sfmt_fill_array32(sfmt_t * sfmt, uint32_t * array, int size);
+void sfmt_fill_array64(sfmt_t * sfmt, uint64_t * array, int size);
+void sfmt_init_gen_rand(sfmt_t * sfmt, uint32_t seed);
+void sfmt_init_by_array(sfmt_t * sfmt, uint32_t * init_key, int key_length);
+const char * sfmt_get_idstring(sfmt_t * sfmt);
+int sfmt_get_min_array_size32(sfmt_t * sfmt);
+int sfmt_get_min_array_size64(sfmt_t * sfmt);
+void sfmt_gen_rand_all(sfmt_t * sfmt);
+
+#ifndef ONLY64
+/**
+ * This function generates and returns 32-bit pseudorandom number.
+ * init_gen_rand or init_by_array must be called before this function.
+ * @param sfmt SFMT internal state
+ * @return 32-bit pseudorandom number
+ */
+inline static uint32_t sfmt_genrand_uint32(sfmt_t * sfmt) {
+ uint32_t r;
+ uint32_t * psfmt32 = &sfmt->state[0].u[0];
+
+ if (sfmt->idx >= SFMT_N32) {
+ sfmt_gen_rand_all(sfmt);
+ sfmt->idx = 0;
+ }
+ r = psfmt32[sfmt->idx++];
+ return r;
+}
+#endif
+/**
+ * This function generates and returns 64-bit pseudorandom number.
+ * init_gen_rand or init_by_array must be called before this function.
+ * The function gen_rand64 should not be called after gen_rand32,
+ * unless an initialization is again executed.
+ * @param sfmt SFMT internal state
+ * @return 64-bit pseudorandom number
+ */
+inline static uint64_t sfmt_genrand_uint64(sfmt_t * sfmt) {
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+ uint32_t * psfmt32 = &sfmt->state[0].u[0];
+ uint32_t r1, r2;
+#else
+ uint64_t r;
+#endif
+ uint64_t * psfmt64 = &sfmt->state[0].u64[0];
+ assert(sfmt->idx % 2 == 0);
+
+ if (sfmt->idx >= SFMT_N32) {
+ sfmt_gen_rand_all(sfmt);
+ sfmt->idx = 0;
+ }
+#if defined(BIG_ENDIAN64) && !defined(ONLY64)
+ r1 = psfmt32[sfmt->idx];
+ r2 = psfmt32[sfmt->idx + 1];
+ sfmt->idx += 2;
+ return ((uint64_t)r2 << 32) | r1;
+#else
+ r = psfmt64[sfmt->idx / 2];
+ sfmt->idx += 2;
+ return r;
+#endif
+}
+
+/* =================================================
+ The following real versions are due to Isaku Wada
+ ================================================= */
+/**
+ * converts an unsigned 32-bit number to a double on [0,1]-real-interval.
+ * @param v 32-bit unsigned integer
+ * @return double on [0,1]-real-interval
+ */
+inline static double sfmt_to_real1(uint32_t v)
+{
+ return v * (1.0/4294967295.0);
+ /* divided by 2^32-1 */
+}
+
+/**
+ * generates a random number on [0,1]-real-interval
+ * @param sfmt SFMT internal state
+ * @return double on [0,1]-real-interval
+ */
+inline static double sfmt_genrand_real1(sfmt_t * sfmt)
+{
+ return sfmt_to_real1(sfmt_genrand_uint32(sfmt));
+}
+
+/**
+ * converts an unsigned 32-bit integer to a double on [0,1)-real-interval.
+ * @param v 32-bit unsigned integer
+ * @return double on [0,1)-real-interval
+ */
+inline static double sfmt_to_real2(uint32_t v)
+{
+ return v * (1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+
+/**
+ * generates a random number on [0,1)-real-interval
+ * @param sfmt SFMT internal state
+ * @return double on [0,1)-real-interval
+ */
+inline static double sfmt_genrand_real2(sfmt_t * sfmt)
+{
+ return sfmt_to_real2(sfmt_genrand_uint32(sfmt));
+}
+
+/**
+ * converts an unsigned 32-bit integer to a double on (0,1)-real-interval.
+ * @param v 32-bit unsigned integer
+ * @return double on (0,1)-real-interval
+ */
+inline static double sfmt_to_real3(uint32_t v)
+{
+ return (((double)v) + 0.5)*(1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+
+/**
+ * generates a random number on (0,1)-real-interval
+ * @param sfmt SFMT internal state
+ * @return double on (0,1)-real-interval
+ */
+inline static double sfmt_genrand_real3(sfmt_t * sfmt)
+{
+ return sfmt_to_real3(sfmt_genrand_uint32(sfmt));
+}
+
+/**
+ * converts an unsigned 32-bit integer to double on [0,1)
+ * with 53-bit resolution.
+ * @param v 32-bit unsigned integer
+ * @return double on [0,1)-real-interval with 53-bit resolution.
+ */
+inline static double sfmt_to_res53(uint64_t v)
+{
+ return v * (1.0/18446744073709551616.0L);
+}
+
+/**
+ * generates a random number on [0,1) with 53-bit resolution
+ * @param sfmt SFMT internal state
+ * @return double on [0,1) with 53-bit resolution
+ */
+inline static double sfmt_genrand_res53(sfmt_t * sfmt)
+{
+ return sfmt_to_res53(sfmt_genrand_uint64(sfmt));
+}
+
+
+/* =================================================
+ The following function are added by Saito.
+ ================================================= */
+/**
+ * generates a random number on [0,1) with 53-bit resolution from two
+ * 32 bit integers
+ */
+inline static double sfmt_to_res53_mix(uint32_t x, uint32_t y)
+{
+ return sfmt_to_res53(x | ((uint64_t)y << 32));
+}
+
+/**
+ * generates a random number on [0,1) with 53-bit resolution
+ * using two 32bit integers.
+ * @param sfmt SFMT internal state
+ * @return double on [0,1) with 53-bit resolution
+ */
+inline static double sfmt_genrand_res53_mix(sfmt_t * sfmt)
+{
+ uint32_t x, y;
+
+ x = sfmt_genrand_uint32(sfmt);
+ y = sfmt_genrand_uint32(sfmt);
+ return sfmt_to_res53_mix(x, y);
+}
+
+#if defined(__cplusplus)
+}
+#endif
+
+#endif
View
17 lib/ccv.h
@@ -825,15 +825,21 @@ void ccv_bbf_classifier_cascade_free(ccv_bbf_classifier_cascade_t* cascade);
* see: http://cvlab.epfl.ch/alumni/oezuysal/ferns.html for more about ferns */
typedef struct {
- int ferns;
+ int structs;
int features;
int scales;
+ int posteriors;
+ float threshold;
+ int* cnum;
+ float* posterior;
// decided to go flat organized fern so that we can profiling different memory layout impacts the performance
ccv_point_t fern[1];
} ccv_ferns_t;
-ccv_ferns_t* __attribute__((warn_unused_result)) ccv_ferns_new(int nferns, int features, int scale, ccv_size_t* sizes);
-float ccv_ferns_predict(ccv_ferns_t* ferns, ccv_dense_matrix_t* a);
+ccv_ferns_t* __attribute__((warn_unused_result)) ccv_ferns_new(int structs, int features, int scales, ccv_size_t* sizes);
+void ccv_ferns_feature(ccv_ferns_t* ferns, ccv_dense_matrix_t* a, int scale, uint32_t* fern);
+void ccv_ferns_correct(ccv_ferns_t* ferns, uint32_t* fern, int c, int repeat);
+float ccv_ferns_predict(ccv_ferns_t* ferns, uint32_t* fern);
void ccv_ferns_free(ccv_ferns_t* ferns);
/* TLD: Track-Learn-Detection is a long-term object tracking framework, which achieved very high
@@ -845,11 +851,16 @@ typedef struct {
int level;
double min_eigen;
double min_forward_backward_error;
+ /* samples generation parameters */
+ int min_win;
+ double include_overlap;
+ double exclude_overlap;
} ccv_tld_param_t;
typedef struct {
ccv_tld_param_t params;
ccv_comp_t box;
+ ccv_ferns_t* ferns;
} ccv_tld_t;
ccv_tld_t* __attribute__((warn_unused_result)) ccv_tld_new(ccv_dense_matrix_t* a, ccv_rect_t box, ccv_tld_param_t params);
View
83 lib/ccv_ferns.c
@@ -2,16 +2,23 @@
#include "ccv_internal.h"
#include "3rdparty/dsfmt/dSFMT.h"
-ccv_ferns_t* ccv_ferns_new(int nferns, int features, int scale, ccv_size_t* sizes)
+ccv_ferns_t* ccv_ferns_new(int structs, int features, int scales, ccv_size_t* sizes)
{
- ccv_ferns_t* ferns = (ccv_ferns_t*)ccmalloc(sizeof(ccv_ferns_t) + sizeof(ccv_point_t) * (nferns * features * scale * 2 - 1));
- ferns->ferns = nferns;
+ assert(structs > 0 && features > 0 && scales > 0);
+ int posteriors = (int)(powf(2.0, features) + 0.5);
+ ccv_ferns_t* ferns = (ccv_ferns_t*)ccmalloc(sizeof(ccv_ferns_t) + sizeof(ccv_point_t) * (structs * features * scales * 2 - 1) + sizeof(float) * structs * posteriors + sizeof(int) * structs * posteriors * 2);
+ ferns->structs = structs;
ferns->features = features;
- ferns->scales = scale;
+ ferns->scales = scales;
+ ferns->posteriors = posteriors;
+ ferns->posterior = (float*)((uint8_t*)(ferns + 1) + sizeof(ccv_point_t) * (structs * features * scales * 2 - 1));
+ // now only for 2 classes
+ ferns->cnum = (int*)(ferns->posterior + structs * posteriors);
+ memset(ferns->posterior, 0, sizeof(float) * structs * posteriors + sizeof(int) * structs * posteriors * 2);
int i, j, k;
dsfmt_t dsfmt;
dsfmt_init_gen_rand(&dsfmt, (uint32_t)ferns);
- for (i = 0; i < nferns; i++)
+ for (i = 0; i < structs; i++)
{
for (k = 0; k < features; k++)
{
@@ -19,16 +26,76 @@ ccv_ferns_t* ccv_ferns_new(int nferns, int features, int scale, ccv_size_t* size
double y1f = dsfmt_genrand_close_open(&dsfmt);
double x2f = dsfmt_genrand_close_open(&dsfmt);
double y2f = dsfmt_genrand_close_open(&dsfmt);
- for (j = 0; j < scale; j++)
+ for (j = 0; j < scales; j++)
{
- ferns->fern[(j * nferns * features + i * features + k) * 2] = ccv_point((int)(x1f * sizes[j].width), (int)(y1f * sizes[j].height));
- ferns->fern[(j * nferns * features + i * features + k) * 2 + 1] = ccv_point((int)(x2f * sizes[j].width), (int)(y2f * sizes[j].height));
+ ferns->fern[(j * structs * features + i * features + k) * 2] = ccv_point((int)(x1f * sizes[j].width), (int)(y1f * sizes[j].height));
+ ferns->fern[(j * structs * features + i * features + k) * 2 + 1] = ccv_point((int)(x2f * sizes[j].width), (int)(y2f * sizes[j].height));
}
}
}
+ ferns->threshold = 0.5 * structs;
return ferns;
}
+void ccv_ferns_feature(ccv_ferns_t* ferns, ccv_dense_matrix_t* a, int scale, uint32_t* fern)
+{
+ ccv_point_t* fern_feature = ferns->fern + scale * ferns->structs * ferns->features;
+ int i, j;
+ unsigned char* a_ptr = a->data.u8;
+ assert(CCV_GET_CHANNEL(a->type) == CCV_C1);
+#define for_block(_, _for_get) \
+ for (i = 0; i < ferns->structs; i++) \
+ { \
+ uint32_t leaf = 0; \
+ for (j = 0; j < ferns->features; j++) \
+ { \
+ if (_for_get(a_ptr + fern_feature[0].y * a->step, fern_feature[0].x, 0) > _for_get(a_ptr + fern_feature[0].y * a->step, fern_feature[0].x, 0)) \
+ leaf = (leaf << 1) | 1; \
+ else \
+ leaf = leaf << 1; \
+ fern_feature += 2; \
+ } \
+ fern[i] = leaf; \
+ }
+ ccv_matrix_getter(a->type, for_block);
+#undef for_block
+}
+
+void ccv_ferns_correct(ccv_ferns_t* ferns, uint32_t* fern, int c, int repeat)
+{
+ assert(c == 0 || c == 1);
+ assert(repeat >= 0);
+ repeat += 1;
+ int i;
+ int* cnum = ferns->cnum;
+ float* post = ferns->posterior;
+ for (i = 0; i < ferns->structs; i++)
+ {
+ uint32_t k = fern[i];
+ cnum[k * 2 + c] += repeat;
+ if (cnum[k * 2] == 0)
+ post[k] = 0;
+ else
+ // needs to compute the log of it
+ post[k] = logf((float)cnum[k * 2] / (cnum[k * 2] + cnum[k * 2 + 1]));
+ cnum += ferns->posteriors * 2;
+ post += ferns->posteriors;
+ }
+}
+
+float ccv_ferns_predict(ccv_ferns_t* ferns, uint32_t* fern)
+{
+ float votes = 0;
+ int i;
+ float* post = ferns->posterior;
+ for (i = 0; i < ferns->structs; i++)
+ {
+ votes += post[fern[i]];
+ post += ferns->posteriors;
+ }
+ return votes;
+}
+
void ccv_ferns_free(ccv_ferns_t* ferns)
{
ccfree(ferns);
View
65 lib/ccv_tld.c
@@ -172,8 +172,64 @@ static ccv_rect_t _ccv_tld_short_term_track(ccv_dense_matrix_t* a, ccv_dense_mat
return newbox;
}
-void _ccv_tld_box(ccv_rect_t box, ccv_array_t** good, ccv_array_t** bad, double include_overlap, double exclude_overlap)
+ccv_comp_t _ccv_tld_generate_box_for(ccv_size_t image_size, ccv_rect_t box, ccv_array_t** scale, ccv_array_t** good, ccv_array_t** bad, ccv_tld_param_t params)
{
+ const float SHIFT = 0.1;
+ const float SCALES[] = {
+ 0.16151, 0.19381, 0.23257, 0.27908, 0.33490, 0.40188, 0.48225,
+ 0.57870, 0.69444, 0.83333, 1, 1.20000, 1.44000, 1.72800,
+ 2.07360, 2.48832, 2.98598, 3.58318, 4.29982, 5.15978, 6.19174
+ };
+ ccv_array_t* ascale = *scale = ccv_array_new(sizeof(ccv_size_t), 21, 0);
+ ccv_array_t* agood = *good = ccv_array_new(sizeof(ccv_comp_t), 64, 0);
+ ccv_array_t* abad = *bad = ccv_array_new(sizeof(ccv_comp_t), 64, 0);
+ int s;
+ double max_overlap = -DBL_MAX;
+ ccv_comp_t best_box = {
+ .id = 0,
+ .rect = ccv_rect(0, 0, 0, 0),
+ };
+ for (s = 0; s < 21; s++)
+ {
+ int width = (int)(box.width * SCALES[s] + 0.5);
+ int height = (int)(box.height * SCALES[s] + 0.5);
+ int min_side = ccv_min(width, height);
+ if (min_side < params.min_win || width > image_size.width || height > image_size.height)
+ continue;
+ ccv_size_t scale = ccv_size(width, height);
+ ccv_array_push(ascale, &scale);
+ for (float y = 0; y < image_size.height - height; y += SHIFT * min_side)
+ {
+ for (float x = 0; x < image_size.width - width; x += SHIFT * min_side)
+ {
+ ccv_comp_t comp;
+ comp.rect = ccv_rect((int)(x + 0.5), (int)(y + 0.5), width, height);
+ comp.id = ascale->rnum - 1;
+ double overlap = ccv_max(ccv_min(comp.rect.x + comp.rect.width, box.x + box.width) - ccv_max(comp.rect.x, box.x), 0) *
+ ccv_max(ccv_min(comp.rect.y + comp.rect.height, box.y + box.height) - ccv_max(comp.rect.y, box.y), 0);
+ overlap = overlap / (comp.rect.width * comp.rect.height + box.width * box.height - overlap);
+ if (overlap > params.include_overlap)
+ {
+ if (overlap > max_overlap)
+ {
+ max_overlap = overlap;
+ best_box = comp;
+ }
+ ccv_array_push(agood, &comp);
+ } else if (overlap < params.exclude_overlap) {
+ ccv_array_push(abad, &comp);
+ }
+ }
+ }
+ }
+ return best_box;
+}
+
+void _ccv_tld_ferns_feature_for(ccv_ferns_t* ferns, ccv_dense_matrix_t* a, ccv_comp_t box, uint32_t* fern)
+{
+ ccv_dense_matrix_t roi = ccv_dense_matrix(box.rect.height, box.rect.width, CCV_GET_DATA_TYPE(a->type) | CCV_GET_CHANNEL(a->type), ccv_get_dense_matrix_cell(a, box.rect.y, box.rect.x, 0), 0);
+ roi.step = a->step;
+ ccv_ferns_feature(ferns, &roi, box.id, fern);
}
ccv_tld_t* ccv_tld_new(ccv_dense_matrix_t* a, ccv_rect_t box, ccv_tld_param_t params)
@@ -181,6 +237,12 @@ ccv_tld_t* ccv_tld_new(ccv_dense_matrix_t* a, ccv_rect_t box, ccv_tld_param_t pa
ccv_tld_t* tld = (ccv_tld_t*)ccmalloc(sizeof(ccv_tld_t));
tld->params = params;
tld->box.rect = box;
+ ccv_size_t image_size = ccv_size(a->cols, a->rows);
+ ccv_array_t* scales = 0;
+ ccv_array_t* good = 0;
+ ccv_array_t* bad = 0;
+ _ccv_tld_generate_box_for(image_size, box, &scales, &good, &bad, params);
+ tld->ferns = ccv_ferns_new(1, 1, scales->rnum, (ccv_size_t*)ccv_array_get(scales, 0));
return tld;
}
@@ -195,5 +257,6 @@ ccv_comp_t ccv_tld_track_object(ccv_tld_t* tld, ccv_dense_matrix_t* a, ccv_dense
void ccv_tld_free(ccv_tld_t* tld)
{
+ ccv_ferns_free(tld->ferns);
ccfree(tld);
}
View
2  lib/makefile
@@ -14,7 +14,7 @@ all: libccv.a
clean:
rm -f *.o 3rdparty/sha1/*.o 3rdparty/kissfft/*.o 3rdparty/dsfmt/*.o libccv.a
-libccv.a: ccv_cache.o ccv_memory.o 3rdparty/sha1/sha1.o 3rdparty/kissfft/kiss_fft.o 3rdparty/kissfft/kiss_fftnd.o 3rdparty/kissfft/kiss_fftr.o 3rdparty/kissfft/kiss_fftndr.o 3rdparty/kissfft/kissf_fft.o 3rdparty/kissfft/kissf_fftnd.o 3rdparty/kissfft/kissf_fftr.o 3rdparty/kissfft/kissf_fftndr.o 3rdparty/dsfmt/dSFMT.o ccv_io.o ccv_numeric.o ccv_algebra.o ccv_util.o ccv_basic.o ccv_resample.o ccv_classic.o ccv_daisy.o ccv_sift.o ccv_bbf.o ccv_mser.o ccv_swt.o ccv_dpm.o ccv_tld.o ccv_ferns.o
+libccv.a: ccv_cache.o ccv_memory.o 3rdparty/sha1/sha1.o 3rdparty/kissfft/kiss_fft.o 3rdparty/kissfft/kiss_fftnd.o 3rdparty/kissfft/kiss_fftr.o 3rdparty/kissfft/kiss_fftndr.o 3rdparty/kissfft/kissf_fft.o 3rdparty/kissfft/kissf_fftnd.o 3rdparty/kissfft/kissf_fftr.o 3rdparty/kissfft/kissf_fftndr.o 3rdparty/dsfmt/dSFMT.o 3rdparty/sfmt/SFMT.o ccv_io.o ccv_numeric.o ccv_algebra.o ccv_util.o ccv_basic.o ccv_resample.o ccv_classic.o ccv_daisy.o ccv_sift.o ccv_bbf.o ccv_mser.o ccv_swt.o ccv_dpm.o ccv_tld.o ccv_ferns.o
ar rcs $@ $^
ccv_io.o: ccv_io.c ccv.h ccv_internal.h io/*.c
Please sign in to comment.
Something went wrong with that request. Please try again.