commit fbac6f95d40a18f4b722cc8251f7e155667fa439 Author: Peter Pettersson Date: Sat Apr 20 19:56:34 2013 +0200 Initial version. diff --git a/RandomWELL512a.cpp b/RandomWELL512a.cpp new file mode 100644 index 0000000..11653c8 --- /dev/null +++ b/RandomWELL512a.cpp @@ -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 // Needed for rand(). +#include // 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 +} diff --git a/RandomWELL512a.h b/RandomWELL512a.h new file mode 100644 index 0000000..89535f4 --- /dev/null +++ b/RandomWELL512a.h @@ -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 diff --git a/RandomWELL512a_SSE2.cpp b/RandomWELL512a_SSE2.cpp new file mode 100644 index 0000000..a101528 --- /dev/null +++ b/RandomWELL512a_SSE2.cpp @@ -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 // Needed for rand(). +#include // 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; +} diff --git a/RandomWELL512a_SSE2.h b/RandomWELL512a_SSE2.h new file mode 100644 index 0000000..bc4fe28 --- /dev/null +++ b/RandomWELL512a_SSE2.h @@ -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 // 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 diff --git a/Timer.h b/Timer.h new file mode 100644 index 0000000..41722bc --- /dev/null +++ b/Timer.h @@ -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 +#include + +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 diff --git a/WELL512a.cpp b/WELL512a.cpp new file mode 100644 index 0000000..88b18ad --- /dev/null +++ b/WELL512a.cpp @@ -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; +} diff --git a/WELL512a.h b/WELL512a.h new file mode 100644 index 0000000..102992d --- /dev/null +++ b/WELL512a.h @@ -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); diff --git a/license.txt b/license.txt new file mode 100644 index 0000000..3bad03b --- /dev/null +++ b/license.txt @@ -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. diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..7e7f59f --- /dev/null +++ b/main.cpp @@ -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; +} diff --git a/readme.txt b/readme.txt new file mode 100644 index 0000000..03f0e72 --- /dev/null +++ b/readme.txt @@ -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.