SoftSUSY is hosted by Hepforge, IPPP Durham
SOFTSUSY  4.1
numerics.cpp File Reference
#include "numerics.h"

Macros

#define JMAX   20
 
#define FUNC(x)   ((*func)(x))
 
#define JMAX   20
 
#define JMAXP   (JMAX+1)
 
#define K   5
 

Functions

double dgauss (double(*f)(double x), double a, double b, double eps)
 
double dgauss (double(*f)(double x, const DoubleVector &v), const DoubleVector &v, double a, double b, double eps)
 
double accurateSqrt1Plusx (double x)
 calculate root(1+x), where x<<1 accurately More...
 
double log1minusx (double x)
 This function DOESNT WORK YET!!! More...
 
int integrateOdes (DoubleVector &ystart, double from, double to, 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.
 
double lnLPoisson (unsigned k, double lambda)
 
double LPoisson (unsigned k, double lambda)
 
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
 
void rungeKuttaStep (const DoubleVector &y, const DoubleVector &dydx, double x, double h, DoubleVector &yout, DoubleVector &yerr, DoubleVector(*derivs)(double, const DoubleVector &))
 
double calcDerivative (double(*func)(double), double x, double h, double *err)
 
void shft2 (double &a, double &b, double c)
 
void shft3 (double &a, double &b, double &c, double d)
 
double findMinimum (double ax, double bx, double cx, double(*f)(double), double tol, double &xmin)
 
DoubleVector dd (double x, const DoubleVector &)
 
double integrandThreshbnr (double x)
 
Complex fnfn (double x)
 
double bIntegral (int n1, double p, double m1, double m2, double mt)
 Returns real part of b function, less accurate than analytic expressions.
 
double fB (const Complex &a)
 
Complex fBc (const Complex &a)
 
double b0 (double p, double m1, double m2, double q)
 Passarino-Veltman function definition: real part. More...
 
Complex b0c (double p, double m1, double m2, double q)
 
double b1 (double p, double m1, double m2, double q)
 Note that b1 is NOT symmetric in m1 <-> m2!!! More...
 
Complex b1c (double p, double m1, double m2, double q)
 Note that b1 is NOT symmetric in m1 <-> m2!!! More...
 
double b22 (double p, double m1, double m2, double q)
 Passarino-Veltman function definition: real part. More...
 
Complex b22c (double p, double m1, double m2, double q)
 Passarino-Veltman function definition: complex answer. More...
 
double d0 (double m1, double m2, double m3, double m4)
 Passarino-Veltman function definition: real part.
 
double d27 (double m1, double m2, double m3, double m4)
 Passarino-Veltman function definition: real part.
 
double c0 (double m1, double m2, double m3)
 Passarino-Veltman function definition: real part.
 
double truncGaussWidthHalf (long &idum)
 
double gasdev (long &idum)
 
double ran1 (long &idum)
 
double cauchyRan (long &idum)
 
int bin (double data, double start, double end, int numBins)
 
double logOfSum (double a, double b)
 Adds logs of two numbers in a more careful way that avoids underflow.
 
double sumOfExp (double a, double b)
 
double dilog (double x)
 
Complex dilog (const Complex &x)
 
std::complex< long double > dilogcl (const std::complex< long double > &z)
 Complex dilogarithm $\mathrm{Li}_2(z)$. More...
 
DoubleVector getRandomDirection (int n, int &numChanged, long &idum)
 returns a random direction: total number of dimensions=n
 
double calcCL (double cl, const DoubleVector &l)
 
int * ivector (long nl, long nh)
 
void free_ivector (int *v, long nl, long)
 
double calc1dFraction (double y, const DoubleVector &l)
 
double fps (double z)
 
double fs (double z)
 
double ffbar (double z)
 
double trapzd (double(*func)(double), double a, double b, int n, double)
 Trapezoidal function integration to accuracy tol.
 
double qtrap (double(*func)(double), double a, double b, double EPS)
 Driver for integration.
 
double midpnt (double(*func)(double), double a, double b, int n)
 Mid point function integration.
 
double edgefn (double x, double y, double z)
 
double mllMax (double mChi1, double mSlep, double mChi2)
 
double mllq (double mSq, double mChi1, double mSlep, double mChi2)
 
double llqThresh (double mSq, double mChi1, double mSlep, double mChi2)
 
double lqnear (double mSq, double, double mSlep, double mChi2)
 
double lqfar (double mSq, double mChi1, double mSlep, double mChi2)
 
double lqhigh (double mSq, double mChi1, double mSlep, double mChi2)
 
double lqlow (double mSq, double mChi1, double mSlep, double mChi2)
 
double qromb (double(*func)(double), double a, double b, double EPS)
 
double qromb2 (double(*func)(double), double a, double b, double EPS)
 Identical copy to facilitate double integral.
 
void polint (const DoubleVector &xa, const DoubleVector &ya, double x, double &y, double &dy)
 
DoubleMatrix display3x3RealMixing (double theta12, double theta13, double theta23, double d)
 
void getAngles (const DoubleMatrix &v, double &t12, double &t13, double &t23, double &d)
 Given a matrix v, determines which angles gave it.
 
bool midPtStep (DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tStep)
 
bool integrateReversibly (DoubleVector &xi, DoubleVector(*derivs)(double t, const DoubleVector &v), double tInitial, double tFinal, int numSteps)
 
double den (double a, int b)
 
double fin (double mm1, double mm2)
 useful for 2-loop mb/mt corrections
 
double zriddr (double(*func)(double), double x1, double x2, double xacc)
 
DoubleMatrix fdjac (int n, DoubleVector x, const DoubleVector &fvec, void(*vecfunc)(const DoubleVector &, DoubleVector &))
 You will need to clear this lot up....
 
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)
 
