* Bug 16209 fixed: grand(n,m,exp,1) may freeze Scilab with default generator

[scilab.git] / scilab / modules / randlib / src / c / sexpo.c

diff --git a/scilab/modules/randlib/src/c/sexpo.c b/scilab/modules/randlib/src/c/sexpo.c
index d95ba28..a385839 100644 (file)
--- a/scilab/modules/randlib/src/c/sexpo.c
+++ b/scilab/modules/randlib/src/c/sexpo.c
@@ -1,4 +1,3 @@
-#include "stack-c.h"
 #include <math.h>
 #include "grand.h"
 #include "core_math.h"
@@ -15,7 +14,7 @@ double C2F(sexpo)(void)
 **********************************************************************
 
 This source code was taken in the project "freemat"(BSD license)
-This source code was modified by Gaüzère Sabine according to the 
+This source code was modified by Gaüzère Sabine according to the
 modifications done by JJV
 
      FOR DETAILS SEE:
@@ -36,40 +35,48 @@ modifications done by JJV
      (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
 */
 {
-/* ** OLD VALUES
-static double q[8] = {
-    0.69314718055995, 0.93337368751905, 0.98887779618387, 0.99849592529150,
-    0.99982928110614, 0.99998331641007, 0.99999856914388, 0.99999989069256
-//0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,0.9999999
-};
-*/
-static double q[8] = {
- 0.69314718246459960938,
- 0.93337368965148925781,
- 0.98887777328491210938,
- 0.99849587678909301758,
- 0.99982929229736328125,
- 0.99998331069946289062,
- 0.99999862909317016602,
- 0.99999988079071044922
-};
-static int i;
-static double sexpo,a,u,ustar,umin;
+    /* ** OLD VALUES
+    static double q[8] = {
+        0.69314718055995, 0.93337368751905, 0.98887779618387, 0.99849592529150,
+        0.99982928110614, 0.99998331641007, 0.99999856914388, 0.99999989069256
+    //0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,0.9999999
+    };
+    */
+    static double q[8] =
+    {
+        0.69314718246459960938,
+        0.93337368965148925781,
+        0.98887777328491210938,
+        0.99849587678909301758,
+        0.99982929229736328125,
+        0.99998331069946289062,
+        0.99999862909317016602,
+        0.99999988079071044922
+    };
+    static int i;
+    static double sexpo, a, u, ustar, umin;
 
     a = 0.0;
     u = C2F(ranf)();
+    u = (u == 0.0) ? C2F(ranf)() : u;
     goto S30;
 S20:
     a += q[0];
 S30:
     u += u;
-//  JJV changed the following to reflect the true algorithm and
-//  JJV prevent unpredictable behavior if U is initially 0.5.
-//  IF (u.LE.1.0) GO TO 20
-    if(u < 1.0) goto S20;
+    //  JJV changed the following to reflect the true algorithm and
+    //  JJV prevent unpredictable behavior if U is initially 0.5.
+    //  IF (u.LE.1.0) GO TO 20
+    if (u < 1.0)
+    {
+        goto S20;
+    }
     u -= 1.0;
-    if(u > q[0]) goto S60;
-    sexpo = a+u;
+    if (u > q[0])
+    {
+        goto S60;
+    }
+    sexpo = a + u;
     return sexpo;
 S60:
     i = 1;
@@ -77,9 +84,15 @@ S60:
     umin = ustar;
 S70:
     ustar = C2F(ranf)();
-    if(ustar < umin) umin = ustar;
+    if (ustar < umin)
+    {
+        umin = ustar;
+    }
     i += 1;
-    if(u > q[i-1]) goto S70;
-    sexpo = a+umin*q[0];
+    if (u > q[i - 1])
+    {
+        goto S70;
+    }
+    sexpo = a + umin * q[0];
     return sexpo;
 }