1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Application Programming Interface //
9 // //
10 // Library: SAGA_API //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // mat_formula.cpp //
15 // //
16 // Copyright (C) 2002 by Andre Ringeler //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. //
22 // //
23 // This library is free software; you can redistribute //
24 // it and/or modify it under the terms of the GNU Lesser //
25 // General Public License as published by the Free //
26 // Software Foundation, either version 2.1 of the //
27 // License, or (at your option) any later version. //
28 // //
29 // This library is distributed in the hope that it will //
30 // be useful, but WITHOUT ANY WARRANTY; without even the //
31 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
32 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
33 // License for more details. //
34 // //
35 // You should have received a copy of the GNU Lesser //
36 // General Public License along with this program; if //
37 // not, see <http://www.gnu.org/licenses/>. //
38 // //
39 //-------------------------------------------------------//
40 // //
41 // e-mail: aringel@saga-gis.org //
42 // //
43 // contact: Andre Ringeler //
44 // Institute of Geography //
45 // University of Goettingen //
46 // Goldschmidtstr. 5 //
47 // 37077 Goettingen //
48 // Germany //
49 // //
50 //-------------------------------------------------------//
51 // //
52 // Based on //
53 // FORMULC.C 2.22 as of 2/19/98 //
54 // Copyright (c) 1993 - 98 by Harald Helfgott //
55 // //
56 // Modified for Grid Data by Andre Ringeler 2001 //
57 // Modified for Function-Fitting by Andre Ringeler 2002 //
58 // Converted to C++ by Andre Ringeler 2002 //
59 // //
60 // Modified to fit SAGA needs by Olaf Conrad 2002 //
61 // //
62 ///////////////////////////////////////////////////////////
63
64 //---------------------------------------------------------
65 #include <string.h>
66 #include <stdlib.h>
67 #include <stdio.h>
68 #include <stdarg.h>
69 #include <errno.h>
70 #include <ctype.h>
71
72 #include "mat_tools.h"
73 #include "grid.h"
74
75
76 ///////////////////////////////////////////////////////////
77 // //
78 // //
79 // //
80 ///////////////////////////////////////////////////////////
81
82 //---------------------------------------------------------
83 #define MAX_PARMS 3
84 #define MAX_CTABLE 255
85
86 #define STD_FNC_NUM 19
87
88 //---------------------------------------------------------
89 #define GET_VALUE_BUFSIZE 500
90
91 //---------------------------------------------------------
92 #define EPSILON 1e-9
93
94
95 ///////////////////////////////////////////////////////////
96 // //
97 // //
98 // //
99 ///////////////////////////////////////////////////////////
100
101 //---------------------------------------------------------
f_atan2(double x,double val)102 static double f_atan2(double x, double val)
103 {
104 return( atan2(x, val) );
105 }
106
107 //---------------------------------------------------------
f_pow(double x,double val)108 static double f_pow(double x, double val)
109 {
110 return( pow(x, val) );
111 }
112
113 //---------------------------------------------------------
f_gt(double x,double val)114 static double f_gt(double x, double val)
115 {
116 return( x > val ? 1.0 : 0.0 );
117 }
118
119 //---------------------------------------------------------
f_lt(double x,double val)120 static double f_lt(double x, double val)
121 {
122 return( x < val ? 1.0 : 0.0 );
123 }
124
125 //---------------------------------------------------------
f_eq(double x,double val)126 static double f_eq(double x, double val)
127 {
128 return( fabs(x - val) < EPSILON ? 1.0 : 0.0 );
129 }
130
131 //---------------------------------------------------------
f_min(double a,double b)132 static double f_min(double a, double b)
133 {
134 return( a < b ? a : b );
135 }
136
137 //---------------------------------------------------------
f_max(double a,double b)138 static double f_max(double a, double b)
139 {
140 return( a > b ? a : b );
141 }
142
143 //---------------------------------------------------------
f_pi(void)144 static double f_pi(void)
145 {
146 return( M_PI );
147 }
148
149 //---------------------------------------------------------
f_int(double x)150 static double f_int(double x)
151 {
152 return( (int)(x) );
153 }
154
155 //---------------------------------------------------------
f_sqr(double x)156 static double f_sqr(double x)
157 {
158 return( x*x );
159 }
160
161 //---------------------------------------------------------
f_fmod(double x,double val)162 static double f_fmod(double x, double val)
163 {
164 return( fmod(x, val) );
165 }
166
167 //---------------------------------------------------------
f_rand_u(double min,double max)168 static double f_rand_u(double min, double max)
169 {
170 return( CSG_Random::Get_Uniform(min, max) );
171 }
172
173 //---------------------------------------------------------
f_rand_g(double mean,double stdv)174 static double f_rand_g(double mean, double stdv)
175 {
176 return( CSG_Random::Get_Gaussian(mean, stdv) );
177 }
178
179 //---------------------------------------------------------
f_and(double x,double y)180 static double f_and(double x, double y)
181 {
182 return( x != 0.0 && y != 0.0 ? 1.0 : 0.0 );
183 }
184
185 //---------------------------------------------------------
f_or(double x,double y)186 static double f_or(double x, double y)
187 {
188 return( x != 0.0 || y != 0.0 ? 1.0 : 0.0 );
189 }
190
191 //---------------------------------------------------------
f_ifelse(double condition,double x,double y)192 static double f_ifelse(double condition, double x, double y)
193 {
194 return( condition ? x : y );
195 // return( fabs(condition) >= EPSILON ? x : y );
196 }
197
198
199 ///////////////////////////////////////////////////////////
200 // //
201 // //
202 // //
203 ///////////////////////////////////////////////////////////
204
205 //---------------------------------------------------------
206 static CSG_Formula::TSG_Function gSG_Functions[MAX_CTABLE] =
207 {
208 { "exp" , (TSG_Formula_Function_1)exp , 1, false }, // 1
209 { "ln" , (TSG_Formula_Function_1)log , 1, false }, // 2
210 { "sin" , (TSG_Formula_Function_1)sin , 1, false }, // 3
211 { "cos" , (TSG_Formula_Function_1)cos , 1, false }, // 4
212 { "tan" , (TSG_Formula_Function_1)tan , 1, false }, // 5
213 { "asin" , (TSG_Formula_Function_1)asin , 1, false }, // 6
214 { "acos" , (TSG_Formula_Function_1)acos , 1, false }, // 7
215 { "atan" , (TSG_Formula_Function_1)atan , 1, false }, // 8
216 { "atan2" , (TSG_Formula_Function_1)f_atan2 , 2, false }, // 9
217 { "abs" , (TSG_Formula_Function_1)fabs , 1, false }, // 10
218 { "sqrt" , (TSG_Formula_Function_1)sqrt , 1, false }, // 11
219 { "gt" , (TSG_Formula_Function_1)f_gt , 2, false }, // 12
220 { "lt" , (TSG_Formula_Function_1)f_lt , 2, false }, // 13
221 { "eq" , (TSG_Formula_Function_1)f_eq , 2, false }, // 14
222 { "pi" , (TSG_Formula_Function_1)f_pi , 0, false }, // 15
223 { "int" , (TSG_Formula_Function_1)f_int , 1, false }, // 16
224 { "mod" , (TSG_Formula_Function_1)f_fmod , 2, false }, // 17
225 { "ifelse", (TSG_Formula_Function_1)f_ifelse, 3, false }, // 18
226 { "log" , (TSG_Formula_Function_1)log10 , 1, false }, // 19
227 { "pow" , (TSG_Formula_Function_1)f_pow , 2, false }, // 20
228 { "sqr" , (TSG_Formula_Function_1)f_sqr , 1, false }, // 21
229 { "rand_u", (TSG_Formula_Function_1)f_rand_u, 2, true }, // 22
230 { "rand_g", (TSG_Formula_Function_1)f_rand_g, 2, true }, // 23
231 { "and" , (TSG_Formula_Function_1)f_and , 2, false }, // 24
232 { "or" , (TSG_Formula_Function_1)f_or , 2, false }, // 25
233 { "min" , (TSG_Formula_Function_1)f_min , 2, false }, // 26
234 { "max" , (TSG_Formula_Function_1)f_max , 2, false }, // 27
235 { NULL , (TSG_Formula_Function_1) NULL , 0, false }
236 };
237
238
239 ///////////////////////////////////////////////////////////
240 // //
241 // //
242 // //
243 ///////////////////////////////////////////////////////////
244
245 //---------------------------------------------------------
CSG_Formula(void)246 CSG_Formula::CSG_Formula(void)
247 {
248 m_Formula.code = NULL;
249 m_Formula.ctable = NULL;
250
251 m_bError = false;
252
253 m_ctable = NULL;
254 m_error = NULL;
255
256 //-----------------------------------------------------
257 m_Functions = (TSG_Function *)SG_Calloc(MAX_CTABLE, sizeof(TSG_Function));
258
259 for(int i=0; i<MAX_CTABLE; i++)
260 {
261 m_Functions[i].Name = gSG_Functions[i].Name;
262 m_Functions[i].Function = gSG_Functions[i].Function;
263 m_Functions[i].nParameters = gSG_Functions[i].nParameters;
264 m_Functions[i].bVarying = gSG_Functions[i].bVarying;
265 }
266 }
267
268 //---------------------------------------------------------
~CSG_Formula(void)269 CSG_Formula::~CSG_Formula(void)
270 {
271 Destroy();
272
273 SG_Free(m_Functions);
274 }
275
276 //---------------------------------------------------------
Destroy(void)277 bool CSG_Formula::Destroy(void)
278 {
279 SG_FREE_SAFE(m_Formula.code);
280 SG_FREE_SAFE(m_Formula.ctable);
281
282 m_bError = false;
283
284 return( true );
285 }
286
287
288 ///////////////////////////////////////////////////////////
289 // //
290 // //
291 // //
292 ///////////////////////////////////////////////////////////
293
294 //---------------------------------------------------------
Get_Help_Operators(bool bHTML,const CSG_String Additional[][2])295 CSG_String CSG_Formula::Get_Help_Operators(bool bHTML, const CSG_String Additional[][2])
296 {
297 const int nOperators = 35;
298
299 CSG_String Operators[nOperators][2] =
300 {
301 { "+" , _TL("Addition") },
302 { "-" , _TL("Subtraction") },
303 { "*" , _TL("Multiplication") },
304 { "/" , _TL("Division") },
305 { "abs(x)" , _TL("Absolute Value") },
306 { "mod(x, y)" , _TL("Returns the floating point remainder of x/y") },
307 { "int(x)" , _TL("Returns the integer part of floating point value x") },
308 { "sqr(x)" , _TL("Square") },
309 { "sqrt(x)" , _TL("Square Root") },
310 { "exp(x)" , _TL("Exponential") },
311 { "pow(x, y)" , _TL("Returns x raised to the power of y") },
312 { "x ^ y" , _TL("Returns x raised to the power of y") },
313 { "ln(x)" , _TL("Natural Logarithm") },
314 { "log(x)" , _TL("Base 10 Logarithm") },
315 { "pi()" , _TL("Returns the value of Pi") },
316 { "sin(x)" , _TL("Sine") },
317 { "cos(x)" , _TL("Cosine") },
318 { "tan(x)" , _TL("Tangent") },
319 { "asin(x)" , _TL("Arcsine") },
320 { "acos(x)" , _TL("Arccosine") },
321 { "atan(x)" , _TL("Arctangent") },
322 { "atan2(x, y)" , _TL("Arctangent of x/y") },
323 { "min(x, y)" , _TL("Returns the minimum of values x and y") },
324 { "max(x, y)" , _TL("Returns the maximum of values x and y") },
325 { "gt(x, y)" , _TL("Returns true (1), if x is greater than y, else false (0)") },
326 { "x > y" , _TL("Returns true (1), if x is greater than y, else false (0)") },
327 { "lt(x, y)" , _TL("Returns true (1), if x is less than y, else false (0)") },
328 { "x < y" , _TL("Returns true (1), if x is less than y, else false (0)") },
329 { "eq(x, y)" , _TL("Returns true (1), if x equals y, else false (0)") },
330 { "x = y" , _TL("Returns true (1), if x equals y, else false (0)") },
331 { "and(x, y)" , _TL("Returns true (1), if both x and y are true (i.e. not 0)") },
332 { "or(x, y)" , _TL("Returns true (1), if at least one of both x and y is true (i.e. not 0)") },
333 { "ifelse(c, x, y)", _TL("Returns x, if condition c is true (i.e. not 0), else y") },
334 { "rand_u(x, y)" , _TL("Random number, uniform distribution with minimum x and maximum y") },
335 { "rand_g(x, y)" , _TL("Random number, Gaussian distribution with mean x and standard deviation y") }
336 };
337
338 //-----------------------------------------------------
339 int i;
340 CSG_String s;
341
342 if( bHTML )
343 {
344 s += "<table border=\"0\">";
345
346 for(i=0; i<nOperators; i++)
347 {
348 CSG_String op = Operators[i][0]; op.Replace("<", "<");
349
350 s += "<tr><td><b>" + op + "</b></td><td>" + Operators[i][1] + "</td></tr>";
351 }
352
353 if( Additional )
354 {
355 for(i=0; !Additional[i][0].is_Empty(); i++)
356 {
357 CSG_String op = Additional[i][0]; op.Replace("<", "<");
358
359 s += "<tr><td><b>" + op + "</b></td><td>" + Additional[i][1] + "</td></tr>";
360 }
361 }
362
363 s += "</table>";
364 }
365 else
366 {
367 for(i=0; i<nOperators; i++)
368 {
369 s += Operators[i][0] + " - " + Operators[i][1] + "\n";
370 }
371
372 if( Additional )
373 {
374 for(i=0; !Additional[i][0].is_Empty(); i++)
375 {
376 s += Additional[i][0] + " - " + Additional[i][1] + "\n";
377 }
378 }
379 }
380
381 return( s );
382 }
383
384 //---------------------------------------------------------
Get_Error(CSG_String & Message)385 bool CSG_Formula::Get_Error(CSG_String &Message)
386 {
387 if( m_bError )
388 {
389 Message = CSG_String::Format("%s %s %d\n", _TL("Error in formula"), _TL("at position"), m_Error_Position);
390
391 if( m_Error_Position < 0 || m_Error_Position >= (int)m_sFormula.Length() )
392 {
393 Message += m_sFormula;
394 }
395 else
396 {
397 Message += m_sFormula.Left (m_Error_Position) + " ["
398 + m_sFormula [m_Error_Position] + "] "
399 + m_sFormula.Right(m_sFormula.Length() - (m_Error_Position + 1));
400 }
401
402 Message += "\n";
403 Message += m_sError;
404 Message += "\n";
405
406 return( true );
407 }
408
409 return( false );
410 }
411
412 //---------------------------------------------------------
_Set_Error(const CSG_String & Error)413 void CSG_Formula::_Set_Error(const CSG_String &Error)
414 {
415 if( Error.is_Empty() )
416 {
417 m_bError = false;
418 m_sError .Clear();
419 }
420 else
421 {
422 m_bError = true;
423 m_sError = Error;
424 }
425 }
426
427
428 ///////////////////////////////////////////////////////////
429 // //
430 ///////////////////////////////////////////////////////////
431
432 //---------------------------------------------------------
Set_Variable(char Var,double Value)433 void CSG_Formula::Set_Variable(char Var, double Value)
434 {
435 m_Parameters[Var - 'a'] = Value;
436 }
437
438 //---------------------------------------------------------
Set_Formula(const CSG_String & Formula)439 bool CSG_Formula::Set_Formula(const CSG_String &Formula)
440 {
441 if( Formula.Length() > 0 )
442 {
443 Destroy();
444
445 m_sFormula = Formula;
446 m_Formula = _Translate(Formula, "abcdefghijklmnopqrstuvwxyz", &m_Length, &m_Error_Position);
447
448 if( m_Formula.code != NULL )
449 {
450 return( true );
451 }
452 }
453
454 Destroy();
455
456 return( false );
457 }
458
459
460 ///////////////////////////////////////////////////////////
461 // //
462 ///////////////////////////////////////////////////////////
463
464 //---------------------------------------------------------
Get_Value(void) const465 double CSG_Formula::Get_Value(void) const
466 {
467 return( _Get_Value(m_Parameters, m_Formula) );
468 }
469
470 //---------------------------------------------------------
Get_Value(double x) const471 double CSG_Formula::Get_Value(double x) const
472 {
473 double Parameters[32];
474
475 memcpy(Parameters, m_Parameters, 32 * sizeof(double));
476
477 Parameters['x'-'a'] = x;
478
479 return( _Get_Value(Parameters, m_Formula) );
480 }
481
482 //---------------------------------------------------------
Get_Value(const CSG_Vector & Values) const483 double CSG_Formula::Get_Value(const CSG_Vector &Values) const
484 {
485 return( Get_Value(Values.Get_Data(), Values.Get_N()) );
486 }
487
488 //---------------------------------------------------------
Get_Value(double * Values,int nValues) const489 double CSG_Formula::Get_Value(double *Values, int nValues) const
490 {
491 double Parameters[32];
492
493 for(int i=0; i<nValues; i++)
494 {
495 Parameters[i] = Values[i];
496 }
497
498 return( _Get_Value(Parameters, m_Formula) );
499 }
500
501 //---------------------------------------------------------
Get_Value(const char * Args,...) const502 double CSG_Formula::Get_Value(const char *Args, ...) const
503 {
504 double Parameters[32];
505
506 va_list ap;
507
508 va_start(ap, Args);
509
510 while( *Args )
511 {
512 Parameters[(*Args++) - 'a'] = va_arg(ap, double);
513 }
514
515 va_end(ap);
516
517 return( _Get_Value(Parameters, m_Formula) );
518 }
519
520 //---------------------------------------------------------
_Get_Value(const double * Parameters,TSG_Formula func) const521 double CSG_Formula::_Get_Value(const double *Parameters, TSG_Formula func) const
522 {
523 double x, y, z, buffer[GET_VALUE_BUFSIZE];
524
525 register double *bufp = buffer; // points to the first free space in the buffer
526 register char *function = func.code;
527 register double *ctable = func.ctable;
528 register double result;
529
530 if( !function )
531 {
532 return( 0 ); // _Set_Error(_TL("empty coded function"));
533 }
534
535 for( ; ; )
536 {
537 switch( *function++ )
538 {
539 case '\0':
540 return( buffer[0] );
541
542 case 'D':
543 *bufp++ = ctable[*function++];
544 break;
545
546 case 'V':
547 *bufp++ = Parameters[(*function++) - 'a'];
548 break;
549
550 case 'M':
551 result = -(*--bufp);
552 *bufp++ = result;
553 break;
554
555 case '+':
556 y = *(--bufp);
557 result = y + *(--bufp);
558 *bufp++ = result;
559 break;
560
561 case '-':
562 y = *--bufp;
563 result = *(--bufp) - y;
564 *bufp++ = result;
565 break;
566
567 case '*':
568 y = *(--bufp);
569 result = *(--bufp) * y;
570 *bufp++ = result;
571 break;
572
573 case '/':
574 y = *--bufp;
575 result = *(--bufp) / y;
576 *bufp++ = result;
577 break;
578
579 case '^':
580 y = *--bufp;
581 result = pow(*(--bufp), y);
582 *bufp++ = result;
583 break;
584
585 case '=':
586 y = *--bufp;
587 result = y == *(--bufp) ? 1.0 : 0.0;
588 *bufp++ = result;
589 break;
590
591 case '>':
592 y = *--bufp;
593 result = y < *(--bufp) ? 1.0 : 0.0;
594 *bufp++ = result;
595 break;
596
597 case '<':
598 y = *--bufp;
599 result = y > *(--bufp) ? 1.0 : 0.0;
600 *bufp++ = result;
601 break;
602
603 case '&':
604 y = *--bufp;
605 result = y && *(--bufp) ? 1.0 : 0.0;
606 *bufp++ = result;
607 break;
608
609 case '|':
610 y = *--bufp;
611 result = y || *(--bufp) ? 1.0 : 0.0;
612 *bufp++ = result;
613 break;
614
615 case 'F':
616 switch (m_Functions[*function].nParameters)
617 {
618 case 0:
619 *bufp++ = ((TSG_Formula_Function_0)m_Functions[*function++].Function)();
620 break;
621
622 case 1:
623 x = *--bufp;
624 *bufp++ = ((TSG_Formula_Function_1)m_Functions[*function++].Function)(x);
625 break;
626
627 case 2:
628 y = *--bufp;
629 x = *--bufp;
630 *bufp++ = ((TSG_Formula_Function_2)m_Functions[*function++].Function)(x, y);
631 break;
632
633 case 3:
634 z = *--bufp;
635 y = *--bufp;
636 x = *--bufp;
637 *bufp++ = ((TSG_Formula_Function_3)m_Functions[*function++].Function)(x, y, z);
638 break;
639
640 default:
641 return( 0 ); // _Set_Error(_TL("I2: too many parameters"));
642 }
643 break;
644
645 default:
646 return( 0 ); // _Set_Error(_TL("I1: unrecognizable operator"));
647 }
648 }
649
650 // if( (bufp - buffer) != 1 ) // _Set_Error(_TL("I3: corrupted buffer"));
651
652 return( buffer[0] );
653 }
654
655
656 ///////////////////////////////////////////////////////////
657 // //
658 // //
659 // //
660 ///////////////////////////////////////////////////////////
661
662 //---------------------------------------------------------
Get_Used_Variables(void)663 const char * CSG_Formula::Get_Used_Variables(void)
664 {
665 static CSG_String ret;
666
667 ret.Clear();
668
669 for(int i=0; i<'z'-'a'; i++)
670 {
671 if( m_Vars_Used[i] == true )
672 {
673 ret.Append((char)(i + 'a'));
674 }
675 }
676
677 return( ret );
678 }
679
680
681 ///////////////////////////////////////////////////////////
682 // //
683 // //
684 // //
685 ///////////////////////////////////////////////////////////
686
687 //---------------------------------------------------------
688 // int varying; Does the result of the function vary
689 // even when the parameters stay the same?
690 // varying = 1 for e.g. random - number generators.
691 // Result: 0 is rendered if there is an error
692 // 1 is rendered otherwise
693 //
Add_Function(const char * Name,TSG_Formula_Function_1 Function,int nParameters,bool bVarying)694 bool CSG_Formula::Add_Function(const char *Name, TSG_Formula_Function_1 Function, int nParameters, bool bVarying)
695 {
696 if( nParameters < 0 || nParameters > 3 )
697 {
698 _Set_Error(_TL("invalid number of parameters"));
699
700 return( false );
701 }
702
703 TSG_Function *pFunction;
704
705 for(pFunction=m_Functions; pFunction->Function && strcmp(Name, pFunction->Name); pFunction++)
706 {}
707
708 if( pFunction->Function != NULL ) // old function is superseded
709 {
710 pFunction->Function = Function;
711 pFunction->nParameters = nParameters;
712 pFunction->bVarying = bVarying;
713
714 _Set_Error();
715
716 return( true );
717 }
718
719 if( (pFunction - m_Functions) >= MAX_CTABLE - 1 )
720 {
721 _Set_Error(_TL("function table full"));
722
723 return( false );
724 }
725
726 pFunction->Name = Name;
727 pFunction->Function = Function;
728 pFunction->nParameters = nParameters;
729 pFunction->bVarying = bVarying;
730
731 _Set_Error();
732
733 return( true );
734 }
735
736
737 ///////////////////////////////////////////////////////////
738 // //
739 ///////////////////////////////////////////////////////////
740
741 //---------------------------------------------------------
_Get_Function(int i,char * Name,int * nParameters,int * bVarying)742 int CSG_Formula::_Get_Function(int i, char *Name, int *nParameters, int *bVarying)
743 {
744 if( !m_Functions[i].Function )
745 {
746 _Set_Error(_TL("index out of bounds"));
747
748 return( 0 );
749 }
750
751 strcpy(Name , m_Functions[i].Name);
752 *nParameters = m_Functions[i].nParameters;
753 *bVarying = m_Functions[i].bVarying ? 1 : 0;
754
755 _Set_Error();
756
757 return( 1 );
758 }
759
760 //---------------------------------------------------------
761 // If the function exists, _Get_Function() returns the index
762 // of its name in the table. Otherwise, it returns -1.
763 //
_Get_Function(const char * Name)764 int CSG_Formula::_Get_Function(const char *Name)
765 {
766 TSG_Function *pFunction;
767
768 for(pFunction=m_Functions; pFunction->Function && strcmp(Name, pFunction->Name); pFunction++)
769 {}
770
771 if( pFunction->Function == NULL )
772 {
773 _Set_Error(_TL("function not found"));
774
775 return( -1 );
776 }
777
778 _Set_Error();
779
780 return( (int)(pFunction - m_Functions) );
781 }
782
783
784 ///////////////////////////////////////////////////////////
785 // //
786 // //
787 // //
788 ///////////////////////////////////////////////////////////
789
790 //---------------------------------------------------------
_is_Operand(char c)791 inline int CSG_Formula::_is_Operand(char c)
792 {
793 return( (c == '+')
794 || (c == '-')
795 || (c == '*')
796 || (c == '/')
797 || (c == '^')
798 || (c == '=')
799 || (c == '<')
800 || (c == '>')
801 || (c == '&')
802 || (c == '|')
803 );
804 }
805
806 //---------------------------------------------------------
_is_Operand_Code(char c)807 inline int CSG_Formula::_is_Operand_Code(char c)
808 {
809 return( (c == '+')
810 || (c == '-')
811 || (c == '*')
812 || (c == '/')
813 || (c == '^')
814 || (c == '=')
815 || (c == '<')
816 || (c == '>')
817 || (c == '&')
818 || (c == '|')
819 || (c == 'M')
820 );
821 }
822
823 //---------------------------------------------------------
_is_Number(char c)824 inline int CSG_Formula::_is_Number(char c)
825 {
826 return( isdigit(c) || c == '.' || c == 'E' );
827 }
828
829
830 ///////////////////////////////////////////////////////////
831 // //
832 // //
833 // //
834 ///////////////////////////////////////////////////////////
835
836 ///////////////////////////////////////////////////////////
837 // Interpreting functions //
838 ///////////////////////////////////////////////////////////
839
840 //---------------------------------------------------------
_Translate(const char * sourc,const char * args,int * leng,int * error)841 CSG_Formula::TSG_Formula CSG_Formula::_Translate(const char *sourc, const char *args, int *leng, int *error)
842 {
843 const char *scan, *scarg;
844 char *result, *source, *code, *nfunc;
845 size_t size_estim;
846 double *ctable;
847 TSG_Formula returned;
848
849 //-----------------------------------------------------
850 *leng = 0;
851 *error = 0;
852 m_error = NULL;
853 returned.code = NULL;
854 returned.ctable = NULL;
855
856 //-----------------------------------------------------
857 source = (char *)SG_Malloc((strlen(sourc) + 1) * sizeof(char));
858
859 if( source == NULL )
860 {
861 _Set_Error(_TL("no memory"));
862
863 return( returned );
864 }
865
866 strcpy(source, sourc);
867
868 for(scan=source; *scan!='\0'; scan++)
869 {
870 if( islower(*scan) && !isalpha(*(scan + 1)) && (scan == source || !isalpha(*(scan - 1))) )
871 {
872 for(scarg=args; *scarg!='\0' && *scarg != *scan; scarg++)
873 {}
874
875 if( *scarg == '\0' )
876 {
877 _Set_Error(_TL("undeclared parameter"));
878
879 m_error = scan;
880 *error = (int)(m_error - source);
881
882 SG_Free(source);
883
884 return (returned);
885 }
886 }
887 }
888
889 //-----------------------------------------------------
890 size_estim = _max_size(source);
891
892 if( !(code = (char *)SG_Malloc(size_estim)) )
893 {
894 _Set_Error(_TL("no memory"));
895
896 *error = -1;
897
898 SG_Free(source);
899
900 return (returned);
901 }
902
903
904 //-----------------------------------------------------
905 m_pctable = 0;
906
907 if( !(m_ctable = (double *)SG_Malloc(MAX_CTABLE * sizeof(double))) )
908 {
909 _Set_Error(_TL("no memory"));
910
911 *error = -1;
912
913 SG_Free(source);
914 SG_Free(code);
915
916 return (returned);
917 }
918
919 ctable = m_ctable;
920
921 //-----------------------------------------------------
922 _Set_Error();
923
924 result = _i_trans(code, (char *)source, (char *)source + strlen(source));
925
926 if( !result || m_bError )
927 {
928 *error = (int)(m_error ? m_error - source : -1);
929
930 SG_Free(source);
931 SG_Free(code);
932 SG_Free(m_ctable);
933
934 return (returned);
935 }
936 else
937 {
938 *result = '\0';
939 *error = -1;
940 *leng = (int)(result - code);
941
942 if( ((*leng) + 1) * sizeof(char) > size_estim )
943 {
944 _Set_Error(_TL("I4: size estimate too small"));
945
946 SG_Free(source);
947
948 return( returned );
949 }
950 else if( ((*leng) + 1) * sizeof(char) < size_estim )
951 {
952 nfunc = (char *)SG_Malloc(((*leng) + 1) * sizeof(char));
953
954 if (nfunc)
955 {
956 memcpy(nfunc, code, ((*leng) + 1) * sizeof(char));
957 SG_Free(code);
958 code = nfunc;
959 }
960 }
961
962 if( m_pctable < MAX_CTABLE )
963 {
964 ctable = (double *)SG_Malloc(m_pctable * sizeof(double));
965
966 if( ctable )
967 {
968 memcpy(ctable, m_ctable, m_pctable * sizeof(double));
969
970 SG_Free(m_ctable);
971 }
972 else
973 {
974 ctable = m_ctable;
975 }
976 }
977 else
978 {
979 ctable = m_ctable;
980 }
981
982 returned.code = code;
983 returned.ctable = ctable;
984
985 _Set_Error();
986
987 SG_Free(source);
988
989 return (returned);
990 }
991 }
992
993
994 ///////////////////////////////////////////////////////////
995 // //
996 // //
997 // //
998 ///////////////////////////////////////////////////////////
999
1000 //---------------------------------------------------------
_i_trans(char * function,char * begin,char * end)1001 char * CSG_Formula::_i_trans(char *function, char *begin, char *end)
1002 {
1003 char tempch, *scan, *endf, *tempu, *temp3, *temps = NULL, *par_buf, *paramstr[MAX_PARMS];
1004 int i, pars, space, n_function;
1005 double tempd;
1006
1007 if( begin >= end )
1008 {
1009 _Set_Error(_TL("missing operand"));
1010 m_error = begin;
1011 return NULL;
1012 }
1013
1014 for(pars = 0, scan = begin; scan < end && pars >= 0; scan++)
1015 {
1016 if( *scan == '(' ) { pars++; } else
1017 if( *scan == ')' ) { pars--; }
1018 }
1019
1020 if( pars < 0 || pars > 0 )
1021 {
1022 _Set_Error(_TL("unmatched parentheses"));
1023 m_error = scan - 1;
1024 return NULL;
1025 }
1026
1027 for(pars = 0, scan = end - 1; scan >= begin; scan--)
1028 {
1029 if( *scan == '(' ) pars++; else
1030 if( *scan == ')' ) pars--; else
1031 if( !pars && (*scan == '+' || ((*scan == '-') && scan != begin)) && (scan == begin || *(scan - 1) != 'E') )
1032 break;
1033 }
1034
1035 if( scan >= begin )
1036 {
1037 if ((tempu = _i_trans(function, begin, scan)) &&
1038 (temp3 = _i_trans(tempu, scan + 1, end)))
1039 {
1040 *temp3++ = *scan;
1041 temp3 = _comp_time(function, temp3, 2);
1042 if( m_bError )
1043 return NULL;
1044 else
1045 return temp3;
1046 }
1047 else
1048 return NULL;
1049 }
1050
1051 for (pars = 0, scan = end - 1; scan >= begin; scan--)
1052 {
1053 if (*scan == '(') pars++;
1054 else if (*scan == ')')
1055 pars--;
1056 else if (!pars &&(*scan == '*' || *scan == '/'))
1057 break;
1058 }
1059 if (scan >= begin)
1060 {
1061 if ((tempu = _i_trans(function, begin, scan)) &&
1062 (temp3 = _i_trans(tempu, scan + 1, end)))
1063 {
1064 *temp3++ = *scan;
1065 temp3 = _comp_time(function, temp3, 2);
1066 if (m_bError)
1067 return NULL;
1068 else
1069 return temp3;
1070 }
1071 else
1072 return NULL;
1073 }
1074
1075 /* unary minus */
1076 if (*begin == '-')
1077 {
1078 tempu = _i_trans(function, begin + 1, end);
1079 if (tempu)
1080 {
1081 *tempu++ = 'M';
1082 tempu = _comp_time(function, tempu, 1);
1083 if (m_bError)
1084 return NULL;
1085 else
1086 return tempu;
1087 }
1088 else
1089 return NULL;
1090 }
1091
1092 for (pars = 0, scan = end - 1; scan >= begin; scan--)
1093 {
1094 if (*scan == '(') pars++;
1095 else if (*scan == ')')
1096 pars--;
1097 else if (!pars &&(*scan == '^'))
1098 break;
1099 else if (!pars &&(*scan == '='))
1100 break;
1101 else if (!pars &&(*scan == '>'))
1102 break;
1103 else if (!pars &&(*scan == '<'))
1104 break;
1105 else if (!pars &&(*scan == '&'))
1106 break;
1107 else if (!pars &&(*scan == '|'))
1108 break;
1109 }
1110
1111 if (scan >= begin)
1112 {
1113 if ((tempu = _i_trans(function, begin, scan)) &&
1114 (temp3 = _i_trans(tempu, scan + 1, end)))
1115 {
1116 *temp3++ = *scan;
1117 temp3 = _comp_time(function, temp3, 2);
1118 if (m_bError)
1119 return NULL;
1120 else
1121 return temp3;
1122 }
1123 else
1124 return NULL;
1125 }
1126
1127 /* erase white space */
1128 while (isspace(*begin))
1129 begin++;
1130 while (isspace(*(end - 1)))
1131 end--;
1132
1133 if (*begin == '(' && *(end - 1) == ')')
1134 return _i_trans(function, begin + 1, end - 1);
1135
1136 if (end == begin + 1 && islower(*begin))
1137 {
1138 *function++ = 'V';
1139 *function++ = *begin;
1140 return function;
1141 }
1142
1143 tempch = *end;
1144 *end = '\0';
1145 tempd = strtod(begin, (char **)&tempu);
1146 *end = tempch;
1147
1148 if( (char *)tempu == end )
1149 {
1150 *function++ = 'D';
1151 if (m_pctable < MAX_CTABLE)
1152 {
1153 m_ctable[m_pctable] = tempd;
1154 *function++ = (char)m_pctable++;
1155 }
1156 else
1157 {
1158 _Set_Error(_TL("too many constants"));
1159 m_error = begin;
1160 return NULL;
1161 }
1162 return function;
1163 }
1164
1165 /*function*/
1166 if (!isalpha(*begin) && *begin != '_')
1167 {
1168 _Set_Error(_TL("syntax error"));
1169 m_error = begin;
1170 return NULL;
1171 }
1172
1173 for(endf = begin + 1; endf < end &&(isalnum(*endf) || *endf == '_'); endf++)
1174 {}
1175
1176 tempch = *endf;
1177 *endf = '\0';
1178 if ((n_function = _Get_Function(begin)) == -1)
1179 {
1180 *endf = tempch;
1181 m_error = begin;
1182 return NULL;
1183 }
1184 *endf = tempch;
1185 if (*endf != '(' || *(end - 1) != ')')
1186 {
1187 _Set_Error(_TL("improper function syntax"));
1188 m_error = endf;
1189 return NULL;
1190 }
1191
1192 if( m_Functions[n_function].nParameters == 0 )
1193 {
1194 /*function without parameters(e.g. pi()) */
1195 space = 1;
1196 for (scan = endf + 1; scan <(end - 1); scan++)
1197 if (!isspace(*scan))
1198 space = 0;
1199 if (space)
1200 {
1201 *function++ = 'F';
1202 *function++ = n_function;
1203 function = _comp_time(function - 2, function, 0);
1204 if (m_bError)
1205 return NULL; /* internal error in _comp_time */
1206 else
1207 return function;
1208 }
1209 else
1210 {
1211 m_error = endf + 1;
1212 _Set_Error(_TL("too many parameters"));
1213 return NULL;
1214 }
1215 }
1216 else
1217 { /*function with parameters*/
1218 tempch = *(end - 1);
1219 *(end - 1) = '\0';
1220 par_buf = (char *)SG_Malloc(sizeof(char) * (strlen(endf + 1) + 1));
1221
1222 if (!par_buf)
1223 {
1224 _Set_Error(_TL("no memory"));
1225 m_error = NULL;
1226 return NULL;
1227 }
1228
1229 strcpy(par_buf, endf + 1);
1230 *(end - 1) = tempch;
1231
1232 for (i = 0; i < m_Functions[n_function].nParameters; i++)
1233 {
1234 if ((temps = _my_strtok((i == 0) ? par_buf : NULL)) == NULL)
1235 break;
1236 paramstr[i] = temps;
1237 }
1238
1239 if (temps == NULL)
1240 {
1241 SG_Free(par_buf);
1242 m_error = end - 2;
1243 _Set_Error(_TL("too few parameters"));
1244 return NULL;
1245 }
1246
1247 if ((temps = _my_strtok(NULL)) != NULL)
1248 {
1249 SG_Free(par_buf);
1250 m_error =(temps - par_buf) +(endf + 1);
1251 _Set_Error(_TL("too many parameters"));
1252 return NULL;
1253 }
1254
1255 tempu = function;
1256 for (i = 0; i < m_Functions[n_function].nParameters; i++)
1257 if( !(tempu = _i_trans(tempu, paramstr[i], paramstr[i] + strlen(paramstr[i]))) )
1258 {
1259 m_error =(m_error - par_buf) +(endf + 1);
1260 SG_Free(par_buf);
1261
1262 return NULL;
1263 }
1264
1265 /* OK */
1266 SG_Free(par_buf);
1267 *tempu++ = 'F';
1268 *tempu++ = n_function;
1269 tempu = _comp_time(function, tempu, m_Functions[n_function].nParameters);
1270 if (m_bError)
1271 return NULL; /* internal error in _comp_time */
1272 else
1273 return tempu;
1274 }
1275 }
1276
1277 //---------------------------------------------------------
_comp_time(char * function,char * fend,int npars)1278 char * CSG_Formula::_comp_time(char *function, char *fend, int npars)
1279 {
1280 char *scan, temp;
1281 int i;
1282 double tempd;
1283 TSG_Formula trans;
1284
1285 scan = function;
1286 for (i = 0; i < npars; i++)
1287 {
1288 if (*scan++ != 'D')
1289 return fend;
1290 scan++;
1291 }
1292
1293 if (!((scan == fend -(sizeof((char) 'F') + sizeof(char))
1294 && *(fend - 2) == 'F' && m_Functions[*(fend - 1)].bVarying == 0) ||
1295 (scan == fend - sizeof(char)
1296 && _is_Operand_Code(*(fend - 1))))
1297 )
1298 return fend;
1299
1300 temp = *fend;
1301 *fend = '\0';
1302
1303 trans.code = function;
1304 trans.ctable = m_ctable;
1305 tempd = _Get_Value(m_Parameters, trans);
1306 *fend = temp;
1307 *function++ = 'D';
1308 m_pctable -= npars;
1309 *function++ =(char) m_pctable;
1310 m_ctable[m_pctable++] = tempd;
1311
1312 return function;
1313 }
1314
1315 //---------------------------------------------------------
_max_size(const char * source)1316 int CSG_Formula::_max_size(const char *source)
1317 {
1318 int numbers = 0;
1319 int functions = 0;
1320 int operators = 0;
1321 int variables = 0;
1322
1323 const size_t var_size = 2 * sizeof(char);
1324 const size_t num_size = sizeof(char) + sizeof(double);
1325 const size_t op_size = sizeof(char);
1326 const size_t end_size = sizeof('\0');
1327
1328 for(int i=0; i<'z'-'a'; i++)
1329 {
1330 m_Vars_Used[i] = false;
1331 }
1332
1333 for(const char *scan=source; *scan; scan++)
1334 {
1335 if( isalpha(*scan) && (*scan != 'E') )
1336 {
1337 if( isalpha(*(scan + 1)) || isdigit(*(scan + 1)) )
1338 {
1339 // must be a function name (combination of letters and digits, e.g. sin(..), atan2(..))
1340 }
1341 else if( *(scan + 1) == '(' )
1342 {
1343 functions++;
1344 }
1345 else
1346 {
1347 variables++;
1348 m_Vars_Used[(int)(*scan - 'a')] = true;
1349 }
1350 }
1351 }
1352
1353 if( _is_Operand(*source) )
1354 {
1355 operators++;
1356 }
1357
1358 if( *source != '\0' )
1359 {
1360 for(const char *scan=source+1; *scan; scan++)
1361 {
1362 if( _is_Operand(*scan) && *(scan - 1) != 'E' )
1363 {
1364 operators++;
1365 }
1366 }
1367 }
1368
1369 const char *scan = source;
1370
1371 while( *scan )
1372 {
1373 if( _is_Number(*scan) || ((*scan == '+' || *scan == '-') && scan > source && *(scan - 1) == 'E') )
1374 {
1375 numbers++;
1376 scan++;
1377
1378 while( _is_Number(*scan) ||((*scan == '+' || *scan == '-') && scan > source && *(scan - 1) == 'E') )
1379 scan++;
1380 }
1381 else
1382 {
1383 scan++;
1384 }
1385 }
1386
1387 return( numbers*num_size + operators*op_size + functions*num_size + variables*var_size + end_size );
1388 }
1389
1390 //---------------------------------------------------------
_my_strtok(char * s)1391 char * CSG_Formula::_my_strtok(char *s)
1392 {
1393 static char *token = NULL;
1394
1395 if( s != NULL )
1396 token = s;
1397 else if( token != NULL )
1398 s = token;
1399 else
1400 return( NULL );
1401
1402 for(int pars=0; *s != '\0' && (*s != ',' || pars != 0); s++)
1403 {
1404 if (*s == '(') ++pars;
1405 if (*s == ')') --pars;
1406 }
1407
1408 if( *s == '\0' )
1409 {
1410 s = token;
1411 token = NULL;
1412 }
1413 else
1414 {
1415 *s = '\0';
1416 char *next_token = s + 1;
1417 s = token;
1418 token = next_token;
1419 }
1420
1421 return( s );
1422 }
1423
1424
1425 ///////////////////////////////////////////////////////////
1426 // //
1427 // //
1428 // //
1429 ///////////////////////////////////////////////////////////
1430
1431 //---------------------------------------------------------
1432