\p 67;             /* default precision 38 decimal digits */
/* for higher precision, use values of \p 133 or \p 105 */

/* L, B, scnt, strmat, rtrmat, strm, rtrm, sfuncprec, Period idelta are global variables set by init(base), rtrmat set by riemup */
centerat  =     0; /* optional centerat value for sexp   */
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;
idelta    =     0; /* I*imag(exp(Pi*I*(1/(strm*2)))), 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) */
lastradius = 0.84; /* used by loop(t), use 0.84 for last iteration, and don't update precision */

/* dummy global variables initializations, actually set by init and loop */
strm      =    12; /* terms in the sexp Taylor series approximation    */
rtrm      =    18; /* terms in the Riemann mapping Taylor series       */
sfuncprec = 1E-32; /* precision for L fixed point approximation */
Period    =     1; /* period calculated by init */
scnt      =    90; /* updated by init based on sfuncprec        */
strmat    =  strm; /* chaos/overflow prevention */
rtrmat    =  rtrm; /* chaos/overflow prevention */
xterminc  =     0;

rer       = vector (rtrm*2,i,0);
xterms    = vector (rtrm*2,i,0);
rcrc      = vector (rtrm*2,i,0);
rie       = vector (rtrm,i,0);
scrc      = vector (strm,i,0);
rieest    = vector (strm,i,0);
sie       = vector (strm,i,0);

