4848#define  n_NEWTON_FREE  4
4949//                   --- based on some testing; had = 10 
5050
51- #define  MLOGICAL_NA  -1
52- // an "NA_LOGICAL" substitute for Mathlib {only used here, for now} 
53- 
54- //attribute_hidden 
5551static  void 
5652qbeta_raw (double  alpha , double  p , double  q ,
5753	  Rboolean  lower_tail , Rboolean  log_p ,
58- 	  int   swap_01 ,  double  log_q_cut , int  n_N , double *  qb );
54+ 	  double  log_q_cut , int  n_N , double *  qb );
5955
6056/* This is the public interface in Rmath.h, so has to use 'int' */ 
6157double  qbeta (double  alpha , double  p , double  q , int  lower_tail , int  log_p )
@@ -71,8 +67,8 @@ double qbeta(double alpha, double p, double q, int lower_tail, int log_p)
7167
7268    double  qbet [2 ];// = { qbeta(), 1 - qbeta() } 
7369    qbeta_raw (alpha , p , q , lower_tail , log_p ,
74- 	      // swap_01 ,   log_q_cut ,      n_N 
75- 	      MLOGICAL_NA ,  USE_LOG_X_CUTOFF , n_NEWTON_FREE , qbet );
70+ 	      // log_q_cut ,      n_N 
71+ 	      USE_LOG_X_CUTOFF , n_NEWTON_FREE , qbet );
7672    return  qbet [0 ];
7773}
7874
@@ -111,21 +107,17 @@ static const double
111107#define  const4  0.04481
112108
113109// Returns both qbeta() and its "mirror" 1-qbeta(). Useful notably when qbeta() ~= 1 
114- // This was hidden, but it is only used in this file so could be static 
115- // and be simplifed as swap_01 is always NA. 
116110// attribute_hidden void 
117111static  void 
118112qbeta_raw (double  alpha , double  p , double  q ,
119113	  Rboolean  lower_tail , Rboolean  log_p ,
120- 	  int  swap_01 , // {TRUE, NA, FALSE}: if NA, algorithm decides swap_tail 
121114	  double  log_q_cut , /* if == Inf: return log(qbeta(..)); 
122115			       otherwise, if finite: the bound for 
123116			       switching to log(x)-scale; see use_log_x */ 
124117	  int  n_N ,  // number of "unconstrained" Newton steps before switching to constrained 
125118	  double  * qb ) // = qb[0:1] = { qbeta(), 1 - qbeta() } 
126119{
127120    Rboolean 
128- 	swap_choose  =  (swap_01  ==  MLOGICAL_NA ),
129121	swap_tail ,
130122	log_ , give_log_q  =  (log_q_cut  ==  ML_POSINF ),
131123	use_log_x  =  give_log_q , // or u < log_q_cut  below 
@@ -195,21 +187,20 @@ qbeta_raw(double alpha, double p, double q,
195187    // Conceptually,  0 < p_ < 1  (but can be 0 or 1 because of cancellation!) 
196188    logbeta  =  lbeta (p , q );
197189
198-     // this is only ever called with swap_choose = TRUE 
199-     swap_tail  =  (swap_choose ) ? (p_  >  0.5 ) : (Rboolean ) swap_01 ;
190+     swap_tail  =  (p_  >  0.5 );
200191
201192    R_ifDEBUG_printf (
202- 	"qbeta(%g, %g, %g, lower_t=%d, log_p=%d, swap_01=%d,  log_q_cut=%g, n_N=%d):%s\n" 
193+ 	"qbeta(%g, %g, %g, lower_t=%d, log_p=%d, log_q_cut=%g, n_N=%d):%s\n" 
203194	"  swap_tail=%s :" ,
204- 	alpha , p ,q , lower_tail , log_p , swap_01 ,  log_q_cut , n_N ,
195+ 	alpha , p ,q , lower_tail , log_p , log_q_cut , n_N ,
205196	(log_p  &&  (p_  ==  0.  ||  p_  ==  1. )) ? (p_ == 0. ?" p_=0" :" p_=1" ) : "" ,
206197	(swap_tail  ? "TRUE" : "F" ));
207198
208199    int  n_maybe_swaps  =  0 ;
209200maybe_swap :
210-     // change tail; default (swap_01 = NA):  afterwards 0 < a <= 1/2 
201+     // change tail; afterwards 0 < a <= 1/2 
211202    if (swap_tail ) { /* change tail, swap copies of {p,q}:  p <-> q :*/ 
212- 	a  =  R_DT_CIv (alpha ); // = 1 - p_ , is < 1/2 if(swap_choose)  
203+ 	a  =  R_DT_CIv (alpha ); // = 1 - p_ , is < 1/2 
213204	/* la := log(a), but without numerical cancellation: */ 
214205	la  =  R_DT_Clog (alpha );
215206	pp  =  q ; qq  =  p ;
@@ -346,7 +337,7 @@ qbeta_raw(double alpha, double p, double q,
346337    }
347338
348339    // Problem: If initial u is completely wrong, we make a wrong decision here 
349-     if (swap_choose   &&    //   vvvv/ why -exp(*)? u  on log-x scale! Swapping can be very good, but needs smart if(..) 
340+     if (
350341       (( swap_tail  &&  u  >= - exp (  log_q_cut )) ||  // ==> "swap back" 
351342	(!swap_tail  &&  u  >= - exp (4 * log_q_cut ) &&  pp  / qq  <  1000. ) // ==> "swap now" 
352343	   )) {
0 commit comments