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 calcDerivative(double (*func)(double, void*),
77  double x, double h, double* err, void* params = NULL);
78 
81 double findMinimum(double ax, double bx, double cx, double (*f)(double),
82  double tol, double & xmin);
83 
84 void shft2(double & a, double & b, double & c);
85 void shft3(double & a, double & b, double & c, double & d);
87 
89 double integrandThreshbn(double x);
91 double bIntegral(int n, double p, double m1, double m2, double mt);
92 DoubleVector dd(double x, const DoubleVector & y);
93 
95 double b0(double p, double m1, double m2, double q);
97 double b1(double p, double m1, double m2, double q);
99 double b22(double p, double m1, double m2, double q);
101 double c0(double m1, double m2, double m3);
103 double d27(double m1, double m2, double m3, double m4);
105 double d0(double m1, double m2, double m3, double m4);
107 inline double a0(double m, double q) {
108  if (fabs(m) < EPSTOL) return 0.;
109  return sqr(m) * (1.0 - 2. * log(abs(m / q)));
110 }
112 inline double ffn(double p, double m1, double m2, double q) {
113  return a0(m1, q) - 2.0 * a0(m2, q) -
114  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
115  b0(p, m1, m2, q);
116 }
118 inline double gfn(double p, double m1, double m2, double q) {
119  return (sqr(p) - sqr(m1) - sqr(m2)) * b0(p, m1, m2, q) - a0(m1, q)
120  - a0(m2, q);
121 }
123 inline double hfn(double p, double m1, double m2, double q) {
124  return 4.0 * b22(p, m1, m2, q) + gfn(p, m1, m2, q);
125 }
127 inline double b22bar(double p, double m1, double m2, double q) {
128  return b22(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
129 }
130 
131 Complex b0c(double p, double m1, double m2, double q);
133 Complex b1c(double p, double m1, double m2, double q);
135 Complex b22c(double p, double m1, double m2, double q);
138 inline Complex ffnc(double p, double m1, double m2, double q) {
139  return a0(m1, q) - 2.0 * a0(m2, q) -
140  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
141  b0c(p, m1, m2, q);
142 }
144 inline Complex gfnc(double p, double m1, double m2, double q) {
145  return (sqr(p) - sqr(m1) - sqr(m2)) * b0c(p, m1, m2, q) - a0(m1, q)
146  - a0(m2, q);
147 }
149 inline Complex hfnc(double p, double m1, double m2, double q) {
150  return 4.0 * b22c(p, m1, m2, q) + gfnc(p, m1, m2, q);
151 }
153 inline Complex b22barc(double p, double m1, double m2, double q) {
154  return b22c(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
155 }
156 
157 /*inline double fB(const Complex & x) {
159  // if (x.real() < EPSTOL * 10. && x.imag() < EPSTOL * 10.) return -1.;
160  return (log(Complex(1.0 - x)) - x * log(Complex(1.0 - 1.0 / x))
161  - Complex(1.0)).real();
162  } */
163 double fB(const Complex & x);
164 
165 Complex fBc(const Complex & x);
166 double dilogarg(double t);
167 double dilog(double x);
168 
169 double integrandThreshbnr(double x);
170 Complex fnfn(double x);
171 
174 double gasdev(long & idum);
177 double truncGaussWidthHalf(long & idum);
180 double ran1(long & idum);
183 double cauchyRan(long & idum);
184 
188 int bin(double data, double start, double end, int numBins);
189 
191 double logOfSum(double a, double b);
192 
194 DoubleVector getRandomDirection(int n, int & numChanged, long & idum);
195 
199 double calcCL(double cl, const DoubleVector & l);
200 
204 double calc1dFraction(double y, const DoubleVector & l);
205 
208 double fps(double z);
209 double fs(double z);
210 double ffbar(double z);
211 
212 Complex dilog(const Complex & x);
213 std::complex<long double> dilogcl(const std::complex<long double>&);
214 
216 double trapzd(double (*func)(double), double a, double b, int n,
217  double tol = 1.0e-3);
219 double qtrap(double (*func)(double), double a, double b, double tol);
221 double midpnt(double (*func)(double), double a, double b, int n);
224 double qromb(double (*func)(double), double a, double b, double EPS);
226 double qromb2(double (*func)(double), double a, double b, double EPS);
227 
230 double sumOfExp(double a, double b);
231 
234 void polint(const DoubleVector& xa, const DoubleVector & ya, double x,
235  double & y, double & dy);
236 
241 DoubleMatrix display3x3RealMixing(double theta12, double theta13,
242  double theta23, double d);
243 
245 void getAngles(const DoubleMatrix & v, double & t12, double & t13,
246  double & t23, double & d);
249 bool midPtStep(DoubleVector & xi,
250  DoubleVector (*derivs)(double t, const DoubleVector & v),
251  double tInitial, double tStep);
255  DoubleVector (*derivs)(double t,
256  const DoubleVector & v),
257  double tInitial, double tFinal, int numSteps);
258 
260 double fin(double mm1, double mm2);
261 double den(double a, int b);
262 double log1minusx(double x);
265 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
266 double zriddr(double (*func)(double), double x1, double x2, double xacc);
272 void fdjac(int n, DoubleVector x, const DoubleVector & fvec, DoubleMatrix & df,
273  int (*vecfunc)(const DoubleVector &, void*, DoubleVector &),
274  void* params = NULL);
277 bool lnsrch(const DoubleVector & xold, double fold, const DoubleVector & g,
278  DoubleVector & p,
279  DoubleVector & x, double & f, double stpmax,
280  int (*vecfunc)(const DoubleVector &, void*, DoubleVector &),
281  DoubleVector & fvec, void* params = NULL);
282 /* allocate an int vector with subscript range v[nl..nh] */
283 int *ivector(long nl, long nh);
284 /* free an int vector allocated with ivector() */
285 void free_ivector(int *v, long nl, long nh);
286 void lubksb(const DoubleMatrix & a, int n, int *indx, DoubleVector & b);
293 bool newt(DoubleVector & x,
294  int (*vecfunc)(const DoubleVector &, void*, DoubleVector &),
295  void* params = NULL);
297 void load(float x, const DoubleVector & v, DoubleVector & y);
300 void score(float x, const DoubleVector & y, DoubleVector & f);
306 //void shoot(const DoubleVector & v, DoubleVector & f);
313 void qrdcmp(DoubleMatrix & a, int n, DoubleVector & c, DoubleVector & d,
314  int & sing);
315 void qrupdt(DoubleMatrix & r, DoubleMatrix & qt, int n,
316  DoubleVector & u, DoubleVector & v);
320 void rotate(DoubleMatrix & r, DoubleMatrix & qt, int n, int i, float a,
321  float b);
322 void rsolv(const DoubleMatrix & a, int n, const DoubleVector & d,
323  DoubleVector & b);
324 
326  // returns sqrt(f) for f>0
327  double ccbSqrt(double f);
329  // returns sqrt(f) for f>0 or -sqrt(|f|) for f<0
330  double signedSqrt(double f);
332  double signedSqr(double f);
333 
337 double kinFn(double m1, double m2, double m3);
338 double lambda(double a, double b, double c);
339 #endif
340 
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
bool newt(DoubleVector &x, int(*vecfunc)(const DoubleVector &, void *, DoubleVector &), void *params=NULL)
More work can be done on this: get rid of int n and in subfunctions too.
Definition: numerics.cpp:1946
int bin(double data, double start, double end, int numBins)
Definition: numerics.cpp:1073
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:107
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:2255
double b22bar(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:127
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:1374
double zriddr(double(*func)(double), double x1, double x2, double xacc)
Definition: numerics.cpp:1733
void score(float x, const DoubleVector &y, DoubleVector &f)
std::complex< long double > dilogcl(const std::complex< long double > &)
Complex dilogarithm .
Definition: numerics.cpp:1191
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:1525
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:425
bool lnsrch(const DoubleVector &xold, double fold, const DoubleVector &g, DoubleVector &p, DoubleVector &x, double &f, double stpmax, int(*vecfunc)(const DoubleVector &, void *, DoubleVector &), DoubleVector &fvec, void *params=NULL)
Definition: numerics.cpp:1805
Dilog function.
double b0(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:529
void polint(const DoubleVector &xa, const DoubleVector &ya, double x, double &y, double &dy)
Definition: numerics.cpp:1576
double b1(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:668
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:1007
double gfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:118
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:1079
double calcCL(double cl, const DoubleVector &l)
Definition: numerics.cpp:1302
double calc1dFraction(double y, const DoubleVector &l)
Definition: numerics.cpp:1332
double b22(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:770
double signedSqr(double f)
returns f * f * sign(f)
Definition: numerics.cpp:2302
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:149
Complex fBc(const Complex &x)
Definition: numerics.cpp:510
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:1641
bool integrateReversibly(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tFinal, int numSteps)
Definition: numerics.cpp:1685
Complex b22barc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:153
double ccbSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2298
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:1654
double signedSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2300
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
void fdjac(int n, DoubleVector x, const DoubleVector &fvec, DoubleMatrix &df, int(*vecfunc)(const DoubleVector &, void *, DoubleVector &), void *params=NULL)
double gasdev(long &idum)
Definition: numerics.cpp:1015
void shft3(double &a, double &b, double &c, double &d)
a=b, b=c and c=d
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:1067
double sumOfExp(double a, double b)
Definition: numerics.cpp:1090
Complex b1c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:727
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:2203
double midpnt(double(*func)(double), double a, double b, int n)
Mid point function integration.
Definition: numerics.cpp:1411
double d27(double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:921
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:1392
double fin(double mm1, double mm2)
useful for 2-loop mb/mt corrections
Definition: numerics.cpp:1707
double c0(double m1, double m2, double m3)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:933
double fps(double z)
Definition: numerics.cpp:1343
double qromb2(double(*func)(double), double a, double b, double EPS)
Identical copy to facilitate double integral.
Definition: numerics.cpp:1549
double kinFn(double m1, double m2, double m3)
Definition: numerics.cpp:2314
void lubksb(const DoubleMatrix &a, int n, int *indx, DoubleVector &b)
Get rid of int n.
Definition: numerics.cpp:1867
Complex ffnc(double p, double m1, double m2, double q)
Definition: numerics.h:138
DoubleMatrix display3x3RealMixing(double theta12, double theta13, double theta23, double d)
Definition: numerics.cpp:1609
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:144
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:487
double hfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:123
double ran1(long &idum)
Definition: numerics.cpp:1036
double fB(const Complex &x)
Definition: numerics.cpp:496
Complex b22c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:825
double ffn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:112
DoubleVector getRandomDirection(int n, int &numChanged, long &idum)
returns a random direction: total number of dimensions=n
Definition: numerics.cpp:1280
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:886
Complex b0c(double p, double m1, double m2, double q)
Definition: numerics.cpp:598