SoftSUSY is hosted by Hepforge, IPPP Durham
SOFTSUSY  4.0
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 <iostream>
21 #include "def.h"
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 
34 double accurateSqrt1Plusx(double x);
35 
39 void rungeKuttaStep(const DoubleVector & y, const DoubleVector & dydx,
40  double x, double h, DoubleVector & yout, DoubleVector & yerr,
41  DoubleVector (*derivs)(double, const DoubleVector &));
42 
44 int odeStepper(DoubleVector & y, const DoubleVector & dydx, double *x, double
45  htry, double eps, DoubleVector & yscal, double *hdid,
46  double *hnext,
47  DoubleVector (*derivs)(double, const DoubleVector &));
48 
51 double lnLPoisson(unsigned k, double lambda);
52 
55 double LPoisson(unsigned k, double lambda);
56 
58 int integrateOdes(DoubleVector & ystart, double x1, double x2, double eps,
59  double h1, double hmin,
60  DoubleVector (*derivs)(double, const DoubleVector &),
61  int (*rkqs)
62  (DoubleVector & y, const DoubleVector & dydx, double *x,
63  double htry, double eps, DoubleVector & yscal, double
64  *hdid, double *hnext,
65  DoubleVector (*derivs)(double, const DoubleVector &)));
66 
69 double calcDerivative(double (*func)(double),
70  double x, double h, double *err);
71 
74 double findMinimum(double ax, double bx, double cx, double (*f)(double),
75  double tol, double *xmin);
76 
77 void shft2(double & a, double & b, double & c);
78 void shft3(double & a, double & b, double & c, double & d);
80 
82 double integrandThreshbn(double x);
84 double bIntegral(int n, double p, double m1, double m2, double mt);
85 DoubleVector dd(double x, const DoubleVector & y);
86 
88 double b0(double p, double m1, double m2, double q);
90 double b1(double p, double m1, double m2, double q);
92 double b22(double p, double m1, double m2, double q);
94 double c0(double m1, double m2, double m3);
96 double d27(double m1, double m2, double m3, double m4);
98 double d0(double m1, double m2, double m3, double m4);
100 inline double a0(double m, double q) {
101  if (fabs(m) < EPSTOL) return 0.;
102  return sqr(m) * (1.0 - 2. * log(abs(m / q)));
103 }
105 inline double ffn(double p, double m1, double m2, double q) {
106  return a0(m1, q) - 2.0 * a0(m2, q) -
107  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
108  b0(p, m1, m2, q);
109 }
111 inline double gfn(double p, double m1, double m2, double q) {
112  return (sqr(p) - sqr(m1) - sqr(m2)) * b0(p, m1, m2, q) - a0(m1, q)
113  - a0(m2, q);
114 }
116 inline double hfn(double p, double m1, double m2, double q) {
117  return 4.0 * b22(p, m1, m2, q) + gfn(p, m1, m2, q);
118 }
120 inline double b22bar(double p, double m1, double m2, double q) {
121  return b22(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
122 }
123 
124 Complex b0c(double p, double m1, double m2, double q);
126 Complex b1c(double p, double m1, double m2, double q);
128 Complex b22c(double p, double m1, double m2, double q);
131 inline Complex ffnc(double p, double m1, double m2, double q) {
132  return a0(m1, q) - 2.0 * a0(m2, q) -
133  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
134  b0c(p, m1, m2, q);
135 }
137 inline Complex gfnc(double p, double m1, double m2, double q) {
138  return (sqr(p) - sqr(m1) - sqr(m2)) * b0c(p, m1, m2, q) - a0(m1, q)
139  - a0(m2, q);
140 }
142 inline Complex hfnc(double p, double m1, double m2, double q) {
143  return 4.0 * b22c(p, m1, m2, q) + gfnc(p, m1, m2, q);
144 }
146 inline Complex b22barc(double p, double m1, double m2, double q) {
147  return b22c(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
148 }
149 
150 /*inline double fB(const Complex & x) {
152  // if (x.real() < EPSTOL * 10. && x.imag() < EPSTOL * 10.) return -1.;
153  return (log(Complex(1.0 - x)) - x * log(Complex(1.0 - 1.0 / x))
154  - Complex(1.0)).real();
155  } */
156 double fB(const Complex & x);
157 
158 Complex fBc(const Complex & x);
159 double dilogarg(double t);
160 double dilog(double x);
161 
162 double integrandThreshbnr(double x);
163 Complex fnfn(double x);
164 
167 double gasdev(long & idum);
170 double truncGaussWidthHalf(long & idum);
173 double ran1(long & idum);
176 double cauchyRan(long & idum);
177 
181 int bin(double data, double start, double end, int numBins);
182 
184 double logOfSum(double a, double b);
185 
187 DoubleVector getRandomDirection(int n, int & numChanged, long & idum);
188 
192 double calcCL(double cl, const DoubleVector & l);
193 
197 double calc1dFraction(double y, const DoubleVector & l);
198 
201 double fps(double z);
202 double fs(double z);
203 double ffbar(double z);
204 
205 Complex dilog(const Complex & x);
206 
208 double trapzd(double (*func)(double), double a, double b, int n,
209  double tol = 1.0e-3);
211 double qtrap(double (*func)(double), double a, double b, double tol);
213 double midpnt(double (*func)(double), double a, double b, int n);
216 double qromb(double (*func)(double), double a, double b, double EPS);
218 double qromb2(double (*func)(double), double a, double b, double EPS);
219 
222 double sumOfExp(double a, double b);
223 
226 void polint(const DoubleVector& xa, const DoubleVector & ya, double x,
227  double & y, double & dy);
228 
233 DoubleMatrix display3x3RealMixing(double theta12, double theta13,
234  double theta23, double d);
235 
237 void getAngles(const DoubleMatrix & v, double & t12, double & t13,
238  double & t23, double & d);
241 bool midPtStep(DoubleVector & xi,
242  DoubleVector (*derivs)(double t, const DoubleVector & v),
243  double tInitial, double tStep);
247  DoubleVector (*derivs)(double t,
248  const DoubleVector & v),
249  double tInitial, double tFinal, int numSteps);
250 
252 double fin(double mm1, double mm2);
253 double den(double a, int b);
254 double log1minusx(double x);
257 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
258 double zriddr(double (*func)(double), double x1, double x2, double xacc);
264 void fdjac(int n, DoubleVector x, const DoubleVector & fvec, DoubleMatrix & df,
265  void (*vecfunc)(const DoubleVector &, DoubleVector &));
268 bool lnsrch(const DoubleVector & xold, double fold, const DoubleVector & g,
269  DoubleVector & p,
270  DoubleVector & x, double & f, double stpmax,
271  void (*vecfunc)(const DoubleVector &, DoubleVector &),
272  DoubleVector & fvec);
273 /* allocate an int vector with subscript range v[nl..nh] */
274 int *ivector(long nl, long nh);
275 /* free an int vector allocated with ivector() */
276 void free_ivector(int *v, long nl, long nh);
277 void lubksb(const DoubleMatrix & a, int n, int *indx, DoubleVector & b);
284 bool newt(DoubleVector & x,
285  void (*vecfunc)(const DoubleVector &, DoubleVector &));
287 void load(float x, const DoubleVector & v, DoubleVector & y);
290 void score(float x, const DoubleVector & y, DoubleVector & f);
296 //void shoot(const DoubleVector & v, DoubleVector & f);
303 void qrdcmp(DoubleMatrix & a, int n, DoubleVector & c, DoubleVector & d,
304  int & sing);
305 void qrupdt(DoubleMatrix & r, DoubleMatrix & qt, int n,
306  DoubleVector & u, DoubleVector & v);
310 void rotate(DoubleMatrix & r, DoubleMatrix & qt, int n, int i, float a,
311  float b);
312 void rsolv(const DoubleMatrix & a, int n, const DoubleVector & d,
313  DoubleVector & b);
314 
316  // returns sqrt(f) for f>0
317  double ccbSqrt(double f);
319  // returns sqrt(f) for f>0 or -sqrt(|f|) for f<0
320  double signedSqrt(double f);
322  double signedSqr(double f);
323 
324 #endif
325 
void rungeKuttaStep(const DoubleVector &y, const DoubleVector &dydx, double x, double h, DoubleVector &yout, DoubleVector &yerr, DoubleVector(*derivs)(double, const DoubleVector &))
Definition: numerics.cpp:230
global variable declaration
Definition: def.cpp:13
int bin(double data, double start, double end, int numBins)
Definition: numerics.cpp:966
Complex log(const Complex &a)
Principal log.
Definition: mycomplex.h:67
double log1minusx(double x)
This function DOESNT WORK YET!!!
Definition: numerics.cpp:95
double lnLPoisson(unsigned k, double lambda)
Definition: numerics.cpp:176
double a0(double m, double q)
inlined PV functions: real part
Definition: numerics.h:100
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:2046
double b22bar(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:120
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:1168
double zriddr(double(*func)(double), double x1, double x2, double xacc)
Definition: numerics.cpp:1527
void score(float x, const DoubleVector &y, DoubleVector &f)
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:1319
DoubleVector and DoubleMatrix classes of doubles and operations between them, complexified copies als...
Dilog function.
double b0(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:422
void polint(const DoubleVector &xa, const DoubleVector &ya, double x, double &y, double &dy)
Definition: numerics.cpp:1370
double b1(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:561
double accurateSqrt1Plusx(double x)
calculate root(1+x), where x<<1 accurately
Definition: numerics.cpp:74
void shft2(double &a, double &b, double &c)
double truncGaussWidthHalf(long &idum)
Definition: numerics.cpp:900
double gfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:111
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:117
double logOfSum(double a, double b)
Adds logs of two numbers in a more careful way that avoids underflow.
Definition: numerics.cpp:972
double calcCL(double cl, const DoubleVector &l)
Definition: numerics.cpp:1096
double calc1dFraction(double y, const DoubleVector &l)
Definition: numerics.cpp:1126
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:1598
double b22(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:663
double signedSqr(double f)
returns f * f * sign(f)
Definition: numerics.cpp:2093
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:142
Complex fBc(const Complex &x)
Definition: numerics.cpp:403
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:1435
bool integrateReversibly(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tFinal, int numSteps)
Definition: numerics.cpp:1479
Complex b22barc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:146
double ccbSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2089
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:1448
double signedSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2091
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:192
double gasdev(long &idum)
Definition: numerics.cpp:908
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:1739
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:277
double cauchyRan(long &idum)
Definition: numerics.cpp:960
double sumOfExp(double a, double b)
Definition: numerics.cpp:983
Complex b1c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:620
double dgauss(double(*f)(double x), double a, double b, double eps)
Definition: numerics.cpp:13
void qrdcmp(DoubleMatrix &a, int n, DoubleVector &c, DoubleVector &d, int &sing)
Definition: numerics.cpp:1994
double midpnt(double(*func)(double), double a, double b, int n)
Mid point function integration.
Definition: numerics.cpp:1205
double d27(double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:814
complex numbers and operators between them
double qtrap(double(*func)(double), double a, double b, double tol)
Driver for integration.
Definition: numerics.cpp:1186
double findMinimum(double ax, double bx, double cx, double(*f)(double), double tol, double *xmin)
Definition: numerics.cpp:318
double fin(double mm1, double mm2)
useful for 2-loop mb/mt corrections
Definition: numerics.cpp:1501
double c0(double m1, double m2, double m3)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:826
double fps(double z)
Definition: numerics.cpp:1137
double qromb2(double(*func)(double), double a, double b, double EPS)
Identical copy to facilitate double integral.
Definition: numerics.cpp:1343
void lubksb(const DoubleMatrix &a, int n, int *indx, DoubleVector &b)
Get rid of int n.
Definition: numerics.cpp:1660
Complex ffnc(double p, double m1, double m2, double q)
Definition: numerics.h:131
DoubleMatrix display3x3RealMixing(double theta12, double theta13, double theta23, double d)
Definition: numerics.cpp:1403
const double EPSTOL
underflow accuracy
Definition: def.h:60
Complex gfnc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:137
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:380
double hfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:116
double ran1(long &idum)
Definition: numerics.cpp:929
double fB(const Complex &x)
Definition: numerics.cpp:389
Complex b22c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:718
double ffn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:105
DoubleVector getRandomDirection(int n, int &numChanged, long &idum)
returns a random direction: total number of dimensions=n
Definition: numerics.cpp:1074
double LPoisson(unsigned k, double lambda)
Definition: numerics.cpp:184
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:779
Complex b0c(double p, double m1, double m2, double q)
Definition: numerics.cpp:491