EXITPOINT to delete malloced var in case of check failure

src/maths/cmaths/cmath1.c

@@ -16,6 +16,8 @@ Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
  *
  */
 
+#include <errno.h>
+
 #include "ngspice/ngspice.h"
 #include "ngspice/memory.h"
 #include "ngspice/cpdefs.h"
@@ -215,7 +217,7 @@ void *cx_conj(void *data, short int type, int length,
     }
 
     /* Else real, so just copy */
-    return memcpy(alloc_d(length), data, length * sizeof(double));
+    return memcpy(alloc_d(length), data, (unsigned int) length * sizeof(double));
 } /* end of function cx_conj */
 
 
@@ -241,21 +243,20 @@ cx_pos(void *data, short int type, int length, int *newlength, short int *newtyp
     return ((void *) d);
 }
 
-void *
-cx_db(void *data, short int type, int length, int *newlength, short int *newtype)
+void *cx_db(void *data, short int type, int length,
+        int *newlength, short int *newtype)
 {
+    int xrc = 0;
     double *d = alloc_d(length);
-    double *dd = (double *) data;
-    ngcomplex_t *cc = (ngcomplex_t *) data;
-    double tt;
-    int i;
 
     *newlength = length;
     *newtype = VF_REAL;
-    if (type == VF_COMPLEX)
+    if (type == VF_COMPLEX) {
+        ngcomplex_t *cc = (ngcomplex_t *) data;
+        int i;
         for (i = 0; i < length; i++) {
-            tt = cmag(cc[i]);
-            rcheckn(tt > 0, "db", d);
+            const double tt = cmag(cc[i]);
+            rcheck(tt > 0, "db");
             /*
                 if (tt == 0.0)
                     d[i] = 20.0 * - log(HUGE);
@@ -263,29 +264,44 @@ cx_db(void *data, short int type, int length, int *newlength, short int *newtype
             */
             d[i] = 20.0 * log10(tt);
         }
-    else
+    }
+    else {
+        double *dd = (double *) data;
+        int i;
         for (i = 0; i < length; i++) {
-            rcheckn(dd[i] > 0, "db", d);
+            const double tt = dd[i];
+            rcheck(tt > 0, "db");
             /*
                 if (dd[i] == 0.0)
                     d[i] = 20.0 * - log(HUGE);
                 else
             */
-            d[i] = 20.0 * log10(dd[i]);
+            d[i] = 20.0 * log10(tt);
         }
-    return ((void *) d);
     }
 
-void *
-cx_log10(void *data, short int type, int length, int *newlength, short int *newtype)
+EXITPOINT:
+    if (xrc != 0) {
+        txfree(d);
+        d = (double *) NULL;
+    }
+    return ((void *) d);
+} /* end of function cx_db */
+
+
+
+void *cx_log10(void *data, short int type, int length,
+        int *newlength, short int *newtype)
 {
-    *newlength = length;
+    int xrc = 0;
+    void *rv;
+
     if (type == VF_COMPLEX) {
         ngcomplex_t *c;
         ngcomplex_t *cc = (ngcomplex_t *) data;
         int i;
 
-        c = alloc_c(length);
+        rv = c = alloc_c(length);
         *newtype = VF_COMPLEX;
         for (i = 0; i < length; i++) {
             double td;
@@ -294,79 +310,107 @@ cx_log10(void *data, short int type, int length, int *newlength, short int *newt
             /* Perhaps we should trap when td = 0.0, but Ken wants
              * this to be possible...
              */
-            rcheckn(td >= 0, "log10", c);
+            rcheck(td >= 0, "log10");
             if (td == 0.0) {
                 realpart(c[i]) = - log10(HUGE);
                 imagpart(c[i]) = 0.0;
-            } else {
+            }
+            else {
                 realpart(c[i]) = log10(td);
-                imagpart(c[i]) = atan2(imagpart(cc[i]),
-                                        realpart(cc[i]));
+                imagpart(c[i]) = atan2(imagpart(cc[i]), realpart(cc[i]));
             }
         }
-        return ((void *) c);
-    } else {
+    }
+    else {
         double *d;
         double *dd = (double *) data;
         int i;
 
-        d = alloc_d(length);
+        rv = d = alloc_d(length);
         *newtype = VF_REAL;
         for (i = 0; i < length; i++) {
-            rcheckn(dd[i] >= 0, "log10", d);
-            if (dd[i] == 0.0)
+            rcheck(dd[i] >= 0, "log10");
+            if (dd[i] == 0.0) {
                 d[i] = - log10(HUGE);
-            else
+            }
+            else {
                 d[i] = log10(dd[i]);
             }
-        return ((void *) d);
         }
     }
 
-void *
-cx_log(void *data, short int type, int length, int *newlength, short int *newtype)
-{
     *newlength = length;
+
+EXITPOINT:
+    if (xrc != 0) { /* Free resources on error */
+        txfree(rv);
+        rv = NULL;
+    }
+
+    return rv;
+} /* end of function cx_log10 */
+
+
+
+void *cx_log(void *data, short int type, int length,
+        int *newlength, short int *newtype)
+{
+    int xrc = 0;
+    void *rv;
+
     if (type == VF_COMPLEX) {
         ngcomplex_t *c;
         ngcomplex_t *cc = (ngcomplex_t *) data;
-        int i;
 
-        c = alloc_c(length);
+        rv = c = alloc_c(length);
         *newtype = VF_COMPLEX;
+
+        int i;
         for (i = 0; i < length; i++) {
             double td;
 
             td = cmag(cc[i]);
-            rcheckn(td >= 0, "log", c);
+            rcheck(td >= 0, "log");
             if (td == 0.0) {
                 realpart(c[i]) = - log(HUGE);
                 imagpart(c[i]) = 0.0;
-            } else {
+            }
+            else {
                 realpart(c[i]) = log(td);
-                imagpart(c[i]) = atan2(imagpart(cc[i]),
-                                        realpart(cc[i]));
+                imagpart(c[i]) = atan2(imagpart(cc[i]), realpart(cc[i]));
             }
         }
-        return ((void *) c);
-    } else {
+    }
+    else {
         double *d;
         double *dd = (double *) data;
-        int i;
 
-        d = alloc_d(length);
+        rv = d = alloc_d(length);
         *newtype = VF_REAL;
+
+        int i;
         for (i = 0; i < length; i++) {
-            rcheckn(dd[i] >= 0, "log", d);
+            rcheck(dd[i] >= 0, "log");
             if (dd[i] == 0.0)
                 d[i] = - log(HUGE);
             else
                 d[i] = log(dd[i]);
         }
-        return ((void *) d);
     }
+
+    *newlength = length;
+
+EXITPOINT:
+    if (xrc != 0) { /* Free resources on error */
+        txfree(rv);
+        rv = NULL;
     }
 
+    return rv;
+} /* end of function cx_log */
+
+
+
 void *
 cx_exp(void *data, short int type, int length, int *newlength, short int *newtype)
 {
@@ -614,20 +658,30 @@ cx_cosh(void *data, short int type, int length, int *newlength, short int *newty
     }
 }
 
-static double *
-d_tan(double *dd, int length)
+
+
+static double *d_tan(double *dd, int length)
 {
-    double *d;
-    int i;
+    int xrc = 0;
+    double *d = alloc_d(length);
 
-    d = alloc_d(length);
+    int i;
     for (i = 0; i < length; i++) {
-        rcheckn(tan(degtorad(dd[i])) != 0, "tan", d);
+        rcheck(tan(degtorad(dd[i])) != 0, "tan");
         d[i] = tan(degtorad(dd[i]));
     }
-    return d;
+
+EXITPOINT:
+    if (xrc != 0) { /* Free resources on error */
+        txfree(d);
+        d = (double *) NULL;
     }
 
+    return d;
+} /* end of function d_tan */
+
+
+
 static double *
 d_tanh(double *dd, int length)
 {
@@ -641,64 +695,93 @@ d_tanh(double *dd, int length)
     return d;
 }
 
-static ngcomplex_t *
-c_tan(ngcomplex_t *cc, int length)
+/* See https://proofwiki.org/wiki/Tangent_of_Complex_Number (formulation 4) among
+ * others for the tangent formula:
+
+ * sin z = sin(x + iy) = sin x cos(iy) + cos x sin(iy) = sin x cosh y + i cos x sinh y
+ * cos z = cos(x + iy) = cos x cos(iy) + sin x sin(iy) = cos x cosh y - i sin x sinh y
+ * tan z = ((sin x cosh y + i cos x sinh y) / (cos x cosh y - i sin x sinh y)) *
+            (cos x cosh y + isin x sinh y) / (cos x cosh y + i sin x sinh y)
+      = ...
+ *
+ *
+ * tan(a + bi) = (sin(2a) + i * sinh(2b)) / (cos(2a) + cosh(2b))
+ */
+static ngcomplex_t *c_tan(ngcomplex_t *cc, int length)
 {
-    ngcomplex_t *c;
-    int i;
+    ngcomplex_t * const c = alloc_c(length);
 
-    c = alloc_c(length);
+    int i;
     for (i = 0; i < length; i++) {
-        double u, v;
+        errno = 0;
+        ngcomplex_t *p_dst = c + i;
+        ngcomplex_t *p_src = cc + i;
+        const double a = p_src->cx_real;
+        const double b = p_src->cx_imag;
+        const double u = 2 * degtorad(a);
+        const double v = 2 * degtorad(b);
+        const double n_r = sin(u);
+        const double n_i = sinh(v);
+        const double d1 = cos(u);
+        const double d2 = cosh(v);
+        const double d = d1 + d2;
+        if (errno != 0 || d == 0.0) {
+            (void) fprintf(cp_err,
+                    "Invalid argument %lf + %lf i for compex tangent", a, b);
+            txfree(c);
+            return (ngcomplex_t *) NULL;
+        }
+        p_dst->cx_real = n_r / d;
+        p_dst->cx_imag = n_i / d;
+    } /* end of loop over elements in array */
+    return c;
+} /* end of function c_tan */
 
-        rcheckn(cos(degtorad(realpart(cc[i]))) *
-               cosh(degtorad(imagpart(cc[i]))), "tan", c);
-        rcheckn(sin(degtorad(realpart(cc[i]))) *
-               sinh(degtorad(imagpart(cc[i]))), "tan", c);
-        u = degtorad(realpart(cc[i]));
-        v = degtorad(imagpart(cc[i]));
-        /* The Lattice C compiler won't take multi-line macros, and
-         * CMS won't take >80 column lines....
-         */
-#define xx1 sin(u) * cosh(v)
-#define xx2 cos(u) * sinh(v)
-#define xx3 cos(u) * cosh(v)
-#define xx4 -sin(u) * sinh(v)
-        cdiv(xx1, xx2, xx3, xx4, realpart(c[i]), imagpart(c[i]));
+
+
+/* complex tanh function, uses tanh(z) = -i * tan(i * z).
+ * Unfortunately, it did so very inefficiently (a memory allocations per
+ * value calculated!) and leaked memory badly (all of the aforementioned
+ * allocations.) These issues were fixed 2020-03-04 */
+static ngcomplex_t *c_tanh(ngcomplex_t *cc, int length)
+{
+    ngcomplex_t * const tmp = alloc_c(length); /* i * z */
+
+    /* Build the i * z array to allow tan() to be called */
+    {
+        int i;
+        for (i = 0; i < length; ++i) {
+            ngcomplex_t *p_dst = tmp + i;
+            ngcomplex_t *p_src = cc + i;
+
+            /* multiply by i */
+            p_dst->cx_real = -p_src->cx_imag;
+            p_dst->cx_imag = p_src->cx_real;
         }
-    return c;
     }
 
-/* complex tanh function, uses tanh(z) = -i * tan (i * z) */
-static ngcomplex_t *
-c_tanh(ngcomplex_t *cc, int length)
+   /* Calculat tan(i * z), exiting on failure */
+    ngcomplex_t *const c = c_tan(tmp, length);
+    if (c == (ngcomplex_t *) NULL) {
+        txfree(tmp);
+        return (ngcomplex_t *) NULL;
+    }
+
+    /* Multiply by -i to find final result */
     {
-    ngcomplex_t *c, *s, *t;
         int i;
+        for (i = 0; i < length; ++i) {
+            ngcomplex_t *p_cur = c + i;
+            const double cx_real = p_cur->cx_real;
+            p_cur->cx_real = p_cur->cx_imag;
+            p_cur->cx_imag = -cx_real;
+        }
+    }
 
-    c = alloc_c(length);
-    s = alloc_c(1);
-    t = alloc_c(1);
-
-    for (i = 0; i < length; i++) {
-        /* multiply by i */
-        t[0].cx_real = -1. * imagpart(cc[i]);
-        t[0].cx_imag = realpart(cc[i]);
-        /* get complex tangent */
-        s = c_tan(t, 1);
-        /* if check in c_tan fails */
-        if (s == NULL) {
-            tfree(t);
-            return (NULL);
-        }
-        /* multiply by -i */
-        realpart(c[i]) = imagpart(s[0]);
-        imagpart(c[i]) = -1. * realpart(s[0]);
-    }
-    tfree(s);
-    tfree(t);
     return c;
-}
+} /* end of function c_tanh */
+
+
 
 void *
 cx_tan(void *data, short int type, int length, int *newlength, short int *newtype)
@@ -770,8 +853,8 @@ cx_sortorder(void *data, short int type, int length, int *newlength, short int *
     double *dd = (double *) data;
     int i;
 
-    amplitude_index_t *array_amplitudes;
-    array_amplitudes = (amplitude_index_t *) tmalloc(sizeof(amplitude_index_t) * (size_t) length);
+    amplitude_index_t * const array_amplitudes = (amplitude_index_t *)
+            tmalloc(sizeof(amplitude_index_t) * (size_t) length);
 
     *newlength = length;
     *newtype = VF_REAL;
@@ -788,7 +871,7 @@ cx_sortorder(void *data, short int type, int length, int *newlength, short int *
             d[i] = array_amplitudes[i].index;
     }
 
-    tfree(array_amplitudes);
+    txfree(array_amplitudes);
 
     /* Otherwise it is 0, but tmalloc zeros the stuff already. */
     return ((void *) d);

src/maths/cmaths/cmath2.c

@@ -361,11 +361,17 @@ cx_avg(void *data, short int type, int length, int *newlength, short int *newtyp
 
 /* Compute the mean of a vector. */
 
-void *
-cx_mean(void *data, short int type, int length, int *newlength, short int *newtype)
+void *cx_mean(void *data, short int type, int length,
+        int *newlength, short int *newtype)
 {
+    if (length == 0) {
+        (void) fprintf(cp_err, "mean calculation requires "
+                "at least one element.\n");
+        return NULL;
+    }
+
     *newlength = 1;
-    rcheck(length > 0, "mean");
+
     if (type == VF_REAL) {
         double *d;
         double *dd = (double *) data;
@@ -377,7 +383,8 @@ cx_mean(void *data, short int type, int length, int *newlength, short int *newty
             *d += dd[i];
         *d /= length;
         return ((void *) d);
-    } else {
+    }
+    else {
         ngcomplex_t *c;
         ngcomplex_t *cc = (ngcomplex_t *) data;
         int i;
@@ -397,43 +404,59 @@ cx_mean(void *data, short int type, int length, int *newlength, short int *newty
 
 /* Compute the standard deviation of all elements of a vector. */
 
-void *
-cx_stddev(void *data, short int type, int length, int *newlength, short int *newtype)
+void *cx_stddev(void *data, short int type, int length,
+        int *newlength, short int *newtype)
 {
+    if (length == 0) {
+        (void) fprintf(cp_err, "standard deviation calculation requires "
+                "at least one element.\n");
+        return NULL;
+    }
+
     *newlength = 1;
-    rcheck(length > 1, "stddev");
+
     if (type == VF_REAL) {
-        double *mean = (double *)cx_mean(data, type, length, newlength, newtype);
-        double *d, sum = 0.;
+        double *p_mean = (double *) cx_mean(data, type, length,
+                newlength, newtype);
+        double mean = *p_mean;
+        double *d, sum = 0.0;
         double *dd = (double *) data;
-        int i;
 
         d = alloc_d(1);
         *newtype = VF_REAL;
-        for (i = 0; i < length; i++)
-            sum += (dd[i] - *mean) * (dd[i] - *mean);
+        int i;
+        for (i = 0; i < length; i++) {
+            const double tmp = dd[i] - mean;
+            sum += tmp * tmp;
+        }
+
         *d = sqrt(sum / ((double) length - 1.0));
-        tfree(mean);
-        return ((void *)d);
+        txfree(p_mean);
+        return (void *) d;
     }
     else {
-        ngcomplex_t *cmean = (ngcomplex_t *)cx_mean(data, type, length, newlength, newtype);
-        double *d, sum = 0., a, b;
+        ngcomplex_t * const p_cmean = (ngcomplex_t *) cx_mean(data, type, length,
+                newlength, newtype);
+        const ngcomplex_t cmean = *p_cmean;
+        const double rmean = realpart(cmean);
+        const double imean = imagpart(cmean);
+        double *d, sum = 0.0;
         ngcomplex_t *cc = (ngcomplex_t *) data;
-        int i;
 
         d = alloc_d(1);
         *newtype = VF_REAL;
+        int i;
         for (i = 0; i < length; i++) {
-            a = realpart(cc[i]) - realpart(*cmean);
-            b = imagpart(cc[i]) - imagpart(*cmean);
+            const double a = realpart(cc[i]) - rmean;
+            const double b = imagpart(cc[i]) - imean;
             sum += a * a + b * b;
         }
         *d = sqrt(sum / ((double) length - 1.0));
-        tfree(cmean);
-        return ((void *)d);
-    }
+        txfree(p_cmean);
+        return (void *) d;
     }
+} /* end of function cx_stddev */
+
 
 
 void *
@@ -635,35 +658,41 @@ cx_times(void *data1, void *data2, short int datatype1, short int datatype2, int
 }
 
 
-void *
-cx_mod(void *data1, void *data2, short int datatype1, short int datatype2, int length)
+void *cx_mod(void *data1, void *data2, short int datatype1, short int datatype2,
+        int length)
 {
+    int xrc = 0;
+    void *rv;
     double *dd1 = (double *) data1;
     double *dd2 = (double *) data2;
-    double *d;
-    ngcomplex_t *cc1 = (ngcomplex_t *) data1;
-    ngcomplex_t *cc2 = (ngcomplex_t *) data2;
-    ngcomplex_t *c, c1, c2;
-    int i, r1, r2, i1, i2, r3, i3;
 
     if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
-        d = alloc_d(length);
+        double *d;
+        rv = d = alloc_d(length);
+
+        int i;
         for (i = 0; i < length; i++) {
-            r1 = (int)floor(fabs(dd1[i]));
+            const int r1 = (int) floor(fabs(dd1[i]));
             rcheck(r1 > 0, "mod");
-            r2 = (int)floor(fabs(dd2[i]));
+            const int r2 = (int)floor(fabs(dd2[i]));
             rcheck(r2 > 0, "mod");
-            r3 = r1 % r2;
+            const int r3 = r1 % r2;
             d[i] = (double) r3;
         }
-        return ((void *) d);
-    } else {
-        c = alloc_c(length);
+    }
+    else {
+        ngcomplex_t *c, c1, c2;
+        ngcomplex_t *cc1 = (ngcomplex_t *) data1;
+        ngcomplex_t *cc2 = (ngcomplex_t *) data2;
+        rv = c = alloc_c(length);
+
+        int i;
         for (i = 0; i < length; i++) {
             if (datatype1 == VF_REAL) {
                 realpart(c1) = dd1[i];
                 imagpart(c1) = 0.0;
-            } else {
+            }
+            else {
                 c1 = cc1[i];
             }
             if (datatype2 == VF_REAL) {
@@ -672,32 +701,45 @@ cx_mod(void *data1, void *data2, short int datatype1, short int datatype2, int l
             } else {
                 c2 = cc2[i];
             }
-            r1 = (int)floor(fabs(realpart(c1)));
+            const int r1 = (int) floor(fabs(realpart(c1)));
             rcheck(r1 > 0, "mod");
-            r2 = (int)floor(fabs(realpart(c2)));
+            const int r2 = (int) floor(fabs(realpart(c2)));
             rcheck(r2 > 0, "mod");
-            i1 = (int)floor(fabs(imagpart(c1)));
+            const int i1 = (int) floor(fabs(imagpart(c1)));
             rcheck(i1 > 0, "mod");
-            i2 = (int)floor(fabs(imagpart(c2)));
+            const int i2 = (int) floor(fabs(imagpart(c2)));
             rcheck(i2 > 0, "mod");
-            r3 = r1 % r2;
-            i3 = i1 % i2;
+            const int r3 = r1 % r2;
+            const int i3 = i1 % i2;
             realpart(c[i]) = (double) r3;
             imagpart(c[i]) = (double) i3;
         }
-        return ((void *) c);
     }
+
+EXITPOINT:
+    if (xrc != 0) { /* Free resources on error */
+        txfree(rv);
+        rv = NULL;
     }
 
+    return rv;
+} /* end of function cx_mod */
+
+
 
 /* Routoure JM : Compute the max of a vector. */
 
-void *
-cx_max(void *data, short int type, int length, int *newlength, short int *newtype)
+void *cx_max(void *data, short int type, int length,
+        int *newlength, short int *newtype)
 {
+    if (length == 0) {
+        (void) fprintf(cp_err, "maximum calculation requires "
+                "at least one element.\n");
+        return NULL;
+    }
+
     *newlength = 1;
-    /* test if length >0 et affiche un message d'erreur */
-    rcheck(length > 0, "mean");
+
     if (type == VF_REAL) {
         double largest = 0.0;
         double *d;
@@ -707,12 +749,16 @@ cx_max(void *data, short int type, int length, int *newlength, short int *newtyp
         d = alloc_d(1);
         *newtype = VF_REAL;
         largest = dd[0];
-        for (i = 1; i < length; i++)
-            if (largest < dd[i])
-                largest = dd[i];
+        for (i = 1; i < length; i++) {
+            const double tmp = dd[i];
+            if (largest < tmp) {
+                largest = tmp;
+            }
+        }
         *d = largest;
-        return ((void *) d);
-    } else {
+        return (void *) d;
+    }
+    else {
         double largest_real = 0.0;
         double largest_complex = 0.0;
         ngcomplex_t *c;
@@ -723,27 +769,37 @@ cx_max(void *data, short int type, int length, int *newlength, short int *newtyp
         *newtype = VF_COMPLEX;
         largest_real = realpart(*cc);
         largest_complex = imagpart(*cc);
-        for (i = 0; i < length; i++) {
-            if (largest_real < realpart(cc[i]))
-                largest_real = realpart(cc[i]);
-            if (largest_complex < imagpart(cc[i]))
-                largest_complex = imagpart(cc[i]);
+        for (i = 1; i < length; i++) {
+            const double tmpr = realpart(cc[i]);
+            if (largest_real < tmpr) {
+                largest_real = tmpr;
+            }
+            const double tmpi = imagpart(cc[i]);
+            if (largest_complex < tmpi) {
+                largest_complex = tmpi;
+            }
         }
         realpart(*c) = largest_real;
         imagpart(*c) = largest_complex;
-        return ((void *) c);
-    }
+        return (void *) c;
     }
+} /* end of function cx_max */
+
 
 
 /* Routoure JM : Compute the min of a vector. */
 
-void *
-cx_min(void *data, short int type, int length, int *newlength, short int *newtype)
+void *cx_min(void *data, short int type, int length,
+        int *newlength, short int *newtype)
 {
+    if (length == 0) {
+        (void) fprintf(cp_err, "minimum calculation requires "
+                "at least one element.\n");
+        return NULL;
+    }
+
     *newlength = 1;
-    /* test if length >0 et affiche un message d'erreur */
-    rcheck(length > 0, "mean");
+
     if (type == VF_REAL) {
         double smallest;
         double *d;
@@ -753,12 +809,16 @@ cx_min(void *data, short int type, int length, int *newlength, short int *newtyp
         d = alloc_d(1);
         *newtype = VF_REAL;
         smallest = dd[0];
-        for (i = 1; i < length; i++)
-            if (smallest > dd[i])
-                smallest = dd[i];
+        for (i = 1; i < length; i++) {
+            const double tmp = dd[i];
+            if (smallest > tmp) {
+                smallest = tmp;
+            }
+        }
         *d = smallest;
-        return ((void *) d);
-    } else {
+        return (void *) d;
+    }
+    else {
         double smallest_real;
         double smallest_complex;
         ngcomplex_t *c;
@@ -770,26 +830,35 @@ cx_min(void *data, short int type, int length, int *newlength, short int *newtyp
         smallest_real = realpart(*cc);
         smallest_complex = imagpart(*cc);
         for (i = 1; i < length; i++) {
-            if (smallest_real > realpart(cc[i]))
-                smallest_real = realpart(cc[i]);
-            if (smallest_complex > imagpart(cc[i]))
-                smallest_complex = imagpart(cc[i]);
+            const double tmpr = realpart(cc[i]);
+            if (smallest_real > tmpr) {
+                smallest_real = tmpr;
+            }
+            const double tmpi = imagpart(cc[i]);
+            if (smallest_complex > tmpi) {
+                smallest_complex = tmpi;
+            }
         }
         realpart(*c) = smallest_real;
         imagpart(*c) = smallest_complex;
-        return ((void *) c);
-    }
+        return (void *) c;
     }
+} /* end of function cx_min */
 
 
 /* Routoure JM : Compute the differential  of a vector. */
 
-void *
-cx_d(void *data, short int type, int length, int *newlength, short int *newtype)
+void *cx_d(void *data, short int type, int length,
+        int *newlength, short int *newtype)
 {
+    if (length == 0) {
+        (void) fprintf(cp_err, "differential calculation requires "
+                "at least one element.\n");
+        return NULL;
+    }
+
     *newlength = length;
-    /* test if length >0 et affiche un message d'erreur */
-    rcheck(length > 0, "deriv");
+
     if (type == VF_REAL) {
         double *d;
         double *dd = (double *) data;
@@ -799,12 +868,13 @@ cx_d(void *data, short int type, int length, int *newlength, short int *newtype)
         *newtype = VF_REAL;
         d[0] = dd[1] - dd[0];
         d[length-1] = dd[length-1] - dd[length-2];
-        for (i = 1; i < length - 1; i++)
+        for (i = 1; i < length - 1; i++) {
             d[i] = dd[i+1] - dd[i-1];
+        }
 
-        return ((void *) d);
-    } else {
-
+        return (void *) d;
+    }
+    else {
         ngcomplex_t *c;
         ngcomplex_t *cc = (ngcomplex_t *) data;
         int i;
@@ -820,11 +890,11 @@ cx_d(void *data, short int type, int length, int *newlength, short int *newtype)
         for (i = 1; i < length - 1; i++) {
             realpart(c[i]) = realpart(cc[i+1]) - realpart(cc[i-1]);
             imagpart(c[i]) = imagpart(cc[i+1]) - imagpart(cc[i-1]);
-
-        }
-        return ((void *) c);
         }
+
+        return (void *) c;
     }
+} /* end of function cx_d */
 
 
 void *

src/maths/cmaths/cmath3.c

@@ -24,51 +24,65 @@ Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
 
 
 static ngcomplex_t *cexp_sp3(ngcomplex_t *c); /* cexp exist's in some newer compiler */
-static ngcomplex_t *cln(ngcomplex_t *c);
-static ngcomplex_t *ctimes(ngcomplex_t *c1, ngcomplex_t *c2);
+static int cln(ngcomplex_t *c, ngcomplex_t *rv);
+static void ctimes(ngcomplex_t *c1, ngcomplex_t *c2, ngcomplex_t *rv);
 
-void *
-cx_divide(void *data1, void *data2, short int datatype1, short int datatype2, int length)
+void *cx_divide(void *data1, void *data2,
+        short int datatype1, short int datatype2, int length)
 {
+    int xrc = 0;
+    void *rv;
     double *dd1 = (double *) data1;
     double *dd2 = (double *) data2;
-    double *d;
     ngcomplex_t *cc1 = (ngcomplex_t *) data1;
     ngcomplex_t *cc2 = (ngcomplex_t *) data2;
     ngcomplex_t *c, c1, c2;
     int i;
 
     if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
-        d = alloc_d(length);
+        double *d;
+        rv = d = alloc_d(length);
         for (i = 0; i < length; i++) {
-            rcheckn(dd2[i] != 0, "divide", d);
+            rcheck(dd2[i] != 0, "divide");
             d[i] = dd1[i] / dd2[i];
         }
-        return ((void *) d);
-    } else {
-        c = alloc_c(length);
+    }
+    else {
+        rv = c = alloc_c(length);
         for (i = 0; i < length; i++) {
             if (datatype1 == VF_REAL) {
                 realpart(c1) = dd1[i];
                 imagpart(c1) = 0.0;
-            } else {
+            }
+            else {
                 c1 = cc1[i];
             }
             if (datatype2 == VF_REAL) {
                 realpart(c2) = dd2[i];
                 imagpart(c2) = 0.0;
-            } else {
+            }
+            else {
                 c2 = cc2[i];
             }
-        rcheckn((realpart(c2) != 0) || (imagpart(c2) != 0), "divide", c);
+        rcheck((realpart(c2) != 0) || (imagpart(c2) != 0), "divide");
 #define xx5 realpart(c1)
 #define xx6 imagpart(c1)
-cdiv(xx5, xx6, realpart(c2), imagpart(c2), realpart(c[i]), imagpart(c[i]));
+        cdiv(xx5, xx6, realpart(c2), imagpart(c2), realpart(c[i]),
+                imagpart(c[i]));
         }
-        return ((void *) c);
     }
+
+EXITPOINT:
+    if (xrc != 0) { /* Free resources on error */
+        tfree(rv);
+        rv = NULL;
     }
 
+    return rv;
+} /* end of function cx_divide */
+
+
+
 /* Should just use "j( )" */
 /* The comma operator. What this does (unless it is part of the argument
  * list of a user-defined function) is arg1 + j(arg2).
@@ -105,46 +119,61 @@ cx_comma(void *data1, void *data2, short int datatype1, short int datatype2, int
     return ((void *) c);
 }
 
-void *
-cx_power(void *data1, void *data2, short int datatype1, short int datatype2, int length)
+void *cx_power(void *data1, void *data2,
+        short int datatype1, short int datatype2, int length)
 {
+    int xrc = 0;
+    void *rv;
     double *dd1 = (double *) data1;
     double *dd2 = (double *) data2;
-    double *d;
-    ngcomplex_t *cc1 = (ngcomplex_t *) data1;
-    ngcomplex_t *cc2 = (ngcomplex_t *) data2;
-    ngcomplex_t *c, c1, c2, *t;
-    int i;
 
     if ((datatype1 == VF_REAL) && (datatype2 == VF_REAL)) {
-        d = alloc_d(length);
+        double *d;
+        rv = d = alloc_d(length);
+
+        int i;
         for (i = 0; i < length; i++) {
             rcheck((dd1[i] >= 0) || (floor(dd2[i]) == ceil(dd2[i])), "power");
             d[i] = pow(dd1[i], dd2[i]);
         }
-        return ((void *) d);
-    } else {
-        c = alloc_c(length);
+    }
+    else {
+        ngcomplex_t *cc1 = (ngcomplex_t *) data1;
+        ngcomplex_t *cc2 = (ngcomplex_t *) data2;
+        ngcomplex_t *c, c1, c2, *t;
+        rv = c = alloc_c(length);
+
+        int i;
         for (i = 0; i < length; i++) {
             if (datatype1 == VF_REAL) {
                 realpart(c1) = dd1[i];
                 imagpart(c1) = 0.0;
-            } else {
+            }
+            else {
                 c1 = cc1[i];
             }
             if (datatype2 == VF_REAL) {
                 realpart(c2) = dd2[i];
                 imagpart(c2) = 0.0;
-            } else {
+            }
+            else {
                 c2 = cc2[i];
             }
 
             if ((realpart(c1) == 0.0) && (imagpart(c1) == 0.0)) {
                 realpart(c[i]) = 0.0;
                 imagpart(c[i]) = 0.0;
-            } else { /* if ((imagpart(c1) != 0.0) && 
+            }
+            else { /* if ((imagpart(c1) != 0.0) && 
                         (imagpart(c2) != 0.0)) */
-                t = cexp_sp3(ctimes(&c2, cln(&c1)));
+                ngcomplex_t tmp, tmp2;
+                if (cln(&c1, &tmp) != 0) {
+                    (void) fprintf(cp_err, "power of 0 + i 0 not allowed.\n");
+                    xrc = -1;
+                    goto EXITPOINT;
+                }
+                ctimes(&c2, &tmp, &tmp2);
+                t = cexp_sp3(&tmp2);
                 c[i] = *t;
                 /*
                 } else {
@@ -154,10 +183,19 @@ cx_power(void *data1, void *data2, short int datatype1, short int datatype2, int
                 */
             }
         }
-        return ((void *) c);
     }
+
+EXITPOINT:
+    if (xrc != 0) { /* Free resources on error */
+        txfree(rv);
+        rv = NULL;
     }
 
+    return rv;
+} /* end of function cx_power */
+
+
+
 /* These are unnecessary... Only cx_power uses them... */
 
 static ngcomplex_t *
@@ -175,30 +213,35 @@ cexp_sp3(ngcomplex_t *c)
     return (&r);
 }
 
-static ngcomplex_t *
-cln(ngcomplex_t *c)
+static int cln(ngcomplex_t *c, ngcomplex_t *rv)
 {
-    static ngcomplex_t r;
+    double c_r = c->cx_real;
+    double c_i = c->cx_imag;
 
-    rcheck(cmag(*c) != 0, "ln");
-    realpart(r) = log(cmag(*c));
-    if (imagpart(*c) != 0.0)
-        imagpart(r) = atan2(imagpart(*c), realpart(*c));
-    else
-        imagpart(r) = 0.0;
-    return (&r);
+    if (c_r == 0 && c_i == 0) {
+        (void) fprintf(cp_err, "Complex log of 0 + i0 is undefined.\n");
+        return -1;
     }
 
-static ngcomplex_t *
-ctimes(ngcomplex_t *c1, ngcomplex_t *c2)
-{
-    static ngcomplex_t r;
+    rv->cx_real = log(cmag(*c));
+    if (c_i != 0.0) {
+        rv->cx_imag = atan2(c_i, c_r);
+    }
+    else {
+        rv->cx_imag = 0.0;
+    }
+    return 0;
+} /* end of functon cln */
+
 
-    realpart(r) = realpart(*c1) * realpart(*c2) - 
+
+static void ctimes(ngcomplex_t *c1, ngcomplex_t *c2, ngcomplex_t *rv)
+{
+    rv->cx_real = realpart(*c1) * realpart(*c2) -
                imagpart(*c1) * imagpart(*c2);
-    imagpart(r) = imagpart(*c1) * realpart(*c2) +
+    rv->cx_imag = imagpart(*c1) * realpart(*c2) +
                realpart(*c1) * imagpart(*c2);
-    return (&r);
+    return;
 }