> Exercise for the reader: What's a "use" of -\ ? (Published in June 1981.)
Sorry, -\ in APL, equivalent to -/\ in J for vectors. On Tue, Dec 22, 2020 at 11:48 AM Roger Hui <rogerhui.can...@gmail.com> wrote: > There are simpler demonstrations of the use of -/ . See > https://code.jsoftware.com/wiki/Essays/Extended_Precision_Functions#Sine_and_Cosine > > Exercise for the reader: What's a "use" of -\ ? (Published in June 1981.) > > On Tue, Dec 22, 2020 at 10:44 AM 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
