Chris Angelico wrote:
Wow, someone else who knows REXX and OS/2! REXX was the first bignum
language I met, and it was really cool after working in BASIC and
80x86 assembly to suddenly be able to work with arbitrary-precision
numbers!
Yes, my "big num" research stuff was initially done in REXX, on VM/CMS.
I later ported my libraries over to OS/2 and continued with that well
into the '90s, when I discovered Unix and 'bc'. Many folks are not
aware that 'bc' also has arbitrary precision floating point math and a
standard math library. It is much faster than REXX because the libraries
are written in C, and unlike REXX they do not have to be implemented in
the interpreter. The syntax of 'bc' is C-like, which is its only
down-side for new students who have never had any language training.
REXX was a great business language, particularly when CMS Pipelines was
introduced.
Just for fun, and a trip down memory lane, you can take a look at some
of my ancient REXX code... this code is a REXX math library designed to
be used from the command line, written in REXX. The primary scientific
functions were implemented, as well as some code to calculate PI...
most of the algorithms can be ported to 'bc' easily, but the 'bc'
algorithms will run much faster, of course.
Back in the day, REXX was the 'new BASIC'
... now I use Python.
/********************************************************************/
/********************************************************************/
/** COMPUTE COMMAND LINE SCIENTIFIC CALCULATOR **/
/********************************************************************/
/********************************************************************/
/** **/
/** IBM Corporation **/
/** HARRISMH at RCHVMP2 **/
/** EL8/663-1 B627 **/
/** Rochester, MN 55901 **/
/** **/
/** AUTHOR: Marcus Harris **/
/** DATE: 93/10/22 ( port from VM/CMS ) **/
/** UPDATE: 98/03/07 **/
/** VER: 2.0a **/
/** **/
/********************************************************************/
/********************************************************************/
/** syntax **/
/** **/
/** COMPUTE expression<;expression><;expression><@digits> **/
/** **/
/** **/
/** The expression(s) will be computed and displayed in terminal **/
/** line mode to arbitrary <@digits> of accuracy. **/
/** **/
/** The expression(s) may be any REXX math clause and may include **/
/** a variable assignment, for use in a subsequent expression: **/
/** ie., COMPUTE A=3;B=7;A/B;A*B@30 **/
/** (will compute both the quotient and product of A & B) **/
/** **/
/********************************************************************/
/********************************************************************/
/** NOTES **/
/** **/
/** COMPUTE expression<;expression><;expression><@digits> **/
/** **/
/** This program exploits Rexx INTERPRET and NUMERIC DIGITS **/
/** **/
/** This exec is portable to OS/2, NT and w95 Rexx/objectRexx **/
/** without changes. **/
/** **/
/********************************************************************/
/********************************************************************/
/** updates **/
/** 93/10/21 mhh New Program (port from VM/CMS) **/
/** 98/03/07 mhh New Version 2.0a **/
/** This version includes trigonometric, logarithmic, **/
/** exponential and combinatorial functions. **/
/** **/
/** See the REXROOTS ABOUT file for a discussion **/
/** of the included functions and some examples, **/
/** using the COMPUTE program. **/
/** **/
/********************************************************************/
/********************************************************************/
/** dependencies **/
/** **/
/** Numeric Fuzz: **/
/** Most of these routines bump-up the numeric digits **/
/** by a fuzzFactor. Rexx fuzz is used to terminate **/
/** the do-forever iteration at the correct number of **/
/** required places. **/
/** The routine then sets the numeric digits back to **/
/** the application setting and returns result-0. **/
/** The (result-0) forces truncation/rounding to give **/
/** the expected result. **/
/** Some implementations of REXX do not support the **/
/** fuzz option; the iterative do-forever methods **/
/** would need modification for correct termination. **/
/** **/
/********************************************************************/
/********************************************************************/
/** mainline **/
/********************************************************************/
arg expression '@' max_digits .
select
when expression='' | expression='?' then say syntax_diagram()
when datatype(max_digits,'W') & max_digits>0 then,
numeric digits max_digits
otherwise
numeric digits 18
end
call evaluate expression
exit 0
/********************************************************************/
/** procedures **/
/********************************************************************/
SYNTAX_DIAGRAM: procedure
return'Syntax: COMPUTE expression <;expression><@precision>'
EVALUATE: procedure expose max_digits
arg expression
signal on error; signal on syntax
/* 502 digits pi/2 */
hpi="1.57079632679489661923132169163975144209858469968755291048747229615390820314310449931401741267105853"
hpi=hpi"3991074043256641153323546922304775291115862679704064240558725142051350969260552779822311474477465190"
hpi=hpi"9822144054878329667230642378241168933915826356009545728242834617301743052271633241066968036301245706"
hpi=hpi"3686229350330315779408744076046048141462704585768218394629518000566526527441023326069207347597075580"
hpi=hpi"471652863518287979597654609305869096630589655255927403723118998137478367594287636244561396909150597456"
/* 502 digits pi */
pi="3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679"
pi=pi"8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"
pi=pi"4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273"
pi=pi"7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094"
pi=pi"3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912"
resetTime=time('R')
do while expression <> ''
parse value expression with exp1';'expression
select
when pos('=',exp1)<>0 then do
interpret 'do; 'exp1'; end;'
end
otherwise
interpret 'do; answer='exp1'; end;'
select
when answer='' then say syntax_diagram()
when datatype(answer,'N') then say answer
otherwise signal error
end
end
end
say "elapsed time:" time('E')
return
/********************************************************************/
/** routines **/
/********************************************************************/
ERROR:; SYNTAX:
say 'Does Not Compute'
exit 28
/********************************************************************/
/** functions **/
/********************************************************************/
fx: procedure /* template
*/
arg x
maxDigits=digits()
numeric digits maxDigits*2
numeric fuzz maxDigits
/* f(x) goes here */
numeric fuzz 0
numeric digits maxDigits
return x-0
/********************************************************************/
/** Exponential - Logarithmic **/
/********************************************************************/
exp: procedure /* (a>=0) (n>=0)
*/
arg a,n
select
when datatype(n, "W")=1 then x=a**n
otherwise x=e(ln(a)*n)
end
return x
ln: procedure /* series methods ( x > 0 ) */
arg x
if x<=0 then return "ERROR"
maxDigits=digits()
numeric digits maxDigits+7
numeric fuzz 7
select
when x>80 then do
i=0
do while x>=5
i=i+1
x=x/10
end
x=ln(x)+i*ln(10)
end
when x<.5 then do
i=0
do while x<.5
i=i+1
x=x*10
end
x=ln(x)-i*ln(10)
end
when x>1.6 then do
u=(x-1)/(x+1); x=u; c=1; s=u**2
do forever
c=c+1
u=s*u
x1=x+u/(2*c-1)
if x1=x then leave
x=x1
end
x=x*2
end
otherwise
u=(x-1); x=u; c=1; s=u; sign=1
do forever
c=c+1
u=u*s
sign=sign*(-1)
x1=x+sign*u/c
if x1=x then leave
x=x1
end
end
numeric fuzz 0
numeric digits maxDigits
return x-0
e: procedure /* (all real x)
e^x = 1+x+x^2/2!+x^3/3!-x^4/4!... */
arg x
maxDigits=digits()
select
when x<(-5) then do
fuzzFactor=abs(x)%5
numeric digits maxDigits*fuzzFactor
numeric fuzz maxDigits*(fuzzFactor-1)
end
otherwise
numeric digits maxDigits+7
numeric fuzz 7
end
n=x; s=x; x=x+1; d=1; c=1
do forever
c=c+1
n=n*s
d=d*c
x1=x+n/d
if x1=x then leave
x=x1
end
numeric fuzz 0
numeric digits maxDigits
return x-0
/********************************************************************/
/** Trigonometric **/
/********************************************************************/
tan: procedure /* (all real x)
*/
arg x,hyper
maxDigits=digits()
numeric digits maxDigits*2
numeric fuzz maxDigits
if hyper="" then x=sin(x)/cos(x)
else x=sinh(x)/cosh(x)
numeric fuzz 0
numeric digits maxDigits
return x-0
tanh: procedure /* (all real x)
*/
arg x
return tan(x,"h")
cos: procedure /* (all real x)
cos(x)=1-x^2/2!+x^4/4!-x^6/6!...
cosh(x)=1+x^2/2!+x^4/4!+x^6/6!... */
arg x,hyper
maxDigits=digits()
numeric digits maxDigits*2
numeric fuzz maxDigits
n=1; s=x**2; x=1; sign=1; d=1; c=1
select
when hyper="" then do forever /* cosine function */
c=c+1
n=n*s
sign=(-1)*sign
do i=2 to 3 by 1
d=d*(2*c-i)
end
x1=x+sign*n/d
if x1=x then leave
x=x1
end
otherwise
do forever /* hyperbolic cosine function */
c=c+1
n=n*s
do i=2 to 3 by 1
d=d*(2*c-i)
end
x1=x+n/d
if x1=x then leave
x=x1
end
end
numeric fuzz 0
numeric digits maxDigits
return x-0
cosh: procedure /* (all real x)
*/
arg x
return cos(x,"h")
sin: procedure /* (all real x)
sin(x)=x/1!-x^3/3!+x^5/5!-x^7/7!...
sinh(x)=x/1!+x^3/3!+x^5/5!+x^7/7!... */
arg x,hyper
maxDigits=digits()
numeric digits maxDigits*2
numeric fuzz maxDigits
n=x; s=x**2; sign=1; d=1; c=1
select
when hyper="" then do forever /* sine function */
c=c+1
n=n*s
sign=(-1)*sign
do i=1 to 2 by 1
d=d*(2*c-i)
end
x1=x+sign*n/d
if X1=x then leave
x=x1
end
otherwise
do forever /* hyperbolic sine Function */
c=c+1
n=n*s
do i=1 to 2 by 1
d=d*(2*c-i)
end
x1=x+n/d
if X1=x then leave
x=x1
end
end /*select*/
numeric fuzz 0
numeric digits maxDigits
return x-0
sinh: procedure /* (all real x)
*/
arg x
return sin(x,"h")
sin_1: procedure expose hpi /* ( x^2<1 -pi/2 < sin_1(x) < pi/2 )
sin_1(x)=x+(x^3/2*3)+(x^5*1*3/2*4*5)+(x^7*1*3*5/2*4*6*7)... */
arg x
if x**2 >1 then return "ERROR"
maxDigits=digits()
numeric digits maxDigits*2
numeric fuzz maxDigits
n=x; d=1; s=x**2; c=1; p=0
if abs(x)>(.7854) then do
p=1; if x<0 then p=-1
s=1-s; x=sqrt(s); n=x
end
do forever
c=c+1; cc=2*c
n=n*(cc-3)*s
d=d*(cc-2)
x1=x+n/(d*(cc-1))
if x1=x then leave
x=x1
end
numeric fuzz 0
numeric digits maxDigits
select
when p>0 then return hpi-x
when p<0 then return x-hpi
otherwise nop
end
return x-0
cos_1: procedure expose hpi /* ( x^2<1 0 < cos_1(x) < pi )
*/
arg x
return hpi-sin_1(x)
tan_1: procedure expose hpi /* ( all real x )
*/
arg x
return hpi-cos_1(x/sqrt(x**2+1))
/********************************************************************/
/** Combinatorial **/
/********************************************************************/
fact: procedure /* n! n=(0,1,2,3...) */
arg n
numeric digits 500
m=1
do i=2 to n by 1
m=i*m
end
return m
part: procedure /* Ordered Partition */
arg m parts
numeric digits 500
mFact=fact(m)
do while parts<>''
parse var parts part parts
if datatype(part,'W') & m-part>=0 then do
m=m-part
mFact=mFact/fact(part)
end
end
return mFact/fact(m)
/********************************************************************/
/** PI routines **/
/********************************************************************/
PI: procedure /* AGM Gauss-Legendre Method
Crandall: Projects In Scientific Computation */
arg n
maxDigits=digits()
numeric digits maxDigits+7
numeric fuzz 7
x=sqrt(2); p=2+x; y=sqrt(x)
do forever
s=sqrt(x)
x=(s+1/s)/2
p1=p*((x+1)/(y+1))
if p1=p then leave
p=p1
s=sqrt(x)
y=((y*s)+(1/s))/(y+1)
end
numeric fuzz 0
numeric digits maxDigits
return p-0
PIa: procedure /* AGM Gauss-Legendre Method */
arg n
maxDigits=digits()
numeric digits maxDigits+7
numeric fuzz 7
a=1; b=1/sqrt(2); t=1/4; x=1
do i=1 to 100
y=a
a=(a+b)/2
b=sqrt(b*y)
if a=b then leave
t=t-x*(y-a)**2
x=2*x
end
PI=((a+b)**2)/(4*t)
numeric fuzz 0
numeric digits maxDigits
return PI-0
/********************************************************************/
/** Roots **/
/********************************************************************/
sqrt: procedure /* (all real x) */
arg x
return nrt(x,2)
cbrt: procedure /* (all real x) */
arg x
return nrt(x,3)
nrt: procedure /* ( a>=0 ) ( n>=0 )
*/
arg a,n
maxDigits=digits()
numeric digits maxDigits+7
numeric fuzz 7
if a=1 | a=0 then return a
if n=0 then return 1
if a<0 | n<0 then return "ERROR"
select
when datatype(n, "W" )=1 then do
m=n-1
x=bisect(a,n)
do forever /* Newton Method x:= x - ( f(x)/f'(x) ) */
x1=(m*x**n+a)/(n*x**m)
if x1=x then leave
x=x1
end
end
otherwise
x=e(ln(a)/n)
end
numeric fuzz 0
numeric digits maxDigits
return x-0
bisect: procedure /* Bisect Method used to seed nrt() */
arg n,root /* (n>=0) (root: 2,3,4,5,6...) */
numeric digits 12
numeric fuzz 3
a=0; b=n
if n<1 then b=1
else a=1
do forever
m=(a+b)/2
Fm=m**root-n
select
when Fm>0 then b=m
when Fm<0 then a=m
otherwise
a=b
end
if a=b then leave
end
return b
--
http://mail.python.org/mailman/listinfo/python-list
