X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=lib%2Fgsl-extras%2Fbetadistinv.c;h=490914056d691e7fc7dd06dca94d299bde8db80f;hb=4fb6f8878aa476e411f7d490a514e0892d3e51e2;hp=0b721cac759c1158601a872d35ab1ca6cbe09761;hpb=4f71ffe74a94ea89787705f1ba12a3dd810adda5;p=pspp-builds.git

diff --git a/lib/gsl-extras/betadistinv.c b/lib/gsl-extras/betadistinv.c
index 0b721cac..49091405 100644
--- a/lib/gsl-extras/betadistinv.c
+++ b/lib/gsl-extras/betadistinv.c
@@ -15,7 +15,7 @@
  *
  * You should have received a copy of the GNU General Public License
  * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  */
 
 /*
@@ -32,7 +32,7 @@
  * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8,
  * August 1968, pages 1264-1273.
  */
-#include <config.h>
+
 #include <math.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_errno.h>
@@ -41,7 +41,7 @@
 #include <gsl/gsl_randist.h>
 #include "gsl-extras.h"
 
-#define BETAINV_INIT_ERR .01
+#define BETAINV_INIT_ERR .001
 #define BETADISTINV_N_TERMS 3
 #define BETADISTINV_MAXITER 20
 
@@ -390,7 +390,7 @@ gslextras_cdf_beta_Pinv ( const double p, const double a, const double b)
   err = beta_result - p;
   abserr = fabs(err);
   relerr = abserr / p;
-  while ( relerr > BETAINV_INIT_ERR && n_iter < 100)
+  while ( relerr > BETAINV_INIT_ERR)
     {
       tmp = new_guess_P ( state, lower, upper, 
 			  p, a, b);
@@ -431,6 +431,19 @@ gslextras_cdf_beta_Pinv ( const double p, const double a, const double b)
 	  result = state;
 	  min_err = relerr;
 	}
+      else
+	{
+	  /*
+	   * Lagrange polynomial failed to reduce the error.
+	   * This will happen with a very skewed beta density. 
+	   * Undo previous steps.
+	   */
+	  state = result;
+	  beta_result = gsl_cdf_beta_P(state,a,b);
+	  err = p - beta_result;
+	  abserr = fabs(err);
+	  relerr = abserr / p;
+	}
     }
 
   n_iter = 0;
@@ -581,7 +594,7 @@ gslextras_cdf_beta_Qinv (double q, double a, double b)
   err = beta_result - q;
   abserr = fabs(err);
   relerr = abserr / q;
-  while ( relerr > BETAINV_INIT_ERR && n_iter < 100)
+  while ( relerr > BETAINV_INIT_ERR)
     {
       n_iter++;
       tmp = new_guess_Q ( state, lower, upper, 
@@ -622,6 +635,19 @@ gslextras_cdf_beta_Qinv (double q, double a, double b)
 	  result = state;
 	  min_err = relerr;
 	}
+      else
+	{
+	  /*
+	   * Lagrange polynomial failed to reduce the error.
+	   * This will happen with a very skewed beta density. 
+	   * Undo previous steps.
+	   */
+	  state = result;
+	  beta_result = gsl_cdf_beta_P(state,a,b);
+	  err = q - beta_result;
+	  abserr = fabs(err);
+	  relerr = abserr / q;
+	}
     }
 
   /*