Discussion:
Happy π Day!
(too old to reply)
Paul Gilmartin
2017-03-16 02:30:59 UTC
Permalink
Raw Message
After a thread in the z390 forum:

/* Rexx */ signal on novalue; /*
# Nilakantha approximation to π.
After: https://helloacm.com/two-simple-equations-to-compute-pi/
See: https://en.wikipedia.org/wiki/Alternating_series#Alternating_series_test
*/
parse value arg( 1 ) 5 with ndig . /* Digits precision. */
numeric digits ndig * 2 /* Be generous with guard digits. */
EPSILON = 3 / 10 ** ndig
ans = 3
do j = 2 by 2 until term <= EPSILON
term = 4.0 / (j * (j + 1) * (j + 2))
ans = ans + ( j // 4 - 1 ) * term
say right( j % 2, 10 ) format( ans, 2, ndig + 2 ) format( term, 2, 6, 3, 0 )
end j

-- gil

----------------------------------------------------------------------
For TSO-REXX subscribe / signoff / archive access instructions,
send email to ***@VM.MARIST.EDU with the message: INFO TSO-REXX
Gerard Schildberger
2017-03-16 03:30:19 UTC
Permalink
Raw Message
Post by Paul Gilmartin
/* Rexx */ signal on novalue; /*
# Nilakantha approximation to π.
After: https://helloacm.com/two-simple-equations-to-compute-pi/
See: https://en.wikipedia.org/wiki/Alternating_series#Alternating_series_test
*/
parse value arg( 1 ) 5 with ndig . /* Digits precision. */
numeric digits ndig * 2 /* Be generous with guard digits. */
EPSILON = 3 / 10 ** ndig
ans = 3
do j = 2 by 2 until term <= EPSILON
term = 4.0 / (j * (j + 1) * (j + 2))
ans = ans + ( j // 4 - 1 ) * term
say right( j % 2, 10 ) format( ans, 2, ndig + 2 ) format( term, 2, 6, 3, 0 )
end j
-- gil
Here is another version to calculate pi using the AGM algorithm:


@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
/*REXX program calculates the value of pi using the AGM algorithm. */
parse arg d .; if d=='' | d=="," then d=500 /*D not specified? Then use default. */
numeric digits d+5 /*set the numeric decimal digits to D+5*/
z=1/4; a=1; g=sqrt(1/2) /*calculate some initial values. */
n=1
do j=1 until a==old; old=a /*keep calculating until no more noise.*/
x=(a+g)*.5; g=sqrt(a*g) /*calculate the next set of terms. */
z=z - n*(x-a)**2; n=n+n; a=x /*Z is used in the final calculation. */
end /*j*/ /* [↑] stop if A equals OLD. */

pi=a**2 / z /*compute the finished value of pi. */
numeric digits d /*set the numeric decimal digits to D.*/
say pi / 1 /*display the computed value of pi. */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2
do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/
numeric digits d; return g/1
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

It's pretty fast. BTW, the sqrt function if one of my own and written for large number of decimal digits (numeric digits).


Here is another version (of the above) which shows (more or less) a running
precision of the calculation:


$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
/*REXX program calculates value of pi using the AGM algorithm (with running digits).*/
parse arg d .; if d=='' | d=="," then d=500 /*D not specified? Then use default. */
numeric digits d+5 /*set the numeric decimal digits to D+5*/
z=1/4; a=1; g=sqrt(1/2) /*calculate some initial values. */
n=1
do j=1 until a==old; old=a /*keep calculating until no more noise.*/
x=(a+g)*.5; g=sqrt(a*g) /*calculate the next set of terms. */
z=z-n*(x-a)**2; n=n+n; a=x /*Z is used in the final calculation. */
many=compare(a,old) /*how many accurate digits computed? */
if many==0 then many=d /*adjust for the very last time. */
say right('iteration' j, 20) right(many, 9) "digits" /*show digits.*/
end /*j*/ /* [↑] stop if A equals OLD. */
say /*display a blank line for a separator.*/
pi=a**2 / z /*compute the finished value of pi. */
numeric digits d /*set the numeric decimal digits to D.*/
say pi / 1 /*display the computed value of pi. */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2
do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/
numeric digits d; return g/1
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

For 250 decimal digits, it (only) takes 9 iterates.
For 500 decimal digits, it (only) takes 10 iterates.
For 1000 decimal digits, it (only) takes 11 iterates.
For 2000 decimal digits, it (only) takes 12 iterates.
For 4000 decimal digits, it (only) takes 13 iterates.
For 8000 decimal digits, it (only) takes 14 iterates.

Needless to say, the larger the required decimal digits, the
longer the iterations take.



Here is another REXX program that "spits out" decimal digits
of pi as they are calculated (based on the number of
decimal digits specified; the number of digits used for
computation is specifiable. I suggest trying a "small"
number (on the command line) such as 1,000 or thereabouts.


############################################################
/*REXX program spits out decimal digits of pi (one digit at a time) until Ctrl-Break.*/
parse arg digs oFID . /*obtain optional argument from the CL.*/
if digs=='' | digs=="," then digs=1e6 /*Not specified? Then use the default.*/
if oFID=='' | oFID=="," then oFID='PI_SPIT.OUT' /* " " " " " " */
numeric digits digs /*with bigger digs, spitting is slower.*/
call time 'Reset' /*reset the wall─clock (elapsed) timer.*/
signal on halt /*───► HALT when Ctrl─Break is pressed.*/
pi=0; v=5; vv=v*v; g=239; gg=g*g; spit=0 /*assign some values to some variables.*/
s=16 /*calculate π with increasing accuracy */
r=4; do n=1 by 2 until old=pi; old=pi /*just calculate pi with odd integers*/
pi=pi + s/(n*v) - r/(n*g) /* ··· using John Machin's formula.*/
if pi==old then leave /*have we exceeded the DIGITS accuracy?*/
s=-s; r=-r; v=v*vv; g=g*gg /*compute some variables for shortcuts.*/
do j=spit+1 to compare(pi,old) /*spit out some (new) digits of π (pi)*/
parse var pi =(j) spit +1 /*equivalent to: spit=substr(pi,j,1) */
call charout ,spit /*display one (new) decimal digit of π.*/
call charout oFID,spit /*··· and also write π digit to a file.*/
end /*j*/ /* [↑] 0, 1, or 2 decimal dig are spit*/
spit=j-1 /*adjust for DO loop index increment.*/
end /*n*/
say /*stick a fork in it, we're all done. */
exit: say; say n%2+1 'iterations took' format(time("Elapsed"),,2) 'seconds.'; exit
halt: say; say 'PI_SPIT halted via use of Ctrl-Break.'; signal exit /*show iterations.*/
############################################################

It uses the charout BIF, so you may want to use a REXX such
as Regina on your PeeCee.
________________________________________ Gerard Schildberger

Loading...