/* * 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 */ #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> #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. */ #define a 1664525lu #define mod 4294967296.0 /* is 2**32 */ #if 1 void randominit(unsigned seed_init) { int i; // this need to be integrated with the existing rng, like taus: navid printf("Initializing random number generator, seed %x\n",seed_init); if (seed_init == 0) { srand((unsigned)time(NULL)); seed = (unsigned int) rand(); } else { seed = seed_init; } 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 */ } iy=1; } #endif /*!\brief Uniform linear congruential random number generator on \f$[0,1)\f$. Returns a double-precision floating-point number.*/ 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 ); } /* @defgroup _gaussdouble Gaussian random number generator based on modified Box-Muller transformation. @ingroup numerical */ /*!\brief Gaussian random number generator based on modified Box-Muller transformation.Returns a double-precision floating-point number. */ double gaussdouble(double mean, double variance) { static int iset=0; static double gset; double fac,r,v1,v2; if (iset == 0) { do { v1 = 2.0*uniformrandom()-1.0; v2 = 2.0*uniformrandom()-1.0; r = v1*v1+v2*v2; } while (r >= 1.0); fac = sqrt(-2.0*log(r)/r); gset= v1*fac; iset=1; return(sqrt(variance)*v2*fac + mean); } else { iset=0; return(sqrt(variance)*gset + mean); } } #ifdef MAIN main(int argc,char **argv) { int i; randominit(); for (i=0; i<10; i++) { printf("%f\n",gaussdouble(0.0,1.0)); } } #endif