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 #ifdef USE_LOOPTOOLS
25 #include "clooptools.h"
26 #endif
27 using namespace softsusy;
28 
31 double dgauss(double (*f)(double x), double a, double b, double eps);
32 
35 double dgauss(double (*f)(double x, const DoubleVector & v),
36  const DoubleVector & v, double a, double b, double eps);
37 
39 double accurateSqrt1Plusx(double x);
40 
44 void rungeKuttaStep(const DoubleVector & y, const DoubleVector & dydx,
45  double x, double h, DoubleVector & yout, DoubleVector & yerr,
46  DoubleVector (*derivs)(double, const DoubleVector &));
47 
49 int odeStepper(DoubleVector & y, const DoubleVector & dydx, double *x, double
50  htry, double eps, DoubleVector & yscal, double *hdid,
51  double *hnext,
52  DoubleVector (*derivs)(double, const DoubleVector &));
53 
56 double lnLPoisson(unsigned k, double lambda);
57 
60 double LPoisson(unsigned k, double lambda);
61 
63 int integrateOdes(DoubleVector & ystart, double x1, double x2, double eps,
64  double h1, double hmin,
65  DoubleVector (*derivs)(double, const DoubleVector &),
66  int (*rkqs)
67  (DoubleVector & y, const DoubleVector & dydx, double *x,
68  double htry, double eps, DoubleVector & yscal, double
69  *hdid, double *hnext,
70  DoubleVector (*derivs)(double, const DoubleVector &)));
71 
74 double calcDerivative(double (*func)(double),
75  double x, double h, double *err);
76 
79 double findMinimum(double ax, double bx, double cx, double (*f)(double),
80  double tol, double & xmin);
81 
82 void shft2(double & a, double & b, double & c);
83 void shft3(double & a, double & b, double & c, double & d);
85 
87 double integrandThreshbn(double x);
89 double bIntegral(int n, double p, double m1, double m2, double mt);
90 DoubleVector dd(double x, const DoubleVector & y);
91 
93 double b0(double p, double m1, double m2, double q);
95 double b1(double p, double m1, double m2, double q);
97 double b22(double p, double m1, double m2, double q);
99 double c0(double m1, double m2, double m3);
101 double d27(double m1, double m2, double m3, double m4);
103 double d0(double m1, double m2, double m3, double m4);
105 inline double a0(double m, double q) {
106  if (fabs(m) < EPSTOL) return 0.;
107  return sqr(m) * (1.0 - 2. * log(abs(m / q)));
108 }
110 inline double ffn(double p, double m1, double m2, double q) {
111  return a0(m1, q) - 2.0 * a0(m2, q) -
112  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
113  b0(p, m1, m2, q);
114 }
116 inline double gfn(double p, double m1, double m2, double q) {
117  return (sqr(p) - sqr(m1) - sqr(m2)) * b0(p, m1, m2, q) - a0(m1, q)
118  - a0(m2, q);
119 }
121 inline double hfn(double p, double m1, double m2, double q) {
122  return 4.0 * b22(p, m1, m2, q) + gfn(p, m1, m2, q);
123 }
125 inline double b22bar(double p, double m1, double m2, double q) {
126  return b22(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
127 }
128 
129 Complex b0c(double p, double m1, double m2, double q);
131 Complex b1c(double p, double m1, double m2, double q);
133 Complex b22c(double p, double m1, double m2, double q);
136 inline Complex ffnc(double p, double m1, double m2, double q) {
137  return a0(m1, q) - 2.0 * a0(m2, q) -
138  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
139  b0c(p, m1, m2, q);
140 }
142 inline Complex gfnc(double p, double m1, double m2, double q) {
143  return (sqr(p) - sqr(m1) - sqr(m2)) * b0c(p, m1, m2, q) - a0(m1, q)
144  - a0(m2, q);
145 }
147 inline Complex hfnc(double p, double m1, double m2, double q) {
148  return 4.0 * b22c(p, m1, m2, q) + gfnc(p, m1, m2, q);
149 }
151 inline Complex b22barc(double p, double m1, double m2, double q) {
152  return b22c(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
153 }
154 
155 /*inline double fB(const Complex & x) {
157  // if (x.real() < EPSTOL * 10. && x.imag() < EPSTOL * 10.) return -1.;
158  return (log(Complex(1.0 - x)) - x * log(Complex(1.0 - 1.0 / x))
159  - Complex(1.0)).real();
160  } */
161 double fB(const Complex & x);
162 
163 Complex fBc(const Complex & x);
164 double dilogarg(double t);
165 double dilog(double x);
166 
167 double integrandThreshbnr(double x);
168 Complex fnfn(double x);
169 
172 double gasdev(long & idum);
175 double truncGaussWidthHalf(long & idum);
178 double ran1(long & idum);
181 double cauchyRan(long & idum);
182 
186 int bin(double data, double start, double end, int numBins);
187 
189 double logOfSum(double a, double b);
190 
192 DoubleVector getRandomDirection(int n, int & numChanged, long & idum);
193 
197 double calcCL(double cl, const DoubleVector & l);
198 
202 double calc1dFraction(double y, const DoubleVector & l);
203 
206 double fps(double z);
207 double fs(double z);
208 double ffbar(double z);
209 
210 Complex dilog(const Complex & x);
211 std::complex<long double> dilogcl(const std::complex<long double>&);
212 
214 double trapzd(double (*func)(double), double a, double b, int n,
215  double tol = 1.0e-3);
217 double qtrap(double (*func)(double), double a, double b, double tol);
219 double midpnt(double (*func)(double), double a, double b, int n);
222 double qromb(double (*func)(double), double a, double b, double EPS);
224 double qromb2(double (*func)(double), double a, double b, double EPS);
225 
228 double sumOfExp(double a, double b);
229 
232 void polint(const DoubleVector& xa, const DoubleVector & ya, double x,
233  double & y, double & dy);
234 
239 DoubleMatrix display3x3RealMixing(double theta12, double theta13,
240  double theta23, double d);
241 
243 void getAngles(const DoubleMatrix & v, double & t12, double & t13,
244  double & t23, double & d);
247 bool midPtStep(DoubleVector & xi,
248  DoubleVector (*derivs)(double t, const DoubleVector & v),
249  double tInitial, double tStep);
253  DoubleVector (*derivs)(double t,
254  const DoubleVector & v),
255  double tInitial, double tFinal, int numSteps);
256 
258 double fin(double mm1, double mm2);
259 double den(double a, int b);
260 double log1minusx(double x);
263 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
264 double zriddr(double (*func)(double), double x1, double x2, double xacc);
270 void fdjac(int n, DoubleVector x, const DoubleVector & fvec, DoubleMatrix & df,
271  void (*vecfunc)(const DoubleVector &, DoubleVector &));
274 bool lnsrch(const DoubleVector & xold, double fold, const DoubleVector & g,
275  DoubleVector & p,
276  DoubleVector & x, double & f, double stpmax,
277  void (*vecfunc)(const DoubleVector &, DoubleVector &),
278  DoubleVector & fvec);
279 /* allocate an int vector with subscript range v[nl..nh] */
280 int *ivector(long nl, long nh);
281 /* free an int vector allocated with ivector() */
282 void free_ivector(int *v, long nl, long nh);
283 void lubksb(const DoubleMatrix & a, int n, int *indx, DoubleVector & b);
290 bool newt(DoubleVector & x,
291  void (*vecfunc)(const DoubleVector &, DoubleVector &));
293 void load(float x, const DoubleVector & v, DoubleVector & y);
296 void score(float x, const DoubleVector & y, DoubleVector & f);
302 //void shoot(const DoubleVector & v, DoubleVector & f);
309 void qrdcmp(DoubleMatrix & a, int n, DoubleVector & c, DoubleVector & d,
310  int & sing);
311 void qrupdt(DoubleMatrix & r, DoubleMatrix & qt, int n,
312  DoubleVector & u, DoubleVector & v);
316 void rotate(DoubleMatrix & r, DoubleMatrix & qt, int n, int i, float a,
317  float b);
318 void rsolv(const DoubleMatrix & a, int n, const DoubleVector & d,
319  DoubleVector & b);
320 
322  // returns sqrt(f) for f>0
323  double ccbSqrt(double f);
325  // returns sqrt(f) for f>0 or -sqrt(|f|) for f<0
326  double signedSqrt(double f);
328  double signedSqr(double f);
329 
333 double kinFn(double m1, double m2, double m3);
334 double lambda(double a, double b, double c);
335 #endif
336 
void rungeKuttaStep(const DoubleVector &y, const DoubleVector &dydx, double x, double h, DoubleVector &yout, DoubleVector &yerr, DoubleVector(*derivs)(double, const DoubleVector &))
Definition: numerics.cpp:294
global variable declaration
Definition: def.cpp:13
int bin(double data, double start, double end, int numBins)
Definition: numerics.cpp:1030
Complex log(const Complex &a)
Principal log.
Definition: mycomplex.h:67
double log1minusx(double x)
This function DOESNT WORK YET!!!
Definition: numerics.cpp:159
double lnLPoisson(unsigned k, double lambda)
Definition: numerics.cpp:240
double a0(double m, double q)
inlined PV functions: real part
Definition: numerics.h:105
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:2209
double b22bar(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:125
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:1331
double zriddr(double(*func)(double), double x1, double x2, double xacc)
Definition: numerics.cpp:1690
void score(float x, const DoubleVector &y, DoubleVector &f)
std::complex< long double > dilogcl(const std::complex< long double > &)
Complex dilogarithm .
Definition: numerics.cpp:1148
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:1482
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:382
Dilog function.
double b0(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:486
void polint(const DoubleVector &xa, const DoubleVector &ya, double x, double &y, double &dy)
Definition: numerics.cpp:1533
double b1(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:625
double accurateSqrt1Plusx(double x)
calculate root(1+x), where x<<1 accurately
Definition: numerics.cpp:138
void shft2(double &a, double &b, double &c)
double truncGaussWidthHalf(long &idum)
Definition: numerics.cpp:964
double gfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:116
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:181
double logOfSum(double a, double b)
Adds logs of two numbers in a more careful way that avoids underflow.
Definition: numerics.cpp:1036
double calcCL(double cl, const DoubleVector &l)
Definition: numerics.cpp:1259
double calc1dFraction(double y, const DoubleVector &l)
Definition: numerics.cpp:1289
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:1761
double b22(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:727
double signedSqr(double f)
returns f * f * sign(f)
Definition: numerics.cpp:2256
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:147
Complex fBc(const Complex &x)
Definition: numerics.cpp:467
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:1598
bool integrateReversibly(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tFinal, int numSteps)
Definition: numerics.cpp:1642
Complex b22barc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:151
double ccbSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2252
Matrix from 1..rows, 1..cols of double values.
Definition: linalg.h:211
double sqr(double a)
square of a number
Definition: utils.h:53
bool midPtStep(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tStep)
Definition: numerics.cpp:1611
double signedSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2254
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:256
double gasdev(long &idum)
Definition: numerics.cpp:972
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:1902
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:341
double cauchyRan(long &idum)
Definition: numerics.cpp:1024
double sumOfExp(double a, double b)
Definition: numerics.cpp:1047
Complex b1c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:684
double dgauss(double(*f)(double x), double a, double b, double eps)
Definition: numerics.cpp:12
void qrdcmp(DoubleMatrix &a, int n, DoubleVector &c, DoubleVector &d, int &sing)
Definition: numerics.cpp:2157
double midpnt(double(*func)(double), double a, double b, int n)
Mid point function integration.
Definition: numerics.cpp:1368
double d27(double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:878
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:1349
double fin(double mm1, double mm2)
useful for 2-loop mb/mt corrections
Definition: numerics.cpp:1664
double c0(double m1, double m2, double m3)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:890
double fps(double z)
Definition: numerics.cpp:1300
double qromb2(double(*func)(double), double a, double b, double EPS)
Identical copy to facilitate double integral.
Definition: numerics.cpp:1506
double kinFn(double m1, double m2, double m3)
Definition: numerics.cpp:2268
void lubksb(const DoubleMatrix &a, int n, int *indx, DoubleVector &b)
Get rid of int n.
Definition: numerics.cpp:1823
Complex ffnc(double p, double m1, double m2, double q)
Definition: numerics.h:136
DoubleMatrix display3x3RealMixing(double theta12, double theta13, double theta23, double d)
Definition: numerics.cpp:1566
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:142
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:444
double hfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:121
double ran1(long &idum)
Definition: numerics.cpp:993
double fB(const Complex &x)
Definition: numerics.cpp:453
Complex b22c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:782
double ffn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:110
DoubleVector getRandomDirection(int n, int &numChanged, long &idum)
returns a random direction: total number of dimensions=n
Definition: numerics.cpp:1237
double LPoisson(unsigned k, double lambda)
Definition: numerics.cpp:248
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:843
Complex b0c(double p, double m1, double m2, double q)
Definition: numerics.cpp:555