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

diff --git a/lib/gsl-extras/betadistinv.c b/lib/gsl-extras/betadistinv.c
index a0838d27..49091405 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
@@ -14,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.
  */
 
 /*
@@ -31,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>
@@ -40,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
 
@@ -389,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);
@@ -430,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;
@@ -580,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, 
@@ -621,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;
+	}
     }
 
   /*