From 78322cc701786d27138a41b18e5122506cff5fb7 Mon Sep 17 00:00:00 2001 From: Marcel Stimberg Date: Mon, 9 Sep 2024 18:19:37 +0200 Subject: [PATCH] Replace randomkit by standard library implementation Closes #1483 --- LICENSE | 71 ---- brian2/devices/cpp_standalone/codeobject.py | 6 +- brian2/devices/cpp_standalone/device.py | 38 +- .../devices/cpp_standalone/templates/main.cpp | 2 +- .../cpp_standalone/templates/objects.cpp | 15 +- .../devices/cpp_standalone/templates/run.cpp | 4 +- brian2/random/randomkit/randomkit.c | 402 ------------------ brian2/random/randomkit/randomkit.h | 189 -------- 8 files changed, 20 insertions(+), 707 deletions(-) delete mode 100644 brian2/random/randomkit/randomkit.c delete mode 100644 brian2/random/randomkit/randomkit.h diff --git a/LICENSE b/LICENSE index c06e706f5..ac81f18a1 100644 --- a/LICENSE +++ b/LICENSE @@ -661,77 +661,6 @@ THE SOFTWARE. ------------------------------------------------------------------------------- -The Mersenne-Twister implementation "Randomkit" in brian2.random.randomkit is -based on code by Jean-Sebastien Roy, provided under the MIT license: - -Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org) - -Permission is hereby granted, free of charge, to any person obtaining a -copy of this software and associated documentation files (the -"Software"), to deal in the Software without restriction, including -without limitation the rights to use, copy, modify, merge, publish, -distribute, sublicense, and/or sell copies of the Software, and to -permit persons to whom the Software is furnished to do so, subject to -the following conditions: - -The above copyright notice and this permission notice shall be included -in all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS -OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - --- - -The rk_random and rk_seed functions algorithms and the original design of -the Mersenne Twister RNG: - -Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, -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. - -3. The names of its contributors may not 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. - --- - -Original algorithm for the implementation of rk_interval function from -Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by -Magnus Jonsson. - --- - -Constants used in the rk_double implementation by Isaku Wada. - --- - The floor division and modulo implementations for integer numbers in C++ standalone are based on the implementation in Cython (https://github.com/cython/cython), published under the Apache License 2.0: diff --git a/brian2/devices/cpp_standalone/codeobject.py b/brian2/devices/cpp_standalone/codeobject.py index 9778c538d..f0d581335 100644 --- a/brian2/devices/cpp_standalone/codeobject.py +++ b/brian2/devices/cpp_standalone/codeobject.py @@ -139,14 +139,14 @@ def generate_rand_code(rand_func, owner): else: thread_number = "omp_get_thread_num()" if rand_func == "rand": - rk_call = "rk_double" + rk_call = "_uniform_random" elif rand_func == "randn": - rk_call = "rk_gauss" + rk_call = "_normal_random" else: raise AssertionError(rand_func) code = """ double _%RAND_FUNC%(const int _vectorisation_idx) { - return %RK_CALL%(brian::_mersenne_twister_states[%THREAD_NUMBER%]); + return brian::%RK_CALL%(brian::_mersenne_twister_generators[%THREAD_NUMBER%]); } """ code = replace( diff --git a/brian2/devices/cpp_standalone/device.py b/brian2/devices/cpp_standalone/device.py index a4050a512..72eef38f4 100644 --- a/brian2/devices/cpp_standalone/device.py +++ b/brian2/devices/cpp_standalone/device.py @@ -19,7 +19,6 @@ import numpy as np -import brian2 from brian2.codegen.codeobject import check_compiler_kwds from brian2.codegen.cpp_prefs import get_compiler_and_args, get_msvc_env from brian2.codegen.generators.cpp_generator import c_data_type @@ -228,8 +227,8 @@ def __init__(self): self.extra_compile_args = [] self.define_macros = [] self.headers = [] - self.include_dirs = ["brianlib/randomkit"] - self.library_dirs = ["brianlib/randomkit"] + self.include_dirs = [] + self.library_dirs = [] self.runtime_library_dirs = [] self.run_environment_variables = {} if sys.platform.startswith("darwin"): @@ -874,12 +873,11 @@ def generate_main_source(self, writer): main_lines.append(f"for (int _i=0; _i<{nb_threads}; _i++)") if seed is None: # random main_lines.append( - " rk_randomseed(brian::_mersenne_twister_states[_i]);" + " brian::_mersenne_twister_generators[_i] = std::mt19937(_rd());" ) else: main_lines.append( - f" rk_seed({seed!r}L + _i," - " brian::_mersenne_twister_states[_i]);" + f"brian::_mersenne_twister_generators[_i].seed({seed!r}L + _i);" ) else: raise NotImplementedError(f"Unknown main queue function type {func}") @@ -1084,28 +1082,6 @@ def copy_source_files(self, writer, directory): os.path.join(directory, "brianlib", "stdint_compat.h"), ) - # Copy the RandomKit implementation - if not os.path.exists(os.path.join(directory, "brianlib", "randomkit")): - os.mkdir(os.path.join(directory, "brianlib", "randomkit")) - shutil.copy2( - os.path.join( - os.path.split(inspect.getsourcefile(brian2))[0], - "random", - "randomkit", - "randomkit.c", - ), - os.path.join(directory, "brianlib", "randomkit", "randomkit.c"), - ) - shutil.copy2( - os.path.join( - os.path.split(inspect.getsourcefile(brian2))[0], - "random", - "randomkit", - "randomkit.h", - ), - os.path.join(directory, "brianlib", "randomkit", "randomkit.h"), - ) - def _insert_func_namespace(self, func, code_object, namespace): impl = func.implementations[self.code_object_class()] func_namespace = impl.get_namespace(code_object.owner) @@ -1583,9 +1559,7 @@ def build( for codeobj in self.code_objects.values() for source_file in codeobj.compiler_kwds.get("sources", []) ] - additional_source_files += codeobj_source_files + [ - "brianlib/randomkit/randomkit.c" - ] + additional_source_files += codeobj_source_files for d in ["code_objects", "results", "static_arrays"]: ensure_directory(os.path.join(directory, d)) @@ -1703,7 +1677,6 @@ def delete(self, code=True, data=True, directory=True, force=False): fnames.extend( [ os.path.join("brianlib", "spikequeue.h"), - os.path.join("brianlib", "randomkit", "randomkit.h"), os.path.join("brianlib", "stdint_compat.h"), ] ) @@ -1731,7 +1704,6 @@ def delete(self, code=True, data=True, directory=True, force=False): if directory: directories = [ - os.path.join("brianlib", "randomkit"), "brianlib", "code_objects", "results", diff --git a/brian2/devices/cpp_standalone/templates/main.cpp b/brian2/devices/cpp_standalone/templates/main.cpp index 5758fff70..ba57df813 100644 --- a/brian2/devices/cpp_standalone/templates/main.cpp +++ b/brian2/devices/cpp_standalone/templates/main.cpp @@ -5,7 +5,6 @@ {{ openmp_pragma('include') }} #include "run.h" #include "brianlib/common_math.h" -#include "randomkit.h" {% for codeobj in code_objects | sort(attribute='name') %} #include "code_objects/{{codeobj.name}}.h" @@ -36,6 +35,7 @@ void set_from_command_line(const std::vector args) } int main(int argc, char **argv) { + std::random_device _rd; std::vector args(argv + 1, argv + argc); if (args.size() >=2 && args[0] == "--results_dir") { diff --git a/brian2/devices/cpp_standalone/templates/objects.cpp b/brian2/devices/cpp_standalone/templates/objects.cpp index 166caab29..0008db62f 100644 --- a/brian2/devices/cpp_standalone/templates/objects.cpp +++ b/brian2/devices/cpp_standalone/templates/objects.cpp @@ -20,7 +20,7 @@ set_variable_from_value(name, {{array_name}}, var_size, (char)atoi(s_value.c_str #include "brianlib/dynamic_array.h" #include "brianlib/stdint_compat.h" #include "network.h" -#include "randomkit.h" +#include #include #include #include @@ -32,7 +32,9 @@ set_variable_from_value(name, {{array_name}}, var_size, (char)atoi(s_value.c_str namespace brian { std::string results_dir = "results/"; // can be overwritten by --results_dir command line arg -std::vector< rk_state* > _mersenne_twister_states; +std::vector< std::mt19937 > _mersenne_twister_generators; +std::uniform_real_distribution _uniform_random; +std::normal_distribution _normal_random; //////////////// networks ///////////////// {% for net in networks | sort(attribute='name') %} @@ -223,8 +225,9 @@ void _init_arrays() {% endfor %} // Random number generator states + std::random_device rd; for (int i=0; i<{{openmp_pragma('get_num_threads')}}; i++) - _mersenne_twister_states.push_back(new rk_state()); + _mersenne_twister_generators.push_back(std::mt19937(rd())); } void _load_arrays() @@ -369,7 +372,7 @@ void _dealloc_arrays() #include "brianlib/dynamic_array.h" #include "brianlib/stdint_compat.h" #include "network.h" -#include "randomkit.h" +#include #include {{ openmp_pragma('include') }} @@ -377,7 +380,9 @@ namespace brian { extern std::string results_dir; // In OpenMP we need one state per thread -extern std::vector< rk_state* > _mersenne_twister_states; +extern std::vector< std::mt19937 > _mersenne_twister_generators; +extern std::uniform_real_distribution _uniform_random; +extern std::normal_distribution _normal_random; //////////////// clocks /////////////////// {% for clock in clocks | sort(attribute='name') %} diff --git a/brian2/devices/cpp_standalone/templates/run.cpp b/brian2/devices/cpp_standalone/templates/run.cpp index 20a429bbf..798493dd0 100644 --- a/brian2/devices/cpp_standalone/templates/run.cpp +++ b/brian2/devices/cpp_standalone/templates/run.cpp @@ -2,7 +2,7 @@ #include #include "objects.h" #include -#include "randomkit.h" +#include {% for codeobj in code_objects | sort(attribute='name') %} #include "code_objects/{{codeobj.name}}.h" @@ -22,8 +22,6 @@ void brian_start() brian::{{clock.name}}.dt = brian::{{array_specs[clock.variables['dt']]}}; brian::{{clock.name}}.t = brian::{{array_specs[clock.variables['t']]}}; {% endfor %} - for (int i=0; i<{{openmp_pragma('get_num_threads')}}; i++) - rk_randomseed(brian::_mersenne_twister_states[i]); // Note that this seed can be potentially replaced in main.cpp } void brian_end() diff --git a/brian2/random/randomkit/randomkit.c b/brian2/random/randomkit/randomkit.c deleted file mode 100644 index 169c288dc..000000000 --- a/brian2/random/randomkit/randomkit.c +++ /dev/null @@ -1,402 +0,0 @@ -/* Random kit 1.3 */ - -/* - * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org) - * - * The rk_random and rk_seed functions algorithms and the original design of - * the Mersenne Twister RNG: - * - * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - * 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. - * - * 3. The names of its contributors may not 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. - * - * Original algorithm for the implementation of rk_interval function from - * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by - * Magnus Jonsson. - * - * Constants used in the rk_double implementation by Isaku Wada. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* static char const rcsid[] = - "@(#) $Jeannot: randomkit.c,v 1.28 2005/07/21 22:14:09 js Exp $"; */ -#include -#include -#include -#include -#include -#include - -#ifdef _WIN32 -/* - * Windows - * XXX: we have to use this ugly defined(__GNUC__) because it is not easy to - * detect the compiler used in distutils itself - */ -#if (defined(__GNUC__) && defined(NPY_NEEDS_MINGW_TIME_WORKAROUND)) - -/* - * FIXME: ideally, we should set this to the real version of MSVCRT. We need - * something higher than 0x601 to enable _ftime64 and co - */ -#define __MSVCRT_VERSION__ 0x0700 -#include -#include - -/* - * mingw msvcr lib import wrongly export _ftime, which does not exist in the - * actual msvc runtime for version >= 8; we make it an alias to _ftime64, which - * is available in those versions of the runtime - */ -#define _FTIME(x) _ftime64((x)) -#else -#include -#include -#define _FTIME(x) _ftime((x)) -#endif - -#ifndef RK_NO_WINCRYPT -/* Windows crypto */ -#ifndef _WIN32_WINNT -#define _WIN32_WINNT 0x0400 -#endif -#include -#include -#endif - -#else -/* Unix */ -#include -#include -#include -#endif - -#include "randomkit.h" - -#ifndef RK_DEV_URANDOM -#define RK_DEV_URANDOM "/dev/urandom" -#endif - -#ifndef RK_DEV_RANDOM -#define RK_DEV_RANDOM "/dev/random" -#endif - -char *rk_strerror[RK_ERR_MAX] = -{ - "no error", - "random device unvavailable" -}; - -/* static functions */ -static unsigned long rk_hash(unsigned long key); - -void -rk_seed(unsigned long seed, rk_state *state) -{ - int pos; - seed &= 0xffffffffUL; - - /* Knuth's PRNG as used in the Mersenne Twister reference implementation */ - for (pos = 0; pos < RK_STATE_LEN; pos++) { - state->key[pos] = seed; - seed = (1812433253UL * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffUL; - } - state->pos = RK_STATE_LEN; - state->gauss = 0; - state->has_gauss = 0; - state->has_binomial = 0; -} - -/* Thomas Wang 32 bits integer hash function */ -unsigned long -rk_hash(unsigned long key) -{ - key += ~(key << 15); - key ^= (key >> 10); - key += (key << 3); - key ^= (key >> 6); - key += ~(key << 11); - key ^= (key >> 16); - return key; -} - -rk_error -rk_randomseed(rk_state *state) -{ -#ifndef _WIN32 - struct timeval tv; -#else - struct _timeb tv; -#endif - int i; - - if (rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) { - /* ensures non-zero key */ - state->key[0] |= 0x80000000UL; - state->pos = RK_STATE_LEN; - state->gauss = 0; - state->has_gauss = 0; - state->has_binomial = 0; - - for (i = 0; i < 624; i++) { - state->key[i] &= 0xffffffffUL; - } - return RK_NOERR; - } - -#ifndef _WIN32 - gettimeofday(&tv, NULL); - rk_seed(rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec) - ^ rk_hash(clock()), state); -#else - _FTIME(&tv); - rk_seed(rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()), state); -#endif - - return RK_ENODEV; -} - -/* Magic Mersenne Twister constants */ -#define N 624 -#define M 397 -#define MATRIX_A 0x9908b0dfUL -#define UPPER_MASK 0x80000000UL -#define LOWER_MASK 0x7fffffffUL - -/* Slightly optimised reference implementation of the Mersenne Twister */ -unsigned long -rk_random(rk_state *state) -{ - unsigned long y; - - if (state->pos == RK_STATE_LEN) { - int i; - - for (i = 0; i < N - M; i++) { - y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); - state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - } - for (; i < N - 1; i++) { - y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); - state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - } - y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK); - state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A); - - state->pos = 0; - } - y = state->key[state->pos++]; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - -long -rk_long(rk_state *state) -{ - return rk_ulong(state) >> 1; -} - -unsigned long -rk_ulong(rk_state *state) -{ -#if ULONG_MAX <= 0xffffffffUL - return rk_random(state); -#else - return (rk_random(state) << 32) | (rk_random(state)); -#endif -} - -unsigned long -rk_interval(unsigned long max, rk_state *state) -{ - unsigned long mask = max, value; - - if (max == 0) { - return 0; - } - /* Smallest bit mask >= max */ - mask |= mask >> 1; - mask |= mask >> 2; - mask |= mask >> 4; - mask |= mask >> 8; - mask |= mask >> 16; -#if ULONG_MAX > 0xffffffffUL - mask |= mask >> 32; -#endif - - /* Search a random value in [0..mask] <= max */ -#if ULONG_MAX > 0xffffffffUL - if (max <= 0xffffffffUL) { - while ((value = (rk_random(state) & mask)) > max); - } - else { - while ((value = (rk_ulong(state) & mask)) > max); - } -#else - while ((value = (rk_ulong(state) & mask)) > max); -#endif - return value; -} - -double -rk_double(rk_state *state) -{ - /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ - long a = rk_random(state) >> 5, b = rk_random(state) >> 6; - return (a * 67108864.0 + b) / 9007199254740992.0; -} - -void -rk_fill(void *buffer, size_t size, rk_state *state) -{ - unsigned long r; - unsigned char *buf = (unsigned char *)buffer; - - for (; size >= 4; size -= 4) { - r = rk_random(state); - *(buf++) = r & 0xFF; - *(buf++) = (r >> 8) & 0xFF; - *(buf++) = (r >> 16) & 0xFF; - *(buf++) = (r >> 24) & 0xFF; - } - - if (!size) { - return; - } - r = rk_random(state); - for (; size; r >>= 8, size --) { - *(buf++) = (unsigned char)(r & 0xFF); - } -} - -rk_error -rk_devfill(void *buffer, size_t size, int strong) -{ -#ifndef _WIN32 - FILE *rfile; - int done; - - if (strong) { - rfile = fopen(RK_DEV_RANDOM, "rb"); - } - else { - rfile = fopen(RK_DEV_URANDOM, "rb"); - } - if (rfile == NULL) { - return RK_ENODEV; - } - done = fread(buffer, size, 1, rfile); - fclose(rfile); - if (done) { - return RK_NOERR; - } -#else - -#ifndef RK_NO_WINCRYPT - HCRYPTPROV hCryptProv; - BOOL done; - - if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL, - CRYPT_VERIFYCONTEXT) || !hCryptProv) { - return RK_ENODEV; - } - done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer); - CryptReleaseContext(hCryptProv, 0); - if (done) { - return RK_NOERR; - } -#endif - -#endif - return RK_ENODEV; -} - -rk_error -rk_altfill(void *buffer, size_t size, int strong, rk_state *state) -{ - rk_error err; - - err = rk_devfill(buffer, size, strong); - if (err) { - rk_fill(buffer, size, state); - } - return err; -} - -double -rk_gauss(rk_state *state) -{ - if (state->has_gauss) { - const double tmp = state->gauss; - state->gauss = 0; - state->has_gauss = 0; - return tmp; - } - else { - double f, x1, x2, r2; - - do { - x1 = 2.0*rk_double(state) - 1.0; - x2 = 2.0*rk_double(state) - 1.0; - r2 = x1*x1 + x2*x2; - } - while (r2 >= 1.0 || r2 == 0.0); - - /* Box-Muller transform */ - f = sqrt(-2.0*log(r2)/r2); - /* Keep for next call */ - state->gauss = f*x1; - state->has_gauss = 1; - return f*x2; - } -} diff --git a/brian2/random/randomkit/randomkit.h b/brian2/random/randomkit/randomkit.h deleted file mode 100644 index e049488ee..000000000 --- a/brian2/random/randomkit/randomkit.h +++ /dev/null @@ -1,189 +0,0 @@ -/* Random kit 1.3 */ - -/* - * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org) - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -/* @(#) $Jeannot: randomkit.h,v 1.24 2005/07/21 22:14:09 js Exp $ */ - -/* - * Typical use: - * - * { - * rk_state state; - * unsigned long seed = 1, random_value; - * - * rk_seed(seed, &state); // Initialize the RNG - * ... - * random_value = rk_random(&state); // Generate random values in [0..RK_MAX] - * } - * - * Instead of rk_seed, you can use rk_randomseed which will get a random seed - * from /dev/urandom (or the clock, if /dev/urandom is unavailable): - * - * { - * rk_state state; - * unsigned long random_value; - * - * rk_randomseed(&state); // Initialize the RNG with a random seed - * ... - * random_value = rk_random(&state); // Generate random values in [0..RK_MAX] - * } - */ - -/* - * Useful macro: - * RK_DEV_RANDOM: the device used for random seeding. - * defaults to "/dev/urandom" - */ - -#include - -#ifndef _RANDOMKIT_ -#define _RANDOMKIT_ - -#define RK_STATE_LEN 624 - -typedef struct rk_state_ -{ - unsigned long key[RK_STATE_LEN]; - int pos; - int has_gauss; /* !=0: gauss contains a gaussian deviate */ - double gauss; - - /* The rk_state structure has been extended to store the following - * information for the binomial generator. If the input values of n or p - * are different than nsave and psave, then the other parameters will be - * recomputed. RTK 2005-09-02 */ - - int has_binomial; /* !=0: following parameters initialized for - binomial */ - double psave; - long nsave; - double r; - double q; - double fm; - long m; - double p1; - double xm; - double xl; - double xr; - double c; - double laml; - double lamr; - double p2; - double p3; - double p4; - -} -rk_state; - -typedef enum { - RK_NOERR = 0, /* no error */ - RK_ENODEV = 1, /* no RK_DEV_RANDOM device */ - RK_ERR_MAX = 2 -} rk_error; - -/* error strings */ -extern char *rk_strerror[RK_ERR_MAX]; - -/* Maximum generated random value */ -#define RK_MAX 0xFFFFFFFFUL - -#ifdef __cplusplus -extern "C" { -#endif - -/* - * Initialize the RNG state using the given seed. - */ -extern void rk_seed(unsigned long seed, rk_state *state); - -/* - * Initialize the RNG state using a random seed. - * Uses /dev/random or, when unavailable, the clock (see randomkit.c). - * Returns RK_NOERR when no errors occurs. - * Returns RK_ENODEV when the use of RK_DEV_RANDOM failed (for example because - * there is no such device). In this case, the RNG was initialized using the - * clock. - */ -extern rk_error rk_randomseed(rk_state *state); - -/* - * Returns a random unsigned long between 0 and RK_MAX inclusive - */ -extern unsigned long rk_random(rk_state *state); - -/* - * Returns a random long between 0 and LONG_MAX inclusive - */ -extern long rk_long(rk_state *state); - -/* - * Returns a random unsigned long between 0 and ULONG_MAX inclusive - */ -extern unsigned long rk_ulong(rk_state *state); - -/* - * Returns a random unsigned long between 0 and max inclusive. - */ -extern unsigned long rk_interval(unsigned long max, rk_state *state); - -/* - * Returns a random double between 0.0 and 1.0, 1.0 excluded. - */ -extern double rk_double(rk_state *state); - -/* - * fill the buffer with size random bytes - */ -extern void rk_fill(void *buffer, size_t size, rk_state *state); - -/* - * fill the buffer with randombytes from the random device - * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is - * On Unix, if strong is defined, RK_DEV_RANDOM is used. If not, RK_DEV_URANDOM - * is used instead. This parameter has no effect on Windows. - * Warning: on most unixes RK_DEV_RANDOM will wait for enough entropy to answer - * which can take a very long time on quiet systems. - */ -extern rk_error rk_devfill(void *buffer, size_t size, int strong); - -/* - * fill the buffer using rk_devfill if the random device is available and using - * rk_fill if is is not - * parameters have the same meaning as rk_fill and rk_devfill - * Returns RK_ENODEV if the device is unavailable, or RK_NOERR if it is - */ -extern rk_error rk_altfill(void *buffer, size_t size, int strong, - rk_state *state); - -/* - * return a random gaussian deviate with variance unity and zero mean. - */ -extern double rk_gauss(rk_state *state); - -#ifdef __cplusplus -} -#endif - -#endif /* _RANDOMKIT_ */