void lubksb (const DoubleMatrix &a, int n, int *indx, DoubleVector &b)
 Get rid of int n.
 
void ludcmp (DoubleMatrix &a, int n, int *indx, double &d)
 Get rid of int n.
 
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. More...
 
DoubleVector testDerivs (double, const DoubleVector &y)
 
void broydn (DoubleVector x, int &check, void(*vecfunc)(const DoubleVector &, DoubleVector &))
 
void qrdcmp (DoubleMatrix &a, int n, DoubleVector &c, DoubleVector &d, int &sing)
 
void qrupdt (DoubleMatrix &r, DoubleMatrix &qt, int n, DoubleVector &u, DoubleVector &v)
 
void rotate (DoubleMatrix &r, DoubleMatrix &qt, int n, int i, float a, float b)
 
void rsolv (const DoubleMatrix &a, int n, const DoubleVector &d, DoubleVector &b)
 
double ccbSqrt (double f)
 returns the square root of the absolute value of the argument
 
double signedSqrt (double f)
 returns the square root of the absolute value of the argument
 
double signedSqr (double f)
 returns f * f * sign(f)
 
double kinFn (double m1, double m2, double m3)
 
double lambda (double a, double b, double c)
 

Detailed Description

Function Documentation

double accurateSqrt1Plusx ( double  x)

calculate root(1+x), where x<<1 accurately

calculate approximate number of terms we're going to need to get 16 digit precision

double b0 ( double  p,
double  m1,
double  m2,
double  q 
)

Passarino-Veltman function definition: real part.

Avoids IR infinities

Try to increase the accuracy of s

Decides level at which one switches to p=0 limit of calculations

p is not 0

Complex b0c ( double  p,
double  m1,
double  m2,
double  q 
)

Avoids IR infinities

Try to increase the accuracy of s

Decides level at which one switches to p=0 limit of calculations

p is not 0

double b1 ( double  p,
double  m1,
double  m2,
double  q 
)

Note that b1 is NOT symmetric in m1 <-> m2!!!

Passarino-Veltman function definition: real part.

Decides level at which one switches to p=0 limit of calculations

< checked

Complex b1c ( double  p,
double  m1,
double  m2,
double  q 
)

Note that b1 is NOT symmetric in m1 <-> m2!!!

Passarino-Veltman function definition: complex answer.

Decides level at which one switches to p=0 limit of calculations

< checked

< checked

double b22 ( double  p,
double  m1,
double  m2,
double  q 
)

