Subversion Repositories wimsdev

Rev

Rev 7798 | Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
10 reyssat 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