Added undocumented LN function.

diff --git a/cupl.doc b/cupl.doc
index b2700aa067d819f1f5790a9c86e40a08437e7fc6..46a3c3b0911f24c2de53235265b5071a2a26849b 100644 (file)
--- a/cupl.doc
+++ b/cupl.doc
@@ -8,6 +8,12 @@ references (by section-page number in the original manual) to notes which
 follow the appendix transcription.  These notes are also excerpts from the
 manual.
 
+There is one correction in the text.  The original manual listed both LN
+and LOG as built-in functions and wrote "LOG(a)" is "natural log of a".  We
+believe this is incorrect; we have changed the "LOG(a)" in the original to
+LN(a) and inserted a new "LOG(a)" entry.  This implementation's LOG function
+is, accordingly, log_10() and not log_e().
+
 There are a few typographical changes to fit it into the ASCII character set.
 The differences:
 
@@ -160,7 +166,9 @@ ELEMENTS OF THE LANGUAGE
                                                    expression
           FLOOR(a)      greatest integer not     numerical valued
                           exceeding a              expression
-          LOG(a)        natural og of a          numerical valued
+          LN(a)         natural log of a         numerical valued
+                                                   expression
+          LOG(a)        log to base 10 of a      numerical valued
                                                    expression
           SQRT(a)       positive square root     numerical valued
                           of a                     expression

diff --git a/cupl.h b/cupl.h
index 98d570f9eda8da9498db1c4931786f2f89279870..5b063628f49aa4659689f07f39c31d2ffa476023 100644 (file)
--- a/cupl.h
+++ b/cupl.h
@@ -94,6 +94,11 @@ extern lvar *idlist;
 
 #define for_symbols(s)    for (s = idlist; s; s = s->next)
 
+/* subscripting operations */
+#define SUB(v, i, j)   (v.elements + i * v.width + j)
+#define SUBI(v, n)     (n / v.width)
+#define SUBJ(v, n)     (n % v.width)
+
 #define NOMEM  "out of memory\n"
 
 /* miscellaneous */
@@ -111,6 +116,7 @@ extern void warn(char *msg, ...);
 
 extern void make_scalar(value *v, scalar n);
 extern value copy_value(value);
+extern value allocate_value(int rank, int i, int j);
 extern void deallocate_value(value *);
 
 void cupl_reset_write();
@@ -130,13 +136,21 @@ extern value cupl_atan(value);
 extern value cupl_cos(value);
 extern value cupl_exp(value);
 extern value cupl_floor(value);
+extern value cupl_ln(value);
 extern value cupl_log(value);
 extern value cupl_sqrt(value);
 extern value cupl_max(value, value);
 extern value cupl_min(value, value);
 extern value cupl_rand(value);
 
-extern value cupl_sgm(value right);
+extern value cupl_det(value);
+extern value cupl_dot(value, value);
+extern value cupl_inv(value);
+extern value cupl_posmax(value);
+extern value cupl_posmin(value);
+extern value cupl_sgm(value);
+extern value cupl_trc(value);
+extern value cupl_trn(value);
 
 extern bool cupl_eq(value, value);
 extern bool cupl_lt(value, value);

diff --git a/execute.c b/execute.c
index f0a846c54cb5e983d6b02670d7025aa4e172ff65..d338512befdb7e7a6d87712467983ac3694f8801 100644 (file)
--- a/execute.c
+++ b/execute.c
@@ -234,6 +234,13 @@ static value cupl_eval(node *tree)
        RETURN_WRAP(tree, tree->car, tree->cdr, result)
        return(result);
 
