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