fast accurate Kneser sexp algorithm

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