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