1 //==============================================================================
2 //
3 // This file is part of GPSTk, the GPS Toolkit.
4 //
5 // The GPSTk is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published
7 // by the Free Software Foundation; either version 3.0 of the License, or
8 // any later version.
9 //
10 // The GPSTk is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with GPSTk; if not, write to the Free Software Foundation,
17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
18 //
19 // This software was developed by Applied Research Laboratories at the
20 // University of Texas at Austin.
21 // Copyright 2004-2020, The Board of Regents of The University of Texas System
22 //
23 //==============================================================================
24
25 //==============================================================================
26 //
27 // This software was developed by Applied Research Laboratories at the
28 // University of Texas at Austin, under contract to an agency or agencies
29 // within the U.S. Department of Defense. The U.S. Government retains all
30 // rights to use, duplicate, distribute, disclose, or release this software.
31 //
32 // Pursuant to DoD Directive 523024
33 //
34 // DISTRIBUTION STATEMENT A: This software has been approved for public
35 // release, distribution is unlimited.
36 //
37 //==============================================================================
38
39 /// @file random.cpp
40 /// Simple random number generator.
41
42 // -----------------------------------------------------------------------------
43 #include <iostream>
44 #include <iomanip>
45 #include <cmath>
46 #include <ctime>
47
48 #include "random.hpp"
49
50 using namespace std;
51
52 // -----------------------------------------------------------------------------
53 // Generate random numbers uniformly distributed from 0.0 to 1.0. Mbig and
54 // Mseed are large but arbitrary, but Mbig > Mseed. The 55 is not arbitrary.
Rand(long seed)55 double Rand(long seed)
56 {
57 #define Mbig 1000000000.
58 #define Mseed 161803398.
59 #define imod(x,y) ((x)-((x)/(y))*(y))
60 static short iff=0,inext,inextp;
61 static double Ma[55];
62 double mj,mk;
63 short i,ii,k;
64 if(!iff) {
65 if(seed < 0) seed=-seed;
66 mj = Mseed-seed;
67 mj = fmod(mj,Mbig);
68 Ma[54] = mj;
69 mk = 1.0;
70 for(i=0; i<55; i++) {
71 ii = imod(21*(i+1),55);
72 Ma[ii] = mk;
73 mk = mj - mk;
74 if(mk < 0.0) mk += Mbig;
75 mj = Ma[ii];
76 }
77 for(k=0; k<4; k++) {
78 for(i=0; i<55; i++) {
79 Ma[i] -= Ma[imod(i+30,55)];
80 if(Ma[i] < 0.0) Ma[i] += Mbig;
81 }
82 }
83 inext = -1;
84 inextp = 30;
85 iff=1;
86 }
87 inext++; if(inext == 55) inext=0;
88 inextp++; if(inextp == 55) inextp=0;
89 mj = Ma[inext]-Ma[inextp];
90 if(mj < 0.0) mj += Mbig;
91 Ma[inext] = mj;
92 return mj/Mbig;
93 #undef Mbig
94 #undef Mseed
95 #undef imod
96 }
97
98 // -----------------------------------------------------------------------------
99 // Generate normally distributed random numbers, zero mean and
100 // sqrt of variance sigma. Uses Box-Muller and Rand() above.
RandNorm(double sigma)101 double RandNorm(double sigma)
102 {
103 #ifdef RAND_NORM_SAVE
104 static short iset=0;
105 static double saved;
106 double r, v1, v2, fact;
107 if(!iset) {
108 do {
109 v1 = 2.0*Rand(1)-1.0;
110 v2 = 2.0*Rand(1)-1.0;
111 r = v1*v1 + v2*v2;
112 } while( r >= 1.0 || r == 0.0);
113 fact = sigma*sqrt(-2.*log(r)/r);
114 saved = v1*fact;
115 iset = 1;
116 return v2*fact;
117 }
118 iset = 0;
119 return saved;
120 #else
121 double r, v1, v2, fact;
122 do {
123 v1 = 2.0*Rand(1)-1.0;
124 v2 = 2.0*Rand(1)-1.0;
125 r = v1*v1 + v2*v2;
126 } while( r >= 1.0 || r == 0.0);
127 fact = sigma*sqrt(-2.*log(r)/r);
128 return v2*fact;
129 #endif
130 }
131
132 // -----------------------------------------------------------------------------
133 // Return random integers between low and hi. If you want a different seed,
134 // call Rand(seed) before you call this.
ARand(int low,int hi)135 int ARand(int low, int hi)
136 {
137 double r=Rand(),d=(double)(hi-low);
138 if(d < 0.0) d = -d;
139 d = r*d;
140 int i=(int)(d+0.5) + low;
141 return i;
142 }
143
144 // -----------------------------------------------------------------------------
145 // Return random doubles between low and hi. If you want a different seed,
146 // call Rand(seed) before you call this.
ARand(double low,double hi)147 double ARand(double low, double hi)
148 {
149 double r=Rand(),d=(hi-low);
150 if(d < 0.0) d = -d;
151 d = r*d;
152 return (low+d);
153 }
154
155 // -----------------------------------------------------------------------------
156 // Generate a random walk sequence, given sqrt variance sigma, time step dt
157 // and previous point xlast.
RandomWalk(double dt,double sigma,double xlast)158 double RandomWalk(double dt, double sigma, double xlast)
159 {
160 return xlast+RandNorm(sigma)*dt;
161 }
162
163 // -----------------------------------------------------------------------------
164 // Generate an exponentially correlated random sequence, given time step dt,
165 // sqrt variance sigma, time constant T and previous point xlast.
RandExpCor(double dt,double sigma,double T,double xlast)166 double RandExpCor(double dt, double sigma, double T, double xlast)
167 {
168 return exp(-dt/T)*xlast+RandNorm(sigma);
169 }
170
171 // -----------------------------------------------------------------------------
172 // integer mod function. assume arguments positive.
173 //static int imod(int x, int y)
174 //{
175 //if(x == 0 || y == 0) return 0;
176 //return (x-(x/y)*y);
177 //}
178
179 // -----------------------------------------------------------------------------
180 // -----------------------------------------------------------------------------
181