ashikk
/
random_well512a_simd
mirror of https://github.com/ppettersson/random_well512a_simd.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Initial version.

This commit is contained in:
Peter Pettersson 2013-04-20 19:56:34 +02:00
commit fbac6f95d4
10 changed files with 600 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
RandomWELL512a.cpp Normal file
View File

@ -0,0 +1,82 @@
// Copyright (c) 2013, Peter Pettersson
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
// 2. 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.
//
// 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.
//
// The views and conclusions contained in the software and documentation are those
// of the authors and should not be interpreted as representing official policies,
// either expressed or implied, of the FreeBSD Project.
#include "RandomWELL512a.h"
#include <stdlib.h> // Needed for rand().
#include <memory.h> // Needed for memcpy().
RandomWELL512a::RandomWELL512a(int seed)
: index(0)
{
srand(seed);
for (int i = 0; i < 16; ++i)
state[i] = rand();
}
RandomWELL512a::RandomWELL512a(unsigned *seed)
: index(0)
{
memcpy(state, seed, 16 * sizeof(unsigned));
}
unsigned RandomWELL512a::GetUnsigned()
{
#define MUTATE_LEFT(value, shift) value ^ (value << shift)
#define MUTATE_RIGHT(value, shift) value ^ (value >> shift)
#define MUTATE_LEFT_MIX(value, shift, mix) value ^ ((value << shift) & mix)
unsigned index_9 = (index + 9) & 15;
unsigned index_13 = (index + 13) & 15;
unsigned index_15 = (index + 15) & 15;
unsigned state_index = state[index];
unsigned state_index_9 = state[index_9];
unsigned state_index_13 = state[index_13];
unsigned state_index_15 = state[index_15];
unsigned z1 = MUTATE_LEFT(state_index, 16);
z1 ^= MUTATE_LEFT(state_index_13, 15);
unsigned z2 = MUTATE_RIGHT(state_index_9, 11);
unsigned result0 = z1 ^ z2;
state[index] = result0;
unsigned result1 = MUTATE_LEFT(state_index_15, 2);
result1 ^= MUTATE_LEFT(z1, 18);;
result1 ^= z2 << 28;
result1 ^= MUTATE_LEFT_MIX(result0, 5, 0xda442d24U);
index = index_15;
state[index] = result1;
return result1;
#undef MUTATE_LEFT
#undef MUTATE_RIGHT
#undef MUTATE_LEFT_MIX
}

50
RandomWELL512a.h Normal file
View File

@ -0,0 +1,50 @@
// Copyright (c) 2013, Peter Pettersson
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
// 2. 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.
//
// 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.
//
// The views and conclusions contained in the software and documentation are those
// of the authors and should not be interpreted as representing official policies,
// either expressed or implied, of the FreeBSD Project.
#ifndef RANDOM_WELL512A_H
#define RANDOM_WELL512A_H
class RandomWELL512a
{
public:
RandomWELL512a(int seed);
RandomWELL512a(unsigned *seed);
unsigned GetUnsigned();
double GetDouble()
{
const double kToFloat = 2.32830643653869628906e-10;
return GetUnsigned() * kToFloat;
}
private:
unsigned state[16];
unsigned index;
};
#endif // RANDOM_WELL512A_H

88
RandomWELL512a_SSE2.cpp Normal file
View File

