[comp.lang.rexx] On-line calculator

HARPER@rs1.vuw.ac.nz (John Harper) (03/27/90)

I have written CALC EXEC and CALC MANUAL as described below and will be pleased
to e-mail them to anybody who asks. (Our computer has no ftp facilities.)

CALC EXEC (496 lines) is an on-line calculator in REXX for the VM/CMS operating
system. It knows the usual real elementary functions and some others.
Design principle: simplicity and brevity of program rather than fast running.
Example: CALC  3*2  LN(2)  gives  6 0.6931...  to whatever accuracy you want
(default 10 figures; over 50 is rather slow - our IBM4381 took 0.2sec cpu time
for the above to 10 figures, 2.3sec to 50, for example.)
CALC by itself waits for expressions or commands to be entered.
More detail: see CALC MANUAL file (268 lines) or enter CALC ?
By J F Harper Mathematics Dept. Victoria University Wellington NZ
e-mail  <harper@rs1.vuw.ac.nz>
    or  <harper%rs1.vuw.ac.nz@uunet.uu.net>
    or  <harper%rs1.vuw.ac.nz@relay.cs.net>
    or  uunet!vuwcomp!rs1!harper
    or  HARPER%nz.ac.vuw.rs1@ean-relay.ac.uk
    or  postmaster%harper%nz.ac.vuw.rs1@uk.ac.ean-relay
phone +64 4 721 000  or + 64 4 715 341

HARPER@rs1.vuw.ac.nz (HARPER) (03/28/90)

In article <506@rs1.vuw.ac.nz>, HARPER@rs1.vuw.ac.nz (John Harper) writes:

>I have written CALC EXEC and CALC MANUAL as described below and will be please
>to e-mail them to anybody who asks. (Our computer has no ftp facilities.)
>
I tried to. Some people who asked have e-addresses our e-mailer can't recognise
and I am spending too much of my time e-mailing things. And one person said
it was brief and relevant enough to post here in full. So here is a copy
of CALC MANUAL (268 lines); stay tuned for CALC EXEC . Note: I have REXX
access only to a VM/CMS operating system so I can't implement any changes
that might be required for Amiga or any other REXX-capable system.

                                 CALC MANUAL
                                     by
        J F Harper Mathematics Dept Victoria University Wellington NZ

        e-mail  <harper@rs1.vuw.ac.nz>
            or  <harper%rs1.vuw.ac.nz@uunet.uu.net>
            or  <harper%rs1.vuw.ac.nz@relay.cs.net>
            or  uunet!vuwcomp!rs1!harper
            or  HARPER%nz.ac.vuw.rs1@ean-relay.ac.uk
            or  postmaster%harper%nz.ac.vuw.rs1@uk.ac.ean-relay
        phone +64 4 721 000  or + 64 4 715 341

                     For CALC EXEC version of 27 Mar 1990
CONTENTS
    1. INTRODUCTION
    2. WHAT CALC CAN DO
    3. ACCURACY
    4. METHODS
    5. ENVIRONMENT-SETTING COMMANDS
    6. FILES

1. INTRODUCTION

CALC is an on-line calculator which can do many simple calculations, and some
not so simple. It is written in REXX for the IBM VM/CMS operating system.
If it seems wrong or incomprehensible, please tell me (address above).

Example 1: enter (using either lower or UPPER case)
   CALC 2+3
and the reply on the screen will be 5

If you enter CALC on a line by itself it will go on doing sums for you until
you enter a blank line (just push the RETURN key).

Example 2: enter
   CALC
   E**10
   EXP(10)
   EXP(5)**2
and you will see three answers that should be the same.

2. WHAT CALC CAN DO

It knows about REXX built-in operators e.g. +, -, *, /, % (divide and discard
remainder), // (remainder), || (concatenate strings), **integer, REXX built-in
functions e.g ABS, SIGN, (lots more in the REXX manual), also E (2.71828...),
PI (3.14159...), ANSWER (to your previous calculation) and these mathematical
functions: SQRT, FAC (integer factorial), LN (natural log), LOG10 (log to base
10), EXP, POWER, SIN, COS, TAN, ASIN, ACOS, ATAN, ARCSIN, ARCCOS, ARCTAN, SINH,
COSH, TANH, SH, CH, TH, ERF, ERFC, GAMMA, LGAMMA, BETA, VSET, VANG, DIST,
ANSWER, X, Y, Z, LAT, LON, SIZE.
Here ASIN = ARCSIN, ACOS = ARCCOS, ATAN = ARCTAN are inverse trig functions;
all trigonometry  is in radians until you enter DEG to change to degrees;
RAD changes back to radians;
POWER(a,b) is EXP(b*LN(a)) where a must be +ve but b need not be integer,
(if b is in fact an integer a**b will evaluate faster);
SINH = SH, COSH = CH, TANH = TH are hyperbolic functions;
ERF,ERFC are error functions (ERFC(0)=1,ERF(x)->1 as x->infinity);
GAMMA(n)=FAC(n-1) if n is a positive integer, LGAMMA(x) is LN(GAMMA(x)), and
BETA is the beta function: BETA(x,y)=GAMMA(x)*GAMMA(y)/GAMMA(x+y);
For GAMMA(x) and BETA(x,y) x and y must be positive.
VSET, VANG and DIST are useful for geophysical and geographical work:
VSET(lat,lon,size) gives the 3 components (x,y,z) of the vector with length
size pointing from the centre of the Earth to latitude lat and longitude lon,
VANG(x,y,z) gives lat,lon, and size (i.e. the inverse function of vset),
DIST(lat1,lon1,lat2,lon2) is the great-circle angular distance between two
points with given latitudes and longitudes. After using VSET the variables
X, Y, Z are the vector components found; after using VANG the variables
LAT, LON, SIZE are the latitude, longitude, size found.
These variables may then be used in further work. VANG and VSET leave ANSWER
unchanged. VANG, VSET and DIST all assume a spherical Earth.

