ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/common/calexpr.c
Revision: 2.43
Committed: Fri Feb 23 03:47:57 2024 UTC (2 months, 3 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.42: +65 -37 lines
Log Message:
perf: Added array index optimization to calcomp routines

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: calexpr.c,v 2.42 2022/04/08 23:32:25 greg Exp $";
3 #endif
4 /*
5 * Compute data values using expression parser
6 *
7 * 7/1/85 Greg Ward
8 *
9 * 11/11/85 Made channel input conditional with (INCHAN) compiles.
10 *
11 * 4/2/86 Added conditional compiles for function definitions (FUNCTION).
12 *
13 * 1/29/87 Made variables conditional (VARIABLE)
14 *
15 * 5/19/88 Added constant subexpression elimination (RCONST)
16 *
17 * 2/19/03 Eliminated conditional compiles in favor of esupport extern.
18 */
19
20 #include "copyright.h"
21
22 #include <ctype.h>
23 #include <errno.h>
24 #include <math.h>
25 #include <stdlib.h>
26
27 #include "rtmisc.h"
28 #include "rtio.h"
29 #include "rterror.h"
30 #include "calcomp.h"
31
32 #define MAXLINE 256 /* maximum line length */
33
34 #define newnode() (EPNODE *)ecalloc(1, sizeof(EPNODE))
35
36 #define isdecimal(c) (isdigit(c) | ((c) == '.'))
37
38 static double euminus(EPNODE *), eargument(EPNODE *), enumber(EPNODE *);
39 static double echannel(EPNODE *);
40 static double eadd(EPNODE *), esubtr(EPNODE *),
41 emult(EPNODE *), edivi(EPNODE *),
42 epow(EPNODE *);
43 static double ebotch(EPNODE *);
44
45 unsigned int esupport = /* what to support */
46 E_VARIABLE | E_FUNCTION ;
47
48 int eofc = 0; /* optional end-of-file character */
49 int nextc; /* lookahead character */
50
51 double (*eoper[])(EPNODE *) = { /* expression operations */
52 ebotch,
53 evariable,
54 enumber,
55 euminus,
56 echannel,
57 efunc,
58 eargument,
59 ebotch,
60 ebotch,
61 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
62 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
63 emult,
64 eadd,
65 0,
66 esubtr,
67 0,
68 edivi,
69 0,0,0,0,0,0,0,0,0,0,
70 ebotch,
71 0,0,
72 ebotch,
73 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
74 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
75 epow,
76 };
77
78 static FILE *infp; /* input file pointer */
79 static char *linbuf; /* line buffer */
80 static char *infile; /* input file name */
81 static int lineno; /* input line number */
82 static int linepos; /* position in buffer */
83
84
85 EPNODE *
86 eparse( /* parse an expression string */
87 char *expr
88 )
89 {
90 EPNODE *ep;
91
92 initstr(expr, NULL, 0);
93 curfunc = NULL;
94 ep = getE1();
95 if (nextc != EOF)
96 syntax("unexpected character");
97 return(ep);
98 }
99
100
101 double
102 eval( /* evaluate an expression string */
103 char *expr
104 )
105 {
106 int prev_support = esupport;
107 EPNODE *ep;
108 double rval;
109
110 esupport &= ~E_RCONST; /* don't bother reducing constant expr */
111 ep = eparse(expr);
112 esupport = prev_support; /* as you were */
113 rval = evalue(ep);
114 epfree(ep,1);
115 return(rval);
116 }
117
118
119 int
120 epcmp( /* compare two expressions for equivalence */
121 EPNODE *ep1,
122 EPNODE *ep2
123 )
124 {
125 double d;
126
127 if (ep1->type != ep2->type)
128 return(1);
129
130 switch (ep1->type) {
131
132 case VAR:
133 return(ep1->v.ln != ep2->v.ln);
134
135 case NUM:
136 if (ep2->v.num == 0)
137 return(ep1->v.num != 0);
138 d = ep1->v.num / ep2->v.num;
139 return((d > 1.000000000001) | (d < 0.999999999999));
140
141 case CHAN:
142 case ARG:
143 return(ep1->v.chan != ep2->v.chan);
144
145 case '=':
146 case ':':
147 return(epcmp(ep1->v.kid->sibling, ep2->v.kid->sibling));
148
149 case CLKT:
150 case SYM: /* should never get this one */
151 return(0);
152
153 default:
154 ep1 = ep1->v.kid;
155 ep2 = ep2->v.kid;
156 while (ep1 != NULL) {
157 if (ep2 == NULL)
158 return(1);
159 if (epcmp(ep1, ep2))
160 return(1);
161 ep1 = ep1->sibling;
162 ep2 = ep2->sibling;
163 }
164 return(ep2 != NULL);
165 }
166 }
167
168
169 void
170 epfree( /* free a parse tree */
171 EPNODE *epar,
172 int frep
173 )
174 {
175 EPNODE *ep;
176
177 switch (epar->type) {
178
179 case VAR:
180 varfree(epar->v.ln);
181 break;
182
183 case SYM:
184 freestr(epar->v.name);
185 break;
186
187 case NUM:
188 case CHAN:
189 case ARG:
190 case CLKT:
191 break;
192
193 default:
194 if (epar->nkids < 0) {
195 ep = epar->v.kid - epar->nkids;
196 while (ep > epar->v.kid)
197 epfree(--ep, 0);
198 efree(ep); /* free array space */
199 } else
200 while ((ep = epar->v.kid) != NULL) {
201 epar->v.kid = ep->sibling;
202 epfree(ep, 1);
203 }
204 break;
205
206 }
207 if (frep)
208 efree(epar);
209 }
210
211
212 void
213 epoptimize( /* realloc lists as arrays if > length 3 */
214 EPNODE *epar
215 )
216 {
217 EPNODE *ep;
218
219 if (epar->nkids > 3) { /* do this node if > 3 kids */
220 int n = 1;
221 epar->v.kid = (EPNODE *)erealloc(epar->v.kid,
222 sizeof(EPNODE)*epar->nkids);
223 while (n < epar->nkids) {
224 ep = epar->v.kid[n-1].sibling;
225 epar->v.kid[n] = *ep;
226 efree(ep); /* not epfree()! */
227 epar->v.kid[n-1].sibling = epar->v.kid + n;
228 n++;
229 }
230 epar->nkids = -epar->nkids;
231 }
232 if (epar->nkids) /* do children if any */
233 for (ep = epar->v.kid; ep != NULL; ep = ep->sibling)
234 epoptimize(ep);
235 }
236
237 /* the following used to be a switch */
238 static double
239 eargument(
240 EPNODE *ep
241 )
242 {
243 return(argument(ep->v.chan));
244 }
245
246 static double
247 enumber(
248 EPNODE *ep
249 )
250 {
251 return(ep->v.num);
252 }
253
254 static double
255 euminus(
256 EPNODE *ep
257 )
258 {
259 EPNODE *ep1 = ep->v.kid;
260
261 return(-evalue(ep1));
262 }
263
264 static double
265 echannel(
266 EPNODE *ep
267 )
268 {
269 return(chanvalue(ep->v.chan));
270 }
271
272 static double
273 eadd(
274 EPNODE *ep
275 )
276 {
277 EPNODE *ep1 = ep->v.kid;
278
279 return(evalue(ep1) + evalue(ep1->sibling));
280 }
281
282 static double
283 esubtr(
284 EPNODE *ep
285 )
286 {
287 EPNODE *ep1 = ep->v.kid;
288
289 return(evalue(ep1) - evalue(ep1->sibling));
290 }
291
292 static double
293 emult(
294 EPNODE *ep
295 )
296 {
297 EPNODE *ep1 = ep->v.kid;
298
299 return(evalue(ep1) * evalue(ep1->sibling));
300 }
301
302 static double
303 edivi(
304 EPNODE *ep
305 )
306 {
307 EPNODE *ep1 = ep->v.kid;
308 double d;
309
310 d = evalue(ep1->sibling);
311 if (d == 0.0) {
312 wputs("Division by zero\n");
313 errno = ERANGE;
314 return(0.0);
315 }
316 return(evalue(ep1) / d);
317 }
318
319 static double
320 epow(
321 EPNODE *ep
322 )
323 {
324 EPNODE *ep1 = ep->v.kid;
325 double d;
326 int lasterrno;
327
328 lasterrno = errno;
329 errno = 0;
330 d = pow(evalue(ep1), evalue(ep1->sibling));
331 #ifdef isnan
332 if (errno == 0) {
333 if (isnan(d))
334 errno = EDOM;
335 else if (isinf(d))
336 errno = ERANGE;
337 }
338 #endif
339 if ((errno == EDOM) | (errno == ERANGE)) {
340 wputs("Illegal power\n");
341 return(0.0);
342 }
343 errno = lasterrno;
344 return(d);
345 }
346
347 static double
348 ebotch(
349 EPNODE *ep
350 )
351 {
352 eputs("Bad expression!\n");
353 quit(1);
354 return 0.0; /* pro forma return */
355 }
356
357
358 EPNODE *
359 ekid( /* return pointer to a node's nth kid */
360 EPNODE *ep,
361 int n
362 )
363 {
364 if (ep->nkids < 0) { /* allocated array? */
365 if (n >= -ep->nkids)
366 return(NULL);
367 return(ep->v.kid + n);
368 }
369 ep = ep->v.kid; /* else get from list */
370 while (n-- > 0)
371 if ((ep = ep->sibling) == NULL)
372 break;
373 return(ep);
374 }
375
376
377 void
378 initfile( /* prepare input file */
379 FILE *fp,
380 char *fn,
381 int ln
382 )
383 {
384 static char inpbuf[MAXLINE];
385
386 infp = fp;
387 linbuf = inpbuf;
388 infile = fn;
389 lineno = ln;
390 linepos = 0;
391 inpbuf[0] = '\0';
392 scan();
393 }
394
395
396 void
397 initstr( /* prepare input string */
398 char *s,
399 char *fn,
400 int ln
401 )
402 {
403 infp = NULL;
404 infile = fn;
405 lineno = ln;
406 linbuf = s;
407 linepos = 0;
408 scan();
409 }
410
411
412 void
413 getscanpos( /* return current scan position */
414 char **fnp,
415 int *lnp,
416 char **spp,
417 FILE **fpp
418 )
419 {
420 if (fnp != NULL) *fnp = infile;
421 if (lnp != NULL) *lnp = lineno;
422 if (spp != NULL) *spp = linbuf+linepos;
423 if (fpp != NULL) *fpp = infp;
424 }
425
426
427 int
428 scan(void) /* scan next character, return literal next */
429 {
430 int lnext = 0;
431
432 do {
433 if (linbuf[linepos] == '\0')
434 if (infp == NULL || fgets(linbuf, MAXLINE, infp) == NULL)
435 nextc = EOF;
436 else {
437 nextc = linbuf[0];
438 lineno++;
439 linepos = 1;
440 }
441 else
442 nextc = linbuf[linepos++];
443 if (!lnext)
444 lnext = nextc;
445 if (nextc == eofc) {
446 nextc = EOF;
447 break;
448 }
449 if (nextc == '{') {
450 scan();
451 while (nextc != '}')
452 if (nextc == EOF)
453 syntax("'}' expected");
454 else
455 scan();
456 scan();
457 }
458 } while (isspace(nextc));
459 return(lnext);
460 }
461
462
463 char *
464 long2ascii( /* convert long to ascii */
465 long l
466 )
467 {
468 static char buf[16];
469 char *cp;
470 int neg = 0;
471
472 if (l == 0)
473 return("0");
474 if (l < 0) {
475 l = -l;
476 neg++;
477 }
478 cp = buf + sizeof(buf);
479 *--cp = '\0';
480 while (l) {
481 *--cp = l % 10 + '0';
482 l /= 10;
483 }
484 if (neg)
485 *--cp = '-';
486 return(cp);
487 }
488
489
490 void
491 syntax( /* report syntax error and quit */
492 char *err
493 )
494 {
495 int i;
496
497 if ((infile != NULL) | (lineno != 0)) {
498 if (infile != NULL) eputs(infile);
499 if (lineno != 0) {
500 eputs(infile != NULL ? ", line " : "line ");
501 eputs(long2ascii((long)lineno));
502 }
503 eputs(":\n");
504 }
505 eputs(linbuf);
506 if (linbuf[strlen(linbuf)-1] != '\n')
507 eputs("\n");
508 for (i = 0; i < linepos-1; i++)
509 eputs(linbuf[i] == '\t' ? "\t" : " ");
510 eputs("^ ");
511 eputs(err);
512 eputs("\n");
513 quit(1);
514 }
515
516
517 void
518 addekid( /* add a child to ep */
519 EPNODE *ep,
520 EPNODE *ek
521 )
522 {
523 if (ep->nkids < 0) {
524 eputs("Cannot add child after optimization\n");
525 quit(1);
526 }
527 ep->nkids++;
528 if (ep->v.kid == NULL)
529 ep->v.kid = ek;
530 else {
531 for (ep = ep->v.kid; ep->sibling != NULL; ep = ep->sibling)
532 ;
533 ep->sibling = ek;
534 }
535 ek->sibling = NULL;
536 }
537
538
539 char *
540 getname(void) /* scan an identifier */
541 {
542 static char str[RMAXWORD+1];
543 int i, lnext;
544
545 lnext = nextc;
546 for (i = 0; i < RMAXWORD && isid(lnext); i++, lnext = scan())
547 str[i] = lnext;
548 str[i] = '\0';
549 while (isid(lnext)) /* skip rest of name */
550 lnext = scan();
551
552 return(str);
553 }
554
555
556 int
557 getinum(void) /* scan a positive integer */
558 {
559 int n, lnext;
560
561 n = 0;
562 lnext = nextc;
563 while (isdigit(lnext)) {
564 n = n * 10 + lnext - '0';
565 lnext = scan();
566 }
567 return(n);
568 }
569
570
571 double
572 getnum(void) /* scan a positive float */
573 {
574 int i, lnext;
575 char str[RMAXWORD+1];
576
577 i = 0;
578 lnext = nextc;
579 while (isdigit(lnext) && i < RMAXWORD) {
580 str[i++] = lnext;
581 lnext = scan();
582 }
583 if ((lnext == '.') & (i < RMAXWORD)) {
584 str[i++] = lnext;
585 lnext = scan();
586 if (i == 1 && !isdigit(lnext))
587 syntax("badly formed number");
588 while (isdigit(lnext) && i < RMAXWORD) {
589 str[i++] = lnext;
590 lnext = scan();
591 }
592 }
593 if ((lnext == 'e') | (lnext == 'E') && i < RMAXWORD) {
594 str[i++] = lnext;
595 lnext = scan();
596 if ((lnext == '-') | (lnext == '+') && i < RMAXWORD) {
597 str[i++] = lnext;
598 lnext = scan();
599 }
600 if (!isdigit(lnext))
601 syntax("missing exponent");
602 while (isdigit(lnext) && i < RMAXWORD) {
603 str[i++] = lnext;
604 lnext = scan();
605 }
606 }
607 str[i] = '\0';
608
609 return(atof(str));
610 }
611
612
613 EPNODE *
614 getE1(void) /* E1 -> E1 ADDOP E2 */
615 /* E2 */
616 {
617 EPNODE *ep1, *ep2;
618
619 ep1 = getE2();
620 while ((nextc == '+') | (nextc == '-')) {
621 ep2 = newnode();
622 ep2->type = nextc;
623 scan();
624 addekid(ep2, ep1);
625 addekid(ep2, getE2());
626 if (esupport&E_RCONST &&
627 (ep1->type == NUM) & (ep1->sibling->type == NUM))
628 ep2 = rconst(ep2);
629 ep1 = ep2;
630 }
631 return(ep1);
632 }
633
634
635 EPNODE *
636 getE2(void) /* E2 -> E2 MULOP E3 */
637 /* E3 */
638 {
639 EPNODE *ep1, *ep2;
640
641 ep1 = getE3();
642 while ((nextc == '*') | (nextc == '/')) {
643 ep2 = newnode();
644 ep2->type = nextc;
645 scan();
646 addekid(ep2, ep1);
647 addekid(ep2, getE3());
648 if (esupport&E_RCONST) {
649 EPNODE *ep3 = ep1->sibling;
650 if ((ep1->type == NUM) & (ep3->type == NUM)) {
651 ep2 = rconst(ep2);
652 } else if (ep3->type == NUM) {
653 if (ep2->type == '/') {
654 if (ep3->v.num == 0)
655 syntax("divide by zero constant");
656 ep2->type = '*'; /* for speed */
657 ep3->v.num = 1./ep3->v.num;
658 } else if (ep3->v.num == 0) {
659 ep1->sibling = NULL; /* (E2 * 0) */
660 epfree(ep2,1);
661 ep2 = ep3;
662 }
663 } else if (ep1->type == NUM && ep1->v.num == 0) {
664 epfree(ep3,1); /* (0 * E3) or (0 / E3) */
665 ep1->sibling = NULL;
666 efree(ep2);
667 ep2 = ep1;
668 }
669 }
670 ep1 = ep2;
671 }
672 return(ep1);
673 }
674
675
676 EPNODE *
677 getE3(void) /* E3 -> E4 ^ E3 */
678 /* E4 */
679 {
680 EPNODE *ep1, *ep2;
681
682 ep1 = getE4();
683 if (nextc != '^')
684 return(ep1);
685 ep2 = newnode();
686 ep2->type = nextc;
687 scan();
688 addekid(ep2, ep1);
689 addekid(ep2, getE3());
690 if (esupport&E_RCONST) {
691 EPNODE *ep3 = ep1->sibling;
692 if ((ep1->type == NUM) & (ep3->type == NUM)) {
693 ep2 = rconst(ep2);
694 } else if (ep1->type == NUM && ep1->v.num == 0) {
695 epfree(ep3,1); /* (0 ^ E3) */
696 ep1->sibling = NULL;
697 efree(ep2);
698 ep2 = ep1;
699 } else if ((ep3->type == NUM && ep3->v.num == 0) |
700 (ep1->type == NUM && ep1->v.num == 1)) {
701 epfree(ep2,1); /* (E4 ^ 0) or (1 ^ E3) */
702 ep2 = newnode();
703 ep2->type = NUM;
704 ep2->v.num = 1;
705 } else if (ep3->type == NUM && ep3->v.num == 1) {
706 efree(ep3); /* (E4 ^ 1) */
707 ep1->sibling = NULL;
708 efree(ep2);
709 ep2 = ep1;
710 }
711 }
712 return(ep2);
713 }
714
715
716 EPNODE *
717 getE4(void) /* E4 -> ADDOP E5 */
718 /* E5 */
719 {
720 EPNODE *ep1, *ep2;
721
722 if (nextc == '-') {
723 scan();
724 ep2 = getE5();
725 if (ep2->type == NUM) {
726 ep2->v.num = -ep2->v.num;
727 return(ep2);
728 }
729 if (ep2->type == UMINUS) { /* don't generate -(-E5) */
730 ep1 = ep2->v.kid;
731 efree(ep2);
732 return(ep1);
733 }
734 ep1 = newnode();
735 ep1->type = UMINUS;
736 addekid(ep1, ep2);
737 return(ep1);
738 }
739 if (nextc == '+')
740 scan();
741 return(getE5());
742 }
743
744
745 EPNODE *
746 getE5(void) /* E5 -> (E1) */
747 /* VAR */
748 /* NUM */
749 /* $N */
750 /* FUNC(E1,..) */
751 /* ARG */
752 {
753 int i;
754 char *nam;
755 EPNODE *ep1, *ep2;
756
757 if (nextc == '(') {
758 scan();
759 ep1 = getE1();
760 if (nextc != ')')
761 syntax("')' expected");
762 scan();
763 return(ep1);
764 }
765
766 if (esupport&E_INCHAN && nextc == '$') {
767 scan();
768 ep1 = newnode();
769 ep1->type = CHAN;
770 ep1->v.chan = getinum();
771 return(ep1);
772 }
773
774 if (esupport&(E_VARIABLE|E_FUNCTION) &&
775 (isalpha(nextc) | (nextc == CNTXMARK))) {
776 nam = getname();
777 ep1 = NULL;
778 if ((esupport&(E_VARIABLE|E_FUNCTION)) == (E_VARIABLE|E_FUNCTION)
779 && curfunc != NULL)
780 for (i = 1, ep2 = curfunc->v.kid->sibling;
781 ep2 != NULL; i++, ep2 = ep2->sibling)
782 if (!strcmp(ep2->v.name, nam)) {
783 ep1 = newnode();
784 ep1->type = ARG;
785 ep1->v.chan = i;
786 break;
787 }
788 if (ep1 == NULL) {
789 ep1 = newnode();
790 ep1->type = VAR;
791 ep1->v.ln = varinsert(nam);
792 }
793 if (esupport&E_FUNCTION && nextc == '(') {
794 ep2 = newnode();
795 ep2->type = FUNC;
796 addekid(ep2, ep1);
797 ep1 = ep2;
798 do {
799 scan();
800 addekid(ep1, getE1());
801 } while (nextc == ',');
802 if (nextc != ')')
803 syntax("')' expected");
804 scan();
805 } else if (!(esupport&E_VARIABLE))
806 syntax("'(' expected");
807 if (esupport&E_RCONST && isconstvar(ep1))
808 ep1 = rconst(ep1);
809 return(ep1);
810 }
811
812 if (isdecimal(nextc)) {
813 ep1 = newnode();
814 ep1->type = NUM;
815 ep1->v.num = getnum();
816 return(ep1);
817 }
818 syntax("unexpected character");
819 return NULL; /* pro forma return */
820 }
821
822
823 EPNODE *
824 rconst( /* reduce a constant expression */
825 EPNODE *epar
826 )
827 {
828 EPNODE *ep;
829
830 ep = newnode();
831 ep->type = NUM;
832 errno = 0;
833 ep->v.num = evalue(epar);
834 if ((errno == EDOM) | (errno == ERANGE))
835 syntax("bad constant expression");
836 epfree(epar,1);
837
838 return(ep);
839 }
840
841
842 int
843 isconstvar( /* is ep linked to a constant expression? */
844 EPNODE *ep
845 )
846 {
847 EPNODE *ep1;
848
849 if (esupport&E_FUNCTION && ep->type == FUNC) {
850 if (!isconstfun(ep->v.kid))
851 return(0);
852 for (ep1 = ep->v.kid->sibling; ep1 != NULL; ep1 = ep1->sibling)
853 if (ep1->type != NUM && !isconstfun(ep1))
854 return(0);
855 return(1);
856 }
857 if (ep->type != VAR)
858 return(0);
859 ep1 = ep->v.ln->def;
860 if (ep1 == NULL || ep1->type != ':')
861 return(0);
862 if (esupport&E_FUNCTION && ep1->v.kid->type != SYM)
863 return(0);
864 return(1);
865 }
866
867
868 int
869 isconstfun( /* is ep linked to a constant function? */
870 EPNODE *ep
871 )
872 {
873 EPNODE *dp;
874 LIBR *lp;
875
876 if (ep->type != VAR)
877 return(0);
878 if ((dp = ep->v.ln->def) != NULL) {
879 if (dp->v.kid->type == FUNC)
880 return(dp->type == ':');
881 else
882 return(0); /* don't identify masked library functions */
883 }
884 if ((lp = ep->v.ln->lib) != NULL)
885 return(lp->atyp == ':');
886 return(0);
887 }