From: jbrock@csfb1.fir.fbc.com (John Brock)
Subject: RXMATH math package for REXX
Reply-To: uunet!csfb1!jbrock
Organization: First Boston Corporation
Date: Thu, 8 Oct 1992 06:40:38 GMT

Here is the code for my RXMATH math package.  Any comments
or questions will be welcome.  I am particularly interested
in UNIX or OS/2 ports.  Please let me know if you find any
hidden CMS dependencies!  (I wonder, for example, about my
use of "parse source" to determine whether RXMATH is running
as a command or a function call).


-------- Beginning of RXMATH math package for REXX ----------------------

This package contains three files:  1) RXMATH EXEC
                                    2) RXMATH$ EXEC
                                    3) RXMATH HELPREXX

The first two files are REXX execs; the third is a CMS style help file.
The execs were developed under CMS, and have run without modification
under TSO.  I have not tried to run them using OS/2 or UNIX or PC REXX
interpreters, and would be very interested to hear what was involved.
My installation's CMS level is VM/XA CMS 5.6, and my REXX version is
REXX370 3.40.

RXMATH is the main exec.  It handles all calculations, and may be
compiled.  RXMATH$ exists only to support the "calc" function, and may
not be compiled, as it contains "interpret" statements.

RXMATH provides arbitrary precision math functions, and must be run under
a version of REXX which supports arbitrary precision arithmetic.  RXMATH
and RXMATH$ use the REXX stack to communicate, so this is required if
you are going to use the "calc" function.

RXMATH is offered as unsupported freeware.  I only ask, if you modify it,
that you leave my name in it somewhere.

And while I really can't promise any kind of support, I would certainly
be inclined to fix any bugs that were pointed out to me, or answer any
questions.  I can currently be reached at the addresses below.  Please
note that while I can send mail to Usenet address from my installation
I cannot send mail to Internet addresses.


John Brock
uunet!csfb1!jbrock
jbrock@csfb1.fir.fbc.com
October 7, 1992


-------- Cut here -------- Beginning of RXMATH EXEC ---------------------
/* REXX -- Written by John Brock, 09/08/89, CMS/TSO version 09/07/90.
 
RXMATH:     Arbitrary precision mathematical functions for REXX
 
Syntax:     Command:        EXEC RXMATH <precision> function <x <y>>
            REXX function:  RXMATH(<precision,> function <,x <,y>>)
 
Where:      precision  is the desired REXX precision of the result
                         (if omitted or zero precision defaults to 9)
            function   is the name of a mathematical function
            x, y       are arguments to the function
 
Functions:  FACT(x)    --  Factorial of x
            PERM(x,y)  --  Permutations of x by y
            COMB(x,y)  --  Combinations of x by y
            SQRT(x)    --  Square root of x
            POW(x,y)   --  x to the power of y
            LOG(x,y)   --  Log of y base x
            EXP(x)     --  e to the power of x
            LN(x)      --  Natural log of x
            PI()       --  Value of pi
            SIN(x)     --  Sine of x
            COS(x)     --  Cosine of x
            TAN(x)     --  Tangent of x
            COT(x)     --  Cotangent of x
            SEC(x)     --  Secant of x
            CSC(x)     --  Cosecant of x
            ARCSIN(x)  --  Inverse sine of x
            ARCCOS(x)  --  Inverse cosine of x
            ARCTAN(x)  --  Inverse tangent of x
            ARCCOT(x)  --  Inverse cotangent of x
            ARCSEC(x)  --  Inverse secant of x
            ARCCSC(x)  --  Inverse cosecant of x
            CALC(x)    --  Calculate the value of an expression
                             (e.g. '1 + exp(pi() / 2)')
 
Notes:      All angles for the trigonometric functions are in radians.
            Calling RXMATH as a command without arguments prints help.
            When called as a command result is printed to screen,
              and on error RXMATH returns 1.
            When called as a function result is returned to caller,
              and on error RXMATH exits without returning a value.
            When the CALC function is called without an argument
              RXMATH goes into an interpretive REXX calculator mode
              (change or drop REXX variable "$prompt" to alter prompt).
 
*/
 
signal on halt /* First instruction, because we could halt anytime! */
signal on syntax /* Mainly in case we hit internal numeric limits. */
signal on novalue /* Shouldn't ever happen. */
 
parse source . $cmd? $pname . /* Get program name. */
$cmd? = $cmd? = "COMMAND" /* Were we called as a command? */
parse numeric $dgt . 1 $.0dgt . /* Default for precision argument. */
$extra_dgt? = 1 /* Increase internal precision above requested value? */
 