Example 3: enter
   CALC
   2+3
   ANSWER*ANSWER
   ANSWER+2
and you should get 5 then 25 then 27

If you enter the filename of an EXEC followed by () or by its parameters in
parentheses, and the filename is not the name of any CALC function listed
above, and the EXEC is written in REXX, then you will find yourself running
the EXEC. If the EXEC does not return any value on exit it will run but you
will get an error message afterwards. If the EXEC is in EXEC or EXEC 2 dialect
instead of REXX it will not even run. To run an EXEC with any name in any
dialect see under CMD EXEC below.

When using non-numerical input, remember that CALC uses REXX which always
puts it in UPPER case unless you enclose the lower case letters in quotes:
c2d(a) = C2D(A) = 193 = EBCDIC code of A ; c2d('a') = 129 = EBCDIC code of a.
Input parameters which are neither letters nor numbers and are not possible
expressions must be in quotes: c2x(/) gives an error while c2x('/') = 61.
(However 3/2 is an expression, which = 1.5) For quotes, use ' ' not " " unless
you want to explore IBM oddities, e.g. c2x(""/"") = 61  but c2x("/") gives an
error if your escape character is the usual " (see TERMINAL in the CP manual).

Non-numerical output works: C2X('A') = C1 = EBCDIC hex code of A. As a result,
some "wrong" input is interpreted simply as a string of characters, and its
output is just the input with lower case letters changed to UPPER case.

To make CALC say how long each calculation took (to 0.01 sec total CPU time)
when it has finished, enter TIME (NOTIME stops this).

CALC works even if called from XEDIT at the ====> line, or from CMS SUBSET.

