home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 18 REXX / 18-REXX.zip / EVX.ZIP / EVX.CMD
OS/2 REXX Batch file  |  1990-09-17  |  17KB  |  546 lines
  1. /* Copyright 1988, B. E. Chi (BEC@ALBNYVM1.BITNET), Computing Services, */
  2. /* University of NY at Albany.  power(x,y) added by MikeA.              */
  3. /* EVX has been given away free on BITNET for years.  It is freeware,   */
  4. /* and is not for sale.  You are free to alter this program, as long    */
  5. /* as you note that in the comments.                                    */
  6. /* Ported to OS/2 by Michael Antonio (713221.1742@CompuServ.COM), 1990
  7. EVX  expression   evaluates any numeric (or, more generally,
  8. any REXX) expression.  Examples:
  9.  
  10. EVX 24*(35/2)-.5        returns 419.5
  11. EVX max(3.5,3/2,17/3)   returns 5.6666667
  12. EVX 3*sqrt(2.7)         returns 4.9295030
  13. EVX reverse(abcxyz)     returns ZYXCBA
  14. EVX date()              returns 19 Feb 1987  (for instance)
  15.  
  16. The expression may contain certain mathematical functions if
  17. required.  Among the standard REXX functions, the ones using
  18. numeric arguments and returning a numeric value are:
  19.  
  20.       ABS(x), MAX(x,y,...,z), MIN(x,y,...,z), SIGN(x).
  21.  
  22. Two utility functions are rad(x) and deg(x),  which  convert
  23. degrees to radians and radians to degrees, respectively.
  24.  
  25. Finally,  the  following  algebraic and transcendental func-
  26. tions are provided:
  27.  
  28.   sqrt(x)     sin(x)      asin(x)     fctrl(x)
  29.   cbrt(x)     cos(x)      acos(x)     gamma(x)
  30.   ln(x)       tan(x)      atan(x)     lngamma(x)
  31.   exp(x)      cot(x)      acot(x)     erf(x)
  32.   sinh(x)     sec(x)      asec(x)     erfc(x)
  33.   cosh(x)     csc(x)      acsc(x)     power(x,y)
  34.  
  35. The symbols "pi" and "e" are defined to have their customary
  36. values.  Also, the value of the  expression  is  saved  with
  37. symbol  "x"  which  may be used in an expression in a subse-
  38. quent call of EVX.  Example:
  39.  
  40. EVX 3*4  returns 12.
  41. EVX x/6  returns 2.
  42.  
  43. Calling EVX with no argument invokes nonstop mode.  For each
  44. subsquent expression typed, the corresponding value is  dis-
  45. played.    To return to the operating system,  enter a blank 
  46. line.  The symbol "x" may be used as before.
  47.  
  48. EVX may be called from within another exec.  To retrieve the
  49. result, execute the statement:
  50.    "GLOBALV  STACK  evx_p;  PULL x"   
  51.       in CMS and:
  52.     "x = VALUE('evx_p', 'OS2ENVIRONMENT')"
  53.       under OS/2.
  54. Alternately,  EVX  expression "(STACK" in CMS or "-STACK" in 
  55. OS/2  suppresses printing and  PUSHes  the result  onto  the 
  56. current stack in the usual manner.
  57. */
  58. SIGNAL ON syntax
  59. PARSE ARG evx_expr
  60. PARSE SOURCE env .
  61. IF env = 'OS/2' THEN
  62.    delim = '-'
  63. ELSE
  64.    delim = '('
  65.  
  66. IF evx_expr = "?" THEN DO i = 7
  67.   IF SOURCELINE(i) = "*/" THEN EXIT
  68.   SAY SOURCELINE(i)
  69. END
  70.  
  71. x = LASTPOS(delim,evx_expr)
  72. evx_stack = x > 0 & TRANSLATE(SUBSTR(evx_expr,x+1)) = "STACK"
  73. IF evx_stack THEN evx_expr = SUBSTR(evx_expr,1,x-1)
  74.                  /* evx_stack = 0 or 1; evx_expr = actual expression. */
  75.  
  76. evx_stop = evx_expr \= ""
  77. IF evx_stop THEN PUSH evx_expr
  78. ELSE IF evx_stack THEN DO
  79.   SAY "EVX nonstop and stack modes may not be used together."
  80.   EXIT 1
  81. END
  82. ELSE SAY "Enter EVX nonstop mode:"
  83.  
  84. NUMERIC DIGITS 12
  85. pi = 3.14159265359
  86. e = 2.71828182846
  87. x = getvar("EVX_P")
  88. IF x = "" THEN x = 0
  89.  
  90. DO UNTIL evx_stop
  91.   PARSE PULL evx_expr
  92.   IF evx_expr = "" THEN evx_stop = 1
  93.   ELSE DO
  94.     NUMERIC DIGITS 12
  95.     INTERPRET "x = "evx_expr
  96.     NUMERIC DIGITS 8
  97.     IF DATATYPE(x,N) & \DATATYPE(x,W) THEN x = x + 0
  98.     IF evx_stack THEN PUSH x
  99.     ELSE SAY evx_expr" =" x
  100.   END
  101. END
  102.  
  103. rt=setvar("EVX_P", x)
  104. NUMERIC DIGITS
  105. EXIT 0
  106.  
  107. /*====================================================================*/
  108. /* rad(x).  Convert degrees to radians.                               */
  109. rad: PROCEDURE
  110. ARG x
  111. CALL argtest "rad",x
  112. RETURN .0174532925199*x
  113. /*====================================================================*/
  114. /* deg(x).  Convert radians to degrees.                               */
  115. deg: PROCEDURE
  116. ARG x
  117. CALL argtest "deg",x
  118. RETURN 57.2957795131*x
  119. /*====================================================================*/
  120. /* sqrt(x).  Method:  assuming x >= 0, set x = a * 4**n, choosing n so
  121. .25 < a <= 1.  Then sqrt(x) = 2**n sqrt(a).  Evaluate sqrt(a) by using
  122. 0.707 as an initial approximation and then iterating four times, using
  123. Newton's method.                                                      */
  124. sqrt: PROCEDURE
  125. ARG x
  126. CALL argtest "sqrt",x
  127. IF x < 0 THEN DO
  128.   SAY "sqrt argument" x " out of range."
  129.   EXIT
  130. END
  131. IF x = 0 THEN RETURN 0
  132. n = 0
  133. DO WHILE x > 1
  134.   x = x/4
  135.   n = n + 1
  136. END
  137. DO WHILE x <= .25
  138.   x = 4*x
  139.   n = n - 1
  140. END
  141. y = .707
  142. DO 4
  143.   y = .5*(y + x/y)
  144. END
  145. RETURN 2**n*y
  146. /*====================================================================*/
  147. /* cbrt(x).  Method:  set x = a * 8**n, choosing n so .125 < a <= 1.
  148. Then cbrt(x) = 2**n cbrt(a).  Evaluate cbrt(a) by using 0.737 as an
  149. initial approximation and then iterating five times, using Newton's
  150. method.                                                               */
  151. cbrt: PROCEDURE
  152. ARG x
  153. CALL argtest "cbrt",x
  154. IF x = 0 THEN RETURN 0
  155. s = SIGN(x)
  156. x = ABS(x)
  157. n = 0
  158. DO WHILE x > 1
  159.   x = x/8
  160.   n = n + 1
  161. END
  162. DO WHILE x <= .125
  163.   x = 8*x
  164.   n = n - 1
  165. END
  166. y = .737
  167. DO 5
  168.   y = (2*y + x/(y**2))/3
  169. END
  170. RETURN s*2**n*y
  171. /*====================================================================*/
  172. /* sin(x).  Method:  x is first reduced so that ABS(x) <= pi/2.  Then a
  173. polynomial approximation is used for evaluation.                      */
  174. sin: PROCEDURE
  175. ARG x
  176. CALL argtest "sin",x
  177. s=SIGN(x)
  178. x=ABS(x)
  179. x = x//6.28318530718
  180. IF x > 3.14159265359 THEN DO
  181.   x = x - 3.14159265359
  182.   s = -s
  183. END
  184. IF x > 1.57079632679 THEN
  185.   x = 3.14159265359 - x
  186. y = x*x
  187. x = s*x*(((((-.0000000239*y + .0000027526)*y - .0001984090)*y + ,
  188.     .0083333315)*y - .1666666664)*y + 1)
  189. RETURN x
  190. /*====================================================================*/
  191. /* cos(x).  Method:  x is first reduced so that ABS(x) <= pi/2.  Then a
  192. polynomial approximation is used for evaluation.                      */
  193. cos: PROCEDURE
  194. ARG x
  195. CALL argtest "cos",x
  196. s=1
  197. x=ABS(x)
  198. x = x//6.28318530718
  199. IF x > 3.14159265359 THEN DO
  200.   x = x - 3.14159265359
  201.   s = -s
  202. END
  203. IF x > 1.57079632679 THEN DO
  204.   x = 3.14159265359 - x
  205.   s = -s
  206. END
  207. y = x*x
  208. x = s*(((((-.0000002605*y + .0000247609)*y - .0013888397)*y + ,
  209.     .0416666418)*y - .4999999963)*y + 1)
  210. RETURN x
  211. /*====================================================================*/
  212. /* tan(x).  Method:  x is first reduced so that ABS(x) <= pi/2.  Then if
  213. ABS(x) > pi/4, tan(x) = cot(pi/2 - x).  Otherwise, a polynomial approxi-
  214. mation is used for evaluation.                                        */
  215. tan: PROCEDURE
  216. ARG x
  217. CALL argtest "tan",x
  218. s=SIGN(x)
  219. x=ABS(x)
  220. x = x//6.28318530718
  221. IF x > 3.14159265359 THEN DO
  222.   x = x - 3.14159265359
  223.   s = -s
  224. END
  225. IF x > 1.57079632679 THEN DO
  226.   x = 3.14159265359 - x
  227.   s = -s
  228. END
  229. IF x > .785398163397 THEN
  230.   x = cot(1.57079632679 - x)
  231. ELSE DO
  232.   y = x*x
  233.   x = x*((((((.0095168091*y + .0029005250)*y + .0245650893)*y + ,
  234.     .0533740603)*y + .1333923995)*y + .3333314036)*y + 1)
  235. END
  236. RETURN s*x
  237. /*====================================================================*/
  238. /* cot(x).  Method:  x is first reduced so that ABS(x) <= pi/2.  Then if
  239. ABS(x) > pi/4, cot(x) = tan(pi/2 - x).  Otherwise, a polynomial approxi-
  240. mation is used for evaluation.                                        */
  241. cot: PROCEDURE
  242. ARG x
  243. CALL argtest "cot",x
  244. s=SIGN(x)
  245. x=ABS(x)
  246. x = x//6.28318530718
  247. IF x > 3.14159265359 THEN DO
  248.   x = x - 3.14159265359
  249.   s = -s
  250. END
  251. IF x > 1.57079632679 THEN DO
  252.   x = 3.14159265359 - x
  253.   s = -s
  254. END
  255. IF x > .785398163397 THEN
  256.   x = tan(1.57079632679 - x)
  257. ELSE DO
  258.   y = x*x
  259.   x = (((((-.0000262619*y - .0002078504)*y - .0021177168)*y - ,
  260.     .0222220287)*y - .3333333410)*y + 1)/x
  261. END
  262. RETURN s*x
  263. /*====================================================================*/
  264. /* sec(x) .                                                           */
  265. sec: PROCEDURE
  266. ARG x
  267. CALL argtest "sec",x
  268. RETURN 1/cos(x)
  269. /*====================================================================*/
  270. /* csc(x) .                                                           */
  271. csc: PROCEDURE
  272. ARG x
  273. CALL argtest "csc",x
  274. RETURN 1/sin(x)
  275. /*====================================================================*/
  276. /* exp(x).  Method:  Let n = integer part of x/ln 2, f = fractional part
  277. of x/ln 2.  Then exp(x) = exp((n + f) ln 2) = exp(ln 2**n) exp (f ln 2)
  278. = 2**n exp (f ln 2), where f ln 2 <= ln 2.  (See A&S p. 90, ex. 11.)
  279. The remaining exponential is evaluated using a polynomial approximation
  280. (see A&S 4.2.45).                                                     */
  281. exp: PROCEDURE
  282. ARG x
  283. CALL argtest "exp",x
  284. s = SIGN(x)
  285. x = ABS(x/.69314718056)
  286. n = x % 1
  287. d = (x - n)*.69314718056
  288. x = ((((((-.0001413161*d + .0013298820)*d - .0083013598)*d + ,
  289.    .0416573475)*d - .1666653019)*d + .4999999206)*d - .9999999995)*d + 1
  290. x = 2**n/x
  291. IF s = -1 THEN x = 1/x
  292. RETURN x
  293. /*====================================================================*/
  294. /* ln(x).  Method:  assuming x > 0, set x = a * 2**n, choosing n so that
  295. 1 <= a <= 2.  Then ln(x) = ln(a) + n*ln(2).  Evaluate ln(a) using a
  296. polynomial approximation (A&S 4.1.44).                                */
  297. ln: PROCEDURE
  298. ARG x
  299. CALL argtest "ln",x
  300. IF x <= 0 THEN DO
  301.   SAY "ln argument" x " out of range."
  302.   EXIT
  303. END
  304. n = 0
  305. DO WHILE x > 2
  306.   x = x/2
  307.   n = n + 1
  308. END
  309. DO WHILE x < 1
  310.   x = 2*x
  311.   n = n - 1
  312. END
  313. x = x - 1
  314. x = (((((((-.0064535442*x + .0360884937)*x - .0953293897)*x + ,
  315.    .1676540711)*x - .2407338084)*x + .3317990258)*x - .4998741238)*x + ,
  316.    .9999964239)*x
  317. x = x + .6931471806*n
  318. RETURN x
  319. /*====================================================================*/
  320. /* power(x,y).  Method:  x(y) = e^(y*ln(x)).                          */
  321. power: PROCEDURE
  322. ARG x,y
  323. CALL argtest "power",x
  324. CALL argtest "power",y
  325. IF datatype(x,W) & datatype(y,W) THEN
  326.    pow = x**y
  327. ELSE
  328.    DO
  329.    IF x<0 THEN 
  330.       DO
  331.       SAY "power argument" y "out of range."
  332.       EXIT
  333.       END
  334.    pow = exp(y*ln(x))
  335.    END
  336. RETURN pow
  337. /*====================================================================*/
  338. /* asin(x).  Method:  a polynomial approximation is used.             */
  339. asin: PROCEDURE
  340. ARG x
  341. CALL argtest "asin",x
  342. IF x<-1 | x>1 THEN DO
  343.   SAY "asin argument" x " out of range."
  344.   EXIT
  345. END
  346. s = SIGN(x)
  347. x = ABS(x)
  348. x = (((((((-.0012624911*x + .0066700901)*x - .0170881256)*x + ,
  349.            .0308918810)*x - .0501743046)*x + .0889789874)*x - ,
  350.            .2145988016)*x + 1.570796305)*sqrt(1 - x)
  351. RETURN s*(1.57079632679 - x)
  352. /*====================================================================*/
  353. /* acos(x).  Method:  use acos(x) = pi/2 - asin(x)                    */
  354. acos: PROCEDURE
  355. ARG x
  356. CALL argtest "acos",x
  357. IF x<-1 | x>1 THEN DO
  358.   SAY "acos argument" x " out of range."
  359.   EXIT
  360. END
  361. RETURN 1.57079632679 - asin(x)
  362. /*====================================================================*/
  363. /* atan(x).  Method:  use atan(x) = pi/4 - atan((1-x)/(1+x)).  Evaluate
  364. the latter function using a polynomial approximation, taking advantage
  365. of the fact that its argument is less than one as long as x > -1.     */
  366. atan: PROCEDURE
  367. ARG x
  368. CALL argtest "atan",x
  369. IF x = 0 THEN RETURN 0
  370. s=SIGN(x)
  371. x=ABS(x)
  372. x = (x - 1)/(x + 1)
  373. y = x*x
  374. x = ((((((((.0028662257*y - .0161657367)*y + .0429096138)*y - ,
  375.            .0752896400)*y + .1065626393)*y - .1420889944)*y + ,
  376.            .1999355085)*y - .3333314528)*y + 1)*x
  377. RETURN .785398163397 + s*x
  378. /*====================================================================*/
  379. /* acot(x).  Method:  use acot(x) = pi/2 - atan(x)                    */
  380. acot: PROCEDURE
  381. ARG x
  382. CALL argtest "acot",x
  383. RETURN 1.57079632679 - atan(x)
  384. /*====================================================================*/
  385. /* asec(x).  Method:  use asec(x) = acos(1/x)                         */
  386. asec: PROCEDURE
  387. ARG x
  388. CALL argtest "asec",x
  389. IF x>-1 & x<1 THEN DO
  390.   SAY "asec argument" x " out of range."
  391.   EXIT
  392. END
  393. RETURN acos(1/x)
  394. /*====================================================================*/
  395. /* acsc(x).  Method:  use acsc(x) = asin(1/x)                         */
  396. acsc: PROCEDURE
  397. ARG x
  398. CALL argtest "acsc",x
  399. IF x>-1 & x<1 THEN DO
  400.   SAY "acsc argument" x " out of range."
  401.   EXIT
  402. END
  403. RETURN asin(1/x)
  404. /*====================================================================*/
  405. /* sinh(x) .                                                          */
  406. sinh: PROCEDURE
  407. ARG x
  408. CALL argtest "sinh",x
  409. RETURN .5*(exp(x) - exp(-x))
  410. /*====================================================================*/
  411. /* cosh(x) .                                                          */
  412. cosh: PROCEDURE
  413. ARG x
  414. CALL argtest "cosh",x
  415. RETURN .5*(exp(x) + exp(-x))
  416. /*====================================================================*/
  417. /* tanh(x) .                                                          */
  418. tanh: PROCEDURE
  419. ARG x
  420. CALL argtest "tanh",x
  421. RETURN (exp(x) - exp(-x))/(exp(x) + exp(-x))
  422. /*====================================================================*/
  423. /* gamma(x).  Method:  use
  424.    gamma(x) = (x - 1)(x - 2)...(x - n)gamma(x - n + 1)  if x >= 1, or
  425.    gamma(x) = gamma(x + n + 1)/x(x + 1)(x + 2)...(x + n)  if x < 1,
  426. in either case, choosing n so that the argument of the gamma function
  427. on the right hand side is in the range 0 through 1.  Evaluate this
  428. gamma function using a polynomial approximation.                      */
  429. gamma: PROCEDURE
  430. ARG x
  431. CALL argtest "gamma",x
  432. p = 1
  433. SELECT
  434. WHEN x >= 1 THEN
  435.   DO WHILE x >= 2
  436.     x = x - 1
  437.     p = x*p
  438.   END
  439. WHEN x%1 = x THEN DO
  440.   SAY "gamma argument" x "out of range."
  441.   EXIT
  442.   END
  443. OTHERWISE
  444.   DO UNTIL x >= 1
  445.     IF x <= -1 & ABS(p) >= ABS(1E999999999/x) THEN RETURN 0
  446.     p = x*p
  447.     x = x + 1
  448.   END
  449.   p = 1/p
  450. END
  451. IF x \= 1 THEN DO
  452.   x = x - 1
  453.   x = (((((((.035868343*x - .193527818)*x + .482199394)*x - ,
  454.             .756704078)*x + .918206857)*x - .897056937)*x + ,
  455.             .988205891)*x - .577191652)*x + 1
  456.   p = x*p
  457. END
  458. RETURN p
  459. /*====================================================================*/
  460. /* lngamma(x), x > 0.  Method:  If x < 5, calculate ln(gamma(x)).
  461. Otherwise, use Stirling's formula thorugh x**-7.                      */
  462. lngamma: PROCEDURE
  463. ARG x
  464. CALL argtest "lngamma",x
  465. IF x < 0 THEN DO
  466.   SAY "lngamma argument" x "out of range."
  467.   EXIT
  468. END
  469. IF x < 5 THEN RETURN ln(gamma(x))
  470. RETURN (x - .5)*ln(x) - x + .91893853320 + 1/(12*x) - 1/(360*x**3) + ,
  471.        1/(1260*x**5) - 1/(1680*x**7)
  472.  
  473. /* fctrl(x) = x! = gamma(x + 1):                                      */
  474. fctrl: PROCEDURE
  475. ARG x
  476. CALL argtest "fctrl",x
  477. IF x%1 = x & x < 0 THEN DO
  478.   SAY "fctrl argument" x "out of range."
  479.   EXIT
  480. END
  481. RETURN gamma(x + 1)
  482. /*====================================================================*/
  483. /* erf(x).  One or another polynomial approximation is used, the choice
  484. depending upon the magnitude of x.  (CACM Algorithm 209.)             */
  485. erf: PROCEDURE
  486. ARG x
  487. CALL argtest "erf",x
  488. s = SIGN(x)
  489. x = .707106781187*ABS(x)
  490. IF x >= 3 | x = 0 THEN RETURN s
  491. IF x < 1 THEN DO
  492.   w = x*x
  493.   x = ((((((((.000124818987*w - .001075204047)*w + .005198775019)*w - ,
  494.              .019198292004)*w + .059054035642)*w - .151968751364)*w + ,
  495.              .319152932694)*w - .531923007300)*w + .797884560593)*2*x
  496.   END
  497. ELSE DO
  498.   x = x - 2
  499.   x = (((((((((((((-.000045255659*x + .000152529290)*x - ,
  500.              .000019538132)*x - .000676904986)*x + .001390604284)*x - ,
  501.              .000794620820)*x - .002034254874)*x + .006549791214)*x - ,
  502.              .010557625006)*x + .011630447319)*x - .009279453341)*x + ,
  503.              .005353579108)*x - .002141268741)*x + .000535310849)*x + ,
  504.              .999936657524
  505.   END
  506. RETURN s*x
  507. /*====================================================================*/
  508. /* erfc(x) = 1 - erf(x).                                              */
  509. erfc: PROCEDURE
  510. ARG x
  511. CALL argtest "erfc",x
  512. RETURN 1 - erf(x)
  513. /*====================================================================*/
  514. argtest: PROCEDURE
  515. PARSE ARG name,x
  516. IF DATATYPE(x,"N") THEN RETURN
  517. SAY name "argument" x "not a number."
  518. EXIT
  519. /*====================================================================*/
  520. getvar: PROCEDURE
  521. PARSE UPPER ARG varName
  522.    PARSE UPPER SOURCE env .
  523.    rt = ""
  524.    IF env = 'CMS' THEN
  525.       DO
  526.       'GLOBALV STACK' varName
  527.       PULL rt
  528.       END
  529.    ELSE IF env = 'OS/2' THEN
  530.       rt = value(varName,,'OS2ENVIRONMENT')
  531. RETURN rt
  532. /*====================================================================*/
  533. setvar: PROCEDURE
  534. PARSE UPPER ARG varName,value
  535.    PARSE UPPER SOURCE env .
  536.    rt = ""
  537.    IF env = 'CMS' THEN
  538.       'GLOBALV SET' varName value
  539.    ELSE IF env = 'OS/2' THEN
  540.       rt = value(varName,value,'OS2ENVIRONMENT')
  541. RETURN rt
  542. SYNTAX:
  543.    say '"'evx_expr'" is not quite valid syntax.  Try again!'
  544. exit
  545. /* END */
  546.