/* Get precision, function name, and arguments from argument list. */
/* Note that the "$extra_dgt?" variable may be taken from argument */
/* list.  This is an "undocumented interface", which is used to    */
/* avoid a second and redundant increase in internal precision     */
/* when RXMATH calls itself recursively through the CALC function. */
if $cmd? /* Do basic argument parsing. */
  then parse arg $func $1 $2 . 1 . $str /* ($str is for CALC). */
  else parse arg $func 1 $$extra_dgt? $$dgt, $1 1 $str, $2
select /* Deal with possible non-basic argument parsing. */
  when $cmd?
    then if datatype($func, "N") /* First argument is precision. */
      then parse arg $dgt $func $1 $2 . 1 . . $str
      else if $func = ""
        then signal help /* No arguments, so print help message. */
  when ($$extra_dgt? == 0 | $$extra_dgt? == 1) & datatype($$dgt, "N")
    then parse arg $extra_dgt? $dgt, $func, $1 1 $str, $2
  otherwise if datatype($func, "N") /* First argument is precision. */
    then parse arg $dgt, $func, $1 1 $str, $2
  end
$func = space(translate($func))
$dgt = $dgt / 1
 
/* Validate precision argument. */
if ^datatype($dgt, "W") | $dgt < 0
  then call error "Illegal value for precision:" $dgt
 
/* List of all functions supported by RXMATH (and available in CALC). */
$funcs = "FACT PERM COMB SQRT POW LOG EXP LN PI SIN COS TAN",
  "COT SEC CSC ARCSIN ARCCOS ARCTAN ARCCOT ARCSEC ARCCSC CALC"
 
/* "SCAN" and "PULLSCAN" are special functions which support CALC. */
/* They are "undocumented", and not listed in the function list    */
/* (which means that they are not available as CALC functions).    */
/* Note that $dgt = 0 has a special meaning for these functions.   */
if ^$cmd? & ($func = "SCAN" | $func = "PULLSCAN") then do
  if $func = "PULLSCAN" then do
    if $2 ^= "LIT" & queued() + externals() = 0 then say $1
    parse pull $1
    end
  exit scan($pname, $funcs, $extra_dgt?, $dgt, $1)
  end
 
/* Validate function name. */
if $func = "" then call error "Function name is missing."
if words($func) ^= 1 | find($funcs, $func) = 0
  then call error "Invalid function name:" $func
$.0func = $func /* Save function name for use in error messages. */
 
/* If precision is zero then use default value. */
if $dgt = 0 then $dgt = $.0dgt
 
/* Remove excess precision from arguments. */
numeric digits $dgt
if datatype($1, "N")
  then $1 = $1 / 1
  else $1 = space($1)
if datatype($2, "N")
  then $2 = $2 / 1
  else $2 = space($2)
 
numeric digits intdgt($dgt, $extra_dgt?) /* Set internal precision. */
 
select /* Call appropriate function. */
  when $func == "FACT" then $0 = fact($1)
  when $func == "PERM" then $0 = perm($1, $2)
  when $func == "COMB" then $0 = comb($1, $2)
  when $func == "SQRT" then $0 = sqrt($1)
  when $func == "POW" then $0 = pow($1, $2)
  when $func == "LOG" then $0 = log($1, $2)
  when $func == "EXP" then $0 = exp($1)
  when $func == "LN" then $0 = ln($1)
  when $func == "PI" then $0 = pi()
  when $func == "SIN" then $0 = sin($1)
  when $func == "COS" then $0 = cos($1)
  when $func == "TAN" then $0 = tan($1)
  when $func == "COT" then $0 = cot($1)
  when $func == "SEC" then $0 = sec($1)
  when $func == "CSC" then $0 = csc($1)
  when $func == "ARCSIN" then $0 = arcsin($1)
  when $func == "ARCCOS" then $0 = arccos($1)
  when $func == "ARCTAN" then $0 = arctan($1)
  when $func == "ARCCOT" then $0 = arccot($1)
  when $func == "ARCSEC" then $0 = arcsec($1)
  when $func == "ARCCSC" then $0 = arccsc($1)
  /* The interpretive CALC function is supported by RXMATH$,      */
  /* leaving RXMATH eligible for REXX compilation (if available). */
  when $func == "CALC" then do
    call "RXMATH$" $str, $extra_dgt? & $str = "",,
      intdgt($dgt, $extra_dgt? & $str ^= ""), errmsg(),,
      $pname "interactive mode -- enter any valid REXX instruction:"
    if symbol("result") = "VAR" /* Process output from CALC function. */
      then if $str = "" /* If we were in interactive CALC mode... */
        then exit result /* ...then return result immediately. */
        else $0 = result /* Process non-interactive result further. */
      else if $str = ""
        then exit /* A simple exit is OK from interactive mode... */
        else call error /* ...but otherwise means there was an error. */
    end
  otherwise call error "Function name is valid but not implemented."
  end
 
