--- ray/src/common/calexpr.c	2008/05/01 16:42:05	2.32
+++ ray/src/common/calexpr.c	2024/02/25 04:11:10	2.47
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: calexpr.c,v 2.32 2008/05/01 16:42:05 greg Exp $";
+static const char	RCSid[] = "$Id: calexpr.c,v 2.47 2024/02/25 04:11:10 greg Exp $";
 #endif
 /*
  *  Compute data values using expression parser
@@ -19,8 +19,6 @@ static const char	RCSid[] = "$Id: calexpr.c,v 2.32 200
 
 #include "copyright.h"
 
-#include  <stdio.h>
-#include  <string.h>
 #include  <ctype.h>
 #include  <errno.h>
 #include  <math.h>
@@ -35,8 +33,10 @@ static const char	RCSid[] = "$Id: calexpr.c,v 2.32 200
 
 #define	 newnode()	(EPNODE *)ecalloc(1, sizeof(EPNODE))
 
-#define	 isdecimal(c)	(isdigit(c) || (c) == '.')
+#define	 isdecimal(c)	(isdigit(c) | ((c) == '.'))
 
+#define  envalue(ep)	((ep)->type==NUM ? (ep)->v.num : evalue(ep))
+
 static double  euminus(EPNODE *), eargument(EPNODE *), enumber(EPNODE *);
 static double  echannel(EPNODE *);
 static double  eadd(EPNODE *), esubtr(EPNODE *),
@@ -105,20 +105,23 @@ eval(			/* evaluate an expression string */
     char  *expr
 )
 {
-    register EPNODE  *ep;
+    int  prev_support = esupport;
+    EPNODE  *ep;
     double  rval;
 
+    esupport &= ~E_RCONST;	/* don't bother reducing constant expr */
     ep = eparse(expr);
+    esupport = prev_support;	/* as you were */
     rval = evalue(ep);
-    epfree(ep);
+    epfree(ep,1);
     return(rval);
 }
 
 
 int
 epcmp(			/* compare two expressions for equivalence */
-    register EPNODE  *ep1,
-    register EPNODE  *ep2
+    EPNODE  *ep1,
+    EPNODE  *ep2
 )
 {
 	double  d;
@@ -167,10 +170,11 @@ epcmp(			/* compare two expressions for equivalence */
 
 void
 epfree(			/* free a parse tree */
-    register EPNODE	 *epar
+    EPNODE	 *epar,
+    int		frep
 )
 {
-    register EPNODE  *ep;
+    EPNODE	*ep;
 
     switch (epar->type) {
 
@@ -189,17 +193,78 @@ epfree(			/* free a parse tree */
 	    break;
 
 	default:
-	    while ((ep = epar->v.kid) != NULL) {
-		epar->v.kid = ep->sibling;
-		epfree(ep);
-	    }
+	    if (epar->nkids < 0) {
+	    	ep = epar->v.kid - epar->nkids;
+	    	while (ep > epar->v.kid)
+	    		epfree(--ep, 0);
+		efree(ep);	/* free array space */
+	    } else
+	    	while ((ep = epar->v.kid) != NULL) {
+		    epar->v.kid = ep->sibling;
+		    epfree(ep, 1);
+		}
 	    break;
 
     }
+    if (frep)
+    	efree(epar);
+}
 
-    efree((char *)epar);
+
+static void
+epflatten(			/* flatten hierarchies for '+', '*' */
+	EPNODE *epar
+)
+{
+    EPNODE	*ep;
+
+    if (epar->nkids < 0) {
+    	eputs("Cannot flatten EPNODE array\n");
+    	quit(1);
+    }
+    for (ep = epar->v.kid; ep != NULL; ep = ep->sibling)
+    	while (ep->type == epar->type) {
+	    EPNODE	*ep1 = ep->v.kid;
+	    while (ep1->sibling != NULL)
+	    	ep1 = ep1->sibling;
+	    ep1->sibling = ep->sibling;
+	    epar->nkids += nekids(ep) - 1;
+	    ep1 = ep->v.kid;
+	    *ep = *ep1;
+	    efree(ep1);		/* not epfree()! */
+	}
 }
 
+
+void
+epoptimize(			/* flatten operations, lists -> arrays */
+	EPNODE	*epar
+)
+{
+    EPNODE	*ep;
+
+    if ((epar->type == '+') | (epar->type == '*'))
+    	epflatten(epar);	/* flatten associative operations */
+
+    if (epar->nkids)		/* do children if any */
+    	for (ep = epar->v.kid; ep != NULL; ep = ep->sibling)
+	    epoptimize(ep);
+
+    if (epar->nkids > 4) {	/* make list into array if > 4 kids */
+        int	n = 1;
+    	epar->v.kid = (EPNODE *)erealloc(epar->v.kid,
+    					sizeof(EPNODE)*epar->nkids);
+    	while (n < epar->nkids) {
+	    ep = epar->v.kid[n-1].sibling;
+	    epar->v.kid[n] = *ep;
+	    efree(ep);		/* not epfree()! */
+	    epar->v.kid[n-1].sibling = epar->v.kid + n;
+	    n++;
+	}
+	epar->nkids = -n;
+    }
+}
+
 				/* the following used to be a switch */
 static double
 eargument(
@@ -222,7 +287,7 @@ euminus(
     EPNODE	*ep
 )
 {
-    register EPNODE  *ep1 = ep->v.kid;
+    EPNODE  *ep1 = ep->v.kid;
 
     return(-evalue(ep1));
 }
@@ -240,9 +305,14 @@ eadd(
     EPNODE	*ep
 )
 {
-    register EPNODE  *ep1 = ep->v.kid;
+    double  sum = 0;
+    EPNODE  *ep1 = ep->v.kid;
 
-    return(evalue(ep1) + evalue(ep1->sibling));
+    do
+    	sum += envalue(ep1);
+    while ((ep1 = ep1->sibling) != NULL);
+
+    return(sum);
 }
 
 static double
@@ -250,9 +320,10 @@ esubtr(
     EPNODE	*ep
 )
 {
-    register EPNODE  *ep1 = ep->v.kid;
+    EPNODE  *ep1 = ep->v.kid;
+    EPNODE  *ep2 = ep1->sibling;
 
-    return(evalue(ep1) - evalue(ep1->sibling));
+    return(envalue(ep1) - envalue(ep2));
 }
 
 static double
@@ -260,9 +331,14 @@ emult(
     EPNODE	*ep
 )
 {
-    register EPNODE  *ep1 = ep->v.kid;
+    double  prod = 1;
+    EPNODE  *ep1 = ep->v.kid;
 
-    return(evalue(ep1) * evalue(ep1->sibling));
+    do
+    	prod *= envalue(ep1);
+    while ((ep1 = ep1->sibling) != NULL);
+
+    return(prod);
 }
 
 static double
@@ -270,16 +346,17 @@ edivi(
     EPNODE	*ep
 )
 {
-    register EPNODE  *ep1 = ep->v.kid;
+    EPNODE  *ep1 = ep->v.kid;
+    EPNODE  *ep2 = ep1->sibling;
     double  d;
 
-    d = evalue(ep1->sibling);
+    d = evalue(ep2);
     if (d == 0.0) {
 	wputs("Division by zero\n");
 	errno = ERANGE;
 	return(0.0);
     }
-    return(evalue(ep1) / d);
+    return(envalue(ep1) / d);
 }
 
 static double
@@ -287,7 +364,7 @@ epow(
     EPNODE	*ep
 )
 {
-    register EPNODE  *ep1 = ep->v.kid;
+    EPNODE  *ep1 = ep->v.kid;
     double  d;
     int	 lasterrno;
 
@@ -295,13 +372,14 @@ epow(
     errno = 0;
     d = pow(evalue(ep1), evalue(ep1->sibling));
 #ifdef  isnan
-    if (errno == 0)
+    if (errno == 0) {
 	if (isnan(d))
 	    errno = EDOM;
 	else if (isinf(d))
 	    errno = ERANGE;
+    }
 #endif
-    if (errno == EDOM || errno == ERANGE) {
+    if ((errno == EDOM) | (errno == ERANGE)) {
 	wputs("Illegal power\n");
 	return(0.0);
     }
@@ -322,33 +400,23 @@ ebotch(
 
 EPNODE *
 ekid(			/* return pointer to a node's nth kid */
-    register EPNODE	 *ep,
-    register int  n
+    EPNODE	 *ep,
+    int  n
 )
 {
-
-    for (ep = ep->v.kid; ep != NULL; ep = ep->sibling)
-	if (--n < 0)
-	    break;
-
+    if (ep->nkids < 0) {	/* allocated array? */
+    	if (n >= -ep->nkids)
+	    return(NULL);
+    	return(ep->v.kid + n);
+    }
+    ep = ep->v.kid;		/* else get from list */
+    while (n-- > 0)
+    	if ((ep = ep->sibling) == NULL)
+		break;
     return(ep);
 }
 
 
-int
-nekids(			/* return # of kids for node ep */
-    register EPNODE	 *ep
-)
-{
-    register int  n = 0;
-
-    for (ep = ep->v.kid; ep != NULL; ep = ep->sibling)
-	n++;
-
-    return(n);
-}
-
-
 void
 initfile(		/* prepare input file */
     FILE  *fp,
@@ -402,7 +470,7 @@ getscanpos(	/* return current scan position */
 int
 scan(void)		/* scan next character, return literal next */
 {
-    register int  lnext = 0;
+    int  lnext = 0;
 
     do {
 	if (linbuf[linepos] == '\0')
@@ -441,7 +509,7 @@ long2ascii(			      /* convert long to ascii */
 )
 {
     static char	 buf[16];
-    register char  *cp;
+    char  *cp;
     int	 neg = 0;
 
     if (l == 0)
@@ -467,9 +535,9 @@ syntax(			/* report syntax error and quit */
     char  *err
 )
 {
-    register int  i;
+    int  i;
 
-    if (infile != NULL || lineno != 0) {
+    if ((infile != NULL) | (lineno != 0)) {
 	if (infile != NULL) eputs(infile);
 	if (lineno != 0) {
 	    eputs(infile != NULL ? ", line " : "line ");
@@ -491,18 +559,23 @@ syntax(			/* report syntax error and quit */
 
 void
 addekid(			/* add a child to ep */
-    register EPNODE	 *ep,
-    EPNODE	*ekid
+    EPNODE	 *ep,
+    EPNODE	*ek
 )
 {
+    if (ep->nkids < 0) {
+    	eputs("Cannot add kid to EPNODE array\n");
+    	quit(1);
+    }
+    ep->nkids++;
     if (ep->v.kid == NULL)
-	ep->v.kid = ekid;
+	ep->v.kid = ek;
     else {
 	for (ep = ep->v.kid; ep->sibling != NULL; ep = ep->sibling)
 	    ;
-	ep->sibling = ekid;
+	ep->sibling = ek;
     }
-    ekid->sibling = NULL;
+    ek->sibling = NULL;		/* shouldn't be necessary */
 }
 
 
@@ -510,7 +583,7 @@ char *
 getname(void)			/* scan an identifier */
 {
     static char	 str[RMAXWORD+1];
-    register int  i, lnext;
+    int  i, lnext;
 
     lnext = nextc;
     for (i = 0; i < RMAXWORD && isid(lnext); i++, lnext = scan())
@@ -526,7 +599,7 @@ getname(void)			/* scan an identifier */
 int
 getinum(void)			/* scan a positive integer */
 {
-    register int  n, lnext;
+    int  n, lnext;
 
     n = 0;
     lnext = nextc;
@@ -541,7 +614,7 @@ getinum(void)			/* scan a positive integer */
 double
 getnum(void)			/* scan a positive float */
 {
-    register int  i, lnext;
+    int  i, lnext;
     char  str[RMAXWORD+1];
 
     i = 0;
@@ -550,7 +623,7 @@ getnum(void)			/* scan a positive float */
 	str[i++] = lnext;
 	lnext = scan();
     }
-    if (lnext == '.' && i < RMAXWORD) {
+    if ((lnext == '.') & (i < RMAXWORD)) {
 	str[i++] = lnext;
 	lnext = scan();
 	if (i == 1 && !isdigit(lnext))
@@ -581,20 +654,20 @@ getnum(void)			/* scan a positive float */
 
 
 EPNODE *
-getE1(void)				/* E1 -> E1 ADDOP E2 */
+getE1(void)			/* E1 -> E1 ADDOP E2 */
 				/*	 E2 */
 {
-    register EPNODE  *ep1, *ep2;
+    EPNODE  *ep1, *ep2;
 
     ep1 = getE2();
-    while (nextc == '+' || nextc == '-') {
+    while ((nextc == '+') | (nextc == '-')) {
 	ep2 = newnode();
 	ep2->type = nextc;
 	scan();
 	addekid(ep2, ep1);
 	addekid(ep2, getE2());
 	if (esupport&E_RCONST &&
-			ep1->type == NUM && ep1->sibling->type == NUM)
+			(ep1->type == NUM) & (ep1->sibling->type == NUM))
 		ep2 = rconst(ep2);
 	ep1 = ep2;
     }
@@ -603,21 +676,40 @@ getE1(void)				/* E1 -> E1 ADDOP E2 */
 
 
 EPNODE *
-getE2(void)				/* E2 -> E2 MULOP E3 */
+getE2(void)			/* E2 -> E2 MULOP E3 */
 				/*	 E3 */
 {
-    register EPNODE  *ep1, *ep2;
+    EPNODE  *ep1, *ep2;
 
     ep1 = getE3();
-    while (nextc == '*' || nextc == '/') {
+    while ((nextc == '*') | (nextc == '/')) {
 	ep2 = newnode();
 	ep2->type = nextc;
 	scan();
 	addekid(ep2, ep1);
 	addekid(ep2, getE3());
-	if (esupport&E_RCONST &&
-			ep1->type == NUM && ep1->sibling->type == NUM)
-		ep2 = rconst(ep2);
+	if (esupport&E_RCONST) {
+		EPNODE	*ep3 = ep1->sibling;
+		if ((ep1->type == NUM) & (ep3->type == NUM)) {
+			ep2 = rconst(ep2);
+		} else if (ep3->type == NUM) {
+			if (ep2->type == '/') {
+				if (ep3->v.num == 0)
+					syntax("divide by zero constant");
+				ep2->type = '*';	/* for speed */
+				ep3->v.num = 1./ep3->v.num;
+			} else if (ep3->v.num == 0) {
+				ep1->sibling = NULL;	/* (E2 * 0) */
+				epfree(ep2,1);
+				ep2 = ep3;
+			}
+		} else if (ep1->type == NUM && ep1->v.num == 0) {
+			epfree(ep3,1);		/* (0 * E3) or (0 / E3) */
+			ep1->sibling = NULL;
+			efree(ep2);
+			ep2 = ep1;
+		}
+	}
 	ep1 = ep2;
     }
     return(ep1);
@@ -625,32 +717,50 @@ getE2(void)				/* E2 -> E2 MULOP E3 */
 
 
 EPNODE *
-getE3(void)				/* E3 -> E4 ^ E3 */
+getE3(void)			/* E3 -> E4 ^ E3 */
 				/*	 E4 */
 {
-    register EPNODE  *ep1, *ep2;
+	EPNODE  *ep1, *ep2;
 
-    ep1 = getE4();
-    if (nextc == '^') {
+	ep1 = getE4();
+	if (nextc != '^')
+		return(ep1);
 	ep2 = newnode();
 	ep2->type = nextc;
 	scan();
 	addekid(ep2, ep1);
 	addekid(ep2, getE3());
-	if (esupport&E_RCONST &&
-			ep1->type == NUM && ep1->sibling->type == NUM)
-		ep2 = rconst(ep2);
+	if (esupport&E_RCONST) {
+		EPNODE	*ep3 = ep1->sibling;
+		if ((ep1->type == NUM) & (ep3->type == NUM)) {
+			ep2 = rconst(ep2);
+		} else if (ep1->type == NUM && ep1->v.num == 0) {
+			epfree(ep3,1);		/* (0 ^ E3) */
+			ep1->sibling = NULL;
+			efree(ep2);
+			ep2 = ep1;
+		} else if ((ep3->type == NUM && ep3->v.num == 0) |
+				(ep1->type == NUM && ep1->v.num == 1)) {
+			epfree(ep2,1);		/* (E4 ^ 0) or (1 ^ E3) */
+			ep2 = newnode();
+			ep2->type = NUM;
+			ep2->v.num = 1;
+		} else if (ep3->type == NUM && ep3->v.num == 1) {
+			efree(ep3);	/* (E4 ^ 1) */
+			ep1->sibling = NULL;
+			efree(ep2);
+			ep2 = ep1;
+		}
+	}
 	return(ep2);
-    }
-    return(ep1);
 }
 
 
 EPNODE *
-getE4(void)				/* E4 -> ADDOP E5 */
+getE4(void)			/* E4 -> ADDOP E5 */
 				/*	 E5 */
 {
-    register EPNODE  *ep1, *ep2;
+    EPNODE  *ep1, *ep2;
 
     if (nextc == '-') {
 	scan();
@@ -660,8 +770,9 @@ getE4(void)				/* E4 -> ADDOP E5 */
 		return(ep2);
 	}
 	if (ep2->type == UMINUS) {	/* don't generate -(-E5) */
-	    efree((char *)ep2);
-	    return(ep2->v.kid);
+	    ep1 = ep2->v.kid;
+	    efree(ep2);
+	    return(ep1);
 	}
 	ep1 = newnode();
 	ep1->type = UMINUS;
@@ -684,7 +795,7 @@ getE5(void)			/* E5 -> (E1) */
 {
 	int	 i;
 	char  *nam;
-	register EPNODE  *ep1, *ep2;
+	EPNODE  *ep1, *ep2;
 
 	if (nextc == '(') {
 		scan();
@@ -694,7 +805,6 @@ getE5(void)			/* E5 -> (E1) */
 		scan();
 		return(ep1);
 	}
-
 	if (esupport&E_INCHAN && nextc == '$') {
 		scan();
 		ep1 = newnode();
@@ -702,9 +812,8 @@ getE5(void)			/* E5 -> (E1) */
 		ep1->v.chan = getinum();
 		return(ep1);
 	}
-
 	if (esupport&(E_VARIABLE|E_FUNCTION) &&
-			(isalpha(nextc) || nextc == CNTXMARK)) {
+			(isalpha(nextc) | (nextc == CNTXMARK))) {
 		nam = getname();
 		ep1 = NULL;
 		if ((esupport&(E_VARIABLE|E_FUNCTION)) == (E_VARIABLE|E_FUNCTION)
@@ -740,7 +849,6 @@ getE5(void)			/* E5 -> (E1) */
 			ep1 = rconst(ep1);
 		return(ep1);
 	}
-
 	if (isdecimal(nextc)) {
 		ep1 = newnode();
 		ep1->type = NUM;
@@ -754,18 +862,18 @@ getE5(void)			/* E5 -> (E1) */
 
 EPNODE *
 rconst(			/* reduce a constant expression */
-    register EPNODE	 *epar
+    EPNODE	 *epar
 )
 {
-    register EPNODE  *ep;
+    EPNODE  *ep;
 
     ep = newnode();
     ep->type = NUM;
     errno = 0;
     ep->v.num = evalue(epar);
-    if (errno == EDOM || errno == ERANGE)
+    if ((errno == EDOM) | (errno == ERANGE))
 	syntax("bad constant expression");
-    epfree(epar);
+    epfree(epar,1);
  
     return(ep);
 }
@@ -773,10 +881,10 @@ rconst(			/* reduce a constant expression */
 
 int
 isconstvar(			/* is ep linked to a constant expression? */
-    register EPNODE	 *ep
+    EPNODE	 *ep
 )
 {
-    register EPNODE  *ep1;
+    EPNODE  *ep1;
 
     if (esupport&E_FUNCTION && ep->type == FUNC) {
 	if (!isconstfun(ep->v.kid))
@@ -799,11 +907,11 @@ isconstvar(			/* is ep linked to a constant expression
 
 int
 isconstfun(			/* is ep linked to a constant function? */
-    register EPNODE	 *ep
+    EPNODE	 *ep
 )
 {
-    register EPNODE  *dp;
-    register LIBR  *lp;
+    EPNODE  *dp;
+    LIBR  *lp;
 
     if (ep->type != VAR)
 	return(0);