numerics.h

Go to the documentation of this file.

 
#ifndef NUMERICS_H
#define NUMERICS_H
 
 
#include "utils.h"
#include "mycomplex.h"
#include "def.h"
#include <iostream>
#include <cmath>
#include <cfloat>
#include <limits>
#include "linalg.h"
#include "dilogwrap.h"
using namespace softsusy;
 
int cubicRoots(double a, double b, double c, double d,
               DoubleVector & ans);
 
int cubicRootsInside(double a, double b, double c, DoubleVector & ans);
 
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);
 
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 &));
 
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 &)));
 
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); 
 
double integrandThreshbn(double x);
double bIntegral(int n, double p, double m1, double m2, double mt);
DoubleVector dd(double x, const DoubleVector & y);
 
double b0(double p, double m1, double m2, double q);
double b1(double p, double m1, double m2, double q);
double b22(double p,  double m1, double m2, double q);
double c0(double m1, double m2, double m3);
double d27(double m1, double m2, double m3, double m4);
double d0(double m1, double m2, double m3, double m4);
inline double a0(double m, double q) {
  if (fabs(m) < EPSTOL) return 0.;
  return sqr(m) * (1.0 - 2. * log(abs(m / q)));
}
inline double ffn(double p, double m1, double m2, double q) {
  return a0(m1, q) - 2.0 * a0(m2, q) - 
    (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) * 
    b0(p, m1, m2, q);
}
inline double gfn(double p, double m1, double m2, double q) {
  return (sqr(p) - sqr(m1) - sqr(m2)) * b0(p, m1, m2, q) - a0(m1, q) 
    - a0(m2, q); 
}
inline double hfn(double p, double m1, double m2, double q) {
  return 4.0 * b22(p, m1, m2, q) + gfn(p, m1, m2, q);
}
inline double b22bar(double p, double m1, double m2, double q) {
  return b22(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
}
 
Complex b0c(double p, double m1, double m2, double q);
Complex b1c(double p, double m1, double m2, double q);
Complex b22c(double p,  double m1, double m2, double q);
inline Complex ffnc(double p, double m1, double m2, double q) {
  return a0(m1, q) - 2.0 * a0(m2, q) - 
    (2.0 * sqr(p) + 2.0 * sqr(m1) - sqr(m2)) * 
    b0c(p, m1, m2, q);
}
inline Complex gfnc(double p, double m1, double m2, double q) {
  return (sqr(p) - sqr(m1) - sqr(m2)) * b0c(p, m1, m2, q) - a0(m1, q) 
    - a0(m2, q); 
}
inline Complex hfnc(double p, double m1, double m2, double q) {
  return 4.0 * b22c(p, m1, m2, q) + gfnc(p, m1, m2, q);
}
inline Complex b22barc(double p, double m1, double m2, double q) {
  return b22c(p, m1, m2, q) - 0.25 * a0(m1, q) - 0.25 * a0(m2, q);
}
 
/*inline double fB(const Complex & x) {
  //  if (x.real() < EPSTOL * 10. && x.imag() < EPSTOL * 10.) return -1.;
  return (log(Complex(1.0 - x)) - x * log(Complex(1.0 - 1.0 / x)) 
         - Complex(1.0)).real();
         } */
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);
 
DoubleVector getRandomDirection(int n, int & numChanged, long & idum); 
 
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>&);
 
double trapzd(double (*func)(double), double a, double b, int n, 
              double tol = 1.0e-3);
double qtrap(double (*func)(double), double a, double b, double tol);
double midpnt(double (*func)(double), double a, double b, int n);
double qromb(double (*func)(double), double a, double b, double EPS);
double qromb2(double (*func)(double), double a, double b, double EPS);
 
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);
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);
double den(double a, int b); 
double log1minusx(double x);
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
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);
/* allocate an int vector with subscript range v[nl..nh] */
int *ivector(long nl, long nh);
/* free an int vector allocated with ivector() */
void free_ivector(int *v, long nl, long nh);
void lubksb(const DoubleMatrix & a, int n, int *indx, DoubleVector & b);
bool newt(DoubleVector & x, 
          int (*vecfunc)(const DoubleVector &, void*, DoubleVector &),
          void* params = NULL);
void load(float x, const DoubleVector & v, DoubleVector & y);
void score(float x, const DoubleVector & y, DoubleVector & f);
//void shoot(const DoubleVector & v, 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);
 
  // returns sqrt(f) for f>0 
  double ccbSqrt(double f);
  // returns sqrt(f) for f>0 or -sqrt(|f|) for f<0
  double signedSqrt(double f);
  double signedSqr(double f);
 
double kinFn(double m1, double m2, double m3);
double lambda(double a, double b, double c); 
#endif