Passarino-Veltman function definition: real part.

Decides level at which one switches to p=0 limit of calculations

This zero p limit is good

Complex b22c ( double  p,
double  m1,
double  m2,
double  q 
)

Passarino-Veltman function definition: complex answer.

Decides level at which one switches to p=0 limit of calculations

This zero p limit is good

int bin ( double  data,
double  start,
double  end,
int  numBins 
)

Returns the number of a bin that the data is in: from 1 to numBins in the range (bins other than this range are also possible - you must deal with them outside the function...)

double calc1dFraction ( double  y,
const DoubleVector l 
)

given a normalised binned likelihood vector l, calculates the fraction of bins with likelihood less than or equal to y as a % of maximum: approximates area OUTSIDE confidence level y

double calcCL ( double  cl,
const DoubleVector l 
)

Calculates the vertical level of confidence level required to contain a fraction cl of area of the histogrammed binned likelihood l. If err is true, a satisfactory answer could not be found for some reason.

double calcDerivative ( double(*)(double)  func,
double  x,
double  h,
double *  err 
)

func is user-supplied, h is an estimate of what step-size to start with and err returns error flags

double cauchyRan ( long &  idum)

Cauchy distribution ie 1 / [ pi gamma (1 + x^2/gamma^2) ]. For a width, you must multiply the x coming out by the width.

double dgauss ( double(*)(double x)  f,
double  a,
double  b,
double  eps 
)

adaptive Gaussian one dimensional integration of f(x) between a and b to precision eps

set the initial values if they are uninitialised

double dgauss ( double(*)(double x, const DoubleVector &v)  f,
const DoubleVector v,
double  a,
double  b,
double  eps 
)

adaptive Gaussian one dimensional integration of f(x) between a and b to precision eps with user - given parameters in DoubleVector & v

set the initial values if they are uninitialised

std::complex<long double> dilogcl ( const std::complex< long double > &  z)

Complex dilogarithm $\mathrm{Li}_2(z)$.

Parameters
zcomplex argument
Note
Implementation translated from SPheno to C++
Returns
$\mathrm{Li}_2(z)$
DoubleMatrix display3x3RealMixing ( double  theta12,
double  theta13,
double  theta23,
double  d 
)

Returns a 3 by 3 real mixing matrix. Input angles are standard CKM parameterisation. If the phase d is not zero, the result is only an approximation to the full complex matrix: see SOFTSUSY manual for details.

phase factor e^i delta: we'll set it to + or - 1 depending on the sign of s13

double fB ( const Complex a)

First, special cases at problematic points

Complex fBc ( const Complex a)

First, special cases at problematic points

double findMinimum ( double  ax,
double  bx,
double  cx,
double(*)(double)  f,
double  tol,
double &  xmin 
)

f is user-defined function, minimum value returned in xmin. Based on a golden section search

double fps ( double  z)

These three functions are for the calculation of 2-loop log pieces of g-2 of the muon

answer should always be real

double gasdev ( long &  idum)

Gaussian deviated random number, mean 0 variance 1. Don't re-set idum once you've initially set it. Initialise with a NEGATIVE integer

bool integrateReversibly ( DoubleVector xi,
DoubleVector(*)(double t, const DoubleVector &v)  derivs,
double  tInitial,
double  tFinal,
int  numSteps 
)

Do a fixed number of steps of approximately reversible integration: hoping it will help convergence in difficult cases. It is really slow though.

double kinFn ( double  m1,
double  m2,
double  m3 
)

Kinematic mass function - differs from the pure one in PDG by a root and factor $ f(m_1,m_2,m_3)=\frac{\sqrt{(m_1^2-(m_2^2+m_3^2))(m_1^2-(m_2^2-m_3^2))}}{2 m_1}$.

double lnLPoisson ( unsigned  k,
double  lambda 
)

Calculates log likelihood of a Poisson with k observed events, expecting lambda>0.

bool lnsrch ( const DoubleVector xold,
double  fold,
const DoubleVector g,
DoubleVector p,
DoubleVector x,
double &  f,
double  stpmax,
void(*)(const DoubleVector &, DoubleVector &)  vecfunc,
DoubleVector fvec 
)

