From: imc@comlab.ox.ac.uk (Ian Collier)
Newsgroups: comp.lang.rexx
Subject: Re: REXX math functions
Date: 13 Oct 92 15:16:45 GMT
Organization: Oxford University Computing Laboratory, UK
X-Local-Date: Tuesday, 13th October 1992 at 4:16pm BST
Originator: imc@msc3.comlab

In article <dej.01rz@qpoint.ocunix.on.ca>, dej@qpoint.ocunix.on.ca (David Jones) wrote:
>>     10 digit precision -- 6.76 seconds-- 3.141592656
>>     20 digit precision -- 40.15 seconds-- 3.14159265354941579892
>>     25 digit precision -- 63.0 seconds-- 3.141592653594157988802306

>Unfortunately, if these numbers are supposed to be approximations of pi,
>then they're off.

>I'll stick to my assembly routines (which agree with the value of pi
>generated by the Maple symbolic algebra package).

For the record, the above numbers seem to have been generated by the
algorithm which I posted earlier (copied below), and are probably due to the
use of a sqrt() routine which gives only about 10 digits of precision.  I
have said before and will say again that "my" algorithm is known to be
correct.  I didn't invent it - I heard it at a recreational maths lecture.
It has quadratic convergence: that is, the number of correct figures doubles
at each iteration.  It is therefore much better than any series summation.
I strongly recommend that Mr Jones tries it out.  I didn't post a sqrt
function last time.  I will this time, and thus all implementations of REXX
are more-or-less guaranteed to return the correct value of pi.

After 8 iterations, pi = 3.141592653589793238462643383279502884197169399
375105820974944592307816406286208998628034825342117069
16 seconds elapsed

After 10 iterations, pi = 3.14159265358979323846264338327950288419716939
937510582097494459230781640628620899862803482534211706798214808651328230
664709384460955058223172535940812848111745028410270193852110555964462294
895493038196442881097566593344612847564823378678316527120190914564856692
346034861045432664821339360726024914127372458700660631558817488152092096
2829254091715364367892590360011330530548820466521384146951941511610
214 seconds elapsed

    parse arg d       /* required number of digits */
    if d='' then d=20
    d=d+1
    numeric digits d
    numeric fuzz 1
    parse value 1 1/sqrt(2) 1/4 1 0 with a b c x p
    
    do j=1 until p=p1
       parse value (a+b)/2 sqrt(a*b) c-x*(b-a)*(b-a)/4 x+x with a b c x
       parse value p (a+b)*(a+b)/(4*c) with p1 p
    end
    
    numeric digits d-1
    say "After" j "iterations, pi =" format(p)
    exit

    sqrt:  /* You need this if and only if there is no builtin
              arbitrary-precision sqrt() function */
    procedure expose d
    parse arg a
    if a=0 then b=0
    else if a<0 then do
       say "sqrt("a"): Incorrect call to routine"
       return
    end
    else do
       numeric digits d+1
       numeric form scientific
       parse value format(a,1,,,0) with before '.' after 'E' e
       if pos('E',before)>0 then parse var before before 'E' e
       a=before||after  /* a is now pointless */
       if e='' then e=0 /* e is the exponent  */
       if e//2 then e=e-1
       else a='0'a      /* exponent is odd, and decimal comes after 2 digits */
       b='';b1=0;c=0
       numeric digits d%2+2
       do i=1 by 2 until length(b)=d+1
          b1=trunc(b1)||substr(a,i,2,0)
          /* find least x s.t. (c||x)*x>b1. Find (c||x-1)*(x-1). */
          d2=c||1
          d1=0
          x=1
          do while d2<=b1
             d1=d2
             d2=d2+(c||x)
             x=x+1
             d2=d2+x
          end
          x=x-1
          b=b||x
          b1=b1-d1
          c=trunc((c||x)+x) /* never in exponential form */
       end
       b=left(b,1)'.'substr(b,2)'E'e/2
    end
    numeric digits d
    return format(b)

Ian Collier
Ian.Collier@prg.ox.ac.uk | imc@ecs.ox.ac.uk
Newsgroups: comp.lang.rexx