The R Project SVN R

Rev

Rev 60118 | Rev 64656 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 60118 Rev 60268
Line 17... Line 17...
17 
 */
 17 
 */
18 
 
 18 
 
19 
#include "nmath.h"
 19 
#include "nmath.h"
20 
#include "dpq.h"
 20 
#include "dpq.h"
21 
 
 21 
 
22 
long double attribute_hidden
 22 
LDOUBLE attribute_hidden
23 
pnbeta_raw(double x, double o_x, double a, double b, double ncp)
 23 
pnbeta_raw(double x, double o_x, double a, double b, double ncp)
24 
{
 24 
{
25 
    /* o_x  == 1 - x  but maybe more accurate */
 25 
    /* o_x  == 1 - x  but maybe more accurate */
26 
 
 26 
 
27 
    /* change errmax and itrmax if desired;
 27 
    /* change errmax and itrmax if desired;
Line 31... Line 31...
31 
				     see PR#11277 */
 31 
				     see PR#11277 */
32 
 
 32 
 
33 
    double a0, lbeta, c, errbd, x0, temp, tmp_c;
 33 
    double a0, lbeta, c, errbd, x0, temp, tmp_c;
34 
    int j, ierr;
 34 
    int j, ierr;
35 
 
 35 
 
36 
    long double ans, ax, gx, q, sumq;
 36 
    LDOUBLE ans, ax, gx, q, sumq;
37 
 
 37 
 
38 
    if (ncp < 0. || a <= 0. || b <= 0.) ML_ERR_return_NAN;
 38 
    if (ncp < 0. || a <= 0. || b <= 0.) ML_ERR_return_NAN;
39 
 
 39 
 
40 
    if(x < 0. || o_x > 1. || (x == 0. && o_x == 1.)) return 0.;
 40 
    if(x < 0. || o_x > 1. || (x == 0. && o_x == 1.)) return 0.;
41 
    if(x > 1. || o_x < 0. || (x == 1. && o_x == 0.)) return 1.;
 41 
    if(x > 1. || o_x < 0. || (x == 1. && o_x == 0.)) return 1.;
Line 85... Line 85...
85 
double attribute_hidden
 85 
double attribute_hidden
86 
pnbeta2(double x, double o_x, double a, double b, double ncp,
 86 
pnbeta2(double x, double o_x, double a, double b, double ncp,
87 
	/* o_x  == 1 - x  but maybe more accurate */
 87 
	/* o_x  == 1 - x  but maybe more accurate */
88 
	int lower_tail, int log_p)
 88 
	int lower_tail, int log_p)
89 
{
 89 
{
90 
    long double ans = pnbeta_raw(x, o_x, a,b, ncp);
 90 
    LDOUBLE ans = pnbeta_raw(x, o_x, a,b, ncp);
- 
 
 91 
 
91 
 
 92 
 
92 
    /* return R_DT_val(ans), but we want to warn about cancellation here */
 93 
    /* return R_DT_val(ans), but we want to warn about cancellation here */
- 
 
 94 
    if (lower_tail)
- 
 
 95 
#ifdef HAVE_LONG_DOUBLE
93 
    if (lower_tail) return (double) (log_p ? logl(ans) : ans);
 96 
	return (double) (log_p ? logl(ans) : ans);
- 
 
 97 
#else
- 
 
 98 
	return log_p ? log(ans) : ans;
- 
 
 99 
#endif
94 
    else {
 100 
    else {
95 
	if (ans > 1. - 1e-10) ML_ERROR(ME_PRECISION, "pnbeta");
 101 
	if (ans > 1. - 1e-10) ML_ERROR(ME_PRECISION, "pnbeta");
96 
	if (ans > 1.0) ans = 1.0;  /* Precaution */
 102 
	if (ans > 1.0) ans = 1.0;  /* Precaution */
97 
#ifdef HAVE_LOG1PL
 103 
#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_LOG1PL)
98 
	return (double) (log_p ? log1pl(-ans) : (1. - ans));
 104 
	return (double) (log_p ? log1pl(-ans) : (1. - ans));
99 
#else
 105 
#else
100 
	/* include standalone case */
 106 
	/* include standalone case */
101 
	return (double) (log_p ? log1p((double)-ans) : (1. - ans));
 107 
	return (double) (log_p ? log1p((double)-ans) : (1. - ans));
102 
#endif
 108 
#endif