@ -0,0 +1,88 @@
// Copyright (c) 2013, Peter Pettersson
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
// 2. 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.
//
// 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.
//
// The views and conclusions contained in the software and documentation are those
// of the authors and should not be interpreted as representing official policies,
// either expressed or implied, of the FreeBSD Project.
#include "RandomWELL512a_SSE2.h"
#include <stdlib.h> // Needed for rand().
#include <memory.h> // Needed for memcpy().
RandomWELL512a_SSE2::RandomWELL512a_SSE2(int seed)
: index(0)
{
srand(seed);
for (int i = 0; i < 16; ++i)
xmm_state[i] = _mm_set_epi32(rand(), rand(), rand(), rand());
}
RandomWELL512a_SSE2::RandomWELL512a_SSE2(unsigned *seed)
: index(0)
{
for (int i = 0; i < 16; ++i)
xmm_state[i] = _mm_set_epi32(seed[i + 48], seed[i + 32], seed[i + 16], seed[i]);
}
void RandomWELL512a_SSE2::GetUnsigned4(unsigned *result4)
{
// Calculations.
#define MUTATE_LEFT(value, shift) _mm_xor_si128(value, _mm_slli_epi32(value, shift))
#define MUTATE_RIGHT(value, shift) _mm_xor_si128(value, _mm_srli_epi32(value, shift))
#define MUTATE_LEFT_MIX(value, shift, mix) _mm_xor_si128(value, _mm_and_si128(_mm_slli_epi32(value, shift), mix))
unsigned index_15 = (index + 15) & 15;
__m128i state_index = xmm_state[index];
__m128i state_index_9 = xmm_state[(index + 9) & 15];
__m128i state_index_13 = xmm_state[(index + 13) & 15];
__m128i state_index_15 = xmm_state[index_15];
const __m128i kMix = _mm_set1_epi32(0xda442d24);
__m128i z1 = _mm_xor_si128(MUTATE_LEFT(state_index, 16), MUTATE_LEFT(state_index_13, 15));
__m128i z2 = MUTATE_RIGHT(state_index_9, 11);
__m128i result0 = _mm_xor_si128(z1, z2);
xmm_state[index] = result0;
__m128i result1 = MUTATE_LEFT(state_index_15, 2);
result1 = _mm_xor_si128(result1, MUTATE_LEFT(z1, 18));
result1 = _mm_xor_si128(result1, _mm_slli_epi32(z2, 28));
result1 = _mm_xor_si128(result1, MUTATE_LEFT_MIX(result0, 5, kMix));
index = index_15;
xmm_state[index] = result1;
_mm_storeu_si128((__m128i *)result4, result1);
#undef MUTATE_LEFT
#undef MUTATE_RIGHT
#undef MUTATE_LEFT_MIX
}
void RandomWELL512a_SSE2::GetDouble4(double *result4)
{
unsigned unsignedResult[4];
GetUnsigned4(unsignedResult);
const double kToFloat = 2.32830643653869628906e-10;
for (unsigned loop = 0; loop < 4; ++loop)
result4[loop] = unsignedResult[loop] * kToFloat;
}

68
RandomWELL512a_SSE2.h Normal file
View File

@ -0,0 +1,68 @@
// Copyright (c) 2013, Peter Pettersson
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
// 2. 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.
//
// 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.
//
// The views and conclusions contained in the software and documentation are those
// of the authors and should not be interpreted as representing official policies,
// either expressed or implied, of the FreeBSD Project.
#ifndef RANDOM_WELL512A_SSE2_H
#define RANDOM_WELL512A_SSE2_H
#include <emmintrin.h> // SSE2 instrinsics.
class RandomWELL512a_SSE2
{
public:
RandomWELL512a_SSE2(int seed);
RandomWELL512a_SSE2(unsigned *seed);
unsigned GetUnsigned()
{
if (resultIndex < 4)
return result[resultIndex++];
GetUnsigned4(result);
resultIndex = 1;
return result[0];
}
double GetDouble()
{
const double kToFloat = 2.32830643653869628906e-10;
return GetUnsigned() * kToFloat;
}
void GetUnsigned4(unsigned *result4);
void GetDouble4(double *result4);
private:
__m128i xmm_state[16];
unsigned index;
// Helper to allow us to return one number per call.
unsigned result[4];
unsigned resultIndex;
};
#endif // RANDOM_WELL512A_SSE2_H

55
Timer.h Normal file
View File

@ -0,0 +1,55 @@
// Copyright (c) 2013, Peter Pettersson
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
// 2. 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.
//
// 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.
//
// The views and conclusions contained in the software and documentation are those
// of the authors and should not be interpreted as representing official policies,
// either expressed or implied, of the FreeBSD Project.
#ifndef TIMER_H
#define TIMER_H
#include <time.h>
#include <iostream>
class Timer
{
public:
Timer()
{
startTime = clock();
}
void Report(const char *msg = NULL)
{
clock_t totalTime = clock() - startTime;
float seconds = (float)totalTime / CLOCKS_PER_SEC;
if (msg)
std::cout << msg << " ";
std::cout << "took " << seconds << " seconds." << std::endl;
}
private:
clock_t startTime;
};
#endif // TIMER_H

