/* fie.c Copyright (c) Kapteyn Laboratorium Groningen 1990 All Rights Reserved. #> fie.dc2 Document: FIE Purpose: Describes the routines FIEINI (parses an input string which contains mathematical expression), FIEDO (which does the actual calculations) and FIEDUMP (which dumps the code generated by FIEINI). Category: MATH File: fie.c Author: K.G. Begeman Description: FIEINI parses an input string which contains a mathematical formula. The input string may contain: 1) functions; syntax ff(..); where ff is one of the following available functions: abs(x) absolute value of x sqrt(x) square root of x sin(x) sine of x asin(x) inverse sine of x cos(x) cosine of x acos(x) inverse cosine of x tan(x) tangent of x sec(x) secans of x csc(x) cosecans of x cot(x) cotangent of x atan(x) inverse tan of x atan2(x,y) inverse tan (mod 2pi) x = sin, y = cos exp(x) exponential of x ln(x) natural log of x log(x) log (base 10) of x sinh(x) hyperbolic sine of x asinh(x) inverse hyperbolic sine of x cosh(x) hyperbolic cosine of x acosh(x) inverse hyperbolic cosine of x tanh(x) hyperbolic tangent of x sech(x) hyperbolic secans of x csch(x) hyperbolic cosecans of x coth(x) hyperbolic cotangent of x bj(n,x) Bessel function of 1st kind, integer order by(n,x) Bessel function of 2nd kind, integer order bi(n,x) modified Bessel fn. of 1st kind, integer order bk(n,x) modified Bessel fn. of 2nd kind, integer order rad(x) convert x to radians deg(x) convert x to degrees erf(x) error function of x erfc(x) 1-error function sinc(x) sin(x)/x (sinc function) sign(x) sign of x (-1,0,1) step(x) returns 1 if x > 0, else 0 rect(x) returns 1 if |x| < 0.5, else 0 mod(x,y) gives remainder of x/y int(x) truncates to integer nint(x) nearest integer ranu(x,y) generates uniform noise between x and y rang(x,y) generates gaussian noise with mean x and dispersion y ranp(x) generates poisson noise with mean x ifeq(x,y,a,b) returns a if x == y, else b ifne(x,y,a,b) returns a if x != y, else b ifgt(x,y,a,b) returns a if x > y, else b ifge(x,y,a,b) returns a if x >= y, else b iflt(x,y,a,b) returns a if x < y, else b ifle(x,y,a,b) returns a if x <= y, else b ifblank(x,a,b) returns a if x == BLANK, else b Some (statistical) functions have a variable number of arguments: sum(x1,..,xn) sum of elements x1 .. xn mean(x1,..,xn) mean var(x1,..,xn) variance sdev(x1,..,xn) standard deviation adev(x1,..,xn) absolute deviation skew(x1,..,xn) skewness kurt(x1,..,xn) kurtosis median(x1,..,xn) median nblanks(x1,..,xn) number of blanks max(x1,..,xn) maximum of elements x1 .. xn min(x1,..,xn) minimum Note that n <= 128 2) constants; syntax cc; where cc is one of the following available constants: PI 3.14159.... C speed of light (SI) H Planck (SI) K Boltzmann (SI) G gravitation (SI) S Stefan-Boltzman (SI) M mass of sun (SI) P parsec (SI) BLANK Universal undefined value Note: the Hubble constant is not included. 3) operators; syntax op; where op is one of the following available operators: + addition - subtraction * multiplication / division ** power 4) parameters; syntax $n; where 1 <= n <= 128. A parameter is a value taken from a real array inserted at the position of $n in the expression. FIEDO inserts the parameters in the expression. Note: With the subroutine fiepar the programmer can define parameter names. The working of FIEINI is best explained by an example; the expression '1.0+exp($1)*sinc($3)' is decoded into the following commands: LDC 1.0 ! load constant LDP 1 ! load parameter 1 FIE EXP ! exp() LDP 3 ! load parameter 3 FIE SINC ! sinc() MUL ! * ADD ! + HLT ! end of code This code is stored internally for later processing by FIEDO. The procedure FIEDUMP displays the internally stored instructions as shown above on stdin. There are at the moment 16 different buffers for storage of instructions. Notes: 1) The calculations are all done in double precision. 2) FIEDO recognizes BLANK values. Any operation on a BLANK causes the result to be set to BLANK (except the function IFBLANK). Also all illegal operations, like dividing by zero or the logarithm of a negative value, will cause the result to be set to BLANK. 3) If you cannot find your favorite constant or function in the list, please contact the author. He might be persuaded to put it in. Updates: Mar 11, 1987: RPK original code Mar 12, 1987: KGB document created May 28, 1987: KGB RPK bug removed Oct 27, 1987: RPK KGB bug removed Aug 15, 1988: KGB RPK bug in decoding reals fixed Aug 1, 1989: KGB converted to GPS Mar 28, 1991: KGB parameter names allowed Feb 21, 1995: KGB/VOG multi-parameter functions added Feb 7, 1997: JPT named parameters supersede built-in constants. Jun 30, 1998: JPT implemented Bessel functions. Aug 11, 1998: JPT fixed bug in nint; added sec, csc, cot, sech, csch, coth, asinh, acosh, atanh, step and rect. Nov 11, 1999: VOG Avoid r1 == 0.0 in log in rang function. #< */ #include "stdio.h" /* */ #include "stdlib.h" /* */ #include "string.h" /* */ #include "ctype.h" /* */ #include "math.h" /* */ #include "gipsyc.h" /* GIPSY symbols and definitions */ #include "fblank.h" /* defines fblank_c */ #include "setfblank.h" /* defines setfblank_c */ #include "setdblank.h" /* defines setdblank_c */ #include "sortda.h" /* defines sortda_c */ #include "momsd.h" /* defines momsd_c */ #include "bessel.h" /* Bessel functions */ #include "float.h" /* Define DBL_EPSILON */ /* * While evaluating, each constant on the stack has a flag whether it * is undefined (result of an illegal operation) or blank (from data set) */ typedef struct { double d; bool u; } fiedbl; #define byte char /* * definitions/declarations for the function code */ #define MAXFIECODE 256 #define bid sizeof(double) /* the number of bytes in a double real */ #define ERR -1 /* an error occurred, too bad! */ #define HLT 0 /* we're done with evaluating */ #define ADD 1 /* + operator */ #define SUB 2 /* - operator */ #define MUL 3 /* * operator */ #define DIV 4 /* / operator */ #define NEG 5 /* negate */ #define PWR 6 /* ** operator */ #define LDP 7 /* load parameter */ #define LDC 8 /* load constant */ #define FIE 9 /* function call (i) has opcode fie + i */ static char *mnem[] = { "HLT","ADD","SUB","MUL","DIV","NEG","PWR","LDP","LDC","FIE" }; static int saveifgen[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; #define MAXID (sizeof(saveifgen)/sizeof(int)) static union { byte opcode[bid]; double c; } fiecode[MAXFIECODE], savefiecode[MAXID][MAXFIECODE]; static int savenpars[MAXID]; static int codeptr; static int opcodeptr; static float BLANK; static void fie_gencode(opc) int opc; { fiecode[codeptr].opcode[opcodeptr++] = opc; if (opcodeptr == bid) { codeptr++; opcodeptr = 0; } } static void fie_genconst(cst) double cst; { fie_gencode(LDC); if (opcodeptr != 0) codeptr++; fiecode[codeptr++].c = cst; opcodeptr = 0; } /* * functions we know */ static char *functs[] = { "SIN" , "ASIN" , "SINH" , "COS" , "ACOS" , "COSH" , "TAN" , "ATAN" , "TANH" , "ATAN2" , "RAD" , "DEG" , "PI" , "EXP" , "LN" , "LOG" , "SQRT" , "ABS" , "SINC" , "C" , "G" , "M" , "ERF" , "ERFC" , "K" , "H" , "P" , "S" , "MAX" , "MIN" , "MOD" , "INT" , "NINT" , "SIGN" , "BLANK" , "IFGT" , "IFLT" , "IFGE" , "IFLE" , "IFEQ" , "IFNE" , "RANU" , "RANG" , "RANP" , "IFBLANK", "SUM" , "MEAN" , "VAR" , "SDEV" , "ADEV" , "SKEW" , "KURT" , "MEDIAN" , "NBLANKS", "BJ" , "BY" , "BK" , "BI" , "SEC" , "CSC" , "COT" , "SECH" , "CSCH" , "COTH" , "ASINH" , "ACOSH" , "ATANH" , "STEP" , "RECT" }; static int nargs[] = { 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 2 , 1 , 1 , 0 , 1 , 1 , 1 , 1 , 1 , 1 , 0 , 0 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , -2 , -2 , 2 , 1 , 1 , 1 , 0 , 4 , 4 , 4 , 4 , 4 , 4 , 2 , 2 , 1 , 3 , -1 , -1 , -2 , -2 , -1 , -2 , -2 , -2 , -1 , 2 , 2 , 2 , 2 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 }; #define MAXFUNCTS (sizeof(nargs)/sizeof(int)) #define MAXFUNLEN 10 #define MAXARG 128 /* * definitions/declarations for the scanner and parser */ #define END 0 /* end of string */ #define PLUS 1 /* '+' */ #define MINUS 2 /* '-' */ #define TIMES 3 /* '*' */ #define DIVIDE 4 /* '/' */ #define PARAM 5 /* parameter */ #define CONST 6 /* constant */ #define FUNCT 7 /* function */ #define LPAR 8 /* '(' */ #define RPAR 9 /* ')' */ #define COMMA 10 /* ',' */ #define POWER 11 /* '**' */ #define MAXPAR MAXARG /* number of parameters */ #define MAXPARNAMELEN 16 /* length of par. names */ static int nparname = 0; static char parname[MAXPAR][MAXPARNAMELEN+1]; static char *in_line; static int pos, maxpos, errorpos; static int ch; static int sym; static int curpar, curfun, npar; static double curconst; static int oddran, error; /* * define the procedures */ static void fie_nextch(); static void fie_nextsym(); static void fie_expression(); static void fie_term(); static void fie_factor(); static void fie_function(); static void fie_error(); static void fie_push(); static fiedbl fie_erf(); static fiedbl fie_erfc(); static fiedbl fie_mod(); static fiedbl fie_int(); static fiedbl fie_nint(); static fiedbl fie_sign(); static double fie_ran(); static fiedbl fie_ranu(); static fiedbl fie_rang(); static fiedbl fie_ranp(); static fiedbl fie_pop(); static fiedbl fie_max(); static fiedbl fie_min(); static fiedbl fie_sum(); static fiedbl fie_mean(); static fiedbl fie_var(); static fiedbl fie_sdev(); static fiedbl fie_adev(); static fiedbl fie_skew(); static fiedbl fie_kurt(); static fiedbl fie_median(); static fiedbl fie_nblanks(); static void fie_nextch() { if (pos < maxpos) ch = in_line[pos++]; else ch = 0; } static void fie_nextsym() { int tenp; double f, frac; int i; char fun[MAXFUNLEN]; while (ch == ' ') fie_nextch(); if (isdigit(ch)||(ch == '.')) { curconst = 0.0; while (isdigit(ch)) { curconst = 10.0 * curconst + (ch - '0'); fie_nextch(); } if (ch == '.') { fie_nextch(); f = 1; frac = 0; while (isdigit(ch)) { frac = 10 * frac + (ch - '0'); f = 10 * f; fie_nextch(); }; curconst = curconst + ((float) frac)/((float) f); }; if ((ch == 'E')||(ch == 'D')||(ch =='e')||(ch == 'd')) { fie_nextch(); tenp = 1; frac = 0; if (ch == '+') { fie_nextch(); } else if (ch == '-') { tenp = -tenp; fie_nextch(); } while (isdigit(ch)) { frac = 10 * frac + (ch - '0'); fie_nextch(); } frac = frac * tenp; curconst = curconst * pow(10.0,(double) frac); } sym = CONST; } else if (isalpha(ch)) { i = 0; while ((isalpha(ch)||isdigit(ch)) && (i < MAXFUNLEN)) { fun[i++] = toupper(ch); fie_nextch(); } fun[i] = 0; sym = FUNCT; if (nparname) { curpar = 0; while ((curpar < nparname) && strcmp(fun,parname[curpar])) curpar++; if (curpar == nparname) curpar = MAXPAR + 1; else curpar++; if (curpar < MAXPAR) { if (curpar > npar) npar = curpar; sym = PARAM; } } if (sym == FUNCT) { curfun = 0; while ((curfun < MAXFUNCTS) && strcmp(fun,functs[curfun])) curfun++; if (curfun == MAXFUNCTS) { fie_error(); } else { sym = FUNCT; } } } else switch(ch) { case '\0' : sym = END; fie_nextch(); break; case '+' : sym = PLUS; fie_nextch(); break; case '-' : sym = MINUS; fie_nextch(); break; case '*' : { sym = TIMES; fie_nextch(); if (ch == '*') { sym = POWER; fie_nextch(); } break; } case '/' : sym = DIVIDE; fie_nextch(); break; case '(' : sym = LPAR; fie_nextch(); break; case ')' : sym = RPAR; fie_nextch(); break; case ',' : sym = COMMA; fie_nextch(); break; case '$' : { if (nparname) { char par[MAXPARNAMELEN+1]; fie_nextch(); i = 0; while ((isalpha(ch)||isdigit(ch)) && (i < MAXPARNAMELEN)) { par[i++] = toupper(ch); fie_nextch(); } par[i] = 0; curpar = 0; while ((curpar < nparname) && strcmp(par,parname[curpar])) curpar++; if (curpar == nparname) curpar = MAXPAR + 1; else curpar++; } else { fie_nextch(); if (!isdigit(ch)) fie_error(); curpar = 0; while (isdigit(ch)) { curpar = curpar * 10 + ch - '0'; fie_nextch(); } } if (curpar > MAXPAR || curpar < 1) fie_error(); if (curpar > npar) npar = curpar; sym = PARAM; break; } default: fie_error(); fie_nextch(); break; } } static void fie_expression() { int s; fie_term(); while ((sym == PLUS) || (sym == MINUS)) { s = sym; fie_nextsym(); fie_term(); if (s == PLUS) fie_gencode(ADD); else fie_gencode(SUB); } } static void fie_term() { int s; fie_factor(); while ((sym == TIMES) || (sym == DIVIDE)) { s = sym; fie_nextsym(); fie_factor(); if (s == TIMES) fie_gencode(MUL); else fie_gencode(DIV); } } static void fie_factor() { switch (sym) { case LPAR : { fie_nextsym(); fie_expression(); if (sym == RPAR) fie_nextsym(); else fie_error(); break; } case PLUS : fie_nextsym(); fie_factor(); break; case MINUS : fie_nextsym(); fie_factor(); fie_gencode(NEG); break; case PARAM : fie_gencode(LDP); fie_gencode(curpar); fie_nextsym(); break; case CONST : fie_genconst(curconst); fie_nextsym(); break; case FUNCT : fie_function(); break; default : fie_error(); break; } if (sym == POWER) { fie_nextsym(); fie_factor(); fie_gencode(PWR); } } static void fie_function() { int f = curfun; int n = nargs[curfun]; int m = 0; fie_nextsym(); if (n > 0) { if (sym == LPAR) fie_nextsym(); else fie_error(); while (n > 0) { fie_expression(); if (--n > 0) { if (sym == COMMA) fie_nextsym(); else fie_error(); } } if (sym == RPAR) fie_nextsym(); else fie_error(); } else if ( n < 0 ) { if (sym == LPAR) fie_nextsym(); else fie_error(); while (sym != RPAR) { m += 1; if ( m > MAXARG ) { fie_error( ); break; } fie_expression(); if (sym == COMMA) { fie_nextsym(); } else if (sym != RPAR) { fie_error( ); break; } } if (sym == RPAR) { fie_nextsym(); if ( m < -n ) fie_error( ); if ( m == 0 ) fie_error( ); } } fie_gencode( FIE + f ); if ( n < 0 ) fie_gencode(m); } static void fie_error() { if (errorpos == 0) errorpos = pos; sym = END; } #define STACKMAX (20*MAXARG) static fiedbl stack[STACKMAX]; static int sp; static double x[MAXARG]; static double work[MAXARG]; static fiedbl getval( fiedbl arg[], fint narg, fint opt ) { fiedbl r; fint nout; double dblank; int i; setdblank_c( &dblank ); for ( i = 0; i < narg; i++ ) { if (!arg[i].u) x[i] = arg[i].d; /* Copy element */ else x[i] = dblank; } r.u = 0; r.d = momsd_c( &opt, x, work, &narg, &nout ); if (r.d == dblank) r.u = 1; return( r ); } static fiedbl fie_sum( fiedbl arg[], int narg ) { return( getval(arg, narg, 1) ); } static fiedbl fie_mean( fiedbl arg[], int narg ) { return( getval(arg, narg, 2) ); } static fiedbl fie_var( fiedbl arg[], int narg ) { return( getval(arg, narg, 3) ); } static fiedbl fie_sdev( fiedbl arg[], int narg ) { return( getval(arg, narg, 4) ); } static fiedbl fie_adev( fiedbl arg[], int narg ) { return( getval(arg, narg, 5) ); } static fiedbl fie_skew( fiedbl arg[], int narg ) { return( getval(arg, narg, 6) ); } static fiedbl fie_kurt( fiedbl arg[], int narg ) { return( getval(arg, narg, 7) ); } static fiedbl fie_median( fiedbl arg[], int narg ) { return( getval(arg, narg, 8) ); } static fiedbl fie_max( fiedbl arg[], int narg ) { return( getval(arg, narg, 9) ); } static fiedbl fie_min( fiedbl arg[], int narg ) { return( getval(arg, narg, 10) ); } static fiedbl fie_nblanks( fiedbl arg[], int narg ) { return( getval(arg, narg, 11) ); } static void fie_push( fiedbl r ) { stack[++sp] = r; } static fiedbl fie_erf( fiedbl arg ) { if (!arg.u) { double p = 0.327591100; double a1 = 0.254829592; double a2 = -0.284496736; double a3 = 1.421413741; double a4 = -1.453152027; double a5 = 1.061405429; double t1 = 1.0 / ( 1.0 + p * fabs(arg.d)); double t2 = t1*t1, t3 = t1*t2, t4 = t1*t3, t5 = t4*t1; if (arg.d > 0.0) { arg.d = (1.0 - (a1*t1+a2*t2+a3*t3+a4*t4+a5*t5)*exp(-arg.d*arg.d)); } else { arg.d = ((a1*t1+a2*t2+a3*t3+a4*t4+a5*t5)*exp(-arg.d*arg.d) - 1.0); } } return( arg ); } static fiedbl fie_erfc( fiedbl arg ) { arg = fie_erf(arg); if (!arg.u) arg.d = 1.0 - arg.d; return(arg); } static fiedbl fie_mod( fiedbl arg1, fiedbl arg2 ) { fiedbl val; if (arg1.u || arg2.u) { val.d = 0.0; val.u = 1; } else if (arg2.d == 0.0) { val.d = 0.0; val.u = 1; } else { int xxx = arg1.d/arg2.d; val.d = arg1.d - xxx * arg2.d; val.u = 0; } return(val); } static fiedbl fie_int( fiedbl arg ) { if (!arg.u) { int xxx = arg.d; arg.d = xxx; } return(arg); } static fiedbl fie_nint( fiedbl arg ) { if (!arg.u) { int xxx = (fabs(arg.d) + 0.5); arg.d = arg.d<0?-xxx:xxx; } return(arg); } static fiedbl fie_sign( fiedbl arg ) { if (!arg.u) { if (arg.d > 0.0) arg.d = 1.0; else if (arg.d < 0.0) arg.d = -1.0; } return(arg); } static double fie_ran() { int xxx = rand(); return( (double) xxx / (double) RAND_MAX ); } static fiedbl fie_ranu( fiedbl arg1, fiedbl arg2 ) { fiedbl val; if (arg1.u || arg2.u) { val.d = 0.0; val.u = 1; } else { val.d = arg1.d + fie_ran() * (arg2.d - arg1.d); val.u = 0; } return(val); } static fiedbl fie_rang( fiedbl arg1, fiedbl arg2 ) { fiedbl val; if (arg1.u || arg2.u) { val.d = 0.0; val.u = 1; } else { double v, r1, r2; r1 = fie_ran(); if (r1 == 0.0) r1 = DBL_EPSILON; /* Avoid log(0) */ r2 = fie_ran(); if (oddran == 0) { v = sqrt(-2*log(r1))*cos(6.283185307179586476925286*r2); oddran = 1; } else { v = sqrt(-2*log(r1))*cos(6.283185307179586476925286*r2); oddran = 0; } val.d = arg1.d + fabs(arg2.d) * v; val.u = 0; } return(val); } static fiedbl fie_ranp( fiedbl arg ) { if (!arg.u) { double val, cum, p, f; if (arg.d < 0.0) { val = 0.0; arg.u = 1; } else if (arg.d < 40) { fiedbl a0, a1, a2; int xxx; a1 = a2 = arg; a2.d = sqrt(a2.d); a0 = fie_rang(a1,a2); xxx = a0.d + 0.5; val = xxx; } else { cum = exp(-arg.d); p = cum; val = 0.0; f = fie_ran(); while ( f >= cum) { val = val + 1.0; p = p * arg.d / val; cum = cum + p; } } arg.d = val; } return(arg); } static fiedbl fie_pop() { return( stack[sp--] ); } /* #> fiepar.dc2 Function: FIEPAR Purpose: Defines the parameter names for the next call to FIEINI. Category: MATH File: fie.c Author: K.G. Begeman Use: INTEGER FIEPAR( PARNAMES , Input CHARACTER*(*) ARRAY NPARS ) Input INTEGER FIEPAR Returns 0 on success, otherwise: -1: Too many parameters, max = 32. PARNAMES The parameter names. Each name has a maximum of 16 characters. NPARS The number of parameter names in PARNAMES. Notes: FIEPAR does not check whether the parameter names are unique. If a parameter name with the same name as a built-in constant or function is used, the former will supersede the latter. Updates: Mar 19, 1991: KGB Document created. Feb 7, 1997: JPT Named parameters supersede built-in constants. jun 26, 1998: VOG Increased MAXFIECODE from 100 to 256 #< Fortran to C interface: @ integer function fiepar( character, integer ) */ fint fiepar_c( fchar parnames , fint *nparnames ) { char *ptr = parnames.a; fint k; fint l = parnames.l; fint m; fint n; if ((*nparnames) > MAXPAR) return( -1 ); nparname = (*nparnames); for (n = 0; n < nparname; n++) { for (k = 0, m = 0; m < l; m++, ptr++) { if (k < MAXPARNAMELEN && (isalpha(*ptr) || isdigit(*ptr))) { parname[n][k++] = toupper(*ptr); } } parname[n][k] = 0; } return( 0 ); } /* #> fieini.dc2 Function: FIEINI Purpose: Decodes a string containing a mathematical expression for FIEDO. Category: MATH File: fie.c Author: K.G. Begeman Use: INTEGER FIEINI( STRING, Input CHARACTER*(*) FIEID, Output INTEGER ERRPOS ) Output INTEGER FIEINI Returns: >0 : largest parameter number in expression. This number of parameters must be present in FIEDO. 0 : no parameter in expression. -1 : syntax error in expression at ERRPOS. -2 : no storage space left. STRING String containing the mathematical expression to be decoded for evaluation by FIEDO. See FIE.DC2 for a detailed description on the syntax and possibilities. FIEID An non negative number which is unique for the evaluated expression. FIEDO needs this id to evaluate the expression. ERRPOS If a syntax error was detected by FIEINI, this number gives the approximate position in the string where the error occurred. Updates: Aug 1, 1989: KGB, original document. #< @ integer function fieini( character, integer, integer ) */ fint fieini_c( fchar expr, fint *id, fint *errpos ) /* * expr string containing the expression * id pointer to codebuffer * errpos position of error in expression */ { int ret; in_line = expr.a; maxpos = expr.l; oddran = 0; error = 0; pos = 0; errorpos = 0; codeptr = 0; npar = 0; opcodeptr = 0; ch = ' '; fie_nextsym(); fie_expression(); if (sym != END) fie_error(); fie_gencode(HLT); if (errorpos == 0) { int nid = 0; while ((nid < MAXID) && (saveifgen[nid])) nid++; if (nid == MAXID) ret = -2; else { int n; for (n = 0; n < MAXFIECODE; n++) savefiecode[nid][n] = fiecode[n]; saveifgen[nid] = 1; savenpars[nid] = npar; *id = nid; ret = npar; } } else { ret = -1; *errpos = errorpos; } return(ret); } /* #> fiedo.dc2 Function: FIEDO Purpose: FIEDO evaluates the code generated by FIEINI. Category: MATH File: fie.c Author: K.G. Begeman Use: INTEGER FIEDO( PARS, Input REAL ARRAY NOPS, Input INTEGER RESULT, Output REAL ARRAY FIEID ) Input INTEGER FIEDO Returns: 0 : successful -1 : id out of range -2 : no code generated for id by FIEINI PARS Parameters as specified in string processed by FIEINI. The first NOPS elements contain parameter 1, the second NOPS elements contain parameter 2 and so on. The number of elements in PARS is equal to NOPS * (the number of parameters as returned by FIEINI). NOPS Number of operations, that is, how many times to repeat the evaluation for different sets of parameters. RESULT The result of each evaluation. The number of elements in RESULT is NOPS. FIEID Id of code generated by FIEINI. Updates: Aug 1, 1989: KGB, original document. #< @ integer function fiedo( real, integer, real, integer ) */ fint fiedo_c( float *data, fint *nop, float *results, fint *id ) /* * data dimensions: * * nop number of operations * results dim. * id id of generated code */ { int iop, i, c, o, opc; fiedbl a, b, r; fiedbl params[MAXPAR], arg[MAXARG]; setfblank_c( &BLANK ); if ((*id >= 0) && (*id < MAXID)) { if (saveifgen[*id]) { int n = 0; for (n = 0; n < MAXFIECODE; n++) fiecode[n] = savefiecode[*id][n]; npar = savenpars[*id]; } else return(-1); } else return(-2); for (iop = 0; iop < *nop; iop++) { c = o = sp = 0; for (i = 0; i < npar; i++) { int q = iop + (*nop) * (i); if (data[q] == BLANK) { params[i].d = 0.0; params[i].u = 1; } else { params[i].d = data[q]; params[i].u = 0; } } do { int narg = 0; opc = fiecode[c].opcode[o++]; if (o == bid) { c++ ; o = 0; } if (opc >= FIE) { int n; narg = nargs[opc-FIE]; if ( narg < 0 ) { narg = fiecode[c].opcode[o++]; if (o == bid) { c++ ; o = 0; } } for (n = 1; n <= narg; n++) arg[narg-n] = fie_pop(); } switch(opc) { case HLT: break; case ADD: { b = fie_pop(); a = fie_pop(); if (a.u || b.u) { r.d = 0.0; r.u = 1; } else { r.d = a.d + b.d; r.u = 0; } fie_push(r); break; } case SUB: { b = fie_pop(); a = fie_pop(); if (a.u || b.u) { r.d = 0.0; r.u = 1; } else { r.d = a.d - b.d; r.u = 0; } fie_push(r); break; } case MUL: { b = fie_pop(); a = fie_pop(); if (a.u || b.u) { r.d = 0.0; r.u = 1; } else { r.d = a.d * b.d; r.u = 0; } fie_push(r); break; } case DIV: { b = fie_pop(); a = fie_pop(); if (a.u || b.u) { r.d = 0.0; r.u = 1; } else if (b.d == 0.0) { r.d = 0.0; r.u = 1; } else { r.d = a.d / b.d; r.u = 0; } fie_push(r); break; } case NEG: { a = fie_pop(); if (!a.u) a.d = -a.d; fie_push(a); break; } case PWR: { b = fie_pop(); a = fie_pop(); if (a.u || b.u) { r.d = 0.0; r.u = 1; } else if ( a.d == 0.0 && b.d < 0 ) { r.d = 0.0; r.u = 1; } else if (a.d >= 0) { r.d = pow(a.d,b.d); r.u = 0; } else { int p = b.d, t; double epsilon = 0.000001; if (fabs(b.d - p) <= epsilon) { t = (p % 2 == 0) ? 1 : -1; r.d = t * pow(fabs(a.d),b.d); r.u = 0; } else { r.d = 0.0; r.u = 1; } } fie_push(r); break; } case LDP: { opc = fiecode[c].opcode[o++]; if (o == bid) { c++; o = 0; } fie_push( params[opc-1] ); break; } case LDC: { if (o != 0) { c++; o = 0; } r.d = fiecode[c++].c; r.u = 0; fie_push(r); break; } default: switch(opc-FIE) { case 0: { a = arg[0]; if (!a.u) a.d = sin(a.d); fie_push(a); break; } case 1: { a = arg[0]; if (!a.u) { if (fabs(a.d) > 1) { a.d = 0.0; a.u = 1; } else { a.d = asin(a.d); } } fie_push(a); break; } case 2: { a = arg[0]; if (!a.u) { if (fabs(a.d) > 70) { a.d = 0.0; a.u = 1; } else { a.d = sinh(a.d); } } fie_push(a); break; } case 3: { a = arg[0]; if (!a.u) a.d = cos(a.d); fie_push(a); break; } case 4: { a = arg[0]; if (!a.u) { if (fabs(a.d) > 1) { a.d = 0.0; a.u = 1; } else { a.d = acos(a.d); } } fie_push(a); break; } case 5: { a = arg[0]; if (!a.u) { if (fabs(a.d) > 70) { a.d = 0.0; a.u = 1; } else { a.d = cosh(a.d); } } fie_push(a); break; } case 6: { a = arg[0]; if (!a.u) a.d = tan(a.d); fie_push(a); break; } case 7: { a = arg[0]; if (!a.u) a.d = atan(a.d); fie_push(a); break; } case 8: { a = arg[0]; if (!a.u) { if (fabs(a.d) > 70) { a.d = a.d<0?-1.0:1.0; } else { a.d = tanh(a.d); } } fie_push(a); break; } case 9: { a = arg[0]; b = arg[1]; if (a.u || b.u) { r.d = 0.0; r.u = 1; } else { r.d = atan2(a.d,b.d); r.u = 0; } fie_push(r); break; } case 10: { a = arg[0]; if (!a.u) a.d *= 0.017453292519943295769237; fie_push(a); break; } case 11: { a = arg[0]; if (!a.u) a.d *= 57.295779513082320876798155; fie_push(a); break; } case 12: { r.d = 3.141592653589793238462643; r.u = 0; fie_push(r); break; } case 13: { a = arg[0]; if (!a.u) { if (a.d > 70) { a.d = 0.0; a.u = 1; } else if (a.d < -70) { a.d = 0.0; } else { a.d = exp(a.d); } } fie_push(a); break; } case 14: { a = arg[0]; if (!a.u) { if (a.d > 0) { a.d = log(a.d); } else { a.d = 0.0; a.u = 1; } } fie_push(a); break; } case 15: { a = arg[0]; if (!a.u) { if (a.d > 0) { a.d = log10(a.d); } else { a.d = 0.0; a.u = 1; } } fie_push(a); break; } case 16: { a = arg[0]; if (!a.u) { if (a.d < 0) { a.d = 0.0; a.u = 1; } else { a.d = sqrt(a.d); } } fie_push(a); break; } case 17: { a = arg[0]; if (!a.u) a.d = fabs(a.d); fie_push(a); break; } case 18: { a = arg[0]; if (!a.u) { if (fabs(a.d) < 1.0e-30) a.d = 1.0; else a.d = sin(a.d)/a.d; } fie_push(a); break; } case 19: { r.d = 2.997925e+8; r.u = 0; fie_push(r); break; } case 20: { r.d = 6.6732e-11; r.u = 0; fie_push(r); break; } case 21: { r.d = 1.99e30; r.u = 0; fie_push(r); break; } case 22: fie_push(fie_erf(arg[0])); break; case 23: fie_push(fie_erfc(arg[0])); break; case 24: { r.d = 1.380622e-23; r.u = 0; fie_push(r); break; } case 25: { r.d = 6.6256196e-34; r.u = 0; fie_push(r); break; } case 26: { r.d = 3.086e16; r.u = 0; fie_push(r); break; } case 27: { r.d = 5.66961e-8; r.u =0; fie_push(r); break; } case 28: fie_push(fie_max(arg,narg)); break; case 29: fie_push(fie_min(arg,narg)); break; case 30: fie_push(fie_mod(arg[0],arg[1])); break; case 31: fie_push(fie_int(arg[0])); break; case 32: fie_push(fie_nint(arg[0])); break; case 33: fie_push(fie_sign(arg[0])); break; case 34: { r.d = 0.0; r.u = 1; fie_push(r); break; } case 35: { a = arg[0]; b = arg[1]; if (a.u || b.u) { r.d = 0.0; r.u = 1; } else if (a.d > b.d) r = arg[2]; else r = arg[3]; fie_push(r); break; } case 36: { a = arg[0]; b = arg[1]; if (a.u || b.u) { r.d = 0.0; r.u = 1; } else if (a.d < b.d) r = arg[2]; else r = arg[3]; fie_push(r); break; } case 37: { a = arg[0]; b = arg[1]; if (a.u || b.u) { r.d = 0.0; r.u = 1; } else if (a.d >= b.d) r = arg[2]; else r = arg[3]; fie_push(r); break; } case 38: { a = arg[0]; b = arg[1]; if (a.u || b.u) { r.d = 0.0; r.u = 1; } else if (a.d <= b.d) r = arg[2]; else r = arg[3]; fie_push(r); break; } case 39: { a = arg[0]; b = arg[1]; if (a.u && b.u) { r = arg[2]; } else if (a.u || b.u) { r = arg[3]; } else if (a.d == b.d) { r = arg[2]; } else { r = arg[3]; } fie_push(r); break; } case 40: { a = arg[0]; b = arg[1]; if (a.u && b.u) { r = arg[3]; } else if (a.u || b.u) { r = arg[2]; } else if (a.d != b.d) { r = arg[2]; } else { r = arg[3]; } fie_push(r); break; } case 41: fie_push(fie_ranu(arg[0],arg[1])); break; case 42: fie_push(fie_rang(arg[0],arg[1])); break; case 43: fie_push(fie_ranp(arg[0])); break; case 44: { a = arg[0]; if (a.u) r = arg[1]; else r = arg[2]; fie_push(r); break; } case 45: fie_push(fie_sum(arg,narg)); break; case 46: fie_push(fie_mean(arg,narg)); break; case 47: fie_push(fie_var(arg,narg)); break; case 48: fie_push(fie_sdev(arg,narg)); break; case 49: fie_push(fie_adev(arg,narg)); break; case 50: fie_push(fie_skew(arg,narg)); break; case 51: fie_push(fie_kurt(arg,narg)); break; case 52: fie_push(fie_median(arg,narg)); break; case 53: fie_push(fie_nblanks(arg,narg)); break; case 54: { a = arg[0]; b = arg[1]; if (a.u || b.u) { r.d = 0.0; r.u = 1; } else { r.d = bessj((int)(a.d+0.5),b.d); r.u = 0; } fie_push(r); break; } case 55: { a = arg[0]; b = arg[1]; if (a.u || b.u || b.d<1.0e-300) { r.d = 0.0; r.u = 1; } else { r.d = bessy((int)(a.d+0.5),b.d); r.u = 0; } fie_push(r); break; } case 56: { a = arg[0]; b = arg[1]; if (a.u || b.u || b.d<1.0e-300) { r.d = 0.0; r.u = 1; } else { r.d = bessk((int)(a.d+0.5),b.d); r.u = 0; } fie_push(r); break; } case 57: { a = arg[0]; b = arg[1]; if (a.u || b.u || b.d>90) { r.d = 0.0; r.u = 1; } else { r.d = bessi((int)(a.d+0.5),b.d); r.u = 0; } fie_push(r); break; } case 58: { /* SEC */ a = arg[0]; if (!a.u) { double q=cos(a.d); if (q != 0.0) { a.d = 1.0/q; } else { a.d = 0.0; a.u = 1; } } fie_push(a); break; } case 59: { /* CSC */ a = arg[0]; if (!a.u) { double q=sin(a.d); if (q != 0.0) { a.d = 1.0/q; } else { a.d = 0.0; a.u = 1; } } fie_push(a); break; } case 60: { /* COT */ a = arg[0]; if (!a.u) { double q=tan(a.d); if (q != 0.0) { a.d = 1.0/q; } else { a.d = 0.0; a.u = 1; } } fie_push(a); break; } case 61: { /* SECH */ a = arg[0]; if (!a.u) { if (fabs(a.d) > 70) { a.d = 0.0; } else { a.d = 1.0/cosh(a.d); } } fie_push(a); break; } case 62: { /* CSCH */ a = arg[0]; if (!a.u) { if (fabs(a.d) > 70) { a.d = 0.0; } else { double q=sinh(a.d); if (q != 0.0) { a.d = 1.0/q; } else { a.d = 0.0; a.u = 1; } } } fie_push(a); break; } case 63: { /* COTH */ a = arg[0]; if (!a.u) { if (fabs(a.d) > 70) { a.d = a.d<0.0?-1.0:1.0; } else { double q=tanh(a.d); if (q != 0.0) { a.d = 1.0/q; } else { a.d = 0.0; a.u = 1; } } } fie_push(a); break; } case 64: { /* ASINH */ a = arg[0]; if (!a.u) a.d = log(a.d+sqrt(a.d*a.d+1.0)); fie_push(a); break; } case 65: { /* ACOSH */ a = arg[0]; if (!a.u) { if (a.d < 1.0) { a.d = 0.0; a.u = 1; } else { a.d = log(a.d+sqrt(a.d*a.d-1.0)); } } fie_push(a); break; } case 66: { /* ATANH */ a = arg[0]; if (!a.u) { if ((a.d < 1.0) && (a.d > -1.0) && (a.d != 0.0)) { a.d = 0.5*log((1+a.d)/(1-a.d)); } else { a.d = 0.0; a.u = 1; } } fie_push(a); break; } case 67: { /* STEP */ a = arg[0]; if (!a.u) a.d = a.d>0.0?1.0:0.0; fie_push(a); break; } case 68: { /* RECT */ a = arg[0]; if (!a.u) a.d = fabs(a.d)<0.5?1.0:0.0; fie_push(a); break; } default: break; } break; } } while (opc != HLT); r = fie_pop(); if (r.u) results[iop] = BLANK; else results[iop] = r.d; } return(0); } /* #> fieclr.dc2 Function: FIECLR Purpose: Clears code previous generated by FIEINI. Category: MATH File: fie.c Author: K.G. Begeman Use: INTEGER FIECLR( FIEID ) Input INTEGER FIECLR Returns: 0 : successful -1 : id out of range -2 : no code generated for this id FIEID Id of code generated by FIEINI which will be cleared. Updates: Mar 7, 1994: KGB, original document. #< @ integer function fieclr( integer ) */ fint fieclr_c( fint *id ) /* * id id of generated code */ { if ( (*id >= 0) && (*id < MAXID)) { if (saveifgen[*id]) { int n; saveifgen[*id] = 0; savenpars[*id] = 0; for (n = 0; n < MAXFIECODE; n++) { int m; for (m = 0; m < bid; m++ ) { savefiecode[*id][n].opcode[m] = 0; } } } else return( -2 ); } else return( -1 ); return( 0 ); } /* #> fiedump.dc2 Function: FIEDUMP Purpose: Display the code generated by FININI on stdin. Category: MATH File: fie.c Author: K.G. Begeman Use: INTEGER FIEDUMP( FIEID ) Input INTEGER FIEDUMP Returns: 0 : successful -1 : id out of range. -2 : no code generated for this id. FIEID Id of code generated by FIEINI. Updates: Aug 1, 1989: KGB, original document. #< @ integer function fiedump( integer ) */ fint fiedump_c( fint *id ) /* * id pointer to codebuffer */ { int c, o, opc, op; if ((*id >= 0) && (*id < MAXID)) { if (saveifgen[*id]) { int n = 0; for (n = 0; n < MAXFIECODE; n++) fiecode[n] = savefiecode[*id][n]; npar = savenpars[*id]; } else return( (fint) -1 ); } else return( (fint) -2 ); c = o = 0; do { op = fiecode[c].opcode[o++]; if (o == bid) { c++ ; o = 0; }; opc = op>FIE ? FIE : op; printf(" %s",mnem[opc]); if (opc == LDP) { opc = fiecode[c].opcode[o++]; if (o == bid) { c++ ; o = 0; }; printf(" %d",opc++); } else if (opc == FIE) { printf(" %s",functs[op-opc]); if ( nargs[op-opc] < 0 ) { printf( " [%d]", fiecode[c].opcode[o++] ); if (o == bid) { c++ ; o = 0; }; } } else if (opc == LDC) { if (o != 0) c++; printf(" %f",fiecode[c++].c); o = 0; } printf("\n"); } while (opc != HLT); return( (fint) 0 ); } #if defined(TESTBED) int main() { #if 1 char *parb = "x y z "; fchar pars; fint three = 3; char *fie[] = {"sdev(x,y,z)", "x+y/z" }; #else char *fie[] = { "$1+$2/$3", "1.0+exp($1)*sinc($3)" }; #endif fint id[2]; int j; float par[30], res[10], *p1, *p2, *p3, *r; fint i, n = 10, nret, epos; #if 1 pars.a = parb; pars.l = 2; nret = fiepar_c( pars, &three ); printf( "fiepar = %d\n", nret ); #endif for (i = 0; i < 30; i++) { par[i] = (float) i; }; setfblank_c(&par[3]); setfblank_c(&par[14]); setfblank_c(&par[25]); setfblank_c(&par[6]); setfblank_c(&par[16]); setfblank_c(&par[26]); for ( j = 0; j < 2; j++ ) { if (j>0) { nret = fieclr_c( &id[0] ); printf("fieclr : %d\n", nret ); } nret = fieini_c( tofchar(fie[j]), &id[j], &epos ); printf("fieini : %d id=%d\n",nret,id[j]); if (nret < 0) { if (nret == -1) { printf("string : %s\n", fie[j] ); printf("error at: %d\n", epos ); } exit(1); } nret = fiedump_c(&id[j]); printf("fiedump: %d\n", nret ); if (nret < 0) { exit(1); } nret = fiedo_c( par, &n, res, &id[j] ); printf("fiedo : %d\n", nret ); if (nret < 0) { exit(1); } p1 = &par[0]; p2 = &par[10]; p3 = &par[20]; r = &res[0]; for (i = 0; i < n; i++) { if (fblank_c(p1)) printf(" BLANK "); else printf("%10f ", *p1); p1++; if (fblank_c(p2)) printf(" BLANK "); else printf("%10f ", *p2); p2++; if (fblank_c(p3)) printf(" BLANK "); else printf("%10f ", *p3); p3++; if (fblank_c(r)) printf(" BLANK \n"); else printf("%10f\n", *r); r++; } } return( 0 ); } #endif