The R Project SVN R

Rev

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

Rev 68947 Rev 77685
Line 59... Line 59...
59 
    /* note - itrmax and errmax may be changed to suit one's needs. */
 59 
    /* note - itrmax and errmax may be changed to suit one's needs. */
60 
 
 60 
 
61 
    const int itrmax = 1000;
 61 
    const int itrmax = 1000;
62 
    const static double errmax = 1.e-12;
 62 
    const static double errmax = 1.e-12;
63 
 
 63 
 
64 
    if (df <= 0.0) ML_ERR_return_NAN;
 64 
    if (df <= 0.0) ML_WARN_return_NAN;
65 
    if(ncp == 0.0) return pt(t, df, lower_tail, log_p);
 65 
    if(ncp == 0.0) return pt(t, df, lower_tail, log_p);
66 
 
 66 
 
67 
    if(!R_FINITE(t))
 67 
    if(!R_FINITE(t))
68 
	return (t < 0) ? R_DT_0 : R_DT_1;
 68 
	return (t < 0) ? R_DT_0 : R_DT_1;
69 
    if (t >= 0.) {
 69 
    if (t >= 0.) {
Line 103... Line 103...
103 
	REprintf("\t p=%10Lg\n",p);
 103 
	REprintf("\t p=%10Lg\n",p);
104 
#endif
 104 
#endif
105 
	if(p == 0.) { /* underflow! */
 105 
	if(p == 0.) { /* underflow! */
106 
 
 106 
 
107 
	    /*========== really use an other algorithm for this case !!! */
 107 
	    /*========== really use an other algorithm for this case !!! */
108 
	    ML_ERROR(ME_UNDERFLOW, "pnt");
 108 
	    ML_WARNING(ME_UNDERFLOW, "pnt");
109 
	    ML_ERROR(ME_RANGE, "pnt"); /* |ncp| too large */
 109 
	    ML_WARNING(ME_RANGE, "pnt"); /* |ncp| too large */
110 
	    return R_DT_0;
 110 
	    return R_DT_0;
111 
	}
 111 
	}
112 
#ifdef DEBUG_pnt
 112 
#ifdef DEBUG_pnt
113 
	REprintf("it  1e5*(godd,   geven)|          p           q           s"
 113 
	REprintf("it  1e5*(godd,   geven)|          p           q           s"
114 
	       /* 1.3 1..4..7.9 1..4..7.9|1..4..7.901 1..4..7.901 1..4..7.901 */
 114 
	       /* 1.3 1..4..7.9 1..4..7.9|1..4..7.901 1..4..7.901 1..4..7.901 */
Line 144... Line 144...
144 
	    q *= lambda / (2 * it + 1);
 144 
	    q *= lambda / (2 * it + 1);
145 
	    tnc += p * xodd + q * xeven;
 145 
	    tnc += p * xodd + q * xeven;
146 
	    s -= p;
 146 
	    s -= p;
147 
	    /* R 2.4.0 added test for rounding error here. */
 147 
	    /* R 2.4.0 added test for rounding error here. */
148 
	    if(s < -1.e-10) { /* happens e.g. for (t,df,ncp)=(40,10,38.5), after 799 it.*/
 148 
	    if(s < -1.e-10) { /* happens e.g. for (t,df,ncp)=(40,10,38.5), after 799 it.*/
149 
		ML_ERROR(ME_PRECISION, "pnt");
 149 
		ML_WARNING(ME_PRECISION, "pnt");
150 
#ifdef DEBUG_pnt
 150 
#ifdef DEBUG_pnt
151 
		REprintf("s = %#14.7Lg < 0 !!! ---> non-convergence!!\n", s);
 151 
		REprintf("s = %#14.7Lg < 0 !!! ---> non-convergence!!\n", s);
152 
#endif
 152 
#endif
153 
		goto finis;
 153 
		goto finis;
154 
	    }
 154 
	    }
Line 159... Line 159...
159 
		     it, 1e5*(double)godd, 1e5*(double)geven, p, q, s, tnc, errbd);
 159 
		     it, 1e5*(double)godd, 1e5*(double)geven, p, q, s, tnc, errbd);
160 
#endif
 160 
#endif
161 
	    if(fabs(errbd) < errmax) goto finis;/*convergence*/
 161 
	    if(fabs(errbd) < errmax) goto finis;/*convergence*/
162 
	}
 162 
	}
163 
	/* non-convergence:*/
 163 
	/* non-convergence:*/
164 
	ML_ERROR(ME_NOCONV, "pnt");
 164 
	ML_WARNING(ME_NOCONV, "pnt");
165 
    }
 165 
    }
166 
    else { /* x = t = 0 */
 166 
    else { /* x = t = 0 */
167 
	tnc = 0.;
 167 
	tnc = 0.;
168 
    }
 168 
    }
169 
 finis:
 169 
 finis:
170 
    tnc += pnorm(- del, 0., 1., /*lower*/TRUE, /*log_p*/FALSE);
 170 
    tnc += pnorm(- del, 0., 1., /*lower*/TRUE, /*log_p*/FALSE);
171 
 
 171 
 
172 
    lower_tail = lower_tail != negdel; /* xor */
 172 
    lower_tail = lower_tail != negdel; /* xor */
173 
    if(tnc > 1 - 1e-10 && lower_tail)
 173 
    if(tnc > 1 - 1e-10 && lower_tail)
174 
	ML_ERROR(ME_PRECISION, "pnt{final}");
 174 
	ML_WARNING(ME_PRECISION, "pnt{final}");
175 
 
 175 
 
176 
    return R_DT_val(fmin2((double)tnc, 1.) /* Precaution */);
 176 
    return R_DT_val(fmin2((double)tnc, 1.) /* Precaution */);
177 
}
 177 
}