#include "ideal.h" #include "y.tab.h" void nonlinerr (funcname) char *funcname; { fprintf (stderr, "ideal: %s() of unknown\n >>>Returning 1.0\n", funcname); } static DEPPTR depvarlist = NULL; static boolean incon_warn = TRUE; static boolean nl_warn; boolean nl_fail; static EQNPTR nl_eqns = NULL; INTLPTR expreval (exprn, givennoad) EXPR exprn; NOADPTR givennoad; { /* This routine returns an INTLPTR whose operator is ';'--a promoted commanode containing the dependency list representing the real part in its left field, the imag part in its right */ register INTLPTR intl; register EXTLPTR extl; if (!exprn) return (commagen (0.0, 0.0)); if (((EXTLPTR)exprn)->leaf) { extl = (EXTLPTR) exprn; dprintf "At a leaf of kind %d\n", extl->kind); switch (extl->kind) { case PATH: return (pathfind (extl->info.path, givennoad)); break; case CONST: return (commagen (extl->info.const, 0.0)); break; } } intl = (INTLPTR) exprn; if (intl->oper == NAME) { dprintf "Looking for a function named %s\n", idprint ((int) intl->left)); } else { dprintf "At an internal node with operator %c\n", intl->oper); } switch (intl->oper) { INTLPTR lefttemp, righttemp, temp, temp2; DEPPTR drek, drek2; float repart, impart, modulus; case NAME: if (((int) intl->left) == lookup ("re")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); depfree ((DEPPTR) temp->right); temp->right = (EXPR) depgen ((VARPTR) NULL, 0.0); return (temp); } else if (((int) intl->left) == lookup ("im")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); depfree ((DEPPTR) temp->left); temp->left = temp->right; temp->right = (EXPR) depgen ((VARPTR) NULL, 0.0); return (temp); } else if (((int) intl->left) == lookup ("conj")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); temp2 = intlgen ( ';', (EXPR) depadd ( (DEPPTR) NULL, 0.0, (DEPPTR) temp->left, 1.0 ), (EXPR) depadd ( (DEPPTR) NULL, 0.0, (DEPPTR) temp->right, -1.0 ) ); intlfree (temp); return (temp2); } else if (((int) intl->left) == lookup ("abs")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); if (!known(temp)) { if (nl_warn) nonlinerr ("abs"); nl_fail = !nl_warn; intlfree (temp); return (commagen (1.0, 0.0)); } else { repart = Re(temp); impart = Im(temp); intlfree (temp); return (commagen (hypot (repart, impart), 0.0)); } } else if (((int) intl->left) == lookup ("cis")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); if (!known(temp)) { if (nl_warn) nonlinerr ("cis"); nl_fail = !nl_warn; intlfree (temp); return (commagen (1.0, 0.0)); } else { repart = Re(temp); impart = Im(temp); if (!radflag) { dtor(repart); dtor(impart); } intlfree (temp); if (impart > EPSILON) fprintf (stderr, "ideal: cis of complex value\n >>>Ignoring imaginary part\n"); return (commagen (cos (repart), sin (repart))); } } else if (((int) intl->left) == lookup ("int")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); if (!known(temp)) { if (nl_warn) nonlinerr ("int"); nl_fail = !nl_warn; intlfree (temp); return (commagen (1.0,0.0)); } else { double intpart; repart = Re(temp); impart = Im(temp); intlfree (temp); if (impart > EPSILON) fprintf (stderr, "ideal: int of complex value\n >>>Ignoring imaginary part\n"); modf (repart, &intpart); return (commagen ((float) intpart, 0.0)); } } else if (((int) intl->left) == lookup ("atan2") || ((int) intl->left) == lookup ("angle")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); if (!known(temp)) { if (nl_warn) nonlinerr ("angle"); nl_fail = !nl_warn; intlfree (temp); return (commagen (1.0,0.0)); } else { repart = Re(temp); impart = Im(temp); intlfree (temp); repart = atan2 (impart, repart); if (!radflag) rtod(repart); return (commagen (repart, 0.0)); } } else if (((int) intl->left) == lookup ("E")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); if (!known(temp)) { if (nl_warn) nonlinerr ("E"); nl_fail = !nl_warn; intlfree (temp); return (commagen (1.0, 0.0)); } else { repart = Re(temp); impart = Im(temp); if (impart > EPSILON) fprintf (stderr, "ideal: E of complex value\n >>>Ignoring imaginary part\n"); repart *= 2*PI; return (commagen (cos (repart), sin (repart))); } } else if (((int) intl->left) == lookup ("unit")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); if (!known(temp)) { if (nl_warn) nonlinerr ("unit"); nl_fail = !nl_warn; intlfree (temp); return (commagen (1.0, 0.0)); } else { repart = Re(temp); impart = Im(temp); intlfree (temp); if ((modulus = hypot (repart, impart)) < EPSILON) return (commagen (0.0, 0.0)); else return ( commagen ( repart/modulus, impart/modulus ) ); } } else if (((int) intl->left) == lookup ("sqrt")) { temp = expreval (((EXPRPTR) intl->right)->expr, givennoad); if (!known(temp)) { if (nl_warn) nonlinerr ("sqrt"); nl_fail = !nl_warn; intlfree (temp); return (commagen (1.0, 0.0)); } else { repart = Re(temp); impart = Im(temp); intlfree (temp); if ((modulus = hypot (repart, impart)) < EPSILON) return (commagen (0.0, 0.0)); else { float theta; modulus = sqrt (modulus); theta = 0.5*atan2 (impart,repart); return ( commagen ( modulus*cos(theta), modulus*sin(theta) ) ); }; } } else { fprintf (stderr, "ideal: unknown function name: %s\n >>>Returning 1.0\n", idprint ((int) intl->left)); return (commagen (1.0, 0.0)); } break; case '~': incon_warn = FALSE; /* FALL THROUGH TO '=' case */ case '=': { DEPPTR depvarwalk; lefttemp = expreval (intl->left, givennoad); righttemp = expreval (intl->right, givennoad); if (nl_fail) { dprintf "Non-linear equation: failure\n"); intlfree (lefttemp); intlfree (righttemp); return (commagen (0.0,0.0)); } for (depvarwalk = depvarlist; depvarwalk; depvarwalk = depvarwalk->next) { lefttemp->left = (EXPR) depsubst ( (DEPPTR) lefttemp->left, (DEPPTR) depvarwalk->var->deplist, depvarwalk->var ); lefttemp->right = (EXPR) depsubst ( (DEPPTR) lefttemp->right, (DEPPTR) depvarwalk->var->deplist, depvarwalk->var ); righttemp->left = (EXPR) depsubst ( (DEPPTR) righttemp->left, (DEPPTR) depvarwalk->var->deplist, depvarwalk->var ); righttemp->right = (EXPR) depsubst ( (DEPPTR) righttemp->right, (DEPPTR) depvarwalk->var->deplist, depvarwalk->var ); } dprintf "equating real parts...\n"); drek = depadd ( (DEPPTR) lefttemp->left, 1.0, (DEPPTR) righttemp->left, -1.0 ); eqndo (drek, exprn, givennoad); depfree (drek); if (depvarlist) { /* trick: at most one variable became /* dependent by the above processing, /* so only it must be replaced in the /* equation on the imaginary parts */ lefttemp->right = (EXPR) depsubst ( (DEPPTR) lefttemp->right, (DEPPTR) depvarlist->var->deplist, depvarlist->var ); righttemp->right = (EXPR) depsubst ( (DEPPTR) righttemp->right, (DEPPTR) depvarlist->var->deplist, depvarlist->var ); } dprintf "equating imag parts...\n"); drek = depadd ( (DEPPTR) lefttemp->right, 1.0, (DEPPTR) righttemp->right, -1.0 ); eqndo (drek, exprn, givennoad); depfree (drek); intlfree (lefttemp); return (righttemp); } break; case '+': lefttemp = expreval (intl->left, givennoad); righttemp = expreval (intl->right, givennoad); drek = depadd ( (DEPPTR) lefttemp->left, 1.0, (DEPPTR) righttemp->left, 1.0 ); drek2 = depadd ( (DEPPTR) lefttemp->right, 1.0, (DEPPTR) righttemp->right, 1.0 ); intlfree (lefttemp); intlfree (righttemp); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); break; case '-': lefttemp = expreval (intl->left, givennoad); righttemp = expreval (intl->right, givennoad); drek = depadd ( (DEPPTR) lefttemp->left, 1.0, (DEPPTR) righttemp->left, -1.0 ); drek2 = depadd ( (DEPPTR) lefttemp->right, 1.0, (DEPPTR) righttemp->right, -1.0 ); intlfree (lefttemp); intlfree (righttemp); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); break; case '*': lefttemp = expreval (intl->left, givennoad); righttemp = expreval (intl->right, givennoad); if (known(lefttemp)) { repart = ((DEPPTR) lefttemp->left)->coeff; impart = ((DEPPTR) lefttemp->right)->coeff; intlfree (lefttemp); drek = depadd ( (DEPPTR) righttemp->left, repart, (DEPPTR) righttemp->right, -impart ); drek2 = depadd ( (DEPPTR) righttemp->left, impart, (DEPPTR) righttemp->right, repart ); intlfree (righttemp); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); } else if (known(righttemp)) { repart = ((DEPPTR) righttemp->left)->coeff; impart = ((DEPPTR) righttemp->right)->coeff; intlfree (righttemp); drek = depadd ( (DEPPTR) lefttemp->left, repart, (DEPPTR) lefttemp->right, -impart ); drek2 = depadd ( (DEPPTR) lefttemp->left, impart, (DEPPTR) lefttemp->right, repart ); intlfree (lefttemp); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); } else { if (nl_warn) fprintf (stderr, "ideal: multiplication of two unknowns\n >>>Returning 1.0\n"); nl_fail = !nl_warn; intlfree (lefttemp); intlfree (righttemp); return (commagen (1.0, 0.0)); } break; case ELEWISE: lefttemp = expreval (intl->left, givennoad); righttemp = expreval (intl->right, givennoad); if (known(lefttemp)) { repart = ((DEPPTR) lefttemp->left)->coeff; impart = ((DEPPTR) lefttemp->right)->coeff; intlfree (lefttemp); drek = depadd ( (DEPPTR) righttemp->left, repart, (DEPPTR) NULL, 0.0 ); drek2 = depadd ( (DEPPTR) righttemp->right, impart, (DEPPTR) NULL, 0.0 ); intlfree (righttemp); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); } else if (known(righttemp)) { repart = ((DEPPTR) righttemp->left)->coeff; impart = ((DEPPTR) righttemp->right)->coeff; intlfree (righttemp); drek = depadd ( (DEPPTR) lefttemp->left, repart, (DEPPTR) NULL, 0.0 ); drek2 = depadd ( (DEPPTR) lefttemp->right, impart, (DEPPTR) NULL, 0.0 ); intlfree (lefttemp); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); } else { if (nl_warn) fprintf (stderr, "ideal: multiplication of two unknowns\n >>>Returning 1.0\n"); nl_fail = !nl_warn; intlfree (lefttemp); intlfree (righttemp); return (commagen (1.0, 0.0)); } break; case '/': lefttemp = expreval (intl->left, givennoad); righttemp = expreval (intl->right, givennoad); if (!known(righttemp)) { if (nl_warn) fprintf (stderr, "ideal: division by an unknown\n >>>Returning 1.0\n"); nl_fail = !nl_warn; intlfree (lefttemp); intlfree (righttemp); return (commagen (1.0, 0.0)); } else { repart = ((DEPPTR) righttemp->left)->coeff; impart = - ((DEPPTR) righttemp->right)->coeff; modulus = repart*repart + impart*impart; intlfree (righttemp); if (modulus < EPSILON*EPSILON) { fprintf (stderr, "ideal: division by zero\n >>>Returning 1.0\n"); intlfree (lefttemp); return (commagen (1.0, 0.0)); } else { drek = depadd ( (DEPPTR) lefttemp->left, repart/modulus, (DEPPTR) lefttemp->right, -impart/modulus ); drek2 = depadd ( (DEPPTR) lefttemp->left, impart/modulus, (DEPPTR) lefttemp->right, repart/modulus ); intlfree (lefttemp); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); } } break; case ',': lefttemp = expreval (intl->left, givennoad); righttemp = expreval (intl->right, givennoad); depfree((DEPPTR) lefttemp->right); depfree((DEPPTR) righttemp->right); temp = intlgen ( ';', (EXPR) lefttemp->left, (EXPR) righttemp->left ); tryfree(lefttemp); tryfree(righttemp); return (temp); break; case ';': drek = depadd ( (DEPPTR) intl->left, 1.0, (DEPPTR) NULL, 0.0 ); drek2 = depadd ( (DEPPTR) intl->right, 1.0, (DEPPTR) NULL, 0.0 ); return (intlgen (';', (EXPR) drek, (EXPR) drek2)); case '^': return (expreval (intl->right, givennoad)); default: fprintf (stderr, "ideal: unknown operator: %c\n >>>Returning 1.0\n", intl->oper); return (commagen (1.0, 0.0)); break; } } void eqndo (deplist, eqn, givennoad) DEPPTR deplist; EXPR eqn; NOADPTR givennoad; { /* when called, equation system says deplist == 0 */ if (!deplist->next && !deplist->var) { if (fabs (deplist->coeff) > EPSILON) { if (incon_warn) { fprintf (stderr, "ideal: inconsistent equation in %s named %s\n", idprint (givennoad->defnode->parm->name), idprint (givennoad->defnode->name) ); exprprint (((INTLPTR) eqn)->left); fprintf (stderr, "="); exprprint (eqn); fprintf (stderr, "\n"); } dprintf "Inconsistent equation\n"); } else dprintf "Redundant equation\n"); } else { DEPPTR curmax; float maxcoeff; DEPPTR depvarwalk; DEPPTR listwalk; maxcoeff = -1; /* find variable whose coefficient is largest in absolute value */ for (listwalk = deplist; listwalk; listwalk = listwalk->next) if (listwalk->var && (maxcoeff < fabs (listwalk->coeff))) { maxcoeff = fabs (listwalk->coeff); curmax = listwalk; } /* get that variable represented in terms of the others */ listwalk = depadd ( curmax->var->deplist, 1.0, deplist, -1.0/curmax->coeff ); depfree (curmax->var->deplist); curmax->var->deplist = listwalk; /* put it on a list of dependent variables /* replace occurrences of it in other dependent variables */ if (!depvarlist) { depvarlist = depgen (curmax->var, 0.0); } else { DEPPTR newhead; for (depvarwalk = depvarlist; depvarwalk; depvarwalk = depvarwalk->next) { depvarwalk->var->deplist = depsubst ( depvarwalk->var->deplist, curmax->var->deplist, curmax->var ); } newhead = depgen (curmax->var, 0.0); newhead->next = depvarlist; depvarlist = newhead; } } } void depvarclean () { /* clean known variables out of the dependent variable list */ DEPPTR prevdep, depvarwalk; DEPNODE nuhead; prevdep = &nuhead; prevdep->next = depvarwalk = depvarlist; while (depvarwalk) { if (!depvarwalk->var->deplist->var) { dprintf "Removing %s(%s) = %f from dependent variable list\n", ISREAL(depvarwalk->var)?"re":"im", idprint (THENAME(depvarwalk->var)), depvarwalk->var->deplist->coeff); prevdep->next = depvarwalk->next; tryfree(depvarwalk); depvarwalk = prevdep->next; } else { prevdep = depvarwalk; depvarwalk = depvarwalk->next; } } depvarlist = nuhead.next; } void reqneval (noadtree) NOADPTR noadtree; { STMTPTR slist[2]; STMTPTR eqnwalk; int i; if (!noadtree) return; nl_warn = FALSE; slist[0] = noadtree->defnode->parm->stmtlist; slist[1] = findbox (noadtree->defnode->parm->name,FALSE)->stmtlist; for (i = 0; i < 2; i ++) for (eqnwalk = nextstmt ('=', slist[i]); eqnwalk; eqnwalk = nextstmt ('=', eqnwalk->next)) { INTLPTR junk; nl_fail = FALSE; junk = expreval ((EXPR) eqnwalk->stmt, noadtree); intlfree (junk); if (nl_fail) { EQNPTR nueqn; nueqn = eqngen ( (EXPR) eqnwalk->stmt, noadtree ); nueqn->next = nl_eqns; nl_eqns = nueqn; nl_fail = FALSE; } depvarclean (); incon_warn = TRUE; } reqneval (noadtree->son); reqneval (noadtree->brother); } void eqneval (noadtree) NOADPTR noadtree; { if (when_bug & 04) bug_on; reqneval (noadtree); bug_off; } void nl_eval () { static boolean nl_succ; INTLPTR junk; { EQNPTR nl_prev, nl_curr, nl_temp; if (when_bug & 010) bug_on; nl_prev = nl_curr = nl_eqns; nl_temp = NULL; while (nl_curr) { nl_curr = nl_prev->next; nl_prev->next = nl_temp; nl_temp = nl_prev; nl_prev = nl_curr; } nl_eqns = nl_temp; nl_succ = TRUE; } while (nl_eqns && nl_succ) { EQNPTR prev_eqn, nl_walk; EQNNODE dummy_eqn; dprintf "Retrying nonlinear equations\n"); prev_eqn = &dummy_eqn; prev_eqn->next = nl_walk = nl_eqns; nl_succ = FALSE; while (nl_walk) { nl_fail = FALSE; junk = expreval (nl_walk->eqn, nl_walk->noad); intlfree (junk); depvarclean (); if (!nl_fail) { prev_eqn->next = nl_walk->next; tryfree(nl_walk); nl_walk = prev_eqn->next; nl_succ = TRUE; } else { prev_eqn = nl_walk; nl_walk = nl_walk->next; } } nl_eqns = dummy_eqn.next; } if (nl_eqns) { EQNPTR nl_walk, nl_next; dprintf "Nonlinear failure\n"); nl_warn = TRUE; for (nl_walk = nl_eqns; nl_walk; nl_walk = nl_next) { junk = expreval (nl_walk->eqn, nl_walk->noad); intlfree (junk); depvarclean (); nl_next = nl_walk->next; tryfree(nl_walk); } } bug_off; } void depvarkill () { /* remove all unknown variables from depvarlist ... no chance for them to be determined now */ if (!depvarlist) return; if (when_bug & 020) fprintf (stderr, "killing depvarlist\n"); depfree (depvarlist); depvarlist = NULL; }