1 #ifndef lint
2 static char *sccsid ="exprn.c (CWI) 1.1 85/03/01";
3 #endif
4 #include "ideal.h"
5 #include "y.tab.h"
6
nonlinerr(funcname)7 void nonlinerr (funcname)
8 char *funcname;
9 {
10 fprintf (stderr, "ideal: %s() of unknown\n >>>Returning 1.0\n", funcname);
11 }
12
13 static DEPPTR depvarlist = NULL;
14 static boolean incon_warn = TRUE;
15 static boolean nl_warn;
16 boolean nl_fail;
17 static EQNPTR nl_eqns = NULL;
18
expreval(exprn,givennoad)19 INTLPTR expreval (exprn, givennoad)
20 EXPR exprn;
21 NOADPTR givennoad;
22 {
23 /* This routine returns an INTLPTR whose operator
24 is ';'--a promoted commanode containing the
25 dependency list representing the real part in
26 its left field, the imag part in its right */
27 register INTLPTR intl;
28 register EXTLPTR extl;
29 if (!exprn)
30 return (commagen (0.0, 0.0));
31 if (((EXTLPTR)exprn)->leaf) {
32 extl = (EXTLPTR) exprn;
33 dprintf "At a leaf of kind %d\n", extl->kind);
34 switch (extl->kind) {
35 case PATH:
36 return (pathfind (extl->info.path, givennoad));
37 break;
38 case CONST:
39 return (commagen (extl->info.const, 0.0));
40 break;
41 }
42 }
43 intl = (INTLPTR) exprn;
44 if (intl->oper == NAME) {
45 dprintf "Looking for a function named %s\n", idprint ((int) intl->left));
46 } else {
47 dprintf "At an internal node with operator %c\n", intl->oper);
48 }
49 switch (intl->oper) {
50 INTLPTR lefttemp, righttemp, temp, temp2;
51 DEPPTR drek, drek2;
52 float repart, impart, modulus;
53 case NAME:
54 if (((int) intl->left) == lookup ("re")) {
55 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
56 depfree ((DEPPTR) temp->right);
57 temp->right = (EXPR) depgen ((VARPTR) NULL, 0.0);
58 return (temp);
59 } else if (((int) intl->left) == lookup ("im")) {
60 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
61 depfree ((DEPPTR) temp->left);
62 temp->left = temp->right;
63 temp->right = (EXPR) depgen ((VARPTR) NULL, 0.0);
64 return (temp);
65 } else if (((int) intl->left) == lookup ("conj")) {
66 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
67 temp2 = intlgen (
68 ';',
69 (EXPR) depadd (
70 (DEPPTR) NULL, 0.0,
71 (DEPPTR) temp->left, 1.0
72 ),
73 (EXPR) depadd (
74 (DEPPTR) NULL, 0.0,
75 (DEPPTR) temp->right, -1.0
76 )
77 );
78 intlfree (temp);
79 return (temp2);
80 } else if (((int) intl->left) == lookup ("abs")) {
81 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
82 if (!known(temp)) {
83 if (nl_warn)
84 nonlinerr ("abs");
85 nl_fail = !nl_warn;
86 intlfree (temp);
87 return (commagen (1.0, 0.0));
88 } else {
89 repart = Re(temp);
90 impart = Im(temp);
91 intlfree (temp);
92 return (commagen (hypot (repart, impart), 0.0));
93 }
94 } else if (((int) intl->left) == lookup ("cis")) {
95 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
96 if (!known(temp)) {
97 if (nl_warn)
98 nonlinerr ("cis");
99 nl_fail = !nl_warn;
100 intlfree (temp);
101 return (commagen (1.0, 0.0));
102 } else {
103 repart = Re(temp);
104 impart = Im(temp);
105 if (!radflag) {
106 dtor(repart);
107 dtor(impart);
108 }
109 intlfree (temp);
110 if (impart > EPSILON)
111 fprintf (stderr, "ideal: cis of complex value\n >>>Ignoring imaginary part\n");
112 return (commagen (cos (repart), sin (repart)));
113 }
114 } else if (((int) intl->left) == lookup ("int")) {
115 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
116 if (!known(temp)) {
117 if (nl_warn)
118 nonlinerr ("int");
119 nl_fail = !nl_warn;
120 intlfree (temp);
121 return (commagen (1.0,0.0));
122 } else {
123 double intpart;
124 repart = Re(temp);
125 impart = Im(temp);
126 intlfree (temp);
127 if (impart > EPSILON)
128 fprintf (stderr, "ideal: int of complex value\n >>>Ignoring imaginary part\n");
129 modf (repart, &intpart);
130 return (commagen ((float) intpart, 0.0));
131 }
132 } else if (((int) intl->left) == lookup ("atan2")
133 || ((int) intl->left) == lookup ("angle")) {
134 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
135 if (!known(temp)) {
136 if (nl_warn)
137 nonlinerr ("angle");
138 nl_fail = !nl_warn;
139 intlfree (temp);
140 return (commagen (1.0,0.0));
141 } else {
142 repart = Re(temp);
143 impart = Im(temp);
144 intlfree (temp);
145 repart = atan2 (impart, repart);
146 if (!radflag)
147 rtod(repart);
148 return (commagen (repart, 0.0));
149 }
150 } else if (((int) intl->left) == lookup ("E")) {
151 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
152 if (!known(temp)) {
153 if (nl_warn)
154 nonlinerr ("E");
155 nl_fail = !nl_warn;
156 intlfree (temp);
157 return (commagen (1.0, 0.0));
158 } else {
159 repart = Re(temp);
160 impart = Im(temp);
161 if (impart > EPSILON)
162 fprintf (stderr, "ideal: E of complex value\n >>>Ignoring imaginary part\n");
163 repart *= 2*PI;
164 return (commagen (cos (repart), sin (repart)));
165 }
166 } else if (((int) intl->left) == lookup ("unit")) {
167 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
168 if (!known(temp)) {
169 if (nl_warn)
170 nonlinerr ("unit");
171 nl_fail = !nl_warn;
172 intlfree (temp);
173 return (commagen (1.0, 0.0));
174 } else {
175 repart = Re(temp);
176 impart = Im(temp);
177 intlfree (temp);
178 if ((modulus = hypot (repart, impart)) < EPSILON)
179 return (commagen (0.0, 0.0));
180 else return (
181 commagen (
182 repart/modulus,
183 impart/modulus
184 )
185 );
186 }
187 } else if (((int) intl->left) == lookup ("sqrt")) {
188 temp = expreval (((EXPRPTR) intl->right)->expr, givennoad);
189 if (!known(temp)) {
190 if (nl_warn)
191 nonlinerr ("sqrt");
192 nl_fail = !nl_warn;
193 intlfree (temp);
194 return (commagen (1.0, 0.0));
195 } else {
196 repart = Re(temp);
197 impart = Im(temp);
198 intlfree (temp);
199 if ((modulus = hypot (repart, impart)) < EPSILON)
200 return (commagen (0.0, 0.0));
201 else {
202 float theta;
203 modulus = sqrt (modulus);
204 theta = 0.5*atan2 (impart,repart);
205 return (
206 commagen (
207 modulus*cos(theta),
208 modulus*sin(theta)
209 )
210 );
211 };
212 }
213 } else {
214 fprintf (stderr, "ideal: unknown function name: %s\n >>>Returning 1.0\n", idprint ((int) intl->left));
215 return (commagen (1.0, 0.0));
216 }
217 break;
218 case '~':
219 incon_warn = FALSE;
220 /* FALL THROUGH TO '=' case */
221 case '=':
222 {
223 DEPPTR depvarwalk;
224 lefttemp = expreval (intl->left, givennoad);
225 righttemp = expreval (intl->right, givennoad);
226 if (nl_fail) {
227 dprintf "Non-linear equation: failure\n");
228 intlfree (lefttemp);
229 intlfree (righttemp);
230 return (commagen (0.0,0.0));
231 }
232 for (depvarwalk = depvarlist;
233 depvarwalk;
234 depvarwalk = depvarwalk->next) {
235 lefttemp->left = (EXPR) depsubst (
236 (DEPPTR) lefttemp->left,
237 (DEPPTR) depvarwalk->var->deplist,
238 depvarwalk->var
239 );
240 lefttemp->right = (EXPR) depsubst (
241 (DEPPTR) lefttemp->right,
242 (DEPPTR) depvarwalk->var->deplist,
243 depvarwalk->var
244 );
245 righttemp->left = (EXPR) depsubst (
246 (DEPPTR) righttemp->left,
247 (DEPPTR) depvarwalk->var->deplist,
248 depvarwalk->var
249 );
250 righttemp->right = (EXPR) depsubst (
251 (DEPPTR) righttemp->right,
252 (DEPPTR) depvarwalk->var->deplist,
253 depvarwalk->var
254 );
255 }
256 dprintf "equating real parts...\n");
257 drek = depadd (
258 (DEPPTR) lefttemp->left, 1.0,
259 (DEPPTR) righttemp->left, -1.0
260 );
261 eqndo (drek, exprn, givennoad);
262 depfree (drek);
263 if (depvarlist) {
264 /* trick: at most one variable became
265 /* dependent by the above processing,
266 /* so only it must be replaced in the
267 /* equation on the imaginary parts */
268 lefttemp->right = (EXPR) depsubst (
269 (DEPPTR) lefttemp->right,
270 (DEPPTR) depvarlist->var->deplist,
271 depvarlist->var
272 );
273 righttemp->right = (EXPR) depsubst (
274 (DEPPTR) righttemp->right,
275 (DEPPTR) depvarlist->var->deplist,
276 depvarlist->var
277 );
278 }
279 dprintf "equating imag parts...\n");
280 drek = depadd (
281 (DEPPTR) lefttemp->right, 1.0,
282 (DEPPTR) righttemp->right, -1.0
283 );
284 eqndo (drek, exprn, givennoad);
285 depfree (drek);
286 intlfree (lefttemp);
287 return (righttemp);
288 }
289 break;
290 case '+':
291 lefttemp = expreval (intl->left, givennoad);
292 righttemp = expreval (intl->right, givennoad);
293 drek = depadd (
294 (DEPPTR) lefttemp->left, 1.0,
295 (DEPPTR) righttemp->left, 1.0
296 );
297 drek2 = depadd (
298 (DEPPTR) lefttemp->right, 1.0,
299 (DEPPTR) righttemp->right, 1.0
300 );
301 intlfree (lefttemp);
302 intlfree (righttemp);
303 return (intlgen (';', (EXPR) drek, (EXPR) drek2));
304 break;
305 case '-':
306 lefttemp = expreval (intl->left, givennoad);
307 righttemp = expreval (intl->right, givennoad);
308 drek = depadd (
309 (DEPPTR) lefttemp->left, 1.0,
310 (DEPPTR) righttemp->left, -1.0
311 );
312 drek2 = depadd (
313 (DEPPTR) lefttemp->right, 1.0,
314 (DEPPTR) righttemp->right, -1.0
315 );
316 intlfree (lefttemp);
317 intlfree (righttemp);
318 return (intlgen (';', (EXPR) drek, (EXPR) drek2));
319 break;
320 case '*':
321 lefttemp = expreval (intl->left, givennoad);
322 righttemp = expreval (intl->right, givennoad);
323 if (known(lefttemp)) {
324 repart = ((DEPPTR) lefttemp->left)->coeff;
325 impart = ((DEPPTR) lefttemp->right)->coeff;
326 intlfree (lefttemp);
327 drek = depadd (
328 (DEPPTR) righttemp->left, repart,
329 (DEPPTR) righttemp->right, -impart
330 );
331 drek2 = depadd (
332 (DEPPTR) righttemp->left, impart,
333 (DEPPTR) righttemp->right, repart
334 );
335 intlfree (righttemp);
336 return (intlgen (';', (EXPR) drek, (EXPR) drek2));
337 } else if (known(righttemp)) {
338 repart = ((DEPPTR) righttemp->left)->coeff;
339 impart = ((DEPPTR) righttemp->right)->coeff;
340 intlfree (righttemp);
341 drek = depadd (
342 (DEPPTR) lefttemp->left, repart,
343 (DEPPTR) lefttemp->right, -impart
344 );
345 drek2 = depadd (
346 (DEPPTR) lefttemp->left, impart,
347 (DEPPTR) lefttemp->right, repart
348 );
349 intlfree (lefttemp);
350 return (intlgen (';', (EXPR) drek, (EXPR) drek2));
351 } else {
352 if (nl_warn)
353 fprintf (stderr, "ideal: multiplication of two unknowns\n >>>Returning 1.0\n");
354 nl_fail = !nl_warn;
355 intlfree (lefttemp);
356 intlfree (righttemp);
357 return (commagen (1.0, 0.0));
358 }
359 break;
360 case '/':
361 lefttemp = expreval (intl->left, givennoad);
362 righttemp = expreval (intl->right, givennoad);
363 if (!known(righttemp)) {
364 if (nl_warn)
365 fprintf (stderr, "ideal: division by an unknown\n >>>Returning 1.0\n");
366 nl_fail = !nl_warn;
367 intlfree (lefttemp);
368 intlfree (righttemp);
369 return (commagen (1.0, 0.0));
370 } else {
371 repart = ((DEPPTR) righttemp->left)->coeff;
372 impart = - ((DEPPTR) righttemp->right)->coeff;
373 modulus = repart*repart + impart*impart;
374 intlfree (righttemp);
375 if (modulus < EPSILON*EPSILON) {
376 fprintf (stderr, "ideal: division by zero\n >>>Returning 1.0\n");
377 intlfree (lefttemp);
378 return (commagen (1.0, 0.0));
379 } else {
380 drek = depadd (
381 (DEPPTR) lefttemp->left, repart/modulus,
382 (DEPPTR) lefttemp->right, -impart/modulus
383 );
384 drek2 = depadd (
385 (DEPPTR) lefttemp->left, impart/modulus,
386 (DEPPTR) lefttemp->right, repart/modulus
387 );
388 intlfree (lefttemp);
389 return (intlgen (';', (EXPR) drek, (EXPR) drek2));
390 }
391 }
392 break;
393 case ',':
394 lefttemp = expreval (intl->left, givennoad);
395 righttemp = expreval (intl->right, givennoad);
396 depfree((DEPPTR) lefttemp->right);
397 depfree((DEPPTR) righttemp->right);
398 temp = intlgen (
399 ';',
400 (EXPR) lefttemp->left,
401 (EXPR) righttemp->left
402 );
403 tryfree(lefttemp);
404 tryfree(righttemp);
405 return (temp);
406 break;
407 case ';':
408 drek = depadd (
409 (DEPPTR) intl->left, 1.0,
410 (DEPPTR) NULL, 0.0
411 );
412 drek2 = depadd (
413 (DEPPTR) intl->right, 1.0,
414 (DEPPTR) NULL, 0.0
415 );
416 return (intlgen (';', (EXPR) drek, (EXPR) drek2));
417 case '^':
418 return (expreval (intl->right, givennoad));
419 default:
420 fprintf (stderr, "ideal: unknown operator: %c\n >>>Returning 1.0\n", intl->oper);
421 return (commagen (1.0, 0.0));
422 break;
423 }
424 }
425
eqndo(deplist,eqn,givennoad)426 void eqndo (deplist, eqn, givennoad)
427 DEPPTR deplist;
428 EXPR eqn;
429 NOADPTR givennoad;
430 {
431 /* when called, equation system says deplist == 0 */
432 if (!deplist->next && !deplist->var) {
433 if (fabs (deplist->coeff) > EPSILON) {
434 if (incon_warn) {
435 fprintf (stderr, "ideal: inconsistent equation in %s named %s\n",
436 idprint (givennoad->defnode->parm->name),
437 idprint (givennoad->defnode->name)
438 );
439 exprprint (((INTLPTR) eqn)->left);
440 fprintf (stderr, "=");
441 exprprint (eqn);
442 fprintf (stderr, "\n");
443 }
444 dprintf "Inconsistent equation\n");
445 } else
446 dprintf "Redundant equation\n");
447 }
448 else {
449 DEPPTR curmax;
450 float maxcoeff;
451 DEPPTR depvarwalk;
452 DEPPTR listwalk;
453 maxcoeff = -1;
454 /* find variable whose coefficient is largest in absolute value */
455 for (listwalk = deplist;
456 listwalk;
457 listwalk = listwalk->next)
458 if (listwalk->var && (maxcoeff < fabs (listwalk->coeff))) {
459 maxcoeff = fabs (listwalk->coeff);
460 curmax = listwalk;
461 }
462 /* get that variable represented in terms of the others */
463 listwalk = depadd (
464 curmax->var->deplist, 1.0,
465 deplist, -1.0/curmax->coeff
466 );
467 depfree (curmax->var->deplist);
468 curmax->var->deplist = listwalk;
469 /* put it on a list of dependent variables
470 /* replace occurrences of it in other dependent variables */
471 if (!depvarlist) {
472 depvarlist = depgen (curmax->var, 0.0);
473 }
474 else {
475 DEPPTR newhead;
476 for (depvarwalk = depvarlist;
477 depvarwalk;
478 depvarwalk = depvarwalk->next) {
479 depvarwalk->var->deplist = depsubst (
480 depvarwalk->var->deplist,
481 curmax->var->deplist,
482 curmax->var
483 );
484 }
485 newhead = depgen (curmax->var, 0.0);
486 newhead->next = depvarlist;
487 depvarlist = newhead;
488 }
489 }
490 }
491
depvarclean()492 void depvarclean ()
493 {
494 /* clean known variables out of the dependent variable list */
495 DEPPTR prevdep, depvarwalk;
496 DEPNODE nuhead;
497 prevdep = &nuhead;
498 prevdep->next = depvarwalk = depvarlist;
499 while (depvarwalk) {
500 if (!depvarwalk->var->deplist->var) {
501 dprintf "Removing %s(%s) = %f from dependent variable list\n",
502 ISREAL(depvarwalk->var)?"re":"im",
503 idprint (THENAME(depvarwalk->var)),
504 depvarwalk->var->deplist->coeff);
505 prevdep->next = depvarwalk->next;
506 tryfree(depvarwalk);
507 depvarwalk = prevdep->next;
508 } else {
509 prevdep = depvarwalk;
510 depvarwalk = depvarwalk->next;
511 }
512 }
513 depvarlist = nuhead.next;
514 }
515
reqneval(noadtree)516 void reqneval (noadtree)
517 NOADPTR noadtree;
518 {
519 STMTPTR slist[2];
520 STMTPTR eqnwalk;
521 int i;
522 if (!noadtree)
523 return;
524 nl_warn = FALSE;
525 slist[0] = noadtree->defnode->parm->stmtlist;
526 slist[1] = findbox (noadtree->defnode->parm->name,FALSE)->stmtlist;
527 for (i = 0; i < 2; i ++)
528 for (eqnwalk = nextstmt ('=', slist[i]);
529 eqnwalk;
530 eqnwalk = nextstmt ('=', eqnwalk->next)) {
531 INTLPTR junk;
532 nl_fail = FALSE;
533 junk = expreval ((EXPR) eqnwalk->stmt, noadtree);
534 intlfree (junk);
535 if (nl_fail) {
536 EQNPTR nueqn;
537 nueqn = eqngen (
538 (EXPR) eqnwalk->stmt,
539 noadtree
540 );
541 nueqn->next = nl_eqns;
542 nl_eqns = nueqn;
543 nl_fail = FALSE;
544 }
545 depvarclean ();
546 incon_warn = TRUE;
547 }
548 reqneval (noadtree->son);
549 reqneval (noadtree->brother);
550 }
551
eqneval(noadtree)552 void eqneval (noadtree)
553 NOADPTR noadtree;
554 {
555 if (when_bug & 04) bug_on;
556 reqneval (noadtree);
557 bug_off;
558 }
559
nl_eval()560 void nl_eval ()
561 {
562 static boolean nl_succ;
563 INTLPTR junk;
564 {
565 EQNPTR nl_prev, nl_curr, nl_temp;
566 if (when_bug & 010) bug_on;
567 nl_prev = nl_curr = nl_eqns;
568 nl_temp = NULL;
569 while (nl_curr) {
570 nl_curr = nl_prev->next;
571 nl_prev->next = nl_temp;
572 nl_temp = nl_prev;
573 nl_prev = nl_curr;
574 }
575 nl_eqns = nl_temp;
576 nl_succ = TRUE;
577 }
578 while (nl_eqns && nl_succ) {
579 EQNPTR prev_eqn, nl_walk;
580 EQNNODE dummy_eqn;
581 dprintf "Retrying nonlinear equations\n");
582 prev_eqn = &dummy_eqn;
583 prev_eqn->next = nl_walk = nl_eqns;
584 nl_succ = FALSE;
585 while (nl_walk) {
586 nl_fail = FALSE;
587 junk = expreval (nl_walk->eqn, nl_walk->noad);
588 intlfree (junk);
589 depvarclean ();
590 if (!nl_fail) {
591 prev_eqn->next = nl_walk->next;
592 tryfree(nl_walk);
593 nl_walk = prev_eqn->next;
594 nl_succ = TRUE;
595 } else {
596 prev_eqn = nl_walk;
597 nl_walk = nl_walk->next;
598 }
599 }
600 nl_eqns = dummy_eqn.next;
601 }
602 if (nl_eqns) {
603 EQNPTR nl_walk, nl_next;
604 dprintf "Nonlinear failure\n");
605 nl_warn = TRUE;
606 for (nl_walk = nl_eqns;
607 nl_walk;
608 nl_walk = nl_next) {
609 junk = expreval (nl_walk->eqn, nl_walk->noad);
610 intlfree (junk);
611 depvarclean ();
612 nl_next = nl_walk->next;
613 tryfree(nl_walk);
614 }
615 }
616 bug_off;
617 }
618
depvarkill()619 void depvarkill ()
620 {
621 /* remove all unknown variables from depvarlist ...
622 no chance for them to be determined now */
623 if (!depvarlist)
624 return;
625 if (when_bug & 020)
626 fprintf (stderr, "killing depvarlist\n");
627 depfree (depvarlist);
628 depvarlist = NULL;
629 }
630