/*REXX*/ /* ------------------------------------------------------------------ */ /* This source code contains REXX implementations for the following */ /* mathematical functions: */ /* */ /* sqrt(x) sin(x) asin(x) fctrl(x) */ /* cbrt(x) cos(x) acos(x) gamma(x) */ /* ln(x) tan(x) atan(x) lngamma(x) */ /* exp(x) cot(x) acot(x) erf(x) */ /* sinh(x) sec(x) asec(x) erfc(x) */ /* cosh(x) csc(x) acsc(x) power(x,y) */ /* */ /* ------------------------------------------------------------------ */ /* Copyright 1988, B.E. Chi (BEC@ALBNYVM1.BITNET), Computing Services */ /* University of NY at Albany. power(x,y) added by MikeA. */ /* */ /* EVX has been given away free on BITNET for years. It is freeware, */ /* and is not for sale. You are free to alter this program, as long */ /* as you note that in the comments. */ /* Ported to OS/2 by */ /* Michael Antonio (713221.1742@CompuServ.COM), 1990 */ /* */ /* EVX expression evaluates any numeric (or, more generally, */ /* any REXX) expression. Examples: */ /* */ /* EVX 24*(35/2)-.5 returns 419.5 */ /* EVX max(3.5,3/2,17/3) returns 5.6666667 */ /* EVX 3*sqrt(2.7) returns 4.9295030 */ /* EVX reverse(abcxyz) returns ZYXCBA */ /* EVX date() returns 19 Feb 1987 (for instance) */ /* */ /* The expression may contain certain mathematical functions if */ /* required. Among the standard REXX functions, the ones using */ /* numeric arguments and returning a numeric value are: */ /* */ /* ABS(x), MAX(x,y,...,z), MIN(x,y,...,z), SIGN(x). */ /* */ /* Two utility functions are rad(x) and deg(x), which convert */ /* degrees to radians and radians to degrees, respectively. */ /* */ /* Finally, the following algebraic and transcendental functions */ /* are provided: */ /* */ /* sqrt(x) sin(x) asin(x) fctrl(x) */ /* cbrt(x) cos(x) acos(x) gamma(x) */ /* ln(x) tan(x) atan(x) lngamma(x) */ /* exp(x) cot(x) acot(x) erf(x) */ /* sinh(x) sec(x) asec(x) erfc(x) */ /* cosh(x) csc(x) acsc(x) power(x,y) */ /* */ /* The symbols "pi" and "e" are defined to have their customary */ /* values. Also, the value of the expression is saved with */ /* symbol "x" which may be used in an expression in a subsequent */ /* call of EVX. Example: */ /* */ /* EVX 3*4 returns 12. */ /* EVX x/6 returns 2. */ /* */ /* Calling EVX with no argument invokes nonstop mode. For each */ /* subsquent expression typed, the corresponding value is */ /* displayed. To return to the operating system, enter a blank */ /* line. The symbol "x" may be used as before. */ /* */ /* EVX may be called from within another exec. To retrieve the */ /* result, execute the statement */ /* */ /* "GLOBALV STACK evx_p; PULL x" (in CMS) */ /* or */ /* "x = VALUE('evx_p', 'OS2ENVIRONMENT')" (under OS/2). */ /* Alternately, */ /* EVX expression "(STACK" (in CMS) */ /* or */ /* EVX expression "-STACK" (in OS/2) */ /* suppresses printing and PUSHes the result onto the */ /* current stack in the usual manner. */ /*--------------------------------------------------------------------*/ /* install an error handler */ SIGNAL ON syntax /* get the parameter */ PARSE ARG evx_expr /* get the current operating system */ PARSE SOURCE env . IF env = 'OS/2' THEN delim = '-' ELSE delim = '(' IF evx_expr = "?" THEN DO /* show the usage help */ DO i = 22 IF SOURCELINE(i) = "" THEN EXIT SAY SOURCELINE(i) END /* DO i = 22 */ END /* IF */ x = LASTPOS( delim,evx_expr ) evx_stack = x > 0 & TRANSLATE(SUBSTR(evx_expr,x+1)) = "STACK" IF evx_stack THEN evx_expr = SUBSTR(evx_expr,1,x-1) /* evx_stack = 0 or 1; evx_expr = actual expression. */ evx_stop = evx_expr \= "" IF evx_stop THEN PUSH evx_expr ELSE IF evx_stack THEN DO SAY "EVX nonstop and stack modes may not be used together." EXIT 1 END ELSE SAY "Enter EVX nonstop mode:" NUMERIC DIGITS 12 /* define some constants */ pi = 3.14159265359 e = 2.71828182846 /* get the result of the previous calculation */ x = getvar( "EVX_P" ) IF x = "" THEN x = 0 /* default is zero */ DO UNTIL evx_stop PARSE PULL evx_expr IF evx_expr = "" THEN evx_stop = 1 ELSE DO /* increase the accuracy */ NUMERIC DIGITS 12 /* calculate the expression */ INTERPRET "x = "evx_expr NUMERIC DIGITS 8 /* round the result */ IF DATATYPE(x,N) & \DATATYPE(x,W) THEN x = x + 0 IF evx_stack THEN PUSH x ELSE SAY evx_expr" =" x END END /* save the result in the environment variable */ rt=setvar("EVX_P", x) /* reset the NUMERIC DIGITS settings */ NUMERIC DIGITS EXIT 0 /* ================================================================== */ /* function: Convert degrees to radians */ /* */ /* call: frad(x) */ /* */ rad: PROCEDURE ARG x CALL argtest "rad",x RETURN .0174532925199*x /* ================================================================== */ /* function: Convert radians to degrees */ /* */ /* call: deg(x) */ /* */ deg: PROCEDURE ARG x CALL argtest "deg",x RETURN 57.2957795131*x /* ================================================================== */ /* sqrt(x) */ /* */ /* Method: */ /* Assuming x >= 0, set x = a * 4**n, choosing n so .25 < a <= 1. */ /* Then sqrt(x) = 2**n sqrt(a). Evaluate sqrt(a) by using 0.707 as */ /* an initial approximation and then iterating four times, using */ /* Newton's method. */ /* */ sqrt: PROCEDURE ARG x CALL argtest "sqrt",x IF x < 0 THEN DO SAY "sqrt argument" x " out of range." EXIT END IF x = 0 THEN RETURN 0 n = 0 DO WHILE x > 1 x = x/4 n = n + 1 END DO WHILE x <= .25 x = 4*x n = n - 1 END y = .707 DO 4 y = .5*(y + x/y) END RETURN 2**n*y /* ================================================================== */ /* cbrt(x) */ /* */ /* Method: */ /* Set x = a * 8**n, choosing n so .125 < a <= 1. Then */ /* cbrt(x) = 2**n cbrt(a). Evaluate cbrt(a) by using 0.737 as */ /* an initial approximation and then iterating five times, using */ /* Newton's method. */ /* */ cbrt: PROCEDURE ARG x CALL argtest "cbrt",x IF x = 0 THEN RETURN 0 s = SIGN(x) x = ABS(x) n = 0 DO WHILE x > 1 x = x/8 n = n + 1 END DO WHILE x <= .125 x = 8*x n = n - 1 END y = .737 DO 5 y = (2*y + x/(y**2))/3 END RETURN s*2**n*y /* ================================================================== */ /* sin(x) */ /* */ /* Method: */ /* x is first reduced so that ABS(x) <= pi/2. Then a polynomial */ /* approximation is used for evaluation. */ /* */ sin: PROCEDURE ARG x CALL argtest "sin",x s=SIGN(x) x=ABS(x) x = x//6.28318530718 IF x > 3.14159265359 THEN DO x = x - 3.14159265359 s = -s END IF x > 1.57079632679 THEN x = 3.14159265359 - x y = x*x x = s*x*(((((-.0000000239*y + .0000027526)*y - .0001984090)*y + , .0083333315)*y - .1666666664)*y + 1) RETURN x /* ================================================================== */ /* cos(x) */ /* */ /* Method: */ /* x is first reduced so that ABS(x) <= pi/2. Then a */ /* polynomial approximation is used for evaluation. */ /* */ cos: PROCEDURE ARG x CALL argtest "cos",x s=1 x=ABS(x) x = x//6.28318530718 IF x > 3.14159265359 THEN DO x = x - 3.14159265359 s = -s END IF x > 1.57079632679 THEN DO x = 3.14159265359 - x s = -s END y = x*x x = s*(((((-.0000002605*y + .0000247609)*y - .0013888397)*y + , .0416666418)*y - .4999999963)*y + 1) RETURN x /* ================================================================== */ /* tan(x) */ /* */ /* Method: */ /* x is first reduced so that ABS(x) <= pi/2. Then if */ /* ABS(x) > pi/4, tan(x) = cot(pi/2 - x). Otherwise, a polynomial */ /* approximation is used for evaluation. */ /* */ tan: PROCEDURE ARG x CALL argtest "tan",x s=SIGN(x) x=ABS(x) x = x//6.28318530718 IF x > 3.14159265359 THEN DO x = x - 3.14159265359 s = -s END IF x > 1.57079632679 THEN DO x = 3.14159265359 - x s = -s END IF x > .785398163397 THEN x = cot(1.57079632679 - x) ELSE DO y = x*x x = x*((((((.0095168091*y + .0029005250)*y + .0245650893)*y + , .0533740603)*y + .1333923995)*y + .3333314036)*y + 1) END RETURN s*x /* ================================================================== */ /* cot(x) */ /* */ /* Method: */ /* x is first reduced so that ABS(x) <= pi/2. Then if */ /* ABS(x) > pi/4, cot(x) = tan(pi/2 - x). Otherwise, a polynomial */ /* approximation is used for evaluation. */ /* */ cot: PROCEDURE ARG x CALL argtest "cot",x s=SIGN(x) x=ABS(x) x = x//6.28318530718 IF x > 3.14159265359 THEN DO x = x - 3.14159265359 s = -s END IF x > 1.57079632679 THEN DO x = 3.14159265359 - x s = -s END IF x > .785398163397 THEN x = tan(1.57079632679 - x) ELSE DO y = x*x x = (((((-.0000262619*y - .0002078504)*y - .0021177168)*y - , .0222220287)*y - .3333333410)*y + 1)/x END RETURN s*x /* ================================================================== */ /* sec(x) */ /* */ /* Method: */ /* sec(x) = 1/cos(x) */ /* */ sec: PROCEDURE ARG x CALL argtest "sec",x RETURN 1/cos(x) /* ================================================================== */ /* csc(x) */ /* */ /* Method: */ /* csc(x) = 1/sin(x) */ /* */ csc: PROCEDURE ARG x CALL argtest "csc",x RETURN 1/sin(x) /* ================================================================== */ /* exp(x) */ /* */ /* Method: */ /* Let n = integer part of x/ln 2, f = fractional part of x/ln 2. */ /* Then */ /* exp(x) = exp((n + f) ln 2) = exp(ln 2**n) exp (f ln 2) */ /* = 2**n exp (f ln 2), */ /* where f ln 2 <= ln 2. (See A&S p. 90, ex. 11.) */ /* The remaining exponential is evaluated using a polynomial */ /* approximation (see A&S 4.2.45). */ /* */ exp: PROCEDURE ARG x CALL argtest "exp",x s = SIGN(x) x = ABS(x/.69314718056) n = x % 1 d = (x - n)*.69314718056 x = ((((((-.0001413161*d + .0013298820)*d - .0083013598)*d + , .0416573475)*d - .1666653019)*d + .4999999206)*d - .9999999995)*d + 1 x = 2**n/x IF s = -1 THEN x = 1/x RETURN x /* ================================================================== */ /* ln(x) */ /* */ /* Method: */ /* Assuming x > 0, set x = a * 2**n, choosing n so that 1 <= a <= 2.*/ /* Then ln(x) = ln(a) + n*ln(2). Evaluate ln(a) using a polynomial */ /* approximation (A&S 4.1.44). */ /* */ ln: PROCEDURE ARG x CALL argtest "ln",x IF x <= 0 THEN DO SAY "ln argument" x " out of range." EXIT END n = 0 DO WHILE x > 2 x = x/2 n = n + 1 END DO WHILE x < 1 x = 2*x n = n - 1 END x = x - 1 x = (((((((-.0064535442*x + .0360884937)*x - .0953293897)*x + , .1676540711)*x - .2407338084)*x + .3317990258)*x - .4998741238)*x + , .9999964239)*x x = x + .6931471806*n RETURN x /* ================================================================== */ /* power(x,y) */ /* */ /* Method: */ /* x(y) = e^(y*ln(x)) */ /* */ power: PROCEDURE ARG x,y CALL argtest "power",x CALL argtest "power",y IF datatype(x,W) & datatype(y,W) THEN pow = x**y ELSE DO IF x<0 THEN DO SAY "power argument" y "out of range." EXIT END pow = exp(y*ln(x)) END RETURN pow /* ================================================================== */ /* asin(x) */ /* */ /* Method: */ /* A polynomial approximation is used. */ /* */ asin: PROCEDURE ARG x CALL argtest "asin",x IF x<-1 | x>1 THEN DO SAY "asin argument" x " out of range." EXIT END s = SIGN(x) x = ABS(x) x = (((((((-.0012624911*x + .0066700901)*x - .0170881256)*x + , .0308918810)*x - .0501743046)*x + .0889789874)*x - , .2145988016)*x + 1.570796305)*sqrt(1 - x) RETURN s*(1.57079632679 - x) /* ================================================================== */ /* acos(x) */ /* */ /* Method: */ /* Use acos(x) = pi/2 - asin(x) */ /* */ acos: PROCEDURE ARG x CALL argtest "acos",x IF x<-1 | x>1 THEN DO SAY "acos argument" x " out of range." EXIT END RETURN 1.57079632679 - asin(x) /* ================================================================== */ /* atan(x) */ /* */ /* Method: */ /* Use atan(x) = pi/4 - atan((1-x)/(1+x)). Evaluate */ /* the latter function using a polynomial approximation, taking */ /* advantage of the fact that its argument is less than one as */ /* long as x > -1. */ /* */ atan: PROCEDURE ARG x CALL argtest "atan",x IF x = 0 THEN RETURN 0 s=SIGN(x) x=ABS(x) x = (x - 1)/(x + 1) y = x*x x = ((((((((.0028662257*y - .0161657367)*y + .0429096138)*y - , .0752896400)*y + .1065626393)*y - .1420889944)*y + , .1999355085)*y - .3333314528)*y + 1)*x RETURN .785398163397 + s*x /* ================================================================== */ /* acot(x) */ /* */ /* Method: */ /* Use acot(x) = pi/2 - atan(x) */ /* */ acot: PROCEDURE ARG x CALL argtest "acot",x RETURN 1.57079632679 - atan(x) /* ================================================================== */ /* asec(x) */ /* */ /* Method: */ /* Use asec(x) = acos(1/x) */ /* */ asec: PROCEDURE ARG x CALL argtest "asec",x IF x>-1 & x<1 THEN DO SAY "asec argument" x " out of range." EXIT END RETURN acos(1/x) /* ================================================================== */ /* acsc(x) */ /* */ /* Method: */ /* Use acsc(x) = asin(1/x) */ /* */ acsc: PROCEDURE ARG x CALL argtest "acsc",x IF x>-1 & x<1 THEN DO SAY "acsc argument" x " out of range." EXIT END RETURN asin(1/x) /* ================================================================== */ /* sinh(x) */ /* */ sinh: PROCEDURE ARG x CALL argtest "sinh",x RETURN .5*(exp(x) - exp(-x)) /* ================================================================== */ /* cosh(x) */ /* */ cosh: PROCEDURE ARG x CALL argtest "cosh",x RETURN .5*(exp(x) + exp(-x)) /* ================================================================== */ /* tanh(x) */ /* */ tanh: PROCEDURE ARG x CALL argtest "tanh",x RETURN (exp(x) - exp(-x))/(exp(x) + exp(-x)) /* ================================================================== */ /* gamma(x) */ /* */ /* Method: */ /* Use */ /* gamma(x) = (x - 1)(x - 2)...(x - n)gamma(x - n + 1) if x >= 1, */ /* or */ /* gamma(x) = gamma(x + n + 1)/x(x + 1)(x + 2)...(x + n) if x < 1, */ /* in either case, choosing n so that the argument of the gamma */ /* function on the right hand side is in the range 0 through 1. */ /* Evaluate this gamma function using a polynomial approximation. */ /* */ gamma: PROCEDURE ARG x CALL argtest "gamma",x p = 1 SELECT WHEN x >= 1 THEN DO WHILE x >= 2 x = x - 1 p = x*p END WHEN x%1 = x THEN DO SAY "gamma argument" x "out of range." EXIT END OTHERWISE DO UNTIL x >= 1 IF x <= -1 & ABS(p) >= ABS(1E999999999/x) THEN RETURN 0 p = x*p x = x + 1 END p = 1/p END /* SELECT */ IF x \= 1 THEN DO x = x - 1 x = (((((((.035868343*x - .193527818)*x + .482199394)*x - , .756704078)*x + .918206857)*x - .897056937)*x + , .988205891)*x - .577191652)*x + 1 p = x*p END RETURN p /* ================================================================== */ /* lngamma(x), x > 0 */ /* */ /* Method: */ /* If x < 5, calculate ln(gamma(x)). */ /* Otherwise, use Stirling's formula through x**-7. */ /* */ lngamma: PROCEDURE ARG x CALL argtest "lngamma",x IF x < 0 THEN DO SAY "lngamma argument" x "out of range." EXIT END IF x < 5 THEN RETURN ln(gamma(x)) RETURN (x - .5)*ln(x) - x + .91893853320 + 1/(12*x) - 1/(360*x**3) + , 1/(1260*x**5) - 1/(1680*x**7) /* ================================================================== */ /* fctrl(x) */ /* */ /* Method: */ /* fctrl(x) = x! = gamma(x + 1): */ /* */ fctrl: PROCEDURE ARG x CALL argtest "fctrl",x IF x%1 = x & x < 0 THEN DO SAY "fctrl argument" x "out of range." EXIT END RETURN gamma(x + 1) /* ================================================================== */ /* erf(x) */ /* */ /* Method: */ /* One or another polynomial approximation is used, the choice */ /* depending upon the magnitude of x. (CACM Algorithm 209.) */ /* */ erf: PROCEDURE ARG x CALL argtest "erf",x s = SIGN(x) x = .707106781187*ABS(x) IF x >= 3 | x = 0 THEN RETURN s IF x < 1 THEN DO w = x*x x = ((((((((.000124818987*w - .001075204047)*w + .005198775019)*w - , .019198292004)*w + .059054035642)*w - .151968751364)*w + , .319152932694)*w - .531923007300)*w + .797884560593)*2*x END ELSE DO x = x - 2 x = (((((((((((((-.000045255659*x + .000152529290)*x - , .000019538132)*x - .000676904986)*x + .001390604284)*x - , .000794620820)*x - .002034254874)*x + .006549791214)*x - , .010557625006)*x + .011630447319)*x - .009279453341)*x + , .005353579108)*x - .002141268741)*x + .000535310849)*x + , .999936657524 END RETURN s*x /* ================================================================== */ /* erfc(x) */ /* */ /* Method: */ /* erfc(x) = 1 - erf(x) */ /* */ erfc: PROCEDURE ARG x CALL argtest "erfc",x RETURN 1 - erf(x) /* ------------------- misc. sub routines of EVX -------------------- */ /* ------------------------------------------------------------------ */ /* function: check if the argument for a function is numeric */ /* */ /* call: ArgTest functionName, functionArgument */ /* */ /* where: functionName - name of the function */ /* functionArgument - argument for the function */ /* */ /* returns: nothing */ /* */ /* note: exits the program if the argument is not numeric */ /* */ Argtest: PROCEDURE PARSE ARG name,x IF DATATYPE(x,"N") THEN RETURN SAY name "argument" x "not a number." EXIT /* ------------------------------------------------------------------ */ /* function: get the value of an environment variable */ /* */ /* call: GetVar varName */ /* */ /* where: varName - name of the environment variable */ /* */ /* returns: the value of the environment variable (or "") */ /* */ Getvar: PROCEDURE PARSE UPPER ARG varName PARSE UPPER SOURCE env . rt = "" IF env = 'CMS' THEN DO 'GLOBALV STACK' varName PULL rt END ELSE IF env = 'OS/2' THEN rt = value(varName,,'OS2ENVIRONMENT') RETURN rt /* ------------------------------------------------------------------ */ /* function: set the value for an environment variable */ /* */ /* call: SetVar varName, varValue */ /* */ /* where: varName - name of the environment variable */ /* varValue - value for the environment variable */ /* */ /* returns: the old value of the environment variable (or "") */ /* */ Setvar: PROCEDURE PARSE UPPER ARG varName,value PARSE UPPER SOURCE env . rt = "" IF env = 'CMS' THEN 'GLOBALV SET' varName value ELSE IF env = 'OS/2' THEN rt = value(varName,value,'OS2ENVIRONMENT') RETURN rt /* ------------------------------------------------------------------ */ /* syntax error handler */ /* */ SYNTAX: say '"'evx_expr'" is not quite valid syntax. Try again!' exit