Back to index

wims  3.65+svn20090927
evalue.c
Go to the documentation of this file.
00001 /*    Copyright (C) 1998-2003 XIAO, Gang of Universite de Nice - Sophia Antipolis
00002  *
00003  *  This program is free software; you can redistribute it and/or modify
00004  *  it under the terms of the GNU General Public License as published by
00005  *  the Free Software Foundation; either version 2 of the License, or
00006  *  (at your option) any later version.
00007  *
00008  *  This program is distributed in the hope that it will be useful,
00009  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00010  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011  *  GNU General Public License for more details.
00012  *
00013  *  You should have received a copy of the GNU General Public License
00014  *  along with this program; if not, write to the Free Software
00015  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00016  */
00017 
00018 /* log(-1) does not make sense in real */
00019 #ifndef NAN
00020 #define NAN log(-1)
00021 #endif
00022 
00023        /* Only two decimal points, less than 1 million.
00024         * No check of buffer length which should be at least 12.
00025         * returns the end of buffer. */
00026 char *moneyprint(char *p, double s)
00027 {
00028     char *p1, *p2, buf[16];
00029     int t, t1, t2;
00030     if(s<0) {*p++='-'; s=-s;}
00031     if(s>999999) s=999999;
00032     t=floor(s*100+0.5); if(t>99999999) t=99999999; if(t<0) t=0;
00033     if(t==0) {*p++='0'; *p=0; return p;}
00034     t1=t/100; t2=t%100; p1=buf+10;
00035     for(*p1--=t1%10+'0',t1/=10;t1>0;*p1--=t1%10+'0',t1/=10);
00036     p2=buf+11;
00037     if(t2) {
00038        *p2++='.';
00039        *p2++=t2/10+'0'; t2%=10;
00040        if(t2) *p2++=t2+'0';
00041     }
00042     p1++; *p2=0; memmove(p,p1,p2-p1+1); p+=p2-p1;
00043     return p;
00044 }
00045 
00046 /* #define RAND_BUF_SIZE 4096
00047 static char rand_buf[RAND_BUF_SIZE];
00048 */     /* The trouble here is that httpd does not initialize 
00049         * the variable RANDOM. 
00050         * So I use time (microseconds) to get a quick solution. */
00051 void init_random(void)
00052 {
00053     int r;
00054     struct timeval t;
00055 /*    initstate(1,rand_buf,RAND_BUF_SIZE); */
00056     gettimeofday(&t,NULL);
00057     r=t.tv_usec+t.tv_sec*1000;
00058     if(r<0) r=-r; if(r==0) r=1;
00059     srandom(r);
00060 }
00061 
00062        /* gives a double random number between 0 and m */
00063 double drand(double m)
00064 {
00065     double r;
00066     r=((double) random()+(double) random()/(double) RAND_MAX);
00067     return (r/(double) RAND_MAX)*m;
00068 }
00069 
00070        /* gives a random integer between 0 and n.
00071         * n maybe floating, but will be rounded */
00072 double irand(double n)
00073 {
00074     int  end,r;
00075     if(n==0) return 0;
00076     if(n>0) end=n; else end=-n;
00077     r=(double) random()*end/RAND_MAX;
00078     if(r==n) r--;
00079     if(n>0) return r; else return -r;    
00080 }
00081 
00082        /* sign of d */
00083 double sign(double d)
00084 {
00085     if(d==0) return 0;
00086     if(d<0) return -1;
00087     else return 1;
00088 }
00089 
00090        /* rounding to integer: problem with half-way rounding */
00091 double myround(double d)
00092 {
00093     long int t;
00094     if(d<0) t=d-0.5; else t=d+0.5;
00095     return t;
00096 }
00097 
00098        /* log of base 2 */
00099 double mylog2(double d)
00100 {
00101     return log(d)/log(2);
00102 }
00103 
00104        /* sec */
00105 double sec(double d)
00106 {      return 1/cos(d);     }
00107 
00108        /* csc */
00109 double csc(double d)
00110 {      return 1/sin(d);     }
00111 
00112        /* cotangent function */
00113 double cotan(double d)
00114 {
00115     return 1/tan(d);
00116 }
00117 
00118        /* hyperbolic cotangent */
00119 double cotanh(double d)
00120 {
00121     return 1/tanh(d);
00122 }
00123 
00124        /* factorial of an integer */
00125 double factorial(double d)
00126 {
00127     int i,n; double t;
00128     n=d;
00129     if(n<0 || n!=d) return NAN;
00130     if(n>1000) return HUGE_VAL;
00131     t=1; for(i=1;i<=n;i++) t=t*i;
00132     return t;
00133 }
00134 
00135        /* binomial coefficient */
00136 double binomial(double d1,double d2)
00137 {
00138     return factorial(d1)/(factorial(d2)*factorial(d1-d2));
00139 }
00140 
00141        /* max and min */
00142 double max(double d1, double d2)
00143 {
00144     if(!finite(d1) || !finite(d2)) return NAN;
00145     if(d1<d2) return d2; else return d1;
00146 }
00147 double min(double d1, double d2)
00148 {
00149     if(!finite(d1) || !finite(d2)) return NAN;
00150     if(d1<d2) return d1; else return d2;
00151 }
00152 
00153        /* gcd and lcm, not really checking errors. */
00154 double gcd(double n1, double n2)
00155 {
00156     unsigned long long int l1, l2, ll;
00157     n1=abs(n1); n2=abs(n2);
00158     if(!finite(n1) || !finite(n2) || n1<0 || n2<0 ||
00159        n1>1E18 || n2>1E18) return NAN;
00160     l1=n1; l2=n2;
00161     if(l1<l2) {
00162        ll=l1;l1=l2;l2=ll;
00163     }
00164     if(l1==0) return NAN;
00165     while(l2>0) {
00166        ll=l2;l2=l1%l2;l1=ll;
00167     }
00168     return l1;
00169 }
00170 
00171 double lcm(double n1, double n2)
00172 {
00173     return n1*n2/gcd(n1,n2);
00174 }
00175 
00176 struct {
00177     char *name;
00178     int type;
00179     double val;
00180     double (*f1) (double parm);
00181     double (*f2) (double parm1, double parm2);
00182 } evalname[]={
00183       {"Argch",      1,     0,     acosh, NULL},
00184       {"Argsh",      1,     0,     asinh, NULL},
00185       {"Argth",      1,     0,     atanh, NULL},
00186       {"E",   0,     M_E,   NULL,  NULL},
00187       {"EULER",      0,     0.57721566490153286, NULL,  NULL},
00188       {EV_S,  0,     0,     NULL,  NULL},
00189       {EV_T,  0,     0,     NULL,  NULL},
00190       {EV_X,  0,     0,     NULL,  NULL},
00191       {EV_Y,  0,     0,     NULL,  NULL},
00192       {"Euler",      0,     0.57721566490153286, NULL,  NULL},
00193       {"Inf", 0,     1,     log,   NULL},
00194       {"NaN", 0,     0,     log,   NULL},
00195       {"PI",  0,     M_PI,  NULL,  NULL},
00196       {"Pi",  0,     M_PI,  NULL,  NULL},
00197       {"abs", 1,     0,     fabs,  NULL},
00198       {"acos",       1,     0,     acos,  NULL},
00199       {"acosh",      1,     0,     acosh, NULL},
00200       {"arccos",1,   0,     acos,  NULL},
00201       {"arcsin",1,   0,     asin,  NULL},
00202       {"arctan",1,   0,     atan,  NULL},
00203       {"arctg",      1,     0,     atan,  NULL},
00204       {"argch",      1,     0,     acosh, NULL},
00205       {"argsh",      1,     0,     asinh, NULL},
00206       {"argth",      1,     0,     atanh, NULL},
00207       {"asin",       1,     0,     asin,  NULL},
00208       {"asinh",      1,     0,     asinh, NULL},
00209       {"atan",       1,     0,     atan,  NULL},
00210       {"atanh",      1,     0,     atanh, NULL},
00211       {"binomial",2, 0,     NULL,  binomial},
00212       {"ceil",       1,     0,     ceil,  NULL}, /* round-up integer */
00213       {"ch",  1,     0,     cosh,  NULL},
00214       {"cos", 1,     0,     cos,   NULL},
00215       {"cosh",       1,     0,     cosh,  NULL},
00216       {"cot", 1,     0,     cotan, NULL},
00217       {"cotan",      1,     0,     cotan, NULL},
00218       {"cotanh",1,   0,     cotanh,       NULL},
00219       {"coth",       1,     0,     cotanh,       NULL},
00220       {"csc", 1,     0,     csc,   NULL},
00221       {"ctg", 1,     0,     cotan, NULL},
00222       {"cth", 1,     0,     cotanh,       NULL},
00223       {"drand",      1,     0,     drand, NULL},
00224       {"e",   0,     M_E,   NULL,  NULL},
00225       {"erf", 1,     0,     erf,   NULL},
00226       {"erfc",       1,     0,     erfc,  NULL},
00227       {"euler",      0,     0.57721566490153286, NULL,  NULL},
00228       {"exp", 1,     0,     exp,   NULL},
00229       {"factorial",1,       0,     factorial,    NULL},
00230       {"floor",      1,     0,     floor, NULL},
00231       {"gcd", 2,     0,     NULL,  gcd},
00232       {"irand",      1,     0,     irand, NULL},
00233 /*      {"j0",       1,     0,     j0,    NULL}, */ /* Bessel functions */
00234 /*      {"j1",       1,     0,     j1,    NULL}, */
00235       {"lcm", 2,     0,     NULL,  lcm},
00236       {"lg",  1,     0,     log10, NULL},
00237       {"lgamma",1,   0,     lgamma,       NULL}, /* log of Gamma function */
00238       {"ln",  1,     0,     log,   NULL},
00239       {"log", 1,     0,     log,   NULL},
00240       {"log10",      1,     0,     log10, NULL},
00241       {"log2",       1,     0,     mylog2,       NULL},
00242       {"max", 2,     0,     NULL,  max},
00243       {"min", 2,     0,     NULL,  min},
00244       {"pi",  0,     M_PI,  NULL,  NULL},
00245       {"pow", 2,     0,     NULL,  pow},
00246       {"rand",       1,     0,     drand, NULL},
00247       {"randdouble",1,      0,     drand, NULL},
00248       {"randfloat",1,       0,     drand, NULL},
00249       {"randint",1,  0,     irand, NULL},
00250       {"random",1,   0,     drand, NULL},
00251       {"randreal",1, 0,     drand, NULL},
00252       {"rint",       1,     0,     myround,      NULL}, /* closest integer */
00253       {"round",      1,     0,     myround,      NULL}, /* closest integer */
00254       {"sec", 1,     0,     sec,   NULL},
00255       {"sgn", 1,     0,     sign,  NULL}, /* sign of the value */
00256       {"sh",  1,     0,     sinh,  NULL},
00257       {"sign",       1,     0,     sign,  NULL}, /* sign of the value */
00258       {"sin", 1,     0,      sin,  NULL},
00259       {"sinh",       1,     0,     sinh,  NULL},
00260       {"sqrt",       1,     0,     sqrt,  NULL},
00261       {"tan", 1,     0,     tan,   NULL},
00262       {"tanh",       1,     0,     tanh,  NULL},
00263       {"tg",  1,     0,     tan,   NULL},
00264       {"th",  1,     0,     tanh,  NULL},
00265 /*      {"y0",       1,     0,     y0,    NULL}, */
00266 /*      {"y1",       1,     0,     y1,    NULL},  */
00267 };
00268 #define evalname_no (sizeof(evalname)/sizeof(evalname[0]))
00269 
00270 int get_evalcnt(void) {return evalname_no;}
00271 char *get_evalname(int i) {return evalname[i].name;}
00272 int get_evaltype(int i) {return evalname[i].type;}
00273 int evaltab_verify(void) {return verify_order(evalname,evalname_no,sizeof(evalname[0]));}
00274 int search_evaltab(char *p) {
00275     return search_list(evalname,evalname_no,sizeof(evalname[0]),p);
00276 }
00277 
00278 static char *evalue_pt;
00279 int evalue_error;
00280 
00281 int get_evalue_error(void) { return evalue_error; }
00282 void set_evalue_error(int e) {evalue_error=e; return;}
00283 
00284 /* prepare pointer for evaluation */
00285 void set_evalue_pointer(char *p)
00286 {
00287     evalue_pt=p;
00288 }
00289 
00290        /* get position of name in nametable */
00291 int eval_getpos(char *name)
00292 {
00293     return search_list(evalname,evalname_no,sizeof(evalname[0]),name);
00294 }
00295 
00296        /* set value to name */
00297 void eval_setval(int pos, double v)
00298 {
00299     if(pos>=0 && pos<evalname_no) evalname[pos].val=v;
00300 }
00301 
00302        /* get string pointer (after evaluation) */
00303 char *get_evalue_pointer(void)
00304 {
00305     return evalue_pt;
00306 }
00307 
00308 double _evalue(int ord)
00309 {
00310     double d,dd;
00311     int i,k;
00312     char buf[32];
00313 
00314     
00315     if(evalue_error) return NAN;
00316     d=0;
00317     while(*evalue_pt=='+') evalue_pt++;
00318     if(*evalue_pt==0) return 0; /* empty string */
00319     switch(*evalue_pt) {
00320        case '(':
00321        evalue_pt++; d=_evalue(')');goto vld;
00322        case '|':
00323        if(ord=='|') {
00324            evalue_pt++; return 0;
00325        }
00326        evalue_pt++; d=fabs(_evalue('|'));goto vld;
00327        case '-':
00328        evalue_pt++; d=-_evalue(6);goto vld;
00329     }
00330     if((128&*evalue_pt)!=0) {
00331        k=(*evalue_pt)&255; evalue_pt++;
00332        if(k>=130 && k<140) {
00333            i=(k-130)*200; k=(*evalue_pt)&255; evalue_pt++;
00334            if(k<33 || k>=233) goto badeval;
00335            i+=k-33; if(i<0 || i>=evalname_no) goto badeval;
00336            goto ename;
00337        }
00338        if(k>=140 && k<150) {
00339            i=(k-140)*200; k=(*evalue_pt)&255; evalue_pt++;
00340            if(k<33 || k>=233) goto badeval;
00341            if(ev_var==NULL || ev_varcnt==NULL) goto badeval;
00342            i+=k-33; if(i<0 || i>=*ev_varcnt) goto badeval;
00343            goto vname;
00344        }
00345        evalue_pt++; goto badeval;
00346     }
00347     if(*evalue_pt=='.' || myisdigit(*evalue_pt))
00348        {d=strtod(evalue_pt,&evalue_pt);goto binary;}
00349     for(i=0;myisalnum(*(evalue_pt+i)) && i<16; i++)
00350       buf[i]=*(evalue_pt+i);
00351     buf[i]=0; evalue_pt+=i;
00352     if(i==0) goto badeval;
00353     if(ev_varcnt!=NULL && ev_var!=NULL && *ev_varcnt>0)
00354       for(i=0;i<*ev_varcnt;i++) {
00355          if(strcmp(buf,ev_var[i].name)==0) {
00356              vname: d=ev_var[i].value; goto vld;
00357          }
00358       }
00359     i=search_list(evalname,evalname_no,sizeof(evalname[0]),buf);
00360     ename: if(i>=0) switch(evalname[i].type) {
00361        case 0: {
00362            d=evalname[i].val;
00363            if(evalname[i].f1!=NULL) {
00364               if(d==0) d=NAN; if(d==1) d=HUGE_VAL;
00365            }
00366            break;
00367        }
00368        case 1: {
00369            if(*evalue_pt!='(') return NAN;
00370            evalue_pt++;
00371            d=evalname[i].f1(_evalue(')')); break;
00372        }
00373        case 2: {
00374            double parm1,parm2;
00375            if(*evalue_pt!='(') return NAN;
00376            evalue_pt++;
00377            parm1=_evalue(',');parm2=_evalue(')');
00378            d=evalname[i].f2(parm1,parm2); break;
00379        }
00380        default: {    /* This is impossible. */
00381            return NAN;
00382        }
00383     }
00384     else {
00385        badeval: evalue_error=-1; return NAN;
00386     }
00387   vld:
00388     if(evalue_error) return NAN;
00389   binary:
00390     if(*evalue_pt=='!') {
00391        evalue_pt++; d=factorial(d);
00392     }
00393     if(*evalue_pt==ord) {evalue_pt++;goto ok;}
00394     if(*evalue_pt==0 || 
00395        (ord<10 && (*evalue_pt==',' || *evalue_pt==';' || *evalue_pt==')'
00396                  || *evalue_pt=='|')))
00397        goto ok;
00398     switch(*evalue_pt) {
00399               case '+':
00400                      if(ord<=8) break;
00401                      evalue_pt++; d+=_evalue(8);goto vld;
00402               case '-':
00403                      if(ord<=8) break;
00404                      evalue_pt++; d-=_evalue(8);goto vld;
00405               case '*':
00406                      if(ord<=6) break;
00407                      evalue_pt++; d*=_evalue(6);goto vld;
00408               case '/':
00409                      if(ord<=6) break;
00410                      evalue_pt++; dd=_evalue(6);
00411                      if(dd==0) {evalue_error=10;return NAN;}
00412                      d/=dd;goto vld;
00413        case '%': {
00414            int di, ddi;
00415            if(ord<=6) break;
00416            evalue_pt++; dd=_evalue(6);
00417            if(dd==0) {evalue_error=10;return NAN;}
00418            di=d; ddi=dd; d=di%ddi;goto vld;
00419        }
00420        case '^': {
00421            if(ord<5) break;
00422            evalue_pt++; d=pow(d,_evalue(5));goto vld;
00423        }
00424        default : {
00425            /*if(ord<=6) break;
00426            d*=_evalue(6);goto vld;*/
00427            return NAN;
00428        }
00429     }
00430     ok: return d;
00431 }
00432 
00433        /* substitute variable names by their environment strings
00434         * The buffer pointed to by p must have enough space
00435         * (defined by MAX_LINELEN). */
00436 char *_substit(char *p)
00437 {
00438     return p;
00439 }
00440 
00441 char *(*substitute) (char *p)=_substit;
00442 
00443        /* evalue a string to double */
00444 double strevalue(char *p)
00445 {
00446     char buf[MAX_LINELEN+1];
00447     
00448     if(p==NULL) return 0;
00449     mystrncpy(buf,p,sizeof(buf));
00450     substitute(buf); nospace(buf);
00451     if(check_parentheses(buf,0)) return NAN;
00452     set_evalue_error(0);
00453     set_evalue_pointer(buf);
00454     return _evalue(10);
00455 }
00456 
00457        /* compile an expression for faster evaluation
00458         * returns -1 if cannot be compiled.
00459         * else returns the number of compilations. */
00460 int evalue_compile(char *p)
00461 {
00462     char *p1, *p2, *pe, name[256], buf[8];
00463     int i,k;
00464     
00465     k=0;
00466     for(p1=p; *p1; p1++) if((128&*p1)!=0) return -1;
00467     nospace(p);
00468     for(p1=find_mathvar_start(p); *p1; p1=find_mathvar_start(pe)) {
00469        pe=find_mathvar_end(p1);
00470        if(!myisalpha(*p1)) continue;
00471        p2=pe; if(p2-p1>16) continue;
00472        memmove(name,p1,p2-p1); name[p2-p1]=0;
00473        if(ev_varcnt!=NULL && ev_var!=NULL && *ev_varcnt>0) {
00474            for(i=0;i<*ev_varcnt && strcmp(name,ev_var[i].name)!=0;i++);
00475            if(i<*ev_varcnt && i<2000) {
00476               buf[0]=i/200+140; buf[1]=i%200+33; buf[2]=0;
00477               string_modify(p,p1,p2,"%s",buf);
00478               pe=p1+2; k++; continue;
00479            }
00480        }
00481        i=search_list(evalname,evalname_no,sizeof(evalname[0]),name);
00482        if(i>=0 && i<2000) {
00483            buf[0]=i/200+130; buf[1]=i%200+33; buf[2]=0;
00484            string_modify(p,p1,p2,"%s",buf);
00485            pe=p1+2; k++; continue;
00486        }
00487     }
00488     return k;
00489 }
00490