Impala Documentation

template<class DstArrayT, class SrcArrayT>
void Impala::Core::Feature::IntWeibullNgbPnLoop< DstArrayT, SrcArrayT >::dpsifn	(	double	x,
		int	n,
		int	kode,
		int	m,
		double *	ans,
		int *	nz,
		int *	ierr
	)			`[inline]`
Definition at line 415 of file IntWeibullNgbPnLoop.h.
References Impala::Core::Feature::IntWeibullNgbPnLoop< DstArrayT, SrcArrayT >::fmax2(), Impala::Core::Feature::IntWeibullNgbPnLoop< DstArrayT, SrcArrayT >::fmin2(), Impala::Core::Feature::IntWeibullNgbPnLoop< DstArrayT, SrcArrayT >::imin2(), M_PI, ML_NAN, ML_POSINF, n_max, Impala::Core::Feature::IntWeibullNgbPnLoop< DstArrayT, SrcArrayT >::Rf_d1mach(), and Impala::Core::Feature::IntWeibullNgbPnLoop< DstArrayT, SrcArrayT >::Rf_i1mach().
Referenced by Impala::Core::Feature::IntWeibullNgbPnLoop< DstArrayT, SrcArrayT >::TriGamma().
00416 {
00417   const static double bvalues[] = {/* Bernoulli Numbers */
00418           1.00000000000000000e+00,
00419           -5.00000000000000000e-01,
00420           1.66666666666666667e-01,
00421           -3.33333333333333333e-02,
00422           2.38095238095238095e-02,
00423           -3.33333333333333333e-02,
00424           7.57575757575757576e-02,
00425           -2.53113553113553114e-01,
00426           1.16666666666666667e+00,
00427           -7.09215686274509804e+00,
00428           5.49711779448621554e+01,
00429           -5.29124242424242424e+02,
00430           6.19212318840579710e+03,
00431           -8.65802531135531136e+04,
00432           1.42551716666666667e+06,
00433           -2.72982310678160920e+07,
00434           6.01580873900642368e+08,
00435           -1.51163157670921569e+10,
00436           4.29614643061166667e+11,
00437           -1.37116552050883328e+13,
00438           4.88332318973593167e+14,
00439           -1.92965793419400681e+16
00440         };
00441 
00442         int i, j, k, mm, mx, nn, np, nx, fn;
00443         double arg, den, elim, eps, fln, fx, rln, rxsq,
00444           r1m4, r1m5, s, slope, t, ta, tk, tol, tols, tss, tst,
00445           tt, t1, t2, wdtol, xdmln, xdmy, xinc, xln = 0.0 /* -Wall */, 
00446           xm, xmin, xq, yint;
00447         double trm[23], trmr[n_max + 1];
00448 
00449         *ierr = 0;
00450         if (n < 0 || kode < 1 || kode > 2 || m < 1) {
00451           *ierr = 1;
00452           return;
00453         }
00454         if (x <= 0.) {
00455           /* useAbramowitz & Stegun 6.4.7 "Reflection Formula"
00456            *psi(k, x) = (-1)^k psi(k, 1-x)-  pi^{n+1} (d/dx)^n cot(x)
00457            */
00458           if (x == (long)x) {
00459             /* non-positive integer : +Inf or NaN depends on n */
00460             for(j=0; j < m; j++) /* k = j + n : */
00461               ans[j] = ((j+n) % 2) ? ML_POSINF : ML_NAN;
00462             return;
00463           }
00464           dpsifn(1. - x, n, /*kode = */ 1, m, ans, nz, ierr);
00465           /* ans[j] == (-1)^(k+1) / gamma(k+1) * psi(k, 1 - x)
00466            *     for j = 0:(m-1) ,k = n + j
00467            */
00468 
00469           /* Cheat for now: only work for m = 1, n in {0,1,2,3} : */
00470           if(m > 1 || n > 3) {/* doesn't happen for digamma() .. pentagamma() */
00471             /* not yet implemented */
00472             *ierr = 4; return;
00473           }
00474           x *= M_PI; /* pi * x */
00475           if (n == 0)
00476             tt = cos(x)/sin(x);
00477           else if (n == 1)
00478             tt = -1/pow(sin(x),2);
00479           else if (n == 2)
00480             tt = 2*cos(x)/pow(sin(x),3);
00481           else if (n == 3)
00482             tt = -2*(2*pow(cos(x),2) + 1)/pow(sin(x),4);
00483           else /* can not happen! */
00484             tt = ML_NAN;
00485           /* end cheat */
00486 
00487           s = (n % 2) ? -1. : 1.;/* s = (-1)^n */
00488           /* t := pi^(n+1) * d_n(x) / gamma(n+1), where
00489            *   d_n(x) := (d/dx)^n cot(x)*/
00490           t1 = t2 = s = 1.;
00491           for(k=0, j=k-n; j < m; k++, j++, s = -s) {
00492             /* k == n+j , s = (-1)^k */
00493             t1 *= M_PI;/* t1 == pi^(k+1) */
00494             if(k >= 2)
00495               t2 *= k;/* t2 == k! == gamma(k+1) */
00496             if(j >= 0) /* by cheat above,  tt === d_k(x) */
00497               ans[j] = s*(ans[j] + t1/t2 * tt);
00498           }
00499           if (n == 0 && kode == 2)
00500             ans[0] += xln;
00501           return;
00502         } /* x <= 0 */
00503 
00504         *nz = 0;
00505         mm = m;
00506         nx = imin2(-Rf_i1mach(15), Rf_i1mach(16));/* = 1021 */
00507         r1m5 = Rf_d1mach(5);
00508         r1m4 = Rf_d1mach(4) * 0.5;
00509         wdtol = fmax2(r1m4, 0.5e-18); /* 1.11e-16 */
00510 
00511         /* elim = approximate exponential over and underflow limit */
00512 
00513         elim = 2.302 * (nx * r1m5 - 3.0);/* = 700.6174... */
00514         xln = log(x);
00515         for(;;) {
00516           nn = n + mm - 1;
00517           fn = nn;
00518           t = (fn + 1) * xln;
00519 
00520           /* overflow and underflow test for small and large x */
00521 
00522           if (fabs(t) > elim) {
00523             if (t <= 0.0) {
00524               *nz = 0;
00525               *ierr = 2;
00526               return;
00527             }
00528           }
00529           else {
00530             if (x < wdtol) {
00531               ans[0] = pow(x, -n-1.0);
00532               if (mm != 1) {
00533                 for(k = 1; k < mm ; k++)
00534                   ans[k] = ans[k-1] / x;
00535               }
00536               if (n == 0 && kode == 2)
00537                 ans[0] += xln;
00538               return;
00539             }
00540 
00541             /* compute xmin and the number of terms of the series,  fln+1 */
00542 
00543             rln = r1m5 * Rf_i1mach(14);
00544             rln = fmin2(rln, 18.06);
00545             fln = fmax2(rln, 3.0) - 3.0;
00546             yint = 3.50 + 0.40 * fln;
00547             slope = 0.21 + fln * (0.0006038 * fln + 0.008677);
00548             xm = yint + slope * fn;
00549             mx = (int)xm + 1;
00550             xmin = mx;
00551             if (n != 0) {
00552               xm = -2.302 * rln - fmin2(0.0, xln);
00553               arg = xm / n;
00554               arg = fmin2(0.0, arg);
00555               eps = exp(arg);
00556               xm = 1.0 - eps;
00557               if (fabs(arg) < 1.0e-3)
00558                 xm = -arg;
00559               fln = x * xm / eps;
00560               xm = xmin - x;
00561               if (xm > 7.0 && fln < 15.0)
00562                 break;
00563             }
00564             xdmy = x;
00565             xdmln = xln;
00566             xinc = 0.0;
00567             if (x < xmin) {
00568               nx = (int)x;
00569               xinc = xmin - nx;
00570               xdmy = x + xinc;
00571               xdmln = log(xdmy);
00572             }
00573 
00574             /* generate w(n+mm-1, x) by the asymptotic expansion */
00575 
00576             t = fn * xdmln;
00577             t1 = xdmln + xdmln;
00578             t2 = t + xdmln;
00579             tk = fmax2(fabs(t), fmax2(fabs(t1), fabs(t2)));
00580             if (tk <= elim)
00581               goto L10;
00582           }
00583           nz++;
00584           mm--;
00585           ans[mm] = 0.;
00586           if (mm == 0)
00587             return;
00588         }
00589         nn = (int)fln + 1;
00590         np = n + 1;
00591         t1 = (n + 1) * xln;
00592         t = exp(-t1);
00593         s = t;
00594         den = x;
00595         for(i=1; i <= nn; i++) {
00596           den += 1.;
00597           trm[i] = pow(den, (double)-np);
00598           s += trm[i];
00599         }
00600         ans[0] = s;
00601         if (n == 0 && kode == 2)
00602           ans[0] = s + xln;
00603 
00604         if (mm != 1) { /* generate higher derivatives, j > n */
00605 
00606           tol = wdtol / 5.0;
00607           for(j = 1; j < mm; j++) {
00608             t /= x;
00609             s = t;
00610             tols = t * tol;
00611             den = x;
00612             for(i=1; i <= nn; i++) {
00613               den += 1.;
00614               trm[i] /= den;
00615               s += trm[i];
00616               if (trm[i] < tols)
00617                 break;
00618             }
00619             ans[j] = s;
00620           }
00621         }
00622         return;
00623 
00624       L10:
00625         tss = exp(-t);
00626         tt = 0.5 / xdmy;
00627         t1 = tt;
00628         tst = wdtol * tt;
00629         if (nn != 0)
00630           t1 = tt + 1.0 / fn;
00631         rxsq = 1.0 / (xdmy * xdmy);
00632         ta = 0.5 * rxsq;
00633         t = (fn + 1) * ta;
00634         s = t * bvalues[2];
00635         if (fabs(s) >= tst) {
00636           tk = 2.0;
00637           for(k = 4; k <= 22; k++) {
00638             t = t * ((tk + fn + 1)/(tk + 1.0))*((tk + fn)/(tk + 2.0)) * rxsq;
00639             trm[k] = t * bvalues[k-1];
00640             if (fabs(trm[k]) < tst)
00641               break;
00642             s += trm[k];
00643             tk += 2.;
00644           }
00645         }
00646         s = (s + t1) * tss;
00647         if (xinc != 0.0) {
00648 
00649           /* backward recur from xdmy to x */
00650 
00651           nx = (int)xinc;
00652           np = nn + 1;
00653           if (nx > n_max) {
00654             *nz = 0;
00655             *ierr = 3;
00656             return;
00657           }
00658           else {
00659             if (nn==0)
00660               goto L20;
00661             xm = xinc - 1.0;
00662             fx = x + xm;
00663 
00664             /* this loop should not be changed. fx is accurate when x is small */
00665             for(i = 1; i <= nx; i++) {
00666               trmr[i] = pow(fx, (double)-np);
00667               s += trmr[i];
00668               xm -= 1.;
00669               fx = x + xm;
00670             }
00671           }
00672         }
00673         ans[mm-1] = s;
00674         if (fn == 0)
00675           goto L30;
00676 
00677         /* generate lower derivatives,  j < n+mm-1 */
00678 
00679         for(j = 2; j <= mm; j++) {
00680           fn--;
00681           tss *= xdmy;
00682           t1 = tt;
00683           if (fn!=0)
00684             t1 = tt + 1.0 / fn;
00685           t = (fn + 1) * ta;
00686           s = t * bvalues[2];
00687           if (fabs(s) >= tst) {
00688             tk = 4 + fn;
00689             for(k=4; k <= 22; k++) {
00690               trm[k] = trm[k] * (fn + 1) / tk;
00691               if (fabs(trm[k]) < tst)
00692                 break;
00693               s += trm[k];
00694               tk += 2.;
00695             }
00696           }
00697           s = (s + t1) * tss;
00698           if (xinc != 0.0) {
00699             if (fn == 0)
00700               goto L20;
00701             xm = xinc - 1.0;
00702             fx = x + xm;
00703             for(i=1 ; i<=nx ; i++) {
00704               trmr[i] = trmr[i] * fx;
00705               s += trmr[i];
00706               xm -= 1.;
00707               fx = x + xm;
00708             }
00709           }
00710           ans[mm - j] = s;
00711           if (fn == 0)
00712             goto L30;
00713         }
00714         return;
00715 
00716       L20:
00717         for(i = 1; i <= nx; i++)
00718           s += 1. / (x + nx - i);
00719 
00720       L30:
00721         if (kode!=2)
00722           ans[0] = s - xdmln;
00723         else if (xdmy != x) {
00724           xq = xdmy / x;
00725           ans[0] = s - log(xq);
00726         }
00727         return;
00728       } /* dpsifn() */
Here is the call graph for this function: