[Scipy-svn] r4180 - trunk/scipy/special/c_misc

scipy-svn@scip... scipy-svn@scip...
Sat Apr 26 18:33:36 CDT 2008


Author: cookedm
Date: 2008-04-26 18:33:30 -0500 (Sat, 26 Apr 2008)
New Revision: 4180

Modified:
   trunk/scipy/special/c_misc/fsolve.c
   trunk/scipy/special/c_misc/misc.h
Log:
scipy.special: Fix for #655: gammaincinv never ends
This shouldn't happen, as I constructed the root-finding routine to guarantee
convergence. I added a maximum number of iterations nevertheless.


Modified: trunk/scipy/special/c_misc/fsolve.c
===================================================================
--- trunk/scipy/special/c_misc/fsolve.c	2008-04-26 03:08:04 UTC (rev 4179)
+++ trunk/scipy/special/c_misc/fsolve.c	2008-04-26 23:33:30 UTC (rev 4180)
@@ -1,6 +1,10 @@
 #include "misc.h"
 #include <math.h>
 
+#define MAX_ITERATIONS              100
+#define FP_CMP_WITH_BISECT_NITER    4
+#define FP_CMP_WITH_BISECT_WIDTH    4.0
+
 static inline double
 max(double a, double b)
 {
@@ -54,6 +58,7 @@
     double x3, f3;
     double w, last_bisect_width;
     double tol;
+    int niter;
 
     if (f1*f2 >= 0.0) {
         return FSOLVE_NOT_BRACKET;
@@ -65,7 +70,7 @@
     }
     w = fabs(x2 - x1);
     last_bisect_width = w;
-    while (1) {
+    for (niter=0; niter < MAX_ITERATIONS; niter++) {
         switch (state) {
         case bisect: {
             x3 = 0.5 * (x1 + x2);
@@ -130,8 +135,8 @@
                sure that we always end up decreasing the interval width
                with a bisection.
              */
-            if (n_falsep > 4) {
-                if (w*4 > last_bisect_width) {
+            if (n_falsep > FP_CMP_WITH_BISECT_NITER) {
+                if (w*FP_CMP_WITH_BISECT_WIDTH > last_bisect_width) {
                     state = bisect;
                 }
                 n_falsep = 0;
@@ -150,6 +155,9 @@
             goto finish;
         }
     }
+    r = FSOLVE_MAX_ITERATIONS;
+    *best_x = x3; *best_f = f3;
+    goto finish;
 exact_soln:
     *best_x = x3; *best_f = 0.0;
     r = FSOLVE_EXACT;

Modified: trunk/scipy/special/c_misc/misc.h
===================================================================
--- trunk/scipy/special/c_misc/misc.h	2008-04-26 03:08:04 UTC (rev 4179)
+++ trunk/scipy/special/c_misc/misc.h	2008-04-26 23:33:30 UTC (rev 4180)
@@ -9,6 +9,8 @@
   FSOLVE_CONVERGED,
   /* Not a bracket */
   FSOLVE_NOT_BRACKET,
+  /* Root-finding didn't converge in a set number of iterations. */
+  FSOLVE_MAX_ITERATIONS
 } fsolve_result_t;
 
 typedef double (*objective_function)(double, void *);



More information about the Scipy-svn mailing list