Path: utzoo!utgpu!news-server.csri.toronto.edu!clyde.concordia.ca!uunet!munnari.oz.au!kaukau.comp.vuw.ac.nz!vuwst1!vuw4341!HARPER
From: HARPER@rs1.vuw.ac.nz (HARPER)
Newsgroups: comp.lang.rexx
Subject: Re: On-line calculator
Message-ID: <508@rs1.vuw.ac.nz>
Date: 28 Mar 90 03:34:32 GMT
References: <506@rs1.vuw.ac.nz>
Reply-To: HARPER@rs1.vuw.ac.nz
Organization: Victoria University of Wellington, NZ.
Lines: 508
Disclaimer: Author bears full responsibility for contents of this article

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 */