我实现一个随机数生成器。

emacsnw

我做机器学习方面的研究，因此非常需要高质量的（伪）随机数生成器，下面这个是我参考"Numerical Recipes in C (NRIC)"和一些其它资料实现的一个C++ Wrapper for random number generator. 当然网上有很多类似的代码，不过我觉得这个应付常见的需求足够了，至少基本满足自己的需求。
它特点就是简单高效，但是随机数质量高，基于48位伪随机数生成引擎，实现了比较常见的binomial, cauchy, exponential, gamma(erlang), gaussian, uniform integers, uniform, poisson, 还有一个数学函数log(Gamma(x))。
因为参考了一些其它资料，我不知道该怎么发布这个东西，不过我有NRIC的正版授权，因此自己用应该是完全合法的。
使用的时候比较简单，比如需要生成200个正态分布的随机数：

1. 
   
2. double x[200];
   
3. Rand engine; // default to use current time as random seed
   
4. for (int i = 0; i < 200; ++i)
   
5.     x[i] = engine.gasdev();
   

复制代码

// rand.h

1. 
   
2. // -*- C++ -*-
   
3. #ifndef RAND_H
   
4. #define RAND_H
   
5. 
   
6. #include <cmath>
   
7. #include <cstdlib>
   
8. 
   
9. class Rand {
   
10.     static const double PI;
    
11.     unsigned short xsubi[3];
    
12.     void set_seed(long seed = 0);
    
13. 
    
14. public:
    
15.     Rand(long seed = 0) {set_seed(seed);}
    
16.     static double log_gamma(double xx); // returns log(Gamma(xx))
    
17. 
    
18.     // supported random deviates
    
19.     double bnldev(double pp, int n); // binomial
    
20.     double caudev();                 // cauchy
    
21.     double caudev(double location, double scale);
    
22.     double expdev();                 // exponential
    
23.     double expdev(double intensity);
    
24.     double gamdev(int shape);        // gamma(erlang)
    
25.     double gamdev(int shape, double intensity);
    
26.     double gasdev();                 // gaussian
    
27.     double gasdev(double mean, double std);
    
28.     int    intdev(int end_point);    // uniform integeters [0, end_point)
    
29.     int    intdev(int start_point, int end_point);
    
30.     double poidev(double mean);      // poisson
    
31.     double unidev();                 // uniform [0, 1)
    
32.     double unidev(double start_point, double end_point);
    
33. };
    
34. 
    
35. inline double Rand::caudev()
    
36. {
    
37.     return ::tan(PI * (unidev() - 0.5));
    
38. }
    
39. 
    
40. inline double Rand::caudev(double location, double scale)
    
41. {
    
42.     return location + scale * caudev();
    
43. }
    
44. 
    
45. inline double Rand::expdev(double intensity)
    
46. {
    
47.     return expdev() / intensity;
    
48. }
    
49. 
    
50. inline double Rand::gamdev(int shape, double intensity)
    
51. {
    
52.     return gamdev(shape) / intensity;
    
53. }
    
54. 
    
55. inline double Rand::gasdev(double mean, double std)
    
56. {
    
57.     return mean + std * gasdev();
    
58. }
    
59. 
    
60. inline int Rand::intdev(int end_point)
    
61. {
    
62.     return static_cast<int>(::floor(unidev() * end_point));
    
63. }
    
64. 
    
65. inline int Rand::intdev(int start_point, int end_point)
    
66. {
    
67.     return start_point + intdev(end_point - start_point);
    
68. }
    
69. 
    
70. inline double Rand::unidev()
    
71. {
    
72.     return ::erand48(xsubi);
    
73. }
    
74. 
    
75. inline double Rand::unidev(double start_point, double end_point)
    
76. {
    
77.     return start_point + (end_point - start_point) * unidev();
    
78. }
    
79. 
    
80. #endif
    

复制代码

// rand.cc

1. 
   
2. #include "rand.h"
   
3. #include <ctime>
   
4. 
   
5. const double Rand::PI = 3.14159265358979324;
   
6. 
   
7. void Rand::set_seed(long seed)
   
8. {
   
9.     long rand_seed = seed;
   
10. 
    
11.     if (seed == 0)
    
12.         rand_seed = static_cast<long>(::time(0));
    
13.     for (int i = 0; i < 3; ++i) {
    
14.         xsubi[i] = rand_seed & 0xFFFF;
    
15.         rand_seed >>= 16;
    
16.     }
    
17.     // the first few numbers are somewhat related with the seed
    
18.     erand48(xsubi);
    
19.     erand48(xsubi);
    
20.     erand48(xsubi);
    
21.     erand48(xsubi);
    
22. }
    
23. 
    
24. double Rand::log_gamma(double xx)
    
25. {
    
26.     double x, y, tmp, ser;
    
27.     static double cof[6] = {76.18009172947146,
    
28.                             -86.50532032941677,
    
29.                             24.01409824083091,
    
30.                             -1.231739572450155,
    
31.                             0.1208650973866179e-2,
    
32.                             -0.5395239384953e-5};
    
33.     y = x = xx;
    
34.     tmp = x + 5.5;
    
35.     tmp -= (x + 0.5) * ::log(tmp);
    
36.     ser = 1.000000000190015;
    
37.     for (int j = 0; j <= 5; ++j)
    
38.         ser += cof[j] / ++y;
    
39.     return -tmp + ::log(2.5066282746310005 * ser / x);
    
40. }
    
41. 
    
