1 /***************************************************************************
2 * eval.cpp is part of Math Graphic Library
3 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru> *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU Lesser General Public License as *
7 * published by the Free Software Foundation; either version 3 of the *
8 * License, or (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU Lesser General Public *
16 * License along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 #include <time.h>
21 #include <ctype.h>
22
23 #include "mgl2/data.h"
24 #include "mgl2/eval.h"
25
26 #if MGL_HAVE_GSL
27 #include <gsl/gsl_sf.h>
28 #include <gsl/gsl_rng.h>
29 #include <gsl/gsl_errno.h>
30 #include <sys/stat.h>
31 #endif
32 //-----------------------------------------------------------------------------
33 // constants for expression parsing
34 enum{
35 EQ_NUM=0, // a variable substitution
36 EQ_RND, // random number
37 EQ_A, // numeric constant
38 // normal functions of 2 arguments
39 EQ_LT, // comparison x<y !!! MUST BE FIRST 2-PLACE FUNCTION
40 EQ_GT, // comparison x>y
41 EQ_EQ, // comparison x=y
42 EQ_OR, // comparison x|y
43 EQ_AND, // comparison x&y
44 EQ_ADD, // addition x+y
45 EQ_SUB, // substraction x-y
46 EQ_MUL, // multiplication x*y
47 EQ_DIV, // division x/y
48 EQ_IPOW, // power x^n for integer n
49 EQ_POW, // power x^y
50 EQ_MOD, // x modulo y
51 EQ_LOG, // logarithm of x on base a, log_a(x) = ln(x)/ln(a)
52 EQ_ARG, // argument of complex number arg(x,y) = atan2(x,y)
53 EQ_HYPOT, // sqrt(x^2+y^2)=hypot(x,y)
54 EQ_MAX, // maximum of x and y
55 EQ_MIN, // minimum of x and y
56 // special functions of 2 arguments
57 EQ_BESJ, // regular cylindrical Bessel function of fractional order
58 EQ_BESY, // irregular cylindrical Bessel function of fractional order
59 EQ_BESI, // regular modified Bessel function of fractional order
60 EQ_BESK, // irregular modified Bessel function of fractional order
61 EQ_ELE, // elliptic integral E(\phi,k) = \int_0^\phi dt \sqrt((1 - k^2 \sin^2(t)))
62 EQ_ELF, // elliptic integral F(\phi,k) = \int_0^\phi dt 1/\sqrt((1 - k^2 \sin^2(t)))
63 EQ_LP, // Legendre polynomial P_l(x), (|x|<=1, l>=0)
64 EQ_BETA, // beta function B(x,y) = Gamma(x)*Gamma(y)/Gamma(x+y)
65 EQ_GAMMA_INC, // incomplete gamma function Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for x>=0.
66 // normal functions of 1 argument
67 EQ_SIN, // sine function \sin(x). !!! MUST BE FIRST 1-PLACE FUNCTION
68 EQ_COS, // cosine function \cos(x).
69 EQ_TAN, // tangent function \tan(x).
70 EQ_ASIN, // inverse sine function \sin(x).
71 EQ_ACOS, // inverse cosine function \sin(x).
72 EQ_ATAN, // inverse tangent function \tan(x).
73 EQ_SINH, // hyperbolic sine function \sin(x).
74 EQ_COSH, // hyperbolic cosine function \sin(x).
75 EQ_TANH, // hyperbolic tangent function \tan(x).
76 EQ_ASINH, // inverse hyperbolic sine function \sin(x).
77 EQ_ACOSH, // inverse hyperbolic cosine function \sin(x).
78 EQ_ATANH, // inverse hyperbolic tangent function \tan(x).
79 EQ_SQRT, // square root function \sqrt(x)
80 EQ_EXP, // exponential function \exp(x)
81 EQ_LN, // logarithm of x, ln(x)
82 EQ_LG, // decimal logarithm of x, lg(x) = ln(x)/ln(10)
83 EQ_SIGN, // sign of number
84 EQ_STEP, // step function
85 EQ_INT, // integer part [x]
86 EQ_ABS, // absolute value of x
87 // special functions of 1 argument
88 EQ_LI2, // dilogarithm for a real argument Li2(x) = - \Re \int_0^x ds \log(1-s)/s.
89 EQ_ELLE, // complete elliptic integral is denoted by E(k) = E(\pi/2, k).
90 EQ_ELLK, // complete elliptic integral is denoted by K(k) = F(\pi/2, k).
91 EQ_AI, // Airy function Ai(x)
92 EQ_BI, // Airy function Bi(x)
93 EQ_ERF, // error function erf(x) = (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2).
94 EQ_EI3, // exponential integral Ei_3(x) = \int_0^x dt \exp(-t^3) for x >= 0.
95 EQ_EI, // exponential integral Ei(x), Ei(x) := - PV(\int_{-x}^\infty dt \exp(-t)/t), where PV denotes the principal value of the integral.
96 EQ_E1, // exponential integral E_1(x), E_1(x) := Re \int_1^\infty dt \exp(-xt)/t.
97 EQ_E2, // exponential integral E_2(x), E_2(x) := Re \int_1^\infty dt \exp(-xt)/t^2.
98 EQ_SI, // Sine integral Si(x) = \int_0^x dt \sin(t)/t.
99 EQ_CI, // Cosine integral Ci(x) = \int_0^x dt \cos(t)/t.
100 EQ_GAMMA, // Gamma function \Gamma(x) = \int_0^\infty dt t^{x-1} \exp(-t)
101 EQ_PSI, // digamma function \psi(x) = \Gamma'(x)/\Gamma(x) for general x, x \ne 0.
102 EQ_W0, // principal branch of the Lambert W function, W_0(x). Functions W(x), are defined to be solutions of the equation W\exp(W) = x.
103 EQ_W1, // secondary real-valued branch of the Lambert W function, W_{-1}(x). Functions W(x), are defined to be solutions of the equation W\exp(W) = x.
104 EQ_SINC, // compute \sinc(x) = \sin(\pi x) / (\pi x) for any value of x.
105 EQ_ZETA, // Riemann zeta function \zeta(s) = \sum_{k=1}^\infty k^{-s}for arbitrary s, s \ne 1.
106 EQ_ETA, // eta function \eta(s) = (1-2^{1-s}) \zeta(s) for arbitrary s.
107 EQ_AID, // Derivative of Airy function Ai(x)
108 EQ_BID, // Derivative of Airy function Bi(x)
109 EQ_Z, // Dawson function \exp(-x^2) \int_0^x dt \exp(t^2)
110 // Jacoby functions of 2 arguments
111 EQ_SN, // Jacobian elliptic function sn(u|m) // !!! MUST BE FIRST NON 1-PLACE FUNCTION
112 EQ_SC, // Jacobian elliptic function sn(u|m)/cn(u|m)
113 EQ_SD, // Jacobian elliptic function sn(u|m)/dn(u|m)
114 EQ_NS, // Jacobian elliptic function 1/sn(u|m)
115 EQ_NC, // Jacobian elliptic function 1/cn(u|m)
116 EQ_ND, // Jacobian elliptic function 1/dn(u|m)
117 EQ_CN, // Jacobian elliptic function cn(u|m)
118 EQ_CS, // Jacobian elliptic function cn(u|m)/sn(u|m)
119 EQ_CD, // Jacobian elliptic function cn(u|m)/dn(u|m)
120 EQ_DN, // Jacobian elliptic function dn(u|m)
121 EQ_DS, // Jacobian elliptic function dn(u|m)/sn(u|m)
122 EQ_DC, // Jacobian elliptic function dn(u|m)/cn(u|m)
123 // MUST BE LAST ELLIPTIC FUNCTION
124 // not-ready
125 EQ_EN,
126 EQ_CL // Clausen function
127 };
128 //-----------------------------------------------------------------------------
129 #ifndef M_PI
130 #define M_PI 3.14159265358979323846
131 #endif
132 //-----------------------------------------------------------------------------
133 int mglFormula::Error=0;
134 bool MGL_LOCAL_PURE mglCheck(char *str,int n);
135 long MGL_LOCAL_PURE mglFindInText(const char *str, const char *lst);
136 //-----------------------------------------------------------------------------
137 #if MGL_HAVE_GSL
138 MGL_NO_EXPORT gsl_rng *mgl_rng=0; // NOTE: should be deleted by gsl_rng_free() but I don't know where :(
139 #endif
mgl_srnd(long seed)140 void MGL_EXPORT mgl_srnd(long seed)
141 {
142 #if MGL_HAVE_GSL
143 if(mgl_rng==0)
144 {
145 gsl_rng_env_setup();
146 mgl_rng = gsl_rng_alloc(gsl_rng_default);
147 }
148 gsl_rng_set(mgl_rng, seed);
149 #else
150 srand(seed);
151 #endif
152 }
mgl_srnd_(int * seed)153 void MGL_EXPORT mgl_srnd_(int *seed) { mgl_srnd(*seed); }
154 //-----------------------------------------------------------------------------
mgl_hypot(double x,double y)155 double MGL_EXPORT_CONST mgl_hypot(double x, double y) { return hypot(x,y); }
156 //-----------------------------------------------------------------------------
157 #if MGL_HAVE_PTHREAD
158 extern pthread_mutex_t mutexRnd;
159 #endif
mgl_rnd()160 double MGL_EXPORT mgl_rnd()
161 {
162 #if MGL_HAVE_PTHREAD
163 pthread_mutex_lock(&mutexRnd);
164 #endif
165 double res;
166 #pragma omp critical(rnd)
167 {
168 #if MGL_HAVE_GSL
169 if(mgl_rng==0)
170 {
171 gsl_rng_env_setup();
172 mgl_rng = gsl_rng_alloc(gsl_rng_default);
173 gsl_rng_set(mgl_rng, time(0));
174 }
175 res = gsl_rng_uniform(mgl_rng);
176 // gsl_rng_free(r);
177 #else
178 res = rand()/(RAND_MAX-1.);
179 #endif
180 }
181 #if MGL_HAVE_PTHREAD
182 pthread_mutex_unlock(&mutexRnd);
183 #endif
184 return res;
185 }
mgl_rnd_()186 double MGL_EXPORT mgl_rnd_() { return mgl_rnd(); }
187 //-----------------------------------------------------------------------------
~mglFormula()188 mglFormula::~mglFormula()
189 {
190 if(tmp) delete tmp;
191 if(Left) delete Left;
192 if(Right) delete Right;
193 }
194 //-----------------------------------------------------------------------------
195 // Formula constructor (automatically parse and "compile" formula)
mglFormula(const char * string)196 mglFormula::mglFormula(const char *string)
197 {
198 dat = tmp = NULL;
199 dx1=dy1=dz1=0; dx2=dy2=dz2=1;
200 #if MGL_HAVE_GSL
201 gsl_set_error_handler_off();
202 #endif
203 Error=0;
204 Left=Right=0;
205 Res=0; Kod=0;
206 if(!string) { Kod = EQ_NUM; Res = 0; return; }
207 char *str = new char[strlen(string)+1];
208 strcpy(str,string);
209 long n,len;
210 mgl_strtrim(str);
211 mgl_strlwr(str);
212 len=strlen(str);
213 if(str[0]==0) { delete []str; return; }
214 if(str[0]=='(' && mglCheck(str+1,len-2)) // remove braces
215 {
216 memmove(str,str+1,len);
217 len-=2; str[len]=0;
218 }
219 len=strlen(str);
220 if(str[0]==':' && str[1]!=0) // this data file for interpolation
221 {
222 double sx1,sx2,sy1,sy2,sz1,sz2;
223 char *buf = strchr(str+1,':');
224 if(buf && *buf)
225 {
226 *buf = 0;
227 int r = sscanf(buf+1,"%lg:%lg:%lg:%lg:%lg:%lg",&sx1,&sx2,&sy1,&sy2,&sz1,&sz2);
228 if(r>1 && sx1!=sx2) { dx1=sx1; dx2=sx2; }
229 if(r>3 && sy1!=sy2) { dy1=sy1; dy2=sy2; }
230 if(r>5 && sz1!=sz2) { dz1=sz1; dz2=sz2; }
231 }
232 dat = tmp = new mglData(str+1);
233 delete []str; return;
234 }
235 n=mglFindInText(str,"&|"); // lowest priority -- logical
236 if(n>=0)
237 {
238 if(str[n]=='|') Kod=EQ_OR; else Kod=EQ_AND;
239 str[n]=0;
240 Left=new mglFormula(str);
241 Right=new mglFormula(str+n+1);
242 delete []str; return;
243 }
244 n=mglFindInText(str,"<>="); // low priority -- conditions
245 if(n>=0)
246 {
247 if(str[n]=='<') Kod=EQ_LT;
248 else if(str[n]=='>') Kod=EQ_GT;
249 else Kod=EQ_EQ;
250 str[n]=0;
251 Left=new mglFormula(str);
252 Right=new mglFormula(str+n+1);
253 delete []str; return;
254 }
255 n=mglFindInText(str,"+-"); // normal priority -- additions
256 if(n>=0 && (n<2 || !strchr("eE",str[n-1]) || (str[n-2]!='.' && !isdigit(str[n-2]))))
257 {
258 if(str[n]=='+') Kod=EQ_ADD; else Kod=EQ_SUB;
259 str[n]=0;
260 Left=new mglFormula(str);
261 Right=new mglFormula(str+n+1);
262 delete []str; return;
263 }
264 n=mglFindInText(str,"*/%"); // high priority -- multiplications
265 if(n>=0)
266 {
267 if(str[n]=='*') Kod=EQ_MUL;
268 else if(str[n]=='/') Kod=EQ_DIV;
269 else Kod=EQ_MOD;
270 str[n]=0;
271 Left=new mglFormula(str);
272 Right=new mglFormula(str+n+1);
273 delete []str; return;
274 }
275 n=mglFindInText(str,"^"); // highest priority -- power
276 if(n>=0)
277 {
278 Kod=EQ_IPOW; str[n]=0;
279 Left=new mglFormula(str);
280 Right=new mglFormula(str+n+1);
281 delete []str; return;
282 }
283
284 for(n=0;n<len;n++) if(str[n]=='(') break;
285 if(n>=len) // this is number or variable
286 {
287 Kod = EQ_NUM;
288 // Left = Right = 0;
289 if(str[1]==0 && str[0]>='a' && str[0]<='z') // available variables
290 { Kod=EQ_A; Res = str[0]-'a'; }
291 else if(!strcmp(str,"rnd")) Kod=EQ_RND;
292 else if(!strcmp(str,"pi")) Res=M_PI;
293 else if(!strcmp(str,"inf")) Res=INFINITY;
294 else Res=atof(str); // this is number
295 }
296 else
297 {
298 char name[128];
299 mgl_strncpy(name,str,128); name[127]=name[n]=0;
300 memmove(str,str+n+1,len-n);
301 len=strlen(str); str[--len]=0;
302 if(!strncmp(name,"jacobi_",7))
303 memmove(name,name+7,(strlen(name+7)+1)*sizeof(char));
304 if(name[0]=='a')
305 {
306 if(!strcmp(name+1,"sin")) Kod=EQ_ASIN;
307 else if(!strcmp(name+1,"cos")) Kod=EQ_ACOS;
308 else if(!strcmp(name+1,"tan")) Kod=EQ_ATAN;
309 else if(!strcmp(name+1,"sinh")) Kod=EQ_ASINH;
310 else if(!strcmp(name+1,"cosh")) Kod=EQ_ACOSH;
311 else if(!strcmp(name+1,"tanh")) Kod=EQ_ATANH;
312 else if(!strcmp(name+1,"rg")) Kod=EQ_ARG;
313 else if(!strcmp(name+1,"bs")) Kod=EQ_ABS;
314 else if(!strcmp(name+1,"i")) Kod=EQ_AI;
315 else if(!strcmp(name+1,"iry_ai")) Kod=EQ_AI;
316 else if(!strcmp(name+1,"iry_bi")) Kod=EQ_BI;
317 else if(!strcmp(name+1,"iry_dai")) Kod=EQ_AID;
318 else if(!strcmp(name+1,"iry_dbi")) Kod=EQ_BID;
319 }
320 else if(name[0]=='b')
321 {
322 if(!strcmp(name+1,"i")) Kod=EQ_BI;
323 else if(!strcmp(name+1,"essel_j")) Kod=EQ_BESJ;
324 else if(!strcmp(name+1,"essel_i")) Kod=EQ_BESI;
325 else if(!strcmp(name+1,"essel_k")) Kod=EQ_BESK;
326 else if(!strcmp(name+1,"essel_y")) Kod=EQ_BESY;
327 else if(!strcmp(name+1,"eta")) Kod=EQ_BETA;
328 }
329 else if(name[0]=='c')
330 {
331 if(!strcmp(name+1,"os")) Kod=EQ_COS;
332 else if(!strcmp(name+1,"osh")) Kod=EQ_COSH;
333 else if(!strcmp(name+1,"h")) Kod=EQ_COSH;
334 else if(!strcmp(name+1,"i")) Kod=EQ_CI;
335 else if(!strcmp(name+1,"n")) Kod=EQ_CN;
336 else if(!strcmp(name+1,"s")) Kod=EQ_CS;
337 else if(!strcmp(name+1,"d")) Kod=EQ_CD;
338 else if(!strcmp(name+1,"l")) Kod=EQ_CL;
339 }
340 else if(name[0]=='d')
341 {
342 if(!strcmp(name+1,"n")) Kod=EQ_DN;
343 else if(!strcmp(name+1,"s")) Kod=EQ_DS;
344 else if(!strcmp(name+1,"c")) Kod=EQ_DC;
345 else if(!strcmp(name+1,"ilog")) Kod=EQ_LI2;
346 }
347 else if(name[0]=='e')
348 {
349 if(!strcmp(name+1,"xp")) Kod=EQ_EXP;
350 else if(!strcmp(name+1,"rf")) Kod=EQ_ERF;
351 else if(!strcmp(name+1,"n")) Kod=EQ_EN;
352 else if(!strcmp(name+1,"e")) Kod=EQ_ELLE;
353 else if(!strcmp(name+1,"k")) Kod=EQ_ELLK;
354 else if(name[0]==0) Kod=EQ_ELE;
355 else if(!strcmp(name+1,"i")) Kod=EQ_EI;
356 else if(!strcmp(name+1,"1")) Kod=EQ_E1;
357 else if(!strcmp(name+1,"2")) Kod=EQ_E2;
358 else if(!strcmp(name+1,"ta")) Kod=EQ_ETA;
359 else if(!strcmp(name+1,"i3")) Kod=EQ_EI3;
360 else if(!strcmp(name+1,"lliptic_e")) Kod=EQ_ELE;
361 else if(!strcmp(name+1,"lliptic_f")) Kod=EQ_ELF;
362 else if(!strcmp(name+1,"lliptic_ec")) Kod=EQ_ELLE;
363 else if(!strcmp(name+1,"lliptic_kc")) Kod=EQ_ELLK;
364 }
365 else if(name[0]=='l')
366 {
367 if(!strcmp(name+1,"og")) Kod=EQ_LOG;
368 else if(!strcmp(name+1,"g")) Kod=EQ_LG;
369 else if(!strcmp(name+1,"n")) Kod=EQ_LN;
370 else if(!strcmp(name+1,"i2")) Kod=EQ_LI2;
371 else if(!strcmp(name+1,"egendre")) Kod=EQ_LP;
372 }
373 else if(name[0]=='s')
374 {
375 if(!strcmp(name+1,"qrt")) Kod=EQ_SQRT;
376 else if(!strcmp(name+1,"in")) Kod=EQ_SIN;
377 else if(!strcmp(name+1,"tep")) Kod=EQ_STEP;
378 else if(!strcmp(name+1,"ign")) Kod=EQ_SIGN;
379 else if(!strcmp(name+1,"inh")) Kod=EQ_SINH;
380 else if(!strcmp(name+1,"h")) Kod=EQ_SINH;
381 else if(!strcmp(name+1,"i")) Kod=EQ_SI;
382 else if(!strcmp(name+1,"n")) Kod=EQ_SN;
383 else if(!strcmp(name+1,"c")) Kod=EQ_SC;
384 else if(!strcmp(name+1,"d")) Kod=EQ_SD;
385 else if(!strcmp(name+1,"inc")) Kod=EQ_SINC;
386 }
387 else if(name[0]=='t')
388 {
389 if(!strcmp(name+1,"g")) Kod=EQ_TAN;
390 else if(!strcmp(name+1,"an")) Kod=EQ_TAN;
391 else if(!strcmp(name+1,"anh")) Kod=EQ_TANH;
392 else if(!strcmp(name+1,"h")) Kod=EQ_TANH;
393 }
394 else if(name[0]=='m')
395 {
396 if(!strcmp(name+1,"od")) Kod=EQ_MOD;
397 else if(!strcmp(name+1,"ax")) Kod=EQ_MAX;
398 else if(!strcmp(name+1,"in")) Kod=EQ_MIN;
399 }
400 else if(!strcmp(name,"hypot")) Kod=EQ_HYPOT;
401 else if(!strcmp(name,"pow")) Kod=EQ_POW;
402 else if(!strcmp(name,"i")) Kod=EQ_BESI;
403 else if(!strcmp(name,"int")) Kod=EQ_INT;
404 else if(!strcmp(name,"j")) Kod=EQ_BESJ;
405 else if(!strcmp(name,"k")) Kod=EQ_BESK;
406 else if(!strcmp(name,"y")) Kod=EQ_BESY;
407 else if(!strcmp(name,"f")) Kod=EQ_ELF;
408 else if(!strcmp(name,"gamma")) Kod=EQ_GAMMA;
409 else if(!strcmp(name,"gamma_inc")) Kod=EQ_GAMMA_INC;
410 else if(!strcmp(name,"ns")) Kod=EQ_NS;
411 else if(!strcmp(name,"nc")) Kod=EQ_NC;
412 else if(!strcmp(name,"nd")) Kod=EQ_ND;
413 else if(!strcmp(name,"w0")) Kod=EQ_W0;
414 else if(!strcmp(name,"w1")) Kod=EQ_W1;
415 else if(!strcmp(name,"psi")) Kod=EQ_PSI;
416 else if(!strcmp(name,"zeta")) Kod=EQ_ZETA;
417 else if(!strcmp(name,"z")) Kod=EQ_Z;
418 else { delete []str; return; } // unknown function
419 n=mglFindInText(str,",");
420 if(n>=0)
421 {
422 str[n]=0;
423 Left=new mglFormula(str);
424 Right=new mglFormula(str+n+1);
425 }
426 else Left=new mglFormula(str);
427 }
428 delete []str;
429 }
430 //-----------------------------------------------------------------------------
431 // evaluate formula for 'x'='r', 'y'='n'='v', 't'='z', 'u'='a' variables
Calc(mreal x,mreal y,mreal t,mreal u) const432 mreal mglFormula::Calc(mreal x,mreal y,mreal t,mreal u) const
433 {
434 Error=0;
435 mreal a1[MGL_VS]; memset(a1,0,MGL_VS*sizeof(mreal));
436 a1['a'-'a'] = a1['c'-'a'] = a1['u'-'a'] = u;
437 a1['x'-'a'] = a1['r'-'a'] = x;
438 a1['y'-'a'] = a1['n'-'a'] = a1['v'-'a'] = y;
439 a1['z'-'a'] = a1['t'-'a'] = t;
440 mreal b = CalcIn(a1);
441 return mgl_isfin(b) ? b : NAN;
442 }
443 //-----------------------------------------------------------------------------
444 // evaluate formula for 'x'='r', 'y'='n', 't'='z', 'u'='a', 'v'='b', 'w'='c' variables
Calc(mreal x,mreal y,mreal t,mreal u,mreal v,mreal w) const445 mreal mglFormula::Calc(mreal x,mreal y,mreal t,mreal u,mreal v,mreal w) const
446 {
447 Error=0;
448 mreal a1[MGL_VS]; memset(a1,0,MGL_VS*sizeof(mreal));
449 a1['c'-'a'] = a1['w'-'a'] = w;
450 a1['b'-'a'] = a1['v'-'a'] = v;
451 a1['a'-'a'] = a1['u'-'a'] = u;
452 a1['x'-'a'] = a1['r'-'a'] = x;
453 a1['y'-'a'] = a1['n'-'a'] = y;
454 a1['z'-'a'] = a1['t'-'a'] = t;
455 mreal b = CalcIn(a1);
456 return mgl_isfin(b) ? b : NAN;
457 }
458 //-----------------------------------------------------------------------------
459 // evaluate formula for arbitrary set of variables
Calc(const mreal var[MGL_VS]) const460 mreal mglFormula::Calc(const mreal var[MGL_VS]) const
461 {
462 Error=0;
463 mreal b = CalcIn(var);
464 return mgl_isfin(b) ? b : NAN;
465 }
466 //-----------------------------------------------------------------------------
467 // evaluate formula for 'x'='r', 'y'='n'='v', 't'='z', 'u'='a' variables
CalcD(char diff,mreal x,mreal y,mreal t,mreal u) const468 mreal mglFormula::CalcD(char diff,mreal x,mreal y,mreal t,mreal u) const
469 {
470 Error=0;
471 mreal a1[MGL_VS]; memset(a1,0,MGL_VS*sizeof(mreal));
472 a1['a'-'a'] = a1['c'-'a'] = a1['u'-'a'] = u;
473 a1['x'-'a'] = a1['r'-'a'] = x;
474 a1['y'-'a'] = a1['n'-'a'] = a1['v'-'a'] = y;
475 a1['z'-'a'] = a1['t'-'a'] = t;
476 mreal b = CalcDIn(diff-'a', a1);
477 return mgl_isfin(b) ? b : NAN;
478 }
479 //-----------------------------------------------------------------------------
480 // evaluate formula for 'x'='r', 'y'='n', 't'='z', 'u'='a', 'v'='b', 'w'='c' variables
CalcD(char diff,mreal x,mreal y,mreal t,mreal u,mreal v,mreal w) const481 mreal mglFormula::CalcD(char diff,mreal x,mreal y,mreal t,mreal u,mreal v,mreal w) const
482 {
483 Error=0;
484 mreal a1[MGL_VS]; memset(a1,0,MGL_VS*sizeof(mreal));
485 a1['c'-'a'] = a1['w'-'a'] = w;
486 a1['b'-'a'] = a1['v'-'a'] = v;
487 a1['a'-'a'] = a1['u'-'a'] = u;
488 a1['x'-'a'] = a1['r'-'a'] = x;
489 a1['y'-'a'] = a1['n'-'a'] = y;
490 a1['z'-'a'] = a1['t'-'a'] = t;
491 mreal b = CalcDIn(diff-'a', a1);
492 return mgl_isfin(b) ? b : NAN;
493 }
494 //-----------------------------------------------------------------------------
495 // evaluate derivate of formula respect to 'diff' variable for arbitrary set of other variables
CalcD(const mreal var[MGL_VS],char diff) const496 mreal mglFormula::CalcD(const mreal var[MGL_VS], char diff) const
497 {
498 Error=0;
499 mreal b = CalcDIn(diff-'a', var);
500 return mgl_isfin(b) ? b : NAN;
501 }
502 //-----------------------------------------------------------------------------
cand(double a,double b)503 double MGL_LOCAL_CONST cand(double a,double b) {return a&&b?1:0;}
cor(double a,double b)504 double MGL_LOCAL_CONST cor(double a,double b) {return a||b?1:0;}
ceq(double a,double b)505 double MGL_LOCAL_CONST ceq(double a,double b) {return a==b?1:0;}
clt(double a,double b)506 double MGL_LOCAL_CONST clt(double a,double b) {return a<b?1:0;}
cgt(double a,double b)507 double MGL_LOCAL_CONST cgt(double a,double b) {return a>b?1:0;}
add(double a,double b)508 double MGL_LOCAL_CONST add(double a,double b) {return a+b;}
sub(double a,double b)509 double MGL_LOCAL_CONST sub(double a,double b) {return a-b;}
mul(double a,double b)510 double MGL_LOCAL_CONST mul(double a,double b) {return a&&b?a*b:0;}
del(double a,double b)511 double MGL_LOCAL_CONST del(double a,double b) {return b?a/b:NAN;}
ipw(double a,double b)512 double MGL_LOCAL_CONST ipw(double a,double b) {return fabs(b-int(b))<1e-5 ? mgl_ipow(a,int(b)) : pow(a,b);}
llg(double a,double b)513 double MGL_LOCAL_CONST llg(double a,double b) {return log(a)/log(b);}
514 #if MGL_HAVE_GSL
gslEllE(double a,double b)515 double MGL_LOCAL_CONST gslEllE(double a,double b) {return gsl_sf_ellint_E(a,b,GSL_PREC_SINGLE);}
gslEllF(double a,double b)516 double MGL_LOCAL_CONST gslEllF(double a,double b) {return gsl_sf_ellint_F(a,b,GSL_PREC_SINGLE);}
gslLegP(double a,double b)517 double MGL_LOCAL_CONST gslLegP(double a,double b) {return gsl_sf_legendre_Pl(int(a),b);}
gslEllEc(double a)518 double MGL_LOCAL_CONST gslEllEc(double a) {return gsl_sf_ellint_Ecomp(a,GSL_PREC_SINGLE);}
gslEllFc(double a)519 double MGL_LOCAL_CONST gslEllFc(double a) {return gsl_sf_ellint_Kcomp(a,GSL_PREC_SINGLE);}
gslAi(double a)520 double MGL_LOCAL_CONST gslAi(double a) {return gsl_sf_airy_Ai(a,GSL_PREC_SINGLE);}
gslBi(double a)521 double MGL_LOCAL_CONST gslBi(double a) {return gsl_sf_airy_Bi(a,GSL_PREC_SINGLE);}
gslAi_d(double a)522 double MGL_LOCAL_CONST gslAi_d(double a) {return gsl_sf_airy_Ai_deriv(a,GSL_PREC_SINGLE);}
gslBi_d(double a)523 double MGL_LOCAL_CONST gslBi_d(double a) {return gsl_sf_airy_Bi_deriv(a,GSL_PREC_SINGLE);}
524 #endif
sgn(double a)525 double MGL_LOCAL_CONST sgn(double a) {return a<0 ? -1:(a>0?1:0);}
stp(double a)526 double MGL_LOCAL_CONST stp(double a) {return a>0 ? 1:0;}
arg(double a,double b)527 double MGL_LOCAL_CONST arg(double a,double b) { return atan2(b,a); }
mgz1(double)528 double MGL_LOCAL_CONST mgz1(double) {return NAN;} // NOTE I think NAN value is more correct here than 0
mgz2(double,double)529 double MGL_LOCAL_CONST mgz2(double,double) {return NAN;} // NOTE I think NAN value is more correct here than 0
mgl_asinh(double x)530 double MGL_LOCAL_CONST mgl_asinh(double x) { return log(x+sqrt(x*x+1.)); }
mgl_acosh(double x)531 double MGL_LOCAL_CONST mgl_acosh(double x) { return x>1 ? log(x+sqrt(x*x-1.)) : NAN; }
mgl_atanh(double x)532 double MGL_LOCAL_CONST mgl_atanh(double x) { return fabs(x)<1 ? log((1.+x)/(1.-x))/2 : NAN; }
mgl_fmin(double a,double b)533 double MGL_LOCAL_CONST mgl_fmin(double a,double b) { return a > b ? b : a; }
mgl_fmax(double a,double b)534 double MGL_LOCAL_CONST mgl_fmax(double a,double b) { return a > b ? a : b; }
mgl_fmod(double a,double m)535 double MGL_LOCAL_CONST mgl_fmod(double a, double m) { return (a>=0)?fmod(a,m):fmod(a,m)+m; }
536 //-----------------------------------------------------------------------------
537 typedef double (*func_1)(double);
538 typedef double (*func_2)(double, double);
539 //-----------------------------------------------------------------------------
540 static const mreal z2[EQ_SIN-EQ_LT] = {3,3,3,3,0,3,3,0,0,0,0,0,NAN,3,3,3,3
541 #if MGL_HAVE_GSL
542 ,3,NAN, 3,NAN, 0,0,3,1,3
543 #else
544 ,0,0,0,0,0,0,0,0,0
545 #endif
546 };
547 static const func_2 f2[EQ_SIN-EQ_LT] = {clt,cgt,ceq,cor,cand,add,sub,mul,del,ipw,pow,mgl_fmod,llg,arg,hypot,mgl_fmax,mgl_fmin
548 #if MGL_HAVE_GSL
549 ,gsl_sf_bessel_Jnu,gsl_sf_bessel_Ynu,
550 gsl_sf_bessel_Inu,gsl_sf_bessel_Knu,
551 gslEllE,gslEllF,gslLegP,gsl_sf_beta,gsl_sf_gamma_inc
552 #else
553 ,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2
554 #endif
555 };
556 static const func_1 f1[EQ_SN-EQ_SIN] = {sin,cos,tan,asin,acos,atan,sinh,cosh,tanh,
557 mgl_asinh,mgl_acosh,mgl_atanh,sqrt,exp,log,log10,sgn,stp,floor,fabs
558 #if MGL_HAVE_GSL
559 ,gsl_sf_dilog,gslEllEc,gslEllFc,gslAi,gslBi,gsl_sf_erf,
560 gsl_sf_expint_3,gsl_sf_expint_Ei,gsl_sf_expint_E1,gsl_sf_expint_E2,
561 gsl_sf_Si,gsl_sf_Ci,gsl_sf_gamma,gsl_sf_psi,gsl_sf_lambert_W0,
562 gsl_sf_lambert_Wm1,gsl_sf_sinc,gsl_sf_zeta,gsl_sf_eta,gslAi_d,gslBi_d,
563 gsl_sf_dawson
564 #else
565 ,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,
566 mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1
567 #endif
568 };
569 //-----------------------------------------------------------------------------
570 // evaluation of embedded (included) expressions
CalcIn(const mreal * a1) const571 mreal mglFormula::CalcIn(const mreal *a1) const
572 {
573 if(dat)
574 {
575 mreal x = (a1['x'-'a']-dx1)*(dat->GetNx()-1)/(dx2-dx1);
576 mreal y = (a1['y'-'a']-dy1)*(dat->GetNy()-1)/(dy2-dy1);
577 mreal z = (a1['z'-'a']-dz1)*(dat->GetNz()-1)/(dz2-dz1);
578 return mgl_data_spline(dat,x,y,z);
579 }
580 if(Kod<EQ_LT)
581 {
582 if(Kod==EQ_RND) return mgl_rnd();
583 else return (Kod==EQ_A) ? a1[int(Res)] : Res;
584 }
585 double a = Left->CalcIn(a1);
586 if(mgl_isfin(a))
587 {
588 if(Kod<EQ_SIN)
589 {
590 // try to bypass calc b if a==0
591 if(a==0 && z2[Kod-EQ_LT]!=3) return z2[Kod-EQ_LT];
592 return Right?f2[Kod-EQ_LT](a, Right->CalcIn(a1)):NAN;
593 }
594 else if(Kod<EQ_SN) return f1[Kod-EQ_SIN](a);
595 #if MGL_HAVE_GSL
596 else if(Kod<=EQ_DC)
597 {
598 double sn=0, cn=0, dn=0;
599 gsl_sf_elljac_e(a,Right->CalcIn(a1), &sn, &cn, &dn);
600 switch(Kod)
601 {
602 case EQ_SN: return sn;
603 case EQ_SC: return sn/cn;
604 case EQ_SD: return sn/dn;
605 case EQ_CN: return cn;
606 case EQ_CS: return cn/sn;
607 case EQ_CD: return cn/dn;
608 case EQ_DN: return dn;
609 case EQ_DS: return dn/sn;
610 case EQ_DC: return dn/cn;
611 case EQ_NS: return 1./sn;
612 case EQ_NC: return 1./cn;
613 case EQ_ND: return 1./dn;
614 }
615 }
616 #endif
617 }
618 return NAN;
619 }
620 //-----------------------------------------------------------------------------
mgzz(double,double)621 double MGL_LOCAL_CONST mgzz(double,double) {return 0;}
mgp(double,double)622 double MGL_LOCAL_CONST mgp(double ,double ) {return 1;}
mgm(double,double)623 double MGL_LOCAL_CONST mgm(double ,double ) {return -1;}
mul1(double,double b)624 double MGL_LOCAL_CONST mul1(double ,double b) {return b;}
mul2(double a,double)625 double MGL_LOCAL_CONST mul2(double a,double ) {return a;}
div1(double,double b)626 double MGL_LOCAL_CONST div1(double ,double b) {return b?1/b:NAN;}
div2(double a,double b)627 double MGL_LOCAL_CONST div2(double a,double b) {return b?-a/(b*b):NAN;}
ipw1(double a,double b)628 double MGL_LOCAL_CONST ipw1(double a,double b) {return b*(fabs(b-int(b))<1e-5 ? mgl_ipow(a,int(b-1)) : pow(a,b-1));}
pow1(double a,double b)629 double MGL_LOCAL_CONST pow1(double a,double b) {return b*pow(a,b-1);}
pow2(double a,double b)630 double MGL_LOCAL_CONST pow2(double a,double b) {return log(a)*pow(a,b);}
llg1(double a,double b)631 double MGL_LOCAL_CONST llg1(double a,double b) {return 1/(a*log(b));}
llg2(double a,double b)632 double MGL_LOCAL_CONST llg2(double a,double b) {return -log(a)/(b*log(b)*log(b));}
cos_d(double a)633 double MGL_LOCAL_CONST cos_d(double a) {return -sin(a);}
tan_d(double a)634 double MGL_LOCAL_CONST tan_d(double a) {return 1./(cos(a)*cos(a));}
asin_d(double a)635 double MGL_LOCAL_CONST asin_d(double a) {return 1./sqrt(1.-a*a);}
acos_d(double a)636 double MGL_LOCAL_CONST acos_d(double a) {return -1./sqrt(1.-a*a);}
atan_d(double a)637 double MGL_LOCAL_CONST atan_d(double a) {return 1./(1.+a*a);}
tanh_d(double a)638 double MGL_LOCAL_CONST tanh_d(double a) {return 1./(cosh(a)*cosh(a));}
atanh_d(double a)639 double MGL_LOCAL_CONST atanh_d(double a){return 1./(1.-a*a);}
asinh_d(double a)640 double MGL_LOCAL_CONST asinh_d(double a){return 1./sqrt(1.+a*a);}
acosh_d(double a)641 double MGL_LOCAL_CONST acosh_d(double a){return 1./sqrt(a*a-1.);}
sqrt_d(double a)642 double MGL_LOCAL_CONST sqrt_d(double a) {return 0.5/sqrt(a);}
log10_d(double a)643 double MGL_LOCAL_CONST log10_d(double a){return M_LN10/a;}
log_d(double a)644 double MGL_LOCAL_CONST log_d(double a) {return 1./a;}
erf_d(double a)645 double MGL_LOCAL_CONST erf_d(double a) {return 2*exp(-a*a)/sqrt(M_PI);}
dilog_d(double a)646 double MGL_LOCAL_CONST dilog_d(double a){return log(a)/(1.-a);}
ei_d(double a)647 double MGL_LOCAL_CONST ei_d(double a) {return exp(a)/a;}
si_d(double a)648 double MGL_LOCAL_CONST si_d(double a) {return a?sin(a)/a:1;}
ci_d(double a)649 double MGL_LOCAL_CONST ci_d(double a) {return cos(a)/a;}
exp3_d(double a)650 double MGL_LOCAL_CONST exp3_d(double a) {return exp(-a*a*a);}
e1_d(double a)651 double MGL_LOCAL_CONST e1_d(double a) {return exp(-a)/a;}
sinc_d(double a)652 double MGL_LOCAL_CONST sinc_d(double a) {return a ? (cos(M_PI*a)/a-sin(M_PI*a)/(M_PI*a*a)) : 0;}
653 #if MGL_HAVE_GSL
e2_d(double a)654 double MGL_LOCAL_CONST e2_d(double a) {return -gsl_sf_expint_E1(a);}
gslJnuD(double a,double b)655 double MGL_LOCAL_CONST gslJnuD(double a,double b) {return 0.5*(gsl_sf_bessel_Jnu(a-1,b)-gsl_sf_bessel_Jnu(a+1,b));}
gslYnuD(double a,double b)656 double MGL_LOCAL_CONST gslYnuD(double a,double b) {return 0.5*(gsl_sf_bessel_Ynu(a-1,b)-gsl_sf_bessel_Ynu(a+1,b));}
gslKnuD(double a,double b)657 double MGL_LOCAL_CONST gslKnuD(double a,double b) {return -(a*gsl_sf_bessel_Knu(a,b)/b +gsl_sf_bessel_Knu(a-1,b));}
gslInuD(double a,double b)658 double MGL_LOCAL_CONST gslInuD(double a,double b) {return -(a*gsl_sf_bessel_Inu(a,b)/b -gsl_sf_bessel_Inu(a-1,b));}
gslEllE1(double a,double b)659 double MGL_LOCAL_CONST gslEllE1(double a,double b) {return sqrt(1.-sin(a)*sin(a)*b);}
gslEllE2(double a,double b)660 double MGL_LOCAL_CONST gslEllE2(double a,double b) {return (gsl_sf_ellint_E(a,b,GSL_PREC_SINGLE) - gsl_sf_ellint_F(a,b,GSL_PREC_SINGLE))/(2.*b);}
gslEllF1(double a,double b)661 double MGL_LOCAL_CONST gslEllF1(double a,double b) {return 1./sqrt(1.-sin(a)*sin(a)*b);}
gslEllF2(double a,double b)662 double MGL_LOCAL_CONST gslEllF2(double a,double b) {return (gsl_sf_ellint_E(a,b,GSL_PREC_SINGLE) - gsl_sf_ellint_F(a,b,GSL_PREC_SINGLE)*(1.-b))/(2*b*(1.-b)) - sin(2.*a)/(sqrt(1.-sin(a)*sin(a)*b)*2.*(1.-b));}
gslE_d(double a)663 double MGL_LOCAL_CONST gslE_d(double a) {return (gsl_sf_ellint_Ecomp(a,GSL_PREC_SINGLE) - gsl_sf_ellint_Kcomp(a,GSL_PREC_SINGLE))/(2.*a);}
gslK_d(double a)664 double MGL_LOCAL_CONST gslK_d(double a) {return (gsl_sf_ellint_Ecomp(a,GSL_PREC_SINGLE) - (1.-a)*gsl_sf_ellint_Kcomp(a,GSL_PREC_SINGLE))/(2.*a*(1.-a));}
gamma_d(double a)665 double MGL_LOCAL_CONST gamma_d(double a) {return gsl_sf_psi(a)*gsl_sf_gamma(a);}
666 #endif
ginc_d(double a,double x)667 double MGL_LOCAL_CONST ginc_d(double a, double x) {return -exp(-x)*pow(x,a-1);}
668 //-----------------------------------------------------------------------------
669 static const func_2 f21[EQ_SIN-EQ_LT] = {mgzz,mgzz,mgzz, mgzz,mgzz,mgp, mgp,mul1,div1, ipw1,pow1,mgp,llg1, mgz2,mgzz,mgzz
670 #if MGL_HAVE_GSL
671 ,mgz2,mgz2,mgz2, mgz2,gslEllE1,gslEllF1, mgz2,mgz2,mgz2
672 #else
673 ,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2
674 #endif
675 };
676 static const func_2 f22[EQ_SIN-EQ_LT] = {mgzz,mgzz,mgzz,mgzz,mgzz,mgp,mgm,mul2,div2,pow2,pow2,mgz2,llg2, mgz2,mgzz,mgzz
677 #if MGL_HAVE_GSL
678 ,gslJnuD,gslYnuD,gslInuD,gslKnuD,gslEllE2,gslEllF2,mgz2/*gslLegP*/,mgz2,ginc_d
679 #else
680 ,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2,mgz2
681 #endif
682 };
683 static const func_1 f11[EQ_SN-EQ_SIN] = {cos,cos_d,tan_d,asin_d,acos_d,atan_d,cosh,sinh,tanh_d,
684 asinh_d,acosh_d,atanh_d,sqrt_d,exp,log_d,log10_d,mgz1,mgz1,mgz1,sgn
685 #if MGL_HAVE_GSL
686 ,dilog_d,gslE_d,gslK_d,gslAi_d,gslBi_d,erf_d,exp3_d,ei_d,e1_d,e2_d,
687 si_d,ci_d,gamma_d,gsl_sf_psi_1,mgz1,mgz1,sinc_d,mgz1,mgz1,mgz1,mgz1,mgz1
688 #else
689 ,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,
690 mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1,mgz1
691 #endif
692 };
693 //-----------------------------------------------------------------------------
694 // evaluation of derivative of embedded (included) expressions
CalcDIn(int id,const mreal * a1) const695 mreal mglFormula::CalcDIn(int id, const mreal *a1) const
696 {
697 if(dat)
698 {
699 mreal x = (a1['x'-'a']-dx1)*(dat->GetNx()-1)/(dx2-dx1);
700 mreal y = (a1['y'-'a']-dy1)*(dat->GetNy()-1)/(dy2-dy1);
701 mreal z = (a1['z'-'a']-dz1)*(dat->GetNz()-1)/(dz2-dz1);
702 mreal dx,dy,dz, res=0;
703 mgl_data_spline_ext(dat,x,y,z,&dx,&dy,&dz);
704 if(id=='x'-'a') res = dx/(dat->GetNx()-1)*(dx2-dx1);
705 if(id=='y'-'a') res = dy/(dat->GetNy()-1)*(dy2-dy1);
706 if(id=='z'-'a') res = dz/(dat->GetNz()-1)*(dz2-dz1);
707 return res;
708 }
709 if(Kod==EQ_A && id==(int)Res) return 1;
710 else if(Kod<EQ_LT) return 0;
711 double a = Left->CalcIn(a1), d = Left->CalcDIn(id,a1);
712 if(mgl_isfin(a) && mgl_isfin(d))
713 {
714 if(Kod<EQ_SIN)
715 {
716 double b = Right?Right->CalcIn(a1):NAN;
717 double c = Right?Right->CalcDIn(id,a1):NAN;
718 // return mgl_isfin(b) ? (f21[Kod-EQ_LT](a,b)*d + f22[Kod-EQ_LT](a,b)*c) : NAN;
719 return mgl_isfin(b) ? (d?f21[Kod-EQ_LT](a,b)*d:0) + (c?f22[Kod-EQ_LT](a,b)*c:0) : NAN;
720 }
721 // else if(Kod<EQ_SN) return f11[Kod-EQ_SIN](a)*d;
722 else if(Kod<EQ_SN) return d?f11[Kod-EQ_SIN](a)*d:0;
723 #if MGL_HAVE_GSL
724 else if(Kod<=EQ_DC)
725 {
726 double sn=0, cn=0, dn=0, b = Right->CalcIn(a1);
727 if(mgl_isbad(b)) return NAN;
728 gsl_sf_elljac_e(a,b, &sn, &cn, &dn);
729 switch(Kod) // At this moment parse only differentiation or argument NOT mu !!!
730 {
731 case EQ_SN: return cn*dn*d;
732 case EQ_SC: return dn*d/(cn*cn);
733 case EQ_SD: return cn*d/(dn*dn);
734 case EQ_CN: return -dn*sn*d;
735 case EQ_CS: return dn*d/(sn*sn);
736 case EQ_CD: return (b-1)*d*sn/(dn*dn);
737 case EQ_DN: return -b*d*cn*sn;
738 case EQ_DS: return -cn*d/(sn*sn);
739 case EQ_DC: return (1-b)*sn*d/(cn*cn);
740 case EQ_NS: return -cn*dn*d/(sn*sn);
741 case EQ_NC: return dn*sn*d/(cn*cn);
742 case EQ_ND: return b*cn*sn*d/(dn*dn);
743 }
744 }
745 #endif
746 }
747 return NAN;
748 }
749 //-----------------------------------------------------------------------------
750 // Check braces correctness
mglCheck(char * str,int n)751 bool MGL_LOCAL_PURE mglCheck(char *str,int n)
752 {
753 long s = 0;
754 for(long i=0;i<n;i++)
755 {
756 if(str[i]=='(') s++;
757 if(str[i]==')') s--;
758 if(s<0) return false;
759 }
760 return (s==0) ? true : false;
761 }
762 //-----------------------------------------------------------------------------
763 // Try to find one of symbols lst in the string str
mglFindInText(const char * str,const char * lst)764 long MGL_LOCAL_PURE mglFindInText(const char *str, const char *lst)
765 {
766 long l=0,r=0,len=strlen(str);
767 for(long i=len-1;i>=0;i--)
768 {
769 if(str[i]=='(') l++;
770 if(str[i]==')') r++;
771 if(l==r && strchr(lst,str[i])) return i;
772 }
773 return -1;
774 }
775 //-----------------------------------------------------------------------------
mgl_create_expr(const char * expr)776 HMEX MGL_EXPORT mgl_create_expr(const char *expr) { return new mglFormula(expr); }
mgl_delete_expr(HMEX ex)777 void MGL_EXPORT mgl_delete_expr(HMEX ex) { if(ex) delete ex; }
mgl_expr_eval(HMEX ex,double x,double y,double z)778 double MGL_EXPORT_PURE mgl_expr_eval(HMEX ex, double x, double y,double z)
779 { return ex->Calc(x,y,z); }
mgl_expr_eval_v(HMEX ex,mreal * var)780 double MGL_EXPORT_PURE mgl_expr_eval_v(HMEX ex, mreal *var)
781 { return ex->Calc(var); }
mgl_expr_diff(HMEX ex,char dir,double x,double y,double z)782 double MGL_EXPORT_PURE mgl_expr_diff(HMEX ex, char dir, double x, double y,double z)
783 { return ex->CalcD(dir,x,y,z); }
mgl_expr_diff_v(HMEX ex,char dir,mreal * var)784 double MGL_EXPORT_PURE mgl_expr_diff_v(HMEX ex, char dir, mreal *var)
785 { return ex->CalcD(var, dir); }
786 //-----------------------------------------------------------------------------
mgl_create_expr_(const char * expr,int l)787 uintptr_t MGL_EXPORT mgl_create_expr_(const char *expr, int l)
788 { char *s=new char[l+1]; memcpy(s,expr,l); s[l]=0;
789 uintptr_t res = uintptr_t(mgl_create_expr(s));
790 delete []s; return res; }
mgl_delete_expr_(uintptr_t * ex)791 void MGL_EXPORT mgl_delete_expr_(uintptr_t *ex) { mgl_delete_expr((HMEX)ex); }
mgl_expr_eval_(uintptr_t * ex,mreal * x,mreal * y,mreal * z)792 double MGL_EXPORT_PURE mgl_expr_eval_(uintptr_t *ex, mreal *x, mreal *y, mreal *z)
793 { return mgl_expr_eval((HMEX) ex, *x,*y,*z); }
mgl_expr_diff_(uintptr_t * ex,const char * dir,mreal * x,mreal * y,mreal * z,int)794 double MGL_EXPORT_PURE mgl_expr_diff_(uintptr_t *ex, const char *dir, mreal *x, mreal *y, mreal *z, int)
795 { return mgl_expr_diff((HMEX) ex, *dir,*x,*y,*z); }
796 //-----------------------------------------------------------------------------
797