fast accurate Kneser sexp algorithm

[sheldonison]

Here is the most recent code, Nov 12th, 2013:
  kneser.gp (Size: 35.84 KB / Downloads: 1,322)
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

[sheldonison]

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 -I<imag(z)<I.

There's a second routine, that converges in the upper half of the complex plane, for imag(I)>0.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[1]) 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[1] = (rlnB+lnlnB)/2;
  cie[2] = rlnB-lnlnB;
  ctrmat = 2;
  for (s=3, ctrm, cie[s]=0);
  L=0.5;
  curabs = 1;
  lastabs = 1;
  s=120;
  while ((curabs<lastabs) | (s>0),
    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])*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<imag(z)<1.  Real(z) over a very large range   */
  local (i,s,t,y);
  z=z+recenter-centerat;
  t=0;
  while ((real(z)<-0.5), z=z+1; t--; ); /* real z converges over a very large range by mapping back to the unit approximation range */
  while ((real(z)> 0.5), z=z-1; t++; );

  y = z;
  s = cie[1];
  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[1];
  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[1];
  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[1]=rie[1]-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 (((x<y) && (x<z)), renorm=renorm+delta; y=x,
         if (((z<y) && (z<x)), renorm=renorm-delta; y=z));
      delta=delta*2/3;
    );
    z=1; /* renorm update array */
    for (s=2,ctrm,
      z=z*renorm;
      cie[s]=cie[s]*z;  /* renormalize the sexp Taylor series terms */
    );
  );

  delta = 2; y=0; /* recenter update */
  for (s=1,60,
    z=y; y=1-sexp(0);
    if (sign(y) != sign(z), delta=delta/2);
    newdelta = abs(y*recenterK); /* approximation, best for base e */
    if (newdelta<delta,delta=newdelta); /* use the approximation if the approximation is < delta */
    if (y>0, 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[1])));
  );
  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[1]))  );
  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[1]; set by init(b).");
  print ("");
}


{
  init(exp(1));
  help;
}

For the most recent code version: go to the Nov 21st, 2011 thread.

[bo198214]

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

1) generate the Riemann mapping of the unit length segment

Also , his name is "Kneser", not "Knesser"!

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.

[sheldonison]

Wow this sounds to be a great thing!
But before continuing, please explain what you mean by "unit length segment"? In

1) generate the Riemann mapping of the unit length segment

Thanks 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.

\( \theta(z)=\operatorname{isuperf}(\operatorname{sexp}(z))-z \)
- Sheldon

Also , his name is "Kneser", not "Knesser"!

ps. I fixed the pari-GP code to spell Kneser's name correctly!

[nuninho1980]

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?

[sheldonison]

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

.....
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<eta, For example, if you type
init(sqrt(2));
superf(0) You'll get 5.76705325376430. This is the superfunction developed from the upper repelling fixed point, which is "4". For bases<eta, the sexp(z) function doesn't work, nor does the loop(n) function. It would require special unique code to get the function to work for the lower fixed point (=2), since "2" is an attracting fixed point, and the rest of the code is for repelling fixed points. Maybe someday? Also, eta doesn't work at all, and the loop(n) routine overflows for bases<1.7. I haven't debugged to understand why it ges an "exponent (expo) overflow."

[nuninho1980]

Strangely, the "superf"/"isuperf" subroutines sort of work for bases<eta, For example, if you type
init(sqrt(2));
superf(0) You'll get 5.76705325376430. This is the superfunction developed from the upper repelling fixed point, which is "4". For bases<eta, the sexp(z) function doesn't work, nor does the loop(n) function. It would require special unique code to get the function to work for the lower fixed point (=2), since "2" is an attracting fixed point, and the rest of the code is for repelling fixed points. Maybe someday? Also, eta doesn't work at all, and the loop(n) routine overflows for bases<1.7. I haven't debugged to understand why it ges an "exponent (expo) overflow."

ya! I have same error when I tried base=1.7. I try base=sqrt(2) and more bases...

your code already was updated today. wow.

[sheldonison]

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......

Here's an example, settings which I like as the next step up from the term term double precision results, which take less than a minute. Type "init(exp(1));loop(14);" and this next example runs for 7-8 minutes on my laptop, and reports that resulting sexp(z) Taylor series is accurate to 108 binary bits, or just a bit over 32 decimal digits. That precision would be within a radius of abs(z)<=0.5.

\p 67;
ctrm = 75;
cvec = 150;
rtrm = 1000;
rvec = 2000;
sfuncprec = 1E-32;

Here, gp's precision is set to 67 digits. The goal is 32 decimal digits of precision.

The sfuncprec (superfunction precision) is set to 1E-32 digits, to give superf/isuperf results accurate to 32 digits. That is, the isuperf/superf iteration count was set to 235 iterations. If you take the \( \log^{[235]}(0.5) \), the result will be within 10^-32 of the fixed point of L. That is roughly what it takes to get 32 digits of accuracy from the superfunction. We compare exp(L+1E-32) with the fixed point approximation, L+L*1E-32, and we notice that the exponent function has an x^2 term around 1E-64. Thus, the maximum precision we're going to be left with is just a bit over 32 decimal digits.

67-32 leaves 35 digits of precision for gp to work with, or 3 spare bits for rounding errors, which surprisingly seems sufficient! We're asking for a Taylor series with 75 terms. By the Nyquist sampling criteria, we need at least 150 sample points. So cvec (the number of sample terms used to generate the sexp(z) Taylor series) =150. There will be a small sampling error term, due to higher frequencies in the superfunction being sampled. But, for base e sexp(z), the higher frequency terms roll off pretty fast, with each term being roughly half the magnitude of the term before. We want results accurate to 32 decimal digits, so we need 150 sample points, accurate to 32 decimal digits. The sample points will be in the upper half of the complex plane, centered around z=0. The points closest to the real axis are \( \exp(\pi i /300) \) and \( \exp(299\pi i /300) \). This is at \( \Im(z)=0.0105 \)

rtrm = 1000 is driven by the requirement that the Riemann approximation, which is \( \text{superf}(z+\theta(z)) \) be accurate to 32 decimal digits, at \( \Im(z)=0.0105 \). superf is accurate to 32 decimal digits, so the error term is determined by theta(z), which has a radius of convergence at Im(z)=0. The 1000th term in theta(z) has a magnitude of 2.26E-05. The theta(z)=a_n*exp(z*2pi*i). With the known value for Im(z) closest to the real axis, we get magnitude of the a_n term is 6.02*10E-34, which is smaller than the goal of 10E-32. So, using rtrm~=1000 is ok. Sampling for the \( \theta(z) \) function requires 2000 samples by the Nyquist criterion, so we need rvec (sample points for the RiemannCircle/theta(z) function) >= 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

[sheldonison]

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;

[nuninho1980]

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*I
3 strm(s) out of 50 sexp(z) generates 36 Riemann samples
2 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 samples
3 rtrm(s) out of 18 riemaprx(z) generates 12 sexp samples
sexp(-0.5)= 0.00071918820319672353153203228317623
0.3364955007 Riemann/sexp binary precision bits I=0.1305261922*I
1.758198575 recenter/renorm 0
(15:46) gp > loop(4)
1 strm(s) out of 12 sexp(z) generates 36 Riemann samples
2 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.