These are experimental things for trying the shooting method - returns F.F/2 evaluated at x. Boolean value on return is error flag

double log1minusx ( double  x)

This function DOESNT WORK YET!!!

1/a^b

Find largest power that we need from the expansion

double LPoisson ( unsigned  k,
double  lambda 
)

Calculates likelihood of a Poisson with k observed events, expecting lambda>0.

bool midPtStep ( DoubleVector xi,
DoubleVector(*)(double t, const DoubleVector &v)  derivs,
double  tInitial,
double  tStep 
)

Evolves the dependent variables xi by one reversible mid-point step of length tStep. Returns true if there is an error.

initial guess

difference between iterations

bool newt ( DoubleVector x,
void(*)(const DoubleVector &, DoubleVector &)  vecfunc 
)

More work can be done on this: get rid of int n and in subfunctions too.

Multi-dimensional globally convergent multi-dimensional root solver for n variables. adjusts x so that f(x)=0, where f is the last function provided by vecfunc, given vector input x. If false on output, there is no error. If returns 1, a local minimum or saddle-point has been found (df/dx=0). The length of x should be equal to the number of parameters to vary AND the number of constraints ie the length of the vector in vecfunc.

< max iterations

< convergence on function values

< spurious convergence to min of fmin

< maximum dx convergence criterion

< maximum step length allowed in line searches

Check points aren't getting too close together

void polint ( const DoubleVector xa,
const DoubleVector ya,
double  x,
double &  y,
double &  dy 
)

Returns a value of the polynomial (y) and an error estimate (dy), given a bunch of points (xa) and their y values (ya).

void qrdcmp ( DoubleMatrix a,
int  n,
DoubleVector c,
DoubleVector d,
int &  sing 
)

Function that integrates ODEs: to be passed to Newton's method for solving the 2-boundary value problem. Given v[1..n2], it sets all the boundary conditions at the initial scale and calculates f[1..n2] - a score for how well the boundary condition at the high scale is satisfied. Then you should be able to find the solutions for the unknown numbers. QR decomposition of the matrix a. A = Q.R. Upper triangular matrix R is returned in upper triangle of a, except for diagonal elements which are returned in d. Orthog matrix Q is given as a product of n-1 Houselholder matrices $Q_1 Q_{n-1}$ where $Q_j=1-u_j u_j/c_j$. sing returns true (1) if singularity is encountered (but decomposition is still completed) otherwise false (0).

double qromb ( double(*)(double)  func,
double  a,
double  b,
double  EPS 
)

Integrate a function *func which is analytic to precision EPS between a and b

double ran1 ( long &  idum)

Normally distributed random number between 0 and 1. Don't re-set idum once you've initially set it. Initialise with a NEGATIVE integer

void rotate ( DoubleMatrix r,
DoubleMatrix qt,
int  n,
int  i,
float  a,
float  b 
)

Given r and qt, carry out a Jacobi rotation on rows i and i+1 of each matrix. a and b are parameters of the rotation: $ =a/{a^2+b^2}, = b / {a^2 + b^2}$.

void rungeKuttaStep ( const DoubleVector y,
const DoubleVector dydx,
double  x,
double  h,
DoubleVector yout,
DoubleVector yerr,
DoubleVector(*)(double, const DoubleVector &)  derivs 
)

A single step of Runge Kutta (5th order), input: y and dydx (derivative of y), x is independent variable. yout is value after step. derivs is a user-supplied function

double sumOfExp ( double  a,
double  b 
)

This sums two exponentials nof the arguments in a way that tries to avoid underflows.

this underflow determines whether it's actually worth adding the two figures... any more than 15 digits and it's just going to be the maximum one.

double truncGaussWidthHalf ( long &  idum)

Gives 1/2 chance of being between 0 and 1 (flat), otherwise, a decreasing Gaussian

double zriddr ( double(*)(double)  func,
double  x1,
double  x2,
double  xacc 
)

returns a root of function func between x1 and x2 to accuracy xacc (x1 and x2 MUST bracket a root)