numeric digits $dgt /* Reset precision to requested value. */
$0 = $0 / 1 /* Format result, removing trailing zeros. */
 
/* All done, so return result to user! */
if ^$cmd? then exit $0
say $0
exit
 
 
/* Given a string, scan it for function calls that can be executed */
/* by RXMATH, and replace them with direct calls to RXMATH.  For   */
/* example, "exp(2)" might become "RXMATH(0 digits(), 'EXP', 2)".  */
/* Note that the calls to RXMATH that are constructed here supply  */
/* $extra_dgt? as part of the first argument.  This is really the  */
/* only place where this should ever happen.                       */
scan: procedure expose $.
parse arg $pname, $funcs, $extra_dgt?, $dgt, $str
if $dgt = 0 then $dgt = "digits()"
 
do $i = 1 to words($funcs) /* Loop thru list of supported functions. */
  $func = translate(word($funcs, $i)) /* Select a function. */
  $pos = 0 /* Start at the beginning of the string. */
  do forever /* Look for multiple occurrences. */
    $prepstr = prepstr($str, $pos + 1, '"' "'", '"' "'")
    $pos = pos($func"(", $prepstr, $pos + 1) /* Needn't find ")". */
    if $pos = 0 then leave /* Function not found in string. */
    if $pos > 1
      then if symbol(substr($prepstr, $pos - 1, 1)) ^= "BAD"
        then iterate /* (What we found was part of another function!) */
    /* Supported function found, so replace it with call to RXMATH. */
    $str = insert($pname"("$extra_dgt? $dgt", '"$func"', ",,
      delstr($str, $pos, length($func) + 1), $pos - 1)
    end
  end
 
return $str
 
 
/* Take a string, uppercase it, and blank out any quoted       */
/* substrings (to avoid unnecessary function call expansions). */
/* Input is string, starting position, left delimiters, and    */
/* corresponding right delimiters.                             */
/* Example: "ln(2) + calc('pi()')" ==> "ln(2) + calc('    ')". */
prepstr: procedure expose $.
arg $str, $start, $lftdlm, $rgtdlm
 
do while $start <= length($str)
 
  /* Find the nearest left delimiter, if any. */
  parse value (length($str) + 1) with $lftpos $lft $rgt
  do $i = 1 to words($lftdlm)
    $pos = pos(word($lftdlm, $i), $str, $start)
    if $pos > 0 & $pos < $lftpos then do
      $lft = word($lftdlm, $i)
      $rgt = word($rgtdlm, $i)
      $lftpos = $pos
      end
    end
 
  if $lft = "" then leave /* No left delimiter found. */
 
  /* Blank out string between left and right delimiter (if any). */
  $lftpos = $lftpos + length($lft)
  $rgtpos = pos($rgt, $str, $lftpos)
  if $rgtpos = 0 then $rgtpos = length($str) + 1
  $str = overlay(copies(" ",  $rgtpos - $lftpos), $str, $lftpos)
  $start = $rgtpos + length($rgt)
  end
 
return $str
 
 
/* Determine how numeric digits should be set internally.           */
/* If $extra_dgt? = 1 then we double requested precision, using     */
/* the default setting of numeric digits to set a lower limit (so   */
/* that a requested precision of 2, for example, won't result in an */
/* unreasonably low internal precision of 4).  If $extra_dgt? = 0   */
/* then we just use the requested precision, whatever it is.        */
intdgt: procedure expose $.
arg $dgt, $extra_dgt?
if $extra_dgt?
  then return 2 * max($dgt, $.0dgt) /* Extended precision. */
  else return $dgt /* Unextended precision. */
 
 
/* Find value of $x!. */
fact: procedure expose $.
 
$x = arg_ok(arg(1), "W >=0")
 
$y = 1
do $n = 2 to $x
  $y = $y * $n
  end
 
return $y
 
 
/* Find the number of permutations of $x things taken $y at a time. */
perm: procedure expose $.
 
