template<class DstArrayT, class SrcArrayT>
Definition at line 281 of file WeibullNgbPnLoop.h. References M_LN2, M_PI, and MMGAMMA. 00282 { 00283 /* force into the interval 1..3 */ 00284 if( x < 0.0 ) 00285 return DiGamma(1.0-x)+M_PI/tan(M_PI*(1.0-x)) ; /* refection formula */ 00286 else if( x < 1.0 ) 00287 return DiGamma(1.0+x)-1.0/x ; 00288 else if ( x == 1.0) 00289 return -MMGAMMA ; 00290 else if ( x == 2.0) 00291 return 1.0-MMGAMMA ; 00292 else if ( x == 3.0) 00293 return 1.5-MMGAMMA ; 00294 else if ( x > 3.0) 00295 /* duplication formula */ 00296 return 0.5*(DiGamma(x/2.0)+DiGamma((x+1.0)/2.0))+M_LN2 ; 00297 else 00298 { 00299 static double Kncoe[] = 00300 { .30459198558715155634315638246624251, 00301 .72037977439182833573548891941219706, -.12454959243861367729528855995001087, 00302 .27769457331927827002810119567456810e-1, -.67762371439822456447373550186163070e-2, 00303 .17238755142247705209823876688592170e-2, -.44817699064252933515310345718960928e-3, 00304 .11793660000155572716272710617753373e-3, -.31253894280980134452125172274246963e-4, 00305 .83173997012173283398932708991137488e-5, -.22191427643780045431149221890172210e-5, 00306 .59302266729329346291029599913617915e-6, -.15863051191470655433559920279603632e-6, 00307 .42459203983193603241777510648681429e-7, -.11369129616951114238848106591780146e-7, 00308 .304502217295931698401459168423403510e-8, -.81568455080753152802915013641723686e-9, 00309 .21852324749975455125936715817306383e-9, -.58546491441689515680751900276454407e-10, 00310 .15686348450871204869813586459513648e-10, -.42029496273143231373796179302482033e-11, 00311 .11261435719264907097227520956710754e-11, -.30174353636860279765375177200637590e-12, 00312 .80850955256389526647406571868193768e-13, -.21663779809421233144009565199997351e-13, 00313 .58047634271339391495076374966835526e-14, -.15553767189204733561108869588173845e-14, 00314 .41676108598040807753707828039353330e-15, -.11167065064221317094734023242188463e-15 } ; 00315 00316 double Tn_1 = 1.0 ; /* T_{n-1}(x), started at n=1 */ 00317 double Tn = x-2.0 ; /* T_{n}(x) , started at n=1 */ 00318 double resul = Kncoe[0] + Kncoe[1]*Tn ; 00319 00320 x -= 2.0 ; 00321 00322 for(int n = 2 ; n < sizeof(Kncoe)/sizeof(double) ;n++) 00323 { 00324 const double Tn1 = 2.0 * x * Tn - Tn_1 ; /* Chebyshev recursion, Eq. 22.7.4 Abramowitz-Stegun */ 00325 resul += Kncoe[n]*Tn1 ; 00326 Tn_1 = Tn ; 00327 Tn = Tn1 ; 00328 } 00329 return resul ; 00330 } 00331 } // DiGamma(double x)
|