init(initbase) = {
  local(s,x,y,lastabs, curabs, lastl);
  default(format, "g0.24");
  B=initbase;
  for (s=1,strm,
    x=-1/(strm*2)+(s/strm);      /* sexp 1/2 circle, from 0 to 1, 0 to pi, scrc goes from 1 to 0+i to -1 */
    scrc[s]=exp(Pi*I*x);
  );
  /* dumb approximation for base B, first derivative smooth */
  lnB = log(B);
  rlnB = 1/lnB;
  lnlnB = log(rlnB)*rlnB;
  sie[1] = (rlnB+lnlnB)/2;
  sie[2] =  rlnB-lnlnB;
  sie[3] =  0;
  while (sexp(0)>1,recenter=recenter-0.5);
  strmat = 3;
  for (s=4, strm, sie[s]=0);
  L=0.5;
  curabs = 1;
  lastabs = 1;
  s=120; if (B<1.7, s=300);
  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;
  LxlnB = L*lnB;
  rlogLxlnB = 1/log(LxlnB);
  ideltaupd(50); /* works better for initialization renormup */
  xterminc = 0; /* xterminc used to stabilize isuperf for bases<1.7 */
  if ((B<1.7),  xterminc=1);
  if ((B<1.55), xterminc=2);
  if ((B<1.49), xterminc=3);
  centerat=0; if ((B>20), centerat=-0.5);
  for (s=1,rtrm*2,
    x = -1 - 1/(rtrm*4) + (s/(rtrm*2)) + xterminc; /* 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;
  Period = 2*Pi*I/(L*log(B)+log(log(B)));
  print ("   Pseudo Period " Period);
  default(format, "g0.32");
  return(1);
}

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 = sie[1];
  for (i=2,strmat,
    s=s+sie[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 = sie[1];
  for (i=2,strmat,
    s=s+sie[i]*y;
    y=y*z;
  );
  return(s);
}

riemaprx(w) = {
  local(s,y,z,tot);
  z = exp(w*2*Pi*I);
  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 (strmat " strm(s) out of " strm " sexp(z) generates " 2*rtrm " Riemann samples, scnt= " scnt);
  for (s=1,2*rtrm,
    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,2*rtrm,
      tot=tot+rer[s];
      rer[s]=rer[s]*rcrc[s];
    );
    tot=tot/(2*rtrm);
    rie[t]=tot;
    if ((t>1), m=m*imult; rie[t]=rie[t]*m); /* convert rie back to idelta=0 format */
    if ((t>=7),
      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(w) = {
  local (rgt0,argt0,sie1,sie2,sie1s,sie2s);
  rgt0 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
  argt0 = abs(rgt0);
  if (argt0<>0,
    sie1=sie[1];
    sie2=sie[2];

    sie[1]=sie[1]+argt0;
    rgt1 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
    sie1s=(rgt0-rgt1)/argt0;
    sie[1]=sie1;

    sie[2]=sie[2]+argt0;
    rgt2 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
    sie2s=(rgt0-rgt2)/argt0;
    sie[2]=sie2;

    /*   solve for a1, a2, which are the delta values for sie[1], sie[2], respectively */
    /*   a1*imag(sie1s)+a2*imag(sie2s) = imag(rgt0);  */
    /*   a1*real(sie1s)+a2*real(sie2s) = real(rgt0);  */
    sie[1] = sie1 + (imag(rgt0)*real(sie2s) - real(rgt0)*imag(sie2s)) / (real(sie2s)*imag(sie1s) - imag(sie2s)*real(sie1s));
    sie[2] = sie2 + (real(rgt0)*imag(sie1s) - imag(rgt0)*real(sie1s)) / (real(sie2s)*imag(sie1s) - imag(sie2s)*real(sie1s));
  );
  return(rgt0);
}

recenterup(w) = {
  local (delta,newdelta,s,x,y,z,recenterK);
  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 */
    renorm=abs(renormsub);renormsub;
  );

  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,y,z,tot);
/* use the rie array, and regenerate the "sie" approximation */
  print (rtrmat " rtrm(s) out of " rtrm " riemaprx(z) generates " strm " sexp samples");
  for (t=1,strm, rieest[t]=riemaprx(scrc[t]+centerat));
  strmat=strm;
  for (s=1,strm,
    tot=0;
    for (t=1,strm,
      tot=tot+real(rieest[t]);
      rieest[t]=rieest[t]*conj(scrc[t]);
    );
    sie[s]=tot/strm;
    if ((s%2==0), /* odd Taylor series terms should always be positive */
      if ((sie[s]<0), strmat=s-2;s=strm);
    );
    if (((s%2) && (s>1)), /* even Taylor series, sie[s-2]<0 => sie[s]<0 */
      if (((sie[s]>0) && (sie[s-2]<0)), strmat=s-2;s=strm);
    );
  );
}

diffriemc(z) = {
  local(y);
  y=sexp(z)-riemaprx(z);
  return(y);
}

erroravg(w) = {
  local(tot,x);
  x = sexp(-0.5);
  print ("sexp(-0.5)= " x);
  tot=0;
  for (s=1,150,
    z=-0.5-1/300+(s/150);
    tot=tot+abs(diffriemc(z+idelta));
  );
  tot=-log(tot/150)/log(2);
  default(format, "g0.10");
  print (tot " Riemann/sexp binary precision bits I=" idelta);
  print (recenter " recenter/renorm " renorm);
  default(format, "g0.32");
  return(tot);
}

riemcontour(z) = {
  local(y);
  y=sexp(z);
  y = isuperf(y)-z;
  return(y);
}

ideltaupd(c) = {
  local(y);
  y=exp(Pi*I*(1/(2*c)));
  idelta =I*imag(y);
  imult = exp(-idelta*I*2*Pi);
}

changertrmstrm(r,st) = {
  local(s,x,y);
  /* update the rtrm arrays first */
  rtrm=r;
  rer       = vector (rtrm*2,i,0);
  xterms    = vector (rtrm*2,i,0);
  rcrc      = vector (rtrm*2,i,0);
  rie       = vector (rtrm,i,0);

  for (s=1,rtrm*2,
    x = -1 - 1/(rtrm*4) + (s/(rtrm*2)) + xterminc; /* 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);
  );
  /* idelataupd(st) is done outside of this routine */
  lastprecision=curprecision;
  /* now that we have the new array sizes, update the riemaprx rie array */
  riemup;  if (abs(recenter)>0.5,recenter=0);

  /* next, update the sizes of the sie sexp(z) arrays */
  strm=st;
  rieest    = vector (strm,i,0);
  sie       = vector (strm,i,0);
  scrc      = vector (strm,i,0);
  for(s=1, strm,
    x=-1/(strm*2)+(s/strm);
    scrc[s]=exp(Pi*I*x);
  );
  /* and finally update the sie array */
  sexpup;
  /* recenter and renormzlize the sie array */
  recenterup;
  /* calculate the current error term, to be used in the next iteration */
  curprecision = erroravg;
};

loop(t) = {
  local (lp,c,s,x,y,precbits);
  lp=t;
  if ((curprecision<5),
    /* this is the initialization loop, use a reasonable set of initial values */
    ideltaupd(12);
    sfuncprec=1E-8;
    y=0.5;
    s=0;
    while ((abs(y-L)>sfuncprec), y=log(y)*rlnB; s++; ); /* sfuncprec defaults to 1E-8 */
    scnt=s;
    changertrmstrm(18,12);
    lp--;
  );
  precbits=0;
  while ((curprecision>lastprecision) && (lp>0),
    lp--;
    lastprecision = curprecision;
    precbits = curprecision+16;
    if (((lp>0) || (t==1)),
    /* this is not the last iteration, update the precision, and the scnt count used for superf/isuperf */
      x=1/2;
      sfuncprec=2^(-precbits);
      y=0.5;
      s=0;
      while ((abs(y-L)>sfuncprec), y=log(y)*rlnB; s++; );
      scnt=s;
    );
    if ((lp == 0), x=lastradius);  /* use 0.84 for last iteration, and don't update precision */
    /* c will become the new number of terms and sample points in the sexp taylor series array */
    y = precbits*x + 1.5 + log(abs(sie[strmat]))/log(2);
    c = strmat + floor(y);
    /* update idelta to match the value in c */
    ideltaupd(c);
    /* dtermine the value to update for the number of points in the riemaprx taylor series */
    y = (log(2)/log(imult))*(precbits - 2.5 + log(abs(rie[rtrmat]))/log(2));
    rtrm = floor(y);
    /* changertrmstrm actually does the update */
    changertrmstrm(rtrm,c);
  );
}

base_e_2_10(w) = {
/* \p 67; sfuncprec = 1E-32; */
  init(exp(1));
  loop(13);
  dumparray;
  init(2);
  loop(13);
  dumparray;
  init(10);
  loop(16);
  dumparray;
}

dumparray(w) = {
  local (r,s,t);
  default(format, "g0.10");
  print ("sie sexp Taylor series approximation terms");
  if (centerat<>0, print ("sexp Tayler series centerat " centerat));
  default(format, "g0.64");
  write ("kneser.txt", " ");
  write ("kneser.txt", "base          " B);
  write ("kneser.txt", "fixed point   " L);
  write ("kneser.txt", "Pseudo Period " Period);
  if (centerat<>0, write ("kneser.txt", "sexp Tayler series centerat " centerat));
  write ("kneser.txt", " ");
  write ("kneser.txt", "iterations  " scnt " used for superf/isuperf");
  write ("kneser.txt", strmat " strm(s) out of " strm " sexp(z) generates " 2*rtrm " Riemann samples");
  write ("kneser.txt", rtrmat " rtrm(s) out of " rtrm " riemaprx(z) generates " strm " sexp samples");
  write ("kneser.txt", curprecision " Riemann/sexp binary precision bits");
  write ("kneser.txt", recenter " recenter/renorm " renorm);
  write ("kneser.txt", " ");
  write ("kneser.txt", "sie sexp Taylor series approximation terms");
  for (t=1, strmat, write  ("kneser.txt", sie[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]));
  for (t=1, strmat, print  (sie[t]));
  print ("kneser.txt  see file for Riemann and sexp Taylor approximation arrays");

  /* other interesting or semi-interesting plots              */
  /* ploth (t=-0.6,0.6, z=diffriemc(t+idelta*1.2); [real(z), imag(z)]  ); */
  /* print ("plot1, difference Riemann sexp approximation"); */
  /* ploth (t=-1,  1, sexp(t)                          ); */
  /* print ("plot2, sexp, Taylor series approximation"); */
  /* ploth (t= 0, 0.999, s=floor(t*strmat+1);r=floor(t*rtrmat+1);[log(abs(sie[s]))/log(2),log(abs(rie[s]))/log(2)]); */
  /* print ("plot3, log_2() Riemann approximation array, shows precision"); */
  /* ploth (t=-1.21, 1.511, z=riemcontour(t); [real(z), imag(z)]);     */
  /* ploth (t=-1.01, 1.011, z=riemaprx(t); [real(z), imag(z)]); */
  /* 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 7-8 bits precision");
  print ("init(base)   initializes program for base, default init(exp(1))=e");
  print ("base_e_2_10  initializes and loops and generates graphs/data for base e, 2, 10");
  print ("sexp(z)      returns sexp(z), evaluated via the Taylor series approximation");
  print ("riemaprx(z)  Riemann mapping 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 ("dumparray    dumps sexp/Riemann Taylor series arrays to kneser.txt file");
  print ("help         print this menu,  'morestuff'  for more special functions");
  print ("loop(7)      suggestion n=7 iteration loops, ~54 bits precision, ~10 seconds");
  print ("loop(13)     suggestion n=14 iteration loops, ~104 bits precision, ~one minute");
  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 sie 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 ("             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;  default center, other real values are optional, centerat=-0.5");
  print ("             centerat=-0.5 is used for base B>20.  after changing centerat,");
  print ("             type sexpup; to regenerate the Taylor series.");
  print ("sfuncprec  = 1E-32; precision for the super/isuper functions");
  print ("lastradius   radius of conversion for sexp taylor sereies, for last iteration");
  print ("idelta       idelta=imag(scrc[1]); set by init(b).");
  print ("");
}


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