Read /dev/urandom to init RNG for uniformrandom and gaussZiggurat

openair1/SIMULATION/TOOLS/rangen_double.c

@@ -26,42 +26,57 @@
 
 #include  "sim.h"
 
-static unsigned int seed, iy, ir[98];
-/*
-@defgroup _uniformdouble
-@ingroup numerical Uniform linear congruential random number generator.
-*/
-
-/*!\brief Initialization routine for Uniform/Gaussian random number generators. */
+static unsigned int urseed, iy, ir[98]; /// uniformrandom
+static bool tableNordDone = false; /// gaussZiggurat
 
 #define a 1664525lu
 #define mod 4294967296.0                /* is 2**32 */
 
-#if 1
-void randominit(unsigned seed_init)
+/*!\brief Generate a random number form `/dev/urandom`. */
+unsigned long get_random_seed(void)
 {
-  int i;
-  // this need to be integrated with the existing rng, like taus: navid
-
-  if (seed_init == 0) {
-    srand((unsigned)time(NULL));
-
-    seed = (unsigned int) rand();
-  } else {
-    seed = seed_init;
+  unsigned long seed;
+  const char* fn = "/dev/urandom";
+  FILE* f = fopen(fn, "rb");
+  if (f == NULL) {
+    fprintf(stderr, "could not open %s for seed generation: %d %s\n", fn, errno, strerror(errno));
+    abort();
+  }
+  int rc = fread(&seed, sizeof(seed), 1, f);
+  if (rc < 0) {
+    fprintf(stderr, "could not read %s for seed generation: %d %s\n", fn, errno, strerror(errno));
+    abort();
   }
-  printf("Initializing random number generator, seed %x\n",seed);
+  fclose(f);
+  return seed;
+}
 
-  if (seed % 2 == 0) seed += 1; /* seed and mod are relative prime */
 
-  for (i=1; i<=97; i++) {
-    seed = a*seed;                 /* mod 2**32  */
-    ir[i]= seed;                   /* initialize the shuffle table    */
+/*!\brief Initialization routine for Uniform/Gaussian random number generators. */
+void randominit(unsigned long seed_init)
+{
+  unsigned long seed = seed_init != 0 ? seed_init : get_random_seed();
+  printf("Initializing random number generator, seed %ld\n", seed);
+
+  // initialize uniformrandom RNG
+  urseed = (unsigned int) seed;
+  if (urseed % 2 == 0)
+    urseed += 1; /* urseed and mod are relative prime */
+  for (int i = 1; i <= 97; i++) {
+    urseed = a * urseed; /* mod 2**32  */
+    ir[i] = urseed; /* initialize the shuffle table    */
   }
+  iy = 1;
 
- iy=1;
+  // initialize gaussZiggurat RNG
+  tableNor(seed);
+  tableNordDone = true;
 }
-#endif
+
+/*
+@defgroup _uniformdouble
+@ingroup numerical Uniform linear congruential random number generator.
+*/
 
 /*!\brief Uniform linear congruential random number generator on \f$[0,1)\f$.  Returns a double-precision floating-point number.*/
 
@@ -70,13 +85,11 @@ double uniformrandom(void)
 #define a 1664525lu
 #define mod 4294967296.0                /* is 2**32 */
 
-  int j;
-
-  j=1 + 97.0*iy/mod;
-  iy=ir[j];
-  seed = a*seed;                          /* mod 2**32 */
-  ir[j] = seed;
-  return( (double) iy/mod );
+  int j = 1 + 97.0 * iy / mod;
+  iy = ir[j];
+  urseed = a * urseed; /* mod 2**32 */
+  ir[j] = urseed;
+  return (double)iy / mod;
 }
 
 /*
@@ -109,8 +122,12 @@ double __attribute__ ((no_sanitize_address)) gaussdouble(double mean, double var
   }
 }
 
+
+/*
+@defgroup _gaussZiggurat
+@ingroup numerical ziggurat random number generator for exponentially distributed numbers
+*/
 // Ziggurat
-static bool tableNordDone=false;
 static double wn[128], fn[128];
 static uint32_t iz, jz, jsr = 123456789, kn[128];
 static int32_t hz;
@@ -180,9 +197,7 @@ double __attribute__ ((no_sanitize_address)) gaussZiggurat(double mean, double v
 {
   if (!tableNordDone) {
     // let's make reasonnable constant tables
-    struct timespec t;
-    clock_gettime(CLOCK_MONOTONIC, &t);
-    tableNor((long)(t.tv_nsec%INT_MAX));
+    tableNor(get_random_seed());
   }
   hz = SHR3;
   iz = hz & 127;

openair1/SIMULATION/TOOLS/sim.h

@@ -523,8 +523,9 @@ the value \f$\mathrm{sgn}(u)i\f$.  The search requires at most \f$Nbits-1\f$ com
 */
 int gauss(unsigned int *gauss_LUT,unsigned char Nbits);
 
+unsigned long get_random_seed(void);
 double gaussdouble(double,double);
-void randominit(unsigned int seed_init);
+void randominit(unsigned long seed_init);
 double uniformrandom(void);
 double gaussZiggurat(double mean, double variance);
 void tableNor(unsigned long seed);

openair2/UTIL/MATH/random.c

-/*
- * Licensed to the OpenAirInterface (OAI) Software Alliance under one or more
- * contributor license agreements.  See the NOTICE file distributed with
- * this work for additional information regarding copyright ownership.
- * The OpenAirInterface Software Alliance licenses this file to You under
- * the OAI Public License, Version 1.1  (the "License"); you may not use this file
- * except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *      http://www.openairinterface.org/?page_id=698
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- *-------------------------------------------------------------------------------
- * For more information about the OpenAirInterface (OAI) Software Alliance:
- *      contact@openairinterface.org
- */
-
-/*
-                                 random.c
-                             -------------------
-  AUTHOR  : Lionel GAUTHIER
-  COMPANY : EURECOM
-  EMAIL   : Lionel.Gauthier@eurecom.fr
- ***************************************************************************/
-#include "rtos_header.h"
-#include "platform_types.h"
-
-#include <sys/time.h>
-
-
-/* Random generators */
-#define FACTOR       16807
-#define LASTXN       127773
-#define UPTOMOD     -2836
-
-static int      seed;
-
-
-void
-init_uniform (void)
-{
-  struct timeval  tv;
-  struct timezone tz;
-
-  gettimeofday (&tv, &tz);
-  seed = (int) tv.tv_usec;
-#ifdef NODE_MT
-#warning TO DO seed = mobileId
-  //seed += mobileId;
-#endif
-#ifdef NODE_RG
-#warning TO DO seed = rgId
-  //seed += rgId;
-#endif
-}
-
-
-int
-uniform (void)
-{
-  static int      times, rest, prod1, prod2;
-
-  times = seed / LASTXN;
-  rest = seed - times * LASTXN;
-  prod1 = times * UPTOMOD;
-  prod2 = rest * FACTOR;
-  seed = prod1 + prod2;
-
-  return seed;
-}