The R Project SVN R

Rev

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

Rev 24735 Rev 24769
Line 84... Line 84...
84 
	pp = q; qq = p; swap_tail = 1;
 84 
	pp = q; qq = p; swap_tail = 1;
85 
    }
 85 
    }
86 
 
 86 
 
87 
    /* calculate the initial approximation */
 87 
    /* calculate the initial approximation */
88 
 
 88 
 
89 
    r = sqrt(-log(a * a));
 89 
    r = sqrt(-2 * log(a));
90 
    y = r - (const1 + const2 * r) / (1. + (const3 + const4 * r) * r);
 90 
    y = r - (const1 + const2 * r) / (1. + (const3 + const4 * r) * r);
91 
    if (pp > 1 && qq > 1) {
 91 
    if (pp > 1 && qq > 1) {
92 
	r = (y * y - 3.) / 6.;
 92 
	r = (y * y - 3.) / 6.;
93 
	s = 1. / (pp + pp - 1.);
 93 
	s = 1. / (pp + pp - 1.);
94 
	t = 1. / (qq + qq - 1.);
 94 
	t = 1. / (qq + qq - 1.);
Line 98... Line 98...
98 
    } else {
 98 
    } else {
99 
	r = qq + qq;
 99 
	r = qq + qq;
100 
	t = 1. / (9. * qq);
 100 
	t = 1. / (9. * qq);
101 
	t = r * pow(1. - t + y * sqrt(t), 3.0);
 101 
	t = r * pow(1. - t + y * sqrt(t), 3.0);
102 
	if (t <= 0.)
 102 
	if (t <= 0.)
103 
	    xinbta = 1. - exp((log((1. - a) * qq) + logbeta) / qq);
 103 
	    xinbta = 1. - exp((log1p(-a)+ log(qq) + logbeta) / qq);
104 
	else {
 104 
	else {
105 
	    t = (4. * pp + r - 2.) / t;
 105 
	    t = (4. * pp + r - 2.) / t;
106 
	    if (t <= 1.)
 106 
	    if (t <= 1.)
107 
		xinbta = exp((log(a * pp) + logbeta) / pp);
 107 
		xinbta = exp((log(a * pp) + logbeta) / pp);
108 
	    else
 108 
	    else