The R Project SVN R

Rev

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

Rev 77685 Rev 89291
Line 57... Line 57...
57 
	/* Using Abramowitz & Stegun  9.1.2
 57 
	/* Using Abramowitz & Stegun  9.1.2
58 
	 * this may not be quite optimal (CPU and accuracy wise) */
 58 
	 * this may not be quite optimal (CPU and accuracy wise) */
59 
	return(((alpha - na == 0.5) ? 0 : bessel_y(x, -alpha) * cospi(alpha)) -
 59 
	return(((alpha - na == 0.5) ? 0 : bessel_y(x, -alpha) * cospi(alpha)) -
60 
	       ((alpha      == na ) ? 0 : bessel_j(x, -alpha) * sinpi(alpha)));
 60 
	       ((alpha      == na ) ? 0 : bessel_j(x, -alpha) * sinpi(alpha)));
61 
    }
 61 
    }
62 
    else if (alpha > 1e7) {
 62 
    else if (alpha > 1e7) { // NB: same bound 'besselJY_max_nu' in math_2b() and ./bessel_j.c
63 
	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
 63 
	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
64 
			alpha);
 64 
			alpha);
65 
	return ML_NAN;
 65 
	return ML_NAN;
66 
    }
 66 
    }
67 
    nb = 1+ (int)na;/* nb-1 <= alpha < nb */
 67 
    nb = 1+ (int)na;/* nb-1 <= alpha < nb */
Line 97... Line 97...
97 
    vmaxset(vmax);
 97 
    vmaxset(vmax);
98 
#endif
 98 
#endif
99 
    return x;
 99 
    return x;
100 
}
 100 
}
101 
 
 101 
 
102 
/* Called from R: modified version of bessel_y(), accepting a work array
 102 
/* Called from R via math_2b() in ../main/arithmetic.c:
103 
 * instead of allocating one. */
 103 
 * modified version of bessel_y(), accepting a work array instead of allocating one. */
104 
double bessel_y_ex(double x, double alpha, double *by)
 104 
double bessel_y_ex(double x, double alpha, double *by)
105 
{
 105 
{
106 
    int nb, ncalc;
 106 
    int nb, ncalc;
107 
    double na;
 107 
    double na;
108 
 
 108 
 
Line 119... Line 119...
119 
	/* Using Abramowitz & Stegun  9.1.2
 119 
	/* Using Abramowitz & Stegun  9.1.2
120 
	 * this may not be quite optimal (CPU and accuracy wise) */
 120 
	 * this may not be quite optimal (CPU and accuracy wise) */
121 
	return(((alpha - na == 0.5) ? 0 : bessel_y_ex(x, -alpha, by) * cospi(alpha)) -
 121 
	return(((alpha - na == 0.5) ? 0 : bessel_y_ex(x, -alpha, by) * cospi(alpha)) -
122 
	       ((alpha      == na ) ? 0 : bessel_j_ex(x, -alpha, by) * sinpi(alpha)));
 122 
	       ((alpha      == na ) ? 0 : bessel_j_ex(x, -alpha, by) * sinpi(alpha)));
123 
    }
 123 
    }
124 
    else if (alpha > 1e7) {
 124 
    else if (alpha > 1e7) { // NB: same bound 'besselJY_max_nu' in math_2b() and ./bessel_j.c
125 
	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
 125 
	MATHLIB_WARNING(_("besselY(x, nu): nu=%g too large for bessel_y() algorithm"),
126 
			alpha);
 126 
			alpha);
127 
	return ML_NAN;
 127 
	return ML_NAN;
128 
    }
 128 
    }
129 
    nb = 1+ (int)na;/* nb-1 <= alpha < nb */
 129 
    nb = 1+ (int)na;/* nb-1 <= alpha < nb */