SoftSUSY is hosted by Hepforge, IPPP Durham
SOFTSUSY  4.1
numerics.h
Go to the documentation of this file.
1 
12 #ifndef NUMERICS_H
13 #define NUMERICS_H
14 
17 
18 #include "utils.h"
19 #include "mycomplex.h"
20 #include "def.h"
21 #include <iostream>
22 #include "linalg.h"
23 #include "dilogwrap.h"
24 using namespace softsusy;
25 
28 double dgauss(double (*f)(double x), double a, double b, double eps);
29 
32 double dgauss(double (*f)(double x, const DoubleVector & v),
33  const DoubleVector & v, double a, double b, double eps);
34 
36 double accurateSqrt1Plusx(double x);
37 
41 void rungeKuttaStep(const DoubleVector & y, const DoubleVector & dydx,
42  double x, double h, DoubleVector & yout, DoubleVector & yerr,
43  DoubleVector (*derivs)(double, const DoubleVector &));
44 
46 int odeStepper(DoubleVector & y, const DoubleVector & dydx, double *x, double
47  htry, double eps, DoubleVector & yscal, double *hdid,
48  double *hnext,
49  DoubleVector (*derivs)(double, const DoubleVector &));
50 
53 double lnLPoisson(unsigned k, double lambda);
54 
57 double LPoisson(unsigned k, double lambda);
58 
60 int integrateOdes(DoubleVector & ystart, double x1, double x2, double eps,
61  double h1, double hmin,
62  DoubleVector (*derivs)(double, const DoubleVector &),
63  int (*rkqs)
64  (DoubleVector & y, const DoubleVector & dydx, double *x,
65  double htry, double eps, DoubleVector & yscal, double
66  *hdid, double *hnext,
67  DoubleVector (*derivs)(double, const DoubleVector &)));
68 
71 double calcDerivative(double (*func)(double),
72  double x, double h, double *err);
73 
76 double findMinimum(double ax, double bx, double cx, double (*f)(double),
77  double tol, double & xmin);
78 
79 void shft2(double & a, double & b, double & c);
80 void shft3(double & a, double & b, double & c, double & d);
82 
84 double integrandThreshbn(double x);
86 double bIntegral(int n, double p, double m1, double m2, double mt);
87 DoubleVector dd(double x, const DoubleVector & y);
88 
90 double b0(double p, double m1, double m2, double q);
92 double b1(double p, double m1, double m2, double q);
94 double b22(double p, double m1, double m2, double q);
96 double c0(double m1, double m2, double m3);
98 double d27(double m1, double m2, double m3, double m4);
100 double d0(double m1, double m2, double m3, double m4);
102 inline double a0(double m, double q) {
103  if (fabs(m) < EPSTOL) return 0.;
104  return sqr(m) * (1.0 - 2. * log(abs(m / q)));
105 }
107 inline double ffn(double p, double m1, double m2, double q) {
108  return a0(m1, q) - 2.0 * a0(m2, q) -
109  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
110  b0(p, m1, m2, q);
111 }
113 inline double gfn(double p, double m1, double m2, double q) {
114  return (sqr(p) - sqr(m1) - sqr(m2)) * b0(p, m1, m2, q) - a0(m1, q)
115  - a0(m2, q);
116 }
118 inline double hfn(double p, double m1, double m2, double q) {
119  return 4.0 * b22(p, m1, m2, q) + gfn(p, m1, m2, q);
120 }
122 inline double b22bar(double p, double m1, double m2, double q) {
123  return b22(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
124 }
125 
126 Complex b0c(double p, double m1, double m2, double q);
128 Complex b1c(double p, double m1, double m2, double q);
130 Complex b22c(double p, double m1, double m2, double q);
133 inline Complex ffnc(double p, double m1, double m2, double q) {
134  return a0(m1, q) - 2.0 * a0(m2, q) -
135  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
136  b0c(p, m1, m2, q);
137 }
139 inline Complex gfnc(double p, double m1, double m2, double q) {
140  return (sqr(p) - sqr(m1) - sqr(m2)) * b0c(p, m1, m2, q) - a0(m1, q)
141  - a0(m2, q);
142 }
144 inline Complex hfnc(double p, double m1, double m2, double q) {
145  return 4.0 * b22c(p, m1, m2, q) + gfnc(p, m1, m2, q);
146 }
148 inline Complex b22barc(double p, double m1, double m2, double q) {
149  return b22c(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
150 }
151 
152 /*inline double fB(const Complex & x) {
154  // if (x.real() < EPSTOL * 10. && x.imag() < EPSTOL * 10.) return -1.;
155  return (log(Complex(1.0 - x)) - x * log(Complex(1.0 - 1.0 / x))
156  - Complex(1.0)).real();
157  } */
158 double fB(const Complex & x);
159 
160 Complex fBc(const Complex & x);
161 double dilogarg(double t);
162 double dilog(double x);
163 
164 double integrandThreshbnr(double x);
165 Complex fnfn(double x);
166 
169 double gasdev(long & idum);
172 double truncGaussWidthHalf(long & idum);
175 double ran1(long & idum);
178 double cauchyRan(long & idum);
179 
183 int bin(double data, double start, double end, int numBins);
184 
186 double logOfSum(double a, double b);
187 
189 DoubleVector getRandomDirection(int n, int & numChanged, long & idum);
190 
194 double calcCL(double cl, const DoubleVector & l);
195 
199 double calc1dFraction(double y, const DoubleVector & l);
200 
203 double fps(double z);
204 double fs(double z);
205 double ffbar(double z);
206 
207 Complex dilog(const Complex & x);
208 std::complex<long double> dilogcl(const std::complex<long double>&);
209 
211 double trapzd(double (*func)(double), double a, double b, int n,
212  double tol = 1.0e-3);
214 double qtrap(double (*func)(double), double a, double b, double tol);
216 double midpnt(double (*func)(double), double a, double b, int n);
219 double qromb(double (*func)(double), double a, double b, double EPS);
221 double qromb2(double (*func)(double), double a, double b, double EPS);
222 
225 double sumOfExp(double a, double b);
226 
229 void polint(const DoubleVector& xa, const DoubleVector & ya, double x,
230  double & y, double & dy);
231 
236 DoubleMatrix display3x3RealMixing(double theta12, double theta13,
237  double theta23, double d);
238 
240 void getAngles(const DoubleMatrix & v, double & t12, double & t13,
241  double & t23, double & d);
244 bool midPtStep(DoubleVector & xi,
245  DoubleVector (*derivs)(double t, const DoubleVector & v),
246  double tInitial, double tStep);
250  DoubleVector (*derivs)(double t,
251  const DoubleVector & v),
252  double tInitial, double tFinal, int numSteps);
253 
255 double fin(double mm1, double mm2);
256 double den(double a, int b);
257 double log1minusx(double x);
260 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
261 double zriddr(double (*func)(double), double x1, double x2, double xacc);
267 void fdjac(int n, DoubleVector x, const DoubleVector & fvec, DoubleMatrix & df,
268  void (*vecfunc)(const DoubleVector &, DoubleVector &));
271 bool lnsrch(const DoubleVector & xold, double fold, const DoubleVector & g,
272  DoubleVector & p,
273  DoubleVector & x, double & f, double stpmax,
274  void (*vecfunc)(const DoubleVector &, DoubleVector &),
275  DoubleVector & fvec);
276 /* allocate an int vector with subscript range v[nl..nh] */
277 int *ivector(long nl, long nh);
278 /* free an int vector allocated with ivector() */
279 void free_ivector(int *v, long nl, long nh);
280 void lubksb(const DoubleMatrix & a, int n, int *indx, DoubleVector & b);
287 bool newt(DoubleVector & x,
288  void (*vecfunc)(const DoubleVector &, DoubleVector &));
290 void load(float x, const DoubleVector & v, DoubleVector & y);
293 void score(float x, const DoubleVector & y, DoubleVector & f);
299 //void shoot(const DoubleVector & v, DoubleVector & f);
306 void qrdcmp(DoubleMatrix & a, int n, DoubleVector & c, DoubleVector & d,
307  int & sing);
308 void qrupdt(DoubleMatrix & r, DoubleMatrix & qt, int n,
309  DoubleVector & u, DoubleVector & v);
313 void rotate(DoubleMatrix & r, DoubleMatrix & qt, int n, int i, float a,
314  float b);
315 void rsolv(const DoubleMatrix & a, int n, const DoubleVector & d,
316  DoubleVector & b);
317 
319  // returns sqrt(f) for f>0
320  double ccbSqrt(double f);
322  // returns sqrt(f) for f>0 or -sqrt(|f|) for f<0
323  double signedSqrt(double f);
325  double signedSqr(double f);
326 
330 double kinFn(double m1, double m2, double m3);
331 double lambda(double a, double b, double c);
332 #endif
333 
void rungeKuttaStep(const DoubleVector &y, const DoubleVector &dydx, double x, double h, DoubleVector &yout, DoubleVector &yerr, DoubleVector(*derivs)(double, const DoubleVector &))
Definition: numerics.cpp:302
global variable declaration
Definition: def.cpp:13
int bin(double data, double start, double end, int numBins)
Definition: numerics.cpp:1038
Complex log(const Complex &a)
Principal log.
Definition: mycomplex.h:67
double log1minusx(double x)
This function DOESNT WORK YET!!!
Definition: numerics.cpp:167
double lnLPoisson(unsigned k, double lambda)
Definition: numerics.cpp:248
double a0(double m, double q)
inlined PV functions: real part
Definition: numerics.h:102
DoubleVector is of variable length, and contains double precision.
Definition: linalg.h:35
void rotate(DoubleMatrix &r, DoubleMatrix &qt, int n, int i, float a, float b)
Definition: numerics.cpp:2217
double b22bar(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:122
double trapzd(double(*func)(double), double a, double b, int n, double tol=1.0e-3)
Trapezoidal function integration to accuracy tol.
Definition: numerics.cpp:1339
double zriddr(double(*func)(double), double x1, double x2, double xacc)
Definition: numerics.cpp:1698
void score(float x, const DoubleVector &y, DoubleVector &f)
std::complex< long double > dilogcl(const std::complex< long double > &)
Complex dilogarithm .
Definition: numerics.cpp:1156
switches (options) and parameters such as default fermion masses, required accuracy etc ...
double qromb(double(*func)(double), double a, double b, double EPS)
Definition: numerics.cpp:1490
DoubleVector and DoubleMatrix classes of doubles and operations between them, complexified copies als...
double findMinimum(double ax, double bx, double cx, double(*f)(double), double tol, double &xmin)
Definition: numerics.cpp:390
Dilog function.
double b0(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:494
void polint(const DoubleVector &xa, const DoubleVector &ya, double x, double &y, double &dy)
Definition: numerics.cpp:1541
double b1(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:633
double accurateSqrt1Plusx(double x)
calculate root(1+x), where x<<1 accurately
Definition: numerics.cpp:146
void shft2(double &a, double &b, double &c)
double truncGaussWidthHalf(long &idum)
Definition: numerics.cpp:972
double gfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:113
A few handy bits and pieces - little mathematical functions and the like.
int integrateOdes(DoubleVector &ystart, double x1, double x2, double eps, double h1, double hmin, DoubleVector(*derivs)(double, const DoubleVector &), int(*rkqs)(DoubleVector &y, const DoubleVector &dydx, double *x, double htry, double eps, DoubleVector &yscal, double *hdid, double *hnext, DoubleVector(*derivs)(double, const DoubleVector &)))
Organises integration of 1st order system of ODEs.
Definition: numerics.cpp:189
double logOfSum(double a, double b)
Adds logs of two numbers in a more careful way that avoids underflow.
Definition: numerics.cpp:1044
double calcCL(double cl, const DoubleVector &l)
Definition: numerics.cpp:1267
double calc1dFraction(double y, const DoubleVector &l)
Definition: numerics.cpp:1297
bool lnsrch(const DoubleVector &xold, double fold, const DoubleVector &g, DoubleVector &p, DoubleVector &x, double &f, double stpmax, void(*vecfunc)(const DoubleVector &, DoubleVector &), DoubleVector &fvec)
Definition: numerics.cpp:1769
double b22(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:735
double signedSqr(double f)
returns f * f * sign(f)
Definition: numerics.cpp:2264
double integrandThreshbn(double x)
For calculation of PV functions.
Complex hfnc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:144
Complex fBc(const Complex &x)
Definition: numerics.cpp:475
void getAngles(const DoubleMatrix &v, double &t12, double &t13, double &t23, double &d)
Given a matrix v, determines which angles gave it.
Definition: numerics.cpp:1606
bool integrateReversibly(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tFinal, int numSteps)
Definition: numerics.cpp:1650
Complex b22barc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:148
double ccbSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2260
Matrix from 1..rows, 1..cols of double values.
Definition: linalg.h:214
double sqr(double a)
square of a number
Definition: utils.h:49
bool midPtStep(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tStep)
Definition: numerics.cpp:1619
double signedSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2262
int odeStepper(DoubleVector &y, const DoubleVector &dydx, double *x, double htry, double eps, DoubleVector &yscal, double *hdid, double *hnext, DoubleVector(*derivs)(double, const DoubleVector &))
organises the variable step-size for Runge-Kutta evolution
Definition: numerics.cpp:264
double gasdev(long &idum)
Definition: numerics.cpp:980
void fdjac(int n, DoubleVector x, const DoubleVector &fvec, DoubleMatrix &df, void(*vecfunc)(const DoubleVector &, DoubleVector &))
void shft3(double &a, double &b, double &c, double &d)
a=b, b=c and c=d
bool newt(DoubleVector &x, void(*vecfunc)(const DoubleVector &, DoubleVector &))
More work can be done on this: get rid of int n and in subfunctions too.
Definition: numerics.cpp:1910
void load(float x, const DoubleVector &v, DoubleVector &y)
calculates the n-vector y, given freely specifiable values v(1..n2) at x1
double calcDerivative(double(*func)(double), double x, double h, double *err)
Definition: numerics.cpp:349
double cauchyRan(long &idum)
Definition: numerics.cpp:1032
double sumOfExp(double a, double b)
Definition: numerics.cpp:1055
Complex b1c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:692
double dgauss(double(*f)(double x), double a, double b, double eps)
Definition: numerics.cpp:20
void qrdcmp(DoubleMatrix &a, int n, DoubleVector &c, DoubleVector &d, int &sing)
Definition: numerics.cpp:2165
double midpnt(double(*func)(double), double a, double b, int n)
Mid point function integration.
Definition: numerics.cpp:1376
double d27(double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:886
complex numbers and operators between them
double m1
decay global variable declarations
Definition: decays.cpp:14
double qtrap(double(*func)(double), double a, double b, double tol)
Driver for integration.
Definition: numerics.cpp:1357
double fin(double mm1, double mm2)
useful for 2-loop mb/mt corrections
Definition: numerics.cpp:1672
double c0(double m1, double m2, double m3)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:898
double fps(double z)
Definition: numerics.cpp:1308
double qromb2(double(*func)(double), double a, double b, double EPS)
Identical copy to facilitate double integral.
Definition: numerics.cpp:1514
double kinFn(double m1, double m2, double m3)
Definition: numerics.cpp:2276
void lubksb(const DoubleMatrix &a, int n, int *indx, DoubleVector &b)
Get rid of int n.
Definition: numerics.cpp:1831
Complex ffnc(double p, double m1, double m2, double q)
Definition: numerics.h:133
DoubleMatrix display3x3RealMixing(double theta12, double theta13, double theta23, double d)
Definition: numerics.cpp:1574
const double EPSTOL
underflow accuracy
Definition: def.h:33
Complex gfnc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:139
double bIntegral(int n, double p, double m1, double m2, double mt)
Returns real part of b function, less accurate than analytic expressions.
Definition: numerics.cpp:452
double hfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:118
double ran1(long &idum)
Definition: numerics.cpp:1001
double fB(const Complex &x)
Definition: numerics.cpp:461
Complex b22c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:790
double ffn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:107
DoubleVector getRandomDirection(int n, int &numChanged, long &idum)
returns a random direction: total number of dimensions=n
Definition: numerics.cpp:1245
double LPoisson(unsigned k, double lambda)
Definition: numerics.cpp:256
drop-in replacement for the original home-grown Complex class
Definition: mycomplex.h:17
double d0(double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:851
Complex b0c(double p, double m1, double m2, double q)
Definition: numerics.cpp:563