source: JEM-EUSO/esaf_cc_at_lal/packages/common/base/src/utils.cc @ 114

//
// $Id: utils.cc 2955 2011-06-28 18:38:08Z mabl $
//
#include <math.h>
#include "utils.hh"



//_____________________________________________________________________________________________________
pair<Int_t, Double_t*> &findRoots(Double_t a, Double_t b, Double_t c) {
    //
    // Solutions of:
    // a x^2 + b x + c = 0
    // 
    // q = -1/2 (b + sign(b) sqrt(b^2 - 4 a c)
    //
    // x_1 = q / a
    // x_2 = c / q
    //
    // cfr. 
    // Numerical Recipes in Fortran (5.6 Quadratic and Cubic Equations)
    // W. H. Press et al.
    // Cambridge University Press
    // Second Edition
    //
    
    // \Delta=b^2-4ac
    Double_t det = b*b - 4*a*c;
    static Double_t roots[2];
    static pair<Int_t, Double_t*> p;
    p.second = roots;
    

    // det < 0 OR a=b=0: no solution
    if(det < -kTolerance || (fabs(a) < kTolerance && fabs(b) < kTolerance) ) {
        p.first=0;
        p.second=0;
        if((fabs(a) < kTolerance && fabs(b) < kTolerance))
            cout << "\nWARNING : <findRoots> algo called with coeff. a=b=0"<<endl;
    }
    else {
        Double_t q = -0.5 * (b + sign(b)*sqrt(det));

        // if det=0 OR a=0 OR b=c=0: one solution
        if(fabs(det) < kTolerance || fabs(a) < kTolerance || (fabs(b) < kTolerance && fabs(c) < kTolerance) ) {
            if(fabs(a) < kTolerance) p.second[0] = c/q;  // b cannot be null here, so q
            else                     p.second[0] = q/a;
            p.first = 1;
        }
        
        // if det > 0 : two solutions
        else {
            p.second[0] = q/a;
            p.second[1] = c/q;
            p.first = 2;
        }
    }
    
    return p;
}

//_____________________________________________________________________________________________________
Int_t zbrac(Float_t (*func)(Float_t), Float_t *x1, Float_t *x2) {
    //
    // Numerical Recipes in C. Chapter 9. Root Finding and Nonlinear Sets of Equations
    //
    
    if(*x1 == *x2) return -1;
    Float_t f1 = (*func)(*x1);
    Float_t f2 = (*func)(*x2);
    for (Int_t j=1;j<NTRY;j++) {
        if (f1*f2 < 0) return 1;
        if (fabs(f1) < fabs(f2))
            f1 = (*func)(*x1 += FACTOR*(*x1-*x2));
        else
            f2 = (*func)(*x2 += FACTOR*(*x2-*x1));
    }
    return 0;
}

//_____________________________________________________________________________________________________
Float_t rtbis(Float_t (*func)(Float_t), Float_t x1, Float_t x2, Float_t xacc) {
    //
    //
    //
    
    Float_t dx,f,fmid,xmid,rtb;
    f = (*func)(x1);
    fmid = (*func)(x2);
    if(f*fmid >=0) return -1;
    rtb = f<0 ? (dx = x2-x1,x1) : (dx=x1-x2,x2);
    for (Int_t j=1;j<=JMAX;j++) {
        fmid = (*func)(xmid=rtb+(dx*=0.5));
        if (fmid <=0) rtb = xmid;
        if (fabs(dx) < xacc || fmid == 0) return rtb;
    }
    return 0;
}