softsusy is hosted by Hepforge, IPPP Durham
 SOFTSUSY  4.1
numerics.h File Reference

numerical routines - differential equation solver, differentiator and function minimiser for instance More...

#include "utils.h"
#include "mycomplex.h"
#include "def.h"
#include <iostream>
#include <cmath>
#include <cfloat>
#include <limits>
#include "linalg.h"
#include "dilogwrap.h"

Go to the source code of this file.

## Macros

#define SIGN(a, b)   ((b) >= 0.0 ? fabs(a) : -fabs(a))
sign of inputs

## Functions

int cubicRoots (double a, double b, double c, double d, DoubleVector &ans)

int cubicRootsInside (double a, double b, double c, DoubleVector &ans)
cubic roots for x^3 + a*x + b*x^2 + c*x^4 = 0 More...

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...

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

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

double lnLPoisson (unsigned k, double lambda)

double LPoisson (unsigned k, double lambda)

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.

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

double calcDerivative (double(*func)(double, void *), double x, double h, double *err, void *params=NULL)

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

void shft2 (double &a, double &b, double &c)

void shft3 (double &a, double &b, double &c, double &d)
a=b, b=c and c=d

double integrandThreshbn (double x)
For calculation of PV functions.

double bIntegral (int n, double p, double m1, double m2, double mt)
Returns real part of b function, less accurate than analytic expressions.

DoubleVector dd (double x, const DoubleVector &y)

