Add 'xo-randomgen/' from commit '44a31724ec'

git-subtree-dir: xo-randomgen
git-subtree-mainline: 57170366da
git-subtree-split: 44a31724ec
This commit is contained in:
Roland Conybeare 2025-05-10 18:54:18 -05:00
commit ae78de305f
21 changed files with 788 additions and 0 deletions

6
xo-randomgen/.gitignore vendored Normal file
View file

@ -0,0 +1,6 @@
# clangd (run via lsp) keeps state here
.cache
# typical build directory
.build
# compile_commands.json: symlink to build directory, should be created manually
compile_commands.json

View file

@ -0,0 +1,43 @@
# randomgen/CMakeLists.txt
cmake_minimum_required(VERSION 3.10)
project(randomgen VERSION 0.1)
include(GNUInstallDirs)
include(cmake/xo-bootstrap-macros.cmake)
xo_cxx_toplevel_options2()
# ----------------------------------------------------------------
# cmake -DCMAKE_BUILD_TYPE=debug
xo_toplevel_debug_config2()
# ----------------------------------------------------------------
# cmake -DCMAKE_BUILD_TYPE=coverage
xo_toplevel_coverage_config2()
# ----------------------------------------------------------------
# bespoke (usually temporary) c++ settings
set(PROJECT_CXX_FLAGS "")
#set(PROJECT_CXX_FLAGS "-fconcepts-diagnostics-depth=2")
add_definitions(${PROJECT_CXX_FLAGS})
# ----------------------------------------------------------------
add_subdirectory(example)
#add_subdirectory(utest)
# ----------------------------------------------------------------
# output targets
set(SELF_LIB randomgen)
xo_add_headeronly_library(${SELF_LIB})
xo_install_library4(${SELF_LIB} ${PROJECT_NAME}Targets)
xo_export_cmake_config(${PROJECT_NAME} ${PROJECT_VERSION} ${PROJECT_NAME}Targets)
# ----------------------------------------------------------------
# install additional components
install(TARGETS randomgen_ex1 DESTINATION bin/randomgen/example)

18
xo-randomgen/README.md Normal file
View file