63
WELL512a.cpp Normal file
View File

@ -0,0 +1,63 @@
/*************************************************************************
Copyright: Francois Panneton and Pierre L'Ecuyer, Université de Montréal
Makoto Matsumoto, Hiroshima University
Notice: This code can be used freely for personal, academic,
or non-commercial purposes. For commercial purposes,
please contact P. L'Ecuyer at: lecuyer@iro.UMontreal.ca
This code can also be used under the terms of the GNU General Public
License as published by the Free Software Foundation, either version 3
of the License, or any later version. See the GPL licence at URL
http://www.gnu.org/licenses
**************************************************************************/
#include "WELL512a.h"
#define W 32
#define R 16
#define P 0
#define M1 13
#define M2 9
#define M3 5
#define MAT0POS(t,v) (v^(v>>t))
#define MAT0NEG(t,v) (v^(v<<(-(t))))
#define MAT3NEG(t,v) (v<<(-(t)))
#define MAT4NEG(t,b,v) (v ^ ((v<<(-(t))) & b))
#define V0 STATE[state_i ]
#define VM1 STATE[(state_i+M1) & 0x0000000fU]
#define VM2 STATE[(state_i+M2) & 0x0000000fU]
#define VM3 STATE[(state_i+M3) & 0x0000000fU]
#define VRm1 STATE[(state_i+15) & 0x0000000fU]
#define VRm2 STATE[(state_i+14) & 0x0000000fU]
#define newV0 STATE[(state_i+15) & 0x0000000fU]
#define newV1 STATE[state_i ]
#define newVRm1 STATE[(state_i+14) & 0x0000000fU]
#define FACT 2.32830643653869628906e-10
static unsigned int state_i = 0;
static unsigned int STATE[R];
static unsigned int z0, z1, z2;
void InitWELLRNG512a (unsigned int *init)
{
int j;
state_i = 0;
for (j = 0; j < R; j++)
STATE[j] = init[j];
}
double WELLRNG512a (void)
{
z0 = VRm1;
z1 = MAT0NEG(-16, V0) ^ MAT0NEG(-15, VM1);
z2 = MAT0POS(11, VM2);
newV1 = z1 ^ z2;
newV0 = MAT0NEG(-2, z0) ^ MAT0NEG(-18, z1) ^ MAT3NEG(-28, z2) ^ MAT4NEG(-5, 0xda442d24U, newV1);
state_i = (state_i + 15) & 0x0000000fU;
return ((double) STATE[state_i]) * FACT;
}

19
WELL512a.h Normal file
View File

@ -0,0 +1,19 @@
/*************************************************************************
Copyright: Francois Panneton and Pierre L'Ecuyer, Université de Montréal
Makoto Matsumoto, Hiroshima University
Notice: This code can be used freely for personal, academic,
or non-commercial purposes. For commercial purposes,
please contact P. L'Ecuyer at: lecuyer@iro.UMontreal.ca
This code can also be used under the terms of the GNU General Public
License as published by the Free Software Foundation, either version 3
of the License, or any later version. See the GPL licence at URL
http://www.gnu.org/licenses
**************************************************************************/
/* seed is an array of 16 unsigned int's used as seed of the generator */
void InitWELLRNG512a (unsigned int *seed);
double WELLRNG512a (void);

26
license.txt Normal file
View File

@ -0,0 +1,26 @@
Copyright (c) 2013, Peter Pettersson
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. 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.
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.
The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies,
either expressed or implied, of the FreeBSD Project.

135
main.cpp Normal file
View File

