1 // $Id: mathx.h,v 1.32 2011/04/23 02:02:49 bobgian Exp $ 2 3 /* 4 Copyright 2002 Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein 5 6 This software is distributed free of charge for non-commercial use 7 and is copyrighted. Of course, we do not guarantee that the software 8 works, and are not responsible for any damage you may cause or have. 9 */ 10 11 #ifndef MATHX_H 12 #define MATHX_H 13 14 #include <cmath> 15 #include <vector> 16 17 #include "vectorx.h" 18 19 // if we're on an alpha... does this catch that? 20 #ifdef alphaev56 21 #define MYEXP UnderFlowExp 22 #else 23 #define MYEXP exp 24 #endif 25 26 //------------------------------------------------------------------------------------ 27 // The proper version of isnan and the right way to get access to 28 // it depends on the system and compiler. Here we try to get that 29 // all straightened out. 30 31 // Macosx doesn't currently include isnan properly without this 32 // G-d only knows why it needs iostream to compile properly 33 #ifdef LAMARC_COMPILE_MACOSX 34 #if (__GNUC__ < 4) 35 #include <iostream> 36 extern "C" int isnan(double); 37 #endif 38 #endif 39 40 // the default isnan 41 #define systemSpecificIsnan isnan 42 43 // Microsoft VC++ uses a different isnan function. Unfortunately, 44 // Metrowerks defines _MSC_VER sometimes, but we don't want to 45 // do this in that case 46 #ifdef __MWERKS__ 47 #elif defined _MSC_VER 48 #include <float.h> 49 #define systemSpecificIsnan _isnan 50 #endif 51 52 const double LOG2 = 0.69314718055994530942; // log(2.0). 53 const double LOG10 = 2.30258509299404568401799145468; // log(10.0) 54 const double PI = 3.1415926535897932384626433832795028841972; 55 const double SQRT_PI = 1.7724538509055160272981674833411451827975; // sqrt(PI) 56 const double EULERS_CONSTANT = 0.5772156649015328606065120900824024310422; 57 58 using std::vector; 59 double log_gamma(double x); 60 double my_gamma(double x); 61 double erfc(double x); // the "complementary error function" 62 double incompletegamma(double alpha, double x); 63 double probchi(long df, double chi); 64 double find_chi(long df, double prob); 65 DoubleVec1d gamma_rates(double alpha, long ratenum); 66 double alnorm (double x, int up); 67 double SafeDivide(double num, double denom); 68 double logfac(long n); 69 //long factorial(long n); 70 double factorial(double n); 71 double UnderFlowExp(double pow); 72 bool IsEven(long n); 73 DoubleVec2d Invert(const DoubleVec2d& src); // invert square matrix 74 DoubleVec2d Transpose(const DoubleVec2d& src); // transpose square matrix 75 bool CloseEnough(double a, double b); // are values within epsilon? 76 DoubleVec1d SafeExp(const DoubleVec1d& src); 77 double SafeExpAndSum(const DoubleVec1d& src); 78 double SafeExp(double src); 79 double SafeProductWithExp(double a, double x); 80 double OverflowOfProductWithExp(double a, double x); 81 double UnderflowOfProductWithExp(double a, double x); 82 double SafeExpDiff(double x1, double x2); 83 double OverflowOfSafeExpDiff(double x1, double x2); 84 double UnderflowOfSafeExpDiff(double x1, double x2); 85 double LogZero(double x); 86 double SafeLog(double x); 87 DoubleVec1d SafeLog(DoubleVec1d src); 88 long ChooseRandomFromWeights(DoubleVec1d& weights); 89 double ExpE1(const double& x); 90 double BesselK(double v, double x, double& resultFor_vPlusOne); // == K(v,x) 91 double DvBesselK(double v, double x, double Kvx); 92 double psi(double x); // Euler's psi function == d(log_gamma(x))/dx. 93 void ScaleLargestToZero(DoubleVec1d& unscaled); 94 void ScaleToSumToOne(DoubleVec1d& vec); 95 DoubleVec1d AddValsOfLogs(DoubleVec1d vec1, DoubleVec1d vec2); 96 97 #if 0 // Potentially DEAD CODE (bobgian, Feb 2010) 98 class Gamma 99 { 100 double alpha; 101 double logAlpha; 102 double logAlpha1; 103 double invAlpha9; 104 double sqrAlpha3; 105 106 public: 107 Gamma(); 108 Gamma(double); 109 void SetAlpha(double); 110 DoubleVec1d Calculate(long); // number of divisions 111 double IncompleteGamma(double alpha, double x); 112 double LogGamma(double); 113 double Tail(double); 114 double ProbChi2(long df , double chi); 115 double FindChi2(long df, double prob); 116 }; 117 #endif 118 119 class EigenCalculator 120 { 121 private: 122 123 static const double accuracy; 124 125 // Cosine and sine 126 std::pair<double,double> Coeffs(double x, double y); 127 128 // Givens matrix reduction 129 void Givens(DoubleVec2d& a, long x, long y, long size, 130 const std::pair<double,double> cs, bool userow); 131 132 // Matrix tridiagonalization 133 void Tridiag(DoubleVec2d& a, DoubleVec2d& eigvecs); 134 135 // Intermediate work of eigenvalues/vectors 136 void Shiftqr(DoubleVec2d& a, DoubleVec2d& eigvecs); 137 138 // Modified dot-product 139 void DotProduct(const DoubleVec2d& first, const DoubleVec2d& second, 140 DoubleVec2d& answer); 141 142 public: 143 144 // Driver for eigenvalues/vectors 145 std::pair<DoubleVec1d, DoubleVec2d> Eigen(DoubleVec2d a); 146 147 }; 148 149 #endif // MATHX_H 150 151 //____________________________________________________________________________________ 152