learnpython.org - an online interactive Python tutorial

harrismh777 harrismh777 at charter.net
Sat Apr 23 20:37:40 EDT 2011


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 at 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




More information about the Python-list mailing list