@ -0,0 +1,18 @@
# random number generators
## Getting Started
### build + install dependencies
- see [github/Rconybea/cmake](https://github.com/Rconybea/xo-cmake) -- cmake modules
### to build + install locally
```
$ cd randomgen
$ mkdir build
$ cd build
$ PREFIX=/usr/local # for example
$ cmake -DCMAKE_MODULE_PATH=${PREFIX}/share/cmake -DCMAKE_PREFIX_PATH=$(PREFIX) -DCMAKE_INSTALL_PREFIX=${PREFIX} ..
$ make
$ make install
```

View file

@ -0,0 +1,4 @@
@PACKAGE_INIT@
include("${CMAKE_CURRENT_LIST_DIR}/randomgenTargets.cmake")
check_required_components("@PROJECT_NAME@")

View file

@ -0,0 +1,35 @@
# ----------------------------------------------------------------
# for example:
# $ PREFIX=/usr/local
# $ cmake -DCMAKE_MODULE_PATH=prefix -DCMAKE_INSTALL_PREFIX=$PREFIX -B .build
#
# will get
# CMAKE_MODULE_PATH
# from xo-cmake-config --cmake-module-path
#
# and expect .cmake macros in
# CMAKE_MODULE_PATH/xo_macros/xo_cxx.cmake
# ----------------------------------------------------------------
find_program(XO_CMAKE_CONFIG_EXECUTABLE NAMES xo-cmake-config REQUIRED)
if ("${XO_CMAKE_CONFIG_EXECUTABLE}" STREQUAL "XO_CMAKE_CONFIG_EXECUTABLE-NOT_FOUND")
message(FATAL "could not find xo-cmake-config executable")
endif()
message(STATUS "XO_CMAKE_CONFIG_EXECUTABLE=${XO_CMAKE_CONFIG_EXECUTABLE}")
if (NOT XO_SUBMODULE_BUILD)
if (("${CMAKE_MODULE_PATH}" STREQUAL "") OR ("${CMAKE_MODULE_PATH}" STREQUAL prefix))
# default to typical install location for xo-project-macros
execute_process(COMMAND ${XO_CMAKE_CONFIG_EXECUTABLE} --cmake-module-path OUTPUT_VARIABLE CMAKE_MODULE_PATH)
message(STATUS "CMAKE_MODULE_PATH=${CMAKE_MODULE_PATH}")
endif()
endif()
# needs to have been installed somewhere on CMAKE_MODULE_PATH,
# (e.g. from xo-cmake with the same value for CMAKE_INSTALL_PREFIX)
#
include(xo_macros/xo_cxx)
xo_cxx_bootstrap_message()

View file

@ -0,0 +1,2 @@
add_subdirectory(ex1)
add_subdirectory(ex2)

View file

@ -0,0 +1,2 @@
add_executable(randomgen_ex1 ex1.cpp)
xo_include_options2(randomgen_ex1)

View file

@ -0,0 +1,24 @@
/* @file ex1.cpp */
#include "xo/randomgen/xoshiro256.hpp"
#include <algorithm>
#include <iostream>
using namespace xo;
using namespace xo::rng;
int
main(int argc, char ** argv) {
xoshiro256ss rng{123456789};
std::array<std::uint64_t, 20> v;
std::generate(v.begin(), v.end(), rng);
for (std::uint64_t i=0; i<v.size(); ++i)
std::cout << "v[" << i << "]: " << v[i] << std::endl;
return 0;
} /*main*/
/* end ex1.cpp */

View file

@ -0,0 +1,3 @@
add_executable(randomgen_ex2 ex2.cpp)
xo_include_options2(randomgen_ex2)
xo_dependency(randomgen_ex2 indentlog)

View file

@ -0,0 +1,16 @@
/* @file ex2.cpp */
#include "xo/randomgen/xoshiro256.hpp"
#include "xo/randomgen/random_seed.hpp"
using namespace xo;
using namespace xo::rng;
int
main(int argc, char ** argv) {
Seed<xoshiro256ss> seed;
xoshiro256ss eng(seed);
} /*main*/
/* end ex2.cpp */

View file

@ -0,0 +1,29 @@
/* @file bernoulligen.hpp */
#pragma once
#include "generator.hpp"
#include <random>
namespace xo {
namespace rng {
/* Engine: e.g. xo::rng::xoshiro256ss or std::mt19937 */
template <class Engine>
class bernoulligen : public generator<Engine, std::bernoulli_distribution<double>> {
public:
using generator_type = generator<Engine, std::bernoulli_distribution<double>>;
template <class Engine>
static generator_type make(Engine engine, double prob) {
return generator_type::make(std::move(engine),
std::bernoulli_distribution<double>(prob));
}
template <class Engine>
static generator_type conflip(Engine engine) {
return generator_type::make(std::move(engine),
std::bernoulli_distribution<double>(0.5));
}
};
} /*namespace rng*/
} /*namespace xo*/

View file

@ -0,0 +1,40 @@
/* @file distribution_concept.hpp */
#pragma once
#include <concepts>
namespace xo {
namespace rng {
template <class Distribution, class Engine>
concept distribution_concept = requires(Distribution dist, typename Distribution::param_type p) {
typename Distribution::result_type;
typename Distribution::param_type;
{ Distribution() };
{ Distribution(p) };
{ dist.reset() };
{ dist.param() };
{ dist.param(p) };
// { dist(g) }; // generator g satisfying engine_concept
// { dist(g, p) };
{ dist.min() };
{ dist.max() };
{ dist == dist };
{ dist != dist };
// os << dist
// is >> dist
}
#ifdef __clang__
// std::copyable apparently not available in clang11 ?
#else
&& std::copyable<Distribution>
&& std::copyable<typename Distribution::param_type>
&& std::equality_comparable<typename Distribution::param_type>
#endif
;
} /*namespace rng*/
} /*namespace xo*/
/* end distribution_concept.hpp */

View file

@ -0,0 +1,78 @@
/* @file engine_concept.hpp */
#pragma once
#include <concepts>
#include <random>
namespace std {
#ifdef __clang__
# if __clang_major__ <= 11
template < class T >
concept integral = std::is_integral_v<T>;
template < class T >
concept signed_integral = std::integral<T> && std::is_signed_v<T>;
template < class T >
concept unsigned_integral
= std::integral<T> && !std::signed_integral<T>;
template< class F, class... Args >
concept invocable
= requires(F&& f, Args&&... args) {
std::invoke(std::forward<F>(f), std::forward<Args>(args)...);
/* not required to be equality-preserving */
};
template< typename G >
concept uniform_random_bit_generator
= std::invocable<G&>
&& std::unsigned_integral<std::invoke_result_t<G&>>
&& requires { { G::min() } -> std::same_as<std::invoke_result_t<G&>>;
{ G::max() } -> std::same_as<std::invoke_result_t<G&>>;
requires std::bool_constant<(G::min() < G::max())>::value; };
# endif
#else
/* uniform_random_bit_generator provided by gcc 12.3.2 */
/* uniform_random_bit_generator provided by clang 16 */
#endif
} /*namespace std*/
namespace xo {
namespace rng {
/* an engine generates psuedo-random bits.
* given
* RngEngine eng = ...;
*
* RngEngine::result_type x = eng();
*
* puts random bits into x.
*/
template <class RngEngine>
concept engine_concept = requires(RngEngine engine, typename RngEngine::result_type r) {
/* note: the first 4 requirements characterize UniformRandomBitGenerator */
typename RngEngine::result_type;
{ RngEngine(r) };
{ engine.min() } -> std::same_as<typename RngEngine::result_type>;
{ engine.max() } -> std::same_as<typename RngEngine::result_type>;
/* must return value in closed interval [.min(), .max()] */
{ engine() } -> std::same_as<typename RngEngine::result_type>;
{ engine.seed() };
{ engine.seed(r) };
{ engine == engine };
{ engine != engine };
}
#ifdef __clang__
// std::copyable apparently not available in clang11 ?
#else
&& std::copyable<RngEngine>
#endif
&& std::uniform_random_bit_generator<RngEngine>;
} /*namespace rng*/
} /*namespace xo*/
/* end engine_concept.hpp */

View file

@ -0,0 +1,21 @@
/* @file exponentialgen.hpp */
#pragma once
#include "generator.hpp"
#include <random>
namespace xo {
namespace rng {
template <class Engine>
class exponentialgen : public generator<Engine, std::exponential_distribution<double>> {
public:
using generator_type = generator<Engine, std::exponential_distribution<double>>;
template <class Engine>
static generator_type make(Engine eng, double lambda) {
return make_generator(std::move(eng), std::exponential_distribution<double>(lambda));
}
};
} /*namespace rng*/
} /*namespace xo*/

View file

@ -0,0 +1,108 @@
/* @file gaussianpairgen.hpp */
#pragma once
#include "generator.hpp"
#include <random>
#include <array>
namespace xo {
namespace random {
/* editor bait: 2d normal, normal xy
*
* if
* N1 ~ N(0,1)
* N2 ~ N(0,1)
* are two indepenent, normally-distributed r.v's with
* mean 0 and variance 1, then
* let
* A = | 1 0 | X = | N1 |
* | r q | | N2 |
*
* with r^2 + q^2 = 1
*
* and consider
* A.X = | N1 | := | Y1 |
* | r.N1 + q.N2 | | Y2 |
*
* Y1, Y2 both have mean 0,
* since both are linear combination of 0-mean N(0,1) variables
*
* Var(Y1) = 1
* Var(Y2) = r^2.Var(N1) + q^2.Var(N2)
* = r^2 + q^2
* = 1
*
* (since N1,N2 indept, and Var(N1)=Var(N2)=1)
*
* Cov(Y1,Y2) = r.Cov(N1,N1) + q.Cov(N1,N2)
* = r.Var(N1)
* = r
*
* (since Cov(N1,N2)=0)
*
* we have correlation coefficient for Y1,Y2:
*
* Cov(Y1,Y2)
* p(Y1,Y2) = --------------------
* sqrt(Var(Y1).Var(Y2))
*
* = r
*/
template<typename FloatType>
class gaussianpair_dist {
public:
using result_type = std::array<FloatType, 2>;
public:
/* generate pairs of gaussian N(0,1) random numbers,
* with correlation coefficient rho
*
* Require:
* - rho in the interval [-1, +1]
*/
explicit gaussianpair_dist(FloatType rho)
: r_(rho), q_(std::sqrt(1.0 - rho*rho)) {}
template<typename Engine>
result_type operator()(Engine & engine) {
FloatType n1 = this->ndist_(engine);
FloatType n2 = this->ndist_(engine);
FloatType y1 = n1;
FloatType y2 = this->r_ * n1 + this->q_ * n2;
return {y1, y2};
} /*operator()*/
private:
/* correlation coefficient r
* 2nd random variable Y2 in each pair will be constructed by
* r.N1 + sqrt(1-r^2).N2
*/
FloatType r_;
/* q := sqrt(1-r^2) */
FloatType q_;
/* state for generating indept normally-distributed r.v's */
std::normal_distribution<FloatType> ndist_;
}; /*gaussianpair_dist*/
/* generate pairs of correlated gaussian random variables */
template<class Engine>
class gaussianpairgen {
public:
using engine_type = Engine;
using generator_type = generator<Engine, gaussianpair_dist<double>>;
template<typename Engine>
static generator_type make(Engine eng,
double rho) {
return generator_type::make(std::move(eng),
gaussianpair_dist<double>(rho));
}
}; /*GaussianPairGen*/
} /*namespace random*/
} /*namespace xo*/
/* end gaussianpairgen.hpp */

View file

@ -0,0 +1,49 @@
/* @file generator.hpp */
#pragma once
#include "engine_concept.hpp"
#include "distribution_concept.hpp"
#include <utility>
namespace xo {
namespace rng {
/* Engine: uniform integer random number generator, e.g. xoshiro256ss
* Distribution: random number distribution, e.g. std::normal_distribution
*/
template <class Engine, class Distribution> requires engine_concept<Engine> && distribution_concept<Distribution, Engine>
class generator {
public:
using result_type = typename Distribution::result_type;
using engine_type = Engine;
public:
generator(Engine & e, Distribution const & d)
: engine_{e},
distribution_{d} {}
generator(Engine && e, Distribution && d)
: engine_{std::move(e)},
distribution_{std::move(d)} {}
static generator make(Engine e, Distribution d) {
return generator(e, d);
}
result_type operator()() { return this->distribution_(this->engine_); }
private:
/* random number generator; generates uniformly-distributed integers */
Engine engine_;
/* distribution object */
Distribution distribution_;
}; /*generator*/
template <class Engine, class Distribution>
generator<Engine, Distribution> make_generator(Engine e, Distribution d) {
return generator<Engine, Distribution>::make(std::move(e),
std::move(d));
} /*make_generator*/
} /*namespace rng*/
} /*namespace xo*/
/* end generator.hpp */

View file

@ -0,0 +1,14 @@
/* @file normalgen.hpp */
#pragma once
#include "generator.hpp"
#include <random>
namespace xo {
namespace rng {
/* Engine: e.g. xo::rng::xoshiro256 or std::mt19937 */
template <class Engine>
using normalgen = generator<Engine, std::normal_distribution<double>>;
} /*namespace rng*/
} /*namespace xo*/

View file

@ -0,0 +1,7 @@
/* @file print.hpp */
#pragma once
#include "xo/indentlog/print/array.hpp"
/* end print.hpp */

View file

@ -0,0 +1,87 @@
/* @file random_seed.hpp */
#include <iostream>
#include <cstdint>
#include <stdlib.h>
#ifdef _BSD_SOURCE
# include <bsd/stdlib.h>
#else
# include <sys/random.h>
#endif
namespace xo {
namespace rng {
/* generate a 64-bit random seed using /dev/urandom or similar source.
* This is relatively expensive; at least cost of a system call
* + may block if host has rebooted recently
*
* Require:
* - T is null-constructible.
*
* return value will contain a T-instance in which representation
* has been populated with random bits. Expecting T to be something
* like int32_t, or std::array<uint64_t, ..>
*/
template<typename T>
void random_seed(T * p_seed) {
# ifdef __APPLE__
/* NOTE: arc4random_buf() works on darwin/nix;
* probably need to do something else on intel linux
*/
::arc4random_buf(p_seed, sizeof(*p_seed));
# else
/* avail flags: GRND_RANDOM | GRND_NONBLOCK */
while (::getrandom(p_seed, sizeof(*p_seed), 0) == -1) {
if (errno == EINTR) {
/* interrupted by signal, try again */
continue;
} else {
break;
}
}
# endif
} /*random_seed*/
template<typename T>
T random_seed() {
T retval;
random_seed(&retval);
return retval;
} /*random_seed*/
/* RAII-style random-number seed
*
* Usage:
*
* Seed<xoshiro256ss> seed;
*
* auto eng = xoshiro256ss(seed);
* or
* auto rng = UnitIntervalGen<xoshiro256ss>::make(seed);
*/
template<typename Engine>
struct Seed {
using seed_type = typename Engine::seed_type;
Seed() { random_seed(&seed_); }
operator seed_type const & () const { return seed_; }
seed_type seed_;
}; /*Seed*/
template<typename T>
inline std::ostream &
operator<<(std::ostream & os,
Seed<T> const & x)
{
/* NOTE: if compile error here, may want caller to #include [indentlog/print/vector.hpp] */
os << x.seed_;
return os;
} /*operator<<*/
} /*namespace rng*/
} /*namespace xo*/
/* end random_seed.hpp */

View file

@ -0,0 +1,33 @@
/* @file uniformgen.hpp */
#pragma once
#include "generator.hpp"
#include <random>
namespace xo {
namespace rng {
template <class Engine>
class uniformgen : public generator<Engine, std::uniform_real_distribution<double>> {
public:
using generator_type = generator<Engine, std::uniform_real_distribution<double>>;
/* named ctor idiom */
template <class Eng>
static generator_type unit(Eng eng) {
return make_generator(std::move(eng),
std::uniform_real_distribution<double>(0.0, 1.0));
}
/* named ctor idiom */
template <class Eng>
static generator_type interval(Eng eng, double lo, double hi) {
return make_generator(std::move(eng),
std::uniform_real_distribution<double>(lo, hi));
}
};
} /*namespace rng*/
} /*namespace xo*/
/* end uniformgen.hpp */

View file

@ -0,0 +1,169 @@
/* @file xoshiro256.hpp */
#pragma once
#include "engine_concept.hpp"
#include <iostream>
#include <array>
#include <limits>
#include <cstdint>
namespace xo {
namespace rng {
/* engine for producing 64-bit random numbers
*
* see https:/en.wikipedia.org/wiki/Xorshift#xoshiro256**
*
* - satisfies c++ UniformRandomBitGenerator
* - satisfies c++
*
* Note: zero seed --> constant output sequence {0, 0, 0, ...}
*/
class xoshiro256ss {
public:
using result_type = std::uint64_t;
using seed_type = std::array<std::uint64_t, 4>;
public:
/* null state -- generates constant stream of 0 bits */
xoshiro256ss() : xoshiro256ss(0) {}
/* copy ctor */
xoshiro256ss(xoshiro256ss const & x) = default;
xoshiro256ss(seed_type const & seed) : s_(seed) {}
/* fallback version -- deprecated */
xoshiro256ss(std::uint64_t seed)
{
this->s_[0] = 0;
this->s_[1] = seed;
this->s_[2] = 0;
this->s_[3] = 0;
generate();
}
static constexpr std::uint64_t min() { return 0; }
static constexpr std::uint64_t max() { return std::numeric_limits<std::uint64_t>::max(); }
static std::uint64_t rol64(std::uint64_t x, std::int64_t k)
{
return (x << k) | (x >> (64 - k));
}
static bool equal(xoshiro256ss const & x, xoshiro256ss const & y) {
return ((x.s_[0] == y.s_[0])
&& (x.s_[1] == y.s_[1])
&& (x.s_[2] == y.s_[2])
&& (x.s_[3] == y.s_[3]));
}
/* puts generator into null state */
void seed() { *this = xoshiro256ss(); }
void seed(std::uint64_t s) { *this = xoshiro256ss{s}; }
/* e.g. used with std::seed_seq<> */
template <typename SeedSeq>
void seed(SeedSeq & sseq) {
sseq.generate(s_.begin(), s_.end());
}
std::uint64_t generate() {
std::array<std::uint64_t, 4> & s = (this->s_);
std::uint64_t const result = rol64(s[1] * 5, 7) * 9;
std::uint64_t const t = s[1] << 17;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rol64(s[3], 45);
return result;
} /*generate*/
/* advance to same state as obtained from z calls to .generate(). O(z) !
* usually better to use jump().
*
* providing .discard() to satisfy c++ named requirement _RandomNumberEngine_
*/
void discard(std::uint64_t z) {
for (std::uint64_t i=0; i<z; ++i)
this->generate();
}
/* equivalent to .discard(2^128), but uses O(1) time
*
* (may use in multithreaded program to get determinsitic non-overlapping random sequences)
*/
void jump() {
std::array<std::uint64_t, 4> const s_jump_v
= {{0x180ec6d33cfd0aba,
0xd5a61266f0c9392c,
0xa9582618e03fc9aa,
0x39abdc4529b1661c}};
std::array<std::uint64_t, 4> & s = (this->s_);
std::uint64_t s0 = 0;
std::uint64_t s1 = 0;
std::uint64_t s2 = 0;
std::uint64_t s3 = 0;
for (std::uint32_t i = 0; i < s_jump_v.size(); ++i) {
for (std::uint32_t bit = 0; bit < 64; ++bit) {
if (s_jump_v[i] & 1UL << bit) {
s0 ^= s[0];
s1 ^= s[1];
s2 ^= s[2];
s3 ^= s[3];
}
this->generate();
}
}
s[0] = s0;
s[1] = s1;
s[2] = s2;
s[3] = s3;
} /*jump*/
/* inverse of .load() */
void print(std::ostream & os) const {
os << "<xoshiro256ss " << s_[0] << " " << s_[1] << " " << s_[2] << " " << s_[3] << ">";
}
/* inverse of .print() */
void load(std::istream & is) {
std::string header, trailer;
std::array<std::uint64_t, 4> sv;
is >> header >> sv[0] >> sv[1] >> sv[2] >> sv[3] >> trailer;
if ((header != "<xoshiro256ss") || trailer != ">")
throw std::runtime_error("xoshiro256ss.load: bad input format, expecting input like <xoshiro256ss $s0 $s1 $s2 $s3>");
this->s_ = sv;
} /*load*/
std::uint64_t operator()() { return generate(); }
private:
/* state */
std::array<std::uint64_t, 4> s_;
}; /*xoshiro256ss*/
inline bool operator==(xoshiro256ss const & x, xoshiro256ss const & y) {
return xoshiro256ss::equal(x, y);
}
inline bool operator!=(xoshiro256ss const & x, xoshiro256ss const & y) {
return !xoshiro256ss::equal(x, y);
}
static_assert(engine_concept<xoshiro256ss>);
} /*namespace rng*/
} /*namespace xo*/
/* end xoshiro256.hpp */