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