@@ -158,11 +158,11 @@
 #define Q4              xxx[10]
 #define Q5              xxx[11]
 
 double
 expm1(double x) {
-        double y, hi, lo, c, t, e, hxs, hfx, r1;
+        double y, hi, lo, c = 0.0L, t, e, hxs, hfx, r1;
         int k, xsb;
         unsigned hx;
 
         hx = ((unsigned *) &x)[HIWORD];         /* high word of x */
         xsb = hx & 0x80000000;                  /* sign bit of x */

@@ -170,13 +170,14 @@
                 y = x;
         else
                 y = -x;                         /* y = |x| */
         hx &= 0x7fffffff;                       /* high word of |x| */
 
-        /* filter out huge and non-finite arugment */
-        if (hx >= 0x4043687A) {                 /* if |x|>=56*ln2 */
-                if (hx >= 0x40862E42) {         /* if |x|>=709.78... */
+        /* filter out huge and non-finite argument */
+        /* for example exp(38)-1 is approximately 3.1855932e+16 */
+        if (hx >= 0x4043687A) {                 /* if |x|>=56*ln2  (~38.8162...)*/
+                if (hx >= 0x40862E42) {         /* if |x|>=709.78... -> inf */
                         if (hx >= 0x7ff00000) {
                                 if (((hx & 0xfffff) | ((int *) &x)[LOWORD])
                                         != 0)
                                         return x * x;   /* + -> * for Cheetah */
                                 else

@@ -192,55 +193,56 @@
         }
 
         /* argument reduction */
         if (hx > 0x3fd62e42) {                  /* if  |x| > 0.5 ln2 */
                 if (hx < 0x3FF0A2B2) {          /* and |x| < 1.5 ln2 */
-                        if (xsb == 0) {
+                        if (xsb == 0) {         /* positive number */
                                 hi = x - ln2_hi;
                                 lo = ln2_lo;
                                 k = 1;
                         }
-                        else {
+                        else { /* negative number */
                                 hi = x + ln2_hi;
                                 lo = -ln2_lo;
                                 k = -1;
                         }
                 }
-                else {
+                else {  /* |x| > 1.5 ln2 */
                         k = (int) (invln2 * x + (xsb == 0 ? 0.5 : -0.5));
                         t = k;
                         hi = x - t * ln2_hi;    /* t*ln2_hi is exact here */
                         lo = t * ln2_lo;
                 }
                 x = hi - lo;
-                c = (hi - x) - lo;
+                c = (hi - x) - lo; /* still at |x| > 0.5 ln2 */
         }
         else if (hx < 0x3c900000) {             /* when |x|<2**-54, return x */
                 t = huge + x;           /* return x w/inexact when x != 0 */
                 return x - (t - (huge + x));
         }
-        else
+        else    /* |x| <= 0.5 ln2 */
                 k = 0;
 
         /* x is now in primary range */
         hfx = 0.5 * x;
         hxs = x * hfx;
         r1 = one + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
         t = 3.0 - r1 * hfx;
         e = hxs * ((r1 - t) / (6.0 - x * t));
-        if (k == 0)
-                return x - (x * e - hxs);       /* c is 0 */
-        else {
+        if (k == 0) /* |x| <= 0.5 ln2 */
+                return x - (x * e - hxs);
+        else {      /* |x| > 0.5 ln2 */
                 e = (x * (e - c) - c);
                 e -= hxs;
                 if (k == -1)
                         return 0.5 * (x - e) - 0.5;
-                if (k == 1)
+                if (k == 1) {
                         if (x < -0.25)
                                 return -2.0 * (e - (x + 0.5));
                         else
                                 return one + 2.0 * (x - e);
+                }
                 if (k <= -2 || k > 56) {        /* suffice to return exp(x)-1 */
                         y = one - (e - x);
                         ((int *) &y)[HIWORD] += k << 20;
                         return y - one;
                 }