double b0 (double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part. More...

double b1 (double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part. More...

double b22 (double p, double m1, double m2, double q)
Passarino-Veltman function definition: real part. More...

double c0 (double m1, double m2, double m3)
Passarino-Veltman function definition: real part.

double d27 (double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.

double d0 (double m1, double m2, double m3, double m4)
Passarino-Veltman function definition: real part.

double a0 (double m, double q)
inlined PV functions: real part

double ffn (double p, double m1, double m2, double q)
inlined PV functions: real part

double gfn (double p, double m1, double m2, double q)
inlined PV functions: real part

double hfn (double p, double m1, double m2, double q)
inlined PV functions: real part

double b22bar (double p, double m1, double m2, double q)
inlined PV functions: real part

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

Complex b1c (double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer. More...

Complex b22c (double p, double m1, double m2, double q)
Passarino-Veltman function definition: complex answer. More...

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

Complex gfnc (double p, double m1, double m2, double q)
inlined PV functions: complex answer

Complex hfnc (double p, double m1, double m2, double q)
inlined PV functions: complex answer

Complex b22barc (double p, double m1, double m2, double q)
inlined PV functions: complex answer

double fB (const Complex &x)

Complex fBc (const Complex &x)

double dilogarg (double t)

double dilog (double x)

double integrandThreshbnr (double x)

Complex fnfn (double x)

double gasdev (long &idum)

double truncGaussWidthHalf (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.

DoubleVector getRandomDirection (int n, int &numChanged, long &idum)
returns a random direction: total number of dimensions=n

double calcCL (double cl, const DoubleVector &l)

double calc1dFraction (double y, const DoubleVector &l)

double fps (double z)

double fs (double z)

double ffbar (double z)

Complex dilog (const Complex &x)

std::complex< long double > dilogcl (const std::complex< long double > &)
Complex dilogarithm $$\mathrm{Li}_2(z)$$. More...

double trapzd (double(*func)(double), double a, double b, int n, double tol=1.0e-3)
Trapezoidal function integration to accuracy tol.

double qtrap (double(*func)(double), double a, double b, double tol)
Driver for integration.

double midpnt (double(*func)(double), double a, double b, int n)
Mid point function integration.

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.

double sumOfExp (double a, double b)

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 fin (double mm1, double mm2)
useful for 2-loop mb/mt corrections

double den (double a, int b)

double log1minusx (double x)
This function DOESNT WORK YET!!! More...

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

void fdjac (int n, DoubleVector x, const DoubleVector &fvec, DoubleMatrix &df, int(*vecfunc)(const DoubleVector &, void *, DoubleVector &), void *params=NULL)

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)

int * ivector (long nl, long nh)

void free_ivector (int *v, long nl, long nh)

void lubksb (const DoubleMatrix &a, int n, int *indx, DoubleVector &b)
Get rid of int n.

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. More...

void load (float x, const DoubleVector &v, DoubleVector &y)
calculates the n-vector y, given freely specifiable values v(1..n2) at x1

void score (float x, const DoubleVector &y, DoubleVector &f)

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

numerical routines - differential equation solver, differentiator and function minimiser for instance

## ◆ accurateSqrt1Plusx()

 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

## ◆ b0()

 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

alternative form: should be more accurate

## ◆ b0c()

 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

alternative form: should be more accurate

## ◆ b1()

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

Passarino-Veltman function definition: real part.

Passarino-Veltman function definition: real part.

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

< checked

## ◆ b1c()

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

Passarino-Veltman function definition: complex answer.

Passarino-Veltman function definition: complex answer.

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

< checked

< checked

## ◆ b22()

 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

## ◆ b22c()

 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

## ◆ bin()

 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...)

## ◆ calc1dFraction()

 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

## ◆ calcCL()

 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.

## ◆ calcDerivative() [1/2]

 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

## ◆ calcDerivative() [2/2]

 double calcDerivative ( double(*)(double, void *) func, double x, double h, double * err, void * params = NULL )

DH: overloaded version allowing to pass additional parameters to the function

## ◆ cauchyRan()

 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.

## ◆ cubicRoots()

 int cubicRoots ( double a, double b, double c, double d, DoubleVector & ans )

gives DoubleVector ans with all n real roots (returned) in increasing absolute order for: a x^3 + b x^2 + c x + d=0. Needs an upgrade to include complex roots

order the roots in increasing absolute order

It's a linear equation

it's a quadratic really

## ◆ cubicRootsInside()

 int cubicRootsInside ( double a, double b, double c, DoubleVector & ans )

cubic roots for x^3 + a*x + b*x^2 + c*x^4 = 0

Three real roots

only one real root

## ◆ dgauss() [1/2]

 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

## ◆ dgauss() [2/2]

 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

## ◆ dilogcl()

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

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

Parameters
 z complex argument
Note
Implementation translated from SPheno to C++
Returns
$$\mathrm{Li}_2(z)$$

## ◆ display3x3RealMixing()

 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

## ◆ fB()

 double fB ( const Complex & x )

First, special cases at problematic points

## ◆ fBc()

 Complex fBc ( const Complex & x )

First, special cases at problematic points

## ◆ fdjac()

 void fdjac ( int n, DoubleVector x, const DoubleVector & fvec, DoubleMatrix & df, int(*)(const DoubleVector &, void *, DoubleVector &) vecfunc, void * params = NULL )

Forward difference approximation to Jacobian matrix. On input, x is the point to be evaluated, fvec is the vector of function values at the point, and vecfunc(n, x, f) is the Jacobian array

## ◆ ffnc()

 Complex ffnc ( double p, double m1, double m2, double q )
inline

Passarino-Veltman function definition: complex answer inlined PV functions: complex answer

## ◆ findMinimum()

 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

## ◆ fps()

 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

## ◆ gasdev()

 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

## ◆ integrateReversibly()

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

Reversible integrate with numSteps steps between tInitial and tFinal. Returns true if there's an error

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

## ◆ kinFn()

 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}$$.

## ◆ lnLPoisson()

 double lnLPoisson ( unsigned k, double lambda )

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

## ◆ lnsrch()

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

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

## ◆ log1minusx()

 double log1minusx ( double x )

This function DOESNT WORK YET!!!

1/a^b

Find largest power that we need from the expansion

## ◆ LPoisson()

 double LPoisson ( unsigned k, double lambda )

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

## ◆ midPtStep()

 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

## ◆ newt()

 bool newt ( DoubleVector & x, int(*)(const DoubleVector &, void *, DoubleVector &) vecfunc, void * params = NULL )

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

## ◆ polint()

 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).

## ◆ qrdcmp()

 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 \ldots Q_{n-1}$ where $Q_j=1-u_j \otimes u_j/c_j$. sing returns true (1) if singularity is encountered (but decomposition is still completed) otherwise false (0).

## ◆ qromb()

 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

## ◆ ran1()

 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

## ◆ rotate()

 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: $\cos \theta=a/\sqrt{a^2+b^2}, \sin \theta = b / \sqrt{a^2 + b^2}$.

## ◆ rungeKuttaStep()

 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

## ◆ score()

 void score ( float x, const DoubleVector & y, DoubleVector & f )

Gives a discrepancy vector f[1..n2] from ending boundary conditions at endpoint x2 given values for y there

## ◆ shft2()

 void shft2 ( double & a, double & b, double & c )

a=b and b=c

## ◆ sumOfExp()

 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.

## ◆ truncGaussWidthHalf()

 double truncGaussWidthHalf ( long & idum )

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

## ◆ zriddr()

 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)