 fast accurate Kneser sexp algorithm - Printable Version +- Tetration Forum (https://math.eretrandre.org/tetrationforum) +-- Forum: Tetration and Related Topics (https://math.eretrandre.org/tetrationforum/forumdisplay.php?fid=1) +--- Forum: Computation (https://math.eretrandre.org/tetrationforum/forumdisplay.php?fid=8) +--- Thread: fast accurate Kneser sexp algorithm (/showthread.php?tid=486) Pages: 1 2 3 4 5 fast accurate Kneser sexp algorithm - sheldonison - 08/07/2010 Here is the most recent code, Nov 12th, 2013: kneser.gp (Size: 35.84 KB / Downloads: 968) updates Nov 12th 2013 o updated to be compatible with latest pari-gp o removed eval, and kill o enhanced prtpoly routine o added lambertL routine end updates Nov 12th, 2013 This is the first of several posts. update: changed the names of the routines to match current code. I found a fast accurate way to implement the Knesser sexp, and implemented it in Pari-GP, which is a free mathematics program available on the web. The algorithm generates the Knesser sexp accurate to 15 decimal digits in about a minute. This includes a 400 term Riemann mapping approximation, and a 50 term Cauchy unit circle approximation to generate the Taylor series. I've tried the algorithm for base 2, base e, and base 10, and they are all well behaved. I will post the pari-GP code, as well as the Cauchy series it generated. The algorithm is for bases > eta. The algorithm is to start with the linear approximation for the sexp, that has a smooth first derivative. For base e, this is the approximation that the sexp(z) from [-1 to 0]=z+1. This region is extended with the exponent/logarithm to extend the initial unit length estimate. This is the initial trivial Taylor series. In my program, this Taylor series approximation for the sexp is called sexp. Then iterate as follows: 1) generate the Riemann mapping of the unit length segment This involves the Laurent series of the inverse superfunction of the unit length segment, mapped to the unit circle. The complex valued superfunction and its inverse are developed from the complex fixed point. Specifically, over a unit length, generate inverse_superfunction(sexp(z))-z, which is then mapped to the unit circle. Use the Cauchy Integral formula to generate the Taylor series. Strictly speaking, this function has a Laurent series, not a Taylor series, but we ignore the x^-n terms in the Laurent series, and only keep the x^+n terms, where n>=0. This is then trivially mapped back into a complex valued 1-cyclic fourier series, theta(z). 2) riemaprx(z). This 1-cyclic fourier series is used to generate a sexp approximation defined only in the imag(z)>0 portion of the plane. riemaprx(z)=superfunction(z+theta(z)). riemaprx(z) is poorly behaved at the real axis, with high frequency oscillations due to the slowly converging theta function. But the high frequency terms decay very rapidly as imag(z) increases in the complex plane, and if imag(z)>=1*I, the function behaves close to the complex superfunction. 3) sexp(z). Using riemaprx(z), generate the upper half of a unit half circle, conventionally centered on z=0, in the upper half of the complex plane. Treat the lower half of the unit circle as the complex conjugate of the upper half of the unit circle. Use the Cauchy integral to generate a Taylor series. This is a real valued Taylor series. The Taylor series is renormalized and recentered, so that sexp(0)=1, and so that sexp(+0.5)=exponent(sexp(-0.5). Iterate steps 1) 2) and 3). The claim is that this infinite iteration converges to Knessder's sexp(z), which is real valued at the real axis, with sexp(0)=1, and sexp(z+1)=exp(sexp(z)). Seven iterations were required to get 15 digits of accuracy, which each iteration adding approximately 8 binary bits of accuracy. In the Pari-GP program, each iteration takes less than 10 seconds. Implementing this algorithm is highly problematic and requires very large numbers of terms in the 1-cyclic theta function, due to the very slow convergence of the Riemann mapping. I found a trick to dramatically improve the convergence of the Riemann mapping, to handle the singularity, and to allow fast convergence. Here a few more details, before I call it a night. In my implementation, I wanted a 50 term Taylor series for the sexp approximation. This requires 200 points around a unit circle. Since we only generate the upper half of the unit circle, that's 100 points. These points are evenly spaced, with all points having imag(z)>0. The first point has z=e^(I*pi/200), which has imag(z)=0.015. This means that in riemprx(z), the imag(z) will always be >=0.015. That in turns means that the high frequency oscillations decay. For Riemann mapping theta terms>~300 is not computationally significant for any of these 100 points on the unit circle. This is the first of several posts. The next post will deal with the ideas required to improve the Riemann mapping convergence problems due to the singularity. I need to get some sleep, and it may take a few days to post all of my results. I want to clean up the code a little bit, since there's a lot of junk in the code that I needed to debug it and get it all working that I want to remove before posting, to avoid unnecessary confusion. Teaser: This is the output of my program, to generate the Taylor series terms I generated for base e, base 2, and base 10. - Sheldon base 2.71828182845904523536028747135266249775724709369 Cauchy Taylor series approximation terms 0 0.999999999999999920606916968922030701567136645654 1 1.09176735125832094016082238370083804714 2 0.271483212901694628277009772276018058438 3 0.212453248176256341682340133816642561735 4 0.069540376139987455601815633290052837540 5 0.0442919520904733787858836715649205635588 6 0.0147367420963894641781675687873762988951 7 0.0086687818172253171998963830218805136384 8 0.00279647939838551204530926929318637318059 9 0.00161063129058431027448748226836099353371 10 0.00048992723148441335050668409792945192050 11 0.000288181071154071176125207353574308091920 12 0.000080094612538568642532817921403719623613 13 0.0000502911417938227014378026179786713377326 14 0.0000121837903449185625957879474875634578856 15 0.0000086655336673937799848419995967770476324 16 0.00000168778231933152083906008897805134285921 17 0.00000149325324858236409519824063327981138611 18 0.000000198760764215837789507617671297185393027 19 0.000000260867356010921328462837051102110834560 20 0.0000000147099541512649833075949006297264385980 21 0.0000000468344973323925542131416564470746086702 22 -0.00000000154924165843564586518006754343312538792 23 0.0000000087415107851444402981222259166563437222 24 -0.00000000112578730421439340057613203197171894664 25 0.00000000170795927016981912380839845219103887952 26 -0.000000000377858310553104789819392071420497066106 27 0.000000000349577878719151948745500164259927239404 28 -1.05377008157789121410377229095750742668 E-10 29 7.4590973139432053510241775086048122696 E-11 30 -2.71759784858227219934013489755749519135 E-11 31 1.64607673327765883270714032926236921945 E-11 32 -6.7418700690910586159590314585805404498 E-12 33 3.7253295931126574070795526274392527304 E-12 34 -1.63908465536567938669004255576025126318 E-12 35 8.5836440640244217926800645809445203856 E-13 36 -3.94371548069373997097518693401613112984 E-13 37 2.00252641630788814440567087736426238381 E-13 38 -9.4417589560533107254019594050715409934 E-14 39 4.71206677912917615338890731001808112440 E-14 40 -2.25611366289088681618159333858518060425 E-14 41 1.11546306823408543027263377496376191447 E-14 42 -5.3891796465832969287083377311762631343 E-15 43 2.65194707551149745587030410085037067028 E-15 44 -1.28753286810996291447909130936705727867 E-15 45 6.3232296660338767222943776801526973230 E-16 46 -3.07332373088892725649620921357557234127 E-16 47 1.50855996920297588144091655941600050227 E-16 48 -7.2896911298653544932621055870671863118 E-17 49 3.57053526604856938153249432482269325976 E-17 base 2 Cauchy Taylor series approximation terms 0 1.00000000000000007750805374695916476506324315386 1 0.88936495462097636471560015559202608188 2 0.0086765489653699778435866428210914580477 3 0.095238800075181821015052159076811795525 4 -0.0057523485401261058429044155126967586843 5 0.0129665820200371929900639963746185426267 6 -0.00219604962303099108685649788446966214674 7 0.00199674684791144956471972891306784395245 8 -0.00056335481487852293422791227955443116661 9 0.000348242328188164833713322066658348490588 10 -0.000128532441264722069102630177525101178612 11 0.000067081924420528746629841393628806788277 12 -0.0000282987528228002605436935419049427621654 13 0.0000138001319906297777277157382455652644180 14 -0.0000062019093983767935657716936885382308667 15 0.00000295556146480613946456248019931704026334 16 -0.00000136867922453685047021475603425081919190 17 0.00000064905707564831820958280933643834369294 18 -0.000000305166939330940495799879426965798767859 19 0.000000144948206147741006787309611484943294178 20 -0.000000068746643115668186072914364608901276865 21 0.0000000327674451745520233209052672738147657033 22 -0.0000000156310874697440742179774639803902393806 23 0.0000000074781283860060906963656804741326980492 24 -0.00000000358267681486784107606081902426567798118 25 0.00000000171986774278787862816161264617020453732 26 -0.00000000082681596398826102575453499011874838983 27 0.000000000398110465534940913092926440510368851600 28 -1.91942994009799320930467536291965520002 E-10 29 9.2663184086753526923394526399596161631 E-11 30 -4.47869896241156433452311745776051738114 E-11 31 2.16712006481128547956338987778919331607 E-11 32 -1.04969755427424430603006383413954325851 E-11 33 5.0894457415080947575828909582705030799 E-12 34 -2.46987949612974961256130795951078775734 E-12 35 1.19965332790083051126762500373890740311 E-12 36 -5.8316699693946087927290168849662635343 E-13 37 2.83700134396762172106948455139976794164 E-13 38 -1.38119296763538037439839695675052294198 E-13 39 6.7286246573391428225154736432311797155 E-14 40 -3.28040714229361317312282338553709457619 E-14 41 1.59994527980264052587204787980778074081 E-14 42 -7.8111883641343962658543820790964098310 E-15 43 3.8123345996113862398910808416038811110 E-15 44 -1.86468916410329958272266244696629139897 E-15 45 9.0928731462620461748731226835758948880 E-16 46 -4.4652219382963819619261489249575291192 E-16 47 2.16258309701983021090373383869474626354 E-16 48 -1.07562115928439172115561558722739104721 E-16 49 5.05103639405910509662889857211651054178 E-17 base 10 Cauchy Taylor series approximation terms 0 1.00000000003413666821707397610280580973 1 1.6842890280105318114702665361 2 1.4130640093559017784867641066 3 1.3295945895120850984491868131 4 1.1104423522130456306141162585 5 0.90791475960484107525431757972 6 0.70590812180827017012794979492 7 0.53596957102484588479532491192 8 0.39514882149633253685483274509 9 0.28568024972253825344818666302 10 0.20237131278470927324534438107 11 0.14110168269002467755684498200 12 0.096870196895171185635561049276 13 0.065642683799764598157183797364 14 0.043932559903871435514499480872 15 0.029083045516276459146860315440 16 0.019054947684195538710084470697 17 0.0123686846345218426147071658882 18 0.0079582306458762837173385697464 19 0.0050792074578440317507369752633 20 0.0032170388379653970425789785613 21 0.0020231559289094629314623265337 22 0.0012638006922446651516872615141 23 0.00078448455296479211241357653392 24 0.00048404347686816718054617660014 25 0.00029697714889618699433841203816 26 0.00018122446592630282867291500326 27 0.000110023183879937876785387661614 28 0.000066469809662926642710985803857 29 0.000039970210188298586919385209485 30 0.000023927948827218540851314911165 31 0.0000142631537294666506633793631158 32 0.0000084672146369634588559435171282 33 0.0000050067182976642942211046734889 34 0.0000029492933862829666416493809347 35 0.0000017310105129704227353019917697 36 0.0000010124091578563321764991987405 37 0.00000059013055820603257123183162307 38 0.00000034287053679750024168609816182 39 0.00000019859459586612750224033089360 40 0.000000114688324996137110089024467356 41 0.000000066049601906226671780697958378 42 0.000000037941008758273721620067345461 43 0.000000021745359672035340985615756584 44 0.0000000124403916641799613324962201452 45 0.0000000071074328972264520946865200883 46 0.0000000040598591783922772751609136103 47 0.0000000023194687970727546265839419100 48 0.0000000013297332794947389880115713152 49 0.00000000076366914316723959202509556617 The pari-GP code - sheldonison - 08/07/2010 Next, I'd like to post the pari-GP code I have right now. Then I'll translate each of the subroutines in the code back into mathematics language, talk about convergence, and how I handled the singularity in the Riemann mapping, without requiring an infinite number of terms in the very slowly converging Riemann mapping. Lots of posts to come. Download this code, at the end of this post, and save it to a file. There is also an attachment at the end of this post with the "kneser.gp" file. Save it in the pari-GP directory, and start pari-GP. Type "\r kneser.gp". The first thing that comes up is a help menu, with the various routines I wrote, and a short one line description of what they do. The program starts initialized for base=e. By the way, here is the link for pari-GP download. Here's what comes up in the pari-GP window: Code:loop(n)      n iteration loops, each loop adds ~8 bits precision init(base)   initializes program for base, default loop(exp(1))=e" base_e_2_10  initializes and loops and generates graphs/data for base e, 2, 10 sexp(z)      sexp(z), returns Taylor series approximation of sexp(z) riemaprx(z)  riemen approximation of sexp(z), imag(z)>=idelta, 0.016*I superf(z)    superfunction of z, base B isuperf(z)   inverse superfunction of z, base B plotary      plots and dumps sexp/Riemann Taylor series arrays to file help         print this menu loop(6)      suggestion n=6 iteration loops Follow the suggestion and type "loop(6)", to loop through the estimation iteration routine six times. This will take about a minute. The program has now calculated a 50 term Taylor series for the sexp approximation for base=e accurate to 15 significant digits. To access that routine, type "sexp(0)"=1.00000000000000. "sexp(1)"=2.71828182845905. "sexp(-0.5)"=0.498563287941115. This is the half iterate for base=e. The code handles complex arguments just fine. "sexp(-2.5)"=-0.362370072029386 + 3.14159265358979*I. "sexp(0.5*I)"=0.936255672826876 + 0.520646429188355*I. The sexp(z) routine is accurate for -I0.015. For example, type riemaprx(0+4*I)=0.331375339574123 + 1.34164499758753*I. Type "help" at any time to reprint the help menu. Try typing "init(2)". This initializes the program for base=2. Once again "loop(6)" to generate another taylor series. sexp(-0.5)=0.543035666027967. This is the half-iterate for base 2. Cool! To go back to base=e, type "init(exp(1));loop(6)" Two other really useful routines are "superf(z)", and "isuperf(z)". These implement the complex valued superfunction and the complex valued inverse super-function for the current base. "superf(0)"= 1.14934109947623 + 1.12294996859363*I. "isuperf(0.5)" = -0.989970104598650 - 0.856059475588971*I. And here is the code itself (also, see downloadable attachment after the code). For the most recent code version: go to the Nov 21st, 2011 thread. - Sheldon Code:\p 38;             /* default precision 38 decimal digits */ ctrm      =    50; /* terms in the sexp Taylor series approximation    */ cvec      =   100; /* elements to calculate sexp Taylor series, 2*ctrm */ rtrm      =   400; /* terms in the Riemann mapping Taylor series       */ rvec      =   800; /* elements to calculate the Reimann Taylor series, 4*rtrm */ sfuncprec = 1E-17; /* precision for L fixed point approximation */ /* Default values saved for future reference, 56 binary bits, 17 decimal digits */ /* \p 38;             */                                                           /* ctrm      =    50; */                                                           /* cvec      =   100; */                                                           /* rtrm      =   400; */                                                           /* rvec      =   800; */                                                           /* sfuncprec = 1E-17; */  /*  57 binary bits, 17 decimal digits, loop(7) */         /* higher precision 75 term taylor series, 106 binary bits, 32 decimal digits */   /* \p 67;             */  /* loop(14) */                                           /* ctrm      =    75; */                                                           /* cvec      =   150; */                                                           /* rtrm      =  1000; */                                                           /* rvec      =  2000; */                                                           /* sfuncprec = 1E-32; */                                                           /* ultra high precision values, saved for future reference, 213 binary bits */ /* \p 134;            */  /* loop(28), ~6-7minutes/loop */                         /* ctrm      =   110; */                                                           /* cvec      =   220; */                                                           /* rtrm      =  3000; */                                                           /* rvec      =  6000; */                                                           /* sfuncprec = 1E-64; */                                                           /* L, B, scnt, ctrmat, rtrmat, idelta are global variables set by init(base), rtrmat set by riemup */ centerat  =     0; /* for larger bases, try centerat=-0.5 ...   */ renormup  =     1; /* used in the recenterup routine, for experiments, try renormup=0 */ curprecision=0; lastprecision=-2; B         = exp(1); L         = 0.3181315052047641353126542516 + 1.337235701430689408901162143*I; scnt      =    90; /* updated by init based on sfuncprec        */ ctrmat    =    50; /* chaos/overflow prevention */ rtrmat    =   400; /* chaos/overflow prevention */ idelta    =     0; /* I*imag(ccrc) set by init(b) initialization routine   */ imult     =     1; /* also set by init(b), corresponds to e^2pi*idelta */ recenter  =     0; /* updated by recenterup called by init(b) */ renorm    =     1; /* set by init(b) */ rer       = vector (rvec,i,0); xterms    = vector (rvec,i,0); rcrc      = vector (rvec,i,0); rie       = vector (rtrm,i,0); riesave   = vector (rtrm,i,0); ccrc      = vector (cvec,i,0); rieest    = vector (cvec,i,0); cie       = vector (ctrm,i,0); init(initbase) = {   local(s,x,y,lastabs, curabs, lastl);   default(format, "g0.15");   B=initbase;   for (s=1,cvec,     x=-1/(cvec*2)+(s/cvec);      /* sexp 1/2 circle, from 0 to 1, 0 to pi, ccrc goes from 1 to 0+i to -1 */     y=exp(Pi*I*x);     ccrc[s]=y;   );   /* linear approximation for base B, first derivative smooth */   lnB = log(B);   rlnB = 1/lnB;   lnlnB = log(rlnB)*rlnB;   cie = (rlnB+lnlnB)/2;   cie = rlnB-lnlnB;   ctrmat = 2;   for (s=3, ctrm, cie[s]=0);   L=0.5;   curabs = 1;   lastabs = 1;   s=120;   while ((curabs0),     lastabs = curabs;     lastl=L;     L=log(L)*rlnB;     curabs = abs(lastl-L);     s--;   );   print ("   base         " B);   print ("   fixed point  " L); /* calculate value for scnt.  This is the number iterations required for the superfunction, isuperf routines   */   y=0.5;   s=0;   while ((abs(y-L)>sfuncprec), y=log(y)*rlnB; s++; ); /* sfuncprec defaults to 1E-12 */   scnt=s;   print ("   iterations   " scnt " used for superf/isuperf");   LxlnB = L*lnB;   rlogLxlnB = 1/log(LxlnB);   idelta =I*imag(ccrc)*1.00; /* *0.99 is a stability term experiment */   imult = exp(-idelta*I*2*Pi);   for (s=1,rvec,     x = -1 - 1/(rvec*2) + (s/rvec); /* riemann full circle, from -0.5 to +0.5, -pi to +pi, +pi to -pi, in slow mode, avoid the singularity in the horizontal direction */     xterms[s]=x;     rcrc[s]=exp(-2*Pi*I*x);   );   renormup = 1; /* used in the recenterup routine, for experiments, try renormup=0 */   recenterup;  /* recenter sexp(z) algorithm so that sexp(-1)=0, sexp(0)=1 */   for (s=1, rtrm, rie[s]=0);   lastprecision=-2;   curprecision=0; } superf(z) = { /* complex valued superfunction for base B, generated from the fixed point L */   local (y);   y=z-scnt;   y=L+((LxlnB)^y); /* LxlnB=L*log(B) */   for (i=1,scnt, y=exp(y*lnB));   return(y); } isuperf(z) = { /* complex valued inverse superfunction for base B, generated from the fixed point L */   local (y,sc);   y=z;   for (i=1,scnt, y=log(y)*rlnB); /* rlnB=1/log(B) */   y=y-L;   sc=(LxlnB)^scnt; /* LxlnB=log(L*lnB) */   y=y*sc;   y=log(y)*rlogLxlnB; /* rlogLxlnB = 1/log(L*lnB); */   return(y); } sexp(z) = { /* sexp approximation from unit circle Taylor series.  Converges nicely for -1 0.5), z=z-1; t++; );   y = z;   s = cie;   for (i=2,ctrmat,     s=s+cie[i]*y;     y=y*z;   );   while ((t<0), s=log(s)*rlnB; t++; );   while ((t>0), s=exp(s*lnB);  t--; );   return(s); } shortsexp(z) = { /* only used for recenterup generation, needs sexp(z) without recenter */   local (i,s,w,y);   y = z;   s = cie;   for (i=2,ctrmat,     s=s+cie[i]*y;     y=y*z;   );   return(s); } riemaprx(w) = {   local(s,x,y,z,tot); /*w = w - idelta;     skip this if m*imult used */   x = real(w);   x = x - floor(x);   /* function only valid if imag(w)>0, what if imag(w)<0 */   y = imag(w);   z = exp(x*2*Pi*I)*exp(-y*2*Pi);   tot = w+rie;   y = z;   for (s=2,rtrmat,     tot=tot+rie[s]*y;     y=y*z;   );   z = superf(tot);   return(z); } riemup(w) = {   local(s,t,x,m,tot);   print (ctrmat " ctrm(s) out of " ctrm " sexp(z) generates " rvec " Riemann samples");   for (s=1,rvec,     x = xterms[s]; /* riemann full circle, from -0.5 to +0.5, -pi to +pi, +pi to -pi, in slow mode, avoid the singularity in the horizontal direction */     rer[s] = isuperf(sexp(x+idelta))-x;   );   rtrmat=rtrm;   m=1;   for (t=1,rtrm,     tot=0;     for (s=1,rvec,       tot=tot+rer[s];       rer[s]=rer[s]*rcrc[s];     );     tot=tot/rvec;     rie[t]=tot;     if ((t>1), m=m*imult; rie[t]=rie[t]*m); /* convert rie back to idelta=0 format */     if ((t>=3),       if (( abs(rie[t]) > abs(rie[t-1]) ), rtrmat=t-1;t=rtrm); /* stability, chaos prevention */     );   );   /* convert rie back to idelta=0 format, allows easy comparisons between rie arrays for different runs */   rie=rie-idelta; } renormsub(z) = {   local(y);   y=shortsexp(z*(0.5+idelta))-exp(lnB*shortsexp(z*(-0.5+idelta)));   return(abs(y)); } recenterup(w) = {   local (delta,newdelta,s,x,y,z,renormK,recenterK);   renormK  = 1.575;  /* approximations for base e, other bases will also work pretty good, by looping */   recenterK= 1.092;  /* approximations for base e, other bases will also work pretty good, by looping */   if (renormup,     /* default is to renormalize each time recenter */     y = renormsub(1);     delta = 2*y;     renorm=1;     for (s=1,25,       x = renormsub(renorm+delta);       z = renormsub(renorm-delta);       if (((x0, recenter=recenter+delta, recenter=recenter-delta);   ); } sexpup(w) = {   local(s,t,x,y,z,tot); /* use the rie array, and regenerate the "cie" approximation */   print (rtrmat " rtrm(s) out of " rtrm " riemaprx(z) generates " cvec " sexp samples");   for (t=1,cvec, rieest[t]=riemaprx(ccrc[t]+centerat));   ctrmat=ctrm;   for (s=1,ctrm,     tot=0;     for (t=1,cvec,       tot=tot+real(rieest[t]);       rieest[t]=rieest[t]*conj(ccrc[t]);     );     cie[s]=tot/cvec;     if ((s%2==0), /* odd Taylor series terms should always be positive */       if ((cie[s]<0), ctrmat=s-2;s=ctrm);     );     if (((s%2) && (s>1)), /* even Taylor series, cie[s-2]<0 => cie[s]<0 */       if (((cie[s]>0) && (cie[s-2]<0)), ctrmat=s-2;s=ctrm);     );   ); } diffriemc(z) = {   local(y);   y=imag(sexp(z)-riemaprx(z)); /*y=sqrt(real(y)*real(y)+imag(y)*imag(y));*/   return(y); } erroravg(w) = {   local(tot,x);   x = sexp(-0.5);   print ("sexp(-0.5)= " x);   tot=0;   for (s=1,200,     z=-0.5-1/400+(s/200);     tot=tot+abs(diffriemc(z+I*imag(ccrc)));   );   tot=-log(tot/200)/log(2);   print (tot " Riemann/sexp binary precision bits I=" idelta);   print (recenter " recenter/renorm " renorm-1);   return(tot); } dumprie(z) = {   local(s,x,y);   s=floor(z);   if ((s>rtrm), s=rtrm);   if ((s<1),    s=1);   x=real(rie[s]);   y=imag(rie[s]);   x=sqrt(x*x+y*y);   return(x); } dumpcie(z) = {   local(s,x);   s=floor(z);   if ((s>ctrm), s=ctrm);   if ((s<1),    s=1);   x= cie[s];   return(x); } riemcontour(z) = {   local(y);   y=sexp(z);   y = isuperf(y)-z;   return(y); } loop(w) = {   local (i,j);   for (i=1,w,     for (j=1,rtrm,riesave[j]=rie[j]);     rtrmatsave=rtrmat;     lastprecision = curprecision;     riemup;     sexpup;     recenterup;     curprecision = erroravg;     if ( ((curprecision-lastprecision)<0),       print ("backup one iteration, and exit loop with earlier results");       for (j=1,rtrm,rie[j]=riesave[j]);       rtrmat=rtrmatsave;       sexpup;       recenterup;       curprecision = erroravg;       i = w;                             /* exit loop, save results         */     );   ); } base_e_2_10(w) = {   init(exp(1));   loop(7);   plotary;   init(2);   loop(7);   plotary;   init(10);   loop(8);   plotary; } qe210(w) = {   init(exp(1));   loop(16);   plotary;   init(2);   loop(16);   plotary;   /* special logic for base 10, using centerat for improved speed, stability */   centerat=-0.5;   init(10);   loop(16);   centerat=0;   sexpup;   recenterup;   erroravg;   plotary;   /* for comparison, base 10 without special recenterat logic */   init(10);   loop(24);   plotary; } plotary(w) = {   local (t);   ploth (t=-0.6,0.6,  diffriemc(t+I*imag(ccrc))  );   print ("plot1, difference Riemann sexp approximation");   ploth (t=-1,  1, sexp(t)                          );   print ("plot2, sexp, Taylor series approximation");   ploth (t= 1, rtrmat,  log(dumprie(t))/log(2)        );   print ("plot3, log_2() Riemann approximation array, shows precision");   ploth (t= 1, ctrmat,  log(abs(dumpcie(t)))/log(2));   print ("plot4, log_2() sexp taylor series array, shows precision");   default(format, "g0.64");   write ("kneser.txt", " ");   write ("kneser.txt", "base        " B);   write ("kneser.txt", "fixed point " L);   write ("kneser.txt", "iterations  " scnt " used for superf/isuperf");   write ("kneser.txt", " ");   write ("kneser.txt", ctrmat " ctrm(s) out of " ctrm " sexp(z) generates " rvec " Riemann samples");   write ("kneser.txt", rtrmat " rtrm(s) out of " rtrm " riemaprx(z) generates " cvec " sexp samples");   write ("kneser.txt", curprecision " Riemann/sexp binary precision bits");   write ("kneser.txt", recenter " recenter/renorm " renorm-1);   write ("kneser.txt", " ");   write ("kneser.txt", "cie sexp Taylor series approximation terms");   for (t=1, ctrmat, write  ("kneser.txt", cie[t]));   write ("kneser.txt", " ");   write ("kneser.txt", "Riemann approximation imaginary delta");   write ("kneser.txt", idelta);   write ("kneser.txt", "rie RiemannCircle / 1-cyclic fourier series");   for (t=1, rtrmat, write  ("kneser.txt", rie[t]));   default(format, "g0.15");   print ("kneser.txt  file for Riemann and sexp Taylor approximation arrays");   /* other interesting or semi-interesting plots              */   /* ploth (t=-1.21, 1.511, imag(riemcontour(t))       );     */   /* ploth (t=-1.21, 1.511, real(riemcontour(t))       );     */   /* ploth (t=-1.01, 1.011, real(riemaprx(t))          );     */   /* ploth (t=-1.01, 0.011, imag(riemaprx(t))          );     */   /* ploth (t=-1.01, 0.011, imag(riemaprx(t+idelta))   );     */   /* ploth (t=-1.01, 0.011, imag(riemaprx(t)-sexp(t))  );     */   /* ploth (t=-1.01, 0.011, diffriemc(t+0.016*I)       );     */   /* ploth (t=-1.01, 0.011, diffriemc(t)               );     */ } help(w) = {   print ("");   print ("help menu for Kneser.gp program");   print ("");   print ("loop(n)      n iteration loops, each loop adds ~8 bits precision");   print ("init(base)   initializes program for base, default loop(exp(1))=e");   print ("base_e_2_10  initializes and loops and generates graphs/data for base e, 2, 10");   print ("sexp(z)      sexp(z), returns Taylor series approximation of sexp(z)");   print ("riemaprx(z)  riemen approximation of sexp(z), imag(z)>=idelta, 0.016*I");   print ("superf(z)    superfunction of z, base B");   print ("isuperf(z)   inverse superfunction of z, base B");   print ("plotary      plots and dumps sexp/Riemann Taylor series arrays to file");   print ("help         print this menu,  'morestuff'  for more special functions");   print ("loop(6)      suggestion n=6 iteration loops");   print (""); } morestuff(w) = {   print ("");   print ("morestuff    morestuff menu for Kneser.gp program");   print ("functions");   print ("shortsexp(z) used for recenterup generation, sexp(z) without recenter, centerat");   print ("sexpup;      updates the sexp(z) function cie array from riemaprx(z) function");   print ("riemup;      updates the riemaprx(z) function rie array from sexp(z) function");   print ("recenterup;  updates the recenter/renorm values");   print ("plotary;     plots arrays/graphs and writes array data to kneser.txt file");   print ("qe210;       initialies/loops/writes/graphs for quad+ ");   print ("             for base e, base 2, and base 10");   print ("default constants");   print ("renormup=1;  used in the recenterup routine, for experiments, try renormup=0");   print ("centerat=0;  improves convergence for base 10 quad precision. centerat=-0.5;");   print ("             after changing centerat=-0.5; type 'sexpup' to reupdate cie array");   print ("riesave      where the copy of rie is stored");   print ("sfuncprec = 1E-15; precision for the super/isuper functions, increased by quad");   print ("idelta=ccrc; set by init(b).");   print (""); } {   init(exp(1));   help; }For the most recent code version: go to the Nov 21st, 2011 thread. RE: fast accurate Knesser sexp algorithm - bo198214 - 08/08/2010 (08/07/2010, 06:53 AM)sheldonison Wrote: I found a fast accurate way to implement the Knesser sexp, and implemented it in Pari-GP Wow this sounds to be a great thing! But before continuing, please explain what you mean by "unit length segment"? In Quote:1) generate the Riemann mapping of the unit length segment Also , his name is "Kneser", not "Knesser"! (08/07/2010, 09:17 PM)sheldonison Wrote: Henryk, can I save the file, "knesser.gp" itself, to the eretrandre website? I added the .gp extension to the allowed uploads. Just make an attachment. RE: fast accurate Knesser sexp algorithm - sheldonison - 08/08/2010 (08/08/2010, 03:54 PM)bo198214 Wrote: Wow this sounds to be a great thing! But before continuing, please explain what you mean by "unit length segment"? In Quote:1) generate the Riemann mapping of the unit length segmentThanks Henryk. I think this could be an important development that I have stumbled upon. I finally had to expand to something more powerful than an excel spreadsheet (and perl programs) though! Its really neat that pari-GP is available on the web as shareware. Hopefully, my thid post in the Mathematical discussion forum helps a little in explaining the theta(z) funciton and the Riemann mapping. The fourier series for theta(z) can be developed from any arbitrary unit length on the real axis of sexp(z), where z>-2. - Sheldon Quote:Also , his name is "Kneser", not "Knesser"! ps. I fixed the pari-GP code to spell Kneser's name correctly! RE: fast accurate Kneser sexp algorithm - nuninho1980 - 08/15/2010 I like your code that it is better and faster on world! wow! my cpu c2d e6600 benchmarked during 24seconds -> init(exp(1)); loop(6). what's your cpu? same or c2q q6600? if eta > base (> 1 or > 0 or > -oo ?) then you add tetration regular for code, ok? updated kneser.gp code - sheldonison - 08/19/2010 I updated the code today, in the second post,. You can download the attachment there. Major changes:higher precision initialization values stored in comments The default, with loop(7) will get double precision results. There are other values in comments at the top of the program. The highest precision results I've tried so far was to generate a 110 term Taylor values accurate to approximately 213 binary bits (terms accurate to >=64 decimal digits), in 3 hours (28 loops, each loop took 6-7 minutes). The corresponding theta(z) function Fourier series has 3000 terms! Pari-gp doesn't support changing precision in a subroutine, so these higher precision values must be manually edited in the kneser.gp file prior to loading the kneser.gp program. More moderate values, with a Taylor series accurate to 100 binary bits, can be calculated in under 10 minutes. These suggested settings have been optimized for best results, and saved in the comments. As before, Type "loop(7);plotary;" to get a 50 term Taylor series output to the kneser.txt file. Use "plotary" to output the current file, and plot some graphs. For higher precision, more loops are required. Each loop gets 7-8 binary bits additional precision. improved convergence I discovered that pari-gp can get really badly confused in how it handles its precision, which is groups of 32 bit chunks. If you compute a Taylor series from right to left, it can suddenly think there's 9 decimal digits less precision, and go chaotic! It took me 3 or 4 days to figure out it gets completely cured, by executing a Taylor series from left to right! This also seems to make the routine run faster. Before figuring out the culprit, I did a few other things to improve convergence in the Taylor series, by looking for ill behaved terms (non decreasing terms, and Jay's uniqueness lemma). Maybe these weren't necessary, but they also make the routine a little faster. I clock the time for "loop(6)" as 21-22 seconds now. Also, I haven't seen any other chaotic cases since making the change. If the code detects precision decreases from one iteration to the next, it will backup and quit since the optimal point has been reached. Theta/Riemann Taylor/Fourier series is now always stored relative to the real axis, even if idelta>0. This allows comparing results for different ideltas, and changing idelta on the fly. Support for centerat for the sexp(z) Taylor series. The default is to centerat=0, but you can develop the Taylor series for sexp(z) at centerat=-0.5, which improves performance for sexp base 10 "init(10)", although with the new better recenter/renorm algorithm, this is less of a big deal. Either way, you can develop the sexp(z) Taylor series at other centers. For the highest precision results, I would like to modify the code to generate more Taylor series terms, to increase the radius of accurate convergence, to a unit circle, instead of within a circle of radius=0.5. Right now, for example the 110 term Taylor series described above is only accurate to 64 decimal digits when abs(z)<=0.5, but only accurate to 35 decimal digits out to a radius of abs(z)<=1. That's because there's only 110 terms, and the 110th term is ~10E-35. That has to do with the limitations and time optimizations of how the Taylor series is generated from the superf(z+theta(z)), which makes it expensive to add Taylor series terms. Later, I will probably post results showing how the Taylor and Fourier series terms converge, and how many terms it takes to get appropriate precision in the theta(z) fourier series so as to generate sufficiently accurate samples for the sexp(z) series. This will justify the "complex conjugate" trick to generate the Taylor series for sexp(z) from superf(z+theta(z)), which is one very small step along the way to showing convergence. Later, I could also add some of nuninho's requests..... - Sheldon (08/15/2010, 12:09 AM)nuninho1980 Wrote: ..... if eta > base (> 1 or > 0 or > -oo ?) then you add tetration regular for code, ok? Strangely, the "superf"/"isuperf" subroutines sort of work for bases= 2000. The high frequency roll off is much more gradual, so we need a little bit of extra precision here. We also need to show that our 75 term sexp(z) Taylor series supports 32 decimal digits of precision at a radius with abs(z)<=0.5. This is for the unit length we use to generate the RiemannCircle/Theta function from. The unit length will be from -0.5+idelta to +0.5+idelta. We want all terms over the unit length accurate to 32 decimal digits. The 75th term is -7.15 E-25. 7.15E-25*(0.5^75)=1.9E-47, which is much much less then 1E-32. The sexp(z) Taylor series has 32 decimal digits of precision out to a radius of about 0.78, which is greater then 1/2. So, all of the parameters are sufficient to get around 32 bits of precision. Observationally, when iterating the "loop(14)", we get about 108 binary bits equal in the sexp(z) algorithm versus the riemaprx(z) algorithm. This is about 33 decimal digits. Graphing the error gainst a golden copy with higher precision shows similar results. - Sheldon RE: fast accurate Kneser sexp algorithm - sheldonison - 10/14/2010 Below, see the attached file, with a new version of the kneser.gp "pari-gp" code, with minor enhancements to speedup the algorithm. The default "kneser.gp" precision is now "\p 67", which supports results accurate to 32 decimal digits, or 105 binary bits. This takes about 14 or so iterations, and takes about a minute to run. This code version updates the size of the arrays as you go, so it is a lot more efficient. The size of the sexp(z) Taylor series and the Riemann mapping series is now updated to match the algorithm's increases in precision, of about 7-8 bits per iteration. Also, the number of exp/log iterations used in the super/isuper super function routine is now also increased each iteration through the loop. The program initializes the same as before. As an example, to calculate the series for other bases, like base 10, type the following, which will execute in just under a minute. For larger bases, like base10, it actually takes more like "loop(17)" iterations to get the full 32 decimal digits of precision supported by "\p 67". The last iteration always takes extra long, since I chose to also increases the size of both series, so as to allow for more terms in the sexp(z) Taylor series result than are required simply for the next loop iteration, but this allows for a larger radius of convergence in the sexp(z) Taylor series, and is controlled by the variable "lastradius". Code:\r kneser.gp; init(10); loop(17); dumparray; /* dumparray dumps the Taylor series to "kneser.txt" and output console */ Much higher precision results are also made possible by this update. Type the following on the command line. In about 30 minutes or so, after 28 iterations, the result is a Taylor series for the base "e" sexp(z), accurate to 64 decimal digits or so. The Taylor series dumped to the file are printed to 64 decimal digits by default, you'll need to update the code in dumparray to get the Taylor series results with more than 64 decimal digits. I know of several further enhancements that are possible, and I've experimented with implementing some of them. I think my focus is going to shift more towards explaining convergence, and other ideas I've had. - Sheldon Code:\r kneser.gp; \p 134; init(exp(1)); loop(28); dumparray; RE: fast accurate Kneser sexp algorithm - nuninho1980 - 10/15/2010 lastest your code for faster. I tried bases 1.5, 1.68, 1.7 and 100 and worked (not work on your older code). it's better. but I get errors after try init(1.6) and init(200). errors here: base 1.6 PHP Code:(15:51) gp > init(1.6);loop(1)   base          1.60000000000000000000000   fixed point   1.77824328112897128908095 + 1.46912237320727818425047*I   Pseudo Period 8.97674947562306210989089 + 1.04999719021153010066853*I3 strm(s) out of 50 sexp(z) generates 36 Riemann samples2 rtrm(s) out of 18 riemaprx(z) generates 12 sexp samples  *** exp: exponent (expo) overflow  base 200 PHP Code:(15:46) gp > init(200)   base          200   fixed point   -0.167639843463924048976342 + 0.375686229934324598314735*I   Pseudo Period 2.73715856720609094587638 + 1.07145668424439506771849*I%1 = 1(15:46) gp > loop(1)3 strm(s) out of 50 sexp(z) generates 36 Riemann samples3 rtrm(s) out of 18 riemaprx(z) generates 12 sexp samplessexp(-0.5)= 0.000719188203196723531532032283176230.3364955007 Riemann/sexp binary precision bits I=0.1305261922*I1.758198575 recenter/renorm 0(15:46) gp > loop(4)1 strm(s) out of 12 sexp(z) generates 36 Riemann samples2 rtrm(s) out of 18 riemaprx(z) generates 12 sexp samples  *** exp: exponent (expo) overflow  I tried base 1.67 but I get 3 bits precision (last loop (from 2 until 5) is less bits precision than previous loop) may you fix this code. 