#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double bisect2(double hi,double ga) {
	double mulower = 0.0;
	double muupper = 1.0;
	double mu, est, err;
	int i = 1;
	while(true) {
		mu = (mulower+muupper)/2;
		if((muupper-mulower) < 0.0001) {
			break;
		}
		est = Rf_pbeta(hi,mu*ga,(1-mu)*ga,true,false);
		if(est <= 0.995) {
			muupper = mu;
			mulower = mulower;
		} else {
			mulower = mu;
			muupper = muupper;
		}
		err = fabs(est-0.995);
		if(err < 0.0001) {
			break;
		}
		i += 1;
		if(i == 50000) {
			Rcout << "bisect2(" << hi << "," << ga << ") taking way too long\n";
		}
	}
	return mu;
}

// [[Rcpp::export]]
NumericVector bisect3(double lo, double hi) {
	double gaupper = 1.0; 
	double galower = 0.0;
	double ga, mu, est, err;
	int i = 1;
	NumericVector out(2);
	while(true) {
		mu = bisect2(hi,gaupper);
		est = Rf_qbeta(0.005,mu*gaupper,(1-mu)*gaupper,true,false);
		if(est < lo) {
			gaupper = 2*gaupper;
		} else {
			break;
		}
	}
	while(true) {
		ga = (gaupper+galower)/2;
		if((gaupper-galower) < 0.0001) {
			break;
		}
		mu = bisect2(hi,ga);
		est = Rf_qbeta(0.005,mu*ga,(1-mu)*ga,true,false);
		if(est < lo) {
			galower = ga;
			gaupper = gaupper;
		} else {
			galower = galower;
			gaupper = ga;
		}
		err = fabs(est-lo);
		if(err < 0.0001) {
			break;
		}
		i += 1;
		if(i == 50000) {
			Rcout << "bisect3(" << lo << "," << hi << ") taking way too long\n";
		}
	}
	out[0] = mu*ga;
	out[1] = (1-mu)*ga;
	return out;
}

// [[Rcpp::export]]
NumericMatrix nesurb(NumericVector lo, NumericVector hi) {
	int n = lo.size();
	NumericVector ans(2);
	NumericMatrix shape(n,2);
	for(int i = 0; i < n; i++) {
		ans = bisect3(lo[i],hi[i]);
		shape(i,0) = ans[0];
		shape(i,1) = ans[1];
	}
	return shape;
}