@ -0,0 +1,135 @@
// Copyright (c) 2013, Peter Pettersson
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
// 2. 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.
//
// 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.
//
// The views and conclusions contained in the software and documentation are those
// of the authors and should not be interpreted as representing official policies,
// either expressed or implied, of the FreeBSD Project.
#include "Timer.h"
#include "RandomWELL512a.h"
#include "RandomWELL512a_SSE2.h"
#include "WELL512a.h"
void Benchmark()
{
const unsigned kSeed = 123;
const unsigned kNumLoops = 5;
const unsigned kNumIterations = 40000000;
double *result0 = new double [kNumIterations];
double *result1 = new double [kNumIterations];
double *result2 = new double [kNumIterations];
double *result3 = new double [kNumIterations];
srand(kSeed);
unsigned seed[16];
for (unsigned i = 0; i < 16; ++i)
seed[i] = rand();
for (unsigned loop = 0; loop < kNumLoops; ++loop)
{
// RAND implementation.
srand(kSeed);
Timer timer0;
for (unsigned i = 0; i < kNumIterations; ++i)
result0[0] = rand() / (double)RAND_MAX;
timer0.Report("Rand(): ");
// WELL512 C++ implementation.
RandomWELL512a random(kSeed);
Timer timer1;
for (unsigned i = 0; i < kNumIterations; ++i)
result1[i] = random.GetDouble();
timer1.Report("WELL512 C++: ");
// WELL512 SSE2 implementation.
RandomWELL512a_SSE2 randomSSE2(kSeed);
Timer timer2;
for (unsigned i = 0; i < kNumIterations; i += 4)
randomSSE2.GetDouble4(result2 + i);
timer2.Report("WELL512 SSE2: ");
// WELL512 C implementation.
InitWELLRNG512a(seed);
Timer timer3;
for (unsigned i = 0; i < kNumIterations; ++i)
result3[i] = WELLRNG512a();
timer3.Report("WELL512 C: ");
std::cout << "---" << std::endl;
}
delete [] result0;
delete [] result1;
delete [] result2;
delete [] result3;
}
// Verify that the SIMD implementation returns the same values as the original
// algorithm would.
void Test()
{
const int kSeed = 123;
const unsigned kNumIterations = 4 * 1024;
double *result0 = new double [kNumIterations];
double *result1 = new double [kNumIterations];
srand(kSeed);
unsigned seed[4 * 16];
for (unsigned i = 0; i < 4 * 16; ++i)
seed[i] = rand();
RandomWELL512a randomWell0(seed + 0 * 16);
RandomWELL512a randomWell1(seed + 1 * 16);
RandomWELL512a randomWell2(seed + 2 * 16);
RandomWELL512a randomWell3(seed + 3 * 16);
RandomWELL512a_SSE2 randomWellSSE2(seed);
for (unsigned i = 0; i < kNumIterations; i += 4)
{
result0[i + 0] = randomWell0.GetDouble();
result0[i + 1] = randomWell1.GetDouble();
result0[i + 2] = randomWell2.GetDouble();
result0[i + 3] = randomWell3.GetDouble();
result1[i + 0] = randomWellSSE2.GetDouble();
result1[i + 1] = randomWellSSE2.GetDouble();
result1[i + 2] = randomWellSSE2.GetDouble();
result1[i + 3] = randomWellSSE2.GetDouble();
}
if (memcmp(result0, result1, kNumIterations * sizeof(double)))
std::cout << "ERROR: C++ vs SSE2: The results don't match!" << std::endl;
else
std::cout << "C++ vs SSE2: Results match" << std::endl;
delete [] result0;
delete [] result1;
}
int main(int argc, char **argv)
{
Benchmark();
Test();
return 0;
}

14
readme.txt Normal file
View File

@ -0,0 +1,14 @@
A SIMD optimized version of the WELL512a random algorithm.
For a detailed description of how the algorithm works, please see the following pages:
1) http://en.wikipedia.org/wiki/Well_equidistributed_long-period_linear
2) http://www.iro.umontreal.ca/~panneton/WELLRNG.html
Benchmarking shows that the SSE2 implementation is generally about 2-3 times faster than
the C implementation by the original authors.
It requires a SSE2 capable CPU (i.e. Intel Pentium 4 or newer, AMD Athlon 64 or newer).
All the source code except WELL512a.cpp|h that contains the original implementation
are released under the BSD license which allows you to include this in any open
or closed source project.