1 /*
2   Copyright 2009 Andreas Biegert
3 
4   This file is part of the CS-BLAST package.
5 
6   The CS-BLAST package is free software: you can redistribute it and/or modify
7   it under the terms of the GNU General Public License as published by
8   the Free Software Foundation, either version 3 of the License, or
9   (at your option) any later version.
10 
11   The CS-BLAST package is distributed in the hope that it will be useful,
12   but WITHOUT ANY WARRANTY; without even the implied warranty of
13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14   GNU General Public License for more details.
15 
16   You should have received a copy of the GNU General Public License
17   along with this program.  If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 
21 #ifndef RAN_H_
22 #define RAN_H_
23 
24 #include "utils.h"
25 #include <cmath>
26 
27 namespace cs {
28 
29 // Implementation of the highest quality recommended random number generator
30 // from Numerical Recipes. The period of the generator is 3.138e57.
31 struct Ran {
32   typedef unsigned long long int Ullong;
33 
RanRan34   Ran(Ullong j) : v(4101842887655102017LL), w(1) {
35     u = j ^ v; int64();
36     v = u; int64();
37     w = v; int64();
38   }
39 
int64Ran40   inline Ullong int64() {
41     u = u * 2862933555777941757LL + 7046029254386353087LL;
42     v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
43     w = 4294957665U*(w & 0xffffffff) + (w >> 32);
44     Ullong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
45     return (x + v) ^ w;
46   }
47 
doubRan48   inline double doub() { return 5.42101086242752217E-20 * int64(); }
49 
int32Ran50   inline unsigned int int32() { return (unsigned int)int64(); }
51 
operatorRan52   inline unsigned int operator() (unsigned int n) { return int32() % n; }
53 
54   Ullong u,v,w;
55 };
56 
57 
58 // Normal diistribution generator also from Numerical Recipes.
59 struct Gaussian : public Ran {
60   typedef unsigned long long int Ullong;
61 
GaussianGaussian62   Gaussian(double mmu, double ssig, Ullong i) : Ran(i), mu(mmu), sig(ssig) {}
63 
operatorGaussian64   double operator() () {
65     double u,v,x,y,q;
66     do {
67       u = doub();
68       v = 1.7156*(doub()-0.5);
69       x = u - 0.449871;
70       y = std::abs(v) + 0.386595;
71       q = SQR(x) + y*(0.19600*y-0.25472*x);
72     } while (q > 0.27597
73              && (q > 0.27846 || SQR(v) > -4.*log(u)*SQR(u)));
74     return mu + sig*v/u;
75   }
76 
77   double mu,sig;
78 };
79 
80 }  // namespace cs
81 
82 #endif
83