You're taking a square root there, which cannot be represented rationally.
If you want this to work, you'll need to decide on how much precision
you need, and use an approximation which is adequate for that
precision.
You should probably also introduce some words to represent
intermediate results, to avoid email line wrap issues.
(I would also prefer using something along the lines of chudSeries&x:
rather than repeatedly using x: within the definition.)
Good luck,
--
Raul
On Tue, Dec 22, 2020 at 1:44 PM Devon McCormick <devon...@gmail.com> wrote:
>
> 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
