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