$x = arg_ok(arg(1), "W >=0")
$y = arg_ok(arg(2), "W >=0")
if $x < $y then call error "Illegal arguments:" $x "<" $y
 
$z = 1
do $n = ($x - $y) + 1 to $x
  $z = $z * $n
  end
 
return $z
 
 
/* Find the number of combinations of $x things taken $y at a time. */
comb: procedure expose $.
 
$x = arg_ok(arg(1), "W >=0")
$y = arg_ok(arg(2), "W >=0")
if $x < $y then call error "Illegal arguments:" $x "<" $y
 
if $x - $y < $y then $y = $x - $y /* Adjust $y. */
 
$z = 1
do $n = ($x - $y) + 1 to $x
  $z = $z * $n
  end
 
do $n = 2 to $y
  $z = $z / $n
  end
 
return $z
 
 
/* Find square root using Newton's method. */
sqrt: procedure expose $.
 
$x = arg_ok(arg(1), ">=0")
 
if $x = 0 then return 0
if $x = 1 then return 1
 
/* Adjust number so it is between 1 and 100. */
$i = ilog(100, $x, 1, 100)
$x = $x * 100 ** -$i
 
$y = $x / 2
$old = $y
do $n = 1
  $y = (($x / $y) + $y) / 2
  if $y = $old then leave
  $old = $y
  end
 
return $y * 10 ** $i
 
 
/* Find $x to the power of $y. */
pow: procedure expose $.
 
$x = arg_ok(arg(1))
$y = arg_ok(arg(2))
 
if $y = 0 then return 1
if datatype($y, "W") & $x ^= 0 then return $x ** $y
if $x <= 0 then do
  if $x < 0 then call error,
    "Non-integer power of a negative number:" "("$x") **" $y
  if $y < 0 then call error "Zero taken to a negative power:" $y
  return 0
  end
 