+    case LN:
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_ln(rightside);
+       deallocate_value(&rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
     case SQRT:
        rightside = EVAL_WRAP(cupl_eval(tree->cdr));
        result = cupl_sqrt(rightside);
@@ -268,6 +275,41 @@ static value cupl_eval(node *tree)
         * Matrix functions
         */
 
+    case DET:
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_det(rightside);
+       deallocate_value(&rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
+    case DOT:
+       leftside = EVAL_WRAP(cupl_eval(tree->car));
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_dot(leftside, rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
+    case INV:
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_inv(rightside);
+       deallocate_value(&rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
+    case POSMAX:
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_posmax(rightside);
+       deallocate_value(&rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
+    case POSMIN:
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_posmin(rightside);
+       deallocate_value(&rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
     case SGM:
        rightside = EVAL_WRAP(cupl_eval(tree->cdr));
        result = cupl_sgm(rightside);
@@ -275,6 +317,20 @@ static value cupl_eval(node *tree)
        RETURN_WRAP(tree, tree->car, tree->cdr, result)
        return(result);
 
+    case TRC:
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_trc(rightside);
+       deallocate_value(&rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
+    case TRN:
+       rightside = EVAL_WRAP(cupl_eval(tree->cdr));
+       result = cupl_trn(rightside);
+       deallocate_value(&rightside);
+       RETURN_WRAP(tree, tree->car, tree->cdr, result)
+       return(result);
+
        /*
         * Relations.
         *
@@ -345,6 +401,14 @@ static value cupl_eval(node *tree)
        RETURN_WRAP(tree, tree->car, tree->cdr, result)
        return(result);
 
+       /*
+        * Allocation.
+        */
+
+    case ALLOCATE:
+       /* FIXME: implement ALLOCATE */
+       die("ALLOCATE is not implemented\n");
+
        /*
         * Control structures.
         */
@@ -359,6 +423,26 @@ static value cupl_eval(node *tree)
        result.rank = FAIL;
        break;
 
+    case FOR:
+       /* FIXME: implement FOR */
+       die("FOR is not implemented\n");
+
+    case TIMES:
+       /* FIXME: implement TIMES */
+       die("TIMES is not implemented\n");
+
+    case BLOCK:
+       /* FIXME: implement BLOCK */
+       die("BLOCK skip is not implemented\n");
+
+    case GO:
+       /* FIXME: implement GO TO */
+       die("GO TO is not implemented\n");
+
+    case IF:
+       /* FIXME: implement IF */
+       die("GO TO is not implemented\n");
+
     case STOP:
        longjmp(jmpbuf, 1);
        /* no fall through */

diff --git a/monitor.c b/monitor.c
index 879606bc7e043115f4b93d5923529eb756dac3ce..8a81e5517741282676158dad2ba51570518f1ac2 100644 (file)
--- a/monitor.c
+++ b/monitor.c
@@ -9,6 +9,7 @@ SYNOPSIS
 
     void make_scalar(value *v, scalar i)
     void copy_value(value v);
+    value allocate_value(int rank, int i, int j)
     void deallocate_value(value *v)
 
     void cupl_reset_write()
@@ -27,13 +28,27 @@ SYNOPSIS
     value cupl_cos(value)
     value cupl_exp(value)
     value cupl_floor(value)
+    value cupl_ln(value)
     value cupl_log(value)
     value cupl_sqrt(value)
     value cupl_max(value, value)
     value cupl_min(value, value)
     value cupl_rand(value)
 
+    bool cupl_eq(value, value)
+    bool cupl_lt(value, value)
+    bool cupl_gt(value, value)
+    bool cupl_le(value, value)
+    bool cupl_ge(value, value)
+
+    value cupl_det(value);
+    value cupl_dot(value, vlue);
+    value cupl_inv(value);
+    value cupl_posmax(value);
+    value cupl_posmin(value);
     value cupl_sgm(value right)
+    value cupl_trc(value right)
+    value cupl_trn(value right)
 
 DESCRIPTION 
    Runtime support.  This is segregated from the execute() code in case anyone
@@ -107,12 +122,25 @@ value copy_value(value v)
     return(newvalue);
 }
 
+value allocate_value(int rank, int i, int j)
+/* allocate a value of given shape */
+{
+    value v;
+
+    v.rank = rank;
+    v.width = j; v.depth = i;
+    v.elements = (scalar *)calloc(sizeof(scalar), i * j);
+
+    return(v);
+}
+
 void deallocate_value(value *v)
 /* destroy a value copy, only if its reference count is 1 */
 {
     (void) free(v->elements);
 }
 
+
 /****************************************************************************
  *
  * Output support
@@ -144,6 +172,7 @@ static void needspace(int w)
 void cupl_scalar_write(char *name, scalar quant)
 /* write a numeric or skip a field in CUPL style */
 {
+    /* FIXME: cupl_scalar_write screws up line wraps */
     if (name)
     {
        needspace(2 * fieldwidth);
@@ -179,7 +208,7 @@ value cupl_add(value left, value right)
 /* add two CUPL values */
 {
     if (!CONGRUENT(left, right))
-       die("addition failed, operands of different sizes\n");
+       die("addition failed, operands of different sizes or ranks\n");
     else
     {
        value   result;
@@ -196,7 +225,7 @@ value cupl_subtract(value left, value right)
 /* subtract two CUPL values */
 {
     if (!CONGRUENT(left, right))
-       die("subtract failed, operands of different sizes\n");
+       die("subtract failed, operands of different sizes or ranks\n");
     else
     {
        value   result;
@@ -212,6 +241,7 @@ value cupl_subtract(value left, value right)
 value cupl_multiply(value left, value right)
 /* multiply two CUPL values */
 {
+    /* FIXME: implement cupl_multiply for non-scalar operands */
     if (left.rank == 0 && right.rank == 0)
     {
        value   result;
@@ -227,6 +257,7 @@ value cupl_multiply(value left, value right)
 value cupl_divide(value left, value right)
 /* divide two CUPL values */
 {
+    /* FIXME: implement cupl_divide for non-scalar operands */
     if (left.rank == 0 && right.rank == 0)
     {
        value   result;
@@ -343,7 +374,7 @@ value cupl_floor(value right)
     }
 }
 
-value cupl_log(value right)
+value cupl_ln(value right)
 /* apply ln function */
 {
     value      result;
@@ -358,6 +389,21 @@ value cupl_log(value right)
     }
 }
 
+value cupl_log(value right)
+/* apply log10 function */
+{
+    value      result;
+
+    if (right.rank != 0)
+       die("LOG is only defined for scalar arguments\n");
+    else
+    {
+       make_scalar(&result, 0);
+       result.elements[0] = log10(right.elements[0]);
+       return(result);
+    }
+}
+
 value cupl_sqrt(value right)
 /* apply square root */
 {
@@ -383,30 +429,28 @@ value cupl_max(value left, value right)
 /* apply max function */
 {
     value      result;
+    int        n;
 
-    if (right.rank != 0 || left.rank != 0)
-       die("MAX is not yet supported for vectors and matrices\n");
-    else
-    {
-       make_scalar(&result, 0);
-       result.elements[0] = max(right.elements[0], left.elements[0]);
-       return(result);
-    }
+    make_scalar(&result, right.elements[0]);
+    for (n = 0; n < right.width * right.depth; n++)
+       if (result.elements[0] < right.elements[n])
+           result.elements[0] = right.elements[n];
+
+    return(result);
 }
 
 value cupl_min(value left, value right)
 /* apply min functionnot yet supported for vectors and matrices */
 {
     value      result;
+    int        n;
 
-    if (right.rank != 0 || left.rank != 0)
-       die("MIN is \n");
-    else
-    {
-       make_scalar(&result, 0);
-       result.elements[0] = min(right.elements[0], left.elements[0]);
-       return(result);
-    }
+    make_scalar(&result, right.elements[0]);
+    for (n = 0; n < right.width * right.depth; n++)
+       if (result.elements[0] > right.elements[n])
+           result.elements[0] = right.elements[n];
+
+    return(result);
 }
 
 value cupl_rand(value right)
@@ -425,13 +469,200 @@ value cupl_rand(value right)
     }
 }
 
+/****************************************************************************
+ *
+ * Relations
+ *
+ ****************************************************************************/
+
+static bool normalize(scalar *ep1, scalar *ep2)
+/* normalize a pair of scalars for conversion purposes */
+{
+    int        mult;
+
+    if (*ep1 == 0 && *ep2 == 0)
+       return;
+
+    /*
+     * Note: this is not precisely the roundoff method used by original
+     * CUPL, though it should (famous last words) be equivalent.
+     */
+    if (fabs(*ep1) > fabs(*ep2))
+       mult = pow(10, log10(*ep2) + 14);
+    else
+       mult = pow(10, log10(*ep1) + 14);
+
+    *ep1 *= mult;
+    *ep2 *= mult;
+
+    *ep1 = floor(*ep1);
+    *ep2 = floor(*ep2);
+}
+
+bool cupl_eq(value v1, value v2)
+/* test any two CUPL values for pairwise equality */
+{
+    if (!CONGRUENT(v1, v2))
+       die("comparison failed, operands of different sizes or ranks\n");
+    else
+    {
+       int     n;
+
+       for (n = 0; n < v2.width * v2.depth; n++)
+       {
+           scalar e1 = v1.elements[n], e2 = v1.elements[n];
+
+           normalize(&e1, &e2);
+
+           if (e1 != e2)
+               return(FALSE);
+       }
+
+       return(TRUE);
+    }
+}
+
+bool cupl_le(value v1, value v2)
+/* are all elements of v1 less than or equal to their correspondents in v2? */
+{
+    if (!CONGRUENT(v1, v2))
+       die("LE failed, operands of different sizes or ranks\n");
+    else
+    {
+       bool equal = TRUE;
+       int     n;
+
+
+       for (n = 0; n < v2.width * v2.depth; n++)
+       {
+           scalar e1 = v1.elements[n], e2 = v1.elements[n];
+
+           normalize(&e1, &e2);
+
+           if (e1 > e2)
+               return(FALSE);
+       }
+
+       return(TRUE);
+    }
+}
+
+bool cupl_ge(value v1, value v2)
+/* are all els of v1 greater than or equal to their correspondents in v2? */
+{
+    if (!CONGRUENT(v1, v2))
+       die("GE failed, operands of different sizes or ranks\n");
+    else
+    {
+       int     n;
+
+       for (n = 0; n < v2.width * v2.depth; n++)
+       {
+           scalar e1 = v1.elements[n], e2 = v1.elements[n];
+
+           normalize(&e1, &e2);
+
+           if (e1 < e2)
+               return(FALSE);
+       }
+
+       return(TRUE);
+    }
+}
+
+bool cupl_lt(value v1, value v2)
+/* strange CUPL definition #1, see appendix A */
+{
+    return(cupl_le(v1, v2) && !cupl_eq(v1, v2));
+}
+
+bool cupl_gt(value v1, value v2)
+/* strange CUPL definition #2, see appendix A */
+{
+    return(cupl_ge(v1, v2) && !cupl_eq(v1, v2));
+}
+
 /****************************************************************************
  *
  * Matrix functions
  *
  ****************************************************************************/
 
-/* FIXME: most not yet implemented */
+value cupl_det(value right)
+{
+    /* FIXME: implement DET */
+    die("the determinant function is not yet implemented");
+}
+
+value cupl_dot(value left, value right)
+/* compute inner or dot product of two vectors */
+{
+    if (!CONGRUENT(left, right) || right.rank != 1)
+       die("DOT failed, operands of different sizes or ranks\n");
+    else
+    {
+       value   result;
+       int     n;
+
+       make_scalar(&result, 0);
+       for (n = 0; n < left.width * left.depth; n++)
+           result.elements[0] += left.elements[n] * right.elements[n];
+       return(result);
+    }
+}
+
+value cupl_inv(value right)
+/* compute the inverse of a matrix */
+{
+    /* FIXME: implement INV */
+    die("the inverse function is not yet implemented\n");
+}
+
+value cupl_posmax(value right)
+/* row position of maximum element of argument */
+{
+    if (right.rank == 0)
+       die("POSMAX argument is a scalar");
+    else
+    {
+       value   result;
+       int     n;
+       scalar  maxel = right.elements[0];;
+
+       make_scalar(&result, 0);
+       for (n = 0; n < right.width * right.depth; n++)
+           if (maxel < right.elements[n])
+           {
+               result.elements[0] = SUBI(right, n);
+               maxel = right.elements[0];
+           }
+
+       return(result);
+    }
+}
+
+value cupl_posmin(value right)
+/* row position of minimum element of argument */
+{
+    if (right.rank == 0)
+       die("POSMIN argument is a scalar");
+    else
+    {
+       value   result;
+       int     n;
+       scalar  minel = right.elements[0];;
+
+       make_scalar(&result, 0);
+       for (n = 0; n < right.width * right.depth; n++)
+           if (minel > right.elements[n])
+           {
+               result.elements[0] = SUBI(right, n);
+               minel = right.elements[n];
+           }
+
+       return(result);
+    }
+}
 
 value cupl_sgm(value right)
 /* compute the sum of the elements of a matrix */
@@ -445,6 +676,34 @@ value cupl_sgm(value right)
     return(result);
 }
 
-/* monitor.c ends here */
+value cupl_trc(value right)
+/* compute sum of elements on main diagonals */
+{
+    if (right.rank != 2 || right.width != right.depth)
+       die("TRACE failed, operands is not a square matrix\n");
+    else
+    {
+       value   result;
+       int     n;
+
+       make_scalar(&result, 0);
+       for (n = 0; n < right.width; n++)
+           result.elements[0] += SUB(right, n, n)[0];
+       return(result);
+    }
+}
+
+value cupl_trn(value right)
+/* compute the transpose of a matrix */
+{
+    value result = allocate_value(right.rank, right.depth, right.width);
+    int        i, j;
 
+    for (i = 0; i < right.depth; i++)
+       for (j = 0; j < right.width; j++)
+           SUB(result, j, i)[0] = SUB(right, i, j)[0];
 
+    return(result);
+}
+
+/* monitor.c ends here */