I knew I should have noted it when I saw it. I think it was a response on
a blog. I will see if I can dig it up. It was a few months ago.
On Wed, Dec 23, 2020 at 5:03 AM R.E. Boss <r.e.b...@outlook.com> wrote:
> " 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
>
--
Devon McCormick, CFA
Quantitative Consultant
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
