Subversion Repositories wimsdev

Rev

Rev 7847 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. /*    Copyright (C) 1998-2003 XIAO, Gang of Universite de Nice - Sophia Antipolis
  2.  *
  3.  *  This program is free software; you can redistribute it and/or modify
  4.  *  it under the terms of the GNU General Public License as published by
  5.  *  the Free Software Foundation; either version 2 of the License, or
  6.  *  (at your option) any later version.
  7.  *
  8.  *  This program is distributed in the hope that it will be useful,
  9.  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10.  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11.  *  GNU General Public License for more details.
  12.  *
  13.  *  You should have received a copy of the GNU General Public License
  14.  *  along with this program; if not, write to the Free Software
  15.  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  16.  */
  17.  
  18. /* log(-1) does not make sense in real */
  19. #ifndef NAN
  20. #define NAN log(-1)
  21. #endif
  22.  
  23.         /* Only two decimal points, less than 1 million.
  24.          * No check of buffer length which should be at least 12.
  25.          * returns the end of buffer. */
  26. char *moneyprint(char *p, double s)
  27. {
  28.     char *p1, *p2, buf[16];
  29.     int t, t1, t2;
  30.     if(s<0) {*p++='-'; s=-s;}
  31.     if(s>999999) s=999999;
  32.     t=floor(s*100+0.5); if(t>99999999) t=99999999; if(t<0) t=0;
  33.     if(t==0) {*p++='0'; *p=0; return p;}
  34.     t1=t/100; t2=t%100; p1=buf+10;
  35.     for(*p1--=t1%10+'0',t1/=10;t1>0;*p1--=t1%10+'0',t1/=10);
  36.     p2=buf+11;
  37.     if(t2) {
  38.         *p2++='.';
  39.         *p2++=t2/10+'0'; t2%=10;
  40.         if(t2) *p2++=t2+'0';
  41.     }
  42.     p1++; *p2=0; memmove(p,p1,p2-p1+1); p+=p2-p1;
  43.     return p;
  44. }
  45.  
  46. /* #define RAND_BUF_SIZE 4096
  47. static char rand_buf[RAND_BUF_SIZE];
  48. */      /* The trouble here is that httpd does not initialize
  49.          * the variable RANDOM.
  50.          * So I use time (microseconds) to get a quick solution. */
  51. void init_random(void)
  52. {
  53.     int r;
  54.     struct timeval t;
  55. /*    initstate(1,rand_buf,RAND_BUF_SIZE); */
  56.     gettimeofday(&t,NULL);
  57.     r=t.tv_usec+t.tv_sec*1000;
  58.     if(r<0) r=-r; if(r==0) r=1;
  59.     srandom(r);
  60. }
  61.  
  62.         /* gives a double random number between 0 and m */
  63. double drand(double m)
  64. {
  65.     double r;
  66.     r=((double) random()+(double) random()/(double) RAND_MAX);
  67.     return (r/(double) RAND_MAX)*m;
  68. }
  69.  
  70.         /* gives a random integer between 0 and n.
  71.          * n maybe floating, but will be rounded */
  72. double irand(double n)
  73. {
  74.     int  end,r;
  75.     if(n==0) return 0;
  76.     if(n>0) end=n; else end=-n;
  77.     r=(double) random()*end/RAND_MAX;
  78.     if(r==n) r--;
  79.     if(n>0) return r; else return -r;    
  80. }
  81.  
  82.         /* sign of d */
  83. double sign(double d)
  84. {
  85.     if(d==0) return 0;
  86.     if(d<0) return -1;
  87.     else return 1;
  88. }
  89.  
  90.         /* rounding to integer: problem with half-way rounding */
  91. double myround(double d)
  92. {
  93.     long int t;
  94.     if(d<0) t=d-0.5; else t=d+0.5;
  95.     return t;
  96. }
  97.  
  98.         /* log of base 2 */
  99. double mylog2(double d)
  100. {
  101.     return log(d)/log(2);
  102. }
  103.  
  104.         /* sec */
  105. double sec(double d)
  106. {       return 1/cos(d);        }
  107.  
  108.         /* csc */
  109. double csc(double d)
  110. {       return 1/sin(d);        }
  111.  
  112.         /* cotangent function */
  113. double cotan(double d)
  114. {
  115.     return 1/tan(d);
  116. }
  117.  
  118.         /* hyperbolic cotangent */
  119. double cotanh(double d)
  120. {
  121.     return 1/tanh(d);
  122. }
  123.  
  124.         /* factorial of an integer */
  125. double factorial(double d)
  126. {
  127.     int i,n; double t;
  128.     n=d;
  129.     if(n<0 || n!=d) return NAN;
  130.     if(n>1000) return HUGE_VAL;
  131.     t=1; for(i=1;i<=n;i++) t=t*i;
  132.     return t;
  133. }
  134.  
  135.         /* binomial coefficient */
  136. double binomial(double d1,double d2)
  137. {
  138.     return factorial(d1)/(factorial(d2)*factorial(d1-d2));
  139. }
  140.  
  141.         /* max and min */
  142. double max(double d1, double d2)
  143. {
  144.     if(!finite(d1) || !finite(d2)) return NAN;
  145.     if(d1<d2) return d2; else return d1;
  146. }
  147. double min(double d1, double d2)
  148. {
  149.     if(!finite(d1) || !finite(d2)) return NAN;
  150.     if(d1<d2) return d1; else return d2;
  151. }
  152.  
  153.         /* gcd and lcm, not really checking errors. */
  154. double gcd(double n1, double n2)
  155. {
  156.     unsigned long long int l1, l2, ll;
  157.     n1=abs(n1); n2=abs(n2);
  158.     if(!finite(n1) || !finite(n2) || n1<0 || n2<0 ||
  159.        n1>1E18 || n2>1E18) return NAN;
  160.     l1=n1; l2=n2;
  161.     if(l1<l2) {
  162.         ll=l1;l1=l2;l2=ll;
  163.     }
  164.     if(l1==0) return NAN;
  165.     while(l2>0) {
  166.         ll=l2;l2=l1%l2;l1=ll;
  167.     }
  168.     return l1;
  169. }
  170.  
  171. double lcm(double n1, double n2)
  172. {
  173.     return n1*n2/gcd(n1,n2);
  174. }
  175.  
  176. struct {
  177.     char *name;
  178.     int type;
  179.     double val;
  180.     double (*f1) (double parm);
  181.     double (*f2) (double parm1, double parm2);
  182. } evalname[]={
  183.       {"Argch", 1,      0,      acosh,  NULL},
  184.       {"Argsh", 1,      0,      asinh,  NULL},
  185.       {"Argth", 1,      0,      atanh,  NULL},
  186.       {"E",     0,      M_E,    NULL,   NULL},
  187.       {"EULER", 0,      0.57721566490153286,    NULL,   NULL},
  188.       {EV_S,    0,      0,      NULL,   NULL},
  189.       {EV_T,    0,      0,      NULL,   NULL},
  190.       {EV_X,    0,      0,      NULL,   NULL},
  191.       {EV_Y,    0,      0,      NULL,   NULL},
  192.       {"Euler", 0,      0.57721566490153286,    NULL,   NULL},
  193.       {"Inf",   0,      1,      log,    NULL},
  194.       {"NaN",   0,      0,      log,    NULL},
  195.       {"PI",    0,      M_PI,   NULL,   NULL},
  196.       {"Pi",    0,      M_PI,   NULL,   NULL},
  197.       {"abs",   1,      0,      fabs,   NULL},
  198.       {"acos",  1,      0,      acos,   NULL},
  199.       {"acosh", 1,      0,      acosh,  NULL},
  200.       {"arccos",1,      0,      acos,   NULL},
  201.       {"arcsin",1,      0,      asin,   NULL},
  202.       {"arctan",1,      0,      atan,   NULL},
  203.       {"arctg", 1,      0,      atan,   NULL},
  204.       {"argch", 1,      0,      acosh,  NULL},
  205.       {"argsh", 1,      0,      asinh,  NULL},
  206.       {"argth", 1,      0,      atanh,  NULL},
  207.       {"asin",  1,      0,      asin,   NULL},
  208.       {"asinh", 1,      0,      asinh,  NULL},
  209.       {"atan",  1,      0,      atan,   NULL},
  210.       {"atanh", 1,      0,      atanh,  NULL},
  211.       {"binomial",2,    0,      NULL,   binomial},
  212.       {"ceil",  1,      0,      ceil,   NULL}, /* round-up integer */
  213.       {"ch",    1,      0,      cosh,   NULL},
  214.       {"cos",   1,      0,      cos,    NULL},
  215.       {"cosh",  1,      0,      cosh,   NULL},
  216.       {"cot",   1,      0,      cotan,  NULL},
  217.       {"cotan", 1,      0,      cotan,  NULL},
  218.       {"cotanh",1,      0,      cotanh, NULL},
  219.       {"coth",  1,      0,      cotanh, NULL},
  220.       {"csc",   1,      0,      csc,    NULL},
  221.       {"ctg",   1,      0,      cotan,  NULL},
  222.       {"cth",   1,      0,      cotanh, NULL},
  223.       {"drand", 1,      0,      drand,  NULL},
  224.       {"e",     0,      M_E,    NULL,   NULL},
  225.       {"erf",   1,      0,      erf,    NULL},
  226.       {"erfc",  1,      0,      erfc,   NULL},
  227.       {"euler", 0,      0.57721566490153286,    NULL,   NULL},
  228.       {"exp",   1,      0,      exp,    NULL},
  229.       {"factorial",1,   0,      factorial,      NULL},
  230.       {"floor", 1,      0,      floor,  NULL},
  231.       {"gcd",   2,      0,      NULL,   gcd},
  232.       {"irand", 1,      0,      irand,  NULL},
  233. /*      {"j0",  1,      0,      j0,     NULL}, */ /* Bessel functions */
  234. /*      {"j1",  1,      0,      j1,     NULL}, */
  235.       {"lcm",   2,      0,      NULL,   lcm},
  236.       {"lg",    1,      0,      log10,  NULL},
  237.       {"lgamma",1,      0,      lgamma, NULL}, /* log of Gamma function */
  238.       {"ln",    1,      0,      log,    NULL},
  239.       {"log",   1,      0,      log,    NULL},
  240.       {"log10", 1,      0,      log10,  NULL},
  241.       {"log2",  1,      0,      mylog2, NULL},
  242.       {"max",   2,      0,      NULL,   max},
  243.       {"min",   2,      0,      NULL,   min},
  244.       {"pi",    0,      M_PI,   NULL,   NULL},
  245.       {"pow",   2,      0,      NULL,   pow},
  246.       {"rand",  1,      0,      drand,  NULL},
  247.       {"randdouble",1,  0,      drand,  NULL},
  248.       {"randfloat",1,   0,      drand,  NULL},
  249.       {"randint",1,     0,      irand,  NULL},
  250.       {"random",1,      0,      drand,  NULL},
  251.       {"randreal",1,    0,      drand,  NULL},
  252.       {"rint",  1,      0,      myround,        NULL}, /* closest integer */
  253.       {"round", 1,      0,      myround,        NULL}, /* closest integer */
  254.       {"sec",   1,      0,      sec,    NULL},
  255.       {"sgn",   1,      0,      sign,   NULL}, /* sign of the value */
  256.       {"sh",    1,      0,      sinh,   NULL},
  257.       {"sign",  1,      0,      sign,   NULL}, /* sign of the value */
  258.       {"sin",   1,      0,       sin,   NULL},
  259.       {"sinh",  1,      0,      sinh,   NULL},
  260.       {"sqrt",  1,      0,      sqrt,   NULL},
  261.       {"tan",   1,      0,      tan,    NULL},
  262.       {"tanh",  1,      0,      tanh,   NULL},
  263.       {"tg",    1,      0,      tan,    NULL},
  264.       {"th",    1,      0,      tanh,   NULL},
  265. /*      {"y0",  1,      0,      y0,     NULL}, */
  266. /*      {"y1",  1,      0,      y1,     NULL},  */
  267. };
  268. #define evalname_no (sizeof(evalname)/sizeof(evalname[0]))
  269.  
  270. int get_evalcnt(void) {return evalname_no;}
  271. char *get_evalname(int i) {return evalname[i].name;}
  272. int get_evaltype(int i) {return evalname[i].type;}
  273. int evaltab_verify(void) {return verify_order(evalname,evalname_no,sizeof(evalname[0]));}
  274. int search_evaltab(char *p) {
  275.     return search_list(evalname,evalname_no,sizeof(evalname[0]),p);
  276. }
  277.  
  278. static char *evalue_pt;
  279. int evalue_error;
  280.  
  281. int get_evalue_error(void) { return evalue_error; }
  282. void set_evalue_error(int e) {evalue_error=e; return;}
  283.  
  284. /* prepare pointer for evaluation */
  285. void set_evalue_pointer(char *p)
  286. {
  287.     evalue_pt=p;
  288. }
  289.  
  290.         /* get position of name in nametable */
  291. int eval_getpos(char *name)
  292. {
  293.     return search_list(evalname,evalname_no,sizeof(evalname[0]),name);
  294. }
  295.  
  296.         /* set value to name */
  297. void eval_setval(int pos, double v)
  298. {
  299.     if(pos>=0 && pos<evalname_no) evalname[pos].val=v;
  300. }
  301.  
  302.         /* get string pointer (after evaluation) */
  303. char *get_evalue_pointer(void)
  304. {
  305.     return evalue_pt;
  306. }
  307.  
  308. double _evalue(int ord)
  309. {
  310.     double d,dd;
  311.     int i,k;
  312.     char buf[32];
  313.  
  314.    
  315.     if(evalue_error) return NAN;
  316.     d=0;
  317.     while(*evalue_pt=='+') evalue_pt++;
  318.     if(*evalue_pt==0) return 0; /* empty string */
  319.     switch(*evalue_pt) {
  320.         case '(':
  321.         evalue_pt++; d=_evalue(')');goto vld;
  322.         case '|':
  323.         if(ord=='|') {
  324.             evalue_pt++; return 0;
  325.         }
  326.         evalue_pt++; d=fabs(_evalue('|'));goto vld;
  327.         case '-':
  328.         evalue_pt++; d=-_evalue(6);goto vld;
  329.     }
  330.     if((128&*evalue_pt)!=0) {
  331.         k=(*evalue_pt)&255; evalue_pt++;
  332.         if(k>=130 && k<140) {
  333.             i=(k-130)*200; k=(*evalue_pt)&255; evalue_pt++;
  334.             if(k<33 || k>=233) goto badeval;
  335.             i+=k-33; if(i<0 || i>=evalname_no) goto badeval;
  336.             goto ename;
  337.         }
  338.         if(k>=140 && k<150) {
  339.             i=(k-140)*200; k=(*evalue_pt)&255; evalue_pt++;
  340.             if(k<33 || k>=233) goto badeval;
  341.             if(ev_var==NULL || ev_varcnt==NULL) goto badeval;
  342.             i+=k-33; if(i<0 || i>=*ev_varcnt) goto badeval;
  343.             goto vname;
  344.         }
  345.         evalue_pt++; goto badeval;
  346.     }
  347.     if(*evalue_pt=='.' || myisdigit(*evalue_pt))
  348.         {d=strtod(evalue_pt,&evalue_pt);goto binary;}
  349.     for(i=0;myisalnum(*(evalue_pt+i)) && i<16; i++)
  350.       buf[i]=*(evalue_pt+i);
  351.     buf[i]=0; evalue_pt+=i;
  352.     if(i==0) goto badeval;
  353.     if(ev_varcnt!=NULL && ev_var!=NULL && *ev_varcnt>0)
  354.       for(i=0;i<*ev_varcnt;i++) {
  355.           if(strcmp(buf,ev_var[i].name)==0) {
  356.               vname: d=ev_var[i].value; goto vld;
  357.           }
  358.       }
  359.     i=search_list(evalname,evalname_no,sizeof(evalname[0]),buf);
  360.     ename: if(i>=0) switch(evalname[i].type) {
  361.         case 0: {
  362.             d=evalname[i].val;
  363.             if(evalname[i].f1!=NULL) {
  364.                 if(d==0) d=NAN; if(d==1) d=HUGE_VAL;
  365.             }
  366.             break;
  367.         }
  368.         case 1: {
  369.             if(*evalue_pt!='(') return NAN;
  370.             evalue_pt++;
  371.             d=evalname[i].f1(_evalue(')')); break;
  372.         }
  373.         case 2: {
  374.             double parm1,parm2;
  375.             if(*evalue_pt!='(') return NAN;
  376.             evalue_pt++;
  377.             parm1=_evalue(',');parm2=_evalue(')');
  378.             d=evalname[i].f2(parm1,parm2); break;
  379.         }
  380.         default: {      /* This is impossible. */
  381.             return NAN;
  382.         }
  383.     }
  384.     else {
  385.         badeval: evalue_error=-1; return NAN;
  386.     }
  387.   vld:
  388.     if(evalue_error) return NAN;
  389.   binary:
  390.     if(*evalue_pt=='!') {
  391.         evalue_pt++; d=factorial(d);
  392.     }
  393.     if(*evalue_pt==ord) {evalue_pt++;goto ok;}
  394.     if(*evalue_pt==0 ||
  395.        (ord<10 && (*evalue_pt==',' || *evalue_pt==';' || *evalue_pt==')'
  396.                    || *evalue_pt=='|')))
  397.        goto ok;
  398.     switch(*evalue_pt) {
  399.                 case '+':
  400.                         if(ord<=8) break;
  401.                         evalue_pt++; d+=_evalue(8);goto vld;
  402.                 case '-':
  403.                         if(ord<=8) break;
  404.                         evalue_pt++; d-=_evalue(8);goto vld;
  405.                 case '*':
  406.                         if(ord<=6) break;
  407.                         evalue_pt++; d*=_evalue(6);goto vld;
  408.                 case '/':
  409.                         if(ord<=6) break;
  410.                         evalue_pt++; dd=_evalue(6);
  411.                         if(dd==0) {evalue_error=10;return NAN;}
  412.                         d/=dd;goto vld;
  413.         case '%': {
  414.             int di, ddi;
  415.             if(ord<=6) break;
  416.             evalue_pt++; dd=_evalue(6);
  417.             if(dd==0) {evalue_error=10;return NAN;}
  418.             di=d; ddi=dd; d=di%ddi;goto vld;
  419.         }
  420.         case '^': {
  421.             if(ord<5) break;
  422.             evalue_pt++; d=pow(d,_evalue(5));goto vld;
  423.         }
  424.         default : {
  425.             /*if(ord<=6) break;
  426.             d*=_evalue(6);goto vld;*/
  427.             return NAN;
  428.         }
  429.     }
  430.     ok: return d;
  431. }
  432.  
  433.         /* substitute variable names by their environment strings
  434.          * The buffer pointed to by p must have enough space
  435.          * (defined by MAX_LINELEN). */
  436. char *_substit(char *p)
  437. {
  438.     return p;
  439. }
  440.  
  441. char *(*substitute) (char *p)=_substit;
  442.  
  443.         /* evalue a string to double */
  444. double strevalue(char *p)
  445. {
  446.     char buf[MAX_LINELEN+1];
  447.    
  448.     if(p==NULL) return 0;
  449.     mystrncpy(buf,p,sizeof(buf));
  450.     substitute(buf); nospace(buf);
  451.     if(check_parentheses(buf,0)) return NAN;
  452.     set_evalue_error(0);
  453.     set_evalue_pointer(buf);
  454.     return _evalue(10);
  455. }
  456.  
  457.         /* compile an expression for faster evaluation
  458.          * returns -1 if cannot be compiled.
  459.          * else returns the number of compilations. */
  460. int evalue_compile(char *p)
  461. {
  462.     char *p1, *p2, *pe, name[256], buf[8];
  463.     int i,k;
  464.    
  465.     k=0;
  466.     for(p1=p; *p1; p1++) if((128&*p1)!=0) return -1;
  467.     nospace(p);
  468.     for(p1=find_mathvar_start(p); *p1; p1=find_mathvar_start(pe)) {
  469.         pe=find_mathvar_end(p1);
  470.         if(!myisalpha(*p1)) continue;
  471.         p2=pe; if(p2-p1>16) continue;
  472.         memmove(name,p1,p2-p1); name[p2-p1]=0;
  473.         if(ev_varcnt!=NULL && ev_var!=NULL && *ev_varcnt>0) {
  474.             for(i=0;i<*ev_varcnt && strcmp(name,ev_var[i].name)!=0;i++);
  475.             if(i<*ev_varcnt && i<2000) {
  476.                 buf[0]=i/200+140; buf[1]=i%200+33; buf[2]=0;
  477.                 string_modify(p,p1,p2,"%s",buf);
  478.                 pe=p1+2; k++; continue;
  479.             }
  480.         }
  481.         i=search_list(evalname,evalname_no,sizeof(evalname[0]),name);
  482.         if(i>=0 && i<2000) {
  483.             buf[0]=i/200+130; buf[1]=i%200+33; buf[2]=0;
  484.             string_modify(p,p1,p2,"%s",buf);
  485.             pe=p1+2; k++; continue;
  486.         }
  487.     }
  488.     return k;
  489. }
  490.  
  491.