← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to contrib/MathEx/mathex.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-02 18:12:15 UTC
  • mfrom: (1.2.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20090902181215-yla8zvcas2ucvkm9
[Christophe Prud'homme]
* New upstream release
  + fixed surface mesh orientation bug introduced in 2.4.0;
  + mesh and graphics code refactoring;
  + small usability enhancements and bug fixes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/MathEx/mathex.cpp
 
1
///////////////////////////////////////////////////////////////////////////
 
2
// mathex 0.2.3 (beta) - Copyright (C) 2000-2003, by Sadao Massago       //
 
3
// file: mathex.h (math expression evaluator header file)                //
 
4
// requires: none                                                        //
 
5
// project web page: http://sscilib.sourceforge.net/                     //
 
6
// ----------------------------------------------------------------------//
 
7
// The mathex library and related files is licenced under the term of    //
 
8
// GNU LGPL (Lesser General Public License) version 2.1 or latter        //
 
9
// with exceptions that allow for static linking.                        //
 
10
// See license.txt for detail.                                           //
 
11
// For GNU LGPL, see lesser.txt.                                         //
 
12
// For information over GNU or GNU compatible license, visit the site    //
 
13
// http://www.gnu.org.                                                   //
 
14
///////////////////////////////////////////////////////////////////////////
 
15
 
 
16
////////////////////////////////////////////////////////////////////////////
 
17
// references:
 
18
//-------------------------------------------------------------------------
 
19
// title: Algoritms and Data Structures
 
20
// author: Niklaus Wirth
 
21
// publisher: Prentice-Hall
 
22
// year: 1989
 
23
//-------------------------------------------------------------------------
 
24
// title: The C++ Programming Language
 
25
// edition: thrird
 
26
// author: Bjarne Stroustrup
 
27
// publisher: Addilson-Wesley
 
28
// year: 1997
 
29
//-------------------------------------------------------------------------
 
30
// title: building your own C interpreter
 
31
// author: Schildt, Helbert
 
32
// journal: Dr. Dobb's Jornal (http://www.ddj.com)
 
33
// number: August/1989
 
34
// pages: 38-49, 110-122
 
35
///////////////////////////////////////////////////////////////////////////
 
36
 
 
37
// #define _DEBUG_
 
38
 
 
39
#include "mathex.h"
 
40
 
 
41
#include <string.h>
 
42
#include <stdlib.h>
 
43
#include <math.h>
 
44
#include <string>
 
45
#include <vector>
 
46
#include <cctype>
 
47
// debug purpose
 
48
#ifdef _DEBUG_
 
49
  #include <iostream>
 
50
#endif
 
51
 
 
52
   namespace smlib {
 
53
      using namespace std;
 
54
   
 
55
   
 
56
   ////////////////////////////////////
 
57
   // local variables, functions, etc
 
58
   ////////////////////////////////////
 
59
      namespace {
 
60
      
 
61
      //////////////////////////////////
 
62
      // C function with one parameter
 
63
      //--------------------------------
 
64
      
 
65
          double fac(double x)
 
66
         // Function return the factorial of integer belong to range 0..170.
 
67
         // in future, will implement aproximation for gamma function to evaluate for x > 170
 
68
         {
 
69
            double p;
 
70
            unsigned n;
 
71
         
 
72
            if((x <0) || (x>170)) // maximum admited where
 
73
               throw mathex::error("Error [fac()]: range error");
 
74
            else {
 
75
               n = static_cast<unsigned>(x+0.5);
 
76
               if(x < 2)
 
77
                  return 1;
 
78
               else if(n == 2)
 
79
                  return(2);
 
80
               for(p=1;1<n;n--)
 
81
                  p *= n;
 
82
               return p;
 
83
            }
 
84
         } // fac()
 
85
      
 
86
      // convert radian to degree
 
87
          double deg(double x)
 
88
         {
 
89
            return( x * 180.0 / M_PI);
 
90
         } // deg()
 
91
      
 
92
      // convert degree to radian
 
93
          double rad(double x)
 
94
         {
 
95
            return(x * M_PI / 180.0);
 
96
         } // rad()
 
97
      
 
98
      // return square
 
99
      /*    double sqr(double x)
 
100
      {
 
101
         return (x*x);
 
102
      } */
 
103
      // sqr()
 
104
      
 
105
      // return real part
 
106
          double trunc(double x)
 
107
         {
 
108
            return static_cast<long>(x);
 
109
         } // trunc()
 
110
      
 
111
      // return rounded value to int
 
112
          double round(double x)
 
113
         {
 
114
            return static_cast<long>(x+0.5);
 
115
         } // round
 
116
      
 
117
      // return 0 if x negative, 1 otherwise
 
118
          double step(double x)
 
119
         {
 
120
           if ( x < 0 ) return static_cast<long>(0.);
 
121
           return static_cast<long>(1.);
 
122
         } // step
 
123
      
 
124
      // return -1 if x negative, 1 otherwise
 
125
          double sign(double x)
 
126
         {
 
127
           if ( x < 0 ) return static_cast<long>(-1.);
 
128
           return static_cast<long>(1.);
 
129
         } // sign
 
130
      
 
131
      //////////////////////////////////////////////////////////
 
132
      // arithmetic operators
 
133
      //////////////////////////////////////////////////////////
 
134
      // operators was treated in Parsingr time, to check the
 
135
      //   operator level.
 
136
      //   unary operator is the function double op(double) and
 
137
      //   binary operator is function double op(double,double)
 
138
      //--------------------------------------------------------
 
139
      //   unary:
 
140
      //   "+" -> none
 
141
      //   "-" -> unary_minus
 
142
      //   binary:
 
143
      //   "+" -> binary_plus
 
144
      //   "-" -> binary_minus
 
145
      //   "*" -> binary_times
 
146
      //   "/" -> binary_divide
 
147
      //   "^" -> binary_power
 
148
      //   "%" -> fmod (defined in math.h)
 
149
      ///////////////////////////////////////////////////////////
 
150
      
 
151
      //////////////////////
 
152
      // C unary operator
 
153
      
 
154
          double unary_minus(double x)
 
155
         {
 
156
            return -x;
 
157
         } // unary_minus()
 
158
      
 
159
      //////////////////////////////////////////
 
160
      // C binary operators
 
161
      //----------------------------------------
 
162
      
 
163
          double binary_plus(double x, double y)
 
164
         {
 
165
            return x + y;
 
166
         } // binary_plus() 
 
167
      
 
168
          double binary_minus(double x, double y)
 
169
         {
 
170
            return x - y;
 
171
         } // binary_minus()
 
172
      
 
173
          double binary_times(double x, double y)
 
174
         {
 
175
            return x * y;
 
176
         } // binary_timess()
 
177
      
 
178
          double binary_divide(double x, double y)
 
179
         // in future, put more precisery checking, for overflow ?
 
180
         {
 
181
            if (y == 0)
 
182
               throw mathex::error("Error [binary_divide()]: divisin by zero");
 
183
            else
 
184
               return (x/y);
 
185
         } // binary_divide()
 
186
      
 
187
          double binary_power(double x, double y)
 
188
         {
 
189
            return pow(x,y);
 
190
         } // binary_power()
 
191
      
 
192
      /////////////////////////////////////////
 
193
      // pre-defined user defined functions
 
194
      //---------------------------------------
 
195
      
 
196
          double p_rand(vector <double> const &x)
 
197
         // rand() return value between [0,1)
 
198
         {
 
199
            if(x.size() != 0)
 
200
               throw mathex::error("Error [p_rand()]: can not use argument");
 
201
            return rand()/(RAND_MAX+1.0);
 
202
         } // p_rand()
 
203
      
 
204
      // maximum
 
205
          double p_max(vector <double> const & x)
 
206
         {
 
207
         
 
208
            double maxval=0;
 
209
         
 
210
            if(x.size() == 0)
 
211
               throw mathex::error("Error [p_max()]: No arguments");
 
212
            maxval = x[0];
 
213
            for(unsigned i=0; i<x.size(); i++)
 
214
               if(x[i] > maxval)
 
215
                  maxval = x[i];
 
216
         
 
217
            return maxval;
 
218
         } // p_max
 
219
      
 
220
      // minimum
 
221
          double p_min(vector <double> const & x)
 
222
         {
 
223
         
 
224
            double minval=0;
 
225
         
 
226
            if(x.size() == 0)
 
227
               throw mathex::error("Error [p_min()]: No arguments");
 
228
            minval = x[0];
 
229
            for(unsigned i=0; i<x.size(); i++)
 
230
               if(x[i] < minval)
 
231
                  minval = x[i];
 
232
         
 
233
            return minval;
 
234
         } // p_min
 
235
      
 
236
      // sum
 
237
          double p_sum(vector <double> const & x)
 
238
         {
 
239
         
 
240
            double sumval=0;
 
241
         
 
242
            if(x.size() == 0)
 
243
               throw mathex::error("Error [p_sum()]: No arguments");
 
244
            for(unsigned i=0; i<x.size(); i++)
 
245
               sumval += x[i];
 
246
         
 
247
            return sumval;
 
248
         } // p_sum
 
249
      
 
250
      // average
 
251
          double p_med(vector <double> const & x)
 
252
         {
 
253
            if(x.size() == 0)
 
254
               throw mathex::error("Error [p_med()]: No arguments");
 
255
            return p_sum(x)/x.size();
 
256
         } // p_med()
 
257
      
 
258
      ////////////////////////////////
 
259
      // function/constant tables
 
260
      ////////////////////////////////
 
261
      
 
262
      //////////////////////////////////
 
263
      // One parameter C function table
 
264
      
 
265
      // Record for one parameter C functions
 
266
         typedef 
 
267
             struct
 
268
            {
 
269
               const char *name;
 
270
               double  (*f)(double x); // one parameter function
 
271
            } CFUNCREC;
 
272
        
 
273
      ///////////////////////////////////////////////////////////   
 
274
      // One parameter internal C defined functions
 
275
      // the unary operator are stored on first part of table
 
276
      //---------------------------------------------------------
 
277
      // CAUTION:
 
278
      // if add unary operator, adjust the definiftion of NUM_UNARY_OP 
 
279
      #define NUM_UNARY_OP 1
 
280
         CFUNCREC cfunctable[] =
 
281
         {
 
282
         // name and funtion pointer
 
283
         // first part of table is unary operators
 
284
         {"-", unary_minus}, // unary operator
 
285
         // seccond part of table is functions
 
286
         { "abs",     fabs },
 
287
         { "Abs",     fabs },
 
288
         { "fabs",     fabs },
 
289
         { "Fabs",     fabs },
 
290
         { "acos",    acos },
 
291
         { "Acos",    acos },
 
292
         { "asin",    asin },
 
293
         { "Asin",    asin },
 
294
         { "atan",    atan },
 
295
         { "Atan",    atan },
 
296
         { "cos",     cos },
 
297
         { "Cos",     cos },
 
298
         { "cosh",    cosh },
 
299
         { "Cosh",    cosh },
 
300
         { "deg",     deg },   // added
 
301
         { "Deg",     deg },
 
302
         { "exp",     exp },
 
303
         { "Exp",     exp },
 
304
         { "fac",     fac },   // added
 
305
         { "Fac",     fac },
 
306
         { "log",     log },
 
307
         { "Log",     log },
 
308
         { "log10",   log10 },
 
309
         { "Log10",   log10 },
 
310
         // { "pow10",   pow10 } // in future, add it?
 
311
         { "rad",     rad },   // added 
 
312
         { "Rad",     rad },
 
313
         { "round",   round }, // added
 
314
         { "Round",   round },
 
315
         { "sign",    sign },
 
316
         { "Sign",    sign },
 
317
         { "sin",     sin },
 
318
         { "Sin",     sin },
 
319
         { "sinh",    sinh },
 
320
         { "Sinh",    sinh },
 
321
         // { "sqr",     sqr }, // added
 
322
         { "sqrt",    sqrt },
 
323
         { "Sqrt",    sqrt },
 
324
         { "step",    step },
 
325
         { "Step",    step },
 
326
         { "tan",     tan },
 
327
         { "Tan",     tan },
 
328
         { "tanh",    tanh },
 
329
         { "Tanh",    tanh },
 
330
         { "trunc",   trunc }, // added
 
331
         { "Trunc",   trunc },
 
332
         { "floor",   floor }, // largest integer not grather than x
 
333
         { "Floor",   floor },
 
334
         { "ceil",    ceil }, // smallest integer not less than x
 
335
         { "Ceil",    ceil },
 
336
         
 
337
         { 0, 0 } // 0 mark the end of table
 
338
         };
 
339
      
 
340
      /////////////////////
 
341
      // binary operators
 
342
      
 
343
      // binary operator's record
 
344
          class BINOPREC {
 
345
         public:
 
346
            char name;
 
347
            double (*f)(double, double);
 
348
         };
 
349
      
 
350
      ////////////////////
 
351
      // binary operators
 
352
         BINOPREC binoptable[] =
 
353
         {
 
354
         // name, value
 
355
         { '+',  binary_plus},
 
356
         { '-',  binary_minus},
 
357
         { '*',  binary_times},
 
358
         { '/',  binary_divide},
 
359
         { '^',  binary_power},
 
360
         { '%',  fmod},
 
361
         
 
362
         { '\0' , 0 } // '\0' mark the end of table
 
363
         };
 
364
      
 
365
      //////////////
 
366
      // constants
 
367
      
 
368
      // constants record
 
369
          class CONSTREC {
 
370
         public:
 
371
            const char *name;
 
372
            double value;
 
373
         };
 
374
      
 
375
      /////////////
 
376
      // constants
 
377
         CONSTREC consttable[] =
 
378
         {
 
379
         // name, value
 
380
         { "pi",  M_PI },
 
381
         { "e",   M_E },
 
382
         
 
383
         { NULL, 0 } // NULL mark the end of table
 
384
         };
 
385
      
 
386
      } // namespace {
 
387
   
 
388
   
 
389
   ////////////////////////////
 
390
   // implementations
 
391
   ///////////////////////////
 
392
   
 
393
   ///////////////
 
394
   // constatnts
 
395
   /// undefined number of arguments (for user defined functions
 
396
      const int mathex::UNDEFARGS = -1; // for user function arguments   
 
397
   
 
398
   ////////////////
 
399
   // methods
 
400
       void mathex::addstdfunc()
 
401
      {
 
402
         addfunc("rand", p_rand, 0);
 
403
         addfunc("Rand", p_rand, 0);
 
404
         addfunc("sum", p_sum, UNDEFARGS);
 
405
         addfunc("Sum", p_sum, UNDEFARGS);
 
406
         addfunc("max", p_max, UNDEFARGS);
 
407
         addfunc("Max", p_max, UNDEFARGS);
 
408
         addfunc("min", p_min, UNDEFARGS);
 
409
         addfunc("Min", p_min, UNDEFARGS);
 
410
         addfunc("med", p_med, UNDEFARGS);
 
411
         addfunc("Med", p_med, UNDEFARGS);
 
412
      } // addstdfunc()
 
413
   
 
414
       void mathex::reset()
 
415
      {
 
416
         delvar();
 
417
         delfunc();
 
418
         status = notparsed;
 
419
         expr = "";
 
420
         bytecode.clear();
 
421
         pos = 0;
 
422
         addstdfunc();  
 
423
      } // reset
 
424
   
 
425
   /////////////////////////////////
 
426
   // varibles table manipulations
 
427
   
 
428
       bool mathex::addvar(string const &name, double *var)
 
429
      // register the program internal variable
 
430
      {
 
431
         unsigned i;
 
432
      
 
433
         for(i=0; (i<vartable.size()) && (vartable[i].name != name);i++);
 
434
         if(i<vartable.size()) { // found! overwrite
 
435
            vartable[i].var = var;
 
436
            return true;
 
437
         }
 
438
         else if(!isnewvalidname(name))
 
439
            return false;
 
440
         vartable.push_back(VARREC(name, var));
 
441
         return true;
 
442
      } // addvar()
 
443
   
 
444
       bool mathex::delvar(string const &name)
 
445
      // delete one variable
 
446
      {
 
447
         unsigned i;
 
448
      
 
449
         for(i=0; (i<vartable.size()) && (vartable[i].name != name);i++);
 
450
         if(i < vartable.size())  {
 
451
         // how to use erase?
 
452
         // vartable.erase(&vartable[i],&vartable[i+1]);
 
453
            for(unsigned j=0; j<vartable.size()-1; j++)
 
454
               vartable[j] = vartable[j+1];
 
455
            vartable.pop_back(); // delete last
 
456
            status = notparsed;
 
457
            return true;
 
458
         }
 
459
         else
 
460
            return false;
 
461
      } // delvar()
 
462
   
 
463
   //////////////////////////////////////////////
 
464
   // user defined function table manipulation
 
465
   //////////////////////////////////////////////
 
466
   
 
467
       bool mathex::addfunc(string const &name, double (*f)(vector<double> const &), int NumArgs)
 
468
      // register the user defined function
 
469
      {
 
470
         unsigned i;
 
471
      
 
472
         for(i=0; (i<functable.size()) && (functable[i].name != name);i++);
 
473
         if(i<functable.size()) { // found! overwrite
 
474
            functable[i].f = f;
 
475
            functable[i].numargs = NumArgs;
 
476
            return true;
 
477
         }
 
478
         else if(!isnewvalidname(name))
 
479
            return false;
 
480
         functable.push_back(FUNCREC(name, f, NumArgs));
 
481
         return true;
 
482
      } // addfunc()
 
483
   
 
484
       bool mathex::delfunc(string const &name)
 
485
      // delete the user defined function
 
486
      {
 
487
         unsigned i;
 
488
      
 
489
         for(i=0; (i<functable.size()) && (functable[i].name != name);i++);
 
490
         if(i < functable.size())  {
 
491
         // how to use erase?
 
492
         // functable.erase(&vartable[i],&vartable[i+1]);
 
493
            for(unsigned j=0; j<vartable.size()-1; j++)
 
494
               functable[j] = functable[j+1];
 
495
            functable.pop_back(); // delete last
 
496
            return true;
 
497
         }
 
498
         else
 
499
            return false;
 
500
      } // delfunc()
 
501
   
 
502
   ////////////////////////////////////////////////////
 
503
   // get the index of variables/constants/functions/
 
504
   //     binary operator/user defined functions
 
505
   //--------------------------------------------------
 
506
   // return -1 if not found
 
507
   ////////////////////////////////////////////////////
 
508
   
 
509
       int mathex::getconst(string const &name)
 
510
      // get index of const
 
511
      // return -1 if not found
 
512
      {
 
513
         unsigned i;
 
514
      // look up the constants
 
515
         for(i=0;consttable[i].name && strcmp(name.c_str(), consttable[i].name);i++);
 
516
         if(consttable[i].name != NULL) // if found
 
517
            return i;
 
518
         else
 
519
            return -1;
 
520
      } // getconst
 
521
   
 
522
   
 
523
       int mathex::getvar(string const &name)
 
524
      // get index of variable
 
525
      // return -1 if not found
 
526
      {
 
527
         unsigned i;
 
528
      
 
529
      // look up the table
 
530
         for(i=0;(i<vartable.size()) && strcmp(name.c_str(), vartable[i].name.c_str());i++);
 
531
         if(i<vartable.size()) // if found
 
532
            return i;
 
533
         else
 
534
            return -1;
 
535
      } // getvar
 
536
   
 
537
   
 
538
       int mathex::getcfunc(string const &name)
 
539
      // get index of one parameter function
 
540
      // return -1 if not found
 
541
      {
 
542
         unsigned i;
 
543
      // look up the constants
 
544
         for(i=NUM_UNARY_OP;cfunctable[i].name && strcmp(name.c_str(), cfunctable[i].name);i++);
 
545
         if(cfunctable[i].name != NULL) // if found
 
546
            return i;
 
547
         else
 
548
            return -1;
 
549
      } // getcfunc
 
550
   
 
551
       int mathex::getunaryop(string const &name)
 
552
      // get index of unary operator
 
553
      // return -1 if not found
 
554
      {
 
555
         unsigned i;
 
556
      // look up the constants
 
557
         for(i=0; cfunctable[i].name && strcmp(name.c_str(), cfunctable[i].name);i++);
 
558
         if(cfunctable[i].name != NULL) // if found
 
559
            return i;
 
560
         else
 
561
            return -1;
 
562
      } // getunaryop
 
563
   
 
564
       int mathex::getbinop(char name)
 
565
      // get index of one parameter function
 
566
      // return -1 if not found
 
567
      {
 
568
         unsigned i;
 
569
      // look up the constants
 
570
         for(i=0;binoptable[i].name && (name != binoptable[i].name);i++);
 
571
         if(binoptable[i].name != '\0') // if found
 
572
            return i;
 
573
         else
 
574
            return -1;
 
575
      } // getbinop
 
576
   
 
577
   
 
578
       int mathex::getuserfunc(string const &name)
 
579
      // get index of variable
 
580
      // return -1 if not found
 
581
      {
 
582
         unsigned i;
 
583
      // look up the constants
 
584
         for(i=0;(i<functable.size()) && strcmp(name.c_str(), functable[i].name.c_str());i++);
 
585
         if(i<functable.size()) // if found
 
586
            return i;
 
587
         else
 
588
            return -1;
 
589
      } // getuserfunc
 
590
   
 
591
   
 
592
       bool mathex::isnewvalidname(string const &name)
 
593
      // Name validation.
 
594
      {
 
595
         if(name.empty() || (!isalpha(name[0]) && (name[0] != '_')))
 
596
            return false;
 
597
      // check for rest of characters
 
598
         for(unsigned j=0; j<name.size(); j++)
 
599
            if(!isalnum(name[j]) && (name[j] != '-'))
 
600
               return false;
 
601
         return (getcfunc(name) < 0) && (getconst(name) < 0) 
 
602
            && (getuserfunc(name) < 0) && (getvar(name) < 0);
 
603
      } // isnewvalidname
 
604
   
 
605
   
 
606
   //////////////////
 
607
   //  Evaluation 
 
608
   
 
609
   // main evaluator: internal use only
 
610
   
 
611
   // main evaluation function
 
612
       double mathex::eval()
 
613
      //  Eval the parsed stack and return
 
614
      {
 
615
         static vector <double> x; // suppose that eval does not eval
 
616
         evalstack.clear();
 
617
      
 
618
         if(status == notparsed) parse();
 
619
         if(status == invalid) throw error("eval()", "invalid expression");
 
620
      
 
621
         for(unsigned i=0; i<bytecode.size(); i++)
 
622
         {
 
623
            switch(bytecode[i].state) {
 
624
               case CODETOKEN::VALUE: evalstack.push_back(bytecode[i].value);
 
625
                  break;
 
626
               case CODETOKEN::VARIABLE:
 
627
               // get value of variable as value
 
628
                  evalstack.push_back(*vartable[bytecode[i].idx].var);
 
629
                  break;
 
630
               case CODETOKEN::FUNCTION: // Call the C internal functions with one parameter
 
631
               #ifdef _DEBUG_
 
632
               if(evalstack.size()<1) // error: It does not to occur if currect parsed.
 
633
                  throw error("eval()", "stack error");
 
634
               #endif
 
635
                  evalstack.back() = cfunctable[bytecode[i].idx].f(evalstack.back());
 
636
                  break;
 
637
               case CODETOKEN::BINOP: // Call the intern C function with two parameters
 
638
               #ifdef _DEBUG_
 
639
               if(evalstack.size() < 2) // error: It does not to occur if currect parsed.
 
640
                  throw error("eval()", "stack error");
 
641
               #endif
 
642
                  evalstack[evalstack.size()-2] = binoptable[bytecode[i].idx].f
 
643
                     (evalstack[evalstack.size()-2], evalstack.back());
 
644
                  evalstack.pop_back(); // delete last
 
645
                  break;                                
 
646
               case CODETOKEN::USERFUNC: // Call the user defined functions
 
647
               #ifdef _DEBUG_
 
648
               if(bytecode[i].numargs > evalstack.size())
 
649
                  throw error("eval()", "stack error");
 
650
               #endif
 
651
                  if(bytecode[i].numargs > 0) {
 
652
                     x.resize(bytecode[i].numargs);
 
653
                     for(unsigned j=0; j<static_cast<unsigned>(bytecode[i].numargs); j++)
 
654
                        x[bytecode[i].numargs-1-j] = evalstack[evalstack.size()-1-j];
 
655
                     evalstack.resize(evalstack.size()-bytecode[i].numargs+1);
 
656
                                                        
 
657
                     evalstack.back() = functable[bytecode[i].idx].f(x);
 
658
                  }
 
659
                                                else // Fixing bug pointed by  Hugh Denman <denmanh@tcd.ie> November 06, 2003
 
660
                                                   evalstack.push_back(functable[bytecode[i].idx].f(x));
 
661
                  break;         
 
662
               default: // invarid stack. It does not occur if currect parsed
 
663
                  throw  error("eval()", "invalid code token");
 
664
            }
 
665
         } // for(i=0; ByteCode[i].state != EMPTY;i++);
 
666
      
 
667
      #ifdef _DEBUG_
 
668
      if(evalstack.size() != 1)
 
669
         throw error("eval()", "stack error");  
 
670
      #endif
 
671
         return evalstack[0];
 
672
      } // eval()
 
673
   
 
674
   /////////////////
 
675
   // parser
 
676
   //---------------
 
677
   
 
678
   /////////////////
 
679
   // get the token
 
680
   
 
681
   // get the number token
 
682
       bool mathex::getnumber(double &x)
 
683
      {
 
684
         unsigned long i = pos;
 
685
         bool decimal;
 
686
      
 
687
      // is a number?
 
688
         if((i >= expr.size()) || !strchr("0123456789.", expr[i])) // not a number
 
689
            return false;
 
690
      
 
691
      // getting the number
 
692
         for(decimal=false; i<expr.size(); i++) {
 
693
            if(!isdigit(expr[i]) && ((expr[i] != '.') || decimal) )
 
694
               break;
 
695
            if(expr[i] == '.') decimal = true;
 
696
         }
 
697
      
 
698
         if((i==(pos+1)) && (expr[i]=='.')) // is not a number
 
699
            return false;
 
700
      
 
701
      // if scientific notation
 
702
         if((toupper(expr[i])=='E') && (i<expr.size())) { // cientific notation 
 
703
         // decimal = true; // turn on to detect that are double
 
704
            i++; // skip this
 
705
            if((i<expr.size()) && ((expr[i]=='+') || (expr[i]=='-')) ) { // if sign
 
706
               i++;
 
707
            }
 
708
         // get the expoent
 
709
            while((i<expr.size()) && isdigit(expr[i]))
 
710
               i++;
 
711
         }
 
712
      // now, copy the token and conter to number
 
713
      // if decimal is true, the number is double. otherwise, number is integer
 
714
      // for this detection, cientific notation need to enable decimal too.
 
715
      // The integer value are not used for this package
 
716
      
 
717
         x = strtod(expr.substr(pos, i-pos).c_str(), 0);
 
718
         pos = i;
 
719
         return true;
 
720
      } // getnumber()
 
721
   
 
722
       bool mathex::getidentifier(string &name)
 
723
      {
 
724
         unsigned i = pos;
 
725
      
 
726
         name.erase();
 
727
         if((i>=expr.size()) || (!isalpha(expr[i]) && (expr[i] != '_'))) // not a identifier
 
728
            return false;
 
729
      
 
730
      // identifier
 
731
         for(;(i<expr.size()) &&(isalnum(expr[i]) || (expr[i] == '_')); i++);
 
732
      
 
733
         name = expr.substr(pos, i-pos);
 
734
         pos = i; // advance the input
 
735
      
 
736
         return true;
 
737
      } // getidentifier()
 
738
   
 
739
   
 
740
       mathex::PARSERTOKEN::type mathex::nexttoken()
 
741
      // Gets the next token from the expr
 
742
      {
 
743
         string identifier;
 
744
      
 
745
         while((pos<expr.size()) && isspace(expr[pos]) )
 
746
            pos++;
 
747
      
 
748
         if(pos == expr.size()) {
 
749
            curtok.state = PARSERTOKEN::END;
 
750
            return curtok.state;
 
751
         }
 
752
      
 
753
         if(getnumber(curtok.value)) {
 
754
            curtok.state = PARSERTOKEN::VALUE;
 
755
            return curtok.state;
 
756
         }
 
757
         else if(getidentifier(identifier)) { // if identifier
 
758
         // checking the identifier type
 
759
            if((curtok.idx = getcfunc(identifier))>=0) // internal C functions
 
760
               curtok.state = PARSERTOKEN::FUNCTION;
 
761
            else if((curtok.idx=getuserfunc(identifier))>=0) { // user defined functions
 
762
               curtok.numargs = functable[curtok.idx].numargs;
 
763
               curtok.state = PARSERTOKEN::USERFUNC;
 
764
            }
 
765
            else if((curtok.idx=getvar(identifier))>=0) // variable
 
766
               curtok.state = PARSERTOKEN::VARIABLE;
 
767
            else if((curtok.idx=getconst(identifier))>=0) {  // constant are treated as NUMBER
 
768
               curtok.state = PARSERTOKEN::VALUE;
 
769
               curtok.value = consttable[curtok.idx].value;
 
770
            }
 
771
            else { // new identifier: not supported
 
772
               curtok.state = PARSERTOKEN::INVALID;
 
773
            }
 
774
         }
 
775
         else { // will be operators or delimiters
 
776
            switch(expr[pos]) {
 
777
               case '+' : curtok.state = PARSERTOKEN::PLUS;
 
778
                  break;
 
779
               case '-' : curtok.state = PARSERTOKEN::MINUS;
 
780
                  break;
 
781
               case '*' : curtok.state = PARSERTOKEN::TIMES;
 
782
                  break;
 
783
               case '/' : curtok.state = PARSERTOKEN::DIVIDE;
 
784
                  break;
 
785
               case '%' : curtok.state = PARSERTOKEN::MODULE;
 
786
                  break;
 
787
               case '^' : curtok.state = PARSERTOKEN::POWER;
 
788
                  break;
 
789
               case ',' : curtok.state = PARSERTOKEN::COMMA;
 
790
                  break;
 
791
               case '(' : curtok.state = PARSERTOKEN::OPAREN;
 
792
                  break;
 
793
               case ')' : curtok.state = PARSERTOKEN::CPAREN;
 
794
                  break;
 
795
               default  : curtok.state = PARSERTOKEN::INVALID;
 
796
            } // switch 
 
797
            if(curtok.state != PARSERTOKEN::INVALID) {
 
798
               curtok.idx = getbinop(expr[pos]);
 
799
               pos++;
 
800
            }
 
801
         } // else
 
802
      
 
803
         return curtok.state;
 
804
      } // nexttoken()
 
805
   
 
806
   ////////////////////////////
 
807
   // CodeStack operations 
 
808
   ////////////////////////////
 
809
   #ifdef _DEBUG_
 
810
    void mathex::printcoderec(CODETOKEN const &token)
 
811
   // print the Code Stack status
 
812
   {
 
813
      switch(token.state) {
 
814
      // case CODETOKEN::EMPTY: cout << "EMPTY\n";
 
815
      //  break;
 
816
         case CODETOKEN::VALUE:
 
817
            cout << "VALUE: " << token.value << endl;
 
818
            break;
 
819
         case CODETOKEN::VARIABLE:
 
820
            cout << "VARIABLE: " << vartable[token.idx].name << endl;
 
821
            break;
 
822
         case CODETOKEN::FUNCTION:
 
823
            if(token.idx <NUM_UNARY_OP) 
 
824
               cout << "unary operator: " << cfunctable[token.idx].name << endl;
 
825
            else
 
826
               cout << "FUNCTION: " << cfunctable[token.idx].name << endl;
 
827
            break;
 
828
         case CODETOKEN::BINOP:
 
829
            cout << "binary operator: " << binoptable[token.idx].name << endl;
 
830
            break;
 
831
         case CODETOKEN::USERFUNC:
 
832
            cout << "USERFUNC: " << functable[token.idx].name << endl;
 
833
            break;
 
834
         default: printf("INVALID\n");
 
835
      }
 
836
   
 
837
   } // printcoderec()
 
838
   
 
839
    void mathex::printbytecode()
 
840
   {
 
841
      cout << "codesize = "<< bytecode.size() << endl;
 
842
      for(unsigned i=0; i<bytecode.size(); i++)
 
843
         printcoderec(bytecode[i]);
 
844
   }
 
845
   #endif
 
846
   
 
847
       void mathex::parse()
 
848
      // Parse the expression
 
849
      {
 
850
      // parserstatus = true;
 
851
         bytecode.clear();
 
852
         status = invalid;
 
853
         pos=0;
 
854
         nexttoken();
 
855
         parsearithmetic1();
 
856
      #ifdef _DEBUG_
 
857
        printbytecode();
 
858
      #endif
 
859
         if(curtok.state != PARSERTOKEN::END) // if remain token
 
860
            throw error("parse()", "End of expression expected");
 
861
         status = parsed;
 
862
      } // parse()
 
863
   
 
864
       void mathex::parsearithmetic1(void)
 
865
      // level 1 arithmetic operator: binary plus/minus
 
866
      {
 
867
         unsigned savedidx; 
 
868
         parsearithmetic2();
 
869
         while((curtok.state == PARSERTOKEN::PLUS) || (curtok.state == PARSERTOKEN::MINUS)) {
 
870
            savedidx = curtok.idx;
 
871
            nexttoken();
 
872
            if((curtok.state == PARSERTOKEN::PLUS) || (curtok.state == PARSERTOKEN::MINUS)) 
 
873
               throw error("parse()", "Invalid expression");
 
874
            parsearithmetic2();      
 
875
            bytecode.push_back(CODETOKEN(CODETOKEN::BINOP, savedidx));
 
876
         } // while
 
877
      } // parsearithmetic1
 
878
   
 
879
   
 
880
       void mathex::parsearithmetic2(void)
 
881
      // level 2 arithmetic operator: multiplication, division, module
 
882
      {
 
883
         unsigned savedidx;
 
884
         parsearithmetic3();   
 
885
         while( (curtok.state == PARSERTOKEN::TIMES) || (curtok.state == PARSERTOKEN::DIVIDE)
 
886
           || (curtok.state == PARSERTOKEN::MODULE) ) {
 
887
            savedidx = curtok.idx;
 
888
            nexttoken();
 
889
            if((curtok.state == PARSERTOKEN::PLUS) || (curtok.state == PARSERTOKEN::MINUS)) 
 
890
               throw error("parse()", "Invalid expression");
 
891
            parsearithmetic3();
 
892
            bytecode.push_back(CODETOKEN(CODETOKEN::BINOP, savedidx));
 
893
         }
 
894
      } // parsearithmetic3
 
895
   
 
896
       void mathex::parsearithmetic3(void)
 
897
      // level 3 arithmetic operator: power
 
898
      {
 
899
         unsigned savedidx;
 
900
         parsearithmetic4();
 
901
         if(curtok.state == PARSERTOKEN::POWER) {
 
902
            savedidx = curtok.idx;
 
903
            nexttoken();
 
904
            if((curtok.state == PARSERTOKEN::PLUS) || (curtok.state == PARSERTOKEN::MINUS)) 
 
905
               throw error("parse()", "Invalid expression");
 
906
            parsearithmetic4();
 
907
            bytecode.push_back(CODETOKEN(CODETOKEN::BINOP, savedidx));
 
908
         }
 
909
      } // parsearithmetic3()
 
910
   
 
911
       void mathex::parsearithmetic4(void)
 
912
      // level 4 arithmetic operator: unary plus/minus
 
913
      {
 
914
         PARSERTOKEN::type state;
 
915
         if(((state=curtok.state) == PARSERTOKEN::PLUS) || (state == PARSERTOKEN::MINUS))
 
916
            nexttoken();
 
917
         if((curtok.state == PARSERTOKEN::PLUS) || (curtok.state == PARSERTOKEN::MINUS)) 
 
918
            throw error("parse()", "Invalid expression");
 
919
         parseatom();
 
920
      
 
921
         if(state == PARSERTOKEN::MINUS) // stored index are for binary operator. Get correct index
 
922
            bytecode.push_back(CODETOKEN(CODETOKEN::FUNCTION, getunaryop("-"))); 
 
923
        // unary minus are on function table
 
924
      } // parsearithmetic5()
 
925
   
 
926
       void mathex::parseatom(void)
 
927
      // level 6: literal numbers, variables and functions
 
928
      {
 
929
         unsigned i;
 
930
      
 
931
      // parentesis expression
 
932
         if(curtok.state == PARSERTOKEN::OPAREN) {
 
933
            nexttoken();
 
934
            if(curtok.state == PARSERTOKEN::CPAREN)
 
935
               throw error("parseatom()", "No expression inside parentesis");
 
936
            parsearithmetic1();
 
937
         
 
938
            if(curtok.state != PARSERTOKEN::CPAREN)
 
939
               throw error("parseatom()", "\")\" expected");
 
940
            nexttoken();  // Added by Hugh Denman (<denmanh@tcd.ie> or hdenman@cantab.net) Oct/03/2003     
 
941
         }
 
942
         // Number
 
943
         else if(curtok.state == PARSERTOKEN::VALUE) { // numbers
 
944
            bytecode.push_back(CODETOKEN(curtok.value));
 
945
            nexttoken();
 
946
         }
 
947
         // variables
 
948
         else if(curtok.state == PARSERTOKEN::VARIABLE) { // variables
 
949
            bytecode.push_back(CODETOKEN(CODETOKEN::VARIABLE, curtok.idx));
 
950
            nexttoken();
 
951
         }
 
952
         // internal C function with one parameters
 
953
         else if(curtok.state == PARSERTOKEN::FUNCTION) { // one parameter C functions
 
954
            parserstack.push_back(curtok);
 
955
            nexttoken();
 
956
            if(curtok.state != PARSERTOKEN::OPAREN)
 
957
               throw error("parseatom()", "\"(\" expected");
 
958
            nexttoken();
 
959
            if(curtok.state == PARSERTOKEN::CPAREN) 
 
960
               throw error("parseatom()", "invalid number of arguments");
 
961
            parsearithmetic1();
 
962
            if(curtok.state != PARSERTOKEN::CPAREN)
 
963
               throw error("parseatom()", "\")\" expected");
 
964
            curtok = parserstack.back();
 
965
            parserstack.pop_back();
 
966
            bytecode.push_back(CODETOKEN(CODETOKEN::FUNCTION, curtok.idx));
 
967
            nexttoken();
 
968
         }
 
969
         // user defined functions
 
970
         else if(curtok.state == PARSERTOKEN::USERFUNC) { // user defined function
 
971
            parserstack.push_back(curtok);
 
972
            nexttoken();
 
973
            if(curtok.state != PARSERTOKEN::OPAREN)
 
974
               throw error("parseatom()", "\"(\" expected");
 
975
            nexttoken();
 
976
            if(curtok.state == PARSERTOKEN::CPAREN)
 
977
               i = 0;
 
978
            else { // arguments exist
 
979
               parsearithmetic1();
 
980
               i = 1;
 
981
            // while(parserstatus && (curtok.state != PARSERTOKEN::END)
 
982
            //        &&(curtok.state != PARSERTOKEN::CPAREN)) {
 
983
               while((curtok.state != PARSERTOKEN::END) && (curtok.state != PARSERTOKEN::CPAREN)) {
 
984
                  if(curtok.state == PARSERTOKEN::COMMA) {
 
985
                     nexttoken();
 
986
                     i++;
 
987
                  }
 
988
                  else 
 
989
                     throw error("parseatom()", "unknow error");
 
990
                  parsearithmetic1();
 
991
               } // while
 
992
               if(curtok.state != PARSERTOKEN::CPAREN)
 
993
                  throw error("parseatom()", "\")\" expected");
 
994
            } // else
 
995
            curtok = parserstack.back();
 
996
            parserstack.pop_back();
 
997
         
 
998
            if ((curtok.numargs != UNDEFARGS) && (i != static_cast<unsigned>(curtok.numargs)))
 
999
            // i is current number of parameters
 
1000
               throw error("parseatom()", "invalid number of arguments");
 
1001
         
 
1002
         // number of parameters is correct. Now, put the function
 
1003
         // i is number of arguments    
 
1004
            bytecode.push_back(CODETOKEN(CODETOKEN::USERFUNC, curtok.idx, i));
 
1005
         
 
1006
            nexttoken();
 
1007
         } // user defined functions
 
1008
         // End of buffer
 
1009
         else if (curtok.state == PARSERTOKEN::END)
 
1010
            throw error("parseatom()", "unexpected end of expression");
 
1011
         // invalid
 
1012
         else if (curtok.state == PARSERTOKEN::INVALID)
 
1013
            throw error("parseatom()", "invalid token on expression");
 
1014
         // unknow error
 
1015
         else // it not occur
 
1016
            throw error("parseatom()", "unknow error");  
 
1017
      } // parseatom()
 
1018
   
 
1019
   
 
1020
   
 
1021
   } // namespace smlib {
 
1022
 
 
1023
 
 
1024
// end of mathex.c