fast accurate Kneser sexp algorithm

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;

}