" I found it discouraging to read a recent missive from Niklaus Wirth where he
takes a jab at this useful simplification as an example of an unwarranted
complication. It's sad that an eminent luminary like this has failed to grasp,
after all these years, how useful this convention is."
What was that "recent missive from Niklaus Wirth" ?
R.E. Boss
-----Original Message-----
From: Programming <programming-boun...@forums.jsoftware.com> On Behalf Of Devon
McCormick
Sent: dinsdag 22 december 2020 19:44
To: J-programming forum <programm...@jsoftware.com>
Subject: Re: [Jprogramming] Chudnovsky algo for π (Pi) using extended arithmetic
Thanks everyone for the tips!
Bob - your suggestions still stalls:
20j18":%12x*-/chudSeriesa i.2
3.141592653589788641
20j18":%12x*-/chudSeriesa i.3
3.141592653589788641
(Should be
3.14159265358979323846...
as we all know :) )
Roger's suggestion will take me longer to translate.
I got into this diversion from thinking about the array languages FAQ I
mentioned in the KEI Centenary discussion last Thursday. One of the first
objections a lot of people have is the right-to-left order of evaluation and I
wanted a good illustration of how this gives you alternating sum with
-/ . I found it discouraging to read a recent missive from Niklaus Wirth
where he takes a jab at this useful simplification as an example of an
unwarranted complication. It's sad that an eminent luminary like this has
failed to grasp, after all these years, how useful this convention is.
The Chudnovsky formula is a good illustration because it includes a "_1^k"
term to force the alternating sum.
On Tue, Dec 22, 2020 at 12:56 PM Roger Hui <rogerhui.can...@gmail.com>
wrote:
> In the J source, file vx.c, function jtxpi(), you find an
> implementation of the Chudnovsky algorithm using extended precision
> arithmetic. It is in C but it should not be too hard to figure out the J
> equivalent.
>
> static XF1(jtxpi){A e;B p;I i,n,n1,sk;X a,b,c,d,*ev,k,f,m,q,s,s0,t;
> RZ(w); if(!XDIG(w))R xzero; ASSERT(jt->xmode!=XMEXACT,EVDOMAIN);
> RZ(a=xc(545140134L));
> RZ(b=XCUBE(xc(640320L)));
> RZ(c=xc(13591409L));
> RZ(d=xplus(xc(541681608L),xtymes(xc(10L),xc(600000000L))));
> n1=(13+AN(w)*XBASEN)/14; n=1+n1;
> RZ(e=piev(n,b)); ev=XAV(e); m=ev[n1]; f=xzero; s0=xone; sk=1;
> for(i=p=0;;++i,p=!p){
> s=xtymes(s0,xplus(c,xtymes(a,xc(i))));
> t=xdiv(xtymes(s,m),ev[i],XMEXACT);
> f=p?xminus(f,t):xplus(f,t);
> if(i==n1)break;
> DO(6, s0=xtymes(s0,xc(sk++));); RE(s0); /* s0 = ! 6*i */ }
> f=xtymes(d,f); q=xpow(xc(10L),xc(14*n1));
> k=xtymes(xtymes(a,m),xsqrt(xtymes(b,xsq(q))));
> R xdiv(xtymes(k,w),xtymes(f,q),jt->xmode);
> } /* Chudnovsky Bros. */
>
>
> On Tue, Dec 22, 2020 at 9:30 AM Devon McCormick <devon...@gmail.com>
> wrote:
>
> > The Chudnovsky algorithm -
> > https://en.wikipedia.org/wiki/Chudnovsky_algorithm - is supposed to
> > have the fastest convergence for pi (or 1/pi, to be exact). I had
> > tried this one which is OK but only seems to add one digit for each power
> > of 10:
> > 6!:2 'pi=. 4*-/%>:+:i. 1e6'
> > 0.0145579
> > 20j18":pi
> > 3.141591653589793420
> > 6!:2 'pi=. 4*-/%>:+:i. 1e7'
> > 0.140673
> > 20j18":pi
> > 3.141592553589793280
> > 6!:2 'pi=. 4*-/%>:+:i. 1e8'
> > 1.44487
> > 20j18":pi
> > 3.141592643589793177
> > 6!:2 'pi=. 4*-/%>:+:i. 1e9'
> > 16.7988
> > 20j18":pi
> > 3.141592652589793477
> >
> > The Chudnovsky series is fast and gives me 15 correct decimal digits
> > for only 2 iterations but fails to improve after that presumably
> > because of floating point limitations:
> > 20j18":%12*-/chudSeries i.2
> > 3.141592653589795336
> > 20j18":%12*-/chudSeries i.3
> > 3.141592653589795336
> > 20j18":%12*-/chudSeries i.4
> > 3.141592653589795336
> > chudSeries i.4
> > 0.0265258 4.98422e_16 2.59929e_30 1.45271e_44
> >
> > However, when I try to use extended precision, it does not seem to
> > give
> me
> > extended precision results:
> > chudSeries=: 13 : '((!6x*x: y)*13591409x+545140134x*x: y)%(!3x*x:
> > y)*(3x^~!x: y)*640320x^3r2+3x*x: y'
> >
> > I get the same fp values for the "i.4" argument as shown above. Am
> > I
> doing
> > something wrong or overlooking something?
> >
> > Thanks,
> >
> > Devon
> >
> > --
> >
> > Devon McCormick, CFA
> >
> > Quantitative Consultant
> > --------------------------------------------------------------------
> > -- For information about J forums see
> > http://www.jsoftware.com/forums.htm
> >
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
>
--
Devon McCormick, CFA
Quantitative Consultant
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
