@@ -116,13 +116,13 @@ public double ChiSquared(int df) {
116
116
if ( df <= 0 ) {
117
117
throw new ArgumentException ( "df (" + df + ") is not greater than 0" ) ;
118
118
}
119
- return this . Gamma ( df / 2.0 , 2 ) ;
119
+ return this . Gamma ( df * 0.5 , 2 ) ;
120
120
}
121
121
122
122
/// <summary>Not documented yet.</summary>
123
123
/// <returns>A 64-bit floating-point number.</returns>
124
124
public double Exponential ( ) {
125
- return - Math . Log ( this . NonZeroUniform ( ) ) ;
125
+ return - Math . Log ( 1.0 - this . Uniform ( ) ) ;
126
126
}
127
127
128
128
/// <summary>Not documented yet.</summary>
@@ -140,9 +140,8 @@ public double Gamma(double a, double b) {
140
140
/// <param name='a'>Another 64-bit floating-point number.</param>
141
141
/// <returns>A 64-bit floating-point number.</returns>
142
142
public double Gamma ( double a ) {
143
- if ( a <= - 1 ) {
144
- throw new ArgumentException ( "a (" + a + ") is not greater than " +
145
- ( - 1 ) ) ;
143
+ if ( a <= 0 ) {
144
+ throw new ArgumentException ( "a (" + a + ") is not greater than 0" ) ;
146
145
}
147
146
double v , x , u , x2 , d , c ;
148
147
d = ( a < 1 ? 1 + a : a ) - ( 1.0 / 3.0 ) ;
@@ -152,9 +151,9 @@ public double Gamma(double a) {
152
151
x = this . Normal ( ) ;
153
152
v = Math . Pow ( ( c * x ) + 1 , 3 ) ;
154
153
} while ( v <= 0 ) ;
155
- u = this . Uniform ( ) ;
156
- x2 = Math . Pow ( x , 2 ) ;
157
- } while ( u >= 1 - ( 0.0331 * x2 * x2 ) &&
154
+ u = 1.0 - this . Uniform ( ) ;
155
+ x2 = x * x ;
156
+ } while ( u >= 1 - ( 0.0331 * x2 * x2 ) &&
158
157
Math . Log ( u ) >= ( 0.5 * x2 ) + ( d * ( 1 - v + Math . Log ( v ) ) ) ) ;
159
158
if ( a < 1 ) {
160
159
return d * v * Math . Exp ( this . Exponential ( ) / - a ) ;
@@ -226,7 +225,7 @@ public int Hypergeometric(int trials, int ones, int count) {
226
225
/// <param name='sd'>Standard deviation.</param>
227
226
/// <returns>A 64-bit floating-point number.</returns>
228
227
public double LogNormal ( double mean , double sd ) {
229
- return Math . Exp ( ( this . Normal ( ) * sd ) + mean ) ;
228
+ return Math . Exp ( this . Normal ( mean , sd ) ) ;
230
229
}
231
230
232
231
/// <summary>Conceptually, generates either 1 or 0 until the given
@@ -313,7 +312,7 @@ public double Normal() {
313
312
return this . valueLastNormal ;
314
313
}
315
314
}
316
- double x = this . NonZeroUniform ( ) ;
315
+ double x = 1.0 - this . Uniform ( ) ;
317
316
double y = this . Uniform ( ) ;
318
317
double s = Math . Sqrt ( - 2 * Math . Log ( x ) ) ;
319
318
double t = 2 * Math . PI * y ;
@@ -345,13 +344,15 @@ public int Poisson(double mean) {
345
344
") is less than 0" ) ;
346
345
}
347
346
double l = Math . Exp ( - mean ) ;
348
- var k = 0 ;
349
- double p = 0 ;
350
- do {
351
- ++ k ;
347
+ var count = - 1 ;
348
+ double p = 1. 0;
349
+ while ( true ) {
350
+ ++ count ;
352
351
p *= this . Uniform ( ) ;
353
- } while ( p > l ) ;
354
- return k - 1 ;
352
+ if ( p <= l ) {
353
+ return count ;
354
+ }
355
+ }
355
356
}
356
357
357
358
/// <summary>Not documented yet.</summary>
@@ -381,32 +382,14 @@ public double Uniform(double max) {
381
382
/// number from 0 and up, but less than 1.</summary>
382
383
/// <returns>A 64-bit floating-point number.</returns>
383
384
public double Uniform ( ) {
384
- var b = new byte [ 7 ] ;
385
- this . valueIrg . GetBytes ( b , 0 , 7 ) ;
386
- var lb = ( long ) b [ 0 ] & 0xffL ;
387
- lb |= ( ( long ) b [ 1 ] & 0xffL ) << 8 ;
388
- lb |= ( ( long ) b [ 2 ] & 0xffL ) << 16 ;
389
- lb |= ( ( long ) b [ 3 ] & 0xffL ) << 24 ;
390
- lb |= ( ( long ) b [ 4 ] & 0xffL ) << 32 ;
391
- lb |= ( ( long ) b [ 5 ] & 0xffL ) << 40 ;
392
- lb |= ( ( long ) b [ 6 ] & 0xfL ) << 48 ;
393
- lb |= 0x3ffL << 52 ;
394
- return BitConverter . Int64BitsToDouble ( lb ) - 1.0 ;
385
+ return this . UniformLong ( 9007199254740992L ) / 9007199254740992.0 ;
395
386
}
396
387
397
388
/// <summary>Returns a uniformly-distributed 32-bit floating-point
398
389
/// number from 0 and up, but less than 1.</summary>
399
- /// <returns>A 64 -bit floating-point number.</returns>
390
+ /// <returns>A 32 -bit floating-point number.</returns>
400
391
public double UniformSingle ( ) {
401
- var b = new byte [ 3 ] ;
402
- this . valueIrg . GetBytes ( b , 0 , 3 ) ;
403
- var lb = ( int ) b [ 0 ] & 0xff ;
404
- lb |= ( ( int ) b [ 1 ] & 0xff ) << 8 ;
405
- lb |= ( ( int ) b [ 2 ] & 0x7f ) << 16 ;
406
- lb |= 0x7f << 23 ;
407
- return BitConverter . ToSingle (
408
- BitConverter . GetBytes ( lb ) ,
409
- 0 ) - 1.0 ;
392
+ return this . UniformInt ( 16777216 ) / 16777216.0f ;
410
393
}
411
394
412
395
/// <summary>Generates a random 32-bit signed integer within a given
@@ -490,7 +473,7 @@ public int UniformInt(int maxExclusive) {
490
473
throw new ArgumentException ( "maxExclusive (" + maxExclusive +
491
474
") is less than 0" ) ;
492
475
}
493
- if ( maxExclusive <= 0 ) {
476
+ if ( maxExclusive <= 1 ) {
494
477
return 0 ;
495
478
}
496
479
var b = new byte [ 4 ] ;
@@ -520,8 +503,7 @@ public int UniformInt(int maxExclusive) {
520
503
return ib ;
521
504
}
522
505
int maxexc ;
523
- maxexc = ( maxExclusive <= 100 ) ? 2147483600 :
524
- ( ( Int32 . MaxValue / maxExclusive ) * maxExclusive ) ;
506
+ maxexc = ( Int32 . MaxValue / maxExclusive ) * maxExclusive ;
525
507
while ( true ) {
526
508
this . valueIrg . GetBytes ( b , 0 , 4 ) ;
527
509
ib = b [ 0 ] & 0xff ;
@@ -569,22 +551,5 @@ public long UniformLong(long maxExclusive) {
569
551
}
570
552
}
571
553
}
572
-
573
- private double NonZeroUniform ( ) {
574
- var b = new byte [ 7 ] ;
575
- long lb = 0 ;
576
- do {
577
- this . valueIrg . GetBytes ( b , 0 , 7 ) ;
578
- lb = b [ 0 ] & 0xffL ;
579
- lb |= ( b [ 1 ] & 0xffL ) << 8 ;
580
- lb |= ( b [ 2 ] & 0xffL ) << 16 ;
581
- lb |= ( b [ 3 ] & 0xffL ) << 24 ;
582
- lb |= ( b [ 4 ] & 0xffL ) << 32 ;
583
- lb |= ( b [ 5 ] & 0xffL ) << 40 ;
584
- lb |= ( b [ 6 ] & 0xfL ) << 48 ;
585
- } while ( lb == 0 ) ;
586
- lb |= 0x3ffL << 52 ;
587
- return BitConverter . Int64BitsToDouble ( lb ) - 1.0 ;
588
- }
589
554
}
590
555
}
0 commit comments