1 /***************************************************************************
2 * evalp.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 #include <wchar.h>
23 #include "mgl2/base.h"
24 #include "mgl2/parser.h"
25 #if MGL_HAVE_GSL
26 #include <gsl/gsl_sf.h>
27 #include <gsl/gsl_errno.h>
28 #endif
29 //-----------------------------------------------------------------------------
30 std::wstring mgl_trim_ws(const std::wstring &str);
31 HMDT MGL_NO_EXPORT mglFormulaCalc(const std::wstring &string, mglParser *arg, const std::vector<mglDataA*> &head);
32 HADT MGL_NO_EXPORT mglFormulaCalcC(const std::wstring &string, mglParser *arg, const std::vector<mglDataA*> &head);
33 HMDT MGL_NO_EXPORT mglFormulaCalcA(std::wstring string, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns);
34 HADT MGL_NO_EXPORT mglFormulaCalcAC(std::wstring string, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns);
35 //-----------------------------------------------------------------------------
mgl_formula_calc(const char * str,long n,...)36 HMDT MGL_EXPORT mgl_formula_calc(const char *str, long n, ...)
37 {
38 if(n<1) return NULL;
39 std::vector<mglDataA*> head;
40 va_list vl; va_start(vl,n);
41 for(long i=0;i<n;i++) head.push_back(va_arg(vl,mglDataA*));
42 va_end(vl);
43 return mglFormulaCalc(str, head);
44 }
45 //-----------------------------------------------------------------------------
mgl_formula_calc_c(const char * str,long n,...)46 HADT MGL_EXPORT mgl_formula_calc_c(const char *str, long n, ...)
47 {
48 if(n<1) return NULL;
49 std::vector<mglDataA*> head;
50 va_list vl; va_start(vl,n);
51 for(long i=0;i<n;i++) head.push_back(va_arg(vl,mglDataA*));
52 va_end(vl);
53 return mglFormulaCalcC(str, head);
54 }
55 //-----------------------------------------------------------------------------
mglFormulaCalc(const char * str,const std::vector<mglDataA * > & head)56 HMDT MGL_EXPORT mglFormulaCalc(const char *str, const std::vector<mglDataA*> &head)
57 {
58 if(!str || *str==0) return NULL;
59 std::wstring s;
60 for(long i=0;str[i];i++) s.push_back(str[i]);
61 return mglFormulaCalc(s,0,head);
62 }
63 //-----------------------------------------------------------------------------
mglFormulaCalcC(const char * str,const std::vector<mglDataA * > & head)64 HADT MGL_EXPORT mglFormulaCalcC(const char *str, const std::vector<mglDataA*> &head)
65 {
66 if(!str || *str==0) return NULL;
67 std::wstring s;
68 for(long i=0;str[i];i++) s.push_back(str[i]);
69 return mglFormulaCalcC(s,0,head);
70 }
71 //-----------------------------------------------------------------------------
mglApplyFunc(const std::wstring & str,mglParser * arg,const std::vector<mglDataA * > & head,double (* func)(double),const std::vector<std::wstring> & fns)72 HMDT MGL_NO_EXPORT mglApplyFunc(const std::wstring &str, mglParser *arg, const std::vector<mglDataA*> &head, double (*func)(double), const std::vector<std::wstring> &fns)
73 {
74 HMDT d = mglFormulaCalcA(str, arg, head, fns);
75 long n = d->GetNN(); mreal *dd=d->a;
76 #pragma omp parallel for
77 for(long i=0;i<n;i++) dd[i] = func(dd[i]);
78 return d;
79 }
80 //-----------------------------------------------------------------------------
81 #if MGL_HAVE_GSL
mglApplyFuncGSL(const std::wstring & str,mglParser * arg,const std::vector<mglDataA * > & head,double (* func)(double,gsl_mode_t),const std::vector<std::wstring> & fns)82 HMDT MGL_NO_EXPORT mglApplyFuncGSL(const std::wstring &str, mglParser *arg, const std::vector<mglDataA*> &head, double (*func)(double, gsl_mode_t), const std::vector<std::wstring> &fns)
83 {
84 HMDT d = mglFormulaCalcA(str, arg, head, fns);
85 long n = d->GetNN(); mreal *dd=d->a;
86 #pragma omp parallel for
87 for(long i=0;i<n;i++) dd[i] = func(dd[i],GSL_PREC_SINGLE);
88 return d;
89 }
90 #endif
91 //-----------------------------------------------------------------------------
mglApplyOper(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,double (* func)(double,double),const std::vector<std::wstring> & fns)92 HMDT MGL_NO_EXPORT mglApplyOper(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, double (*func)(double,double), const std::vector<std::wstring> &fns)
93 {
94 HMDT a = mglFormulaCalcA(a1,arg,head, fns), b = mglFormulaCalcA(a2,arg,head,fns), r,d;
95 long na = a->GetNN(), nb = b->GetNN(), nn;
96 if(na!=1) { r=a; d=b; nn=na; }
97 else { r=b; d=a; nn=nb; }
98 mreal va=a->a[0], vb=b->a[0], *aa=a->a, *bb=b->a, *cc=r->a;
99 if(na==nb)
100 #pragma omp parallel for
101 for(long i=0;i<nn;i++) cc[i] = func(aa[i], bb[i]);
102 else if(na==1)
103 #pragma omp parallel for
104 for(long i=0;i<nn;i++) cc[i] = func(va, bb[i]);
105 else
106 #pragma omp parallel for
107 for(long i=0;i<nn;i++) cc[i] = func(aa[i], vb);
108 mgl_delete_data(d); return r;
109 }
110 //-----------------------------------------------------------------------------
mglApplyOperAdd(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)111 HMDT MGL_NO_EXPORT mglApplyOperAdd(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
112 {
113 HMDT a = mglFormulaCalcA(a1,arg,head,fns), b = mglFormulaCalcA(a2,arg,head,fns), r,d;
114 long na = a->GetNN(), nb = b->GetNN(), nn;
115 if(na!=1) { r=a; d=b; nn=na; }
116 else { r=b; d=a; nn=nb; }
117 mreal *aa=r->a, *bb=d->a, v=bb[0];
118 if(na==nb)
119 #pragma omp parallel for
120 for(long i=0;i<nn;i++) aa[i] += bb[i];
121 else
122 #pragma omp parallel for
123 for(long i=0;i<nn;i++) aa[i] += v;
124 mgl_delete_data(d); return r;
125 }
126 //-----------------------------------------------------------------------------
mglApplyOperSub(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)127 HMDT MGL_NO_EXPORT mglApplyOperSub(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
128 {
129 HMDT a = mglFormulaCalcA(a1,arg,head,fns), b = mglFormulaCalcA(a2,arg,head,fns), r,d;
130 long na = a->GetNN(), nb = b->GetNN(), nn;
131 if(na!=1) { r=a; d=b; nn=na; }
132 else { r=b; d=a; nn=nb; }
133 mreal va=a->a[0], vb=b->a[0], *aa=a->a, *bb=b->a, *cc=r->a;
134 if(na==nb)
135 #pragma omp parallel for
136 for(long i=0;i<nn;i++) cc[i] = aa[i]-bb[i];
137 else if(na==1)
138 #pragma omp parallel for
139 for(long i=0;i<nn;i++) cc[i] = va-bb[i];
140 else
141 #pragma omp parallel for
142 for(long i=0;i<nn;i++) cc[i] = aa[i]-vb;
143 mgl_delete_data(d); return r;
144 }
145 //-----------------------------------------------------------------------------
mglApplyOperMul(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)146 HMDT MGL_NO_EXPORT mglApplyOperMul(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
147 {
148 HMDT a = mglFormulaCalcA(a1,arg,head,fns), b = mglFormulaCalcA(a2,arg,head,fns), r,d;
149 long na = a->GetNN(), nb = b->GetNN(), nn;
150 if(na!=1) { r=a; d=b; nn=na; }
151 else { r=b; d=a; nn=nb; }
152 mreal *aa=r->a, *bb=d->a, v=bb[0];
153 if(na==nb)
154 #pragma omp parallel for
155 for(long i=0;i<nn;i++) aa[i] *= bb[i];
156 else
157 #pragma omp parallel for
158 for(long i=0;i<nn;i++) aa[i] *= v;
159 mgl_delete_data(d); return r;
160 }
161 //-----------------------------------------------------------------------------
mglApplyOperDiv(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)162 HMDT MGL_NO_EXPORT mglApplyOperDiv(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
163 {
164 HMDT a = mglFormulaCalcA(a1,arg,head,fns), b = mglFormulaCalcA(a2,arg,head,fns), r,d;
165 long na = a->GetNN(), nb = b->GetNN(), nn;
166 if(na!=1) { r=a; d=b; nn=na; }
167 else { r=b; d=a; nn=nb; }
168 mreal va=a->a[0], vb=b->a[0], *aa=a->a, *bb=b->a, *cc=r->a;
169 if(na==nb)
170 #pragma omp parallel for
171 for(long i=0;i<nn;i++) cc[i] = bb[i]!=0?aa[i]/bb[i]:NAN;
172 else if(na==1)
173 #pragma omp parallel for
174 for(long i=0;i<nn;i++) cc[i] = bb[i]!=0?va/bb[i]:NAN;
175 else if(vb!=0)
176 #pragma omp parallel for
177 for(long i=0;i<nn;i++) cc[i] = aa[i]/vb;
178 else
179 #pragma omp parallel for
180 for(long i=0;i<nn;i++) cc[i] = NAN;
181 mgl_delete_data(d); return r;
182 }
183 //-----------------------------------------------------------------------------
mglApplyFuncC(const std::wstring & str,mglParser * arg,const std::vector<mglDataA * > & head,dual (* func)(dual),const std::vector<std::wstring> & fns)184 HADT MGL_NO_EXPORT mglApplyFuncC(const std::wstring &str, mglParser *arg, const std::vector<mglDataA*> &head, dual (*func)(dual), const std::vector<std::wstring> &fns)
185 {
186 HADT d = mglFormulaCalcAC(str, arg, head,fns);
187 long n = d->GetNN(); dual *dd=d->a;
188 #pragma omp parallel for
189 for(long i=0;i<n;i++) dd[i] = func(dd[i]);
190 return d;
191 }
192 //-----------------------------------------------------------------------------
mglApplyOperC(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,dual (* func)(dual,dual),const std::vector<std::wstring> & fns)193 HADT MGL_NO_EXPORT mglApplyOperC(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, dual (*func)(dual,dual), const std::vector<std::wstring> &fns)
194 {
195 HADT a = mglFormulaCalcAC(a1,arg,head,fns), b = mglFormulaCalcAC(a2,arg,head,fns), r,d;
196 long na = a->GetNN(), nb = b->GetNN(), nn;
197 if(na!=1) { r=a; d=b; nn=na; }
198 else { r=b; d=a; nn=nb; }
199 dual va=a->a[0], vb=b->a[0], *aa=a->a, *bb=b->a, *cc=r->a;
200 if(na==nb)
201 #pragma omp parallel for
202 for(long i=0;i<nn;i++) cc[i] = func(aa[i], bb[i]);
203 else if(na==1)
204 #pragma omp parallel for
205 for(long i=0;i<nn;i++) cc[i] = func(va, bb[i]);
206 else
207 #pragma omp parallel for
208 for(long i=0;i<nn;i++) cc[i] = func(aa[i], vb);
209 mgl_delete_datac(d); return r;
210 }
211 //-----------------------------------------------------------------------------
mglApplyOperAddC(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)212 HADT MGL_NO_EXPORT mglApplyOperAddC(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
213 {
214 HADT a = mglFormulaCalcAC(a1,arg,head,fns), b = mglFormulaCalcAC(a2,arg,head,fns), r,d;
215 long na = a->GetNN(), nb = b->GetNN(), nn;
216 if(na!=1) { r=a; d=b; nn=na; }
217 else { r=b; d=a; nn=nb; }
218 dual *aa=r->a, *bb=d->a, v=bb[0];
219 if(na==nb)
220 #pragma omp parallel for
221 for(long i=0;i<nn;i++) aa[i] += bb[i];
222 else
223 #pragma omp parallel for
224 for(long i=0;i<nn;i++) aa[i] += v;
225 mgl_delete_datac(d); return r;
226 }
227 //-----------------------------------------------------------------------------
mglApplyOperSubC(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)228 HADT MGL_NO_EXPORT mglApplyOperSubC(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
229 {
230 HADT a = mglFormulaCalcAC(a1,arg,head,fns), b = mglFormulaCalcAC(a2,arg,head,fns), r,d;
231 long na = a->GetNN(), nb = b->GetNN(), nn;
232 if(na!=1) { r=a; d=b; nn=na; }
233 else { r=b; d=a; nn=nb; }
234 dual va=a->a[0], vb=b->a[0], *aa=a->a, *bb=b->a, *cc=r->a;
235 if(na==nb)
236 #pragma omp parallel for
237 for(long i=0;i<nn;i++) cc[i] = aa[i]-bb[i];
238 else if(na==1)
239 #pragma omp parallel for
240 for(long i=0;i<nn;i++) cc[i] = va-bb[i];
241 else
242 #pragma omp parallel for
243 for(long i=0;i<nn;i++) cc[i] = aa[i]-vb;
244 mgl_delete_datac(d); return r;
245 }
246 //-----------------------------------------------------------------------------
mglApplyOperMulC(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)247 HADT MGL_NO_EXPORT mglApplyOperMulC(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
248 {
249 HADT a = mglFormulaCalcAC(a1,arg,head,fns), b = mglFormulaCalcAC(a2,arg,head,fns), r,d;
250 long na = a->GetNN(), nb = b->GetNN(), nn;
251 if(na!=1) { r=a; d=b; nn=na; }
252 else { r=b; d=a; nn=nb; }
253 dual *aa=r->a, *bb=d->a, v=bb[0];
254 if(na==nb)
255 #pragma omp parallel for
256 for(long i=0;i<nn;i++) aa[i] *= bb[i];
257 else
258 #pragma omp parallel for
259 for(long i=0;i<nn;i++) aa[i] *= v;
260 mgl_delete_datac(d); return r;
261 }
262 //-----------------------------------------------------------------------------
mglApplyOperDivC(const std::wstring & a1,const std::wstring & a2,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)263 HADT MGL_NO_EXPORT mglApplyOperDivC(const std::wstring &a1, const std::wstring &a2, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
264 {
265 HADT a = mglFormulaCalcAC(a1,arg,head,fns), b = mglFormulaCalcAC(a2,arg,head,fns), r,d;
266 long na = a->GetNN(), nb = b->GetNN(), nn;
267 if(na!=1) { r=a; d=b; nn=na; }
268 else { r=b; d=a; nn=nb; }
269 dual va=a->a[0], vb=b->a[0], *aa=a->a, *bb=b->a, *cc=r->a;
270 if(na==nb)
271 #pragma omp parallel for
272 for(long i=0;i<nn;i++) cc[i] = bb[i]!=mreal(0)?aa[i]/bb[i]:NAN;
273 else if(na==1)
274 #pragma omp parallel for
275 for(long i=0;i<nn;i++) cc[i] = bb[i]!=mreal(0)?va/bb[i]:NAN;
276 else if(vb!=mreal(0))
277 #pragma omp parallel for
278 for(long i=0;i<nn;i++) cc[i] = aa[i]/vb;
279 else
280 #pragma omp parallel for
281 for(long i=0;i<nn;i++) cc[i] = NAN;
282 mgl_delete_datac(d); return r;
283 }
284 //-----------------------------------------------------------------------------
mglCheck(std::wstring str)285 bool MGL_LOCAL_PURE mglCheck(std::wstring str)
286 {
287 long s = 0,i,n=str.length();
288 for(i=0;i<n;i++)
289 {
290 if(str[i]=='(') s++;
291 if(str[i]==')') s--;
292 if(s<0) return false;
293 }
294 return (s==0) ? true : false;
295 }
296 //-----------------------------------------------------------------------------
mglFindInText(const std::wstring & str,const char * lst,int num=0)297 long MGL_LOCAL_PURE mglFindInText(const std::wstring &str,const char *lst, int num=0)
298 {
299 long l=0,r=0,ls=0,rs=0;
300 for(long i=str.length()-1;i>=0;i--)
301 {
302 if(str[i]=='(') l++;
303 if(str[i]==')') r++;
304 if(str[i]=='[') ls++;
305 if(str[i]==']') rs++;
306 if(l==r && ls==rs && strchr(lst,str[i]))
307 { num--; if(num<0) return i; }
308 }
309 return -1;
310 }
311 //-----------------------------------------------------------------------------
312 double MGL_LOCAL_CONST cand(double a,double b);// {return a&&b?1:0;}
313 double MGL_LOCAL_CONST cor(double a,double b);// {return a||b?1:0;}
314 double MGL_LOCAL_CONST ceq(double a,double b);// {return a==b?1:0;}
315 double MGL_LOCAL_CONST clt(double a,double b);// {return a<b?1:0;}
316 double MGL_LOCAL_CONST cgt(double a,double b);// {return a>b?1:0;}
317 double MGL_LOCAL_CONST stp(double a);// {return a>0?1:0;}
318 double MGL_LOCAL_CONST sgn(double a);// {return a>0?1:(a<0?-1:0);}
319 double MGL_LOCAL_CONST ipw(double a,double b);// {return mgl_ipow(a,int(b));}
320 double MGL_LOCAL_CONST llg(double a,double b);// {return log(a)/log(b);}
321 //double MGL_LOCAL_CONST asinh(double x);// { return log(x+sqrt(x*x+1)); }
322 //double MGL_LOCAL_CONST acosh(double x);// { return x>1 ? log(x+sqrt(x*x-1)) : NAN; }
323 //double MGL_LOCAL_CONST atanh(double x);// { return fabs(x)<1 ? log((1+x)/(1-x))/2 : NAN; }
324 double MGL_LOCAL_CONST gslEllE(double a,double b);// {return gsl_sf_ellint_E(a,b,GSL_PREC_SINGLE);}
325 double MGL_LOCAL_CONST gslEllF(double a,double b);// {return gsl_sf_ellint_F(a,b,GSL_PREC_SINGLE);}
326 double MGL_LOCAL_CONST gslLegP(double a,double b);// {return gsl_sf_legendre_Pl(int(a),b);}
327 double MGL_LOCAL_CONST mgl_asinh(double x);
328 double MGL_LOCAL_CONST mgl_acosh(double x);
329 double MGL_LOCAL_CONST mgl_atanh(double x);
330 double MGL_LOCAL_CONST mgl_fmin(double a,double b);
331 double MGL_LOCAL_CONST mgl_fmax(double a,double b);
332 double MGL_LOCAL_CONST mgl_fmod(double a, double m);
333 //-----------------------------------------------------------------------------
334 // It seems that standard wcstombs() have a bug. So, I replace by my own.
mgl_wcstombs(char * dst,const wchar_t * src,int size)335 void MGL_EXPORT mgl_wcstombs(char *dst, const wchar_t *src, int size)
336 {
337 int j;
338 for(j=0;j<size-1 && src[j]!=0;j++)
339 dst[j] = src[j]<0x7f ? src[j] : ' ';
340 dst[j] = 0;
341 }
342 //-----------------------------------------------------------------------------
FindVar(const std::vector<mglDataA * > & head,const std::wstring & name)343 MGL_LOCAL_PURE const mglDataA *FindVar(const std::vector<mglDataA*> &head, const std::wstring &name)
344 {
345 for(size_t i=0;i<head.size();i++)
346 if(head[i] && !wcscmp(head[i]->Name(),name.c_str())) // bypass std::string comparison warning
347 return head[i];
348 return 0;
349 }
350 //-----------------------------------------------------------------------------
mgl_wcslwr(wchar_t * str)351 void MGL_EXPORT mgl_wcslwr(wchar_t *str)
352 {
353 size_t l=mgl_wcslen(str);
354 for(size_t k=0;k<l;k++)
355 str[k] = (str[k]>='A' && str[k]<='Z') ? str[k]+'a'-'A' : str[k];
356 }
357 //-----------------------------------------------------------------------------
mgl_gettime(const std::wstring & s)358 mreal MGL_NO_EXPORT mgl_gettime(const std::wstring &s)
359 {
360 mreal t=NAN;
361 tm a; memset(&a,0,sizeof(tm));
362 if(swscanf(s.c_str(),L"%u-%u-%u_%u.%u.%d", &a.tm_hour,&a.tm_min,&a.tm_sec, &a.tm_mday,&a.tm_mon,&a.tm_year)==6)
363 { a.tm_year-=1900; a.tm_mon -= 1;
364 if(a.tm_hour<24 && a.tm_min<60 && a.tm_sec<60 && a.tm_mday>0 && a.tm_mday<32 && a.tm_mon<12)
365 t = mktime(&a);
366 }
367 else if(swscanf(s.c_str(),L"%d.%d.%d", &a.tm_mday,&a.tm_mon,&a.tm_year)==3)
368 { a.tm_year-=1900; a.tm_mon -= 1;
369 if(a.tm_mday>0 && a.tm_mday<32 && a.tm_mon<12)
370 t = mktime(&a);
371 }
372 else if(swscanf(s.c_str(),L"%d-%d-%d", &a.tm_hour,&a.tm_min,&a.tm_sec)==3)
373 { a.tm_mday=1; a.tm_mon=0; a.tm_year=70;
374 if(a.tm_hour<24 && a.tm_min<60 && a.tm_sec<60)
375 t = mktime(&a);
376 }
377 return t;
378 }
379 //-----------------------------------------------------------------------------
mgl_jac_sn(double a,double m)380 double MGL_NO_EXPORT mgl_jac_sn(double a, double m)
381 {
382 double sn=0, cn=0, dn=0;
383 #if MGL_HAVE_GSL
384 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
385 #endif
386 return sn;
387 }
mgl_jac_sc(double a,double m)388 double MGL_NO_EXPORT mgl_jac_sc(double a, double m)
389 {
390 double sn=0, cn=1, dn=0;
391 #if MGL_HAVE_GSL
392 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
393 #endif
394 return sn/cn;
395 }
mgl_jac_sd(double a,double m)396 double MGL_NO_EXPORT mgl_jac_sd(double a, double m)
397 {
398 double sn=0, cn=0, dn=1;
399 #if MGL_HAVE_GSL
400 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
401 #endif
402 return sn/dn;
403 }
404
mgl_jac_cn(double a,double m)405 double MGL_NO_EXPORT mgl_jac_cn(double a, double m)
406 {
407 double sn=1, cn=0, dn=0;
408 #if MGL_HAVE_GSL
409 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
410 #endif
411 return cn;
412 }
mgl_jac_cs(double a,double m)413 double MGL_NO_EXPORT mgl_jac_cs(double a, double m)
414 {
415 double sn=1, cn=0, dn=0;
416 #if MGL_HAVE_GSL
417 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
418 #endif
419 return cn/sn;
420 }
mgl_jac_cd(double a,double m)421 double MGL_NO_EXPORT mgl_jac_cd(double a, double m)
422 {
423 double sn=0, cn=0, dn=1;
424 #if MGL_HAVE_GSL
425 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
426 #endif
427 return cn/dn;
428 }
429
mgl_jac_dn(double a,double m)430 double MGL_NO_EXPORT mgl_jac_dn(double a, double m)
431 {
432 double sn=0, cn=0, dn=0;
433 #if MGL_HAVE_GSL
434 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
435 #endif
436 return dn;
437 }
mgl_jac_ds(double a,double m)438 double MGL_NO_EXPORT mgl_jac_ds(double a, double m)
439 {
440 double sn=1, cn=0, dn=0;
441 #if MGL_HAVE_GSL
442 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
443 #endif
444 return dn/sn;
445 }
mgl_jac_dc(double a,double m)446 double MGL_NO_EXPORT mgl_jac_dc(double a, double m)
447 {
448 double sn=0, cn=1, dn=0;
449 #if MGL_HAVE_GSL
450 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
451 #endif
452 return dn/cn;
453 }
454
mgl_jac_nd(double a,double m)455 double MGL_NO_EXPORT mgl_jac_nd(double a, double m)
456 {
457 double sn=0, cn=0, dn=1;
458 #if MGL_HAVE_GSL
459 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
460 #endif
461 return 1./dn;
462 }
mgl_jac_ns(double a,double m)463 double MGL_NO_EXPORT mgl_jac_ns(double a, double m)
464 {
465 double sn=1, cn=0, dn=0;
466 #if MGL_HAVE_GSL
467 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
468 #endif
469 return 1./sn;
470 }
mgl_jac_nc(double a,double m)471 double MGL_NO_EXPORT mgl_jac_nc(double a, double m)
472 {
473 double sn=0, cn=1, dn=0;
474 #if MGL_HAVE_GSL
475 gsl_sf_elljac_e(a,m, &sn, &cn, &dn);
476 #endif
477 return 1./cn;
478 }
479 //-----------------------------------------------------------------------------
480 /// Parse string and substitute the script argument
481 // All numbers are presented as mglData(1). Do boundary checking.
482 // NOTE: In any case where number is required the mglData::a[0] is used.
483 // String flag is binary 0x1 -> 'x', 0x2 -> 'y', 0x4 -> 'z'
mglFormulaCalc(const std::wstring & str,mglParser * arg,const std::vector<mglDataA * > & head)484 HMDT MGL_NO_EXPORT mglFormulaCalc(const std::wstring &str, mglParser *arg, const std::vector<mglDataA*> &head)
485 {
486 if(str.empty()) return NULL;
487 std::vector<std::wstring> fns = mgl_wcs_args(str,'\\');
488 return mglFormulaCalcA(fns[0],arg,head,fns);
489 }
mglFormulaCalcA(std::wstring str,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)490 HMDT MGL_NO_EXPORT mglFormulaCalcA(std::wstring str, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
491 {
492 #if MGL_HAVE_GSL
493 gsl_set_error_handler_off();
494 #endif
495 if(str.empty()) return new mglData; // nothing to parse
496 str = mgl_trim_ws(str);
497 mreal tval = mgl_gettime(str);
498 if(mgl_isnum(tval))
499 { mglData *r=new mglData; r->a[0] = tval; return r; }
500
501 long n,len=str.length();
502 if(str[0]=='(' && mglCheck(str.substr(1,len-2))) // remove braces
503 { str = str.substr(1,len-2); len-=2; }
504 if(str[0]==':' && str[1]!=0) // this data file
505 {
506 size_t l=str.length()+1;
507 char *buf = new char[l]; memset(buf,0,l);
508 for(size_t i=1;str[i]!=0 && str[i]!=':' && i<l;i++) buf[i-1]=str[i];
509 HMDT res = new mglData(buf); delete []buf;
510 return res;
511 }
512 if(str[0]=='[') // this is manual subdata
513 {
514 long i, j, br=0,k;
515 bool ar=true,mt=false;
516 HMDT res=0;
517 for(i=1,j=1;i<len-1;i++)
518 {
519 if(str[i]=='[') br++;
520 if(str[i]==']' && br>0) br--;
521 if(str[i]==',' && !br)
522 {
523 HMDT a1=mglFormulaCalc(str.substr(j,i-j), arg, head);
524 if(j==1)
525 { res = a1; ar = (a1->nx==1); mt = (a1->nx>1 && a1->ny==1); }
526 else
527 {
528 if(ar) // res 1d array
529 { k = res->nx; res->Insert('x',k); mgl_data_put_dat(res,a1,k,-1,-1); }
530 else if(mt) // res 2d array
531 { k = res->ny; res->Insert('y',k); mgl_data_put_dat(res,a1,-1,k,-1); }
532 else // res 3d array
533 { k = res->nz; res->Insert('z',k); mgl_data_put_dat(res,a1,-1,-1,k); }
534 mgl_delete_data(a1);
535 }
536 j=i+1;
537 }
538 }
539 HMDT a1=mglFormulaCalc(str.substr(j,i-j), arg, head);
540 if(j==1)
541 { res = a1; ar = (a1->nx==1); mt = (a1->nx>1 && a1->ny==1); }
542 else
543 {
544 if(ar) // res 1d array
545 { k = res->nx; res->Insert('x',k); mgl_data_put_dat(res,a1,k,-1,-1); }
546 else if(mt) // res 2d array
547 { k = res->ny; res->Insert('y',k); mgl_data_put_dat(res,a1,-1,k,-1); }
548 else // res 3d array
549 { k = res->nz; res->Insert('z',k); mgl_data_put_dat(res,a1,-1,-1,k); }
550 mgl_delete_data(a1);
551 }
552 return res;
553 }
554
555 n=mglFindInText(str,"&|"); // lowest priority -- logical
556 if(n>=0)
557 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, str[n]=='|'?cor:cand,fns);
558 n=mglFindInText(str,"<>="); // low priority -- conditions
559 if(n>=0)
560 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, str[n]=='<'?clt:(str[n]=='>'?cgt:ceq),fns);
561 n=mglFindInText(str,"+-"); // normal priority -- additions
562 if(n>=0 && (n<2 || !strchr("eE",str[n-1]) || (str[n-2]!='.' && !isdigit(str[n-2])) ))
563 return str[n]=='+'? mglApplyOperAdd(str.substr(0,n),str.substr(n+1),arg, head,fns) :
564 mglApplyOperSub(str.substr(0,n),str.substr(n+1),arg, head,fns);
565 n=mglFindInText(str,"*/%"); // high priority -- multiplications
566 if(n>=0)
567 return str[n]=='*'? mglApplyOperMul(str.substr(0,n),str.substr(n+1),arg, head,fns) :
568 (str[n]=='/'? mglApplyOperDiv(str.substr(0,n),str.substr(n+1),arg, head,fns) :
569 mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_fmod,fns));
570 n=mglFindInText(str,"@"); // high priority -- combine
571 if(n>=0)
572 {
573 HMDT a1 = mglFormulaCalc(str.substr(0,n),arg, head);
574 HMDT a2 = mglFormulaCalc(str.substr(n+1),arg, head);
575 HMDT res = mgl_data_combine(a1,a2);
576 mgl_delete_data(a1); mgl_delete_data(a2);
577 return res?res:new mglData;
578 }
579 n=mglFindInText(str,"^"); // highest priority -- power
580 if(n>=0) return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, ipw,fns);
581 n=mglFindInText(str,":"); // highest priority -- array
582 if(n>=0 && str.compare(L":"))
583 {
584 HMDT a1=mglFormulaCalc(str.substr(0,n), arg, head);
585 HMDT a2=mglFormulaCalc(str.substr(n+1), arg, head);
586 HMDT res = new mglData(abs(int(a2->a[0]+0.5)-int(a1->a[0]+0.5))+1);
587 res->Fill(a1->a[0], a2->a[0]);
588 mgl_delete_data(a1); mgl_delete_data(a2);
589 return res;
590 }
591 n=mglFindInText(str,"."); // highest priority -- suffixes
592 wchar_t c0 = str[n+1];
593 if(n>=0 && c0>='a' && c0!='e' && c0!='E' && !isdigit(c0))
594 {
595 mreal x,y,z,k,v=NAN;
596 HMDT d = mglFormulaCalc(str.substr(0,n), arg, head);
597 long ns[3] = {d->nx-1, d->ny-1, d->nz-1};
598 const std::wstring &p=str.substr(n+1);
599 wchar_t ch = p[1];
600 if(c0=='a')
601 {
602 if(ch==0) v = d->a[0];
603 else
604 {
605 d->Momentum(ch,x,y);
606 if(ch=='a') v = x;
607 else if(ch>='x' && ch<='z') v = x/ns[ch-'x'];
608 }
609 }
610 else if(c0=='n')
611 {
612 if(ch>='x' && ch<='z') v = ns[p[1]-'x']+1;
613 else if(!p.compare(L"nmax")) { v=d->MaximalNeg(); }
614 else if(!p.compare(L"nmin")) { v=d->Minimal(); v = v<0?v:0; }
615 }
616 else if(c0=='k')
617 { d->Momentum(ch,x,y,z,k); v=k; }
618 else if(c0=='w')
619 {
620 d->Momentum(ch,x,y);
621 if(ch=='a') v = y;
622 else if(ch>='x' && ch<='z') v = y/ns[ch-'x'];
623 }
624 else if(c0=='m')
625 {
626 if(ch=='a' && p[2]=='x') v = d->Maximal();
627 else if(ch=='i' && p[2]=='n') v = d->Minimal();
628 else if(ch=='x' && p[2]=='f') v = d->Maximal('x',0)/mreal(ns[0]);
629 else if(ch=='x' && p[2]=='l') v = d->Maximal('x',-1)/mreal(ns[0]);
630 else if(ch=='x') { d->Maximal(x,y,z); v = x/ns[0]; }
631 else if(ch=='y' && p[2]=='f') v = d->Maximal('y',0)/mreal(ns[1]);
632 else if(ch=='y' && p[2]=='l') v = d->Maximal('y',-1)/mreal(ns[1]);
633 else if(ch=='y') { d->Maximal(x,y,z); v = y/ns[1]; }
634 else if(ch=='z' && p[2]=='f') v = d->Maximal('z',0)/mreal(ns[2]);
635 else if(ch=='z' && p[2]=='l') v = d->Maximal('z',-1)/mreal(ns[2]);
636 else if(ch=='z') { d->Maximal(x,y,z); v = z/ns[2]; }
637 }
638 else if(c0=='s')
639 {
640 if(ch=='u' && p[2]=='m') v = d->Momentum('x',x,y);
641 else if(ch=='a' || (ch>='x' && ch<='z'))
642 { d->Momentum(ch,x,y,z,k); v = z; }
643 }
644 else if(!p.compare(L"fst")) { long i=-1,j=-1,l=-1; v = d->Find(0,i,j,l); }
645 else if(!p.compare(L"lst")) { long i=-1,j=-1,l=-1; v = d->Last(0,i,j,l); }
646 else if(!p.compare(L"pmax")) { v=d->Maximal(); v = v>0?v:0; }
647 else if(!p.compare(L"pmin")) { v=d->MinimalPos(); }
648 delete d;
649 // if this is valid suffix when finish parsing (it can be mreal number)
650 if(mgl_isfin(v))
651 { HMDT res = new mglData; res->a[0]=v; return res; }
652 }
653 if(str[0]=='`')
654 {
655 HMDT res = mglFormulaCalc(str.substr(1), arg, head);
656 res->Transpose(); return res;
657 }
658 for(n=0;n<len;n++) if(str[n]=='(') break;
659 if(n>=len) // this is number or variable
660 {
661 HCDT v = (str!=L"#$mgl")?FindVar(head, str):0;
662 if(v) return new mglData(v);
663 const mglNum *f = arg?arg->FindNum(str.c_str()):0;
664 if(f) { HMDT res = new mglData; res->a[0] = f->d; return res; }
665 else if(!str.compare(L"rnd"))
666 {
667 v=FindVar(head, L"#$mgl");
668 HMDT res = v?new mglData(v->GetNx(),v->GetNy(),v->GetNz()) : new mglData;
669 for(long i=0;i<res->GetNN();i++) res->a[i] = mgl_rnd();
670 return res;
671 }
672 else
673 {
674 HMDT res = new mglData;
675 wchar_t ch = str[0];
676 if(ch<':') res->a[0] = wcstod(str.c_str(),0); // this is number
677 else if(!str.compare(L"pi")) res->a[0] = M_PI;
678 else if(ch==':') res->a[0] = -1;
679 else if(!str.compare(L"nan")) res->a[0] = NAN;
680 else if(!str.compare(L"inf")) res->a[0] = INFINITY;
681 return res;
682 }
683 }
684 else
685 {
686 std::wstring nm = str.substr(0,n);
687 str = str.substr(n+1,len-n-2); len -= n+2;
688 HCDT v = FindVar(head, nm);
689 HMDT tmp = 0;
690 // mglVar *v = arg->FindVar(nm.c_str());
691 if(!v && !nm.compare(0,7,L"jacobi_")) nm = nm.substr(7);
692 if(!v && nm.empty())
693 {
694 long m=mglFindInText(str,")");
695 if(m>1)
696 { nm = str.substr(0,m); str = str.substr(m+2);
697 tmp = mglFormulaCalc(nm,arg,head); v = tmp; }
698 }
699 n = mglFindInText(str,",");
700 if(v) // subdata
701 {
702 if(str[0]=='\'' && str[len-1]=='\'') // this is column call
703 {
704 char *buf = new char[len];
705 mgl_wcstombs(buf, str.substr(1).c_str(), len-1); buf[len-1]=0;
706 HMDT res = mgl_data_column(v,buf);
707 if(tmp) mgl_delete_data(tmp);
708 delete []buf; return res?res:new mglData;
709 }
710 else
711 {
712 HMDT a1=0, a2=0, a3=0;
713 if(n>0)
714 {
715 long m=mglFindInText(str.substr(0,n),",");
716 if(m>0)
717 {
718 str[m]=0;
719 a1 = mglFormulaCalc(str.substr(0,m), arg, head);
720 a2 = mglFormulaCalc(str.substr(m+1,n-m-1), arg, head);
721 a3 = mglFormulaCalc(str.substr(n+1), arg, head);
722 }
723 else
724 {
725 a1 = mglFormulaCalc(str.substr(0,n), arg, head);
726 a2 = mglFormulaCalc(str.substr(n+1), arg, head);
727 }
728 }
729 else a1 = mglFormulaCalc(str, arg, head);
730 HMDT res = mgl_data_subdata_ext(v,a1,a2,a3);
731 if(tmp) mgl_delete_data(tmp);
732 mgl_delete_data(a1); mgl_delete_data(a2);
733 mgl_delete_data(a3); return res?res:new mglData;
734 }
735 if(tmp) mgl_delete_data(tmp);
736 }
737 else if(nm[0]=='a') // function
738 {
739 if(!nm.compare(L"asin")) return mglApplyFunc(str, arg, head, asin,fns);
740 else if(!nm.compare(L"acos")) return mglApplyFunc(str, arg, head, acos,fns);
741 else if(!nm.compare(L"atan")) return mglApplyFunc(str, arg, head, atan,fns);
742 else if(!nm.compare(L"asinh")) return mglApplyFunc(str, arg, head, mgl_asinh,fns);
743 else if(!nm.compare(L"acosh")) return mglApplyFunc(str, arg, head, mgl_acosh,fns);
744 else if(!nm.compare(L"atanh")) return mglApplyFunc(str, arg, head, mgl_atanh,fns);
745 else if(!nm.compare(L"arg"))
746 {
747 if(n>0) return mglApplyOper(str.substr(n+1),str.substr(0,n),arg, head, atan2,fns);
748 else
749 {
750 HADT a1 = mglFormulaCalcC(str, arg, head);
751 HMDT res = mgl_datac_arg(a1);
752 mgl_delete_datac(a1); return res;
753 }
754 }
755 else if(!nm.compare(L"abs"))
756 {
757 if(n>0) return mglApplyOper(str.substr(n+1),str.substr(0,n),arg, head, hypot,fns);
758 else
759 {
760 HADT a1 = mglFormulaCalcC(str, arg, head);
761 HMDT res = mgl_datac_abs(a1);
762 mgl_delete_datac(a1); return res;
763 }
764 }
765 #if MGL_HAVE_GSL
766 else if(!nm.compare(L"ai") || !nm.compare(L"airy_ai"))
767 return mglApplyFuncGSL(str, arg, head, gsl_sf_airy_Ai,fns);
768 else if(!nm.compare(L"airy_dai"))
769 return mglApplyFuncGSL(str, arg, head, gsl_sf_airy_Ai_deriv,fns);
770 else if(!nm.compare(L"airy_bi"))
771 return mglApplyFuncGSL(str, arg, head, gsl_sf_airy_Bi,fns);
772 else if(!nm.compare(L"airy_dbi"))
773 return mglApplyFuncGSL(str, arg, head, gsl_sf_airy_Bi_deriv,fns);
774 }
775 else if(nm[0]=='b')
776 {
777 if(!nm.compare(L"beta") && n>0)
778 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_beta,fns);
779 else if(!nm.compare(L"bi"))
780 return mglApplyFuncGSL(str, arg, head, gsl_sf_airy_Bi,fns);
781 else if(!nm.compare(L"bessel_i") && n>0)
782 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Inu,fns);
783 else if(!nm.compare(L"bessel_j") && n>0)
784 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Jnu,fns);
785 else if(!nm.compare(L"bessel_k") && n>0)
786 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Knu,fns);
787 else if(!nm.compare(L"bessel_y") && n>0)
788 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Ynu,fns);
789 #endif
790 }
791 else if(nm[0]=='c')
792 {
793 if(!nm.compare(L"cos")) return mglApplyFunc(str, arg, head, cos,fns);
794 else if(!nm.compare(L"cosh") || !nm.compare(L"ch")) return mglApplyFunc(str, arg, head, cosh,fns);
795 else if(!nm.compare(L"conj"))
796 {
797 HADT a1 = mglFormulaCalcC(str, arg, head);
798 HMDT res = mgl_datac_real(a1);
799 mgl_delete_datac(a1); return res;
800 }
801 #if MGL_HAVE_GSL
802 else if(!nm.compare(L"ci")) return mglApplyFunc(str, arg, head, gsl_sf_Ci,fns);
803 else if(!nm.compare(L"cn") && n>0)
804 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_cn,fns);
805 else if(!nm.compare(L"cs") && n>0)
806 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_cs,fns);
807 else if(!nm.compare(L"cd") && n>0)
808 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_cd,fns);
809 #endif
810 }
811 else if(nm[0]=='e')
812 {
813 if(!nm.compare(L"exp")) return mglApplyFunc(str, arg, head, exp,fns);
814 #if MGL_HAVE_GSL
815 else if(!nm.compare(L"erf")) return mglApplyFunc(str, arg, head, gsl_sf_erf,fns);
816 else if(!nm.compare(L"ee") || !nm.compare(L"elliptic_ec"))
817 return mglApplyFuncGSL(str, arg, head, gsl_sf_ellint_Ecomp,fns);
818 else if(!nm.compare(L"ek") || !nm.compare(L"elliptic_kc"))
819 return mglApplyFuncGSL(str, arg, head, gsl_sf_ellint_Kcomp,fns);
820 else if((!nm.compare(L"e") || !nm.compare(L"elliptic_e")) && n>0)
821 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gslEllE,fns);
822 else if(!nm.compare(L"elliptic_f"))
823 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gslEllF,fns);
824
825 else if(!nm.compare(L"ei")) return mglApplyFunc(str, arg, head, gsl_sf_expint_Ei,fns);
826 else if(!nm.compare(L"e1")) return mglApplyFunc(str, arg, head, gsl_sf_expint_E1,fns);
827 else if(!nm.compare(L"e2")) return mglApplyFunc(str, arg, head, gsl_sf_expint_E2,fns);
828 else if(!nm.compare(L"eta")) return mglApplyFunc(str, arg, head, gsl_sf_eta,fns);
829 else if(!nm.compare(L"ei3")) return mglApplyFunc(str, arg, head, gsl_sf_expint_3,fns);
830 #endif
831 }
832 else if(nm[0]=='l')
833 {
834 if(!nm.compare(L"log") && n>0)
835 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, llg,fns);
836 else if(!nm.compare(L"lg")) return mglApplyFunc(str, arg, head, log10,fns);
837 else if(!nm.compare(L"ln")) return mglApplyFunc(str, arg, head, log,fns);
838 #if MGL_HAVE_GSL
839 else if(!nm.compare(L"li2")) return mglApplyFunc(str, arg, head, gsl_sf_dilog,fns);
840 else if(!nm.compare(L"legendre") && n>0)
841 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gslLegP,fns);
842 #endif
843 }
844 else if(nm[0]=='s')
845 {
846 if(!nm.compare(L"sqrt")) return mglApplyFunc(str, arg, head, sqrt,fns);
847 else if(!nm.compare(L"sin")) return mglApplyFunc(str, arg, head, sin,fns);
848 else if(!nm.compare(L"step")) return mglApplyFunc(str, arg, head, stp,fns);
849 else if(!nm.compare(L"sign")) return mglApplyFunc(str, arg, head, sgn,fns);
850 else if(!nm.compare(L"sinh") || !nm.compare(L"sh")) return mglApplyFunc(str, arg, head, sinh,fns);
851 #if MGL_HAVE_GSL
852 else if(!nm.compare(L"si")) return mglApplyFunc(str, arg, head, gsl_sf_Si,fns);
853 else if(!nm.compare(L"sinc")) return mglApplyFunc(str, arg, head, gsl_sf_sinc,fns);
854 else if(!nm.compare(L"sn") && n>0)
855 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_sn,fns);
856 else if(!nm.compare(L"sc") && n>0)
857 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_sc,fns);
858 else if(!nm.compare(L"sd") && n>0)
859 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_sd,fns);
860 #endif
861 else if(!nm.compare(L"sum") && n>0)
862 {
863 HMDT a=NULL, b=mglFormulaCalcA(str.substr(n+1), arg, head, fns);
864 long m = long(b->a[0]+0.5);
865 const char *s = head.size()>0?head[head.size()-1]->s.s:"";
866 if(m>0)
867 {
868 std::vector<mglDataA*> hh(head);
869 int in=0;
870 if(s[0]=='_' && s[1]>='i' && s[1]<'z') in = s[1]-'i'+1;
871 char name[3] = {'_',char('i'+in),0};
872 b->s = name; b->Create(1,1,1); hh.push_back(b);
873 a = mglFormulaCalcA(str.substr(0,n), arg, hh, fns);
874 long nn = a->GetNN();
875 for(long i=1;i<m;i++)
876 {
877 b->a[0] = i;
878 HMDT c = mglFormulaCalcA(str.substr(0,n), arg, hh, fns);
879 for(long j=0;j<nn;j++) a->a[j] += c->a[j];
880 mgl_delete_data(c);
881 }
882 mgl_delete_data(b);
883 }
884 else
885 { a = new mglData; a->a[0]=NAN; }
886 return a;
887 }
888 else if(!nm.compare(L"spline"))
889 {
890 HMDT a = mglFormulaCalcA(str.substr(0,n), arg, head, fns);
891 long n1 = mglFindInText(str,",",1), n2 = mglFindInText(str,",",2);
892 HMDT ix, iy=NULL, iz=NULL;
893 if(n1>0)
894 {
895 ix = mglFormulaCalcA(str.substr(n+1,n1), arg, head, fns);
896 if(n2>0)
897 {
898 iy = mglFormulaCalcA(str.substr(n1+1,n2), arg, head, fns);
899 iz = mglFormulaCalcA(str.substr(n2+1), arg, head, fns);
900 }
901 else iy = mglFormulaCalcA(str.substr(n1+1), arg, head, fns);
902
903 }
904 else ix = mglFormulaCalcA(str.substr(n+1), arg, head, fns);
905 HMDT res = mgl_data_evaluate(a, ix,iy,iz, false);
906 if(a) mgl_delete_data(a);
907 if(ix) mgl_delete_data(ix);
908 if(iy) mgl_delete_data(iy);
909 if(iz) mgl_delete_data(iz);
910 return res;
911 }
912 }
913 else if(nm[0]=='t')
914 {
915 if(!nm.compare(L"tg") || !nm.compare(L"tan"))
916 return mglApplyFunc(str, arg, head, tan,fns);
917 else if(!nm.compare(L"tanh") || !nm.compare(L"th"))
918 return mglApplyFunc(str, arg, head, tanh,fns);
919 }
920 else if(!nm.compare(L"pow") && n>0)
921 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, pow,fns);
922 else if(nm[0]=='m')
923 {
924 if(!nm.compare(L"mod") && n>0)
925 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_fmod,fns);
926 else if(!nm.compare(L"min") && n>0)
927 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_fmin,fns);
928 else if(!nm.compare(L"max") && n>0)
929 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_fmax,fns);
930 }
931 else if(!nm.compare(L"int")) return mglApplyFunc(str, arg, head, floor,fns);
932 else if(!nm.compare(L"random"))
933 { HMDT res=mglFormulaCalc(str, arg, head); mreal *a = res->a;
934 for(long i=0;i<res->GetNN();i++) a[i] = mgl_rnd();
935 return res; }
936 else if(!nm.compare(L"real"))
937 {
938 HADT a1 = mglFormulaCalcC(str, arg, head);
939 HMDT res = mgl_datac_real(a1);
940 mgl_delete_datac(a1); return res;
941 }
942 else if(!nm.compare(L"imag"))
943 {
944 HADT a1 = mglFormulaCalcC(str, arg, head);
945 HMDT res = mgl_datac_imag(a1);
946 mgl_delete_datac(a1); return res;
947 }
948 else if(!nm.compare(L"norm"))
949 {
950 HADT a1 = mglFormulaCalcC(str, arg, head);
951 HMDT res = mgl_datac_norm(a1);
952 mgl_delete_datac(a1); return res;
953 }
954 #if MGL_HAVE_GSL
955 else if(!nm.compare(L"dn") && n>0)
956 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_dn,fns);
957 else if(!nm.compare(L"ds") && n>0)
958 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_ds,fns);
959 else if(!nm.compare(L"dc") && n>0)
960 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_dc,fns);
961 else if(!nm.compare(L"nc") && n>0)
962 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_nc,fns);
963 else if(!nm.compare(L"ns") && n>0)
964 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_ns,fns);
965 else if(!nm.compare(L"nd") && n>0)
966 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, mgl_jac_nd,fns);
967 else if(!nm.compare(L"i") && n>0)
968 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Inu,fns);
969 else if(!nm.compare(L"j") && n>0)
970 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Jnu,fns);
971 else if(!nm.compare(L"k") && n>0)
972 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Knu,fns);
973 else if(!nm.compare(L"y") && n>0)
974 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_bessel_Ynu,fns);
975 else if(!nm.compare(L"f") && n>0)
976 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gslEllF,fns);
977 else if(!nm.compare(L"hypot") && n>0)
978 return mglApplyOper(str.substr(n+1),str.substr(0,n),arg, head, hypot,fns);
979 else if(!nm.compare(L"gamma")) return mglApplyFunc(str, arg, head, gsl_sf_gamma,fns);
980 else if(!nm.compare(L"gamma_inc") && n>0)
981 return mglApplyOper(str.substr(0,n),str.substr(n+1),arg, head, gsl_sf_gamma_inc,fns);
982 else if(!nm.compare(L"w0")) return mglApplyFunc(str, arg, head, gsl_sf_lambert_W0,fns);
983 else if(!nm.compare(L"w1")) return mglApplyFunc(str, arg, head, gsl_sf_lambert_Wm1,fns);
984 else if(!nm.compare(L"psi")) return mglApplyFunc(str, arg, head, gsl_sf_psi,fns);
985 else if(!nm.compare(L"zeta")) return mglApplyFunc(str, arg, head, gsl_sf_zeta,fns);
986 else if(!nm.compare(L"z")) return mglApplyFunc(str, arg, head, gsl_sf_dawson,fns);
987 #endif
988 else if(!nm.compare(L"value"))
989 {
990 HMDT a = mglFormulaCalcA(str.substr(0,n), arg, head, fns);
991 long n1 = mglFindInText(str,",",1), n2 = mglFindInText(str,",",2);
992 HMDT ix, iy=NULL, iz=NULL;
993 if(n1>0)
994 {
995 ix = mglFormulaCalcA(str.substr(n+1,n1), arg, head, fns);
996 if(n2>0)
997 {
998 iy = mglFormulaCalcA(str.substr(n1+1,n2), arg, head, fns);
999 iz = mglFormulaCalcA(str.substr(n2+1), arg, head, fns);
1000 }
1001 else iy = mglFormulaCalcA(str.substr(n1+1), arg, head, fns);
1002
1003 }
1004 else ix = mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1005 HMDT res = mgl_data_subdata_ext(a, ix,iy,iz);
1006 if(a) mgl_delete_data(a);
1007 if(ix) mgl_delete_data(ix);
1008 if(iy) mgl_delete_data(iy);
1009 if(iz) mgl_delete_data(iz);
1010 return res;
1011 }
1012 else if(!nm.compare(L"dsum") && n>0)
1013 {
1014 HMDT a=NULL, b=mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1015 long m = long(b->a[0]+0.5);
1016 const char *s = head.size()>0?head[head.size()-1]->s.s:"";
1017 if(m>0)
1018 {
1019 std::vector<mglDataA*> hh(head);
1020 int in=0, zn=1;
1021 if(s[0]=='_' && s[1]>='i' && s[1]<'z') in = s[1]-'i'+1;
1022 char name[3] = {'_',char('i'+in),0};
1023 b->s = name; b->Create(1,1,1); hh.push_back(b);
1024 a = mglFormulaCalcA(str.substr(0,n), arg, hh, fns);
1025 long nn = a->GetNN();
1026 for(long i=1;i<m;i++)
1027 {
1028 b->a[0] = i; zn *= -1;
1029 HMDT c = mglFormulaCalcA(str.substr(0,n), arg, hh, fns);
1030 for(long j=0;j<nn;j++) a->a[j] += zn*c->a[j];
1031 mgl_delete_data(c);
1032 }
1033 mgl_delete_data(b);
1034 }
1035 else
1036 { a = new mglData; a->a[0]=NAN; }
1037 return a;
1038 }
1039 else if(!nm.compare(L"prod") && n>0)
1040 {
1041 HMDT a=NULL, b=mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1042 long m = long(b->a[0]+0.5);
1043 const char *s = head.size()>0?head[head.size()-1]->s.s:"";
1044 if(m>0)
1045 {
1046 std::vector<mglDataA*> hh(head);
1047 int in=0;
1048 if(s[0]=='_' && s[1]>='i' && s[1]<'z') in = s[1]-'i'+1;
1049 char name[3] = {'_',char('i'+in),0};
1050 b->s = name; b->Create(1,1,1); hh.push_back(b);
1051 a = mglFormulaCalcA(str.substr(0,n), arg, hh, fns);
1052 long nn = a->GetNN();
1053 for(long i=1;i<m;i++)
1054 {
1055 b->a[0] = i;
1056 HMDT c = mglFormulaCalcA(str.substr(0,n), arg, hh, fns);
1057 for(long j=0;j<nn;j++) a->a[j] *= c->a[j];
1058 mgl_delete_data(c);
1059 }
1060 mgl_delete_data(b);
1061 }
1062 else
1063 { a = new mglData; a->a[0]=NAN; }
1064 return a;
1065 }
1066 else if(nm[0]=='f' && nm[1]=='n' && fns.size()>1)
1067 {
1068 HMDT a=NULL, b=NULL, r=NULL;
1069 std::vector<mglDataA*> hh(head);
1070 std::vector<std::wstring> tmp; // disable recursion
1071 if(n>0)
1072 {
1073 a = mglFormulaCalcA(str.substr(0,n), arg, head, fns);
1074 if(a) { a->s = "_1"; hh.push_back(a); }
1075 b = mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1076 if(b) { b->s = "_2"; hh.push_back(b); }
1077 }
1078 else
1079 {
1080 a = mglFormulaCalcA(str, arg, head, fns);
1081 if(a) { a->s = "_1"; hh.push_back(a); }
1082 }
1083 if(!nm.compare(L"fn1") && fns.size()>1) r = mglFormulaCalcA(fns[1], arg, hh, tmp);
1084 if(!nm.compare(L"fn2") && fns.size()>2) r = mglFormulaCalcA(fns[2], arg, hh, tmp);
1085 if(!nm.compare(L"fn3") && fns.size()>3) r = mglFormulaCalcA(fns[3], arg, hh, tmp);
1086 if(!nm.compare(L"fn4") && fns.size()>4) r = mglFormulaCalcA(fns[4], arg, hh, tmp);
1087 if(!nm.compare(L"fn5") && fns.size()>5) r = mglFormulaCalcA(fns[5], arg, hh, tmp);
1088 if(!nm.compare(L"fn6") && fns.size()>6) r = mglFormulaCalcA(fns[6], arg, hh, tmp);
1089 if(!nm.compare(L"fn7") && fns.size()>7) r = mglFormulaCalcA(fns[7], arg, hh, tmp);
1090 if(!nm.compare(L"fn8") && fns.size()>8) r = mglFormulaCalcA(fns[8], arg, hh, tmp);
1091 if(!nm.compare(L"fn9") && fns.size()>9) r = mglFormulaCalcA(fns[9], arg, hh, tmp);
1092 if(a) mgl_delete_data(a);
1093 if(b) mgl_delete_data(b);
1094 if(!r) { r = new mglData; r->a[0]=NAN; }
1095 return r;
1096 }
1097 }
1098 HMDT res = new mglData; res->a[0]=NAN; return res;
1099 }
1100 //-----------------------------------------------------------------------------
1101 dual MGL_LOCAL_CONST ceqc(dual a,dual b); //{return a==b?1:0;}
1102 dual MGL_LOCAL_CONST cltc(dual a,dual b); //{return real(a-b)<0?1:0;}
1103 dual MGL_LOCAL_CONST cgtc(dual a,dual b); //{return real(a-b)>0?1:0;}
1104 dual MGL_LOCAL_CONST ipwc(dual a,dual b); //{return mgl_ipowc(a,int(b.real()));}
1105 dual MGL_LOCAL_CONST powc(dual a,dual b); //{return exp(b*log(a)); }
1106 dual MGL_LOCAL_CONST llgc(dual a,dual b); //{return log(a)/log(b); }
1107 dual MGL_LOCAL_CONST cmplxc(dual a,dual b); //{return a+dual(0,1)*b; }
1108 dual MGL_LOCAL_CONST expi(dual a); //{ return exp(dual(0,1)*a); }
1109 dual MGL_LOCAL_CONST expi(double a); //{ return dual(cos(a),sin(a)); }
1110 //-----------------------------------------------------------------------------
1111 dual MGL_LOCAL_CONST hypotc(dual x, dual y); //{ return sqrt(x*x+y*y); }
1112 dual MGL_LOCAL_CONST asinhc(dual x); //{ return log(x+sqrt(x*x+mreal(1))); }
1113 dual MGL_LOCAL_CONST acoshc(dual x); //{ return log(x+sqrt(x*x-mreal(1))); }
1114 dual MGL_LOCAL_CONST atanhc(dual x); //{ return log((mreal(1)+x)/(mreal(1)-x))/mreal(2); }
1115 dual MGL_LOCAL_CONST conjc(dual x); //{ return dual(real(x),-imag(x)); }
1116 dual MGL_LOCAL_CONST sinc(dual x); //{ return sin(x); }
1117 dual MGL_LOCAL_CONST cosc(dual x); //{ return cos(x); }
1118 dual MGL_LOCAL_CONST tanc(dual x); //{ return tan(x); }
1119 dual MGL_LOCAL_CONST sinhc(dual x); //{ return sinh(x); }
1120 dual MGL_LOCAL_CONST coshc(dual x); //{ return cosh(x); }
1121 dual MGL_LOCAL_CONST tanhc(dual x); //{ return tanh(x); }
1122 dual MGL_LOCAL_CONST asinc(dual x); //{ return log(ic*x+sqrt(mreal(1)-x*x))/ic; }
1123 dual MGL_LOCAL_CONST acosc(dual x); //{ return log(x+sqrt(x*x-mreal(1)))/ic; }
1124 dual MGL_LOCAL_CONST atanc(dual x); //{ return log((ic-x)/(ic+x))/(mreal(2)*ic); }
1125 dual MGL_LOCAL_CONST expc(dual x); //{ return exp(x); }
1126 dual MGL_LOCAL_CONST sqrtc(dual x); //{ return sqrt(x); }
1127 dual MGL_LOCAL_CONST logc(dual x); //{ return log(x); }
1128 dual MGL_LOCAL_CONST absc(dual x); //{ return abs(x); }
1129 dual MGL_LOCAL_CONST argc(dual x); //{ return arg(x); }
1130 dual MGL_LOCAL_CONST lgc(dual x); //{ return log10(x);}
1131 dual MGL_LOCAL_CONST realc(dual x); //{ return real(x); }
1132 dual MGL_LOCAL_CONST imagc(dual x); //{ return imag(x); }
1133 dual MGL_LOCAL_CONST normc(dual x); //{ return norm(x); }
1134 //-----------------------------------------------------------------------------
1135 /// Parse string and substitute the script argument
1136 // All numbers are presented as mglData(1). Do boundary checking.
1137 // NOTE: In any case where number is required the mglData::a[0] is used.
1138 // String flag is binary 0x1 -> 'x', 0x2 -> 'y', 0x4 -> 'z'
mglFormulaCalcC(const std::wstring & str,mglParser * arg,const std::vector<mglDataA * > & head)1139 HADT MGL_NO_EXPORT mglFormulaCalcC(const std::wstring &str, mglParser *arg, const std::vector<mglDataA*> &head)
1140 {
1141 if(str.empty()) return NULL;
1142 std::vector<std::wstring> fns = mgl_wcs_args(str,'\\');
1143 return mglFormulaCalcAC(fns[0],arg,head,fns);
1144 }
mglFormulaCalcAC(std::wstring str,mglParser * arg,const std::vector<mglDataA * > & head,const std::vector<std::wstring> & fns)1145 HADT MGL_NO_EXPORT mglFormulaCalcAC(std::wstring str, mglParser *arg, const std::vector<mglDataA*> &head, const std::vector<std::wstring> &fns)
1146 {
1147 #if MGL_HAVE_GSL
1148 gsl_set_error_handler_off();
1149 #endif
1150 if(str.empty()) return new mglDataC; // nothing to parse
1151 str = mgl_trim_ws(str);
1152 mreal tval = mgl_gettime(str);
1153 if(mgl_isnum(tval))
1154 { mglDataC *r=new mglDataC; r->a[0] = tval; return r; }
1155
1156 long n,len=str.length();
1157 if(str[0]=='(' && mglCheck(str.substr(1,len-2))) // remove braces
1158 { str = str.substr(1,len-2); len-=2; }
1159 if(str[0]=='[') // this is manual subdata
1160 {
1161 long i, j, br=0,k;
1162 bool ar=true,mt=false;
1163 HADT res=0;
1164 for(i=1,j=1;i<len-1;i++)
1165 {
1166 if(str[i]=='[') br++;
1167 if(str[i]==']' && br>0) br--;
1168 if(str[i]==',' && !br)
1169 {
1170 HADT a1=mglFormulaCalcC(str.substr(j,i-j), arg, head);
1171 if(j==1)
1172 { res = a1; ar = (a1->nx==1); mt = (a1->nx>1 && a1->ny==1); }
1173 else
1174 {
1175 if(ar) // res 1d array
1176 { k = res->nx; res->Insert('x',k); mgl_datac_put_dat(res,a1,k,-1,-1); }
1177 else if(mt) // res 2d array
1178 { k = res->ny; res->Insert('y',k); mgl_datac_put_dat(res,a1,-1,k,-1); }
1179 else // res 3d array
1180 { k = res->nz; res->Insert('z',k); mgl_datac_put_dat(res,a1,-1,-1,k); }
1181 mgl_delete_datac(a1);
1182 }
1183 j=i+1;
1184 }
1185 }
1186 HADT a1=mglFormulaCalcC(str.substr(j,i-j), arg, head);
1187 if(j==1)
1188 { res = a1; ar = (a1->nx==1); mt = (a1->nx>1 && a1->ny==1); }
1189 else
1190 {
1191 if(ar) // res 1d array
1192 { k = res->nx; res->Insert('x',k); mgl_datac_put_dat(res,a1,k,-1,-1); }
1193 else if(mt) // res 2d array
1194 { k = res->ny; res->Insert('y',k); mgl_datac_put_dat(res,a1,-1,k,-1); }
1195 else // res 3d array
1196 { k = res->nz; res->Insert('z',k); mgl_datac_put_dat(res,a1,-1,-1,k); }
1197 mgl_delete_datac(a1);
1198 }
1199 return res;
1200 }
1201
1202 n=mglFindInText(str,"<>="); // low priority -- conditions
1203 if(n>=0)
1204 return mglApplyOperC(str.substr(0,n),str.substr(n+1),arg, head, str[n]=='<'?cltc:(str[n]=='>'?cgtc:ceqc),fns);
1205 n=mglFindInText(str,"+-"); // normal priority -- additions
1206 if(n>=0 && (n<2 || !strchr("eE",str[n-1]) || (str[n-2]!='.' && !isdigit(str[n-2]))))
1207 return str[n]=='+'? mglApplyOperAddC(str.substr(0,n),str.substr(n+1),arg, head,fns) : mglApplyOperSubC(str.substr(0,n),str.substr(n+1),arg, head,fns);
1208 n=mglFindInText(str,"*/"); // high priority -- multiplications
1209 if(n>=0)
1210 return str[n]=='*'? mglApplyOperMulC(str.substr(0,n),str.substr(n+1),arg, head,fns) : mglApplyOperDivC(str.substr(0,n),str.substr(n+1),arg, head,fns);
1211 n=mglFindInText(str,"@"); // high priority -- combine
1212 if(n>=0)
1213 {
1214 HADT a1 = mglFormulaCalcC(str.substr(0,n),arg, head);
1215 HADT a2 = mglFormulaCalcC(str.substr(n+1),arg, head);
1216 HADT res = mgl_datac_combine(a1,a2);
1217 mgl_delete_datac(a1); mgl_delete_datac(a2);
1218 return res?res:new mglDataC;
1219 }
1220 n=mglFindInText(str,"^"); // highest priority -- power
1221 if(n>=0)
1222 return mglApplyOperC(str.substr(0,n),str.substr(n+1),arg, head, ipwc,fns);
1223 n=mglFindInText(str,":"); // highest priority -- array
1224 if(n>=0 && str.compare(L":"))
1225 {
1226 HMDT a1=mglFormulaCalc(str.substr(0,n), arg, head);
1227 HMDT a2=mglFormulaCalc(str.substr(n+1), arg, head);
1228 HADT res = new mglDataC(abs(int(a2->a[0]+0.5)-int(a1->a[0]+0.5))+1);
1229 res->Fill(a1->a[0], a2->a[0]);
1230 mgl_delete_data(a1); mgl_delete_data(a2);
1231 return res;
1232 }
1233 n=mglFindInText(str,"."); // highest priority -- suffixes
1234 wchar_t c0 = str[n+1];
1235 if(n>=0 && c0>='a' && c0!='e' && c0!='E' && !isdigit(c0))
1236 {
1237 dual v=NAN;
1238 HADT d = mglFormulaCalcC(str.substr(0,n), arg, head);
1239 long ns[3] = {d->nx-1, d->ny-1, d->nz-1};
1240 const std::wstring &p=str.substr(n+1);
1241 wchar_t ch = p[1];
1242 if(c0=='a')
1243 {
1244 if(ch==0) v = d->a[0];
1245 else
1246 {
1247 mreal x,y;
1248 d->Momentum(ch,x,y);
1249 if(ch=='a') v = x;
1250 else if(ch>='x' && ch<='z') v = x/ns[ch-'x'];
1251 }
1252 }
1253 else if(c0=='n' && ch>='x' && ch<='z') v = ns[ch-'x']+1;
1254 else if(c0=='k')
1255 {
1256 mreal x,y,z,k;
1257 d->Momentum(ch,x,y,z,k);
1258 if(ch=='a') v = k;
1259 else if(ch>='x' && ch<='z') v = k/ns[ch-'x'];
1260 }
1261 else if(c0=='w')
1262 {
1263 mreal x,y;
1264 d->Momentum(ch,x,y);
1265 if(ch=='a') v = y;
1266 else if(ch>='x' && ch<='z') v = y/ns[ch-'x'];
1267 }
1268 else if(c0=='m')
1269 {
1270 mreal x,y,z;
1271 if(ch=='a' && p[2]=='x') v = d->Maximal();
1272 else if(ch=='i' && p[2]=='n') v = d->Minimal();
1273 else if(ch=='x' && p[2]=='f') v = d->Maximal('x',0)/mreal(ns[0]);
1274 else if(ch=='x' && p[2]=='l') v = d->Maximal('x',-1)/mreal(ns[0]);
1275 else if(ch=='x') { d->Maximal(x,y,z); v = x/ns[0]; }
1276 else if(ch=='y' && p[2]=='f') v = d->Maximal('y',0)/mreal(ns[1]);
1277 else if(ch=='y' && p[2]=='l') v = d->Maximal('y',-1)/mreal(ns[1]);
1278 else if(ch=='y') { d->Maximal(x,y,z); v = y/ns[1]; }
1279 else if(ch=='z' && p[2]=='f') v = d->Maximal('z',0)/mreal(ns[2]);
1280 else if(ch=='z' && p[2]=='l') v = d->Maximal('z',-1)/mreal(ns[2]);
1281 else if(ch=='z') { d->Maximal(x,y,z); v = z/ns[2]; }
1282 }
1283 else if(c0=='s')
1284 {
1285 mreal x,y,z,k;
1286 if(ch=='u' && p[2]=='m') v = d->Momentum('x',x,y);
1287 else if(ch=='a')
1288 { d->Momentum(ch,x,y,z,k); v = z; }
1289 else if(ch>='x' && ch<='z')
1290 { d->Momentum(ch,x,y,z,k); v = z/ns[ch-'x']; }
1291 }
1292 else if(!p.compare(L"fst")) { long i=-1,j=-1,l=-1; v = d->Find(0,i,j,l); }
1293 else if(!p.compare(L"lst")) { long i=-1,j=-1,l=-1; v = d->Last(0,i,j,l); }
1294 delete d;
1295 // if this is valid suffix when finish parsing (it can be mreal number)
1296 if(mgl_isfin(v))
1297 { HADT res = new mglDataC; res->a[0]=v; return res; }
1298 }
1299 if(str[0]=='`')
1300 {
1301 HADT res = mglFormulaCalcC(str.substr(1), arg, head);
1302 res->Transpose(); return res;
1303 }
1304 for(n=0;n<len;n++) if(str[n]=='(') break;
1305 if(n>=len) // this is number or variable
1306 {
1307 HCDT v = (str!=L"#$mgl")?FindVar(head, str):0;
1308 if(v) return new mglDataC(v);
1309 const mglNum *f = arg?arg->FindNum(str.c_str()):0;
1310 if(f) { HADT res = new mglDataC; res->a[0] = f->c; return res; }
1311 else if(!str.compare(L"rnd"))
1312 {
1313 v=FindVar(head, L"#$mgl");
1314 HADT res = v?new mglDataC(v->GetNx(),v->GetNy(),v->GetNz()) : new mglDataC;
1315 for(long i=0;i<res->GetNN();i++)
1316 res->a[i] = dual(mgl_rnd(), mgl_rnd());
1317 return res;
1318 }
1319 else
1320 {
1321 HADT res = new mglDataC;
1322 wchar_t ch = str[0];
1323 if(ch<':') // this is real number
1324 res->a[0] = (str[str.length()-1]=='i') ? dual(0,wcstod(str.c_str(),0)) : mreal(wcstod(str.c_str(),0));
1325 else if(ch=='i') // this is imaginary number
1326 res->a[0] = dual(0,(str[1]>='0' && str[1]<='9')?wcstod(str.c_str()+1,0):1);
1327 else if(!str.compare(L"pi")) res->a[0] = M_PI;
1328 else if(ch==':') res->a[0] = -1;
1329 else if(!str.compare(L"nan")) res->a[0] = NAN;
1330 else if(!str.compare(L"inf")) res->a[0] = INFINITY;
1331 return res;
1332 }
1333 }
1334 else
1335 {
1336 std::wstring nm = str.substr(0,n);
1337 str = str.substr(n+1,len-n-2); len -= n+2;
1338 HCDT v = FindVar(head, nm);
1339 HADT tmp = 0;
1340 // mglVar *v = arg->FindVar(nm.c_str());
1341 if(!v && !nm.compare(0,7,L"jacobi_")) nm = nm.substr(7);
1342 if(!v && nm.empty())
1343 {
1344 long m=mglFindInText(nm,")");
1345 if(m>1)
1346 { nm = str.substr(0,m); str = str.substr(m+2);
1347 tmp = mglFormulaCalcC(nm,arg,head); v = tmp; }
1348 }
1349 n = mglFindInText(str,",");
1350 if(v) // subdata
1351 {
1352 if(str[0]=='\'' && str[len-1]=='\'') // this is column call
1353 {
1354 char *buf = new char[len];
1355 mgl_wcstombs(buf, str.substr(1).c_str(), len-1); buf[len-1]=0;
1356 HADT res = mgl_datac_column(v,buf);
1357 if(tmp) mgl_delete_datac(tmp);
1358 delete []buf; return res?res:new mglDataC;
1359 }
1360 else
1361 {
1362 HMDT a1=0, a2=0, a3=0;
1363 if(n>0)
1364 {
1365 long m=mglFindInText(str.substr(0,n),",");
1366 if(m>0)
1367 {
1368 str[m]=0;
1369 a1 = mglFormulaCalc(str.substr(0,m), arg, head);
1370 a2 = mglFormulaCalc(str.substr(m+1,n-m-1), arg, head);
1371 a3 = mglFormulaCalc(str.substr(n+1), arg, head);
1372 }
1373 else
1374 {
1375 a1 = mglFormulaCalc(str.substr(0,n), arg, head);
1376 a2 = mglFormulaCalc(str.substr(n+1), arg, head);
1377 }
1378 }
1379 else a1 = mglFormulaCalc(str, arg, head);
1380 HADT res = mgl_datac_subdata_ext(v,a1,a2,a3);
1381 if(tmp) mgl_delete_datac(tmp);
1382 mgl_delete_data(a1); mgl_delete_data(a2);
1383 mgl_delete_data(a3); return res?res:new mglDataC;
1384 }
1385 if(tmp) mgl_delete_datac(tmp);
1386 }
1387 else if(nm[0]=='a') // function
1388 {
1389 if(!nm.compare(L"asin")) return mglApplyFuncC(str, arg, head, asinc,fns);
1390 else if(!nm.compare(L"acos")) return mglApplyFuncC(str, arg, head, acosc,fns);
1391 else if(!nm.compare(L"atan")) return mglApplyFuncC(str, arg, head, atanc,fns);
1392 else if(!nm.compare(L"asinh")) return mglApplyFuncC(str, arg, head, asinhc,fns);
1393 else if(!nm.compare(L"acosh")) return mglApplyFuncC(str, arg, head, acoshc,fns);
1394 else if(!nm.compare(L"atanh")) return mglApplyFuncC(str, arg, head, atanhc,fns);
1395 else if(!nm.compare(L"arg")) return mglApplyFuncC(str, arg, head, argc,fns);
1396 else if(!nm.compare(L"abs")) return mglApplyFuncC(str, arg, head, absc,fns);
1397 }
1398 else if(nm[0]=='c')
1399 {
1400 if(!nm.compare(L"cos")) return mglApplyFuncC(str, arg, head, cosc,fns);
1401 else if(!nm.compare(L"cosh") || !nm.compare(L"ch")) return mglApplyFuncC(str, arg, head, coshc,fns);
1402 else if(!nm.compare(L"conj")) return mglApplyFuncC(str, arg, head, conjc,fns);
1403 else if(!nm.compare(L"cmplx") && n>0)
1404 return mglApplyOperC(str.substr(0,n),str.substr(n+1),arg, head, cmplxc,fns);
1405 }
1406 else if(!nm.compare(L"exp")) return mglApplyFuncC(str, arg, head, expc,fns);
1407 else if(nm[0]=='l')
1408 {
1409 if(!nm.compare(L"log") || !nm.compare(L"ln")) return mglApplyFuncC(str, arg, head, logc,fns);
1410 else if(!nm.compare(L"lg")) return mglApplyFuncC(str, arg, head, lgc,fns);
1411 }
1412 else if(nm[0]=='s')
1413 {
1414 if(!nm.compare(L"sqrt")) return mglApplyFuncC(str, arg, head, sqrtc,fns);
1415 else if(!nm.compare(L"sin")) return mglApplyFuncC(str, arg, head, sinc,fns);
1416 else if(!nm.compare(L"sinh") || !nm.compare(L"sh")) return mglApplyFuncC(str, arg, head, sinhc,fns);
1417 else if(!nm.compare(L"sum") && n>0)
1418 {
1419 HADT a=NULL;
1420 HMDT b=mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1421 long m = long(b->a[0]+0.5);
1422 const char *s = head.size()>0?head[head.size()-1]->s.s:"";
1423 if(m>0)
1424 {
1425 std::vector<mglDataA*> hh(head);
1426 int in=0;
1427 if(s[0]=='_' && s[1]>='i' && s[1]<'z') in = s[1]-'i'+1;
1428 char name[3] = {'_',char('i'+in),0};
1429 b->s = name; b->Create(1,1,1); hh.push_back(b);
1430 a = mglFormulaCalcAC(str.substr(0,n), arg, hh, fns);
1431 long nn = a->GetNN();
1432 for(long i=1;i<m;i++)
1433 {
1434 b->a[0] = i;
1435 HADT c = mglFormulaCalcAC(str.substr(0,n), arg, hh, fns);
1436 for(long j=0;j<nn;j++) a->a[j] += c->a[j];
1437 mgl_delete_datac(c);
1438 }
1439 mgl_delete_data(b);
1440 }
1441 else
1442 { a = new mglDataC; a->a[0]=NAN; }
1443 return a;
1444 }
1445 else if(!nm.compare(L"spline"))
1446 {
1447 HADT a = mglFormulaCalcAC(str.substr(0,n), arg, head, fns);
1448 long n1 = mglFindInText(str,",",1), n2 = mglFindInText(str,",",2);
1449 HMDT ix, iy=NULL, iz=NULL;
1450 if(n1>0)
1451 {
1452 ix = mglFormulaCalcA(str.substr(n+1,n1), arg, head, fns);
1453 if(n2>0)
1454 {
1455 iy = mglFormulaCalcA(str.substr(n1+1,n2), arg, head, fns);
1456 iz = mglFormulaCalcA(str.substr(n2+1), arg, head, fns);
1457 }
1458 else iy = mglFormulaCalcA(str.substr(n1+1), arg, head, fns);
1459
1460 }
1461 else ix = mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1462 HADT res = mgl_datac_evaluate(a, ix,iy,iz, false);
1463 if(a) mgl_delete_datac(a);
1464 if(ix) mgl_delete_data(ix);
1465 if(iy) mgl_delete_data(iy);
1466 if(iz) mgl_delete_data(iz);
1467 return res;
1468 }
1469 }
1470 else if(nm[0]=='t')
1471 {
1472 if(!nm.compare(L"tg") || !nm.compare(L"tan")) return mglApplyFuncC(str, arg, head, tanc,fns);
1473 else if(!nm.compare(L"tanh") || !nm.compare(L"th")) return mglApplyFuncC(str, arg, head, tanhc,fns);
1474 }
1475 else if(!nm.compare(L"pow") && n>0)
1476 return mglApplyOperC(str.substr(0,n),str.substr(n+1),arg, head, powc,fns);
1477 else if(!nm.compare(L"random"))
1478 { HADT res=mglFormulaCalcC(str, arg, head); dual *a = res->a;
1479 for(long i=0;i<res->GetNN();i++) a[i] = dual(mgl_rnd(), mgl_rnd());
1480 return res; }
1481 else if(!nm.compare(L"hypot"))
1482 return mglApplyOperC(str.substr(0,n),str.substr(n+1),arg, head, hypotc,fns);
1483 else if(!nm.compare(L"real")) return mglApplyFuncC(str, arg, head, realc,fns);
1484 else if(!nm.compare(L"imag")) return mglApplyFuncC(str, arg, head, imagc,fns);
1485 else if(!nm.compare(L"norm")) return mglApplyFuncC(str, arg, head, normc,fns);
1486 else if(!nm.compare(L"value"))
1487 {
1488 HADT a = mglFormulaCalcAC(str.substr(0,n), arg, head, fns);
1489 long n1 = mglFindInText(str,",",1), n2 = mglFindInText(str,",",2);
1490 HMDT ix, iy=NULL, iz=NULL;
1491 if(n1>0)
1492 {
1493 ix = mglFormulaCalcA(str.substr(n+1,n1), arg, head, fns);
1494 if(n2>0)
1495 {
1496 iy = mglFormulaCalcA(str.substr(n1+1,n2), arg, head, fns);
1497 iz = mglFormulaCalcA(str.substr(n2+1), arg, head, fns);
1498 }
1499 else iy = mglFormulaCalcA(str.substr(n1+1), arg, head, fns);
1500
1501 }
1502 else ix = mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1503 HADT res = mgl_datac_subdata_ext(a, ix,iy,iz);
1504 if(a) mgl_delete_datac(a);
1505 if(ix) mgl_delete_data(ix);
1506 if(iy) mgl_delete_data(iy);
1507 if(iz) mgl_delete_data(iz);
1508 return res;
1509 }
1510 else if(!nm.compare(L"dsum") && n>0)
1511 {
1512 HADT a=NULL;
1513 HMDT b=mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1514 long m = long(b->a[0]+0.5);
1515 const char *s = head.size()>0?head[head.size()-1]->s.s:"";
1516 if(m>0)
1517 {
1518 std::vector<mglDataA*> hh(head);
1519 int in=0, zn=1;
1520 if(s[0]=='_' && s[1]>='i' && s[1]<'z') in = s[1]-'i'+1;
1521 char name[3] = {'_',char('i'+in),0};
1522 b->s = name; b->Create(1,1,1); hh.push_back(b);
1523 a = mglFormulaCalcAC(str.substr(0,n), arg, hh, fns);
1524 long nn = a->GetNN();
1525 for(long i=1;i<m;i++)
1526 {
1527 b->a[0] = i; zn *= -1;
1528 HADT c = mglFormulaCalcAC(str.substr(0,n), arg, hh, fns);
1529 for(long j=0;j<nn;j++) a->a[j] += mreal(zn)*c->a[j];
1530 mgl_delete_datac(c);
1531 }
1532 mgl_delete_data(b);
1533 }
1534 else
1535 { a = new mglDataC; a->a[0]=NAN; }
1536 return a;
1537 }
1538 else if(!nm.compare(L"prod") && n>0)
1539 {
1540 HADT a=NULL;
1541 HMDT b=mglFormulaCalcA(str.substr(n+1), arg, head, fns);
1542 long m = long(b->a[0]+0.5);
1543 const char *s = head.size()>0?head[head.size()-1]->s.s:"";
1544 if(m>0)
1545 {
1546 std::vector<mglDataA*> hh(head);
1547 int in=0;
1548 if(s[0]=='_' && s[1]>='i' && s[1]<'z') in = s[1]-'i'+1;
1549 char name[3] = {'_',char('i'+in),0};
1550 b->s = name; b->Create(1,1,1); hh.push_back(b);
1551 a = mglFormulaCalcAC(str.substr(0,n), arg, hh, fns);
1552 long nn = a->GetNN();
1553 for(long i=1;i<m;i++)
1554 {
1555 b->a[0] = i;
1556 HADT c = mglFormulaCalcAC(str.substr(0,n), arg, hh, fns);
1557 for(long j=0;j<nn;j++) a->a[j] *= c->a[j];
1558 mgl_delete_datac(c);
1559 }
1560 mgl_delete_data(b);
1561 }
1562 else
1563 { a = new mglDataC; a->a[0]=NAN; }
1564 return a;
1565 }
1566 else if(nm[0]=='f' && nm[1]=='n' && fns.size()>1)
1567 {
1568 HADT a=NULL, b=NULL, r=NULL;
1569 std::vector<mglDataA*> hh(head);
1570 std::vector<std::wstring> tmp; // disable recursion
1571 if(n>0)
1572 {
1573 a = mglFormulaCalcAC(str.substr(0,n), arg, head, fns);
1574 if(a) { a->s = "_1"; hh.push_back(a); }
1575 b = mglFormulaCalcAC(str.substr(n+1), arg, head, fns);
1576 if(b) { b->s = "_2"; hh.push_back(b); }
1577 }
1578 else
1579 {
1580 a = mglFormulaCalcAC(str, arg, head, fns);
1581 if(a) { a->s = "_1"; hh.push_back(a); }
1582 }
1583 if(!nm.compare(L"fn1") && fns.size()>1) r = mglFormulaCalcAC(fns[1], arg, hh, tmp);
1584 if(!nm.compare(L"fn2") && fns.size()>3) r = mglFormulaCalcAC(fns[2], arg, hh, tmp);
1585 if(!nm.compare(L"fn3") && fns.size()>3) r = mglFormulaCalcAC(fns[3], arg, hh, tmp);
1586 if(!nm.compare(L"fn4") && fns.size()>4) r = mglFormulaCalcAC(fns[4], arg, hh, tmp);
1587 if(!nm.compare(L"fn5") && fns.size()>5) r = mglFormulaCalcAC(fns[5], arg, hh, tmp);
1588 if(!nm.compare(L"fn6") && fns.size()>6) r = mglFormulaCalcAC(fns[6], arg, hh, tmp);
1589 if(!nm.compare(L"fn7") && fns.size()>7) r = mglFormulaCalcAC(fns[7], arg, hh, tmp);
1590 if(!nm.compare(L"fn8") && fns.size()>8) r = mglFormulaCalcAC(fns[8], arg, hh, tmp);
1591 if(!nm.compare(L"fn9") && fns.size()>9) r = mglFormulaCalcAC(fns[9], arg, hh, tmp);
1592 if(a) mgl_delete_datac(a);
1593 if(b) mgl_delete_datac(b);
1594 if(!r) { r = new mglDataC; r->a[0]=NAN; }
1595 return r;
1596 }
1597 }
1598 HADT res = new mglDataC; res->a[0]=NAN; return res;
1599 }
1600 //-----------------------------------------------------------------------------
1601