Impala Documentation

template<class HistT>
double Impala::Core::Histogram::FitWeibullMarginal< HistT >::DiGamma ( double x ) const [inline, private]
Definition at line 304 of file FitWeibullMarginal.h.
Referenced by Impala::Core::Histogram::FitWeibullMarginal< HistT >::GammaEst().
00305     {
00306         /* force into the interval 1..3 */
00307         if( x < 0.0 )
00308             return DiGamma(1.0-x)+M_PI/tan(M_PI*(1.0-x)) ;  /* refection formula */
00309         else if( x < 1.0 )
00310             return DiGamma(1.0+x)-1.0/x ;
00311         else if ( x == 1.0)
00312             return -MMGAMMA ;
00313         else if ( x == 2.0)
00314             return 1.0-MMGAMMA ;
00315         else if ( x == 3.0)
00316             return 1.5-MMGAMMA ;
00317         else if ( x > 3.0)
00318             /* duplication formula */
00319             return 0.5*(DiGamma(x/2.0)+DiGamma((x+1.0)/2.0))+M_LN2 ;
00320         else
00321         {
00322             static double Kncoe[] =
00323                 { .30459198558715155634315638246624251,
00324                   .72037977439182833573548891941219706, -.12454959243861367729528855995001087,
00325                   .27769457331927827002810119567456810e-1, -.67762371439822456447373550186163070e-2,
00326                   .17238755142247705209823876688592170e-2, -.44817699064252933515310345718960928e-3,
00327                   .11793660000155572716272710617753373e-3, -.31253894280980134452125172274246963e-4,
00328                   .83173997012173283398932708991137488e-5, -.22191427643780045431149221890172210e-5,
00329                   .59302266729329346291029599913617915e-6, -.15863051191470655433559920279603632e-6,
00330                   .42459203983193603241777510648681429e-7, -.11369129616951114238848106591780146e-7,
00331                   .304502217295931698401459168423403510e-8, -.81568455080753152802915013641723686e-9,
00332                   .21852324749975455125936715817306383e-9, -.58546491441689515680751900276454407e-10,
00333                   .15686348450871204869813586459513648e-10, -.42029496273143231373796179302482033e-11,
00334                   .11261435719264907097227520956710754e-11, -.30174353636860279765375177200637590e-12,
00335                   .80850955256389526647406571868193768e-13, -.21663779809421233144009565199997351e-13,
00336                   .58047634271339391495076374966835526e-14, -.15553767189204733561108869588173845e-14,
00337                   .41676108598040807753707828039353330e-15, -.11167065064221317094734023242188463e-15 } ;
00338 
00339             double Tn_1 = 1.0 ; /* T_{n-1}(x), started at n=1 */
00340             double Tn = x-2.0 ; /* T_{n}(x) , started at n=1 */
00341             double resul = Kncoe[0] + Kncoe[1]*Tn ;
00342 
00343             x -= 2.0 ;
00344 
00345             for(int n = 2 ; n < sizeof(Kncoe)/sizeof(double) ;n++)
00346             {
00347                 const double Tn1 = 2.0 * x * Tn - Tn_1 ;    /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */
00348                 resul += Kncoe[n]*Tn1 ;
00349                 Tn_1 = Tn ;
00350                 Tn = Tn1 ;
00351             }
00352             return resul ;
00353         }
00354     }