From 1cdcdfba29474b196c0b7644e12b8aacb990532c Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Thu, 2 Mar 2023 11:11:53 +0100 Subject: [PATCH 1/3] fixed more interface, added test of integer interface --- .../mathcore/inc/Math/MersenneTwisterEngine.h | 9 +- math/mathcore/inc/Math/MixMaxEngine.h | 10 +- math/mathcore/inc/Math/MixMaxEngine.icc | 12 - math/mathcore/inc/Math/RanluxppEngine.h | 13 + math/mathcore/inc/Math/StdEngine.h | 3 +- math/mathcore/src/MersenneTwisterEngine.cxx | 11 +- math/mathcore/test/testMathRandom.cxx | 227 +++++++++++------- 7 files changed, 178 insertions(+), 107 deletions(-) diff --git a/math/mathcore/inc/Math/MersenneTwisterEngine.h b/math/mathcore/inc/Math/MersenneTwisterEngine.h index 5e586a42a50e3..8000aa6f9107c 100644 --- a/math/mathcore/inc/Math/MersenneTwisterEngine.h +++ b/math/mathcore/inc/Math/MersenneTwisterEngine.h @@ -47,7 +47,7 @@ namespace ROOT { @ingroup Random */ - class MersenneTwisterEngine : public TRandomEngine { + class MersenneTwisterEngine final : public TRandomEngine { public: @@ -56,6 +56,9 @@ namespace ROOT { typedef uint32_t Result_t; typedef uint32_t StateInt_t; + // number of random bits generated + static constexpr uint32_t kNumberOfBits = 32; + MersenneTwisterEngine(uint32_t seed=4357) { SetSeed(seed); @@ -75,9 +78,9 @@ namespace ROOT { } /// minimum integer taht can be generated - static unsigned int MinInt() { return 0; } + static constexpr unsigned int MinInt() { return 0; } /// maximum integer that can be generated - static unsigned int MaxInt() { return 0xffffffff; } // 2^32 -1 + static constexpr unsigned int MaxInt() { return 0xffffffff; } // 2^32 -1 static int Size() { return kSize; } diff --git a/math/mathcore/inc/Math/MixMaxEngine.h b/math/mathcore/inc/Math/MixMaxEngine.h index 4d9100d4fd72b..96a09d3c66cf3 100644 --- a/math/mathcore/inc/Math/MixMaxEngine.h +++ b/math/mathcore/inc/Math/MixMaxEngine.h @@ -19,6 +19,7 @@ #include "Math/TRandomEngine.h" + // struct rng_state_st; /// forward declare generator state // typedef struct rng_state_st rng_state_t; @@ -99,7 +100,7 @@ Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33–38 */ template - class MixMaxEngine : public TRandomEngine { + class MixMaxEngine final : public TRandomEngine { public: @@ -113,6 +114,9 @@ Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33–38 #endif typedef uint64_t Result_t; + // number of random bits generated + static constexpr uint32_t kNumberOfBits = 61; + MixMaxEngine(uint64_t seed=1); @@ -123,10 +127,10 @@ Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33–38 static int Size(); /// maximum integer that can be generated. For MIXMAX is 2^61-1 - static uint64_t MaxInt(); + static constexpr uint64_t MaxInt() { return 2305843009213693951ULL;} /// minimum integer that can be generated. For MIXMAX is 0 - static uint64_t MinInt(); + static constexpr uint64_t MinInt() { return 0;} /// set the generator seed void SetSeed(Result_t seed); diff --git a/math/mathcore/inc/Math/MixMaxEngine.icc b/math/mathcore/inc/Math/MixMaxEngine.icc index 45cfdc4810c3e..c5669cb76d989 100644 --- a/math/mathcore/inc/Math/MixMaxEngine.icc +++ b/math/mathcore/inc/Math/MixMaxEngine.icc @@ -107,18 +107,6 @@ namespace Math { return fRng->IntRndm(); } - template - uint64_t MixMaxEngine::MaxInt() { - //return mixmax::mixmax_engine::max(); - return 2305843009213693951ULL; - } - - template - uint64_t MixMaxEngine::MinInt() { - //return mixmax::mixmax_engine::min(); - return 0; - } - template void MixMaxEngine::RndmArray(int n, double *array){ // Return an array of n random numbers uniformly distributed in ]0,1] diff --git a/math/mathcore/inc/Math/RanluxppEngine.h b/math/mathcore/inc/Math/RanluxppEngine.h index 7866a1409b40f..c7dc35d2ce668 100644 --- a/math/mathcore/inc/Math/RanluxppEngine.h +++ b/math/mathcore/inc/Math/RanluxppEngine.h @@ -34,6 +34,19 @@ class RanluxppEngine final : public TRandomEngine { RanluxppEngine(uint64_t seed = 314159265); ~RanluxppEngine() override; + typedef TRandomEngine BaseType; + + typedef uint64_t Result_t; + + // number of random bits generated + static constexpr uint32_t kNumberOfBits = 48; + + /// minimum integer taht can be generated + static constexpr uint64_t MinInt() { return 0; } + + /// maximum integer taht can be generated + static constexpr uint64_t MaxInt() { return (1ULL<<48) - 1; } // 2^48 -1 + /// Generate a double-precision random number with 48 bits of randomness double Rndm() override; /// Generate a double-precision random number (non-virtual method) diff --git a/math/mathcore/inc/Math/StdEngine.h b/math/mathcore/inc/Math/StdEngine.h index 346d81407928b..8b6366143eff0 100644 --- a/math/mathcore/inc/Math/StdEngine.h +++ b/math/mathcore/inc/Math/StdEngine.h @@ -103,7 +103,8 @@ namespace ROOT { return StdEngineType::Name(); } - static uint64_t MaxInt() { return Generator::max(); } + static constexpr uint64_t MinInt() { return Generator::min(); } + static constexpr uint64_t MaxInt() { return Generator::max(); } private: diff --git a/math/mathcore/src/MersenneTwisterEngine.cxx b/math/mathcore/src/MersenneTwisterEngine.cxx index 18476514b45ce..4eee00d4f78dd 100644 --- a/math/mathcore/src/MersenneTwisterEngine.cxx +++ b/math/mathcore/src/MersenneTwisterEngine.cxx @@ -35,7 +35,13 @@ namespace Math { /// generate a random double number double MersenneTwisterEngine::Rndm_impl() { + auto y = IntRndm_impl(); + // 2.3283064365386963e-10 == 1./(max+1) -> then returned value cannot be = 1.0 + if (y) return ( (double) y * 2.3283064365386963e-10); // * Power(2,-32) + return Rndm_impl(); + } + uint32_t MersenneTwisterEngine::IntRndm_impl() { uint32_t y; @@ -71,10 +77,7 @@ namespace Math { y ^= ((y << 15) & kTemperingMaskC ); y ^= (y >> 18); - // 2.3283064365386963e-10 == 1./(max+1) -> then returned value cannot be = 1.0 - if (y) return ( (double) y * 2.3283064365386963e-10); // * Power(2,-32) - return Rndm_impl(); - + return y; } } // namespace Math diff --git a/math/mathcore/test/testMathRandom.cxx b/math/mathcore/test/testMathRandom.cxx index 9ac1655a33a93..894fed6a91b00 100644 --- a/math/mathcore/test/testMathRandom.cxx +++ b/math/mathcore/test/testMathRandom.cxx @@ -2,6 +2,7 @@ #include "Math/TRandomEngine.h" #include "Math/MersenneTwisterEngine.h" #include "Math/MixMaxEngine.h" +#include "Math/RanluxppEngine.h" //#include "Math/MyMixMaxEngine.h" //#include "Math/GSLRndmEngines.h" #include "Math/GoFTest.h" @@ -21,197 +22,255 @@ using namespace ROOT::Math; -const long long NR = 1E6; +const long long NR = 1E6; double minPValue = 1.E-3; -bool showPlots = false; +bool showPlots = false; +inline void fillBits(int64_t *v, uint64_t x) +{ + for (int i = 0; i < 64; ++i) { + v[i] += x & 1; + x >>= 1; + } +} + +template +bool testInt() +{ + Engine engine; + bool ret = true; + static_assert((~0ULL) >> (64 - Engine::kNumberOfBits) == Engine::MaxInt()); + static_assert(0 == Engine::MinInt()); + static_assert(64 - Engine::kNumberOfBits == __builtin_clzll(Engine::MaxInt())); + + std::cout << "testing integer engine " << engine.Name() << std::endl; + int N = 2 * 1000 * 1000; + int64_t vr[64]; + for (int i = 0; i < 64; ++i) + vr[i] = 0; + for (int i = 0; i < N; ++i) { + auto r = engine.IntRndm(); + fillBits(vr, r); + } + int64_t t = 1000 * 1000; + for (uint32_t i = 0; i < 64; ++i) { + if (i >= Engine::kNumberOfBits && vr[i] != 0) { + std::cout << "not zero? " << i << ' ' << vr[i] << std::endl; + ret = false; + } + if (i < Engine::kNumberOfBits && std::abs(vr[i] - t) > 3000) { + std::cout << "off? " << i << ' ' << vr[i] << std::endl; + ret = false; + } + } + return ret; +} -bool testCompatibility(const std::vector & x, const std::vector & y, double xmin = 0, double xmax = 1) { +bool testCompatibility(const std::vector &x, const std::vector &y, double xmin = 0, double xmax = 1) +{ - GoFTest gof(x.size(), x.data(), y.size(), y.data() ); + GoFTest gof(x.size(), x.data(), y.size(), y.data()); - bool ok = true; + bool ok = true; double pvalue = gof.KolmogorovSmirnov2SamplesTest(); - if ( pvalue < minPValue ) { + if (pvalue < minPValue) { std::cout << "KS Test failed with p-value " << pvalue << std::endl; - ok = false; + ok = false; + } else { + std::cout << "KS Test : OK - pvalue = " << pvalue << std::endl; } - else { - std::cout << "KS Test : OK - pvalue = " << pvalue << std::endl; - } - - if (NR < 10000) { - pvalue = gof.AndersonDarling2SamplesTest(); - if ( pvalue < minPValue ) { + if (NR < 10000) { + pvalue = gof.AndersonDarling2SamplesTest(); + if (pvalue < minPValue) { std::cout << "AD Test failed with p-value " << pvalue << std::endl; - ok = false; + ok = false; + } else { + std::cout << "AD Test : OK - pvalue = " << pvalue << std::endl; } - else { - std::cout << "AD Test : OK - pvalue = " << pvalue << std::endl; - } - } + } // do a chi2 binned test - int nbin = TMath::Min(x.size(), y.size() )/1000; - TH1D * h1 = new TH1D("h1","h1", nbin, xmin, xmax); - TH1D * h2 = new TH1D("h2","h2", nbin, xmin, xmax); - h1->FillN(x.size(), x.data(), nullptr ); - h2->FillN(y.size(), y.data(), nullptr ); + int nbin = TMath::Min(x.size(), y.size()) / 1000; + TH1D *h1 = new TH1D("h1", "h1", nbin, xmin, xmax); + TH1D *h2 = new TH1D("h2", "h2", nbin, xmin, xmax); + h1->FillN(x.size(), x.data(), nullptr); + h2->FillN(y.size(), y.data(), nullptr); pvalue = h1->Chi2Test(h2); - if ( pvalue < minPValue ) { + if (pvalue < minPValue) { std::cout << "Chi2 Test failed with p-value " << pvalue << std::endl; - //showPlots = true; - ok = false; - } - else { - std::cout << "Chi2 Test: OK - pvalue = " << pvalue << std::endl; + // showPlots = true; + ok = false; + } else { + std::cout << "Chi2 Test: OK - pvalue = " << pvalue << std::endl; } if (showPlots) { h1->Draw(); h2->SetLineColor(kRed); h2->Draw("SAME"); - if (gPad) gPad->Update(); - } - else { + if (gPad) + gPad->Update(); + } else { delete h1; - delete h2; + delete h2; } - - return ok; + return ok; } -template -bool testUniform(R1 & r1, R2 & r2) { +template +bool testUniform(R1 &r1, R2 &r2) +{ - std::vector x(NR); std::vector y(NR); - TStopwatch w; w.Start(); - for (int i = 0; i < NR; ++i) + TStopwatch w; + w.Start(); + for (int i = 0; i < NR; ++i) x[i] = r1(); w.Stop(); std::cout << "time for uniform filled for " << typeid(r1).name(); - w.Print(); + w.Print(); - w.Start(); - for (int i = 0; i < NR; ++i) + w.Start(); + for (int i = 0; i < NR; ++i) y[i] = r2(); w.Stop(); std::cout << "time for uniform filled for " << typeid(r2).name(); - w.Print(); + w.Print(); - return testCompatibility(x,y); + return testCompatibility(x, y); } -template -bool testGauss(R1 & r1, R2 & r2) { +template +bool testGauss(R1 &r1, R2 &r2) +{ - std::vector x(NR); std::vector y(NR); - TStopwatch w; w.Start(); - for (int i = 0; i < NR; ++i) - x[i] = ROOT::Math::normal_cdf(r1.Gaus(0,1),1); + TStopwatch w; + w.Start(); + for (int i = 0; i < NR; ++i) + x[i] = ROOT::Math::normal_cdf(r1.Gaus(0, 1), 1); w.Stop(); std::cout << "time for GAUS filled for " << typeid(r1).name(); - w.Print(); + w.Print(); - w.Start(); - for (int i = 0; i < NR; ++i) - y[i] = ROOT::Math::normal_cdf(r2.Gaus(0,1),1); + w.Start(); + for (int i = 0; i < NR; ++i) + y[i] = ROOT::Math::normal_cdf(r2.Gaus(0, 1), 1); w.Stop(); std::cout << "time for GAUS filled for " << typeid(r2).name(); - w.Print(); + w.Print(); - return testCompatibility(x,y); + return testCompatibility(x, y); } -bool test1() { +bool test1() +{ - bool ret = true; + bool ret = true; std::cout << "\nTesting MT vs MIXMAX " << std::endl; Random rmt; Random rmx; ret &= testUniform(rmx, rmt); ret &= testGauss(rmx, rmt); - return ret; + return ret; } -bool test2() { +bool test2() +{ - bool ret = true; + bool ret = true; std::cout << "\nTesting MIXMAX240 vs MIXMAX256" << std::endl; Random rmx1(1111); - Random> rmx2(2222); + Random> rmx2(2222); ret &= testUniform(rmx1, rmx2); ret &= testGauss(rmx1, rmx2); - return ret; + return ret; } -bool test3() { +bool test3() +{ - bool ret = true; + bool ret = true; std::cout << "\nTesting MIXMAX240 vs MIXMAX17" << std::endl; - Random rmx1(1111); - Random> rmx2(2222); + Random> rmx2(2222); ret &= testUniform(rmx1, rmx2); ret &= testGauss(rmx1, rmx2); - return ret; + return ret; } -bool test4() { +bool test4() +{ - bool ret = true; + bool ret = true; std::cout << "\nTesting MIXMAX240 vs MIXMAX240 using different seeds" << std::endl; - Random rmx1(1111); Random rmx2(2222); ret &= testUniform(rmx1, rmx2); ret &= testGauss(rmx1, rmx2); - return ret; + return ret; } +bool test5() +{ -bool testMathRandom() { + bool ret = true; + std::cout << "\nTesting RanLun++ vs MIXMAX " << std::endl; + + Random rmt; + Random rmx; + ret &= testUniform(rmx, rmt); + ret &= testGauss(rmx, rmt); + return ret; +} + +bool testMathRandom() +{ - bool ret = true; std::cout << "testing generating " << NR << " numbers " << std::endl; - ret &= test1(); - ret &= test2(); - ret &= test3(); - ret &= test4(); + ret &= test1(); + ret &= test2(); + ret &= test3(); + ret &= test4(); + ret &= test5(); + + ret &= testInt(); + ret &= testInt(); + ret &= testInt(); + ret &= testInt(); - if (!ret) Error("testMathRandom","Test Failed"); + if (!ret) + Error("testMathRandom", "Test Failed"); else std::cout << "\nTestMathRandom: OK \n"; - return ret; + return ret; } -int main() { +int main() +{ bool ret = testMathRandom(); return (ret) ? 0 : -1; } - From 02e7f21e33ec9aa5146584d4d63e9e86bfd4f5ad Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Wed, 8 Mar 2023 12:56:41 +0100 Subject: [PATCH 2/3] fix for WIN, add test for std --- math/mathcore/inc/Math/StdEngine.h | 194 ++++++++++++-------------- math/mathcore/test/testMathRandom.cxx | 29 +++- 2 files changed, 119 insertions(+), 104 deletions(-) diff --git a/math/mathcore/inc/Math/StdEngine.h b/math/mathcore/inc/Math/StdEngine.h index 8b6366143eff0..9500ac5c028e0 100644 --- a/math/mathcore/inc/Math/StdEngine.h +++ b/math/mathcore/inc/Math/StdEngine.h @@ -17,108 +17,100 @@ namespace ROOT { - namespace Math { - - - class StdRandomEngine {}; - - template - struct StdEngineType { - static const char * Name() { return "std_random_eng";} - }; - template<> - struct StdEngineType { - static const char * Name() { return "std_minstd_rand";} - }; - template<> - struct StdEngineType { - static const char * Name() { return "std_mt19937";} - }; - template<> - struct StdEngineType { - static const char * Name() { return "std_mt19937_64";} - }; - template<> - struct StdEngineType { - static const char * Name() { return "std_ranlux24";} - }; - template<> - struct StdEngineType { - static const char * Name() { return "std_ranlux48";} - }; - template<> - struct StdEngineType { - static const char * Name() { return "std_knuth_b";} - }; - template<> - struct StdEngineType { - static const char * Name() { return "std_random_device";} - }; - - - /** - @ingroup Random - Class to wrap engines fron the C++ standard random library in - the ROOT Random interface. - These cases are then used by the generic TRandomGen class - to provide TRandom interrace generators for the C++ random generators. - - See for examples the TRandomMT64 and TRandomRanlux48 generators - which are typede's to TRandomGen instantiated with some - random engine from the C++ standard library. - - */ - - template - class StdEngine { - - - public: - - typedef StdRandomEngine BaseType; - typedef typename Generator::result_type Result_t; - - StdEngine() : fGen() { - fCONS = 1./fGen.max(); - } - - - void SetSeed(Result_t seed) { fGen.seed(seed);} - - double Rndm() { - Result_t rndm = fGen(); // generate integer number according to the type - if (rndm != 0) return fCONS*rndm; - return Rndm(); - } - - Result_t IntRndm() { - return fGen(); - } - - double operator() () { - return Rndm(); - } - - static const char * Name() { - return StdEngineType::Name(); - } - - static constexpr uint64_t MinInt() { return Generator::min(); } - static constexpr uint64_t MaxInt() { return Generator::max(); } - - - private: - Generator fGen; - double fCONS; //! cached value of maximum integer value generated - }; - - - extern template class StdEngine; - extern template class StdEngine; - - } // end namespace Math +namespace Math { + +class StdRandomEngine {}; + +template +struct StdEngineType { + static const char *Name() { return "std_random_eng"; } +}; +template <> +struct StdEngineType { + static const char *Name() { return "std_minstd_rand"; } +}; +template <> +struct StdEngineType { + static const char *Name() { return "std_mt19937"; } +}; +template <> +struct StdEngineType { + static const char *Name() { return "std_mt19937_64"; } +}; +template <> +struct StdEngineType { + static const char *Name() { return "std_ranlux24"; } +}; +template <> +struct StdEngineType { + static const char *Name() { return "std_ranlux48"; } +}; +template <> +struct StdEngineType { + static const char *Name() { return "std_knuth_b"; } +}; +template <> +struct StdEngineType { + static const char *Name() { return "std_random_device"; } +}; + +/** + @ingroup Random + Class to wrap engines fron the C++ standard random library in + the ROOT Random interface. + These cases are then used by the generic TRandomGen class + to provide TRandom interrace generators for the C++ random generators. + + See for examples the TRandomMT64 and TRandomRanlux48 generators + which are typede's to TRandomGen instantiated with some + random engine from the C++ standard library. + +*/ + +template +class StdEngine { + +public: + typedef StdRandomEngine BaseType; + typedef typename Generator::result_type Result_t; + + StdEngine() : fGen() { fCONS = 1. / fGen.max(); } + + void SetSeed(Result_t seed) { fGen.seed(seed); } + + double Rndm() + { + Result_t rndm = fGen(); // generate integer number according to the type + if (rndm != 0) + return fCONS * rndm; + return Rndm(); + } + + Result_t IntRndm() { return fGen(); } + + double operator()() { return Rndm(); } + + static const char *Name() { return StdEngineType::Name(); } + + static constexpr uint64_t MinInt() { return Generator::min(); } + static constexpr uint64_t MaxInt() { return Generator::max(); } + +#ifdef _WIN64 + static constexpr uint64_t kNumberOfBits = _popcount64(Generator::max()); +#else + static constexpr uint64_t kNumberOfBits = __builtin_popcountll(Generator::max()); +#endif + +private: + Generator fGen; + double fCONS; //! cached value of maximum integer value generated +}; + +extern template class StdEngine; +extern template class StdEngine; + +} // end namespace Math } // end namespace ROOT - #endif /* ROOT_Math_StdEngine */ diff --git a/math/mathcore/test/testMathRandom.cxx b/math/mathcore/test/testMathRandom.cxx index 894fed6a91b00..a5d675dbf4018 100644 --- a/math/mathcore/test/testMathRandom.cxx +++ b/math/mathcore/test/testMathRandom.cxx @@ -3,6 +3,7 @@ #include "Math/MersenneTwisterEngine.h" #include "Math/MixMaxEngine.h" #include "Math/RanluxppEngine.h" +#include "Math/StdEngine.h" //#include "Math/MyMixMaxEngine.h" //#include "Math/GSLRndmEngines.h" #include "Math/GoFTest.h" @@ -40,9 +41,16 @@ bool testInt() { Engine engine; bool ret = true; - static_assert((~0ULL) >> (64 - Engine::kNumberOfBits) == Engine::MaxInt()); - static_assert(0 == Engine::MinInt()); - static_assert(64 - Engine::kNumberOfBits == __builtin_clzll(Engine::MaxInt())); + static_assert((~0ULL) >> (64 - Engine::kNumberOfBits) == Engine::MaxInt(), "MaxInt inconsistent with kNumberOfBits"); + static_assert(0 == Engine::MinInt(), "MinInt not 0"); +#ifdef _WIN64 + static_assert(Engine::kNumberOfBits == _popcount64(Engine::MaxInt()), "MaxInt inconsistent with kNumberOfBits"); +#else + static_assert(Engine::kNumberOfBits == __builtin_popcountll(Engine::MaxInt()), + "MaxInt inconsistent with kNumberOfBits"); + static_assert(64 - Engine::kNumberOfBits == __builtin_clzll(Engine::MaxInt()), + "MaxInt inconsistent with kNumberOfBits"); +#endif std::cout << "testing integer engine " << engine.Name() << std::endl; int N = 2 * 1000 * 1000; @@ -245,6 +253,19 @@ bool test5() return ret; } +bool test6() +{ + + bool ret = true; + std::cout << "\nTesting std::mt19937_64 vs MIXMAX " << std::endl; + + Random> rmt; + Random rmx; + ret &= testUniform(rmx, rmt); + // ret &= testGauss(rmx, rmt); + return ret; +} + bool testMathRandom() { @@ -256,11 +277,13 @@ bool testMathRandom() ret &= test3(); ret &= test4(); ret &= test5(); + ret &= test6(); ret &= testInt(); ret &= testInt(); ret &= testInt(); ret &= testInt(); + ret &= testInt>(); if (!ret) Error("testMathRandom", "Test Failed"); From 7c311bb09222e8e44416b3a4baf40a5b56c68a9d Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Wed, 8 Mar 2023 13:05:30 +0100 Subject: [PATCH 3/3] fix return type of gsl wrappers --- math/mathmore/inc/Math/GSLRndmEngines.h | 868 ++++++++++++------------ math/mathmore/src/GSLRngROOTWrapper.h | 137 ++-- 2 files changed, 492 insertions(+), 513 deletions(-) diff --git a/math/mathmore/inc/Math/GSLRndmEngines.h b/math/mathmore/inc/Math/GSLRndmEngines.h index 3cd822c3fecb0..f4d8fae074583 100644 --- a/math/mathmore/inc/Math/GSLRndmEngines.h +++ b/math/mathmore/inc/Math/GSLRndmEngines.h @@ -1,26 +1,26 @@ // @(#)root/mathmore:$Id$ // Author: L. Moneta, A. Zsenei 08/2005 - /********************************************************************** - * * - * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * - * * - * This library is free software; you can redistribute it and/or * - * modify it under the terms of the GNU General Public License * - * as published by the Free Software Foundation; either version 2 * - * of the License, or (at your option) any later version. * - * * - * This library is distributed in the hope that it will be useful, * - * but WITHOUT ANY WARRANTY; without even the implied warranty of * - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * - * General Public License for more details. * - * * - * You should have received a copy of the GNU General Public License * - * along with this library (see file COPYING); if not, write * - * to the Free Software Foundation, Inc., 59 Temple Place, Suite * - * 330, Boston, MA 02111-1307 USA, or contact the author. * - * * - **********************************************************************/ +/********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ // Header file for class GSLRandom // @@ -34,500 +34,486 @@ #include #include - namespace ROOT { namespace Math { +class GSLRngWrapper; +class GSLMCIntegrator; - class GSLRngWrapper; - class GSLMCIntegrator; - - //_________________________________________________________________ - /** - GSLRandomEngine - Base class for all GSL random engines, - normally user instantiate the derived classes - which creates internally the generator. - - The main GSL generators (see - - here) are available as derived classes - In addition to generate uniform numbers it provides method for - generating numbers according to pre-defined distributions - using the GSL functions from - - GSL random number distributions. - - - - @ingroup Random - */ - class GSLRandomEngine { - - friend class GSLMCIntegrator; - - public: - - /** - default constructor. No creation of rng is done. - If then Initialize() is called an engine is created - based on default GSL type (MT) - */ - GSLRandomEngine(); - - /** - create from an existing rng. - User manage the rng pointer which is then deleted olny by calling Terminate() - */ - GSLRandomEngine( GSLRngWrapper * rng); - - /** - Copy constructor : clone the contained GSL generator - */ - GSLRandomEngine(const GSLRandomEngine & eng); - - /** - Assignment operator : make a deep copy of the contained GSL generator - */ - GSLRandomEngine & operator=(const GSLRandomEngine & eng); - - /** - initialize the generator - If no rng is present the default one based on Mersenne and Twister is created - */ - void Initialize(); - - /** - delete pointer to contained rng - */ - void Terminate(); - - /** - call Terminate() - */ - virtual ~GSLRandomEngine(); - - /** - Generate a random number between ]0,1] - 0 is excluded and 1 is included - */ - double operator() () const; - - /** - Generate a random number between ]0,1] - 0 is excluded and 1 is included - */ - double Rndm() const { return (*this)(); } - - /** - Generate an integer number between [0,max-1] (including 0 and max-1) - if max is larger than available range of algorithm - an error message is printed and zero is returned - */ - unsigned long RndmInt(unsigned long max) const; - /** - Generate an integer number between [0,max_generator-1] (including 0 and max-1) - if max is larger than available range of algorithm - an error message is printed and zero is returned - */ - unsigned long IntRndm() const { - return RndmInt(MaxInt()); // max return the largest value the generator can give +1 - } +//_________________________________________________________________ +/** + GSLRandomEngine + Base class for all GSL random engines, + normally user instantiate the derived classes + which creates internally the generator. - /** - Generate an array of random numbers. - The iterators points to the random numbers - */ - template - void RandomArray(Iterator begin, Iterator end) const { - for ( Iterator itr = begin; itr != end; ++itr ) { - *itr = this->operator()(); - } - } + The main GSL generators (see + + here) are available as derived classes + In addition to generate uniform numbers it provides method for + generating numbers according to pre-defined distributions + using the GSL functions from + + GSL random number distributions. - /** - Generate an array of random numbers - The iterators points to the random numbers - */ - void RandomArray(double * begin, double * end) const; - - /** - return name of generator - */ - std::string Name() const; - - /** - return the state size of generator - */ - unsigned int Size() const; - - /** - return the minimum integer a generator can handle - typically this value is 0 - */ - unsigned long MinInt() const; - - /** - return the maximum integer +1 a generator can handle - - */ - unsigned long MaxInt() const; - - /** - set the random generator seed - */ - void SetSeed(unsigned int seed) const; - - - /** @name Random Distributions - Implemented using the - - GSL Random number Distributions - **/ - //@{ - /** - Gaussian distribution - default method is Box-Muller (polar method) - */ - double Gaussian(double sigma) const; - - /** - Gaussian distribution - Ziggurat method - */ - double GaussianZig(double sigma) const; - - /** - Gaussian distribution - Ratio method - */ - double GaussianRatio(double sigma) const; - /** - Gaussian Tail distribution - */ - double GaussianTail(double a, double sigma) const; - - /** - Bivariate Gaussian distribution with correlation - */ - void Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const; - - /** - Multivariate Gaussian distribution - */ - void GaussianND(const int dim, double *pars, double *covmat, double *genpars) const; - - /** - Exponential distribution - */ - double Exponential(double mu) const; - - /** - Cauchy distribution - */ - double Cauchy(double a) const; - - /** - Landau distribution - */ - double Landau() const; - - /** - Gamma distribution - */ - double Gamma(double a, double b) const; - - /** - Beta distribution - */ - double Beta(double a, double b) const; - - /** - Log Normal distribution - */ - double LogNormal(double zeta, double sigma) const; - - /** - Chi square distribution - */ - double ChiSquare(double nu) const; - - /** - F distrbution - */ - double FDist(double nu1, double nu2) const; - - /** - t student distribution - */ - double tDist(double nu) const; - - /** - Rayleigh distribution - */ - double Rayleigh(double sigma) const; - - /** - Logistic distribution - */ - double Logistic(double a) const; - - /** - Pareto distribution - */ - double Pareto(double a, double b) const; - - /** - generate random numbers in a 2D circle of radious 1 - */ - void Dir2D(double &x, double &y) const; - - /** - generate random numbers in a 3D sphere of radious 1 - */ - void Dir3D(double &x, double &y, double &z) const; - - /** - Poisson distribution - */ - unsigned int Poisson(double mu) const; - - /** - Binomial distribution - */ - unsigned int Binomial(double p, unsigned int n) const; - - /** - Negative Binomial distribution - */ - unsigned int NegativeBinomial(double p, double n) const; - - /** - Multinomial distribution - */ - std::vector Multinomial( unsigned int ntot, const std::vector & p ) const; - - //@} - - - - protected: - - /// internal method used by the derived class to set the type of generators - void SetType(GSLRngWrapper * r) { - fRng = r; - } - /// internal method to return the engine - /// Used by class like GSLMCIntegrator to set the engine - GSLRngWrapper * Engine() { - return fRng; - } - - private: - GSLRngWrapper * fRng; // pointer to GSL generator wrapper (managed by the class) - mutable unsigned int fCurTime; // current time used to seed the generator + @ingroup Random +*/ +class GSLRandomEngine { + friend class GSLMCIntegrator; - }; +public: + /** + default constructor. No creation of rng is done. + If then Initialize() is called an engine is created + based on default GSL type (MT) + */ + GSLRandomEngine(); - //_____________________________________________________________________________________ /** - Mersenne-Twister generator - gsl_rng_mt19937 from - here + create from an existing rng. + User manage the rng pointer which is then deleted olny by calling Terminate() + */ + GSLRandomEngine(GSLRngWrapper *rng); + /** + Copy constructor : clone the contained GSL generator + */ + GSLRandomEngine(const GSLRandomEngine &eng); + + /** + Assignment operator : make a deep copy of the contained GSL generator + */ + GSLRandomEngine &operator=(const GSLRandomEngine &eng); + + /** + initialize the generator + If no rng is present the default one based on Mersenne and Twister is created + */ + void Initialize(); + + /** + delete pointer to contained rng + */ + void Terminate(); + + /** + call Terminate() + */ + virtual ~GSLRandomEngine(); - @ingroup Random + /** + Generate a random number between ]0,1] + 0 is excluded and 1 is included */ - class GSLRngMT : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngMT(); - }; + double operator()() const; - //_____________________________________________________________________________________ /** - Old Ranlux generator (James, Luscher) (default luxury level, p = 223) - (This is eequivalent to TRandom1 with default luxury level) - see here + Generate a random number between ]0,1] + 0 is excluded and 1 is included + */ + double Rndm() const { return (*this)(); } - @ingroup Random + /** + Generate an integer number between [0,max-1] (including 0 and max-1) + if max is larger than available range of algorithm + an error message is printed and zero is returned */ - class GSLRngRanLux : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngRanLux(); - }; + unsigned long RndmInt(unsigned long max) const; + /** + Generate an integer number between [0,max_generator-1] (including 0 and max-1) + if max is larger than available range of algorithm + an error message is printed and zero is returned + */ + unsigned long long IntRndm() const + { + return RndmInt(MaxInt()); // max return the largest value the generator can give +1 + } - //_____________________________________________________________________________________ /** - Second generation of Ranlux generator for single precision with luxury level of 1 - (It throws away 202 values for every 12 used) - see here + Generate an array of random numbers. + The iterators points to the random numbers + */ + template + void RandomArray(Iterator begin, Iterator end) const + { + for (Iterator itr = begin; itr != end; ++itr) { + *itr = this->operator()(); + } + } - @ingroup Random + /** + Generate an array of random numbers + The iterators points to the random numbers */ - class GSLRngRanLuxS1 : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngRanLuxS1(); - }; - typedef GSLRngRanLuxS1 GSLRngRanLux1; // for backward compatibility + void RandomArray(double *begin, double *end) const; - //_____________________________________________________________________________________ /** - Second generation of Ranlux generator for Single precision with luxury level of 2 - (It throws away 397 value for every 12 used) - see here + return name of generator + */ + std::string Name() const; - @ingroup Random + /** + return the state size of generator */ - class GSLRngRanLuxS2 : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngRanLuxS2(); - }; - typedef GSLRngRanLuxS2 GSLRngRanLux2; // for backward compatibility + unsigned int Size() const; + + /** + return the minimum integer a generator can handle + typically this value is 0 + */ + unsigned long MinInt() const; - //_____________________________________________________________________________________ /** - Double precision (48 bits) version of Second generation of Ranlux generator with luxury level of 1 - (It throws away 202 value for every 12 used) - see here + return the maximum integer +1 a generator can handle - @ingroup Random + */ + unsigned long MaxInt() const; + + /** + set the random generator seed + */ + void SetSeed(unsigned int seed) const; + + /** @name Random Distributions + Implemented using the + + GSL Random number Distributions + **/ + //@{ + /** + Gaussian distribution - default method is Box-Muller (polar method) */ - class GSLRngRanLuxD1 : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngRanLuxD1(); - }; + double Gaussian(double sigma) const; - //_____________________________________________________________________________________ /** - Double precision (48 bits) version of Second generation of Ranlux generator with luxury level of 2 - (It throws away 397 value for every 12 used) - see here + Gaussian distribution - Ziggurat method + */ + double GaussianZig(double sigma) const; - @ingroup Random + /** + Gaussian distribution - Ratio method */ - class GSLRngRanLuxD2 : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngRanLuxD2(); - }; - typedef GSLRngRanLuxD2 GSLRngRanLux48; // for backward compatibility + double GaussianRatio(double sigma) const; + /** + Gaussian Tail distribution + */ + double GaussianTail(double a, double sigma) const; + /** + Bivariate Gaussian distribution with correlation + */ + void Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const; - //_____________________________________________________________________________________ /** - Tausworthe generator by L'Ecuyer - see here + Multivariate Gaussian distribution + */ + void GaussianND(const int dim, double *pars, double *covmat, double *genpars) const; - @ingroup Random + /** + Exponential distribution */ - class GSLRngTaus : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngTaus(); - }; + double Exponential(double mu) const; - //_____________________________________________________________________________________ /** - Lagged Fibonacci generator by Ziff - see here + Cauchy distribution + */ + double Cauchy(double a) const; - @ingroup Random + /** + Landau distribution */ - class GSLRngGFSR4 : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngGFSR4(); - }; + double Landau() const; - //_____________________________________________________________________________________ /** - Combined multiple recursive generator (L'Ecuyer) - see here + Gamma distribution + */ + double Gamma(double a, double b) const; - @ingroup Random + /** + Beta distribution */ - class GSLRngCMRG : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngCMRG(); - }; + double Beta(double a, double b) const; - //_____________________________________________________________________________________ /** - 5-th order multiple recursive generator (L'Ecuyer, Blouin and Coutre) - see here + Log Normal distribution + */ + double LogNormal(double zeta, double sigma) const; - @ingroup Random + /** + Chi square distribution */ - class GSLRngMRG : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngMRG(); - }; + double ChiSquare(double nu) const; - //_____________________________________________________________________________________ /** - BSD rand() generator - gsl_rmg_rand from - here + F distrbution + */ + double FDist(double nu1, double nu2) const; - @ingroup Random + /** + t student distribution */ - class GSLRngRand : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngRand(); - }; + double tDist(double nu) const; - //_____________________________________________________________________________________ /** - RANMAR generator - see here + Rayleigh distribution + */ + double Rayleigh(double sigma) const; - @ingroup Random + /** + Logistic distribution + */ + double Logistic(double a) const; + + /** + Pareto distribution + */ + double Pareto(double a, double b) const; + + /** + generate random numbers in a 2D circle of radious 1 + */ + void Dir2D(double &x, double &y) const; + + /** + generate random numbers in a 3D sphere of radious 1 */ - class GSLRngRanMar : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngRanMar(); - }; + void Dir3D(double &x, double &y, double &z) const; - //_____________________________________________________________________________________ /** - MINSTD generator (Park and Miller) - see here + Poisson distribution + */ + unsigned int Poisson(double mu) const; - @ingroup Random + /** + Binomial distribution */ - class GSLRngMinStd : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngMinStd(); - }; + unsigned int Binomial(double p, unsigned int n) const; - /** MixMax generator based on ROOT::Math::MixMaxEngine of N=240 + /** + Negative Binomial distribution + */ + unsigned int NegativeBinomial(double p, double n) const; - @ingroup Random + /** + Multinomial distribution */ - class GSLRngMixMax : public GSLRandomEngine { - public: - typedef GSLRandomEngine BaseType; - GSLRngMixMax(); - ~GSLRngMixMax() override; // we need a dtcor since is not a standard GSL engine - }; + std::vector Multinomial(unsigned int ntot, const std::vector &p) const; + + //@} + +protected: + /// internal method used by the derived class to set the type of generators + void SetType(GSLRngWrapper *r) { fRng = r; } + + /// internal method to return the engine + /// Used by class like GSLMCIntegrator to set the engine + GSLRngWrapper *Engine() { return fRng; } + +private: + GSLRngWrapper *fRng; // pointer to GSL generator wrapper (managed by the class) + mutable unsigned int fCurTime; // current time used to seed the generator +}; + +//_____________________________________________________________________________________ +/** + Mersenne-Twister generator + gsl_rng_mt19937 from + here + + + @ingroup Random +*/ +class GSLRngMT : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngMT(); +}; + +//_____________________________________________________________________________________ +/** + Old Ranlux generator (James, Luscher) (default luxury level, p = 223) + (This is eequivalent to TRandom1 with default luxury level) + see here + + @ingroup Random +*/ +class GSLRngRanLux : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngRanLux(); +}; + +//_____________________________________________________________________________________ +/** + Second generation of Ranlux generator for single precision with luxury level of 1 + (It throws away 202 values for every 12 used) + see here + + @ingroup Random +*/ +class GSLRngRanLuxS1 : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngRanLuxS1(); +}; +typedef GSLRngRanLuxS1 GSLRngRanLux1; // for backward compatibility + +//_____________________________________________________________________________________ +/** + Second generation of Ranlux generator for Single precision with luxury level of 2 + (It throws away 397 value for every 12 used) + see here + + @ingroup Random +*/ +class GSLRngRanLuxS2 : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngRanLuxS2(); +}; +typedef GSLRngRanLuxS2 GSLRngRanLux2; // for backward compatibility + +//_____________________________________________________________________________________ +/** + Double precision (48 bits) version of Second generation of Ranlux generator with luxury level of 1 + (It throws away 202 value for every 12 used) + see here + + @ingroup Random +*/ +class GSLRngRanLuxD1 : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngRanLuxD1(); +}; + +//_____________________________________________________________________________________ +/** + Double precision (48 bits) version of Second generation of Ranlux generator with luxury level of 2 + (It throws away 397 value for every 12 used) + see here + + @ingroup Random +*/ +class GSLRngRanLuxD2 : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngRanLuxD2(); +}; +typedef GSLRngRanLuxD2 GSLRngRanLux48; // for backward compatibility + +//_____________________________________________________________________________________ +/** + Tausworthe generator by L'Ecuyer + see here + + @ingroup Random +*/ +class GSLRngTaus : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngTaus(); +}; + +//_____________________________________________________________________________________ +/** + Lagged Fibonacci generator by Ziff + see here + + @ingroup Random +*/ +class GSLRngGFSR4 : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngGFSR4(); +}; + +//_____________________________________________________________________________________ +/** + Combined multiple recursive generator (L'Ecuyer) + see here + + @ingroup Random +*/ +class GSLRngCMRG : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngCMRG(); +}; + +//_____________________________________________________________________________________ +/** + 5-th order multiple recursive generator (L'Ecuyer, Blouin and Coutre) + see here + + @ingroup Random +*/ +class GSLRngMRG : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngMRG(); +}; + +//_____________________________________________________________________________________ +/** + BSD rand() generator + gsl_rmg_rand from + here + + @ingroup Random +*/ +class GSLRngRand : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngRand(); +}; + +//_____________________________________________________________________________________ +/** + RANMAR generator + see here + + @ingroup Random +*/ +class GSLRngRanMar : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngRanMar(); +}; + +//_____________________________________________________________________________________ +/** + MINSTD generator (Park and Miller) + see here + + @ingroup Random +*/ +class GSLRngMinStd : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngMinStd(); +}; + +/** MixMax generator based on ROOT::Math::MixMaxEngine of N=240 + + @ingroup Random +*/ +class GSLRngMixMax : public GSLRandomEngine { +public: + typedef GSLRandomEngine BaseType; + GSLRngMixMax(); + ~GSLRngMixMax() override; // we need a dtcor since is not a standard GSL engine +}; } // namespace Math } // namespace ROOT // random functions specialization for GSL -// needs to be defined after defining GSLRandomEngine class +// needs to be defined after defining GSLRandomEngine class #include "Math/GSLRandomFunctions.h" #endif /* ROOT_Math_GSLRndmEngines */ - diff --git a/math/mathmore/src/GSLRngROOTWrapper.h b/math/mathmore/src/GSLRngROOTWrapper.h index bb6ac48f5fe5f..2a1fbe09dae73 100644 --- a/math/mathmore/src/GSLRngROOTWrapper.h +++ b/math/mathmore/src/GSLRngROOTWrapper.h @@ -21,89 +21,82 @@ namespace ROOT { - namespace Math { - - /** - * class for wrapping ROOT Engines in gsl_rng types which can be used as extra - * GSL random number generators - * For this we need to implment functions which will be called by gsl_rng. - * The functions (Seed, Rndm, IntRndm) are passed in the gsl_rng_type and used to build a gsl_rng object. - * When gsl_rng is alloacated, only the memory state is allocated using calloc(1,size), which gives a memory - * block of the given bytes and it initializes to zero. Therefore no constructor of GSLRngROOTWrapper can be called - * and also we cannot call non-static member funciton of the class. - * The underlined ROOT engine is then built and deleted using the functions CreateEngine() and FreeEngine(), - * called by the specific GSLRandomEngine class that instantiates for the the generator (e.g. GSLRngMixMax) - * - **/ - template - struct GSLRngROOTWrapper { - - Engine *fEngine = nullptr; - - // non need to have specific ctor/dtor since we cannot call them - - // function called by the specific GSLRndmEngine to create the ROOT engine - static void CreateEngine(gsl_rng *r) - { - // the engine pointer is retrieved from r - GSLRngROOTWrapper *wrng = ((GSLRngROOTWrapper *)r->state); - //printf("calling create engine - engine pointer : %p\n", wrng->fEngine); +namespace Math { + +/** + * class for wrapping ROOT Engines in gsl_rng types which can be used as extra + * GSL random number generators + * For this we need to implment functions which will be called by gsl_rng. + * The functions (Seed, Rndm, IntRndm) are passed in the gsl_rng_type and used to build a gsl_rng object. + * When gsl_rng is alloacated, only the memory state is allocated using calloc(1,size), which gives a memory + * block of the given bytes and it initializes to zero. Therefore no constructor of GSLRngROOTWrapper can be called + * and also we cannot call non-static member funciton of the class. + * The underlined ROOT engine is then built and deleted using the functions CreateEngine() and FreeEngine(), + * called by the specific GSLRandomEngine class that instantiates for the the generator (e.g. GSLRngMixMax) + * + **/ +template +struct GSLRngROOTWrapper { + + Engine *fEngine = nullptr; + + // non need to have specific ctor/dtor since we cannot call them + + // function called by the specific GSLRndmEngine to create the ROOT engine + static void CreateEngine(gsl_rng *r) + { + // the engine pointer is retrieved from r + GSLRngROOTWrapper *wrng = ((GSLRngROOTWrapper *)r->state); + // printf("calling create engine - engine pointer : %p\n", wrng->fEngine); + if (!wrng->fEngine) + wrng->fEngine = new Engine(); + // do nothing in case it is already created (e.g. when calling with default_seed) + } + + static double Rndm(void *p) { return ((GSLRngROOTWrapper *)p)->fEngine->operator()(); } + static unsigned long long IntRndm(void *p) { return ((GSLRngROOTWrapper *)p)->fEngine->IntRndm(); } + + static void Seed(void *p, unsigned long seed) + { + auto wrng = ((GSLRngROOTWrapper *)p); + // (GSL calls at the beginning with the defaul seed (typically zero)) + // printf("calling the seed function with %d on %p and engine %p\n", seed, p, wrng->fEngine); + if (seed == gsl_rng_default_seed) { + seed = 111; // avoid using 0 that for ROOT means a specific seed if (!wrng->fEngine) wrng->fEngine = new Engine(); - // do nothing in case it is already created (e.g. when calling with default_seed) } + assert(wrng->fEngine != nullptr); + wrng->fEngine->SetSeed(seed); + } + static void FreeEngine(gsl_rng *r) + { + auto wrng = ((GSLRngROOTWrapper *)r->state); + if (wrng->fEngine) + delete wrng->fEngine; + wrng->fEngine = nullptr; + } + + static unsigned long Max() { return Engine::MaxInt(); } + static unsigned long Min() { return Engine::MinInt(); } + static size_t Size() { return sizeof(GSLRngROOTWrapper); } + static std::string Name() { return std::string("GSL_") + Engine::Name(); } +}; - static double Rndm(void *p) { return ((GSLRngROOTWrapper *)p)->fEngine->operator()(); } - static unsigned long IntRndm(void *p) { return ((GSLRngROOTWrapper *)p)->fEngine->IntRndm(); } - - static void Seed(void *p, unsigned long seed) - { - auto wrng = ((GSLRngROOTWrapper *)p); - // (GSL calls at the beginning with the defaul seed (typically zero)) - //printf("calling the seed function with %d on %p and engine %p\n", seed, p, wrng->fEngine); - if (seed == gsl_rng_default_seed) { - seed = 111; // avoid using 0 that for ROOT means a specific seed - if (!wrng->fEngine) wrng->fEngine = new Engine(); - } - assert(wrng->fEngine != nullptr); - wrng->fEngine->SetSeed(seed); - } - static void FreeEngine(gsl_rng *r) - { - auto wrng = ((GSLRngROOTWrapper *)r->state); - if (wrng->fEngine) - delete wrng->fEngine; - wrng->fEngine = nullptr; - } - - static unsigned long Max() { return Engine::MaxInt(); } - static unsigned long Min() { return Engine::MinInt(); } - static size_t Size() { return sizeof(GSLRngROOTWrapper); } - static std::string Name() { return std::string("GSL_") + Engine::Name(); } - }; - - } // end namespace Math +} // end namespace Math } // end namespace ROOT #include "Math/MixMaxEngine.h" -// now define and implement the specific types +// now define and implement the specific types -typedef ROOT::Math::GSLRngROOTWrapper> GSLMixMaxWrapper; +typedef ROOT::Math::GSLRngROOTWrapper> GSLMixMaxWrapper; -static const std::string gsl_mixmax_name = GSLMixMaxWrapper::Name(); -static const gsl_rng_type mixmax_type = -{ - gsl_mixmax_name.c_str(), - GSLMixMaxWrapper::Max(), - GSLMixMaxWrapper::Min(), - GSLMixMaxWrapper::Size(), - &GSLMixMaxWrapper::Seed, - &GSLMixMaxWrapper::IntRndm, - &GSLMixMaxWrapper::Rndm -}; +static const std::string gsl_mixmax_name = GSLMixMaxWrapper::Name(); +static const gsl_rng_type mixmax_type = {gsl_mixmax_name.c_str(), GSLMixMaxWrapper::Max(), GSLMixMaxWrapper::Min(), + GSLMixMaxWrapper::Size(), &GSLMixMaxWrapper::Seed, &GSLMixMaxWrapper::IntRndm, + &GSLMixMaxWrapper::Rndm}; const gsl_rng_type *gsl_rng_mixmax = &mixmax_type; #endif -