Impala Documentation

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