Entering HELP or ? gives help (by BROWSEing this CALC MANUAL file, unless you
have just made a REXX error, in which case you get that error's help message.)
After errors you are shown the current environment settings (see below).

To start writing results of CALC commands to the file CALC RESULTS A enter FILE
The result will also appear on the terminal. If CALC RESULTS A already exists,
the result will go on a new line at the end. To stop results being filed in
this way either exit from CALC by pressing the <RETURN> key or enter NOFILE
which leaves you in CALC. Every CALC session starts with the setting NOFILE.
If you enter * as the first character on a line CALC will treat it as a comment
and do nothing except copy the line into CALC RESULTS A if FILE is active.
Comments may also appear after environment-setting commands (e.g. Example 4)
except CMD but CALC will ignore them completely. It also ignores ordinary REXX
comments between /* and */
Long output lines will be written in CALC RESULTS A as several shorter lines,
with the same maximum length that your terminal uses, as given by the REXX
linesize() function and set by the CP TERMINAL LINESIZE command. Linesize()
can be 0; if so, the CALC RESULTS maximum is 255 characters.

Example 4: enter the 4 lines
   CALC 37      to force working to 37 significant figures
   FIG  33      to show answers to 33 significant figures
   FILE
   EXP(PI*SQRT(163))
and you should get 262537412640768743.999999999999250 on the screen and in the
file CALC RESULTS A.

NOTE 1: Commands causing errors will get you out of FILE mode and their output
        will appear only at the terminal.
NOTE 2: FILE leaves variables such as ANSWER unchanged.

If the first word of an input line is CMD then the rest of the line is treated
as a REXX instruction (see Examples 5,6,7 below). There may or may not be any
visible output to the terminal. If there is some, it will be filed away by a
previous FILE command only if you change SAY to CALL OUT as in Example 6.
To run a CP command enter CMD CP before it as in Example 8; to run an EXEC
enter CMD EXEC before it, e.g.  CMD EXEC DIRMAINT  or  CMD EXEC DIRMAINT ?
and to run a MODULE put CMD before it, e.g.  CMD BROWSE CALC EXEC

If any CMD command gives a nonzero return code CALC will treat it as an error
(see NOTE 1 to FILE above) even if it isn't really, e.g. CMD MAKEBUF

Example 5:  CMD X = 7*11
gives X that value (77) which can be used in later commands in that CALC
session. Warning: it overrides the meaning given to X by VSET (see above);
other variables which can be altered in this way are ANSWER E LAT LON PI SIZE
X Y and Z , array C. , and all REXX special variables, e.g. RC and SIGL .
Note that the input lines CMD X = 7*11 and X = 7*11 are quite different;
X = 7*11 would evaluate the Boolean expression (X = 77) and give the answer
1 if the previous value of X was 77, 0 if not.

Example 6:  CMD TRACE R
shows in detail how CALC works (I use it for debugging).

Example 7:  CMD DO N=1 TO 10; CALL OUT N SQRT(N); END
writes a table of numbers and their square roots from 1 to 10 (to achieve this
otherwise you would have to issue 10 separate CALC commands). If you change
CALL OUT to SAY in Example 7 your terminal output will be the same but nothing
will appear in CALC RESULTS A even if you have already entered FILE
(OUT is a procedure inside CALC that will SAY an answer, and, if FILE has been
requested, put it in CALC RESULTS also.)

Example 8:  CMD CP Q T
shows the time, day of the week, date, and how much CPU time you have used
since you logged on: CP Q T is a valid abbreviation of CP QUERY TIME

3. ACCURACY

Calculations are done with 15 significant figures and answers displayed to 10
unless you change this by entering CALC or FIG followed by a number, e.g. FIG 5
to display 5 significant figures or CALC 40 to calculate to 40 figures.

CALC uses REXX which copes with VERY large or small numbers, + or -
9.999999999E+999999999 to 1E-999999999, or 0. Your answers won't always be
good to as many digits as you might think: if f is a rapidly varying function,
f(x) can be a long way from f(x*1.00000000000001) and you may get the latter
when hoping for the former! This is inevitable in arithmetic with limited
accuracy at every step. Hence EXP(x),COSH(x),SINH(x),GAMMA(x) lose a few digits
of accuracy if x is large. Similarly, SIN and COS lose accuracy for large x
because of the error in subtracting an appropriate multiple of 2*pi or 360 deg.
The answer may even have the wrong sign if x is close to a zero of SIN or COS.
There can be enormous errors in TAN if x is close to an odd multiple of pi or
90 deg.  Also, doing arithmetic with n figures, even with exact data, is not
the same as guaranteeing function values or answers correct to n figures. As a
result, GAMMA,LGAMMA,PI,EXP,ERF,ERFC and all log, trig and hyperbolic
functions give about n-2 figures, and if the answer is very large (N figures
before the decimal point), replace n - 2 by n - 2 - log10(N).
If you want to be certain of the accuracy of your results, it is often a good
idea to increase CALC and see if anything changes.

4. METHODS:

SQRT uses the Newton-Raphson method, other elementary functions use Taylor's
theorem for small x and analytical properties to get from large to small x.
The algorithms should be mostly clear from reading the CALC EXEC source code.
LGAMMA uses the recurrence relation to get from LGAMMA(x) to LGAMMA(x+N) where
N is a suitably large integer, then the asymptotic expansion. GAMMA(x) is
found as EXP(LGAMMA(x)).
ERF uses its Taylor series, ERFC(x) the recurrence relation backwards on its
repeated integrals, unless x is large enough and the accuracy asked for is low
enough to use the asymptotic series instead, in which case it uses that.
As ERF(x) + ERFC(x) = 1 and ERF is more accurate and faster for small x, ERFC
for large x, ERF has been defined as 1-ERFC for x>3, ERFC as 1-ERF for x<2.
If x is between 2 and 3 you may check them against each other. There is also
ERFCA(x) which is always the asymptotic series and tells you if it can't give
an answer because x is too small.

NOTE: If you ask for many significant figures (e.g. 50 or more) CALC will try
to oblige, but you may run out of patience, computing funds or virtual machine
size, especially if you use gamma or error functions.

5. ENVIRONMENT-SETTING COMMANDS

There are several commands that set the environment within which CALC does its
calculations. CMD applies only to the current input line, FILE applies until
exit from CALC or until a NOFILE command, and all the others stay until you
change them: each new CALC session starts where your last session finished.
The default setting is what you get if you haven't specified otherwise.

  CALC n        Calculate to n significant figures            [default n=15]

  CMD command   Treat command as a REXX instruction, not an expression to be
                evaluated

  DEG           Trigonometry in degrees

  FIG  n        Show answers to n significant figures         [default n=10]

  FILE          Write results to the file CALC RESULTS A and to the terminal

  NOFILE        Write results to terminal only                [default]

  RAD           Trigonometry in radians                       [default]

  TIME          Show total CPU time used by each calculation

  NOTIME        Don't show CPU time used by each calculation  [default]

6. FILES

Files which are needed to use CALC but need not be on your A disk:
CALC    EXEC      which does most of the work
CALC    MANUAL    which you're reading
and the VM/CMS operating system including GLOBALV, EXECIO, BROWSE, STATE and
YHELP or HELP. If YHELP is available it is used, as it works faster.

Files which will appear on your A disk (CALC automatically updates them):
LASTING GLOBALV   which holds your settings for RAD DEG CALC FIG and TIM
CALC    RESULTS   if you use the FILE command above.

You can stack CALC commands in an EXEC and then call CALC to do them. Here
are alternative versions in REXX and EXEC 2 dialects which both evaluate the
curious number exp(pi*sqrt(163)), which is very close to an integer, and save
its value in the file CALC RESULTS A as in Example 4 above:

&TRACE                          |   /* A REXX EXEC using CALC
* An EXEC2 EXEC using CALC      |   */
&STACK CALC 37                  |   queue 'CALC 37'
&STACK FIG  33                  |   queue 'FIG  33'
&STACK FILE                     |   queue 'FILE'
&STACK EXP(PI*SQRT(163))        |   queue 'EXP(PI*SQRT(163))'
&STACK                          |   queue   /* blank line to exit */
EXEC CALC                       |   'EXEC CALC'

HARPER@rs1.vuw.ac.nz (HARPER) (03/28/90)

In article <506@rs1.vuw.ac.nz>, HARPER@rs1.vuw.ac.nz (John Harper) writes:

>I have written CALC EXEC and CALC MANUAL as described below and will be please
>to e-mail them to anybody who asks. (Our computer has no ftp facilities.)
>
As I said a minute ago when posting CALC MANUAL, some of those who asked
have e-addresses I can't send to directly, and I don't want to spend all
my time e-mailing stuff. CALC EXEC (496 lines) follows the /******/ line;
I have access to REXX only via VM/CMS so can't modify it to suit Amiga or
anything else:

/***************/
/* On-line calculator for VM/CMS op. system, knows elementary functions and
   some others, does arithmetic using REXX. Example: CALC  3*2  sqrt(2)
   gives  6 1.4142...  to whatever accuracy you specify.
   CALC by itself waits for expressions to be entered.
   More detail: see CALC MANUAL file or enter CALC ?
   Design principle: simplicity of program rather than fast running.
   By J F Harper Mathematics Dept. Victoria University Wellington NZ
   e-mail  <harper@rs1.vuw.ac.nz>
       or  <harper%rs1.vuw.ac.nz@uunet.uu.net>
       or  <harper%rs1.vuw.ac.nz@relay.cs.net>
       or  uunet!vuwcomp!rs1!harper
       or  HARPER%nz.ac.vuw.rs1@ean-relay.ac.uk
       or  postmaster%harper%nz.ac.vuw.rs1@uk.ac.ean-relay
   phone +64 4 721 000  or + 64 4 715 341
   This edition 27 Mar 1990.
*/
address command
parse arg c.input1 /* = whatever followed the word CALC on entry */
'MAKEBUF'
/* Most global internal variables are c.something */
c.error = 0 ;
call initialize
signal start
syntax:c.error = rc          /* After REXX errors restart here */
       say 'ERROR in CALC EXEC line' sigl':' errortext(c.error)
       say 'For more information on this error enter HELP or ?'
       c.input1 = 'C.INPUT1'
error: if c.error = 0 | datatype(c.error) ^= 'NUM'
       then c.error = 'C' /* After CALC errors restart here*/
       if left(c.tim,2) ^= 'C.' then call showparms
start: c.file = 'NO';
       signal on syntax
do forever /*main loop*/
   if c.tim = 'TIME'
   then c.time = format(c2d(substr(diag('C'),25,8))*1E-6,,2)
    /* After some errors REXX gives 'C.INPUT1' not '' */
   if c.input1 = ''  | c.input1 = 'C.INPUT1'
   then parse pull c.input2; else c.input2 =  c.input1
   c.cmd = translate(word(c.input2,1))
   if c.input2 = ''  | c.input2 = 'C.INPUT1' then call exitcalc
   if c.input2 = '?' | c.cmd  = 'HELP'   then
   do;if c.error = 0 | c.error = 'C' then
         do; 'BROWSE CALC MANUAL *'
             if c.error = 0 then call showparms
         end
      else
         do;'STATE YHELP EXEC *'
            if rc = 0
               then 'EXEC YHELP DMS4'err(c.error)'E'
               else '      HELP DMS4'err(c.error)'E'
         end
   end
   if c.error ^= 0 | c.0 ^= 2 then call initialize
   c.error = 0;
   c.first = left(c.cmd,1);
   if find('DEG RAD TIME NOTIME',c.cmd) > 0
   then call setparms c.cmd 'NEW' c.say
   else if c.first = '*' then call out c.input2
   else if c.cmd = 'FIG' | c.cmd = 'CALC' then
        do;if datatype(word(c.input2,2)) = 'NUM'
                then call setparms c.cmd 'NEW' word(c.input2,2)
                else say 'ERROR:' c.cmd 'not followed by number'
        end
   else if c.cmd = 'FILE' then
        do; say 'Results will also go to file CALC RESULTS A'
            c.file = 'YES'
        end
   else if c.cmd = 'NOFILE' then
        do; say 'No more results to file CALC RESULTS A';
            c.file = 'NO'
        end
   else if c.input2 ^= '?' & c.input2 ^= 'HELP' then
        do; numeric digits c.sig
            parse value c.input2 with keywd '(' rest
            if c.tim = 'NOTIME' then call cputime 'NOTIME'
            rc = 0
            if find('VANG VSET vang vset',keywd) > 0
            then interpret c.input2; else
            do; if c.cmd = 'CMD' then
                do; parse var c.input2 cmd c.rest; interpret c.rest
                end;
                else
                do;interpret 'answer = ' c.input2 /* Does the sums */
                   numeric digits c.say
                   if datatype(answer) = 'NUM'
                   then do; c.num = answer * 1; call out c.num; end
                   else call out answer
            end;end
            if rc <> 0 then
            do; say 'ERROR in CMD command: ' c.input2
                say '  It gave return code ' rc
                signal error
            end
            rc = 0
            if c.tim = 'TIME' then call cputime 'TIME'
        end
   c.error = 0
   if c.input1 ^= '' & c.input1 ^= 'C.INPUT1' then call exitcalc
end /* of main loop starting at  do forever... */
call exitcalc
INITIALIZE: procedure expose c. e pi
   'GLOBALV INIT'
   'GLOBALV SELECT CALC GET C.TRIG C.SIGFIG C.SAYFIG C.TIM'
   c.sig = c.sigfig; c.say = c.sayfig; c.trg = c.trig;
   if datatype(c.sig,'W') = 0 then c.sig = 15
   if datatype(c.say,'W') = 0 then c.say = 10
   if c.trg = '' | c.trg = 'C.TRIG' then c.trg = 'RAD'
   c.maxh   = 2302585000  /* Slightly less than max sh or ch argument */
   call setparms c.trg 'RESET'
   return /* from initialize */
SHOWPARMS: procedure expose c.
   showmsg = 'CALC working to' c.sig 'digits (CALC' c.sig') displaying'
   call out showmsg c.say '(FIG' c.say') angles in' c.trg
   if c.file ^= 'YES'
   then filemsg = 'Results will go to terminal only (NOFILE).'
   else filemsg = 'Results to terminal and CALC RESULTS A (FILE).'
   if c.tim = 'TIME'
   then call out filemsg 'Commands will be timed (TIME).'
   else call out filemsg 'Timing info not given (NOTIME).'
return /* from showparms */
SETPARMS: procedure expose c. e pi ; arg cmd new figs rest
   cmdlist = 'DEG RAD FIG CALC TIME NOTIME'
   if find(cmdlist,cmd) = 0 then
    do;say 'Eh? Must specify one of' cmdlist 'here, but you gave' cmd;
       signal error;
    end
   else if cmd = 'TIME' | cmd = 'NOTIME' then c.tim = cmd
   else if cmd = 'FIG' | cmd = 'CALC' then
    do;if figs < 1 then
       do; say 'ERROR:' cmd figs 'impossible. Minimum 1' ; signal error
       end;
       if cmd = 'FIG' then c.say = figs; else c.sig = figs
       if c.say > c.sig then
       do;if cmd = 'FIG' then c.sig = figs; else c.say = figs
          imposs = "Note: CALC can't work (CALC n) to fewer figures"
          call out imposs "than it shows (FIG n)"
       end;
       numeric digits c.sig
    end;
   c.small = '1.0E-'c.sig
/*find pi c.pi180 e c.hl2p c.ln10 c.1 to c.10 to current accuracy */
   digits  = max(c.sig,52)
   numeric digits digits
   if digits > 53 & (c.sig > c.sigfig | c.0 ^= 2) then
    do; savetrig = c.trigscale; c.trigscale = 1;
        pi = 16*atanser(0.2) - 4*atanser(1/239)
        c.trigscale = savetrig
        e = exp(1)
        c.hl2p = 0.5*ln(2*pi); c.ln10 = ln(10);
    end; else
    do;pi     = 3.141592653589793238462643383279502884197169399375105821
       e      = 2.718281828459045235360287471352662497757247093699959575
       c.hl2p = .9189385332046727417803297364056176398613974736377834128
       c.ln10 = 2.302585092994045684017991454684364207601101488628772976
    end
   c.pi180  = pi/180
   if cmd = 'RAD' then
    do; c.trg = cmd; c.trigscale = 1;     c.triginv = 1; end;
   else if cmd = 'DEG' then
         do; c.trg = cmd; c.trigscale = c.pi180;c.triginv = 1/c.pi180;end
   if (c.sig > c.sigfig | c.0 ^= 2) then /* see LGAMMA */
   do;c.1 = 1/12;   c.2 = -1/360;      c.3 = 1/1260; c.4 = -1/1680;
      c.5 = 1/1188; c.6 = -691/360360; c.7 = 1/156;  c.8 = -3617/122400
      c.9 = 43867/244188;              c.10 = -174611/125400
   end
   call showparms
   c.0 = 2 /* this is to show parms have been set */
   if new = 'RESET' then return
   c.sayfig = c.say; c.sigfig = c.sig; c.trig = c.trg
   'GLOBALV SELECT CALC PUTP C.SAYFIG C.SIGFIG C.TRIG C.TIM'
   return /*from setparms*/
ERR:   procedure; arg code /* find the right REXX error message */
   elist = '0 0 51 52 50 53 54 55 56 57 58 59 60 61 62 63 65 91 82 '
   elist = elist ||  '20 64  0  0 84 85 66 67 86 87 68 69 92 88 70 '
   elist = elist ||  '71 72 73 89 74 75 76 77 78 79 80  0  0 90 81 '
return subword(elist,code,1)
MAGN:  procedure expose c.; arg x
/* For x = + or - x.yyyEzzz  stack powerof10 = 1Ezzz and exponent = zzz */
   w = format(abs(x),2,0,,0)
   epos = pos('E',w)
   if epos = 0 then do; powerof10 = 1; exponent = 0; end
               else do; powerof10 = 1 || delstr(w,1,epos-1)
                        exponent = delstr(w,1,epos);
                   end
   push powerof10 exponent
   return /*from magn*/
SQRT:  procedure expose c.; arg x;
   if x = 0 then return 0
   if x > 0 then
   do; call magn(x); pull . exponent
       t = '1E' || format(exponent/2,,0)
       tlast = -1.0
       do forever while t ^= tlast
          tlast = t
          t = 0.5*(tlast + x/tlast) /* Newton-Raphson */
       end
       return t /*from sqrt*/
   end;
   say "ERROR: SQRT argument ("||x||") negative or missing"; signal error
FAC:   procedure; arg x
   if (datatype(x,'W') = 0) | (x < 0) then
   do;say "ERROR: FACtorial argument ("||x||") must be +ve integer"
      signal error
   end
   if x > 1000 then
   do; warnstart = "Warning: this may take a long time."
       say warnstart "You can escape before the end with HX"
   end
   f = 1
   do n = 1 by 1 while n <= x
      f = f*n
   end
   return f /*from fac*/
LOG:   procedure; arg x
   say "ERROR: use LN or LOG10"; signal error
LN:    procedure expose c. e; arg x
   if x > 0 then
   do;c.sig = c.sig + 2
      numeric digits c.sig
      call magn(x); pull powerof10 exponent
      if datatype(c.ln10) ^= 'NUM' then c.ln10 = lnseri(10)
      ans = lnseri(x/(e*powerof10)) + exponent*c.ln10 + 1
      c.sig = c.sig-2
      numeric digits c.sig
      return ans /*from ln*/
   end
   say "ERROR: LN argument ("||x||") zero, negative or missing"
   signal error
LNSERI:procedure expose c.; arg x /* assumes x>0 */
   if x > 1.25 | x < 0.8 then return 2*lnseri(sqrt(x))
   y = x-1
   c.small = '1.0E-'c.sig
   sum = y; term = y
   do n = 2 by 1 while abs(term) > c.small
      term = -term*y; sum = sum+term/n
   end
   return sum /*from lnseri*/
ATANSER:procedure expose c. pi; arg x
   if     x  < 0 then return -atanser(-x)
   if abs(x) > 1 then return pi*0.5 - atanser(1/x)
   if abs(x) > 0.1 then return 2*atanser(x/(1+sqrt(1+x*x)))
   c.small = '1.0E-'c.sig
   term = x; sum = x; xx = x*x;
   do n = 3 by 2 while abs(term) > c.small
      term = -term*xx; sum = sum+term/n
   end
   return sum /*from atanser */
LOG10: procedure expose c. e; arg x
   if x > 0 then return ln(x)/c.ln10 /*from log10*/
   say "ERROR: LOG10 argument ("||x||") zero, negative or missing"
   signal error
POWER: procedure expose c. e; arg x,y
   return exp(y*ln(x))
EXP:   procedure expose c.; arg x
   sm = -230258058; big = 230258092.888
   if x < sm then
     do;say "ERROR: EXP argument ("||x||") missing or too small (min "sm")"
        signal error;
     end
   if x > big then
     do;say "ERROR: EXP argument ("||x||") too big (max "big")"
        signal error;
     end
   if x < 0 then return 1/exp(-x)
   if x > 1.0E8 then return exp(x*0.125)**8
   xint=trunc(x)
   extra = length(xint)+5
   c.sig = c.sig+extra
   numeric digits c.sig
   xfrac = x-xint
   if xfrac > 0.5 then xcalc = xfrac-0.5; else xcalc = xfrac
   roote = eseries(0.5)
   exint = roote**(xint*2)
   ans   = eseries(xcalc)*exint
   if xfrac > 0.5 then ans = ans*roote
   c.sig = c.sig-extra
   numeric digits c.sig
   c.small = '1.0E-'c.sig
   return ans
ESERIES:procedure expose c.; arg x /* x near 0 */
   c.small = '1.0E-'c.sig
   term = 1; y = 1;
   do n = 1 by 1 while abs(term) > c.small
      term = term*x/n; y = y+term
   end
   return y
SIN:   procedure expose pi c.; arg x
   return cos(pi*0.5*c.triginv - x)
COS:   procedure expose pi c.; arg x
   xabs = abs(x)
   xturns = xabs*c.trigscale/(2*pi)
   xrem   = (xturns//1)*2*pi      /* xrem in radians mod 2*pi*/
   cossign = +1;
   if xrem > pi then xrem = 2*pi-xrem;
   if xrem > pi*0.5 then do; cossign = -1; xrem = pi-xrem; end
   sum = 1; term = 1;
   do n = 2 by 2 while abs(term) > c.small
      term = -term*xrem**2/(n*(n-1)); sum = sum+term
   end;
   return sum*cossign
TAN:   procedure expose pi c.; arg x
   cosx = cos(x)
   if abs(cosx) > c.small then return sin(x)/cos(x); else
   do; say "ERROR: TAN("||x||") seems to be infinite"; signal error
   end;
ARCSIN:procedure expose pi c.; arg x;
   if abs(x) = 1 then return pi*0.5*c.triginv*sign(x)
   else if abs(x) < 1 then return atanser(x/sqrt(1-x**2))*c.triginv
   say "ERROR: ARCSIN argument ("||x||") too big (max 1)"; signal error
ARCCOS:procedure expose pi c.; arg x;
   if x = 1 then return 0
   else if x = -1 then return pi*c.triginv
   else if abs(x) < 1 then return pi*0.5*c.triginv-arcsin(x)
   say "ERROR: ARCCOS argument ("||x||") too big (max 1)"; signal error
ARCTAN:procedure expose c. pi; arg x; return atanser(x)*c.triginv
ASIN:  procedure expose c. pi; arg x
   if abs(x) <= 1 then return arcsin(x)
   say "ERROR: ASIN argument ("||x||") too big (max 1)"; signal error
ACOS:  procedure expose c. pi; arg x
   if abs(x) <= 1 then return arccos(x)
   say "ERROR: ACOS argument too big (max 1):" x; signal error
ATAN:  procedure expose c. pi; arg x
   return arctan(x)
VSET:  procedure expose c. pi x y z; arg lat,lon,size
   small = '1.0E-'c.say
   maxlat = pi * (0.5 + small) * c.triginv
   if abs(lat) > maxlat then
   do; numeric digits c.say
       say 'lat =' lat+0 'but max =' maxlat+0 c.trg
       numeric digits c.sig
       return ' '
   end
   x = size*cos(lat)*cos(lon);y = size*cos(lat)*sin(lon);z = size*sin(lat)
   numeric digits c.say
   call out 'X =' x+0 ' Y =' y+0 ' Z =' z+0
   numeric digits c.sig
   return ' '
VANG  :procedure expose c. pi lat lon long size; arg x,y,z
   if x = 0 then lon = pi*0.5*sign(y)*c.triginv;
            else lon = atan(y/x)
   if x < 0 then lon = lon + pi*sign(y)*c.triginv
   size = sqrt(x*x + y*y + z*z)
   long = lon
   lat  = asin(z/size)
   numeric digits c.say
   call out 'LAT=' lat+0 ' LON=' lon+0 ' SIZE=' size+0
   numeric digits c.sig
   return ' '
DIST:  procedure expose c. pi; arg lat1,lon1,lat2,lon2
   cosdist = cos(lat1)*cos(lat2)*cos(lon1-lon2) + sin(lat1)*sin(lat2)
   return acos(cosdist)
SH:    procedure expose c.; arg x;
   if abs(x) < c.maxh then return sinh(x)
   say "ERROR: SH argument ("||x||") too big (max "c.maxh")"; signal error
SINH:  procedure expose c.; arg x;
   if abs(x) < 1.0 then
   do; sum = x; term = x;
       do n = 3 by 2 while abs(term) > c.small
          term = term*x*x/(n*(n-1));  sum = sum+term
       end;
       return sum
   end; else
   if abs(x) < c.maxh then
   do expx = exp(x);
      return (expx-1/expx)*0.5
   end
   say "ERROR: SINH argument ("||x||") too big (max "c.maxh")"
   signal error
CH:    procedure expose c.; arg x;
   if abs(x) < c.maxh then return cosh(x)
   say "ERROR: CH argument ("||x||") too big (max "c.maxh")"; signal error
COSH:  procedure expose c.; arg x;
   if abs(x) < c.maxh then
   do; y = exp(abs(x)); return (y+1/y)*0.5
   end
   say "ERROR: COSH argument ("||x||") too big (max "c.maxh")"
   signal error
TANH:  procedure expose c.; arg x;if abs(x) < c.maxh
   then return sinh(x)/cosh(x); else return sign(x)
TH:    procedure expose c.;    arg x; return tanh(x)
ERF:   procedure expose c. pi; arg x; /* Taylor series if x <= 3 */
   if x < 0 then return -erf(-x)
   if x > 3 then return 1.0-erfc(x)
   term = 2*x/sqrt(pi); sum = term; xx = -x*x
   do n = 1 by 1 while abs(term*n) > c.small*sum
      term = term*xx/n; sum = sum+term/(2*n+1);
   end
   return sum
ERFCA: procedure expose pi c.;arg x /* asymp series */
   if x < 2 then return 1.0-erf(x)
   if abs(x) > 47900
   then do; say 'Max x for nonzero erfc(x) is 47900'; return 0; end
   exx = exp(-x*x)/(sqrt(pi)*x)
   xx  = 1/(2*x**2)
   sum = 1-xx
   term = -xx
   oldaterm = abs(term)
   do j = 3 by 2 while abs(term) > c.small*sum
      term = term*(-j*xx)
      newaterm = abs(term)
      if oldaterm < newaterm then
      do;say"ERROR: ERFC asymp. series fails for" c.sig "figures, x = "||x
         signal error
      end
      oldaterm = newaterm
      sum = sum+term
   end
   return exx*sum
ERFC  :procedure expose c. pi; arg x;  /* recurrence rel. i^{n}erfc */
   if x < 0 then return 2.0-erfc(-x) /* if asymp or Taylor useless */
   if x < 2 then return 1.0-erf(x)
   if x > max(6,1.52*sqrt(c.sig)) then return erfca(x)
             /* 1.52 > sqrt(ln(10)) */
   nmax = 10
   newans = 0
   do forever until abs(oldans-newans) < 100.0*c.small*newans
      nmax1 = nmax+1
      e.nmax1 = 0; e.nmax = 1
      do j = nmax-1 by -1 to -1
                      j1 = j+1; j2 = j+2
         e.j = 2*(x*e.j1 + j2*e.j2)
      end; minus1 = -1;
      oldans = newans
      newans = e.0*2/(sqrt(pi)*exp(x*x)*e.minus1)
      nmax = nmax*2
   end
   return newans
BETA:  procedure expose c. e; arg x,y
   return gamma(x)*gamma(y)/gamma(x+y)
GAMMA: procedure expose c. e; arg x
   if x > 1.382E-76 then
   do; if x <= 1.302E8 then
       do; elgama = exp(lgamma(x))
           return elgama
   end;end;
   say "ERROR: GAMMA argument ("||x||") out of range 1.382E-76 to 1.302E8"
   signal error
LGAMMA:procedure expose c. e; arg x
   maxl = 4.34294485E999999990;
   p = 1;
   y = x
   ratio = -c.10/c.small
   numeric digits 4
   c.small = 1.0E-4
   nmax = power(ratio,1/19);
   numeric digits c.sig
   c.small = '1.0E-'c.sig
   if x > c.small & x < maxl then
   do;if y < nmax then
      do; do n = 1 by 1 while y < nmax;
             p = p*y; y = y+1;
          end             /* p = x*(x+1)*(x+2)*...*(x+n) */
      end
      ypower = y; sum = 0;
      if c.small*y < 1 then y2 = 1/(y*y); else y2=0
      term = 1
      if y <= maxl then
      do; do n = 1 by 1 while abs(term) > c.small*y & n < 11;
             ypower = ypower*y2; term = ypower*c.n; sum = sum+term
          end
          return (y-0.5)*ln(y)-y-ln(p)+c.hl2p+sum
      end
   end
   say "ERROR: LGAMMA argument ("||x||") out of range "c.small" to "maxl
   signal error
/*end of lgamma*/
CPUTIME:procedure expose c.; arg show
   numeric digits 40
   newtime = format(c2d(substr(diag('C'),25,8))*1E-6,,2)
   numeric digits c.sig
   if datatype(c.time) = 'NUM'
   then do; t = newtime - c.time
            if show = 'TIME' then call out 'That took' t 'sec total CPU.'
   end
   c.time = newtime
   return /*from CPUTIME*/
OUT:   procedure expose c. pi e x y z lat lon long size
   parse arg stuff
   len = linesize()
/* must ask every time: can alter it with CMD CP TERMINAL LINESIZE nnn */
   if len = 0 then len = 255
   stuff1 = substr(stuff,1,len)
   stuff2 = substr(stuff,len+1)
   if c.first ^= '*' then say stuff
   if c.file = 'YES' then
   do; 'MAKEBUF'
       push stuff1
       'EXECIO 1 DISKW CALC RESULTS A'
       'DROPBUF'
   end
   if length(stuff2) > 0 then do; c.first = '*'; call out stuff2; end
   return
EXITCALC: procedure;
   'DROPBUF'
   exit
/* last line of CALC EXEC */