X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;ds=sidebyside;f=lib%2Fgsl-extras%2Fbetadistinv.c;h=76feb0d8d2ca8f82b797685c0ea9269672be210c;hb=084280e4bfd6b94c26320e72553a9a38201610a3;hp=a0838d27eb744b5ce54bddc051c294b264e55741;hpb=e6645d1e3d026dfdcadd1e63a6d6445823c92fa2;p=pspp

diff --git a/lib/gsl-extras/betadistinv.c b/lib/gsl-extras/betadistinv.c
index a0838d27eb..76feb0d8d2 100644
--- a/lib/gsl-extras/betadistinv.c
+++ b/lib/gsl-extras/betadistinv.c
@@ -1,6 +1,7 @@
 /* cdf/betadistinv.c
  *
- * Copyright (C) 2004 Jason H. Stover.
+ * Copyright (C) 2004 Free Software Foundation, Inc.
+ * Written by Jason H. Stover.
  *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
@@ -31,6 +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>
@@ -40,7 +42,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
 
@@ -389,7 +391,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);
@@ -430,6 +432,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;
@@ -580,7 +595,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, 
@@ -621,6 +636,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;
+	}
     }
 
   /*