42. double Rand::bnldev(double pp, int n)
    
43. {
    
44.     int j;
    
45.     static int nold = -1;
    
46.     double am, em, g, angle, p, bnl, sq, t, y;
    
47.     static double pold = -1.0, pc, plog, pclog, en, oldg;
    
48. 
    
49.     p = (pp <= 0.5 ? pp : 1.0 - pp);
    
50.     am = n * p;
    
51.     if (n < 25) {
    
52.         bnl = 0.0;
    
53.         for (j = 1; j <= n; ++j)
    
54.             if (unidev() < p)
    
55.                 ++bnl;
    
56.     } else if (am < 1.0) {
    
57.         g = ::exp(-am);
    
58.         t = 1.0;
    
59.         for (j = 0; j <= n; ++j) {
    
60.             t *= unidev();
    
61.             if (t < g)
    
62.                 break;
    
63.         }
    
64.         bnl = (j <= n ? j : n);
    
65.     } else {
    
66.         if (n != nold) {
    
67.             en = n;
    
68.             oldg = log_gamma(en + 1.0);
    
69.             nold = n;
    
70.         } if (p != pold) {
    
71.             pc = 1.0 - p;
    
72.             plog = ::log(p);
    
73.             pclog = ::log(pc);
    
74.             pold = p;
    
75.         }
    
76.         sq = ::sqrt(2.0 * am * pc);
    
77.         do {
    
78.             do {
    
79.                 angle = PI * unidev();
    
80.                 y = ::tan(angle);
    
81.                 em = sq * y + am;
    
82.             } while (em < 0.0 || em >= (en + 1.0));
    
83.             em = ::floor(em);
    
84.             t = 1.2 * sq * (1.0 + y * y) *
    
85.                 ::exp(oldg - log_gamma(em + 1.0) -
    
86.                       log_gamma(en - em + 1.0) +
    
87.                       em * plog + (en - em) * pclog);
    
88.         } while (unidev() > t);
    
89.         bnl = em;
    
90.     }
    
91.     if (p != pp)
    
92.         bnl = n - bnl;
    
93.     return bnl;
    
94. }
    
95. 
    
96. double Rand::expdev()
    
97. {
    
98.     double dum;
    
99. 
    
100.     do {
     
101.         dum = unidev();
     
102.     } while (dum == 0.0);
     
103.     return -::log(dum);
     
104. }
     
105. 
     
106. double Rand::gamdev(int shape)
     
107. {
     
108.     double am, e, s, v1, v2, x, y;
     
109. 
     
110.     if (shape < 6) {
     
111.         x = 1.0;
     
112.         for (int j = 1; j <= shape; ++j)
     
113.             x *= unidev();
     
114.         x = -::log(x);
     
115.     } else {
     
116.         do {
     
117.             do {
     
118.                 do {
     
119.                     v1 = unidev();
     
120.                     v2 = 2.0 * unidev() - 1.0;
     
121.                 } while (v1 * v1 + v2 * v2 > 1.0);
     
122.                 y = v2 / v1;
     
123.                 am = shape - 1;
     
124.                 s = sqrt(2.0 * am + 1.0);
     
125.                 x = s * y + am;
     
126.             } while (x <= 0.0);
     
127.             e = (1.0 + y * y) * ::exp(am * ::log(x / am) - s * y);
     
128.         } while (unidev() > e);
     
129.     }
     
130.     return x;
     
131. }
     
132. 
     
133. double Rand::gasdev()
     
134. {
     
135.     static int iset = 0;
     
136.     static double gset;
     
137.     double fac, rsq, v1, v2;
     
138. 
     
139.     if (iset == 0) {
     
140.         do {
     
141.             v1 = 2.0 * unidev() - 1.0;
     
142.             v2 = 2.0 * unidev() - 1.0;
     
143.             rsq = v1 * v1 + v2 * v2;
     
144.         } while (rsq >= 1.0 || rsq == 0.0);
     
145.         fac = sqrt(-2.0 * ::log(rsq) / rsq);
     
146.         gset = v1 * fac;
     
147.         iset = 1;
     
148.         return v2 * fac;
     
149.     }
     
150.     iset = 0;
     
151.     return gset;
     
152. }
     
153. 
     
154. double Rand::poidev(double mean)
     
155. {
     
156.     static double sq, alxm, g, oldm = -1.0;
     
157.     double em, t, y;
     
158. 
     
159.     if (mean < 12.0) {
     
160.         if (mean != oldm) {
     
161.             oldm = mean;
     
162.             g = ::exp(-mean);
     
163.         }
     
164.         em = -1;
     
165.         t = 1.0;
     
166.         do {
     
167.             ++em;
     
168.             t *= unidev();
     
169.         } while (t > g);
     
170.     } else {
     
171.         if (mean != oldm) {
     
172.             oldm = mean;
     
173.             sq = sqrt(2.0 * mean);
     
174.             alxm = ::log(mean);
     
175.             g = mean * alxm - log_gamma(mean + 1.0);
     
176.         }
     
177.         do {
     
178.             do {
     
179.                 y = ::tan(PI * unidev());
     
180.                 em = sq * y + mean;
     
181.             } while (em < 0.0);
     
182.             em = ::floor(em);
     
183.             t = 0.9 * (1.0 + y * y) *
     
184.                 ::exp(em * alxm - log_gamma(em + 1.0) - g);
     
185.         } while (unidev() > t);
     
186.     }
     
187.     return em;
     
188. }
     

复制代码