[3839] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c)
|
---|
| 3 | 2007, Sergey Bochkanov (ALGLIB project).
|
---|
| 4 | 1988, Pierre L'Ecuyer
|
---|
| 5 |
|
---|
| 6 | >>> SOURCE LICENSE >>>
|
---|
| 7 | This program is free software; you can redistribute it and/or modify
|
---|
| 8 | it under the terms of the GNU General Public License as published by
|
---|
| 9 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 10 | License, or (at your option) any later version.
|
---|
| 11 |
|
---|
| 12 | This program is distributed in the hope that it will be useful,
|
---|
| 13 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 14 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 15 | GNU General Public License for more details.
|
---|
| 16 |
|
---|
| 17 | A copy of the GNU General Public License is available at
|
---|
| 18 | http://www.fsf.org/licensing/licenses
|
---|
| 19 |
|
---|
| 20 | >>> END OF LICENSE >>>
|
---|
| 21 | *************************************************************************/
|
---|
| 22 |
|
---|
| 23 | using System;
|
---|
| 24 |
|
---|
| 25 | namespace alglib
|
---|
| 26 | {
|
---|
| 27 | public class hqrnd
|
---|
| 28 | {
|
---|
| 29 | /*************************************************************************
|
---|
| 30 | Portable high quality random number generator state.
|
---|
| 31 | Initialized with HQRNDRandomize() or HQRNDSeed().
|
---|
| 32 |
|
---|
| 33 | Fields:
|
---|
| 34 | S1, S2 - seed values
|
---|
| 35 | V - precomputed value
|
---|
| 36 | MagicV - 'magic' value used to determine whether State structure
|
---|
| 37 | was correctly initialized.
|
---|
| 38 | *************************************************************************/
|
---|
| 39 | public struct hqrndstate
|
---|
| 40 | {
|
---|
| 41 | public int s1;
|
---|
| 42 | public int s2;
|
---|
| 43 | public double v;
|
---|
| 44 | public int magicv;
|
---|
| 45 | };
|
---|
| 46 |
|
---|
| 47 |
|
---|
| 48 |
|
---|
| 49 |
|
---|
| 50 | public const int hqrndmax = 2147483563;
|
---|
| 51 | public const int hqrndm1 = 2147483563;
|
---|
| 52 | public const int hqrndm2 = 2147483399;
|
---|
| 53 | public const int hqrndmagic = 1634357784;
|
---|
| 54 |
|
---|
| 55 |
|
---|
| 56 | /*************************************************************************
|
---|
| 57 | HQRNDState initialization with random values which come from standard
|
---|
| 58 | RNG.
|
---|
| 59 |
|
---|
| 60 | -- ALGLIB --
|
---|
| 61 | Copyright 02.12.2009 by Bochkanov Sergey
|
---|
| 62 | *************************************************************************/
|
---|
| 63 | public static void hqrndrandomize(ref hqrndstate state)
|
---|
| 64 | {
|
---|
| 65 | hqrndseed(AP.Math.RandomInteger(hqrndm1), AP.Math.RandomInteger(hqrndm2), ref state);
|
---|
| 66 | }
|
---|
| 67 |
|
---|
| 68 |
|
---|
| 69 | /*************************************************************************
|
---|
| 70 | HQRNDState initialization with seed values
|
---|
| 71 |
|
---|
| 72 | -- ALGLIB --
|
---|
| 73 | Copyright 02.12.2009 by Bochkanov Sergey
|
---|
| 74 | *************************************************************************/
|
---|
| 75 | public static void hqrndseed(int s1,
|
---|
| 76 | int s2,
|
---|
| 77 | ref hqrndstate state)
|
---|
| 78 | {
|
---|
| 79 | state.s1 = s1%(hqrndm1-1)+1;
|
---|
| 80 | state.s2 = s2%(hqrndm2-1)+1;
|
---|
| 81 | state.v = (double)(1)/(double)(hqrndmax);
|
---|
| 82 | state.magicv = hqrndmagic;
|
---|
| 83 | }
|
---|
| 84 |
|
---|
| 85 |
|
---|
| 86 | /*************************************************************************
|
---|
| 87 | This function generates random real number in (0,1),
|
---|
| 88 | not including interval boundaries
|
---|
| 89 |
|
---|
| 90 | State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
---|
| 91 |
|
---|
| 92 | -- ALGLIB --
|
---|
| 93 | Copyright 02.12.2009 by Bochkanov Sergey
|
---|
| 94 | *************************************************************************/
|
---|
| 95 | public static double hqrnduniformr(ref hqrndstate state)
|
---|
| 96 | {
|
---|
| 97 | double result = 0;
|
---|
| 98 |
|
---|
| 99 | result = state.v*hqrndintegerbase(ref state);
|
---|
| 100 | return result;
|
---|
| 101 | }
|
---|
| 102 |
|
---|
| 103 |
|
---|
| 104 | /*************************************************************************
|
---|
| 105 | This function generates random integer number in [0, N)
|
---|
| 106 |
|
---|
| 107 | 1. N must be less than HQRNDMax-1.
|
---|
| 108 | 2. State structure must be initialized with HQRNDRandomize() or HQRNDSeed()
|
---|
| 109 |
|
---|
| 110 | -- ALGLIB --
|
---|
| 111 | Copyright 02.12.2009 by Bochkanov Sergey
|
---|
| 112 | *************************************************************************/
|
---|
| 113 | public static int hqrnduniformi(int n,
|
---|
| 114 | ref hqrndstate state)
|
---|
| 115 | {
|
---|
| 116 | int result = 0;
|
---|
| 117 | int mx = 0;
|
---|
| 118 |
|
---|
| 119 |
|
---|
| 120 | //
|
---|
| 121 | // Correct handling of N's close to RNDBaseMax
|
---|
| 122 | // (avoiding skewed distributions for RNDBaseMax<>K*N)
|
---|
| 123 | //
|
---|
| 124 | System.Diagnostics.Debug.Assert(n>0, "HQRNDUniformI: N<=0!");
|
---|
| 125 | System.Diagnostics.Debug.Assert(n<hqrndmax-1, "HQRNDUniformI: N>=RNDBaseMax-1!");
|
---|
| 126 | mx = hqrndmax-1-(hqrndmax-1)%n;
|
---|
| 127 | do
|
---|
| 128 | {
|
---|
| 129 | result = hqrndintegerbase(ref state)-1;
|
---|
| 130 | }
|
---|
| 131 | while( result>=mx );
|
---|
| 132 | result = result%n;
|
---|
| 133 | return result;
|
---|
| 134 | }
|
---|
| 135 |
|
---|
| 136 |
|
---|
| 137 | /*************************************************************************
|
---|
| 138 | Random number generator: normal numbers
|
---|
| 139 |
|
---|
| 140 | This function generates one random number from normal distribution.
|
---|
| 141 | Its performance is equal to that of HQRNDNormal2()
|
---|
| 142 |
|
---|
| 143 | State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
---|
| 144 |
|
---|
| 145 | -- ALGLIB --
|
---|
| 146 | Copyright 02.12.2009 by Bochkanov Sergey
|
---|
| 147 | *************************************************************************/
|
---|
| 148 | public static double hqrndnormal(ref hqrndstate state)
|
---|
| 149 | {
|
---|
| 150 | double result = 0;
|
---|
| 151 | double v1 = 0;
|
---|
| 152 | double v2 = 0;
|
---|
| 153 |
|
---|
| 154 | hqrndnormal2(ref state, ref v1, ref v2);
|
---|
| 155 | result = v1;
|
---|
| 156 | return result;
|
---|
| 157 | }
|
---|
| 158 |
|
---|
| 159 |
|
---|
| 160 | /*************************************************************************
|
---|
| 161 | Random number generator: random X and Y such that X^2+Y^2=1
|
---|
| 162 |
|
---|
| 163 | State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
---|
| 164 |
|
---|
| 165 | -- ALGLIB --
|
---|
| 166 | Copyright 02.12.2009 by Bochkanov Sergey
|
---|
| 167 | *************************************************************************/
|
---|
| 168 | public static void hqrndunit2(ref hqrndstate state,
|
---|
| 169 | ref double x,
|
---|
| 170 | ref double y)
|
---|
| 171 | {
|
---|
| 172 | double v = 0;
|
---|
| 173 | double mx = 0;
|
---|
| 174 | double mn = 0;
|
---|
| 175 |
|
---|
| 176 | do
|
---|
| 177 | {
|
---|
| 178 | hqrndnormal2(ref state, ref x, ref y);
|
---|
| 179 | }
|
---|
| 180 | while( ! ((double)(x)!=(double)(0) | (double)(y)!=(double)(0)) );
|
---|
| 181 | mx = Math.Max(Math.Abs(x), Math.Abs(y));
|
---|
| 182 | mn = Math.Min(Math.Abs(x), Math.Abs(y));
|
---|
| 183 | v = mx*Math.Sqrt(1+AP.Math.Sqr(mn/mx));
|
---|
| 184 | x = x/v;
|
---|
| 185 | y = y/v;
|
---|
| 186 | }
|
---|
| 187 |
|
---|
| 188 |
|
---|
| 189 | /*************************************************************************
|
---|
| 190 | Random number generator: normal numbers
|
---|
| 191 |
|
---|
| 192 | This function generates two independent random numbers from normal
|
---|
| 193 | distribution. Its performance is equal to that of HQRNDNormal()
|
---|
| 194 |
|
---|
| 195 | State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
---|
| 196 |
|
---|
| 197 | -- ALGLIB --
|
---|
| 198 | Copyright 02.12.2009 by Bochkanov Sergey
|
---|
| 199 | *************************************************************************/
|
---|
| 200 | public static void hqrndnormal2(ref hqrndstate state,
|
---|
| 201 | ref double x1,
|
---|
| 202 | ref double x2)
|
---|
| 203 | {
|
---|
| 204 | double u = 0;
|
---|
| 205 | double v = 0;
|
---|
| 206 | double s = 0;
|
---|
| 207 |
|
---|
| 208 | while( true )
|
---|
| 209 | {
|
---|
| 210 | u = 2*hqrnduniformr(ref state)-1;
|
---|
| 211 | v = 2*hqrnduniformr(ref state)-1;
|
---|
| 212 | s = AP.Math.Sqr(u)+AP.Math.Sqr(v);
|
---|
| 213 | if( (double)(s)>(double)(0) & (double)(s)<(double)(1) )
|
---|
| 214 | {
|
---|
| 215 |
|
---|
| 216 | //
|
---|
| 217 | // two Sqrt's instead of one to
|
---|
| 218 | // avoid overflow when S is too small
|
---|
| 219 | //
|
---|
| 220 | s = Math.Sqrt(-(2*Math.Log(s)))/Math.Sqrt(s);
|
---|
| 221 | x1 = u*s;
|
---|
| 222 | x2 = v*s;
|
---|
| 223 | return;
|
---|
| 224 | }
|
---|
| 225 | }
|
---|
| 226 | }
|
---|
| 227 |
|
---|
| 228 |
|
---|
| 229 | /*************************************************************************
|
---|
| 230 | Random number generator: exponential distribution
|
---|
| 231 |
|
---|
| 232 | State structure must be initialized with HQRNDRandomize() or HQRNDSeed().
|
---|
| 233 |
|
---|
| 234 | -- ALGLIB --
|
---|
| 235 | Copyright 11.08.2007 by Bochkanov Sergey
|
---|
| 236 | *************************************************************************/
|
---|
| 237 | public static double hqrndexponential(double lambda,
|
---|
| 238 | ref hqrndstate state)
|
---|
| 239 | {
|
---|
| 240 | double result = 0;
|
---|
| 241 |
|
---|
| 242 | System.Diagnostics.Debug.Assert((double)(lambda)>(double)(0), "HQRNDExponential: Lambda<=0!");
|
---|
| 243 | result = -(Math.Log(hqrnduniformr(ref state))/lambda);
|
---|
| 244 | return result;
|
---|
| 245 | }
|
---|
| 246 |
|
---|
| 247 |
|
---|
| 248 | /*************************************************************************
|
---|
| 249 |
|
---|
| 250 | L'Ecuyer, Efficient and portable combined random number generators
|
---|
| 251 | *************************************************************************/
|
---|
| 252 | private static int hqrndintegerbase(ref hqrndstate state)
|
---|
| 253 | {
|
---|
| 254 | int result = 0;
|
---|
| 255 | int k = 0;
|
---|
| 256 |
|
---|
| 257 | System.Diagnostics.Debug.Assert(state.magicv==hqrndmagic, "HQRNDIntegerBase: State is not correctly initialized!");
|
---|
| 258 | k = state.s1/53668;
|
---|
| 259 | state.s1 = 40014*(state.s1-k*53668)-k*12211;
|
---|
| 260 | if( state.s1<0 )
|
---|
| 261 | {
|
---|
| 262 | state.s1 = state.s1+2147483563;
|
---|
| 263 | }
|
---|
| 264 | k = state.s2/52774;
|
---|
| 265 | state.s2 = 40692*(state.s2-k*52774)-k*3791;
|
---|
| 266 | if( state.s2<0 )
|
---|
| 267 | {
|
---|
| 268 | state.s2 = state.s2+2147483399;
|
---|
| 269 | }
|
---|
| 270 |
|
---|
| 271 | //
|
---|
| 272 | // Result
|
---|
| 273 | //
|
---|
| 274 | result = state.s1-state.s2;
|
---|
| 275 | if( result<1 )
|
---|
| 276 | {
|
---|
| 277 | result = result+2147483562;
|
---|
| 278 | }
|
---|
| 279 | return result;
|
---|
| 280 | }
|
---|
| 281 | }
|
---|
| 282 | }
|
---|