star-sp

Random number generators and distributions
git clone git://git.meso-star.fr/star-sp.git
Log | Files | Refs | README | LICENSE

ssp_rng.c (18482B)


      1 /* Copyright (C) 2015-2025 |Méso|Star> (contact@meso-star.com)
      2  *
      3  * This program is free software: you can redistribute it and/or modify
      4  * it under the terms of the GNU General Public License as published by
      5  * the Free Software Foundation, either version 3 of the License, or
      6  * (at your option) any later version.
      7  *
      8  * This program is distributed in the hope that it will be useful,
      9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
     10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     11  * GNU General Public License for more details.
     12  *
     13  * You should have received a copy of the GNU General Public License
     14  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
     15 
     16 #include "ssp_rng_c.h"
     17 
     18 #include <rsys/dynamic_array_char.h>
     19 #include <rsys/mem_allocator.h>
     20 
     21 #include <sstream>
     22 #include <cstring>
     23 #include <limits>
     24 
     25 /*******************************************************************************
     26  * KISS PRNG
     27  ******************************************************************************/
     28 struct rng_kiss { uint32_t x, y, z, c; };
     29 
     30 static res_T
     31 rng_kiss_set(void* data, const uint64_t seed)
     32 {
     33   struct rng_kiss* kiss = (struct rng_kiss*)data;
     34   RAN_NAMESPACE::mt19937 rng_mt(uint32_t(seed % UINT32_MAX));
     35   ASSERT(kiss);
     36   kiss->x = (uint32_t)rng_mt();
     37   do {
     38     kiss->y = (uint32_t)rng_mt();
     39   } while(kiss->y == 0); /* Y must be != 0 */
     40   kiss->z = (uint32_t)rng_mt();
     41   /* Offset c by 1 to avoid z=c=0; should be less than 698769069 */
     42   kiss->c = (uint32_t)rng_mt() % 698769068 + 1;
     43   return RES_OK;
     44 }
     45 
     46 static uint64_t
     47 rng_kiss_get(void* data)
     48 {
     49   struct rng_kiss* kiss = (struct rng_kiss*)data;
     50   uint64_t t;
     51   ASSERT(kiss);
     52 
     53   kiss->x = 314527869 * kiss->x + 1234567;
     54   kiss->y ^= kiss->y << 5;
     55   kiss->y ^= kiss->y >> 7;
     56   kiss->y ^= kiss->y << 22;
     57   t = 4294584393ULL * kiss->z + kiss->c;
     58   kiss->c = (uint32_t)(t >> 32);
     59   kiss->z = (uint32_t)(t & 0xFFFFFFFF);
     60   return (uint32_t)(kiss->x + kiss->y + kiss->z);
     61 }
     62 
     63 static res_T
     64 rng_kiss_read(void* data, FILE* file)
     65 {
     66   uint32_t val[4];
     67   struct rng_kiss* rng = (struct rng_kiss*)data;
     68   size_t n;
     69   ASSERT(data && file);
     70   n = fread(val, sizeof(*val), sizeof(val)/sizeof(*val), file);
     71   if(n != 4) return RES_IO_ERR;
     72   rng->x = val[0];
     73   rng->y = val[1];
     74   rng->z = val[2];
     75   rng->c = val[3];
     76   return RES_OK;
     77 }
     78 
     79 static res_T
     80 rng_kiss_read_cstr(void* data, const char* cstr)
     81 {
     82   struct rng_kiss* rng = (struct rng_kiss*)data;
     83   uint32_t val[4];
     84   size_t sz;
     85   ASSERT(data && cstr);
     86   sz = strlen(cstr);
     87   if(sz < sizeof(val)) return RES_IO_ERR;
     88   memcpy(val, cstr, sizeof(val));
     89   rng->x = val[0];
     90   rng->y = val[1];
     91   rng->z = val[2];
     92   rng->c = val[3];
     93   return RES_OK;
     94 }
     95 
     96 static res_T
     97 rng_kiss_write(const void* data, FILE* file)
     98 {
     99   const struct rng_kiss* rng = (const struct rng_kiss*)data;
    100   uint32_t val[4];
    101   size_t n;
    102   ASSERT(data && file);
    103   val[0] = rng->x;
    104   val[1] = rng->y;
    105   val[2] = rng->z;
    106   val[3] = rng->c;
    107   n = fwrite(val, sizeof(*val), sizeof(val)/sizeof(*val), file);
    108   return (n != 4 ? RES_IO_ERR : RES_OK);
    109 }
    110 
    111 static res_T
    112 rng_kiss_write_cstr
    113   (const void* data,
    114    char* buf,
    115    const size_t bufsz,
    116    size_t* out_len)
    117 {
    118   const struct rng_kiss* rng = (const struct rng_kiss*)data;
    119   uint32_t val[4];
    120   ASSERT(data);
    121 
    122   val[0] = rng->x;
    123   val[1] = rng->y;
    124   val[2] = rng->z;
    125   val[3] = rng->c;
    126   if(bufsz >= sizeof(val) + 1/*null char*/) {
    127     memcpy(buf, val, sizeof(val));
    128     buf[sizeof(val)] = '\0';
    129   }
    130   if(out_len) *out_len = sizeof(val);
    131   return RES_OK;
    132 }
    133 
    134 static res_T
    135 rng_kiss_init(struct mem_allocator* allocator, void* data)
    136 {
    137   (void)allocator;
    138   rng_kiss_set(data, 0);
    139   return RES_OK;
    140 }
    141 
    142 static void
    143 rng_kiss_release(void* data)
    144 {
    145   (void)data;
    146 }
    147 
    148 static double
    149 rng_kiss_entropy(const void* data)
    150 {
    151   (void)data;
    152   return 0.;
    153 }
    154 
    155 static res_T
    156 rng_kiss_discard(void* data, uint64_t n)
    157 {
    158   while (n-- > 0) {
    159     rng_kiss_get(data);
    160   }
    161   return RES_OK;
    162 }
    163 
    164 /* Exported type */
    165 static const struct rng_desc rng_kiss = {
    166   rng_kiss_init,
    167   rng_kiss_release,
    168   rng_kiss_set,
    169   rng_kiss_get,
    170   rng_kiss_discard,
    171   rng_kiss_read,
    172   rng_kiss_read_cstr,
    173   rng_kiss_write,
    174   rng_kiss_write_cstr,
    175   rng_kiss_entropy,
    176   0,
    177   UINT32_MAX,
    178   sizeof(struct rng_kiss),
    179   16
    180 };
    181 
    182 /*******************************************************************************
    183  * C++11 PRNG
    184  ******************************************************************************/
    185 template<typename RNG>
    186 static res_T
    187 rng_cxx_set(void* data, const uint64_t seed)
    188 {
    189   RNG* rng = (RNG*)data;
    190   ASSERT(rng);
    191   rng->seed((typename RNG::result_type)seed);
    192   return RES_OK;
    193 }
    194 
    195 template<>
    196 res_T
    197 rng_cxx_set<RAN_NAMESPACE::random_device>(void* data, const uint64_t seed)
    198 {
    199   (void) data; (void) seed;
    200   return RES_BAD_OP;
    201 }
    202 
    203 template<typename RNG>
    204 static uint64_t
    205 rng_cxx_get(void* data)
    206 {
    207   RNG* rng = (RNG*)data;
    208   ASSERT(rng);
    209   return (*rng)();
    210 }
    211 
    212 template<typename RNG>
    213 static res_T
    214 rng_cxx_write(const void* data, FILE* file)
    215 {
    216   size_t i;
    217   std::stringstream stream;
    218   RNG* rng = (RNG*)data;
    219   ASSERT(rng);
    220   stream << *rng << std::endl;
    221   i = fwrite(stream.str().c_str(), 1, stream.str().size(), file);
    222   return i == stream.str().size() ? RES_OK : RES_IO_ERR;
    223 }
    224 
    225 template<>
    226 res_T
    227 rng_cxx_write<RAN_NAMESPACE::random_device>(const void* data, FILE* file)
    228 {
    229   (void) data; (void) file;
    230   return RES_BAD_OP;
    231 }
    232 
    233 template<typename RNG>
    234 static res_T
    235 rng_cxx_write_cstr
    236   (const void* data,
    237    char* buf,
    238    const size_t bufsz,
    239    size_t* out_len)
    240 {
    241   int len = 0;
    242   std::stringstream stream;
    243   RNG* rng = (RNG*)data;
    244   ASSERT(rng);
    245   stream << *rng << std::endl;
    246   len = snprintf(buf, bufsz, "%s", stream.str().c_str());
    247   CHK(len > 0);
    248   if((size_t)len >= (bufsz - 1/*null char*/)) {
    249     buf[bufsz-1] = '\0';
    250   }
    251   if(out_len) *out_len = (size_t)len;
    252   return RES_OK;
    253 }
    254 
    255 template<>
    256 res_T
    257 rng_cxx_write_cstr<RAN_NAMESPACE::random_device>
    258   (const void* data, char* buf, const size_t bufsz, size_t* out_len)
    259 {
    260   (void)data; (void)buf, (void)bufsz, (void)out_len;
    261   return RES_BAD_OP;
    262 }
    263 
    264 template<typename RNG>
    265 static res_T
    266 rng_cxx_read(void* data, FILE* file)
    267 {
    268   std::stringstream stream;
    269   char buf[512];
    270   char* s;
    271   RNG* rng = (RNG*)data;
    272   ASSERT(rng);
    273   while((s = fgets(buf, (int)sizeof(buf), file))) {
    274     stream << std::string(s);
    275     if(s[strlen(s)-1] == '\n') break;
    276   }
    277   stream >> *rng;
    278   return stream.fail() ? RES_IO_ERR : RES_OK;
    279 }
    280 
    281 template<>
    282 res_T
    283 rng_cxx_read<RAN_NAMESPACE::random_device>(void* data, FILE* file)
    284 {
    285   (void) data; (void)file;
    286   return RES_BAD_OP;
    287 }
    288 
    289 template<typename RNG>
    290 static res_T
    291 rng_cxx_read_cstr(void* data, const char* cstr)
    292 {
    293   std::stringstream stream;
    294   RNG* rng = (RNG*)data;
    295   ASSERT(rng && cstr && cstr[strlen(cstr)-1] == '\n');
    296   stream << std::string(cstr);
    297   stream >> *rng;
    298   return stream.fail() ? RES_IO_ERR : RES_OK;
    299 }
    300 
    301 template<>
    302 res_T
    303 rng_cxx_read_cstr<RAN_NAMESPACE::random_device>(void* data, const char* cstr)
    304 {
    305   (void)data, (void)cstr;
    306   return RES_BAD_OP;
    307 }
    308 
    309 template<typename RNG>
    310 static res_T
    311 rng_cxx_init(struct mem_allocator* allocator, void* data)
    312 {
    313   (void)allocator;
    314   ASSERT(data);
    315   new (data) RNG;
    316   return RES_OK;
    317 }
    318 
    319 template<typename RNG>
    320 static void
    321 rng_cxx_release(void* data)
    322 {
    323   RNG* rng = (RNG*)data;
    324   (void)rng;
    325   ASSERT(rng);
    326   rng->~RNG();
    327 }
    328 
    329 template<typename RNG>
    330 static double
    331 rng_cxx_entropy(const void* data)
    332 {
    333   (void)data;
    334   return 0;
    335 }
    336 
    337 template<typename RNG>
    338 res_T
    339 rng_cxx_discard(void* data, uint64_t n)
    340 {
    341   RNG* rng = (RNG*) data;
    342   ASSERT(rng);
    343   rng->discard(n);
    344   return RES_OK;
    345 }
    346 
    347 template<>
    348 res_T
    349 rng_cxx_discard<RAN_NAMESPACE::random_device>(void* data, uint64_t n)
    350 {
    351   (void) data; (void) n;
    352   return RES_BAD_OP;
    353 }
    354 
    355 template<>
    356 double
    357 rng_cxx_entropy<RAN_NAMESPACE::random_device>(const void* data)
    358 {
    359   const RAN_NAMESPACE::random_device* rng =
    360     (const RAN_NAMESPACE::random_device*)data;
    361   ASSERT(rng);
    362   return rng->entropy();
    363 }
    364 
    365 /* 64-bits Mersenne Twister PRNG */
    366 static const struct rng_desc rng_mt19937_64 = {
    367   rng_cxx_init<RAN_NAMESPACE::mt19937_64>,
    368   rng_cxx_release<RAN_NAMESPACE::mt19937_64>,
    369   rng_cxx_set<RAN_NAMESPACE::mt19937_64>,
    370   rng_cxx_get<RAN_NAMESPACE::mt19937_64>,
    371   rng_cxx_discard<RAN_NAMESPACE::mt19937_64>,
    372   rng_cxx_read<RAN_NAMESPACE::mt19937_64>,
    373   rng_cxx_read_cstr<RAN_NAMESPACE::mt19937_64>,
    374   rng_cxx_write<RAN_NAMESPACE::mt19937_64>,
    375   rng_cxx_write_cstr<RAN_NAMESPACE::mt19937_64>,
    376   rng_cxx_entropy<RAN_NAMESPACE::mt19937_64>,
    377   RAN_NAMESPACE::mt19937_64::min(),
    378   RAN_NAMESPACE::mt19937_64::max(),
    379   sizeof(RAN_NAMESPACE::mt19937_64),
    380   16
    381 };
    382 
    383 /* 48-bits RANLUX PRNG */
    384 static const struct rng_desc rng_ranlux48 = {
    385   rng_cxx_init<RAN_NAMESPACE::ranlux48>,
    386   rng_cxx_release<RAN_NAMESPACE::ranlux48>,
    387   rng_cxx_set<RAN_NAMESPACE::ranlux48>,
    388   rng_cxx_get<RAN_NAMESPACE::ranlux48>,
    389   rng_cxx_discard<RAN_NAMESPACE::ranlux48>,
    390   rng_cxx_read<RAN_NAMESPACE::ranlux48>,
    391   rng_cxx_read_cstr<RAN_NAMESPACE::ranlux48>,
    392   rng_cxx_write<RAN_NAMESPACE::ranlux48>,
    393   rng_cxx_write_cstr<RAN_NAMESPACE::ranlux48>,
    394   rng_cxx_entropy<RAN_NAMESPACE::ranlux48>,
    395   RAN_NAMESPACE::ranlux48::min(),
    396   RAN_NAMESPACE::ranlux48::max(),
    397   sizeof(RAN_NAMESPACE::ranlux48),
    398   16
    399 };
    400 
    401 /* random_device generator */
    402 static const struct rng_desc rng_random_device = {
    403   rng_cxx_init<RAN_NAMESPACE::random_device>,
    404   rng_cxx_release<RAN_NAMESPACE::random_device>,
    405   rng_cxx_set<RAN_NAMESPACE::random_device>,
    406   rng_cxx_get<RAN_NAMESPACE::random_device>,
    407   rng_cxx_discard<RAN_NAMESPACE::random_device>,
    408   rng_cxx_read<RAN_NAMESPACE::random_device>,
    409   rng_cxx_read_cstr<RAN_NAMESPACE::random_device>,
    410   rng_cxx_write<RAN_NAMESPACE::random_device>,
    411   rng_cxx_write_cstr<RAN_NAMESPACE::random_device>,
    412   rng_cxx_entropy<RAN_NAMESPACE::random_device>,
    413   RAN_NAMESPACE::random_device::min(),
    414   RAN_NAMESPACE::random_device::max(),
    415   sizeof(RAN_NAMESPACE::random_device),
    416   16
    417 };
    418 
    419 /*******************************************************************************
    420  * Random123 Counter Based RNG
    421  ******************************************************************************/
    422 typedef r123::Engine<r123::Threefry4x64> threefry_T;
    423 
    424 #ifdef WITH_R123_AES
    425 typedef r123::Engine<r123::AESNI4x32> aes_T;
    426 
    427 template<>
    428 res_T
    429 rng_cxx_init<aes_T>(struct mem_allocator* allocator, void* data)
    430 {
    431   (void) allocator;
    432 
    433   if(!haveAESNI()) {
    434     /* AES-NI instructions not available on this hardware */
    435     return RES_BAD_OP;
    436   }
    437   ASSERT(data);
    438   new (data) aes_T;
    439   return RES_OK;
    440 }
    441 #endif
    442 
    443 /* threefry generator */
    444 static const struct rng_desc rng_threefry = {
    445   rng_cxx_init<threefry_T>,
    446   rng_cxx_release<threefry_T>,
    447   rng_cxx_set<threefry_T>,
    448   rng_cxx_get<threefry_T>,
    449   rng_cxx_discard<threefry_T>,
    450   rng_cxx_read<threefry_T>,
    451   rng_cxx_read_cstr<threefry_T>,
    452   rng_cxx_write<threefry_T>,
    453   rng_cxx_write_cstr<threefry_T>,
    454   rng_cxx_entropy<threefry_T>,
    455   threefry_T::min(),
    456   threefry_T::max(),
    457   sizeof(threefry_T),
    458   16
    459 };
    460 
    461 #ifdef WITH_R123_AES
    462 /* aes generator */
    463 static const struct rng_desc rng_aes = {
    464   rng_cxx_init<aes_T>,
    465   rng_cxx_release<aes_T>,
    466   rng_cxx_set<aes_T>,
    467   rng_cxx_get<aes_T>,
    468   rng_cxx_discard<aes_T>,
    469   rng_cxx_read<aes_T>,
    470   rng_cxx_read_cstr<aes_T>,
    471   rng_cxx_write<aes_T>,
    472   rng_cxx_write_cstr<aes_T>,
    473   rng_cxx_entropy<aes_T>,
    474   aes_T::min(),
    475   aes_T::max(),
    476   sizeof(aes_T),
    477   16
    478 };
    479 #endif
    480 
    481 /*******************************************************************************
    482  * Helper functions
    483  ******************************************************************************/
    484 static void
    485 rng_release(ref_T* ref)
    486 {
    487   struct ssp_rng* rng;
    488   ASSERT(ref);
    489   rng = CONTAINER_OF(ref, struct ssp_rng, ref);
    490   if(rng->state) {
    491     if(rng->desc.release) {
    492       rng->desc.release(rng->state);
    493     }
    494     MEM_RM(rng->allocator, rng->state);
    495   }
    496   MEM_RM(rng->allocator, rng);
    497 }
    498 
    499 /*******************************************************************************
    500  * Exported function
    501  ******************************************************************************/
    502 res_T
    503 ssp_rng_create
    504   (struct mem_allocator* mem_allocator,
    505    const enum ssp_rng_type type,
    506    struct ssp_rng** out_rng)
    507 {
    508   res_T res = RES_OK;
    509 
    510   switch(type) {
    511     case SSP_RNG_KISS:
    512       res = rng_create(mem_allocator, &rng_kiss, out_rng);
    513       break;
    514     case SSP_RNG_MT19937_64:
    515       res = rng_create(mem_allocator, &rng_mt19937_64, out_rng);
    516       break;
    517     case SSP_RNG_RANLUX48:
    518       res = rng_create(mem_allocator, &rng_ranlux48, out_rng);
    519       break;
    520     case SSP_RNG_RANDOM_DEVICE:
    521       res = rng_create(mem_allocator, &rng_random_device, out_rng);
    522       break;
    523     case SSP_RNG_THREEFRY:
    524       res = rng_create(mem_allocator, &rng_threefry, out_rng);
    525       break;
    526 #ifdef WITH_R123_AES
    527     case SSP_RNG_AES:
    528       res = rng_create(mem_allocator, &rng_aes, out_rng);
    529       break;
    530 #endif
    531     default: res = RES_BAD_ARG; break;
    532   }
    533   if(res != RES_OK) goto error;
    534   (*out_rng)->type = type;
    535 
    536 exit:
    537   return res;
    538 error:
    539   goto exit;
    540 }
    541 
    542 res_T
    543 ssp_rng_ref_get(struct ssp_rng* rng)
    544 {
    545   if(!rng) return RES_BAD_ARG;
    546   ref_get(&rng->ref);
    547   return RES_OK;
    548 }
    549 
    550 res_T
    551 ssp_rng_ref_put(struct ssp_rng* rng)
    552 {
    553   if(!rng) return RES_BAD_ARG;
    554   ref_put(&rng->ref, rng_release);
    555   return RES_OK;
    556 }
    557 
    558 res_T
    559 ssp_rng_get_type(const struct ssp_rng* rng, enum ssp_rng_type* type)
    560 {
    561   if(!rng || !type) return RES_BAD_ARG;
    562   *type = rng->type;
    563   return RES_OK;
    564 }
    565 
    566 uint64_t
    567 ssp_rng_get(struct ssp_rng* rng)
    568 {
    569   if(!rng) FATAL("The Random Number Generator is NULL\n");
    570   return rng->desc.get(rng->state);
    571 }
    572 
    573 uint64_t
    574 ssp_rng_uniform_uint64
    575   (struct ssp_rng* rng, const uint64_t lower, const uint64_t upper)
    576 {
    577   if(!rng) FATAL("The Random Number Generator is NULL\n");
    578   RAN_NAMESPACE::uniform_int_distribution<uint64_t> distrib(lower, upper);
    579   return wrap_ran(*rng, distrib);
    580 }
    581 
    582 double
    583 ssp_rng_uniform_double
    584   (struct ssp_rng* rng, const double lower, const double upper)
    585 {
    586   if(!rng) FATAL("The Random Number Generator is NULL\n");
    587   RAN_NAMESPACE::uniform_real_distribution<double> distrib(lower, upper);
    588   return wrap_ran(*rng, distrib);
    589 }
    590 
    591 float
    592 ssp_rng_uniform_float
    593   (struct ssp_rng* rng, const float lower, const float upper)
    594 {
    595   if(!rng) FATAL("The Random Number Generator is NULL\n");
    596   RAN_NAMESPACE::uniform_real_distribution<float> distrib(lower, upper);
    597   return wrap_ran(*rng, distrib);
    598 }
    599 
    600 double
    601 ssp_rng_canonical(struct ssp_rng* rng)
    602 {
    603   if(!rng) FATAL("The Random Number Generator is NULL\n");
    604 #if 0
    605   rng_cxx rng_cxx(*rng);
    606   return RAN_NAMESPACE::generate_canonical
    607     <double, std::numeric_limits<double>::digits>(rng_cxx);
    608 #else
    609   /* Optimized version */
    610   size_t k = rng->dbl_k;
    611   double sum = 0;
    612   double tmp = 1;
    613   for(; k !=0; --k) {
    614     sum += (double)(rng->desc.get(rng->state) - rng->desc.min) * tmp;
    615     tmp = (double)(tmp*rng->r);
    616   }
    617   return sum/tmp;
    618 #endif
    619 }
    620 
    621 float
    622 ssp_rng_canonical_float(struct ssp_rng* rng)
    623 {
    624   if(!rng) FATAL("The Random Number Generator is NULL\n");
    625   /* The C++ standard library does not ensure that the generated single
    626    * precision floating point number is effectively canonical, i.e. it may be
    627    * equal to 1. Use a hand made workaround to handle this bug. */
    628 #if 0
    629   rng_cxx rng_cxx(*rng);
    630   return RAN_NAMESPACE::generate_canonical
    631     <float, std::numeric_limits<float>::digits>(rng_cxx);
    632 #else
    633   float r;
    634   do {
    635     /* Optimised version of the generate_canonical function */
    636     size_t k = rng->flt_k;
    637     float sum = 0;
    638     float tmp = 1;
    639     for(; k !=0; --k) {
    640       sum += (float)(rng->desc.get(rng->state) - rng->desc.min) * tmp;
    641       tmp = (float)(tmp*rng->r);
    642     }
    643     r = sum/tmp;
    644   } while(r >= 1);
    645   return r;
    646 #endif
    647 }
    648 
    649 uint64_t
    650 ssp_rng_min(struct ssp_rng* rng)
    651 {
    652   if(!rng) FATAL("The Random Number Generator is NULL\n");
    653   return rng->desc.min;
    654 }
    655 
    656 uint64_t
    657 ssp_rng_max(struct ssp_rng* rng)
    658 {
    659   if(!rng) FATAL("The Random Number Generator is NULL\n");
    660   return rng->desc.max;
    661 }
    662 
    663 res_T
    664 ssp_rng_set(struct ssp_rng* rng, const uint64_t seed)
    665 {
    666   if(!rng) return RES_BAD_ARG;
    667   return rng->desc.set(rng->state, seed);
    668 }
    669 
    670 res_T
    671 ssp_rng_read(struct ssp_rng* rng, FILE* stream)
    672 {
    673   if(!rng || !stream) return RES_BAD_ARG;
    674   return rng->desc.read(rng->state, stream);
    675 }
    676 
    677 res_T
    678 ssp_rng_read_cstr(struct ssp_rng* rng, const char* cstr)
    679 {
    680   if(!rng || !cstr) return RES_BAD_ARG;
    681   return rng->desc.read_cstr(rng->state, cstr);
    682 }
    683 
    684 res_T
    685 ssp_rng_write(const struct ssp_rng* rng, FILE* stream)
    686 {
    687   if(!rng || !stream) return RES_BAD_ARG;
    688   return rng->desc.write(rng->state, stream);
    689 }
    690 
    691 res_T
    692 ssp_rng_write_cstr
    693   (const struct ssp_rng* rng,
    694    char* buf,
    695    const size_t bufsz,
    696    size_t* len)
    697 {
    698   if(!rng) return RES_BAD_ARG;
    699   return rng->desc.write_cstr(rng->state, buf, bufsz, len);
    700 }
    701 
    702 double
    703 ssp_rng_entropy(const struct ssp_rng* rng)
    704 {
    705   if (!rng) FATAL("The Random Number Generator is NULL\n");
    706   return rng->desc.entropy(rng->state);
    707 }
    708 
    709 res_T
    710 ssp_rng_discard(struct ssp_rng* rng, uint64_t n)
    711 {
    712   if (!rng) FATAL("The Random Number Generator is NULL\n");
    713   return rng->desc.discard(rng->state, n);
    714 }
    715 
    716 /*******************************************************************************
    717  * Local function
    718  ******************************************************************************/
    719 res_T
    720 rng_create
    721   (struct mem_allocator* mem_allocator,
    722    const struct rng_desc* desc,
    723    struct ssp_rng** out_rng)
    724 {
    725   struct mem_allocator* allocator;
    726   struct ssp_rng* rng = NULL;
    727   size_t log2r;
    728   res_T res = RES_OK;
    729 
    730   if(!desc || !out_rng) {
    731     res = RES_BAD_ARG;
    732     goto error;
    733   }
    734 
    735   allocator = mem_allocator ? mem_allocator : &mem_default_allocator;
    736 
    737   /* Align the rng on 16 bytes since its 'r' attrib is a long double that on
    738    * Linux 64-bits systems must be aligned on 16 bytes. */
    739   rng = (struct ssp_rng*)MEM_ALLOC_ALIGNED(allocator, sizeof(*rng), 16);
    740   if(!rng) {
    741     res = RES_MEM_ERR;
    742     goto error;
    743   }
    744   memset(rng, 0, sizeof(struct ssp_rng));
    745   rng->allocator = allocator;
    746   ref_init(&rng->ref);
    747   rng->desc = *desc;
    748 
    749   rng->state = MEM_ALLOC_ALIGNED
    750     (rng->allocator, rng->desc.sizeof_state, rng->desc.alignof_state);
    751   if(!rng->state) {
    752     res = RES_MEM_ERR;
    753     goto error;
    754   }
    755   res = rng->desc.init(allocator, rng->state);
    756   if(res != RES_OK) goto error;
    757 
    758   /* Precompute some values for the canonical distribution */
    759   rng->r = (long double)rng->desc.max - (long double)rng->desc.min + 1.0L;
    760   log2r = (size_t)(std::log(rng->r) / std::log(2.0L));
    761   rng->dbl_k = std::max<size_t>(1UL, (53 + log2r - 1UL) / log2r);
    762   rng->flt_k = std::max<size_t>(1UL, (24 + log2r - 1UL) / log2r);
    763 
    764 exit:
    765   if(out_rng) *out_rng = rng;
    766   return res;
    767 error:
    768   if(rng) {
    769     SSP(rng_ref_put(rng));
    770     rng = NULL;
    771   }
    772   goto exit;
    773 }