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 <cmath>
23 #include <cfloat>
24 #include "linalg.h"
25 #include "dilogwrap.h"
26 using namespace softsusy;
27 
32 int cubicRoots(double a, double b, double c, double d,
33  DoubleVector & ans);
34 
36 int cubicRootsInside(double a, double b, double c, DoubleVector & ans);
37 
40 double dgauss(double (*f)(double x), double a, double b, double eps);
41 
44 double dgauss(double (*f)(double x, const DoubleVector & v),
45  const DoubleVector & v, double a, double b, double eps);
46 
48 double accurateSqrt1Plusx(double x);
49 
53 void rungeKuttaStep(const DoubleVector & y, const DoubleVector & dydx,
54  double x, double h, DoubleVector & yout, DoubleVector & yerr,
55  DoubleVector (*derivs)(double, const DoubleVector &));
56 
58 int odeStepper(DoubleVector & y, const DoubleVector & dydx, double *x, double
59  htry, double eps, DoubleVector & yscal, double *hdid,
60  double *hnext,
61  DoubleVector (*derivs)(double, const DoubleVector &));
62 
65 double lnLPoisson(unsigned k, double lambda);
66 
69 double LPoisson(unsigned k, double lambda);
70 
72 int integrateOdes(DoubleVector & ystart, double x1, double x2, double eps,
73  double h1, double hmin,
74  DoubleVector (*derivs)(double, const DoubleVector &),
75  int (*rkqs)
76  (DoubleVector & y, const DoubleVector & dydx, double *x,
77  double htry, double eps, DoubleVector & yscal, double
78  *hdid, double *hnext,
79  DoubleVector (*derivs)(double, const DoubleVector &)));
80 
83 double calcDerivative(double (*func)(double),
84  double x, double h, double *err);
85 
88 double calcDerivative(double (*func)(double, void*),
89  double x, double h, double* err, void* params = NULL);
90 
93 double findMinimum(double ax, double bx, double cx, double (*f)(double),
94  double tol, double & xmin);
95 
96 void shft2(double & a, double & b, double & c);
97 void shft3(double & a, double & b, double & c, double & d);
99 
101 double integrandThreshbn(double x);
103 double bIntegral(int n, double p, double m1, double m2, double mt);
104 DoubleVector dd(double x, const DoubleVector & y);
105 
107 double b0(double p, double m1, double m2, double q);
109 double b1(double p, double m1, double m2, double q);
111 double b22(double p, double m1, double m2, double q);
113 double c0(double m1, double m2, double m3);
115 double d27(double m1, double m2, double m3, double m4);
117 double d0(double m1, double m2, double m3, double m4);
119 inline double a0(double m, double q) {
120  if (fabs(m) < EPSTOL) return 0.;
121  return sqr(m) * (1.0 - 2. * log(abs(m / q)));
122 }
124 inline double ffn(double p, double m1, double m2, double q) {
125  return a0(m1, q) - 2.0 * a0(m2, q) -
126  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
127  b0(p, m1, m2, q);
128 }
130 inline double gfn(double p, double m1, double m2, double q) {
131  return (sqr(p) - sqr(m1) - sqr(m2)) * b0(p, m1, m2, q) - a0(m1, q)
132  - a0(m2, q);
133 }
135 inline double hfn(double p, double m1, double m2, double q) {
136  return 4.0 * b22(p, m1, m2, q) + gfn(p, m1, m2, q);
137 }
139 inline double b22bar(double p, double m1, double m2, double q) {
140  return b22(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
141 }
142 
143 Complex b0c(double p, double m1, double m2, double q);
145 Complex b1c(double p, double m1, double m2, double q);
147 Complex b22c(double p, double m1, double m2, double q);
150 inline Complex ffnc(double p, double m1, double m2, double q) {
151  return a0(m1, q) - 2.0 * a0(m2, q) -
152  (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) *
153  b0c(p, m1, m2, q);
154 }
156 inline Complex gfnc(double p, double m1, double m2, double q) {
157  return (sqr(p) - sqr(m1) - sqr(m2)) * b0c(p, m1, m2, q) - a0(m1, q)
158  - a0(m2, q);
159 }
161 inline Complex hfnc(double p, double m1, double m2, double q) {
162  return 4.0 * b22c(p, m1, m2, q) + gfnc(p, m1, m2, q);
163 }
165 inline Complex b22barc(double p, double m1, double m2, double q) {
166  return b22c(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
167 }
168 
169 /*inline double fB(const Complex & x) {
171  // if (x.real() < EPSTOL * 10. && x.imag() < EPSTOL * 10.) return -1.;
172  return (log(Complex(1.0 - x)) - x * log(Complex(1.0 - 1.0 / x))
173  - Complex(1.0)).real();
174  } */
175 double fB(const Complex & x);
176 
177 Complex fBc(const Complex & x);
178 double dilogarg(double t);
179 double dilog(double x);
180 
181 double integrandThreshbnr(double x);
182 Complex fnfn(double x);
183 
186 double gasdev(long & idum);
189 double truncGaussWidthHalf(long & idum);
192 double ran1(long & idum);
195 double cauchyRan(long & idum);
196 
200 int bin(double data, double start, double end, int numBins);
201 
203 double logOfSum(double a, double b);
204 
206 DoubleVector getRandomDirection(int n, int & numChanged, long & idum);
207 
211 double calcCL(double cl, const DoubleVector & l);
212 
216 double calc1dFraction(double y, const DoubleVector & l);
217 
220 double fps(double z);
221 double fs(double z);
222 double ffbar(double z);
223 
224 Complex dilog(const Complex & x);
225 std::complex<long double> dilogcl(const std::complex<long double>&);
226 
228 double trapzd(double (*func)(double), double a, double b, int n,
229  double tol = 1.0e-3);
231 double qtrap(double (*func)(double), double a, double b, double tol);
233 double midpnt(double (*func)(double), double a, double b, int n);
236 double qromb(double (*func)(double), double a, double b, double EPS);
238 double qromb2(double (*func)(double), double a, double b, double EPS);
239 
242 double sumOfExp(double a, double b);
243 
246 void polint(const DoubleVector& xa, const DoubleVector & ya, double x,
247  double & y, double & dy);
248 
253 DoubleMatrix display3x3RealMixing(double theta12, double theta13,
254  double theta23, double d);
255 
257 void getAngles(const DoubleMatrix & v, double & t12, double & t13,
258  double & t23, double & d);
261 bool midPtStep(DoubleVector & xi,
262  DoubleVector (*derivs)(double t, const DoubleVector & v),
263  double tInitial, double tStep);
267  DoubleVector (*derivs)(double t,
268  const DoubleVector & v),
269  double tInitial, double tFinal, int numSteps);
270 
272 double fin(double mm1, double mm2);
273 double den(double a, int b);
274 double log1minusx(double x);
277 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
278 double zriddr(double (*func)(double), double x1, double x2, double xacc);
284 void fdjac(int n, DoubleVector x, const DoubleVector & fvec, DoubleMatrix & df,
285  int (*vecfunc)(const DoubleVector &, void*, DoubleVector &),
286  void* params = NULL);
289 bool lnsrch(const DoubleVector & xold, double fold, const DoubleVector & g,
290  DoubleVector & p,
291  DoubleVector & x, double & f, double stpmax,
292  int (*vecfunc)(const DoubleVector &, void*, DoubleVector &),
293  DoubleVector & fvec, void* params = NULL);
294 /* allocate an int vector with subscript range v[nl..nh] */
295 int *ivector(long nl, long nh);
296 /* free an int vector allocated with ivector() */
297 void free_ivector(int *v, long nl, long nh);
298 void lubksb(const DoubleMatrix & a, int n, int *indx, DoubleVector & b);
305 bool newt(DoubleVector & x,
306  int (*vecfunc)(const DoubleVector &, void*, DoubleVector &),
307  void* params = NULL);
309 void load(float x, const DoubleVector & v, DoubleVector & y);
312 void score(float x, const DoubleVector & y, DoubleVector & f);
318 //void shoot(const DoubleVector & v, DoubleVector & f);
325 void qrdcmp(DoubleMatrix & a, int n, DoubleVector & c, DoubleVector & d,
326  int & sing);
327 void qrupdt(DoubleMatrix & r, DoubleMatrix & qt, int n,
328  DoubleVector & u, DoubleVector & v);
332 void rotate(DoubleMatrix & r, DoubleMatrix & qt, int n, int i, float a,
333  float b);
334 void rsolv(const DoubleMatrix & a, int n, const DoubleVector & d,
335  DoubleVector & b);
336 
338  // returns sqrt(f) for f>0
339  double ccbSqrt(double f);
341  // returns sqrt(f) for f>0 or -sqrt(|f|) for f<0
342  double signedSqrt(double f);
344  double signedSqr(double f);
345 
349 double kinFn(double m1, double m2, double m3);
350 double lambda(double a, double b, double c);
351 #endif
352 
double b22(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:860
bool integrateReversibly(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tFinal, int numSteps)
Definition: numerics.cpp:1781
const double EPSTOL
underflow accuracy
Definition: def.h:37
global variable declaration
Definition: def.cpp:13
switches (options) and parameters such as default fermion masses, required accuracy etc
double calcDerivative(double(*func)(double), double x, double h, double *err)
Definition: numerics.cpp:412
void qrdcmp(DoubleMatrix &a, int n, DoubleVector &c, DoubleVector &d, int &sing)
Definition: numerics.cpp:2299
double fps(double z)
Definition: numerics.cpp:1439
double ran1(long &idum)
Definition: numerics.cpp:1132
std::complex< long double > dilogcl(const std::complex< long double > &)
Complex dilogarithm .
Definition: numerics.cpp:1287
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:252
void lubksb(const DoubleMatrix &a, int n, int *indx, DoubleVector &b)
Get rid of int n.
Definition: numerics.cpp:1963
double accurateSqrt1Plusx(double x)
calculate root(1+x), where x<<1 accurately
Definition: numerics.cpp:209
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:327
double b22bar(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:139
DoubleMatrix display3x3RealMixing(double theta12, double theta13, double theta23, double d)
Definition: numerics.cpp:1705
double a0(double m, double q)
inlined PV functions: real part
Definition: numerics.h:119
double sqr(double a)
square of a number
Definition: utils.h:49
drop-in replacement for the original home-grown Complex class
Definition: mycomplex.h:17
DoubleVector is of variable length, and contains double precision.
Definition: linalg.h:35
int cubicRootsInside(double a, double b, double c, DoubleVector &ans)
cubic roots for x^3 + a*x + b*x^2 + c*x^4 = 0
Definition: numerics.cpp:53
double fin(double mm1, double mm2)
useful for 2-loop mb/mt corrections
Definition: numerics.cpp:1803
double qromb2(double(*func)(double), double a, double b, double EPS)
Identical copy to facilitate double integral.
Definition: numerics.cpp:1645
double logOfSum(double a, double b)
Adds logs of two numbers in a more careful way that avoids underflow.
Definition: numerics.cpp:1175
double qromb(double(*func)(double), double a, double b, double EPS)
Definition: numerics.cpp:1621
void load(float x, const DoubleVector &v, DoubleVector &y)
calculates the n-vector y, given freely specifiable values v(1..n2) at x1
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:2042
double gasdev(long &idum)
Definition: numerics.cpp:1111
double log1minusx(double x)
This function DOESNT WORK YET!!!
Definition: numerics.cpp:230
double ffn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:124
double integrandThreshbn(double x)
For calculation of PV functions.
Complex b22c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:918
Complex b22barc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:165
Complex ffnc(double p, double m1, double m2, double q)
Definition: numerics.h:150
A few handy bits and pieces - little mathematical functions and the like.
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:1901
void fdjac(int n, DoubleVector x, const DoubleVector &fvec, DoubleMatrix &df, int(*vecfunc)(const DoubleVector &, void *, DoubleVector &), void *params=NULL)
Complex b0c(double p, double m1, double m2, double q)
Definition: numerics.cpp:673
void shft3(double &a, double &b, double &c, double &d)
a=b, b=c and c=d
Matrix from 1..rows, 1..cols of double values.
Definition: linalg.h:214
Complex b1c(double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer.
Definition: numerics.cpp:815
bool midPtStep(DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tStep)
Definition: numerics.cpp:1750
void rungeKuttaStep(const DoubleVector &y, const DoubleVector &dydx, double x, double h, DoubleVector &yout, DoubleVector &yerr, DoubleVector(*derivs)(double, const DoubleVector &))
Definition: numerics.cpp:365
double signedSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2396
double calcCL(double cl, const DoubleVector &l)
Definition: numerics.cpp:1398
DoubleVector getRandomDirection(int n, int &numChanged, long &idum)
returns a random direction: total number of dimensions=n
Definition: numerics.cpp:1376
double gfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:130
double d27(double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:1017
double d0(double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:982
void rotate(DoubleMatrix &r, DoubleMatrix &qt, int n, int i, float a, float b)
Definition: numerics.cpp:2351
double sumOfExp(double a, double b)
Definition: numerics.cpp:1186
double dgauss(double(*f)(double x), double a, double b, double eps)
Definition: numerics.cpp:83
double midpnt(double(*func)(double), double a, double b, int n)
Mid point function integration.
Definition: numerics.cpp:1507
DoubleVector and DoubleMatrix classes of doubles and operations between them, complexified copies als...
complex numbers and operators between them
double fB(const Complex &x)
Definition: numerics.cpp:559
double zriddr(double(*func)(double), double x1, double x2, double xacc)
Definition: numerics.cpp:1829
double hfn(double p, double m1, double m2, double q)
inlined PV functions: real part
Definition: numerics.h:135
double ccbSqrt(double f)
returns the square root of the absolute value of the argument
Definition: numerics.cpp:2394
double b1(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:754
double qtrap(double(*func)(double), double a, double b, double tol)
Driver for integration.
Definition: numerics.cpp:1488
double b0(double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:592
double signedSqr(double f)
returns f * f * sign(f)
Definition: numerics.cpp:2398
double findMinimum(double ax, double bx, double cx, double(*f)(double), double tol, double &xmin)
Definition: numerics.cpp:488
double c0(double m1, double m2, double m3)
Passarino-Veltman function definition: real part.
Definition: numerics.cpp:1029
double kinFn(double m1, double m2, double m3)
Definition: numerics.cpp:2410
Dilog function.
double calc1dFraction(double y, const DoubleVector &l)
Definition: numerics.cpp:1428
double cauchyRan(long &idum)
Definition: numerics.cpp:1163
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:550
Complex hfnc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:161
void shft2(double &a, double &b, double &c)
void polint(const DoubleVector &xa, const DoubleVector &ya, double x, double &y, double &dy)
Definition: numerics.cpp:1672
double lnLPoisson(unsigned k, double lambda)
Definition: numerics.cpp:311
Complex fBc(const Complex &x)
Definition: numerics.cpp:573
void score(float x, const DoubleVector &y, DoubleVector &f)
Complex gfnc(double p, double m1, double m2, double q)
inlined PV functions: complex answer
Definition: numerics.h:156
int cubicRoots(double a, double b, double c, double d, DoubleVector &ans)
Definition: numerics.cpp:20
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:1470
int bin(double data, double start, double end, int numBins)
Definition: numerics.cpp:1169
double truncGaussWidthHalf(long &idum)
Definition: numerics.cpp:1103
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:1737
double LPoisson(unsigned k, double lambda)
Definition: numerics.cpp:319
double m1
decay global variable declarations
Definition: decays.cpp:14