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

• Project: SOFTSUSY
• Author: Ben Allanach
• Manual: hep-ph/0104145, Comp. Phys. Comm. 143 (2002) 305
• Webpage: http://hepforge.cedar.ac.uk/softsusy/

Function Documentation

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<long double> 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)