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 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 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 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 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 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 552 void eqneval (noadtree) 553 NOADPTR noadtree; 554 { 555 if (when_bug & 04) bug_on; 556 reqneval (noadtree); 557 bug_off; 558 } 559 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 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