SUBROUTINE PARSER( C_EXPR , L_PRINT , NUM_CODE , C_CODE ) C C Parse the arithmetic expression in C_EXPR. The code required to C evaluate the expression is returned in the first NUM_CODE entries C of the CHARACTER*8 array C_CODE. If NUM_CODE is returned as zero, C an error occurred. On input, L_PRINT determines whether or not to C print error messages. C C Modified 02/17/89 by RWCox from APEVAL subroutine in APFORT, for PC. C Modified 06/29/89 by RWCox for Sun Fortran. C Modified 04/04/91 by RWCox to fix problem with -x**2 type operations. C Modified 11/20/96 by RWCox to try to control errors in evaluation. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IMPLICIT NONE C CHARACTER*(*) C_EXPR LOGICAL L_PRINT INTEGER NUM_CODE CHARACTER*8 C_CODE(*) C C Compilation, evaluation, and function stacks. C INTEGER NUM_CSTACK , NUM_FSTACK PARAMETER ( NUM_CSTACK = 2048 , NUM_FSTACK = 40 ) C INTEGER N_CODE(NUM_CSTACK) , NCODE , NEXTCODE , X N_FUNC(NUM_FSTACK) , NFUNC C C Random local stuff C INTEGER NLEN , NPOS , IPOS , NTOKEN , NUSED , NF , NARG , X NERR CHARACTER*10000 C_LOCAL CHARACTER*30 C_MESSAGE CHARACTER*1 C_CH REAL*8 VAL_TOKEN , R8_TOKEN CHARACTER*8 C8_TOKEN EQUIVALENCE ( R8_TOKEN , C8_TOKEN ) C INTEGER LC_TO_UC , LAST_NONBLANK C PARAMETER ( LC_TO_UC = ICHAR('A') - ICHAR('a') ) C INCLUDE 'parser.inc' C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C NLEN = LAST_NONBLANK( C_EXPR ) IF( NLEN .LE. 0 .OR. NLEN .GT. 9999 )THEN C!no input, or too much NUM_CODE = 0 GOTO 8000 ENDIF C C Copy input string to local, deleting blanks and converting case. C NPOS = 0 DO 10 IPOS=1,NLEN C_CH = C_EXPR(IPOS:IPOS) IF( C_CH .NE. ' ' )THEN IF( C_CH .GE. 'a' .AND. C_CH .LE. 'z' ) X C_CH = CHAR( ICHAR(C_CH) + LC_TO_UC ) C!convert case NPOS = NPOS+1 C_LOCAL(NPOS:NPOS) = C_CH ENDIF 10 CONTINUE C!tack 1 blank at the end NLEN = NPOS+1 C_LOCAL(NLEN:NLEN) = ' ' C....................................................................... C This routine parses expressions according to the grammar: C C EXPR == E9 E8 E6 E4 $ C C E4 == E9 E8 E6 E4 | C C E6 == E9 E8 E6 | C C E8 == E9 E8 | C C E9 == | ( E9 E8 E6 E4 ARGTL ) C | ( E9 E8 E6 E4 ) | E9 C C ARGTL == , E9 E8 E6 E4 ARGTL | C C is + or - C is * or / C is ** C is a literal number or a DCL variable C is in the list C_FUNCNAME C C The predictive parser described in Aho and Ullman, "Principles of C Compiler Design" on pages 185-191 for LL(1) grammars is used here, C with additions to perform the evaluation as the parsing proceeds. C These consist of adding code (NC_) to the compilation stack when an C expansion is made. When the code is popped off the stack, it is C executed. C C 02/17/89: Now, when code is popped off the stack, it is just C added to the output code list. C....................................................................... C Prepare to process input string. Initialize the stacks, etc. C C!start scan at 1st character NPOS = 1 C!no function calls yet NFUNC = 0 C!initial compile stack is E9 E8 E6 E4 $ N_CODE(1) = NN_END N_CODE(2) = NN_E4 N_CODE(3) = NN_E6 N_CODE(4) = NN_E8 N_CODE(5) = NN_E9 NCODE = 5 NUM_CODE = 0 C....................................................................... C 1000 is the loop back point to process the next token in the input C string. C 1000 CONTINUE CALL GET_TOKEN( C_LOCAL(NPOS:NLEN) , NTOKEN , VAL_TOKEN , NUSED ) C IF( NTOKEN .EQ. NT_ERROR )THEN NERR = 1 C_MESSAGE = 'Can''t interpret symbol' GOTO 9000 C!error exit ENDIF C C At 2000, process the next compile code until the token is used up. C 2000 CONTINUE NEXTCODE = N_CODE(NCODE) C C If next entry on the compile stack is an opcode, then apply it to C the evaluation stack. C 02/17/89: just add it to the output C IF( NEXTCODE .GE. NC_BEG .AND. NEXTCODE .LE. NC_STOP )THEN NUM_CODE = NUM_CODE + 1 CALL EXECUTE( NEXTCODE , C_CODE(NUM_CODE) ) NCODE = NCODE - 1 C!remove opcode from compile stack GOTO 2000 C!loop back for next compile stack entry ENDIF C C If next compile stack entry is a token itself, it must match the C new token from the input. C IF( NEXTCODE .GE. NT_BEG .AND. NEXTCODE .LE. NT_STOP )THEN IF( NEXTCODE .EQ. NTOKEN )THEN C!a match NCODE = NCODE - 1 C!pop token from compile stack GOTO 5000 C!loop back for next token ENDIF NERR = 2 IF( NEXTCODE .EQ. NT_OPEN )THEN C_CH = '(' ELSEIF( NEXTCODE .EQ. NT_CLOSE )THEN C_CH = ')' ELSEIF( NEXTCODE .EQ. NT_COMMA )THEN C_CH = ',' ELSE C_CH = '?' ENDIF C_MESSAGE = 'Expected a "' // C_CH // '"' GOTO 9000 C!error exit ENDIF C C Should have a nonterminal (NN) here. C IF( NEXTCODE .LT. NN_BEG .OR. NEXTCODE .GT. NN_STOP )THEN NERR = 3 C_MESSAGE = 'Internal parser error' GOTO 9000 C!error exit ENDIF C C Expand the nonterminal appropriately, depending on the token. C If no legal expansion, then stop with an error. C C TOKEN = end of string C IF( NTOKEN .EQ. NT_EOS )THEN IF( NEXTCODE .EQ. NN_END )THEN C!end of string = end of expr ==> compilation done GOTO 8000 C ELSEIF( NEXTCODE .EQ. NN_E8 .OR. NEXTCODE .EQ. NN_E6 .OR. X NEXTCODE .EQ. NN_E4 )THEN NCODE = NCODE - 1 C!expand this to nothing GOTO 2000 C!and try this token again ENDIF NERR = 4 C_MESSAGE = 'Unexpected end of input' GOTO 9000 C!error exit ENDIF C C Check if end of input was expected but not encountered. C IF( NEXTCODE .EQ. NN_END )THEN NERR = 15 C_MESSAGE = 'Expected end of input' GOTO 9000 C!error exit ENDIF C C TOKEN = number or symbol C 02/17/89: added NT_SYMBOL token type; now, the code for C pushing the number or symbol onto the stack is C added to the output. C IF( NTOKEN .EQ. NT_NUMBER .OR. NTOKEN .EQ. NT_SYMBOL )THEN IF( NEXTCODE .EQ. NN_E9 )THEN C!only legal time for a number IF( NTOKEN .EQ. NT_NUMBER )THEN C_CODE(NUM_CODE+1) = 'PUSHNUM' ELSE C_CODE(NUM_CODE+1) = 'PUSHSYM' ENDIF R8_TOKEN = VAL_TOKEN C_CODE(NUM_CODE+2) = C8_TOKEN NUM_CODE = NUM_CODE + 2 NCODE = NCODE - 1 C!pop E9 from compile stack GOTO 5000 C!go to next token ENDIF NERR = 5 C_MESSAGE = 'Expected an operator' GOTO 9000 C!error exit ENDIF C C TOKEN = function call C IF( NTOKEN .EQ. NT_FUNC )THEN IF( NEXTCODE .EQ. NN_E9 )THEN C!only legal time for a function C N_CODE(NCODE+7) = NT_OPEN C!expand E9 into ( E9 E8 E6 E4 ARGTL ) N_CODE(NCODE+6) = NN_E9 N_CODE(NCODE+5) = NN_E8 N_CODE(NCODE+4) = NN_E6 N_CODE(NCODE+3) = NN_E4 N_CODE(NCODE+2) = NN_ARGTL N_CODE(NCODE+1) = NT_CLOSE N_CODE(NCODE) = NC_FUNC + INT(VAL_TOKEN) NCODE = NCODE + 7 C NFUNC = NFUNC + 2 C!setup function stack to check # arguments N_FUNC(NFUNC-1) = INT(VAL_TOKEN) N_FUNC(NFUNC) = 0 GOTO 5000 C!process next token ENDIF NERR = 6 C_MESSAGE = 'Expected an operator' GOTO 9000 C!error exit ENDIF C C TOKEN = addition operator C IF( NTOKEN .EQ. NT_ADDOP )THEN IF( NEXTCODE .EQ. NN_E4 )THEN C!expand E4 into E9 E8 E6 E4 N_CODE(NCODE+4) = NN_E9 N_CODE(NCODE+3) = NN_E8 N_CODE(NCODE+2) = NN_E6 IF( VAL_TOKEN .EQ. VT_PLUS )THEN N_CODE(NCODE+1) = NC_ADD ELSE N_CODE(NCODE+1) = NC_SUB ENDIF N_CODE(NCODE) = NN_E4 NCODE = NCODE + 4 GOTO 5000 C!process next token C ELSEIF( NEXTCODE .EQ. NN_E6 .OR. NEXTCODE .EQ. NN_E8 )THEN NCODE = NCODE - 1 C!expand E6 or E8 to null and try again GOTO 2000 C ELSEIF( NEXTCODE .EQ. NN_E9 )THEN C!unary + or - IF( VAL_TOKEN .EQ. VT_MINUS )THEN C!expand E9 into E9 E8 if addop is - otherwise leave E9 alone C [04/04/91 change: C used to expand to E9 , which makes -x**2 become (-x)**2] N_CODE(NCODE+2) = NN_E9 N_CODE(NCODE+1) = NN_E8 N_CODE(NCODE) = NC_MINUS NCODE = NCODE + 2 ENDIF GOTO 5000 C!process next token ENDIF NERR = 7 C_MESSAGE = 'Illegal arithmetic syntax' GOTO 9000 C!error exit ENDIF C C TOKEN = multiplying operator C IF( NTOKEN .EQ. NT_MULOP )THEN IF( NEXTCODE .EQ. NN_E6 )THEN C!expand E6 into E9 E8 E6 N_CODE(NCODE+3) = NN_E9 N_CODE(NCODE+2) = NN_E8 IF( VAL_TOKEN .EQ. VT_STAR )THEN N_CODE(NCODE+1) = NC_MUL ELSE N_CODE(NCODE+1) = NC_DIV ENDIF N_CODE(NCODE) = NN_E6 NCODE = NCODE + 3 GOTO 5000 C!next token C ELSEIF( NEXTCODE .EQ. NN_E8 )THEN C!expand E8 to null and try this token again NCODE = NCODE - 1 GOTO 2000 ENDIF NERR = 8 C_MESSAGE = 'Illegal arithmetic syntax' GOTO 9000 C!error exit ENDIF C C TOKEN = exponentiation operator C IF( NTOKEN .EQ. NT_EXPOP )THEN IF( NEXTCODE .EQ. NN_E8 )THEN C!expand E8 into E9 E8 <**> N_CODE(NCODE+2) = NN_E9 N_CODE(NCODE+1) = NN_E8 N_CODE(NCODE) = NC_POWER NCODE = NCODE + 2 GOTO 5000 C!process next token ENDIF NERR = 9 C_MESSAGE = 'Illegal arithmetic syntax' GOTO 9000 C!error exit ENDIF C C TOKEN = comma C IF( NTOKEN .EQ. NT_COMMA )THEN IF( NEXTCODE .EQ. NN_E4 .OR. NEXTCODE .EQ. NN_E6 .OR. X NEXTCODE .EQ. NN_E8 )THEN C NCODE = NCODE - 1 C!pop this nonterminal and try this token again GOTO 2000 C ELSEIF( NEXTCODE .EQ. NN_ARGTL )THEN C!expand ARGTL into E9 E8 E6 E4 ARGTL N_CODE(NCODE+4) = NN_E9 N_CODE(NCODE+3) = NN_E8 N_CODE(NCODE+2) = NN_E6 N_CODE(NCODE+1) = NN_E4 N_CODE(NCODE) = NN_ARGTL NCODE = NCODE + 4 C!add 1 to no. of args. encountered, and check if there are too many N_FUNC(NFUNC) = N_FUNC(NFUNC) + 1 NF = N_FUNC(NFUNC-1) IF( N_FUNCARGS(NF) .LE. N_FUNC(NFUNC) .AND. X N_FUNCARGS(NF) .GT. 0 )THEN NERR = 12 C_MESSAGE = 'Wrong number of arguments' GOTO 9000 C!error exit ENDIF GOTO 5000 C!process next token ENDIF NERR = 10 C_MESSAGE = 'Expected an expression' GOTO 9000 C!error exit ENDIF C C TOKEN = open parenthesis C IF( NTOKEN .EQ. NT_OPEN )THEN IF( NEXTCODE .EQ. NN_E9 )THEN C!expand E9 into E9 E8 E6 E4 ) N_CODE(NCODE+4) = NN_E9 N_CODE(NCODE+3) = NN_E8 N_CODE(NCODE+2) = NN_E6 N_CODE(NCODE+1) = NN_E4 N_CODE(NCODE) = NT_CLOSE NCODE = NCODE + 4 GOTO 5000 C!process next token ENDIF NERR = 11 C_MESSAGE = 'Expected an operator' GOTO 9000 C!error exit ENDIF C C TOKEN = close parenthesis C IF( NTOKEN .EQ. NT_CLOSE )THEN IF( NEXTCODE .EQ. NN_E4 .OR. NEXTCODE .EQ. NN_E6 .OR. X NEXTCODE .EQ. NN_E8 )THEN C NCODE = NCODE - 1 C!pop this nonterminal and try this token out on the next one GOTO 2000 C ELSEIF( NEXTCODE .EQ. NN_ARGTL )THEN C!end of function call C NARG = N_FUNC(NFUNC) + 1 C!check # arguments NF = N_FUNC(NFUNC-1) NFUNC = NFUNC - 2 IF( N_FUNCARGS(NF) .LE. 0 )THEN C!variable # of args ==> push number of args on stack (Feb 1997) C_CODE(NUM_CODE+1) = 'PUSHNUM' R8_TOKEN = NARG C_CODE(NUM_CODE+2) = C8_TOKEN NUM_CODE = NUM_CODE + 2 ELSEIF( N_FUNCARGS(NF) .NE. NARG )THEN C!illegal # of args NERR = 12 C_MESSAGE = 'Wrong number of arguments' GOTO 9000 C!error exit ENDIF C NCODE = NCODE - 1 C!pop this nonterminal and try to match the ) with the next compile stack entry GOTO 2000 ENDIF NERR = 13 C_MESSAGE = 'Expected an expression' GOTO 9000 C!error exit ENDIF NERR = 14 C_MESSAGE = 'Internal parser error' GOTO 9000 C!error exit C....................................................................... C At 5000, advance to the next token and loop back C 5000 CONTINUE NPOS = NPOS + NUSED GOTO 1000 C....................................................................... C At 8000, exit C 8000 CONTINUE RETURN C....................................................................... C At 9000, error exit C 9000 CONTINUE IF( L_PRINT )THEN IF( NUSED .LT. 1 ) NUSED = 1 WRITE(*,9001) NERR , C_MESSAGE , C_LOCAL(1:NLEN) , X (' ',NF=1,NPOS) , ('#',NF=1,NUSED) 9001 FORMAT(' PARSER error',I4,': ',A/ 1X,A/ 80A1 ) C CCC WRITE(*,9002) (N_CODE(NF),NF=NCODE,1,-1) CCC9002 FORMAT(' Compile stack is (top down)' / 10(1X,I6) ) ENDIF C NUM_CODE = 0 RETURN END C C C SUBROUTINE EXECUTE( N_OPCODE , C_CODE ) C C Execute the opcode on the evaluation stack. Note that no attempt is C made to intercept errors, such as divide by zero, ACOS(2), etc. C IMPLICIT NONE INTEGER N_OPCODE CHARACTER*8 C_CODE C INCLUDE 'parser.inc' C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Branch to special code for function evaluations C IF( N_OPCODE .GE. NC_FUNC )GOTO 5000 C....................................................................... IF( N_OPCODE .EQ. NC_MINUS )THEN C!unary minus, the only unary op. C_CODE = '--' C ELSE C!a binary operation IF( N_OPCODE .EQ. NC_ADD )THEN C!add C_CODE = '+' ELSEIF( N_OPCODE .EQ. NC_SUB )THEN C!subtract C_CODE = '-' ELSEIF( N_OPCODE .EQ. NC_MUL )THEN C!multiply C_CODE = '*' ELSEIF( N_OPCODE .EQ. NC_DIV )THEN C!divide C_CODE = '/' ELSEIF( N_OPCODE .EQ. NC_POWER )THEN C!** C_CODE = '**' ENDIF ENDIF GOTO 8000 C....................................................................... C Function evaluation C 5000 CONTINUE C_CODE = C_FUNCNAME(N_OPCODE-NC_FUNC) C....................................................................... 8000 CONTINUE RETURN END C C C SUBROUTINE GET_TOKEN( C_INPUT , NTYPE , VALUE , NUSED ) C C Return the 1st token in the input stream. C IMPLICIT NONE CHARACTER*(*) C_INPUT INTEGER NTYPE , NUSED REAL*8 VALUE C INTEGER NLEN , NPOS , IFUNC , IPOS , IO_CODE LOGICAL L_ALPHABETIC , L_IDCHAR , L_NUMERIC CHARACTER*32 C_ID , C_VAL CHARACTER*1 C_C , C_FIRST C REAL*8 R8_VAL CHARACTER*8 C8_VAL EQUIVALENCE ( R8_VAL , C8_VAL ) C INCLUDE 'parser.inc' C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Statement function definitions C L_ALPHABETIC(C_C) = C_C .GE. 'A' .AND. C_C .LE. 'Z' C L_IDCHAR(C_C) = ( C_C .GE. 'A' .AND. C_C .LE. 'Z' ) .OR. X ( C_C .GE. '0' .AND. C_C .LE. '9' ) .OR. X C_C .EQ. '_' .OR. X C_C .EQ. '$' C L_NUMERIC(C_C) = C_C .GE. '0' .AND. C_C .LE. '9' C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C NTYPE = NT_EOS NUSED = 0 NLEN = LEN(C_INPUT) IF( NLEN .LE. 0 )GOTO 8000 C C Process the simple cases 1st C C_FIRST = C_INPUT(1:1) C IF( C_FIRST .EQ. ' ' )GOTO 8000 C NUSED = 1 IF( C_FIRST .EQ. '+' )THEN NTYPE = NT_ADDOP VALUE = VT_PLUS ELSEIF( C_FIRST .EQ. '-' )THEN NTYPE = NT_ADDOP VALUE = VT_MINUS ELSEIF( C_FIRST .EQ. '/' )THEN NTYPE = NT_MULOP VALUE = VT_SLASH ELSEIF( C_FIRST .EQ. '*' )THEN IF( C_INPUT(1:2) .EQ. '**' )THEN NTYPE = NT_EXPOP VALUE = VT_STARS NUSED = 2 ELSE NTYPE = NT_MULOP VALUE = VT_STAR ENDIF ELSEIF( C_FIRST .EQ. '^' )THEN NTYPE = NT_EXPOP VALUE = VT_STARS ELSEIF( C_FIRST .EQ. '(' )THEN NTYPE = NT_OPEN ELSEIF( C_FIRST .EQ. ')' )THEN NTYPE = NT_CLOSE ELSEIF( C_FIRST .EQ. ',' )THEN NTYPE = NT_COMMA ENDIF C IF( NTYPE .NE. NT_EOS )GOTO 8000 C!exit if above was successful C....................................................................... C The only possibilities left are a variable name, a function name, C or a number. C IF( L_ALPHABETIC(C_FIRST) )THEN C!a name C NPOS = 2 110 CONTINUE IF( .NOT. L_IDCHAR(C_INPUT(NPOS:NPOS)) )GOTO 120 NPOS = NPOS+1 GOTO 110 120 CONTINUE NPOS = NPOS-1 C_ID = C_INPUT(1:NPOS) C C The name is now in C_ID. Check to see if it is a function name. C IFUNC = 1 C_FUNCNAME(NUM_FUNC+1) = C_ID 210 CONTINUE IF( .NOT. ( C_ID .NE. C_FUNCNAME(IFUNC) ) )GOTO 220 IFUNC = IFUNC + 1 GOTO 210 220 CONTINUE IF( IFUNC .LE. NUM_FUNC )THEN C!it is a function NTYPE = NT_FUNC VALUE = IFUNC NUSED = NPOS ELSEIF( C_ID(1:NPOS) .EQ. 'PI' )THEN C!symbolic pi NTYPE = NT_NUMBER VALUE = 3.1415 92653 58979 32D+00 NUSED = NPOS ELSE C!must be a symbol NTYPE = NT_SYMBOL C8_VAL = C_ID(1:NPOS) VALUE = R8_VAL NUSED = NPOS ENDIF C....................................................................... C try for a number C ELSEIF( L_NUMERIC(C_FIRST) .OR. C_FIRST .EQ. '.' )THEN NPOS = 2 310 CONTINUE IF( .NOT. L_NUMERIC(C_INPUT(NPOS:NPOS)) )GOTO 320 C!skip digits NPOS = NPOS+1 GOTO 310 320 CONTINUE IF( C_FIRST .NE. '.' .AND. C_INPUT(NPOS:NPOS) .EQ. '.' )THEN NPOS = NPOS+1 410 CONTINUE IF( .NOT. L_NUMERIC(C_INPUT(NPOS:NPOS)) )GOTO 420 C!skip digits after decimal pt NPOS = NPOS+1 GOTO 410 420 CONTINUE ENDIF C!allow for exponent IF( C_INPUT(NPOS:NPOS) .EQ. 'E' .OR. X C_INPUT(NPOS:NPOS) .EQ. 'D' )THEN IPOS = NPOS+1 IF( C_INPUT(IPOS:IPOS) .EQ. '+' .OR. X C_INPUT(IPOS:IPOS) .EQ. '-' )IPOS = IPOS+1 IF( L_NUMERIC(C_INPUT(IPOS:IPOS)) )THEN C!only if a digit follows the E can it be legal NPOS = IPOS 510 CONTINUE IF( .NOT. L_NUMERIC(C_INPUT(NPOS:NPOS)) )GOTO 520 NPOS = NPOS+1 GOTO 510 520 CONTINUE ENDIF ENDIF NPOS = NPOS-1 C!number runs from position 1 to NPOS NUSED = NPOS IF( NPOS .LE. 9 )THEN WRITE(UNIT=C_VAL,FMT=5501) NPOS 5501 FORMAT('(F',I1,'.0)') ELSE WRITE(UNIT=C_VAL,FMT=5502) NPOS 5502 FORMAT('(F',I2,'.0)') ENDIF READ(UNIT=C_INPUT(1:NPOS),FMT=C_VAL,IOSTAT=IO_CODE) VALUE C CCC WRITE(*,5509) C_INPUT(1:NPOS) , C_VAL , VALUE CCC5509 FORMAT( CCC X ' scanned text ',A/ CCC X ' using format ',A/ CCC X ' giving VALUE ',1PG14.7) C IF( IO_CODE .EQ. 0 )THEN NTYPE = NT_NUMBER ELSE NTYPE = NT_ERROR ENDIF C....................................................................... C If not a number, an error! C ELSE NTYPE = NT_ERROR NUSED = 1 ENDIF C....................................................................... 8000 CONTINUE RETURN END C C C C(((.................................................................... FUNCTION LAST_NONBLANK( CLINE ) C C Return the position of the last nonblank character in the input C character string. CLINE is CHARACTER*(*). Even if CLINE is all C blanks, LAST_NONBLANK will be returned as 1 so that operations of the C form CLINE(1:LAST_NONBLANK) won't be garbage. C))).................................................................... CHARACTER*(*) CLINE INTEGER LAST_NONBLANK C INTEGER NPOS C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Start at the end and work backwards until a nonblank is found. C Loop back to 100 to check position # NPOS each time. C NPOS = LEN( CLINE ) 100 CONTINUE C quit if at the beginning IF( NPOS .LE. 1 )GOTO 200 C quit if not a blank or a null IF( CLINE(NPOS:NPOS) .NE. ' ' .AND. . CLINE(NPOS:NPOS) .NE. CHAR(0) )GOTO 200 C move back one position and try again NPOS = NPOS - 1 GOTO 100 C....................................................................... 200 CONTINUE LAST_NONBLANK = NPOS RETURN END C C C FUNCTION HASSYM( SYM , NUM_CODE , C_CODE ) IMPLICIT NONE C INTEGER HASSYM CHARACTER*8 SYM INTEGER NUM_CODE CHARACTER*8 C_CODE(*) C INTEGER NCODE CHARACTER*1 SSS C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C HASSYM = 0 IF( NUM_CODE .LE. 0 ) RETURN SSS = SYM(1:1) C DO 1000 NCODE=1,NUM_CODE IF( C_CODE(NCODE) .EQ. 'PUSHSYM' )THEN IF( C_CODE(NCODE+1)(1:1) .EQ. SSS )THEN HASSYM = 1 RETURN ENDIF ENDIF 1000 CONTINUE C RETURN END C C C FUNCTION PAREVAL( NUM_CODE , C_CODE , R8VAL ) IMPLICIT NONE C REAL*8 PAREVAL INTEGER NUM_CODE CHARACTER*8 C_CODE(*) REAL*8 R8VAL(*) C INTEGER NUM_ESTACK PARAMETER ( NUM_ESTACK = 128 ) REAL*8 R8_EVAL(NUM_ESTACK) C INTEGER NEVAL , NCODE , IALPHA , ITM,NTM CHARACTER*8 C8_VAL , CNCODE REAL*8 R8_VAL , X , Y EQUIVALENCE ( C8_VAL , R8_VAL ) C C Internal library functions C REAL*8 QG , QGINV , BELL2 , RECT , STEP , BOOL , X LAND,LOR,LMOFN,MEDIAN , ZTONE , HMODE,LMODE, X GRAN,URAN,IRAN,ERAN,LRAN , ORSTAT , TENT, MAD REAL*8 ARGMAX,ARGNUM C C External library functions C REAL*8 DAI,DBI , DGAMMA , X DBESI0,DBESI1 , DBESJ0,DBESJ1 , DBESK0,DBESK1 , X DBESY0,DBESY1 , X DERF,DERFC C C Statistics functions (01 Mar 1999 - see parser_int.c) C REAL * 8 FICOTP , FICOPT , FICOTZ , X FITTTP , FITTPT , FITTTZ , X FIFTTP , FIFTPT , FIFTTZ , X FIZTTP , FIZTPT , FIZTTZ , X FICTTP , FICTPT , FICTTZ , X FIBTTP , FIBTPT , FIBTTZ , X FIBNTP , FIBNPT , FIBNTZ , X FIGTTP , FIGTPT , FIGTTZ , X FIPTTP , FIPTPT , FIPTTZ C REAL*8 R2D , D2R PARAMETER ( R2D = 57.29577951308232D+00 , X D2R = 0.01745329251994D+00 ) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C IF( NUM_CODE .LE. 0 )THEN PAREVAL = 0.D+0 GOTO 8000 ENDIF C----------------------------------------------------------------------- IALPHA = ICHAR('A') - 1 NEVAL = 0 NCODE = 0 C 1000 CONTINUE NCODE = NCODE + 1 CNCODE = C_CODE(NCODE) C....................................................................... IF( CNCODE .EQ. 'PUSHSYM' )THEN NEVAL = NEVAL + 1 R8_EVAL(NEVAL) = R8VAL( ICHAR(C_CODE(NCODE+1)(1:1))-IALPHA ) NCODE = NCODE + 1 C....................................................................... ELSEIF( CNCODE .EQ. 'PUSHNUM' )THEN NEVAL = NEVAL + 1 C8_VAL = C_CODE(NCODE+1) R8_EVAL(NEVAL) = R8_VAL NCODE = NCODE + 1 C....................................................................... ELSEIF( CNCODE .EQ. '+' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = R8_EVAL(NEVAL) + R8_EVAL(NEVAL+1) C....................................................................... ELSEIF( CNCODE .EQ. '-' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = R8_EVAL(NEVAL) - R8_EVAL(NEVAL+1) C....................................................................... ELSEIF( CNCODE .EQ. '*' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = R8_EVAL(NEVAL) * R8_EVAL(NEVAL+1) C....................................................................... ELSEIF( CNCODE .EQ. '/' )THEN NEVAL = NEVAL - 1 IF( R8_EVAL(NEVAL+1) .NE. 0.0D+0 )THEN R8_EVAL(NEVAL) = R8_EVAL(NEVAL) / R8_EVAL(NEVAL+1) ELSE R8_EVAL(NEVAL) = 0.0D+0 ENDIF C....................................................................... ELSEIF( CNCODE .EQ. '**' )THEN NEVAL = NEVAL - 1 IF( R8_EVAL(NEVAL) .GT. 0.0D+0 .OR. X (R8_EVAL(NEVAL) .NE. 0.D+00 .AND. X R8_EVAL(NEVAL+1) .EQ. DINT(R8_EVAL(NEVAL+1))) ) X X R8_EVAL(NEVAL) = R8_EVAL(NEVAL) ** R8_EVAL(NEVAL+1) C....................................................................... ELSEIF( CNCODE .EQ. '--' )THEN R8_EVAL(NEVAL) = -R8_EVAL(NEVAL) C....................................................................... ELSEIF( CNCODE .EQ. 'SIN' )THEN R8_EVAL(NEVAL) = SIN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'SIND' )THEN R8_EVAL(NEVAL) = SIN( D2R*R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'COS' )THEN R8_EVAL(NEVAL) = COS( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'COSD' )THEN R8_EVAL(NEVAL) = COS( D2R*R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'TAN' )THEN R8_EVAL(NEVAL) = TAN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'TAND' )THEN R8_EVAL(NEVAL) = TAN( D2R*R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'SQRT' )THEN R8_EVAL(NEVAL) = SQRT(ABS(R8_EVAL(NEVAL))) C....................................................................... ELSEIF( CNCODE .EQ. 'ABS' )THEN R8_EVAL(NEVAL) = ABS( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'EXP' )THEN R8_EVAL(NEVAL) = EXP( MIN(87.5D+0,R8_EVAL(NEVAL)) ) C....................................................................... ELSEIF( CNCODE .EQ. 'LOG' )THEN IF( R8_EVAL(NEVAL) .NE. 0.D+0 ) X R8_EVAL(NEVAL) = LOG(ABS(R8_EVAL(NEVAL))) C....................................................................... ELSEIF( CNCODE .EQ. 'LOG10' )THEN IF( R8_EVAL(NEVAL) .NE. 0.D+0 ) X R8_EVAL(NEVAL) = LOG10(ABS(R8_EVAL(NEVAL))) C....................................................................... ELSEIF( CNCODE .EQ. 'INT' )THEN R8_EVAL(NEVAL) = DINT( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'MAX' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = MAX( R8_EVAL(NEVAL) , R8_EVAL(NEVAL+1) ) C....................................................................... ELSEIF( CNCODE .EQ. 'MIN' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = MIN( R8_EVAL(NEVAL) , R8_EVAL(NEVAL+1) ) C....................................................................... ELSEIF( CNCODE .EQ. 'ASIN' )THEN IF( ABS(R8_EVAL(NEVAL)) .LT. 1.D+0 ) X R8_EVAL(NEVAL) = ASIN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'ACOS' )THEN IF( ABS(R8_EVAL(NEVAL)) .LT. 1.D+0 ) X R8_EVAL(NEVAL) = ACOS( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'ATAN' )THEN R8_EVAL(NEVAL) = ATAN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'ATAN2')THEN NEVAL = NEVAL - 1 IF( R8_EVAL(NEVAL).NE.0.D+0.AND.R8_EVAL(NEVAL+1).NE.0.D+0) X R8_EVAL(NEVAL) = ATAN2( R8_EVAL(NEVAL),R8_EVAL(NEVAL+1) ) C....................................................................... ELSEIF( CNCODE .EQ. 'GRAN')THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = GRAN( R8_EVAL(NEVAL),R8_EVAL(NEVAL+1) ) C....................................................................... ELSEIF( CNCODE .EQ. 'URAN' )THEN R8_EVAL(NEVAL) = URAN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'IRAN' )THEN R8_EVAL(NEVAL) = IRAN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'ERAN' )THEN R8_EVAL(NEVAL) = ERAN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'LRAN' )THEN R8_EVAL(NEVAL) = LRAN( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'SINH' )THEN IF( ABS(R8_EVAL(NEVAL)) .LT. 87.5 ) X R8_EVAL(NEVAL) = SINH( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'COSH' )THEN IF( ABS(R8_EVAL(NEVAL)) .LT. 87.5 ) X R8_EVAL(NEVAL) = COSH( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'TANH' )THEN R8_EVAL(NEVAL) = TANH( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'ASINH' )THEN X = ABS(R8_EVAL(NEVAL)) IF( X .LE. 10.D0 )THEN Y = X + SQRT(X**2+1.D0) ELSE Y = X * ( 1.D0 + SQRT(1.D0+(1.D0/X)**2) ) ENDIF Y = LOG(Y) IF( R8_EVAL(NEVAL) .LT. 0.D0 )THEN R8_EVAL(NEVAL) = -Y ELSE R8_EVAL(NEVAL) = Y ENDIF C....................................................................... ELSEIF( CNCODE .EQ. 'ACOSH' )THEN X = R8_EVAL(NEVAL) IF( X .GE. 1.D+0 )THEN IF( X .LE. 10.D0 )THEN Y = X + SQRT(X**2-1.D0) ELSE Y = X * ( 1.D0 + SQRT(1.D0-(1.D0/X)**2) ) ENDIF R8_EVAL(NEVAL) = LOG(Y) ENDIF C....................................................................... ELSEIF( CNCODE .EQ. 'ATANH' )THEN X = R8_EVAL(NEVAL) IF( ABS(X) .LT. 1.D0 ) X R8_EVAL(NEVAL) = 0.5D0 * LOG( (1.D0+X)/(1.D0-X) ) C....................................................................... ELSEIF( CNCODE .EQ. 'AI' )THEN R8_EVAL(NEVAL) = DAI( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'BI' )THEN R8_EVAL(NEVAL) = DBI( R8_EVAL(NEVAL) , 1 ) C....................................................................... ELSEIF( CNCODE .EQ. 'ERF' )THEN R8_EVAL(NEVAL) = DERF( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'ERFC' )THEN R8_EVAL(NEVAL) = DERFC( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'GAMMA' )THEN R8_EVAL(NEVAL) = DGAMMA( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'I0' )THEN R8_EVAL(NEVAL) = DBESI0( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'I1' )THEN R8_EVAL(NEVAL) = DBESI1( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'J0' )THEN R8_EVAL(NEVAL) = DBESJ0( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'J1' )THEN R8_EVAL(NEVAL) = DBESJ1( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'K0' )THEN R8_EVAL(NEVAL) = DBESK0( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'K1' )THEN R8_EVAL(NEVAL) = DBESK1( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'Y0' )THEN R8_EVAL(NEVAL) = DBESY0( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'Y1' )THEN R8_EVAL(NEVAL) = DBESY1( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'QG' )THEN R8_EVAL(NEVAL) = QG( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'QGINV' )THEN R8_EVAL(NEVAL) = QGINV( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'BELL2' )THEN R8_EVAL(NEVAL) = BELL2( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'RECT' )THEN R8_EVAL(NEVAL) = RECT( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'STEP' )THEN R8_EVAL(NEVAL) = STEP( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'TENT' )THEN R8_EVAL(NEVAL) = TENT( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'BOOL' )THEN R8_EVAL(NEVAL) = BOOL( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'ZTONE' )THEN R8_EVAL(NEVAL) = ZTONE( R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'NOTZERO' )THEN R8_EVAL(NEVAL) = BOOL( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'ISZERO' )THEN R8_EVAL(NEVAL) = 1.D+0 - BOOL( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'EQUALS' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = 1.D+0-BOOL(R8_EVAL(NEVAL)-R8_EVAL(NEVAL+1)) ELSEIF( CNCODE .EQ. 'ISPOSITI' )THEN R8_EVAL(NEVAL) = STEP( R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'ISNEGATI' )THEN R8_EVAL(NEVAL) = STEP( -R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'AND' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = LAND( NTM , R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'MEDIAN' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = MEDIAN( NTM , R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'MAD' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = MAD( NTM , R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'ORSTAT' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM NTM = NTM - 1 ITM = R8_EVAL(NEVAL) R8_EVAL(NEVAL) = ORSTAT( ITM,NTM , R8_EVAL(NEVAL+1) ) ELSEIF( CNCODE .EQ. 'HMODE' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = HMODE( NTM , R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'LMODE' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = LMODE( NTM , R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'OR' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = LOR( NTM , R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'MOFN' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM NTM = NTM - 1 ITM = R8_EVAL(NEVAL) R8_EVAL(NEVAL) = LMOFN( ITM,NTM , R8_EVAL(NEVAL+1) ) ELSEIF( CNCODE .EQ. 'ASTEP' )THEN NEVAL = NEVAL - 1 IF( ABS(R8_EVAL(NEVAL)) .GT. R8_EVAL(NEVAL+1) )THEN R8_EVAL(NEVAL) = 1.D+0 ELSE R8_EVAL(NEVAL) = 0.D+0 ENDIF ELSEIF( CNCODE .EQ. 'ARGMAX' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = ARGMAX( NTM , R8_EVAL(NEVAL) ) ELSEIF( CNCODE .EQ. 'ARGNUM' )THEN NTM = R8_EVAL(NEVAL) NEVAL = NEVAL - NTM R8_EVAL(NEVAL) = ARGNUM( NTM , R8_EVAL(NEVAL) ) C....................................................................... ELSEIF( CNCODE .EQ. 'FICO_T2P' )THEN NEVAL = NEVAL - 3 R8_EVAL(NEVAL) = FICOTP(ABS(R8_EVAL(NEVAL)), X R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2),R8_EVAL(NEVAL+3) ) ELSEIF( CNCODE .EQ. 'FICO_P2T' )THEN NEVAL = NEVAL - 3 R8_EVAL(NEVAL) = FICOPT(R8_EVAL(NEVAL) ,R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2),R8_EVAL(NEVAL+3) ) ELSEIF( CNCODE .EQ. 'FICO_T2Z' )THEN NEVAL = NEVAL - 3 R8_EVAL(NEVAL) = FICOTZ(R8_EVAL(NEVAL) ,R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2),R8_EVAL(NEVAL+3) ) C....................................................................... ELSEIF( CNCODE .EQ. 'FITT_T2P' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FITTTP(ABS(R8_EVAL(NEVAL)), X R8_EVAL(NEVAL+1)) ELSEIF( CNCODE .EQ. 'FITT_P2T' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FITTPT(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) ELSEIF( CNCODE .EQ. 'FITT_T2Z' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FITTTZ(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) C....................................................................... ELSEIF( CNCODE .EQ. 'FIFT_T2P' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIFTTP(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIFT_P2T' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIFTPT(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIFT_T2Z' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIFTTZ(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) C....................................................................... ELSEIF( CNCODE .EQ. 'FIZT_T2P' )THEN R8_EVAL(NEVAL) = FIZTTP(ABS(R8_EVAL(NEVAL))) ELSEIF( CNCODE .EQ. 'FIZT_P2T' )THEN R8_EVAL(NEVAL) = FIZTPT(R8_EVAL(NEVAL)) ELSEIF( CNCODE .EQ. 'FIZT_T2Z' )THEN R8_EVAL(NEVAL) = FIZTTZ(R8_EVAL(NEVAL)) C....................................................................... ELSEIF( CNCODE .EQ. 'FICT_T2P' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FICTTP(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) ELSEIF( CNCODE .EQ. 'FICT_P2T' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FICTPT(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) ELSEIF( CNCODE .EQ. 'FICT_T2Z' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FICTTZ(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) C....................................................................... ELSEIF( CNCODE .EQ. 'FIBT_T2P' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIBTTP(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIBT_P2T' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIBTPT(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIBT_T2Z' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIBTTZ(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) C....................................................................... ELSEIF( CNCODE .EQ. 'FIBN_T2P' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIBNTP(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIBN_P2T' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIBNPT(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIBN_T2Z' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIBNTZ(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) C....................................................................... ELSEIF( CNCODE .EQ. 'FIGT_T2P' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIGTTP(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIGT_P2T' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIGTPT(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) ELSEIF( CNCODE .EQ. 'FIGT_T2Z' )THEN NEVAL = NEVAL - 2 R8_EVAL(NEVAL) = FIGTTZ(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1), X R8_EVAL(NEVAL+2) ) C....................................................................... ELSEIF( CNCODE .EQ. 'FIPT_T2P' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FIPTTP(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) ELSEIF( CNCODE .EQ. 'FIPT_P2T' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FIPTPT(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) ELSEIF( CNCODE .EQ. 'FIPT_T2Z' )THEN NEVAL = NEVAL - 1 R8_EVAL(NEVAL) = FIPTTZ(R8_EVAL(NEVAL),R8_EVAL(NEVAL+1)) C....................................................................... ENDIF C....................................................................... IF( NCODE .LT. NUM_CODE )GOTO 1000 PAREVAL = R8_EVAL(NEVAL) C----------------------------------------------------------------------- 8000 CONTINUE RETURN END C C C SUBROUTINE PAREVEC( NUM_CODE , C_CODE , VA, VB, VC, VD, VE, X VF, VG, VH, VI, VJ, VK, VL, VM, VN, VO, X VP, VQ, VR, VS, VT, VU, VV, VW, VX, VY, VZ, X LVEC, VOUT ) IMPLICIT NONE C C Vector version of PAREVAL, where VA..VZ with length LVEC C are supplied as vectors. C [Modified by Raoqiong Tong, August 1997] C INTEGER NUM_CODE , LVEC CHARACTER*8 C_CODE(*) INTEGER NUM_ESTACK , NVMAX PARAMETER ( NUM_ESTACK = 101 , NVMAX = 64 ) REAL*8 VA(LVEC), VB(LVEC), VC(LVEC), VD(LVEC), VE(LVEC), X VF(LVEC), VG(LVEC), VH(LVEC), VI(LVEC), VJ(LVEC), X VK(LVEC), VL(LVEC), VM(LVEC), VN(LVEC), VO(LVEC), X VP(LVEC), VQ(LVEC), VR(LVEC), VS(LVEC), VT(LVEC), X VU(LVEC), VV(LVEC), VW(LVEC), VX(LVEC), VY(LVEC), X VZ(LVEC), VOUT(LVEC) C REAL*8 R8_EVAL(NVMAX,NUM_ESTACK) , R8VAL(NVMAX,26) C INTEGER NEVAL , NCODE , IALPHA , IV,IBV,IVBOT,IVTOP , X JF, NTM,ITM,JTM CHARACTER*8 C8_VAL , CNCODE , C2CODE REAL*8 R8_VAL , X , Y EQUIVALENCE ( C8_VAL , R8_VAL ) C C 14 Jul 1998: add 1D array for stack copy C REAL*8 SCOP(NUM_ESTACK) C C Internal library functions C REAL*8 QG , QGINV , BELL2 , RECT , STEP , BOOL , LAND, X LOR, LMOFN , MEDIAN , ZTONE , HMODE , LMODE , X GRAN,URAN,IRAN,ERAN,LRAN , ORSTAT , TENT, MAD REAL*8 ARGMAX,ARGNUM C C External library functions C REAL*8 DAI,DBI , DGAMMA , X DBESI0,DBESI1 , DBESJ0,DBESJ1 , DBESK0,DBESK1 , X DBESY0,DBESY1 , X DERF,DERFC C C Statistics functions (01 Mar 1999 - see parser_int.c) C REAL * 8 FICOTP , FICOPT , FICOTZ , X FITTTP , FITTPT , FITTTZ , X FIFTTP , FIFTPT , FIFTTZ , X FIZTTP , FIZTPT , FIZTTZ , X FICTTP , FICTPT , FICTTZ , X FIBTTP , FIBTPT , FIBTTZ , X FIBNTP , FIBNPT , FIBNTZ , X FIGTTP , FIGTPT , FIGTTZ , X FIPTTP , FIPTPT , FIPTTZ C REAL*8 R2D , D2R PARAMETER ( R2D = 57.29577951308232D+00 , X D2R = 0.01745329251994D+00 ) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C IF( NUM_CODE.LE.0 .OR. LVEC.LE.0 )GOTO 8000 C IALPHA = ICHAR('A') - 1 C----------------------------------------------------------------------- DO 5000 IBV=0,LVEC-1,NVMAX IVBOT = IBV + 1 IVTOP = IBV + NVMAX IF( IVTOP .GT. LVEC ) IVTOP = LVEC C ccc WRITE(*,9802) IVBOT,IVTOP ccc9802 FORMAT(' .. PAREVEC: loop from',I5,' to',I5) C DO 100 IV=IVBOT,IVTOP R8VAL(IV-IBV, 1) = VA(IV) 100 CONTINUE DO 101 IV=IVBOT,IVTOP R8VAL(IV-IBV, 2) = VB(IV) 101 CONTINUE DO 102 IV=IVBOT,IVTOP R8VAL(IV-IBV, 3) = VC(IV) 102 CONTINUE DO 103 IV=IVBOT,IVTOP R8VAL(IV-IBV, 4) = VD(IV) 103 CONTINUE DO 104 IV=IVBOT,IVTOP R8VAL(IV-IBV, 5) = VE(IV) 104 CONTINUE DO 105 IV=IVBOT,IVTOP R8VAL(IV-IBV, 6) = VF(IV) 105 CONTINUE DO 106 IV=IVBOT,IVTOP R8VAL(IV-IBV, 7) = VG(IV) 106 CONTINUE DO 107 IV=IVBOT,IVTOP R8VAL(IV-IBV, 8) = VH(IV) 107 CONTINUE DO 108 IV=IVBOT,IVTOP R8VAL(IV-IBV, 9) = VI(IV) 108 CONTINUE DO 109 IV=IVBOT,IVTOP R8VAL(IV-IBV,10) = VJ(IV) 109 CONTINUE DO 110 IV=IVBOT,IVTOP R8VAL(IV-IBV,11) = VK(IV) 110 CONTINUE DO 111 IV=IVBOT,IVTOP R8VAL(IV-IBV,12) = VL(IV) 111 CONTINUE DO 112 IV=IVBOT,IVTOP R8VAL(IV-IBV,13) = VM(IV) 112 CONTINUE DO 113 IV=IVBOT,IVTOP R8VAL(IV-IBV,14) = VN(IV) 113 CONTINUE DO 114 IV=IVBOT,IVTOP R8VAL(IV-IBV,15) = VO(IV) 114 CONTINUE DO 115 IV=IVBOT,IVTOP R8VAL(IV-IBV,16) = VP(IV) 115 CONTINUE DO 116 IV=IVBOT,IVTOP R8VAL(IV-IBV,17) = VQ(IV) 116 CONTINUE DO 117 IV=IVBOT,IVTOP R8VAL(IV-IBV,18) = VR(IV) 117 CONTINUE DO 118 IV=IVBOT,IVTOP R8VAL(IV-IBV,19) = VS(IV) 118 CONTINUE DO 119 IV=IVBOT,IVTOP R8VAL(IV-IBV,20) = VT(IV) 119 CONTINUE DO 120 IV=IVBOT,IVTOP R8VAL(IV-IBV,21) = VU(IV) 120 CONTINUE DO 121 IV=IVBOT,IVTOP R8VAL(IV-IBV,22) = VV(IV) 121 CONTINUE DO 122 IV=IVBOT,IVTOP R8VAL(IV-IBV,23) = VW(IV) 122 CONTINUE DO 123 IV=IVBOT,IVTOP R8VAL(IV-IBV,24) = VX(IV) 123 CONTINUE DO 124 IV=IVBOT,IVTOP R8VAL(IV-IBV,25) = VY(IV) 124 CONTINUE DO 125 IV=IVBOT,IVTOP R8VAL(IV-IBV,26) = VZ(IV) 125 CONTINUE C NEVAL = 0 NCODE = 0 C 1000 CONTINUE NCODE = NCODE + 1 CNCODE = C_CODE(NCODE) ccc WRITE(*,9803) CNCODE ccc9803 FORMAT(' .. PAREVEC: opcode=',A) C....................................................................... IF( CNCODE .EQ. 'PUSHSYM' )THEN JF = ICHAR(C_CODE(NCODE+1)(1:1))-IALPHA IF( NCODE+2 .LE. NUM_CODE )THEN C2CODE = C_CODE(NCODE+2) ELSE C2CODE = 'q' ENDIF IF( C2CODE .EQ. '+' )THEN NCODE = NCODE + 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X + R8VAL(IV-IBV,JF) ENDDO ELSEIF( C2CODE .EQ. '-' )THEN NCODE = NCODE + 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X - R8VAL(IV-IBV,JF) ENDDO ELSEIF( C2CODE .EQ. '*' )THEN NCODE = NCODE + 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X * R8VAL(IV-IBV,JF) ENDDO ELSEIF( C2CODE .EQ. '/' )THEN NCODE = NCODE + 2 DO IV=IVBOT,IVTOP IF( R8VAL(IV-IBV,JF) .NE. 0.D+00 )THEN R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X / R8VAL(IV-IBV,JF) ELSE R8_EVAL(IV-IBV,NEVAL) = 0.D+00 ENDIF ENDDO ELSE NEVAL = NEVAL + 1 NCODE = NCODE + 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8VAL(IV-IBV,JF) ENDDO ENDIF C....................................................................... ELSEIF( CNCODE .EQ. 'PUSHNUM' )THEN C8_VAL = C_CODE(NCODE+1) IF( NCODE+2 .LE. NUM_CODE )THEN C2CODE = C_CODE(NCODE+2) ELSE C2CODE = 'q' ENDIF IF( C2CODE .EQ. '+' )THEN NCODE = NCODE + 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X + R8_VAL ENDDO ELSEIF( C2CODE .EQ. '-' )THEN NCODE = NCODE + 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X - R8_VAL ENDDO ELSEIF( C2CODE .EQ. '*' )THEN NCODE = NCODE + 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X * R8_VAL ENDDO ELSEIF( C2CODE .EQ. '/' )THEN NCODE = NCODE + 2 IF( R8_VAL .NE. 0.D+00 )THEN R8_VAL = 1.D+00 / R8_VAL DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X * R8_VAL ENDDO ELSE DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = 0.D+00 ENDDO ENDIF ELSE NCODE = NCODE + 1 NEVAL = NEVAL + 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_VAL ENDDO ENDIF C....................................................................... ELSEIF( CNCODE .EQ. '+' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X + R8_EVAL(IV-IBV,NEVAL+1) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. '-' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X - R8_EVAL(IV-IBV,NEVAL+1) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. '*' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X * R8_EVAL(IV-IBV,NEVAL+1) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. '/' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP IF( R8_EVAL(IV-IBV,NEVAL+1) .NE. 0.0D+0 )THEN R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X / R8_EVAL(IV-IBV,NEVAL+1) ELSE R8_EVAL(IV-IBV,NEVAL) = 0.D+00 ENDIF ENDDO C....................................................................... ELSEIF( CNCODE .EQ. '**' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP IF( R8_EVAL(IV-IBV,NEVAL) .GT. 0.0D+0 .OR. X (R8_EVAL(IV-IBV,NEVAL) .NE. 0.D+00 .AND. X R8_EVAL(IV-IBV,NEVAL+1) .EQ. X DINT(R8_EVAL(IV-IBV,NEVAL+1)) ) ) X X R8_EVAL(IV-IBV,NEVAL) = R8_EVAL(IV-IBV,NEVAL) X ** R8_EVAL(IV-IBV,NEVAL+1) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. '--' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = -R8_EVAL(IV-IBV,NEVAL) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'SIN' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = SIN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'SIND' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = SIN( D2R*R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'COS' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = COS( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'COSD' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = COS( D2R*R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'TAN' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = TAN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'TAND' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = TAN( D2R*R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'SQRT' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = SQRT(ABS(R8_EVAL(IV-IBV,NEVAL))) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ABS' )THEN DO IV=IVBOT,IVTOP ccc WRITE(*,9809) IV ccc9809 FORMAT(' about to ABS #',I5) R8_EVAL(IV-IBV,NEVAL) = ABS( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'EXP' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = X EXP(MIN(87.5,R8_EVAL(IV-IBV,NEVAL))) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'LOG' )THEN DO IV=IVBOT,IVTOP IF( R8_EVAL(IV-IBV,NEVAL) .NE. 0.D+0 ) X R8_EVAL(IV-IBV,NEVAL) = LOG(ABS(R8_EVAL(IV-IBV,NEVAL))) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'LOG10' )THEN DO IV=IVBOT,IVTOP IF( R8_EVAL(IV-IBV,NEVAL) .NE. 0.D+0 ) X R8_EVAL(IV-IBV,NEVAL) = LOG10(ABS(R8_EVAL(IV-IBV,NEVAL))) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'INT' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DINT( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'MAX' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = MAX( R8_EVAL(IV-IBV,NEVAL) , X R8_EVAL(IV-IBV,NEVAL+1) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'MIN' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = MIN( R8_EVAL(IV-IBV,NEVAL) , X R8_EVAL(IV-IBV,NEVAL+1) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ASIN' )THEN DO IV=IVBOT,IVTOP IF( ABS(R8_EVAL(IV-IBV,NEVAL)) .LT. 1.D+0 ) X R8_EVAL(IV-IBV,NEVAL) = ASIN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ACOS' )THEN DO IV=IVBOT,IVTOP IF( ABS(R8_EVAL(IV-IBV,NEVAL)) .LT. 1.D+0 ) X R8_EVAL(IV-IBV,NEVAL) = ACOS( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ATAN' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = ATAN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ATAN2')THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP IF( R8_EVAL(IV-IBV,NEVAL) .NE.0.D+0 .AND. X R8_EVAL(IV-IBV,NEVAL+1).NE.0.D+0 ) X R8_EVAL(IV-IBV,NEVAL) = ATAN2( R8_EVAL(IV-IBV,NEVAL) , X R8_EVAL(IV-IBV,NEVAL+1) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'GRAN')THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = GRAN( R8_EVAL(IV-IBV,NEVAL) , X R8_EVAL(IV-IBV,NEVAL+1) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'URAN')THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = URAN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'IRAN')THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = IRAN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ERAN')THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = ERAN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'LRAN')THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = LRAN( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'SINH' )THEN DO IV=IVBOT,IVTOP IF( ABS(R8_EVAL(IV-IBV,NEVAL)) .LT. 87.5 ) X R8_EVAL(IV-IBV,NEVAL) = SINH( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'COSH' )THEN DO IV=IVBOT,IVTOP IF( ABS(R8_EVAL(IV-IBV,NEVAL)) .LT. 87.5 ) X R8_EVAL(IV-IBV,NEVAL) = COSH( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'TANH' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = TANH( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ASINH' )THEN DO IV=IVBOT,IVTOP X = ABS(R8_EVAL(IV-IBV,NEVAL)) IF( X .LE. 10.D0 )THEN Y = X + SQRT(X**2+1.D0) ELSE Y = X * ( 1.D0 + SQRT(1.D0+(1.D0/X)**2) ) ENDIF Y = LOG(Y) IF( R8_EVAL(IV-IBV,NEVAL) .LT. 0.D0 )THEN R8_EVAL(IV-IBV,NEVAL) = -Y ELSE R8_EVAL(IV-IBV,NEVAL) = Y ENDIF ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ACOSH' )THEN DO IV=IVBOT,IVTOP X = R8_EVAL(IV-IBV,NEVAL) IF( X .GE. 1.D0 )THEN IF( X .LE. 10.D0 )THEN Y = X + SQRT(X**2-1.D0) ELSE Y = X * ( 1.D0 + SQRT(1.D0-(1.D0/X)**2) ) ENDIF R8_EVAL(IV-IBV,NEVAL) = LOG(Y) ENDIF ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ATANH' )THEN DO IV=IVBOT,IVTOP X = R8_EVAL(IV-IBV,NEVAL) IF( ABS(X) .LT. 1.D0 ) X R8_EVAL(IV-IBV,NEVAL) = 0.5D0 * LOG( (1.D0+X)/(1.D0-X) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'AI' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DAI( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'BI' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBI( R8_EVAL(IV-IBV,NEVAL) , 1 ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'ERF' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DERF( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'ERFC' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DERFC( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'GAMMA' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DGAMMA( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'I0' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESI0( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'I1' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESI1( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'J0' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESJ0( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'J1' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESJ1( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'K0' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESK0( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'K1' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESK1( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'Y0' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESY0( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'Y1' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = DBESY1( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'QG' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = QG( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'QGINV' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = QGINV( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'BELL2' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = BELL2( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'RECT' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = RECT( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'STEP' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = STEP( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'TENT' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = TENT( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'BOOL' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = BOOL( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'ZTONE' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = ZTONE( R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'NOTZERO' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = BOOL( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'ISZERO' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = 1.D+0-BOOL(R8_EVAL(IV-IBV,NEVAL)) ENDDO ELSEIF( CNCODE .EQ. 'EQUALS' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = 1.D+0 X -BOOL(R8_EVAL(IV-IBV,NEVAL)-R8_EVAL(IV-IBV,NEVAL+1)) ENDDO ELSEIF( CNCODE .EQ. 'ISPOSITI' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = STEP( R8_EVAL(IV-IBV,NEVAL) ) ENDDO ELSEIF( CNCODE .EQ. 'ISNEGATI' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = STEP( -R8_EVAL(IV-IBV,NEVAL) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'AND' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = LAND( NTM, SCOP ) ENDDO ELSEIF( CNCODE .EQ. 'MEDIAN' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = MEDIAN( NTM, SCOP ) ENDDO ELSEIF( CNCODE .EQ. 'MAD' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = MAD( NTM, SCOP ) ENDDO ELSEIF( CNCODE .EQ. 'ORSTAT' )THEN NTM = R8_EVAL(1,NEVAL) NEVAL = NEVAL - NTM NTM = NTM - 1 DO IV=IVBOT,IVTOP ITM = R8_EVAL(IV-IBV,NEVAL) DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM) ENDDO R8_EVAL(IV-IBV,NEVAL) = ORSTAT(ITM,NTM,SCOP) ENDDO ELSEIF( CNCODE .EQ. 'HMODE' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = HMODE( NTM, SCOP ) ENDDO ELSEIF( CNCODE .EQ. 'LMODE' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = LMODE( NTM, SCOP ) ENDDO ELSEIF( CNCODE .EQ. 'OR' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = LOR( NTM, SCOP ) ENDDO ELSEIF( CNCODE .EQ. 'MOFN' )THEN NTM = R8_EVAL(1,NEVAL) NEVAL = NEVAL - NTM NTM = NTM - 1 DO IV=IVBOT,IVTOP ITM = R8_EVAL(IV-IBV,NEVAL) DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM) ENDDO R8_EVAL(IV-IBV,NEVAL) = LMOFN(ITM,NTM,SCOP) ENDDO ELSEIF( CNCODE .EQ. 'ASTEP' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP IF( ABS(R8_EVAL(IV-IBV,NEVAL)) .GT. R8_EVAL(IV-IBV,NEVAL+1) )THEN R8_EVAL(IV-IBV,NEVAL) = 1.D+0 ELSE R8_EVAL(IV-IBV,NEVAL) = 0.D+0 ENDIF ENDDO ELSEIF( CNCODE .EQ. 'ARGMAX' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = ARGMAX( NTM, SCOP ) ENDDO ELSEIF( CNCODE .EQ. 'ARGNUM' )THEN NTM = R8_EVAL(1, NEVAL) NEVAL = NEVAL - NTM DO IV=IVBOT,IVTOP DO JTM=1,NTM SCOP(JTM) = R8_EVAL(IV-IBV,NEVAL+JTM-1) ENDDO R8_EVAL(IV-IBV,NEVAL) = ARGNUM( NTM, SCOP ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FICO_T2P' )THEN NEVAL = NEVAL - 3 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FICOTP(ABS(R8_EVAL(IV-IBV,NEVAL)), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2), X R8_EVAL(IV-IBV,NEVAL+3) ) ENDDO ELSEIF( CNCODE .EQ. 'FICO_P2T' )THEN NEVAL = NEVAL - 3 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FICOPT(R8_EVAL(IV-IBV,NEVAL) , X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2), X R8_EVAL(IV-IBV,NEVAL+3) ) ENDDO ELSEIF( CNCODE .EQ. 'FICO_T2Z' )THEN NEVAL = NEVAL - 3 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FICOTZ(R8_EVAL(IV-IBV,NEVAL) , X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2), X R8_EVAL(IV-IBV,NEVAL+3) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FITT_T2P' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FITTTP(ABS(R8_EVAL(IV-IBV,NEVAL)), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO ELSEIF( CNCODE .EQ. 'FITT_P2T' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FITTPT(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO ELSEIF( CNCODE .EQ. 'FITT_T2Z' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FITTTZ(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FIFT_T2P' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIFTTP(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIFT_P2T' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIFTPT(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIFT_T2Z' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIFTTZ(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FIZT_T2P' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIZTTP(ABS(R8_EVAL(IV-IBV,NEVAL))) ENDDO ELSEIF( CNCODE .EQ. 'FIZT_P2T' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIZTPT(R8_EVAL(IV-IBV,NEVAL)) ENDDO ELSEIF( CNCODE .EQ. 'FIZT_T2Z' )THEN DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIZTTZ(R8_EVAL(IV-IBV,NEVAL)) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FICT_T2P' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FICTTP(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO ELSEIF( CNCODE .EQ. 'FICT_P2T' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FICTPT(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO ELSEIF( CNCODE .EQ. 'FICT_T2Z' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FICTTZ(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FIBT_T2P' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIBTTP(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIBT_P2T' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIBTPT(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIBT_T2Z' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIBTTZ(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FIBN_T2P' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIBNTP(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIBN_P2T' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIBNPT(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIBN_T2Z' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIBNTZ(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FIGT_T2P' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIGTTP(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIGT_P2T' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIGTPT(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO ELSEIF( CNCODE .EQ. 'FIGT_T2Z' )THEN NEVAL = NEVAL - 2 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIGTTZ(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1), X R8_EVAL(IV-IBV,NEVAL+2) ) ENDDO C....................................................................... ELSEIF( CNCODE .EQ. 'FIPT_T2P' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIPTTP(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO ELSEIF( CNCODE .EQ. 'FIPT_P2T' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIPTPT(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO ELSEIF( CNCODE .EQ. 'FIPT_T2Z' )THEN NEVAL = NEVAL - 1 DO IV=IVBOT,IVTOP R8_EVAL(IV-IBV,NEVAL) = FIPTTZ(R8_EVAL(IV-IBV,NEVAL), X R8_EVAL(IV-IBV,NEVAL+1)) ENDDO C....................................................................... ENDIF C---------------------------------------------------------------------- IF( NCODE .LT. NUM_CODE )GOTO 1000 C DO 4990 IV=IVBOT,IVTOP VOUT(IV) = R8_EVAL(IV-IBV,NEVAL) 4990 CONTINUE C 5000 CONTINUE C----------------------------------------------------------------------- 8000 CONTINUE RETURN END C C C FUNCTION ZTONE( X ) IMPLICIT REAL*8 (A-H,O-Z) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C IF( X .LE. 0.D+00 )THEN ZTONE = 0.D+00 ELSEIF( X .GE. 1.0 )THEN ZTONE = 1.D+00 ELSE Y = (0.5 * 3.1415926535897932D+00) * (1.6D+00 * X - 0.8D+00) ZTONE = 0.50212657D+00 * ( TANH(TAN(Y)) + 0.99576486D+00 ) ENDIF RETURN END C C C FUNCTION QG( X ) C C Compute the reversed cdf of a Gaussian. C IMPLICIT REAL*8 (A-H,O-Z) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C QG = 0.5D+0 * DERFC( X / 1.414213562373095D+0 ) RETURN END C C C CCC The UNIF() function is now in parser_int.c, CCC where it calls CCC upon the C library to do the dirty work. C CCC FUNCTION UNIF( XJUNK ) CCC IMPLICIT REAL*8 (A-H,O-Z) CCC PARAMETER ( IA = 99992 , IB = 12345 , IT = 99991 ) CCC PARAMETER ( F = 1.00009D-05 ) CCC DATA IX / 271 / CCCC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ CCC IX = MOD( IA*IX+IB , IT ) CCC UNIF = F * IX CCC RETURN CCC END C C C CCC FUNCTION UNIF( XJUNK ) CCC IMPLICIT REAL*8 (A-H,O-Z) CCCC CCCC FACTOR - INTEGER OF THE FORM 8*K+5 AS CLOSE AS POSSIBLE CCCC TO 2**26 * (SQRT(5)-1)/2 (GOLDEN SECTION) CCCC TWO28 = 2**28 (I.E. 28 SIGNIFICANT BITS FOR DEVIATES) CCCC CCC PARAMETER ( FACTOR = 41475557.0D+00 , TWO28 = 268435456.0D+00 ) CCCC CCC DATA R / 0.D+00 / CCCC CCCC RETURNS SAMPLE U FROM THE 0,1 -UNIFORM DISTRIBUTION CCCC BY A MULTIPLICATIVE CONGRUENTIAL GENERATOR OF THE FORM CCCC R := R * FACTOR (MOD 1) . CCCC IN THE FIRST CALL R IS INITIALIZED TO CCCC R := IR / 2**28 , CCCC WHERE IR MUST BE OF THE FORM IR = 4*K+1. CCCC THEN R ASSUMES ALL VALUES 0 < (4*K+1)/2**28 < 1 DURING CCCC A FULL PERIOD 2**26 OF SUNIF. CCCC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ CCCC CCC IF( R .EQ. 0.D+00 ) R = 4000001.D+00 / TWO28 CCC R = DMOD(R*FACTOR,1.0D+00) CCC UNIF = R CCC RETURN CCC END C C C FUNCTION IRAN( TOP ) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 IRAN C C Return an integer uniformly distributed among 0..TOP C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ IRAN = DINT( (TOP+1.D+00)*UNIF(0.D+00) ) RETURN END C C C FUNCTION ERAN( TOP ) IMPLICIT REAL*8 (A-H,O-Z) C C Return an exponentially distributed deviate: F(x) = 1-exp(-x/top) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 100 U1 = UNIF(0.D+00) IF( U1 .LE. 0.D+00 ) GOTO 100 ERAN = -TOP*LOG(U1) RETURN END C C C FUNCTION LRAN( TOP ) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 LRAN C C Return a logistically distributed deviate: F(x) = 1/[1+exp(-x/top)] C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 100 U1 = UNIF(0.D+00) IF( U1 .LE. 0.D+00 .OR. U1 .GE. 1.D+00 ) GOTO 100 LRAN = TOP*LOG( 1.D+00/U1 - 1.D+00 ) RETURN END C C C FUNCTION URAN( X ) IMPLICIT REAL*8 (A-H,O-Z) C C Return a U(0,X) random variable. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C URAN = X * UNIF(0.D+00) RETURN END C C C FUNCTION GRAN2( B , S ) IMPLICIT REAL*8 (A-H,O-Z) C C Compute a Gaussian random deviate with mean B and stdev S C INTEGER IP DATA IP / 0 / C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C IF( IP .EQ. 0 )THEN 100 U1 = UNIF(0.D+00) IF( U1 .LE. 0.D+00 ) GOTO 100 U2 = UNIF(0.D+00) GRAN2 = B + S * SQRT(-2.0D+00*LOG(U1)) * SIN(6.2831853D+00*U2) IP = 1 ELSE GRAN2 = B + S * SQRT(-2.0D+00*LOG(U1)) * COS(6.2831853D+00*U2) IP = 0 ENDIF RETURN END C C C FUNCTION GRAN1( B , S ) IMPLICIT REAL*8 (A-H,O-Z) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C G = -6.D+00 + UNIF(1.D+00) + UNIF(2.D+00) + UNIF(3.D+00) X + UNIF(4.D+00) + UNIF(5.D+00) + UNIF(6.D+00) X + UNIF(7.D+00) + UNIF(8.D+00) + UNIF(9.D+00) X + UNIF(10.D+0) + UNIF(11.D+0) + UNIF(12.D+0) GRAN1 = B + S*G RETURN END C C C FUNCTION GRAN( B , S ) IMPLICIT REAL*8 (A-H,O-Z) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C UU = UNIF(0.0D+00) IF( UU .LE. 0.5D+00 )THEN GRAN = GRAN1(B,S) ELSE GRAN = GRAN2(B,S) ENDIF RETURN END C C C FUNCTION QGINV( P ) C C Return x such that Q(x)=P, for 0 < P < 1. Q=reversed Gaussian cdf. C IMPLICIT REAL*8 (A-H,O-Z) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C DP = P IF( DP .GT. 0.5D+0 )DP = 1.D+0 - DP IF( DP .LE. 0.D+0 )THEN DX = 13.D+00 GOTO 8000 ENDIF C C Step 1: use 26.2.23 from Abramowitz and Stegun C DT = SQRT( -2.D+0 * LOG(DP) ) DX = DT X - ((.010328D+0*DT + .802853D+0)*DT + 2.525517D+0) X /(((.001308D+0*DT + .189269D+0)*DT + 1.432788D+0)*DT + 1.D+0) C C Step 2: do 3 Newton steps to improve this C DO 100 NEWT=1,3 DQ = 0.5D+0 * DERFC( DX / 1.414213562373095D+0 ) - DP DDQ = EXP( -0.5D+00 * DX * DX ) / 2.506628274631000D+0 DX = DX + DQ / DDQ 100 CONTINUE C 8000 CONTINUE IF( P .GT. 0.5D+00 )THEN QGINV = -DX ELSE QGINV = DX ENDIF C RETURN END C C C FUNCTION BELL2(X) REAL*8 BELL2 , X , AX C... AX = ABS(X) IF( AX .LE. 0.5D+00 )THEN BELL2 = 1.0D+00 - 1.3333333333333333D+00 * AX*AX ELSEIF( AX .LE. 1.5D+00 )THEN BELL2 = 0.666666666666667D+00 * (1.5D+00 - AX)**2 ELSE BELL2 = 0.0D+00 ENDIF RETURN END C C C FUNCTION RECT(X) REAL*8 RECT , X , AX AX = ABS(X) IF( AX .LE. 0.5D+00 )THEN RECT = 1.0D+0 ELSE RECT = 0.0D+0 ENDIF RETURN END C C C FUNCTION STEP(X) REAL*8 STEP , X IF( X .LE. 0.D+0 )THEN STEP = 0.D+0 ELSE STEP = 1.D+0 ENDIF RETURN END C C C FUNCTION TENT(X) REAL*8 TENT , X , AX AX = ABS(X) IF( AX .GE. 1.D+0 )THEN TENT = 0.D+0 ELSE TENT = 1.D+0 - AX ENDIF RETURN END C C C FUNCTION BOOL(X) REAL*8 BOOL , X IF( X .EQ. 0.D+0 )THEN BOOL = 0.D+0 ELSE BOOL = 1.D+0 ENDIF RETURN END C C C FUNCTION LAND(N,X) REAL*8 LAND,X(N) INTEGER N , I LAND = 0.D+0 DO 100 I=1,N IF( X(I) .EQ. 0.D+0 ) RETURN 100 CONTINUE LAND = 1.D+0 RETURN END C C C SUBROUTINE BSORT(N,X) REAL *8 X(N) , TMP INTEGER N , I , IT C------------------------------------ Bubble sort 50 IT = 0 DO 100 I=2,N IF( X(I-1) .GT. X(I) )THEN TMP = X(I) X(I) = X(I-1) X(I-1) = TMP IT = 1 ENDIF 100 CONTINUE IF( IT .NE. 0 )GOTO 50 RETURN END C C C FUNCTION ORSTAT(M,N,X) REAL*8 ORSTAT,X(N) INTEGER M,N , I C IF( N .LE. 1 )THEN ORSTAT = X(1) RETURN ENDIF C I = M IF( I .LE. 0 )THEN I = 1 ELSEIF( I .GT. N )THEN I = N ENDIF CALL BSORT(N,X) ORSTAT = X(I) RETURN END C C C FUNCTION MEDIAN(N,X) REAL *8 MEDIAN , X(N) , TMP INTEGER N , IT C IF( N .EQ. 1 )THEN MEDIAN = X(1) RETURN ELSEIF( N .EQ. 2 )THEN MEDIAN = 0.5D+00 * (X(1)+X(2)) RETURN ELSEIF( N .EQ. 3 )THEN IF( X(1) .GT. X(2) )THEN TMP = X(2) X(2) = X(1) X(1) = TMP ENDIF IF( X(1) .GT. X(3) )THEN MEDIAN = X(1) ELSEIF( X(2) .GT. X(3) )THEN MEDIAN = X(3) ELSE MEDIAN = X(2) ENDIF RETURN ENDIF C C--- sort it C CALL BSORT(N,X) C C--- Even N --> average of middle 2 C--- Odd N --> middle 1 C IT = N/2 IF( 2*IT .EQ. N )THEN MEDIAN = 0.5D+00 * (X(IT)+X(IT+1)) ELSE MEDIAN = X(IT+1) ENDIF RETURN END C C C FUNCTION MAD(N,X) REAL *8 MAD , X(N) , TMP , MEDIAN INTEGER N , IT C IF( N .EQ. 1 )THEN MAD = 0.D+00 RETURN ELSEIF( N .EQ. 2 )THEN MAD = 0.5D+00*ABS(X(1)-X(2)) RETURN ENDIF C TMP = MEDIAN(N,X) DO 100 IT=1,N X(IT) = ABS(X(IT)-TMP) 100 CONTINUE MAD = MEDIAN(N,X) RETURN END C C C FUNCTION ARGMAX(N,X) REAL*8 ARGMAX , X(N) , TMP INTEGER N, I , IT , NZ C TMP = X(1) IT = 1 NZ = 0 IF( TMP .EQ. 0.D+00 ) NZ = 1 DO 100 I=2,N IF( X(I) .GT. TMP )THEN IT = I TMP = X(I) ENDIF IF( X(I) .EQ. 0.D+00 ) NZ = NZ+1 100 CONTINUE IF( NZ .EQ. N )THEN ARGMAX = 0.D+00 ELSE ARGMAX = IT ENDIF RETURN END C C C FUNCTION ARGNUM(N,X) REAL*8 ARGNUM , X(N) INTEGER N, I, NZ C NZ = 0 DO 100 I=1,N IF( X(I) .NE. 0.D+00 ) NZ = NZ+1 100 CONTINUE ARGNUM = NZ RETURN END C C C FUNCTION HMODE(N,X) REAL *8 HMODE , VAL , VB , X(N) INTEGER N , I , IV , IB C IF( N .EQ. 1 )THEN HMODE = X(1) RETURN ENDIF C CALL BSORT(N,X) C VAL = X(1) IV = 1 IB = 0 DO 100 I=2,N IF( X(I) .NE. VAL )THEN IF( IV .GE. IB )THEN VB = VAL IB = IV ENDIF VAL = X(I) IV = 1 ELSE IV = IV+1 ENDIF 100 CONTINUE IF( IV .GE. IB ) VB = VAL HMODE = VB RETURN END C C C FUNCTION LMODE(N,X) REAL *8 LMODE , VAL , VB , X(N) INTEGER N , I , IV , IB C IF( N .EQ. 1 )THEN LMODE = X(1) RETURN ENDIF C CALL BSORT(N,X) C VAL = X(1) IV = 1 IB = 0 DO 100 I=2,N IF( X(I) .NE. VAL )THEN IF( IV .GT. IB )THEN VB = VAL IB = IV ENDIF VAL = X(I) IV = 1 ELSE IV = IV+1 ENDIF 100 CONTINUE IF( IV .GT. IB ) VB = VAL LMODE = VB RETURN END C C C FUNCTION LOR(N,X) REAL*8 LOR,X(N) INTEGER N , I LOR = 1.D+0 DO 100 I=1,N IF( X(I) .NE. 0.D+0 ) RETURN 100 CONTINUE LOR = 0.D+0 RETURN END C C C FUNCTION LMOFN(M,N,X) REAL*8 LMOFN,X(N) INTEGER M,N , I,C C = 0 DO 100 I=1,N IF( X(I) .NE. 0.D+0 ) C = C + 1 100 CONTINUE IF( C .GE. M )THEN LMOFN = 1.D+0 ELSE LMOFN = 0.D+0 ENDIF RETURN END C C C REAL*8 FUNCTION DAI( X ) REAL*8 X CALL QQQERR DAI = 0.D+0 RETURN END REAL*8 FUNCTION DBI( X , I ) REAL*8 X CALL QQQERR DBI = 0.D+0 RETURN END REAL*8 FUNCTION DGAMMA( X ) REAL*8 X CALL QQQERR DGAMMA = 0.D+0 RETURN END REAL*8 FUNCTION DBESI0( X ) REAL*8 X CALL QQQERR DBESI0 = 0.D+0 RETURN END REAL*8 FUNCTION DBESI1( X ) REAL*8 X CALL QQQERR DBESI1 = 0.D+0 RETURN END ccc REAL*8 FUNCTION DBESJ0( X ) ccc REAL*8 X ccc CALL QQQERR ccc END ccc REAL*8 FUNCTION DBESJ1( X ) ccc REAL*8 X ccc CALL QQQERR ccc END REAL*8 FUNCTION DBESK0( X ) REAL*8 X CALL QQQERR DBESK0 = 0.D+0 RETURN END REAL*8 FUNCTION DBESK1( X ) REAL*8 X CALL QQQERR DBESK1 = 0.D+0 RETURN END ccc REAL*8 FUNCTION DBESY0( X ) ccc REAL*8 X ccc CALL QQQERR ccc END ccc REAL*8 FUNCTION DBESY1( X ) ccc REAL*8 X ccc CALL QQQERR ccc END ccc REAL*8 FUNCTION DERF( X ) ccc REAL*8 X ccc CALL QQQERR ccc END ccc REAL*8 FUNCTION DERFC( X ) ccc REAL*8 X ccc CALL QQQERR ccc END C SUBROUTINE QQQERR WRITE(*,999) 999 FORMAT('*** PARSER: unimplemented function ***') RETURN END