The R Project SVN R

Rev

Rev 89290 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 89290 Rev 89292
Line 215... Line 215...
215 
 
 215 
 
216 
    const static double twopi1 = 6.28125;			// twopi1 = first few significant digits of 2\pi
 216 
    const static double twopi1 = 6.28125;			// twopi1 = first few significant digits of 2\pi
217 
    const static double twopi2 =  .001935307179586476925286767; /* twopi2 = (2*\pi - twopi1) to working precision, i.e.,
 217 
    const static double twopi2 =  .001935307179586476925286767; /* twopi2 = (2*\pi - twopi1) to working precision, i.e.,
218 
								 * twopi1 + twopi2 = 2 \pi to extra precision.
 218 
								 * twopi1 + twopi2 = 2 \pi to extra precision.
219 
 --------------------------------------------------------------------- */
 219 
 --------------------------------------------------------------------- */
- 
 
 220 
#define very_small_nu  0x1p-800 // 2^-800 = 1.4996968....e-241
220 
 
 221 
 
221 
    --b; /* so, we use  b[1] .. b[nb]  in the code below */
 222 
    --b; /* so, we use  b[1] .. b[nb]  in the code below */
222 
 
 223 
 
223 
    double nu = *alpha, // in [0, 1)   {ensured by caller bessel_j*()}
 224 
    double nu = *alpha, // in [0, 1)   {ensured by caller bessel_j*()}
224 
	twonu = ldexp(nu,1); // = 2 nu = nu+nu
 225 
	twonu = ldexp(nu,1); // = 2 nu = nu+nu
Line 243... Line 244...
243 
	/*===================================================================
 244 
	/*===================================================================
244 
	  Branch into  3 cases :
 245 
	  Branch into  3 cases :
245 
	  1) use 2-term ascending series for small X
 246 
	  1) use 2-term ascending series for small X
246 
	  2) use asymptotic form for large X when NB is not too large
 247 
	  2) use asymptotic form for large X when NB is not too large
247 
	  3) use recursion otherwise;
 248 
	  3) use recursion otherwise;
- 
 
 249 
	   3b:  if 0 < |nu| = |alpha| < very_small_nu, use nu = very_small_nu
248 
	  ===================================================================*/
 250 
	  ===================================================================*/
249 
 
 251 
 
250 
	double alpem, alp2em, aa, bb, cc, p, s, en, sum, tover;
 252 
	double alpem, alp2em, aa, bb, cc, p, s, en, sum, tover;
251 
 
 253 
 
252 
	if (*x < rtnsig_BESS) { // x < 1e-4  here
 254 
	if (*x < rtnsig_BESS) { // x < 1e-4  here
Line 356... Line 358...
356 
	       -------------------------------------------------------- ============= branch 3)
 358 
	       -------------------------------------------------------- ============= branch 3)
357 
	       Use recurrence to generate results.
 359 
	       Use recurrence to generate results.
358 
	       First initialize the calculation of P*S.
 360 
	       First initialize the calculation of P*S.
359 
	       -------------------------------------------------------- */
 361 
	       -------------------------------------------------------- */
360 
 
 362 
 
- 
 
 363 
	    if(nu != 0. && fabs(nu) < very_small_nu) {
- 
 
 364 
		nu = (nu < 0.) ? -very_small_nu : very_small_nu; // in R <= 4.5.2  besselJ(2, 2e-16) gave 1.119e+15
- 
 
 365 
		twonu = ldexp(nu, 1);
- 
 
 366 
	    }
- 
 
 367 
 
361 
	    int nbmx = *nb - intx; // = nb - floor(x)
 368 
	    int nbmx = *nb - intx; // = nb - floor(x)
362 
	    n = intx + 1;
 369 
	    n = intx + 1;
363 
	    en = (double)(n + n) + twonu;
 370 
	    en = (double)(n + n) + twonu;
364 
	    p = en / *x;
 371 
	    p = en / *x;
365 
	    /* ---------------------------------------------------
 372 
	    /* ---------------------------------------------------
Line 490... Line 497...
490 
	      Store b[NB].
 497 
	      Store b[NB].
491 
	      --------------------------------------------------*/
 498 
	      --------------------------------------------------*/
492 
	    b[n] = aa;
 499 
	    b[n] = aa;
493 
	    if (nend >= 0) {
 500 
	    if (nend >= 0) {
494 
		if (n <= 1) {
 501 
		if (n <= 1) {
495 
		    sum += b[1] * ((nu + 1. == 1.) ? 1. : nu);
 502 
		    sum += b[1] * ((nu == 0.) ? 1. : nu); // as |nu| >=  very_small_nu
496 
		    goto L250;
 503 
		    goto L250;
497 
		}
 504 
		}
498 
		else {/*-- nb >= 2 : ---------------------------
 505 
		else {/*-- nb >= 2 : ---------------------------
499 
			Calculate and store b[NB-1].
 506 
			Calculate and store b[NB-1].
500 
			----------------------------------------*/
 507 
			----------------------------------------*/
Line 548... Line 555...
548 
 
 555 
 
549 
L250:
 556 
L250:
550 
	    /* ---------------------------------------------------
 557 
	    /* ---------------------------------------------------
551 
	       Normalize.  Divide all b[N] by sum.
 558 
	       Normalize.  Divide all b[N] by sum.
552 
	       ---------------------------------------------------*/
 559 
	       ---------------------------------------------------*/
553 
/*	    if (nu + 1. != 1.) poor test */
 560 
	    // NB. ensured above that |nu| >= very_small_nu
554 
	    if(fabs(nu) > 1e-15)
 561 
	    if(nu != 0.) { /* was if(fabs(nu) > very_small_nu) , was if(nu + 1. != 1.); then '> 1e-15' .. */
555 
		sum *= (Rf_gamma_cody(nu) * pow(.5* *x, -nu));
 562 
		sum *= (Rf_gamma_cody(nu) * pow(.5* *x, -nu));
- 
 
 563 
	    }
556 
 
 564 
 
- 
 
 565 
#ifdef UNDERFLOW_NOT_GOOD_ENOUGH
557 
	    aa = enmten_BESS; // 8.9e-308 (for R in ./bessel.h)
 566 
	    aa = enmten_BESS; // 8.9e-308 (for R in ./bessel.h)
558 
	    if (sum > 1.)
 567 
	    if (sum > 1.)
559 
		aa *= sum;
 568 
		aa *= sum;
- 
 
 569 
#endif
560 
	    for (n = 1; n <= *nb; ++n) {
 570 
	    for (n = 1; n <= *nb; ++n) {
- 
 
 571 
#ifdef UNDERFLOW_NOT_GOOD_ENOUGH
561 
		if (fabs(b[n]) < aa)
 572 
		if (fabs(b[n]) < aa)
562 
		    b[n] = 0.;
 573 
		    b[n] = 0.;
563 
		else
 574 
		else
- 
 
 575 
#endif
564 
		    b[n] /= sum;
 576 
		    b[n] /= sum;
565 
	    }
 577 
	    }
566 
	}
 578 
	}
567 
 
 579 
 
568 
    }
 580 
    }