Corrected a long-standing error in which ending text with a literal
[xcircuit.git] / spiceparser / equations.c
blob1806a705893947c40de7e66a538e30ca0532c146
1 /********************
2 This file is part of the software library CADLIB written by Conrad Ziesler
3 Copyright 2003, Conrad Ziesler, all rights reserved.
5 *************************
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 ******************/
21 /* equations.c, equation evaluation and parsing
23 Conrad Ziesler
27 #include <ctype.h>
28 #include <math.h>
29 #include <stdio.h>
31 #include "debug.h"
32 #include "equations.h"
34 int __debug_eqn__=0;
36 #ifndef STRIP_DEBUGGING
37 #define D(level,a) do { if(__debug_eqn__>(level)) a; } while(0)
38 #else
39 #define D(level,a)
40 #endif
43 #define uchar unsigned char
47 static eqntoken_t *eqntoken_next(eqntoken_t *p);
50 static void debug_eqn(FILE *fp, char *str, eqn_t *eqn)
52 eqntoken_t *eqt;
53 if(eqn==NULL)return;
54 fprintf(fp," %s=",str);
55 if(eqn->eqn!=NULL)
57 char s=' ';
58 for(eqt=eqn->eqn;!OP_ISEND(eqt->op);eqt=eqntoken_next(eqt))
60 switch(eqt->op)
62 case OPlitm:
63 s='-';
64 case OPlit:
65 fprintf(fp,"(%c%s)",s,eqt->z.lit.lit); break;
66 case OPvalm:
67 s='-';
68 case OPval:
69 fprintf(fp,"(%c%g)",s,eqt->z.val.x); break;
70 case OPlitvm:
71 s='-';
72 case OPlitv:
73 fprintf(fp,"(%c%p->%g)",s,eqt->z.pval.xp,eqt->z.pval.xp[0]); break;
75 case OPadd: fprintf(fp,"+"); break;
76 case OPsub: fprintf(fp,"-"); break;
77 case OPmul: fprintf(fp,"*"); break;
78 case OPdiv: fprintf(fp,"/"); break;
79 case OPexp: fprintf(fp,"^"); break;
80 case OPopen: fprintf(fp,"{"); break;
81 case OPclose: fprintf(fp,"}"); break;
82 case OPnorm: fprintf(fp,":norm:"); break; /* A norm B == norm_B(A) */
83 default: break;
86 if(eqt->op==OPeol_const)
87 fprintf(fp,"%g",eqt->z.val.x);
88 if(eqt->op==OPeol_valp)
89 fprintf(fp,"%p->%g",eqt->z.pval.xp, (eqt->z.pval.xp==NULL)?0.0:*(eqt->z.pval.xp));
96 static void *eqn_malloc(int size)
98 void *p=calloc(1,size);
99 return p;
101 static char *eqn_strdup(char *p)
103 return strdup(p);
106 static void eqn_free(void *p)
107 { if(p!=NULL) free(p); }
111 static eqntoken_t *eqntoken_next(eqntoken_t *p)
113 eqntoken_t *r;
114 switch(p->op)
116 case OPeolist: case OPeol_const: case OPeol_valp: case OPeol_qty:
117 r=p; break;
118 case OPval: case OPvalm: r= (void *) (((char *)p)+(sizeof(itemop_t)+sizeof(eqntokenval_t) )) ; break;
119 case OPlit: case OPlitm:
120 case OPlitv: case OPlitvm: r= (void *) (((char *)p)+(sizeof(itemop_t)+sizeof(eqntokenlit_t) )) ; break;
121 default: r= (void *) (((char *)p)+(sizeof(itemop_t))) ; break;
123 return r;
126 #ifdef UNUSED
128 static eqntoken_t *eqntoken_last(eqntoken_t *p)
130 eqntoken_t *r=p;
132 while(!OP_ISEND(p->op))p=eqntoken_next(p);
134 switch(p->op)
136 case OPeol_const:
137 r= (void *) (((char *)p)+(sizeof(itemop_t)+sizeof(eqntokenval_t))) ; break;
138 case OPeol_valp:
139 r= (void *) (((char *)p)+(sizeof(itemop_t)+sizeof(eqntokenpval_t))) ; break;
140 case OPeolist:
141 r= (void *) (((char *)p)+(sizeof(itemop_t))) ; break;
142 case OPeol_qty:
143 assert(0);
144 default:
145 assert(0);
147 return r;
149 #endif
154 static inline char t_peek(eqntop_t *top)
156 if(top->n_unget>0) { return top->unget[top->n_unget-1]; }
157 return top->str[0];
160 static inline char t_eat(eqntop_t *top)
162 char c;
163 c=t_peek(top);
164 if(top->n_unget>0) top->n_unget++;
165 else top->str++;
166 return c;
169 static inline void t_uneat(eqntop_t *top, char a)
171 if(top->n_unget>(sizeof(top->unget)-1)) assert(0);
172 else
174 top->unget[top->n_unget]=a;
175 top->n_unget++;
182 static float float_parse(eqntop_t *top)
184 float v=0.0,sign=1,sig=0,rest=0,digit,place,exp;
185 uchar a;
187 while(isspace(t_peek(top))) t_eat(top);
188 a=t_peek(top);
189 if(a=='-') { sign=-1; t_eat(top); }
190 else if(a=='+') { sign=1; t_eat(top); }
192 if(t_peek(top)=='.') { sig=0; t_eat(top); }
193 else
195 sig=0;
196 while(isdigit(t_peek(top)))
197 { sig*=10; sig+=(t_eat(top)-'0'); }
200 if(t_peek(top)=='.') t_eat(top);
201 if(isdigit(t_peek(top)))
203 place=0.1;
204 while(isdigit(t_peek(top)))
205 { digit=(t_eat(top)-'0'); rest+=(digit*place); place*=0.1; }
208 v=(sig+rest)*sign;
210 a=t_peek(top);
212 switch(a)
214 case 'f': t_eat(top); v*=1e-15; break;
215 case 'p': t_eat(top); v*=1e-12; break;
216 case 'n': t_eat(top); v*=1e-9; break;
217 case 'u': t_eat(top); v*=1e-6; break;
218 case 'm': t_eat(top); v*=1e-3; break;
219 case 'x': t_eat(top); v*=1e6; break;
220 case 'e': t_eat(top);
221 exp=0;
222 while(isdigit(t_peek(top)))
223 { exp*=10; exp+=(t_eat(top)-'0'); }
224 v*=pow(10.0,exp);
225 break;
227 if(isalpha(t_peek(top))) /* eat trailing alphanumerics (ie. 15fF <- F is extra alpha) */
228 t_eat(top);
230 D(10,fprintf(stdout,"NUM[%e]\n",v));
231 return v;
237 static int eqntoken_parse(eqntop_t *top)
239 eqntoken_t *node;
240 int i,inquote=0;
241 char c;
243 c=t_peek(top);
244 if(top->nodep>=MAXEQNSTACK)return 0;
246 node=top->nodes+top->nodep++;
247 node->z.lit.lit=NULL;
248 switch(c)
250 case EQN_DELIM: t_eat(top); /* fall through */
251 case 0: node->op=OPeolist; return 0;
252 case '+': node->op=OPadd; t_eat(top); top->opc++; break;
253 case '-': node->op=OPsub; t_eat(top); top->opc++; break;
254 case '*': node->op=OPmul; t_eat(top); top->opc++; break;
255 case '/': node->op=OPdiv; t_eat(top); top->opc++; break;
256 case '^': node->op=OPexp; t_eat(top); top->opc++; break;
257 case '(': node->op=OPopen; t_eat(top); top->opc++; break;
258 case ')': node->op=OPclose; t_eat(top); top->opc++; break;
259 case '$': node->op=OPnorm ; t_eat(top); top->opc++; break;
261 default:
263 if(isdigit(c))
265 node->z.val.x=float_parse(top);
266 node->op=OPval;
268 else
270 char *p;
271 p=alloca(strlen(top->str)+sizeof(top->unget)+16);
273 if(c=='{') { inquote++; t_eat(top); c=t_peek(top); }
275 for(i=0;c!=0;c=t_peek(top),i++)
277 if(!inquote) { if(strchr("\'+-/*^\n\r()",c)!=NULL)break; }
278 if(inquote && (c=='{')) inquote++;
279 if(inquote && (c=='}')) { inquote--; if(inquote==0) { t_eat(top); break; }}
281 p[i]=t_eat(top);
282 p[i+1]=0;
284 D(10,fprintf(stdout,"LIT[%s]\n",p));
285 node->z.lit.lit=eqn_strdup(p);
286 node->op=OPlit;
288 top->litc++;
289 break;
291 return 1;
295 /* this code assumes dependencies have already been cached, see eqnlist.c */
296 static int token_value(eqntoken_t *token, float *v)
298 int r=0;
299 float a=0.0;
300 switch(token->op)
302 case OPeol_valp:
303 if(token->z.pval.xp!=NULL)
304 a=*(token->z.pval.xp);
305 else r++;
306 D(3,fprintf(stdout,"loaded pval=(%p->%g)\n",token->z.pval.xp,a));
307 break;
308 case OPeol_const: a=token->z.val.x; break;
309 case OPlit: a= 0.0; r++; break;
310 case OPlitm: a=-0.0; r++; break;
311 case OPlitv: a= *(token->z.pval.xp); break;
312 case OPlitvm: a=-(*(token->z.pval.xp)); break;
313 case OPval: a= token->z.val.x; break;
314 case OPvalm: a=-(token->z.val.x); break;
315 default: assert(0);
317 *v=a;
318 return r;
320 #define PUSH(a) do { if(data->stackp<MAXSTACK) data->stack[data->stackp++]=(a); \
321 D(5,fprintf(stdout,"Push(%g) new stackp=%i\n",(a),data->stackp)); } while (0)
322 #define POP(a) do { if(data->stackp>0) (a)=data->stack[--(data->stackp)]; \
323 D(5,fprintf(stdout,"Pop(%g) new stackp=%i\n",(a),data->stackp)); } while (0)
325 static int eqntoken_eval(eqneval_t *data)
327 int loop,toend;
328 int r=0;
329 float a,b;
330 eqntoken_t *token,*lp;
331 int i;
333 static const int precidence[]=OPprecidence;
334 static const int type[]=OPtype;
336 int types[4]={ -10 , -10 , -10 , -10 };
337 int precs[4]={ -10, -10 , -10, -10 };
338 int stack_entry=data->stackp;
340 toend=0;
343 itemop_t ops[4]={OPeolist, OPeolist, OPeolist, OPeolist};
344 eqntoken_t *tokens[4]={ NULL, NULL, NULL, NULL };
346 loop=0;
347 token=data->list;
348 data->list=eqntoken_next(data->list);
349 if(data->list==token)
350 toend++; /* catch endings */
352 for(lp=token,i=0;(i<4);i++,lp=eqntoken_next(lp)) /* lookahead */
354 ops[i]=lp->op;
355 tokens[i]=lp;
356 assert(ops[i]>=0);
357 assert((ops[i]*sizeof(int))<sizeof(precidence));
358 assert((ops[i]*sizeof(int))<sizeof(type));
359 types[i]=type[ops[i]];
360 precs[i]=precidence[ops[i]];
361 if(OP_ISEND(lp->op)) /* we did last token, stop */
362 break;
365 D(8,fprintf(stdout,"Eval: lookahead = "));
366 for(i=0;i<4;i++)
368 char *syms[]=OPsyms;
369 D(8,fprintf(stdout," %s,t%i,p%i" ,syms[ops[i]],types[i],precs[i]));
371 D(8,fprintf(stdout,"\n"));
373 if(OP_ISANYV(token->op))
375 r=token_value(token,&a);
376 D(5,fprintf(stdout,"Value: %g\n",a));
377 PUSH(a);
380 else
382 if(types[0]==2) /* binary op */
384 a=1.0; b=1.0;
385 POP(a);
387 if( /* literal and higher precidence operator after literal not () */
388 (OP_ISV(ops[1]) && (precs[2] > precs[0]) && (types[2]!=4) ) ||
390 /* operator (ie open/close) which has a higher precidence */
391 ((types[1]>=2) && (precs[1] > precs[0]))
394 { /* SHIFT */
395 D(6,fprintf(stdout,"Shift: recursion\n"));
397 eqntoken_eval(data);
400 POP(b);
402 else /* REDUCE */
404 D(6,fprintf(stdout,"Reduce: \n"));
405 if( OP_ISANYV(ops[1]) )
407 data->list=eqntoken_next(data->list);
408 r+=token_value(tokens[1],&b);
410 else { /* error, need literal */ assert(0); }
413 D(6,fprintf(stdout,"Doing Binary Op, %g op %g \n",a,b));
414 switch(token->op)
416 case OPadd: a=a+b; break;
417 case OPsub: a=a-b; break;
418 case OPmul: a=a*b; break;
419 case OPdiv:
420 if(b!=0.0) a=a/b;
421 else a=0.0;
422 break;
424 case OPexp: a=pow(a,b); break;
425 case OPnorm: if(b!=0) a=pow(pow(a,b),1/b); else a=0; break;
427 default: assert(0);
429 PUSH(a);
432 else
433 if(types[0]==4) /* open/close */
435 if(token->op==OPopen)
437 D(6,fprintf(stdout,"open { : \n"));
438 loop=1;
440 else if(token->op==OPclose)
442 D(6,fprintf(stdout,"close } : returning \n"));
443 return 0;
447 else if(token->op == OPeolist);
449 if(toend) { D(6,fprintf(stdout,"Ending: \n")); return 0; }
450 else
452 D(6,fprintf(stdout,"Stack: entry=%i now=%i %s \n",stack_entry,data->stackp,
453 ((stack_entry<data->stackp)||loop)?"looping":"returning"));
456 while((stack_entry<data->stackp)||loop);
457 return r;
460 #undef PUSH
461 #undef POP
464 static eqntoken_t *eqntok_parse(char *str)
466 int i,j,qlit,qval,qvalp;
467 int last_was_op;
468 int type[]=OPtype;
469 eqntoken_t *list,*lp,*lpn;
470 eqntop_t data;
471 eqntoken_t *cp,*np;
472 memset(&data,0,sizeof(data));
473 data.n_unget=0;
474 data.nodep=0;
475 data.litc=0;
476 data.opc=0;
477 data.str=str;
479 D(2,fprintf(stdout,"Parsing [%s]\n",str));
480 data.errors=0;
481 if(data.str[0]==EQN_DELIM) data.str++; /* eat equation prefix delimiter */
482 while(eqntoken_parse(&data));
484 for(i=0,qval=0,qlit=0,qvalp=0,j=0;i<data.nodep;i++)
486 switch(data.nodes[i].op)
488 case OPlit:
489 case OPlitm:
490 qlit++; break;
491 case OPval:
492 case OPvalm:
493 case OPeol_const:
494 qval++; break;
495 case OPeol_valp:
496 qvalp++; break;
497 default: break;
499 j++;
502 /* fixme. undo bitrot */
503 #ifdef BITROT
504 if((qval==1) && (qlit==0) && (qvalp==0) ) /* optimization for single constant, don't store eolist */
506 eq=eqn_const(data.nodes[0].z.val.x);
507 lp=list=eq.eqn;
509 else if((qvalp==1) && (qlit==0) && (qval==0) )
511 eq=eqn_valp(data.nodes[0].z.pval.xp);
512 lp=list=eq.eqn;
514 else if((qvalp==0) && (qlit==0) && (qval==0) )
516 /* it appears we have a parse error, no equation is possible */
517 lp=list=NULL;
518 fprintf(stdout,"Invalid Equation \n");
519 assert(0);
521 else
522 #endif
524 lp=list=eqn_malloc((sizeof(itemop_t)*(j+1)) + (sizeof(eqntokenlit_t)*(qlit)) + (sizeof(eqntokenval_t)*(qval)));
526 for(last_was_op=1,i=j=0;i<data.nodep;i++)
528 cp=data.nodes+i;
529 if((i+1)<data.nodep)
531 np=data.nodes+i+1;
532 if( (last_was_op) && (cp->op==OPsub) && (type[np->op]==1) )
534 np->op++; /* go from OPxxx to OPxxxm */
535 last_was_op=0;
536 continue;
539 if(type[cp->op]>1)
540 last_was_op=1;
541 else last_was_op=0;
543 lp->op=cp->op;
544 lpn=eqntoken_next(lp);
545 memcpy(lp,cp,((char *)lpn)-((char *)lp));
546 lp=lpn;
548 lp->op=OPeolist;
551 D(7,
552 do {
553 char *syms[]=OPsyms;
554 fprintf(stdout, "eqnparse: %i tokens %p ", data.nodep,lp);
555 for(lp=list;!OP_ISEND(lp->op);lp=eqntoken_next(lp))
556 fprintf(stdout, "%s ",syms[lp->op]);
557 fprintf(stdout, "%s ",syms[lp->op]);
558 fprintf(stdout,"\n");
559 } while (0)
562 return list;
569 /* only evaluate after eqntok_depend, and
570 call eqntok_eval IN PROPER ORDER
571 returns 0 if eval ok, 1 if broken dependency, 2 on parse error
574 static int eqntok_eval(float *res, eqntoken_t *list)
576 eqntoken_t *lp;
577 eqneval_t data;
578 memset(&data,0,sizeof(data));
579 data.stack[(data.stackp=0)]=0.0;
580 data.stackp++;
581 data.list=list;
583 eqn_t v;
584 v.eqn=list;
585 D(3,debug_eqn(stdout,"Evaluating \n",&v));
588 /* check that all dependencies have been taken care of */
589 for(lp=list;!OP_ISEND(lp->op);lp=eqntoken_next(lp))
591 if((lp->op==OPlit)||(lp->op==OPlitm))return 1;
594 eqntoken_eval(&data);
596 if(res!=NULL)
597 *res=data.stack[1];
599 if(data.stackp==2)return 0;
600 else return 2;
604 #ifdef UNUSED
605 static eqntoken_t *eqntok_copy(eqntoken_t *p)
607 int q;
608 eqntoken_t *np;
609 np=eqntoken_last(p);
610 q=((char *)np)-((char *)p);
611 np=eqn_malloc(q);
612 assert(np!=NULL);
613 memcpy(np,p,q);
614 return np;
618 static eqntoken_t *eqntok_copy_m(eqntoken_t *p, float m)
620 int q;
621 eqntoken_t *np,*nnp;
623 if(p->op==OPeol_const) /* common case */
625 np=eqn_malloc(sizeof(itemop_t)+sizeof(eqntokenval_t));
626 np->op=OPeol_const;
627 np->z.val.x=m*p->z.val.x;
629 else /* build np = m*(old_equation)EOL */
631 np=eqntoken_last(p);
632 q=((char *)np)-((char *)p);
634 np=eqn_malloc(q+(sizeof(itemop_t)*4)+sizeof(eqntokenval_t));
635 nnp=np;
636 nnp->op=OPval;
637 nnp->z.val.x=m;
638 nnp=eqntoken_next(nnp);
639 nnp->op=OPmul;
640 nnp=eqntoken_next(nnp);
641 nnp->op=OPopen;
642 nnp=eqntoken_next(nnp);
643 memcpy(nnp,p,q);
644 nnp=(void *) (((char *)nnp)+q-sizeof(itemop_t));
645 nnp->op=OPclose;
646 nnp=eqntoken_next(nnp);
647 nnp->op=OPeolist;
649 return np;
653 #endif
661 eqn_t equation_parse(uchar *str)
663 eqn_t eqn;
664 eqn.val=0.0;
665 eqn.eqn=eqntok_parse(str);
666 return eqn;
673 /* returns 0 if equation can be evaluated */
674 int equation_depend(eqn_t eqn, float *(lookup)(void *user, char *str), void *user)
676 int q=0;
677 float *vp;
678 eqntoken_t *lp,*list;
679 list=eqn.eqn;
680 if(list==NULL)return -1;
682 if(list->op==OPeol_const)return 0; /*optimization */
684 for(lp=list;!OP_ISEND(lp->op);lp=eqntoken_next(lp))
686 if((lp->op==OPlit)||(lp->op==OPlitm))
688 D(50,fprintf(stdout,"eqntok_depend: Looking up %s\n",lp->z.lit.lit));
689 if((vp=lookup(user,lp->z.lit.lit))!=NULL)
691 eqn_free(lp->z.lit.lit);
692 lp->z.pval.xp=vp;
693 lp->op=(lp->op==OPlit)?OPlitv:OPlitvm;
695 else q++;
698 return q;
702 int equation_eval(eqn_t *peqn)
704 int r=-1;
705 if(peqn->eqn!=NULL)
707 r=eqntok_eval(&(peqn->val), peqn->eqn);
708 D(1,fprintf(stdout,"Result: %g\n",peqn->val));
710 return r;
713 eqn_t equation_empty(eqn_t eqn)
715 eqntoken_t *eqt;
716 /* should traverse equations first to free any literal strings */
718 for(eqt=eqn.eqn;(eqt!=NULL);eqt=eqntoken_next(eqt))
720 switch(eqt->op)
722 case OPlit:
723 case OPlitm:
724 eqn_free(eqt->z.lit.lit);
725 break;
726 default: break;
728 if(OP_ISEND(eqt->op)) break;
730 eqn_free(eqn.eqn);
731 eqn.eqn=NULL;
732 eqn.val=0.0;
733 return eqn;
737 void equation_debug(eqn_t eqn, void *fp)
739 char buf[64];
740 sprintf(buf,"eqn=%g \n ",eqn.val);
741 debug_eqn(fp, buf, &eqn);