/* If possible use square root function (for speed). */
if abs($y // 1) = 0.5 then return sqrt($x) ** sign($y) * $x ** ($y % 1)
 
return exp($y * ln($x))
 
 
/* Find the log of $y base $x. */
log: procedure expose $.
return ln(arg(2)) / ln(arg_ok(arg(1), "^=1"))
 
 
/* Find e to the power of $x. */
exp: procedure expose $.
 
$x = arg_ok(arg(1))
 
$i = $x % 1
if abs($x - $i) > 0.5 then $i = $i + sign($x)
$x = $x - $i /* Adjust $x for quick convergence. */
 
$y = sum_e($x)
if $i ^= 0 then $y = $y * sum_e(1) ** $i
 
return $y
 
 
/* Find the natural log of $x (using sum of -((-1)**n)*(x-1)**n/n). */
ln: procedure expose $.
 
$x = arg_ok(arg(1), ">0")
 
if $x = 10 then do /* Special assist for ln(10). */
  $m = "2.30258509299404568401799145468436420760110148862877"
  $y = format($m)
  if $y ^== $m then return $y
  end
 
/* Adjust $x between 0.5 and 1.5 (for convergence). */
$i = ilog(sum_e(1), $x, 0.5, 1.5)
$x = $x * sum_e(1) ** -$i
 
$x = $x - 1
$f = -1
$y = 0
$old = $y
do $n = 1
  $f = -$f * $x
  $y = $y + $f / $n
  if $y = $old then leave
  $old = $y
  end
 
return $y + $i
 
 
/* Find e to the power of $x (using the usual series expansion).    */
/* If $x = 1 then save value of e (because it may be needed again). */
sum_e: procedure expose $.
parse arg $x
 
if $x = 1 then do /* Return hard coded value if precision allows. */
  if symbol("$.0e") = "VAR" then return $.0e
  $e = "2.71828182845904523536028747135266249775724709369996"
  $.0e = format($e)
  if $.0e ^== $e /* Does hard coded value exceed internal precision? */
    then return $.0e
  end
 
$y = 1
$f = 1
$old = $y
do $n = 1
  $f = $f * $x / $n
  $y = $y + $f
  if $y = $old then leave
  $old = $y
  end
 
if $x = 1 then $.0e = $y /* Save value of e. */
 
return $y
 
 
/* Return $i such that $bot <= ($y / ($x ** $i)) <= $top.         */
/* Conditions are $x, $y, $top, $bot > 0, $x ^= 1, $top >= $bot   */
/* This procedure works by repeatedly looking for factors of the  */
/* form $x ** (2 ** $j) which are divided out of $y to move it    */
/* toward the intended range (without overshooting, if possible). */
/* Note: if $x > 1 then ilog($x, $y, 1, $x) = trunc(log($x, $y)). */
ilog: procedure expose $.
arg $x, $y, $bot, $top
 
$high = $y > $top
$sign = 1 - 2 * ($high ^= ($x > 1))
 
$i = 0
do while $high & $y > $top | ^$high & $y < $bot
 
  $f = $x
  do $j = -1
    $q = $y * $f ** -$sign
    if $j >= 0 & ($high & $q < $top | ^$high & $q > $bot) then leave
    $f = $f * $f
    $g = $q
    end
 
  $y = $g
  $i = $i + $sign * 2 ** $j
  end
 
return $i
 
 
/* Calculate (and save) pi, using arcsin(1/2) = pi/6. */
/* Saved or hard coded value of pi used if possible.  */
pi: procedure expose $.
 
if symbol("$.0pi") = "VAR" then return $.0pi
 
$pi = "3.14159265358979323846264338327950288419716939937511"
$.0pi = format($pi)
 
if $.0pi == $pi then $.0pi = 6 * arcsin(0.5)
 
return $.0pi
 
 
/* Find the sine of $x (usual series expansion). */
sin: procedure expose $.
 
$x = arg_ok(arg(1), "TRIG")
 
$x = $x // (2 * pi()) /* Adjust between pi and -pi, for convergence. */
if abs($x) > pi() then $x = $x - sign($x) * 2 * pi()
 
$f = $x
$y = $x
$x = $x ** 2
$old = $y
do $n = 2 by 2
  $f = -$f * $x / ($n * ($n + 1))
  $y = $y + $f
  if $y = $old then leave
  $old = $y
  end
 
return $y
 
 
/* Find the cosine of $x (usual series expansion). */
cos: procedure expose $.
 
$x = arg_ok(arg(1), "TRIG")
 
$x = $x // (2 * pi()) /* Adjust between pi and -pi, for convergence. */
if abs($x) > pi() then $x = $x - sign($x) * 2 * pi()
 
$f = 1
$y = 1
$x = $x ** 2
$old = $y
do $n = 2 by 2
  $f = -$f * $x / ($n * ($n - 1))
  $y = $y + $f
  if $y = $old then leave
  $old = $y
  end
 
return $y
 
 
/* Find the tangent of $x. */
tan: procedure expose $.
 
$x = arg_ok(arg(1))
 
$y = cos($x)
if $y = 0 then call error "Result is infinite."
 
return sin($x) / $y
 
 
/* Find the cotangent of $x. */
cot: procedure expose $.
 
$x = arg_ok(arg(1))
 
$y = sin($x)
if $y = 0 then call error "Result is infinite."
 
return cos($x) / $y
 
 
/* Find the secant of $x. */
sec: procedure expose $.
 
$x = arg_ok(arg(1))
 
$y = cos($x)
if $y = 0 then call error "Result is infinite."
 
return 1 / $y
 
 
/* Find the cosecant of $x. */
csc: procedure expose $.
 
$x = arg_ok(arg(1))
 
$y = sin($x)
if $y = 0 then call error "Result is infinite."
 
return 1 / $y
 
 
/* Find the inverse sine of $x (usual series expansion). */
arcsin: procedure expose $.
 
$x = arg_ok(arg(1), "ABS<=1")
 
/* Avoid region where series converges slowly.  (Note recursion!) */
if abs($x) >= 0.75 then return sign($x) * arccos(sqrt(1 - $x ** 2))
 
$f = $x
$y = $x
$x = $x ** 2
$old = $y
do $n = 2 by 2
  $f = $f * $x * ($n - 1) / $n
  $y = $y + $f / ($n + 1)
  if $y = $old then leave
  $old = $y
  end
 
return $y
 
 
/* Find the inverse cosine of $x. */
arccos: procedure expose $.
return pi() / 2 - arcsin(arg(1))
 
 
/* Find the inverse tangent of $x. */
arctan: procedure expose $.
 
$x = arg_ok(arg(1))
 
return arcsin($x / sqrt($x ** 2 + 1))
 
 
/* Find the inverse cotangent of $x. */
arccot: procedure expose $.
 
$x = arg_ok(arg(1))
 
return arccos($x / sqrt($x ** 2 + 1))
 
 
/* Find the inverse secant of $x. */
arcsec: procedure expose $.
 
$x = arg_ok(arg(1), "ABS>=1")
 
return sign($x) * arccos(1 / $x)
 
 
/* Find the inverse cosecant of $x. */
arccsc: procedure expose $.
 
$x = arg_ok(arg(1), "ABS>=1")
 
if $x > 0
  then return arcsin(1 / $x)
  else return -(pi() + arcsin(1 / $x))
 
 
/* Validate function arguments against a given lists of tests. */
arg_ok: procedure expose $.
parse arg $x, $test_list
upper $test_list
 
/* Argument must always exist and be a valid number. */
if $x = "" then call error "Missing argument."
if ^datatype($x, "N") then call error "Illegal argument:" $x
 
if find($test_list, "W") > 0 /* Must be whole number. */
  then if ^datatype($x, "W")
    then call error "Illegal argument:" $x
if find($test_list, ">=0") > 0 /* Must be >= zero. */
  then if $x < 0
    then call error "Argument < 0:" $x
if find($test_list, ">0") > 0 /* Must be > zero. */
  then if $x <= 0
    then call error "Argument <= 0:" $x
if find($test_list, "^=1") > 0 /* Must not equal one. */
  then if $x = 1
    then call error "Illegal argument:" $x
if find($test_list, "ABS<=1") > 0 /* Absolute value must be <= one. */
  then if abs($x) > 1
    then call error "Argument is > 1 or < -1:" $x
if find($test_list, "ABS>=1") > 0 /* Absolute value must be >= one. */
  then if abs($x) < 1
    then call error "Argument is < 1 and > -1:" $x
 
/* To accurately calculate trigonometric functions the absolute value */
/* of the argument must be small enough that we can subtract from it  */
/* a multiple of 2 * pi and get an accurate result between pi and     */
/* -pi.  Actually the test used here is too rigorous.  In most (but   */
/* not all) cases we get a completely accurate result by using        */
/* digits() % 2 (instead of digits() % 4).  But then, who really      */
/* needs sine or cosine of very large numbers anyway?                 */
if find($test_list, "TRIG") > 0
  then if abs($x) >= 10 ** (digits() % 4)
    then call error "Argument is too large:" $x
 
return $x
 
 
halt: call error errortext(4)
 
 
syntax: call error errortext(rc)
 
 
novalue: call error "NOVALUE at line" sigl
 
 
/* Print out leading comment as help. */
help:
do $i = 2 to sourceline() while sourceline($i) ^= "*/"
  say sourceline($i)
  end
exit
 
 
/* Print error message and exit. */
error: procedure expose $.
if arg(1, "E") then say errmsg(arg(1))
parse source . $cmd? .
if $cmd? = "COMMAND"
  then exit 1 /* Error return code is 1 when called as a command. */
  else exit   /* If called as function caller will blow up. */
 
 
/* Create an error message. */
errmsg: procedure expose $.
parse source . . $pname .
$errmsg = "Error in" $pname
if symbol("$.0func") = "VAR" then $errmsg = $errmsg "("$.0func")"
$errmsg = $errmsg "--"
if arg(1, "E") then $errmsg = $errmsg arg(1)
return $errmsg
 
 
/* (Delete when built-in digits() function becomes available). */
digits: procedure
parse numeric $digits .
return $digits
-------- Cut here -------- End of RXMATH EXEC ---------------------------


-------- Cut here -------- Beginning of RXMATH$ EXEC --------------------
/* REXX -- Written by John Brock, 06/01/90, CMS/TSO version 09/07/90.
 
RXMATH$ supports and should only be called by RXMATH.
 
*/
 
signal on halt
signal on syntax
 
numeric digits arg(3)
 
if arg(1) ^= ""
  then interpret "exit ("rxmath(arg(2) arg(3), "SCAN", arg(1))") + 0"
 
if arg(5, "E") then $prompt = arg(5)
 
retry: if arg(1) ^= "" then exit
 
do forever
  interpret rxmath(arg(2) 0, "PULLSCAN", $prompt, symbol("$prompt"))
  end
 
 
halt: signal on halt
say arg(4) errortext(4)
signal retry
 
 
syntax: signal on syntax
if rc ^= 44 then say arg(4) errortext(rc)
signal retry
 
 
/* (Delete when built-in digits() function becomes available). */
digits: procedure
parse numeric $digits .
return $digits
-------- Cut here -------- End of RXMATH$ EXEC --------------------------


-------- Cut here -------- Beginning of RXMATH HELPREXX -----------------
.cm RXMATH HELPREXX -- Written by John Brock, 09/08/89, version 09/07/90.
RXMATH
 
Use the RXMATH exec to do arbitrary precision mathematics.
 
RXMATH may be called as a command from the command line, in which case
the result is printed to the screen, or it may be called as a function
returning a value to REXX.  If an error occurs RXMATH prints an error
message to the screen.  Then if called as a command it exits with a
return code of 1; otherwise it exits without returning a value.
 
When issued as a command, RXMATH has the format:
 
+--------------------------------------------------------------------+
| RXMATH   | <precision> function <x <y>>                            |
+--------------------------------------------------------------------+
 
When called as a function from REXX, RXMATH has the format:
 
+--------------------------------------------------------------------+
| RXMATH(<precision,> function <,x <,y>>)                            |
+--------------------------------------------------------------------+
 
Where:      precision  is the desired REXX precision of the result
                         (if omitted or zero precision defaults to 9)
 
            function   is the name of a mathematical function
 
            x, y       are arguments to the function
 
Functions:  FACT(x)    --  Factorial of x
            PERM(x,y)  --  Permutations of x by y
            COMB(x,y)  --  Combinations of x by y
            SQRT(x)    --  Square root of x
            POW(x,y)   --  x to the power of y
            LOG(x,y)   --  Log of y base x
            EXP(x)     --  e to the power of x
            LN(x)      --  Natural log of x
            PI()       --  Value of pi
            SIN(x)     --  Sine of x
            COS(x)     --  Cosine of x
            TAN(x)     --  Tangent of x
            COT(x)     --  Cotangent of x
            SEC(x)     --  Secant of x
            CSC(x)     --  Cosecant of x
            ARCSIN(x)  --  Inverse sine of x
            ARCCOS(x)  --  Inverse cosine of x
            ARCTAN(x)  --  Inverse tangent of x
            ARCCOT(x)  --  Inverse cotangent of x
            ARCSEC(x)  --  Inverse secant of x
            ARCCSC(x)  --  Inverse cosecant of x
            CALC(x)    --  Calculate the value of an expression
                             (e.g. '1 + exp(pi() / 2)')
 
Examples:
 
1. Find the square root of 2 to 40 places:
 
       EXEC RXMATH 40 SQRT 2  (or EXEC RXMATH 40 POW 2 0.5)
 
   result: 1.41421356237309504880168872420969807857
 
2. Find the square root of 2 to 9 places (the default):
 
       EXEC RXMATH SQRT 2
 
   result: 1.41421356
 
3. From REXX, print the base 10 log of 5 to 40 places:
 
       say "RXMATH"(40, "LOG", 10, 5)
 
   result: 0.6989700043360188047862611052755069732318
 
4. From REXX, print the base 10 log of 5 to 9 places:
 
       say "RXMATH"("LOG", 10, 5)
 
   result: 0.698970004
 
5. From REXX, print 1 plus e to the power of pi/2, to 50 places:
 
       say "RXMATH"(50, "CALC", "1 + exp(pi() / 2)")
 
   result: 5.8104773809653516554730356667038331263901708746645
 
6. Enter interactive REXX calculator mode:
 
       EXEC RXMATH CALC
 
   result: RXMATH interactive mode -- enter any valid REXX instruction:
           a = 1 + exp(pi() / 2)
           RXMATH interactive mode -- enter any valid REXX instruction:
           say a
           5.8104774
           RXMATH interactive mode -- enter any valid REXX instruction:
           numeric digits 50
           RXMATH interactive mode -- enter any valid REXX instruction:
           b = 1 + exp(pi() / 2); say b
           5.8104773809653516554730356667038331263901708746648
           RXMATH interactive mode -- enter any valid REXX instruction:
           exit 321
           Ready(00321);
 
Usage notes:
 
1. All angles for the trigonometric functions are expressed in radians
   (1 radian = 180/pi degrees).
 
2. In general the numeric arguments to an RXMATH function may be any
   legal REXX numbers, of any size and precision.  If an argument has
   more significant digits than the requested precision of the result
   the argument will be rounded before use.  Of course an attempt to use
   a number outside the range of a given function (e.g., SQRT(-2) or
   FACT(3.3)) will result in an error.  Arguments to the trigonometric
   functions may not be of arbitrary size, but must have an absolute
   value which is < 10 ** (precision % 2).  So for example using
   precision = 9 you could calculate sin(9999), but trying to calculate
   sin(10000) would result in an error.  In the unlikely event that you
   really need to calculate the sine or cosine of a very large number,
   it can be done by increasing the size of the of the precision
   argument until the above condition is met (which may however result
   in an unacceptably long execution time).
 
3. RXMATH is extremely accurate.  If RXMATH returns a number with
   500 significant digits, all 500 should be correct, and the last
   digit should be rounded correctly.  The only exception is when the
   CALC function is given an expression, such as "sin(pi())", which ought
   to equal zero exactly, but instead returns a very small non-zero
   number (such as 6.96470695E-24).  This is unavoidable, since
   any rounding error at all will always be larger than and thus
   dominate a target result of zero.  And of course an expression
   like ((exp(99999) + 1) - exp(99999)) will be evaluated as zero
   rather than one, for an obvious and related reason.
 
4. Some performance considerations.  When the precision is <= 25
   RXMATH can use hard coded values for e and pi (which are used
   internally a lot), and ln(10) (which provides a special boost for
   log base 10 calculations).  The POW function knows when it is being
   asked to take a square root (e.g., RXMATH("POW", 0.25, -1.5) = 8),
   and uses the considerably faster SQRT function when appropriate.
   And RXMATH can recognize situations that REXX can deal with directly,
   such as integer powers, so that for example RXMATH("POW", 2.5, 10)
   will execute much faster than RXMATH("POW", 2.5, 10.1).
 
5. The CALC function takes a character string and interprets it as a
   REXX expression, using the RXMATH math functions with names and
   arguments as listed above (e.g., "exp(4)" is replaced by the
   appropriate call to RXMATH itself).  In addition to the math functions
   supplied by RXMATH the expression may contain any REXX built-in
   functions or any external functions (EXECs or MODULEs) which return
   a numeric result.
 
6. If you call the CALC function without any arguments you will go into
   RXMATH interactive mode.  From there you can issue any valid REXX
   instructions, using any of the math functions supported by RXMATH,
   as well as any other functions or commands usable from within REXX.
   When you enter interactive mode you will be prompted for input.  The
   prompt string is contained in the variable "$prompt", and the prompt
   may be altered or suppressed by changing or dropping that variable.
   In addition if there are lines in the system stack the prompt is
   suppressed while those lines are read and executed.  The "return" or
   "exit" instructions will take you out of interactive mode, and may be
   used to return a value if desired.
 
7. The most accurate way to use RXMATH to evaluate complex expressions
   is to use the CALC function non-interactively.  The interactive CALC
   function is slightly less accurate, and gives exactly the same results
   as you would get if instead of using the CALC function you called
   RXMATH to evaluate each math function separately.  Consider the
   following four methods for calculating the expression "pi() - ln(23)",
   using the default precision of 9.
 
   1.  a = rxmath("PI") - rxmath("LN", 23) /* CALC is avoided. */
 
   2.  queue 'a = pi() - ln(23)'
       queue 'return a'
       b = rxmath("CALC") /* Interactive CALC (stacked instructions). */
 
   3.  c = rxmath("CALC", "pi() - ln(23)") /* Non-interactive CALC. */
 
   4.  queue 'c = calc("pi() - ln(23)")' /* CALC from within CALC. */
       queue 'return c'
       d = rxmath("CALC")
 
   Methods 1 and 2 are exactly equivalent, with a = b = 0.00609843.
   Method 3 is more accurate, setting c = 0.00609843766.  In the first
   two cases pi() and ln(23) are both calculated to 9 places, and their
   difference is taken, giving a result with fewer than 9 digits of
   precision.  In case 3 RXMATH calculates both pi() and ln(23) to an
   internal precision that is greater than the precision requested by
   the user, takes the difference (again with greater precision),
   and then rounds the final result to a full 9 digits of accuracy before
   returning it.  Method 4 invokes the non-interactive CALC function
   from within interactive CALC.  This may seem a bit confusing, but
   it is perfectly legitimate (CALC being an RXMATH function like any
   other), and is exactly equivalent to method 3.  Just make sure that
   the argument is a quoted string, so it isn't evaluated before it
   gets into the non-interactive CALC function.  There are also a number
   of other ways to invoke CALC from within CALC, which are amusing but
   not terribly useful.
 
-------- Cut here -------- End of RXMATH HELPREXX -----------------------
-------- End of RXMATH math package for REXX ----------------------------

-- 
John Brock (212) 909-4662
The First Boston Corporation
12 East 49th Street
New York, NY 10017
uunet!csfb1!jbrock
jbrock@csfb1.fir.fbc.com