Commit 38f263cc authored by lfarizav's avatar lfarizav

New Ziggurat Gauss number generation function

parent 57517313
......@@ -442,8 +442,8 @@ clock_t start=clock();*/
else
{
//rx128_gain_lin=mm_mul_set1_ps(rx_gain_lin);
gauss0_sqrt_NOW=sqrt_NOW*gaussdouble(0.0,1.0);
gauss1_sqrt_NOW=sqrt_NOW*gaussdouble(0.0,1.0);
gauss0_sqrt_NOW=sqrt_NOW*ziggurat();
gauss1_sqrt_NOW=sqrt_NOW*ziggurat();
rx128_re = _mm_loadu_pd(&r_re[a][2*i]);//r_re[a][i],r_re[a][i+1]
rx128_im = _mm_loadu_pd(&r_im[a][2*i]);//r_im[a][i],r_im[a][i+1]
rx128_gain_lin = _mm_set1_pd(rx_gain_lin);
......
......@@ -436,9 +436,13 @@ 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);
double zigguratdouble(unsigned int r, double sigma);
double gaussdouble(double,double);
void randominit(unsigned int seed_init);
void setup_nor();
double uniformrandom(void);
void uniformrandomSSE(__m128d *d1,__m128d *d2);
double ziggurat();
int freq_channel(channel_desc_t *desc,uint16_t nb_rb, int16_t n_samples);
int freq_channel_prach(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples,int16_t prach_fmt,int16_t n_ra_prb);
int init_freq_channel(channel_desc_t *desc,uint16_t nb_rb,int16_t n_samples);
......
......@@ -28,120 +28,11 @@
#include "defs.h"
/*Ziggurat tables ytab, ktab, wtab are used in the random number generation routine*/
/*position of right-most step */
#define PARAM_R 3.44428647676
/* Tabulated values for the Ziggurat levels (heigt)*/
static const double ytab[128] = {
1, 0.963598623011, 0.936280813353, 0.913041104253,
0.892278506696, 0.873239356919, 0.855496407634, 0.838778928349,
0.822902083699, 0.807732738234, 0.793171045519, 0.779139726505,
0.765577436082, 0.752434456248, 0.739669787677, 0.727249120285,
0.715143377413, 0.703327646455, 0.691780377035, 0.68048276891,
0.669418297233, 0.65857233912, 0.647931876189, 0.637485254896,
0.62722199145, 0.617132611532, 0.607208517467, 0.597441877296,
0.587825531465, 0.578352913803, 0.569017984198, 0.559815170911,
0.550739320877, 0.541785656682, 0.532949739145, 0.524227434628,
0.515614886373, 0.507108489253, 0.498704867478, 0.490400854812,
0.482193476986, 0.47407993601, 0.466057596125, 0.458123971214,
0.450276713467, 0.442513603171, 0.434832539473, 0.427231532022,
0.419708693379, 0.41226223212, 0.404890446548, 0.397591718955,
0.390364510382, 0.383207355816, 0.376118859788, 0.369097692334,
0.362142585282, 0.355252328834, 0.348425768415, 0.341661801776,
0.334959376311, 0.328317486588, 0.321735172063, 0.31521151497,
0.308745638367, 0.302336704338, 0.29598391232, 0.289686497571,
0.283443729739, 0.27725491156, 0.271119377649, 0.265036493387,
0.259005653912, 0.253026283183, 0.247097833139, 0.241219782932,
0.235391638239, 0.229612930649, 0.223883217122, 0.218202079518,
0.212569124201, 0.206983981709, 0.201446306496, 0.195955776745,
0.190512094256, 0.185114984406, 0.179764196185, 0.174459502324,
0.169200699492, 0.1639876086, 0.158820075195, 0.153697969964,
0.148621189348, 0.143589656295, 0.138603321143, 0.133662162669,
0.128766189309, 0.123915440582, 0.119109988745, 0.114349940703,
0.10963544023, 0.104966670533, 0.100343857232, 0.0957672718266,
0.0912372357329, 0.0867541250127, 0.082318375932, 0.0779304915295,
0.0735910494266, 0.0693007111742, 0.065060233529, 0.0608704821745,
0.056732448584, 0.05264727098, 0.0486162607163, 0.0446409359769,
0.0407230655415, 0.0368647267386, 0.0330683839378, 0.0293369977411,
0.0256741818288, 0.0220844372634, 0.0185735200577, 0.0151490552854,
0.0118216532614, 0.00860719483079, 0.00553245272614, 0.00265435214565
};
/* Tabulated values for 2^24 times x[i]/x[i+1],
* used to accept for U*x[i+1]<=x[i] without any floating point operations */
static const unsigned long ktab[128] = {
0, 12590644, 14272653, 14988939,
15384584, 15635009, 15807561, 15933577,
16029594, 16105155, 16166147, 16216399,
16258508, 16294295, 16325078, 16351831,
16375291, 16396026, 16414479, 16431002,
16445880, 16459343, 16471578, 16482744,
16492970, 16502368, 16511031, 16519039,
16526459, 16533352, 16539769, 16545755,
16551348, 16556584, 16561493, 16566101,
16570433, 16574511, 16578353, 16581977,
16585398, 16588629, 16591685, 16594575,
16597311, 16599901, 16602354, 16604679,
16606881, 16608968, 16610945, 16612818,
16614592, 16616272, 16617861, 16619363,
16620782, 16622121, 16623383, 16624570,
16625685, 16626730, 16627708, 16628619,
16629465, 16630248, 16630969, 16631628,
16632228, 16632768, 16633248, 16633671,
16634034, 16634340, 16634586, 16634774,
16634903, 16634972, 16634980, 16634926,
16634810, 16634628, 16634381, 16634066,
16633680, 16633222, 16632688, 16632075,
16631380, 16630598, 16629726, 16628757,
16627686, 16626507, 16625212, 16623794,
16622243, 16620548, 16618698, 16616679,
16614476, 16612071, 16609444, 16606571,
16603425, 16599973, 16596178, 16591995,
16587369, 16582237, 16576520, 16570120,
16562917, 16554758, 16545450, 16534739,
16522287, 16507638, 16490152, 16468907,
16442518, 16408804, 16364095, 16301683,
16207738, 16047994, 15704248, 15472926
};
/* tabulated values of 2^{-24}*x[i] */
static const double wtab[128] = {
1.62318314817e-08, 2.16291505214e-08, 2.54246305087e-08, 2.84579525938e-08,
3.10340022482e-08, 3.33011726243e-08, 3.53439060345e-08, 3.72152672658e-08,
3.8950989572e-08, 4.05763964764e-08, 4.21101548915e-08, 4.35664624904e-08,
4.49563968336e-08, 4.62887864029e-08, 4.75707945735e-08, 4.88083237257e-08,
5.00063025384e-08, 5.11688950428e-08, 5.22996558616e-08, 5.34016475624e-08,
5.44775307871e-08, 5.55296344581e-08, 5.65600111659e-08, 5.75704813695e-08,
5.85626690412e-08, 5.95380306862e-08, 6.04978791776e-08, 6.14434034901e-08,
6.23756851626e-08, 6.32957121259e-08, 6.42043903937e-08, 6.51025540077e-08,
6.59909735447e-08, 6.68703634341e-08, 6.77413882848e-08, 6.8604668381e-08,
6.94607844804e-08, 7.03102820203e-08, 7.11536748229e-08, 7.1991448372e-08,
7.2824062723e-08, 7.36519550992e-08, 7.44755422158e-08, 7.52952223703e-08,
7.61113773308e-08, 7.69243740467e-08, 7.77345662086e-08, 7.85422956743e-08,
7.93478937793e-08, 8.01516825471e-08, 8.09539758128e-08, 8.17550802699e-08,
8.25552964535e-08, 8.33549196661e-08, 8.41542408569e-08, 8.49535474601e-08,
8.57531242006e-08, 8.65532538723e-08, 8.73542180955e-08, 8.8156298059e-08,
8.89597752521e-08, 8.97649321908e-08, 9.05720531451e-08, 9.138142487e-08,
9.21933373471e-08, 9.30080845407e-08, 9.38259651738e-08, 9.46472835298e-08,
9.54723502847e-08, 9.63014833769e-08, 9.71350089201e-08, 9.79732621669e-08,
9.88165885297e-08, 9.96653446693e-08, 1.00519899658e-07, 1.0138063623e-07,
1.02247952126e-07, 1.03122261554e-07, 1.04003996769e-07, 1.04893609795e-07,
1.05791574313e-07, 1.06698387725e-07, 1.07614573423e-07, 1.08540683296e-07,
1.09477300508e-07, 1.1042504257e-07, 1.11384564771e-07, 1.12356564007e-07,
1.13341783071e-07, 1.14341015475e-07, 1.15355110887e-07, 1.16384981291e-07,
1.17431607977e-07, 1.18496049514e-07, 1.19579450872e-07, 1.20683053909e-07,
1.21808209468e-07, 1.2295639141e-07, 1.24129212952e-07, 1.25328445797e-07,
1.26556042658e-07, 1.27814163916e-07, 1.29105209375e-07, 1.30431856341e-07,
1.31797105598e-07, 1.3320433736e-07, 1.34657379914e-07, 1.36160594606e-07,
1.37718982103e-07, 1.39338316679e-07, 1.41025317971e-07, 1.42787873535e-07,
1.44635331499e-07, 1.4657889173e-07, 1.48632138436e-07, 1.50811780719e-07,
1.53138707402e-07, 1.55639532047e-07, 1.58348931426e-07, 1.61313325908e-07,
1.64596952856e-07, 1.68292495203e-07, 1.72541128694e-07, 1.77574279496e-07,
1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
};
static unsigned int seed, iy, ir[98];
static uint32_t kn[128];static double wn[128],fn[128];
static unsigned long iz,jz,jsr=123456789;
static long hz;
/*
@defgroup _uniformdouble
@ingroup numerical Uniform linear congruential random number generator.
......@@ -153,6 +44,8 @@ static unsigned int seed, iy, ir[98];
#define mod 4294967296.0 /* is 2**32 */
#if 1
//#define randominit_SSE
#ifdef randominit_SSE
void randominit(unsigned seed_init)
{
int i;
......@@ -161,32 +54,39 @@ void randominit(unsigned seed_init)
if (seed_init == 0) {
srand((unsigned)time(NULL));
seed_array[0] = (unsigned int) rand();
seed_array[1] = (unsigned int) rand();
seed_array[2] = (unsigned int) rand();
seed_array[3] = (unsigned int) rand();
seed = (unsigned int) rand();
} else {
seed = seed_init;
seed_array[0] = seed_init;
seed_array[1] = log(seed_init);
seed_array[2] = pow(seed_init,3);
seed_array[3] = sqrt(seed_init);
}
if (seed % 2 == 0) seed += 1; /* seed and mod are relative prime */
for (i=0;i<4;i++)
if (seed_array[i] % 2 == 0) seed_array[i] += 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 */
}
for (i=1; i<4*98; i++) { /*4 times 98 to use in SSE implementations*/
seed_array[i%4] = a*seed_array[i%4]; /* mod 2**32 */
iy=1;
ir[i]= seed_array[i%4]; /* initialize the shuffle table */
}
for ( i = 0; i < 4*98; i++ ){
printf("ir[%d]: %d\n",i,ir[i]);
}
iy=1;
}
#endif
#if 0
#else
void randominit(unsigned seed_init)
{
int i;
// this need to be integrated with the existing rng, like taus: navid
msg("Initializing random number generator, seed %x\n",seed_init);
seed_init = 62110;
if (seed_init == 0) {
srand((unsigned)time(NULL));
......@@ -201,25 +101,120 @@ void randominit(unsigned seed_init)
seed = a*seed; /* mod 2**32 */
ir[i]= seed; /* initialize the shuffle table */
}
iy=1;
iy=1;
}
#endif
#endif
//Procedure to create tables for normal distribution kn,wn and fn
#define SHR3 (jz=jsr, jsr^=(jsr<<13),jsr^=(jsr>>17),jsr^=(jsr<<5),jz+jsr)
#define UNI (0.5+(signed) SHR3 * 0.2328306e-9)
#define NOR (hz=SHR3, iz=hz&127, (abs(hz)<kn[iz])? hz*wn[iz] : nfix())
double nfix()
{
const double r = 3.442620; static double x, y;
for (;;)
{
x=hz*wn[iz];
if (iz==0)
{
do
{
x = - 0.2904764 * log (UNI);
y = - log (UNI);
}
while (y+y < x*x);
return (hz>0)? r+x : -r-x;
}
if (fn[iz]+UNI*(fn[iz-1]-fn[iz])<exp(-0.5*x*x))
return x;
hz = SHR3;
iz = hz&127;
if (abs(hz) < kn[iz])
return (hz*wn[iz]);
}
}
void setup_nor()
{
double dn = 3.442619855899;
int i;
const double m1 = 2147483648.0;
double q;
double tn = 3.442619855899;
const double vn = 9.91256303526217E-03;
q = vn/exp(-0.5*dn*dn);
kn[0] = (int)((dn/q)*m1);
kn[1] = 0;
wn[0] = (double) ( q / m1 );
wn[127] = (double) ( dn / m1 );
fn[0] = 1.0;
fn[127] = (double) ( exp ( - 0.5 * dn * dn ) );
for ( i = 126; 1 <= i; i-- )
{
dn = sqrt (-2.0 * log ( vn/dn + exp(-0.5*dn*dn)));
kn[i+1] = (int) ((dn / tn)*m1);
tn = dn;
fn[i] = (double) (exp (-0.5*dn*dn));
wn[i] = (double) (dn / m1);
}
for ( i = 0; i <= 127; i++ ){
printf("i %d: kn %d, fn %e, wn %e\n",i,kn[i],fn[i],wn[i]);
}
return;
}
/*!\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 );
}
/*void uniformrandomSSE(__m128d *d1,__m128d *d2)
{
int i,j;
__m128d u_rand128;
__m128i j128;
int j_array[]={(int) (1 + 9.103678167e-08*iy_array[0]),(int) (1 + 9.103678167e-08*iy_array[1]),(int) (1 + 9.103678167e-08*iy_array[2]),(int) (1 + 9.103678167e-08*iy_array[3])};
//j128=_mm_setr_epi32(1 + 391.0*iy_array[0]/mod,1 + 391.0*iy_array[1]/mod,1 + 391.0*iy_array[2]/mod,1 + 391.0*iy_array[3]/mod);
//j128=_mm_setr_epi32(391,391,391,391);
//j128=_mm_mul_epu32(j128,_mm_setr_epi32(iy_array[0],iy_array[1],iy_array[2],iy_array[3]));
//j128=_mm_slli_epi32(j128,32);
//j=1 + 97.0*iy/mod;
iy_array[0]=ir[j_array[0]];
iy_array[1]=ir[j_array[1]];
iy_array[2]=ir[j_array[2]];
iy_array[3]=ir[j_array[3]];
seed_array[0] = a*seed_array[0];
seed_array[1] = a*seed_array[1];
seed_array[2] = a*seed_array[2];
seed_array[3] = a*seed_array[3];
ir[j_array[0]] = seed_array[0];
ir[j_array[1]] = seed_array[1];
ir[j_array[2]] = seed_array[2];
ir[j_array[3]] = seed_array[3];
*d1=_mm_setr_pd (2*((double) iy_array[0]/mod)-1, 2*((double) iy_array[1])-1);
*d2=_mm_setr_pd (2*((double) iy_array[2]/mod)-1, 2*((double) iy_array[3])-1);
//return ((double) iy/mod );
return;
}*/
/*!\brief Ziggurat random number generator based on rejection sampling. It returns a pseudorandom normally distributed double number between -4.5 and 4.5*/
double ziggurat()
{
return NOR;
}
/*
@defgroup _gaussdouble Gaussian random number generator based on modified Box-Muller transformation.
......@@ -227,58 +222,103 @@ double uniformrandom(void)
*/
/*!\brief Gaussian random number generator based on modified Box-Muller transformation.Returns a double-precision floating-point number. */
#define Ziggurat_SSE
#ifdef Ziggurat_SSE
double gaussdouble(double mean, double variance)
#define random_SSE
#ifdef random_SSE
double gaussdouble(double mean, double variance)//It is necessary to improve the function. However if we enable SSE the gain in time it is not too much.
{
static int iset=0;
static double gset;
double fac,r,v1,v2;
double fac,rn,r[2],v1[2],v2[2],ones[]={1.0,1.0};
static double max=-1000000;
static double min=1000000;
__m128d v1_128, v2_128, r128, ones128, compge128;
ones128 = _mm_load_pd(ones);
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);
//v2 = 2.0*uniformrandom()-1.0;
v1_128 = _mm_set_pd(2*UNI-1,2*UNI-1);
v2_128 = _mm_set_pd(2*UNI-1,2*UNI-1);
//r = v1*v1+v2*v2;
v1_128= _mm_mul_pd(v1_128,v1_128);
v2_128= _mm_mul_pd(v2_128,v2_128);
r128= _mm_add_pd(v1_128,v2_128);
//compge128=_mm_cmpge_pd(r128,ones128);
_mm_storeu_pd(&r[0],r128);
//printf("Inside do: r[0] %e, r[1] %e\n",r[0],r[1]);
} while (r[0] >= 1.0 && r[1]>=1.0);
//printf("outside do: r[0] %e, r[1] %e\n",r[0],r[1]);
if (r[0]<r[1]){
fac = sqrt(-2.0*log(r[0])/r[0]);
gset= v1[0]*fac;
iset=1;
return(sqrt(variance)*v1[1]*fac + mean);
}
else{
fac = sqrt(-2.0*log(r[1])/r[1]);
gset= v2[0]*fac;
iset=1;
return(sqrt(variance)*v2[1]*fac + mean);
}
} else {
iset=0;
//printf("normal random number %e, max %e, min %e\n",sqrt(variance)*gset + mean, max,min);
if (max<sqrt(variance)*gset + mean)
max=sqrt(variance)*gset + mean;
if (min>sqrt(variance)*gset + mean)
min=sqrt(variance)*gset + mean;
return(sqrt(variance)*gset + mean);
}
}
#else
double gaussdouble(double mean, double variance)
{
/*static int first_run;
static double sum;
static int count;
if (!first_run)
{
first_run=1;
sum=0;
count=0;
} */
static int iset=0;
static double gset;
double fac,r,v1,v2;
static double max=-1000000;
static double min=1000000;
if (iset == 0) {
do {
/*count++;
clock_t start=clock();*/
v1 = 2.0*uniformrandom()-1.0;
/*clock_t stop=clock();
printf("UE_freq_channel time is %f s, AVERAGE time is %f s, count %d, sum %e\n",(float) (stop-start)/CLOCKS_PER_SEC,(float) (sum+stop-start)/(count*CLOCKS_PER_SEC),count,sum+stop-start);
sum=(sum+stop-start);*/
v2 = 2.0*uniformrandom()-1.0;
r = v1*v1+v2*v2;
//printf("Inside do: r %e\n",r);
} while (r >= 1.0);
//printf("outside do: r %e\n",r);
fac = sqrt(-2.0*log(r)/r);
gset= v1*fac;
iset=1;
return(sqrt(variance)*v2*fac + mean);
} else {
iset=0;
//printf("normal random number %e, max %e, min %e\n",sqrt(variance)*gset + mean, max,min);
if (max<sqrt(variance)*gset + mean)
max=sqrt(variance)*gset + mean;
if (min>sqrt(variance)*gset + mean)
min=sqrt(variance)*gset + mean;
return(sqrt(variance)*gset + mean);
}
}
#endif
#ifdef MAIN
main(int argc,char **argv)
{
......
......@@ -139,7 +139,7 @@ eNBs =
////////// MME parameters:
mme_ip_address = ( { ipv4 = "172.24.11.4";
mme_ip_address = ( { ipv4 = "192.168.12.148";
ipv6 = "192:168:30::17";
active = "yes";
preference = "ipv4";
......
......@@ -990,10 +990,12 @@ void init_seed(uint8_t set_seed)
if(set_seed) {
randominit (oai_emulation.info.seed);
setup_nor();//Setup for the normal probability distribution function
set_taus_seed (oai_emulation.info.seed);
} else {
randominit (0);
setup_nor();//Setup for the normal probability distribution function
set_taus_seed (0);
}
}
......@@ -1042,7 +1044,6 @@ extern int subframe_eNB_mask,subframe_UE_mask;
int eNB_trx_read(openair0_device *device, openair0_timestamp *ptimestamp, void **buff, int nsamps, int cc)
{
static int count=0;
int eNB,UE,CC;
if (count==500)
{
//Use ./oaisim -q option to enable the oai performance profiler.
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment