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
