\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, =2 */
curprecision=0;
lastprecision=-2;
slowmode  =     0;     /* old slowmode, idelta changes.  new "fastmode" uses idelta=0.12*I, sexpup uses itself */
fastdelta = 12*I/100;  /* default idelta for use when (slowmode==0); must be a rational number, or else there is precision loss for \p 134 */
debugloop =     0;     /* debug loop plot control */
B         = exp(1);
etaB      = exp(1/exp(1));
L         = 0.318 + 1.337*I;
idelta    = 0.1*I; /* 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) */
gennewsexp =    0; /* use the riemB function for bases<eta */
lastradius = 0.65; /* used by loop(t), use 0.65 for last iteration, and don't update precision */
loopradius = 0.50; /* used by loop(t), use 0.50 for other iterations, experiment */
strmmult   =  3/2;
loopcount =     0; /* loopcount */
noboundchk =    0;

/* 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 */
goalbits  =     0; /* calculated by init(B) */
precis    =     0; /* calculated in init(B) */
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;

/* stuff added for cheta */
chterms   = 51;
chdelta   = 100;
invchetr  = 0.006;
usematrix = 1;
chetadlt  = 0;
sxpetadlt = 0;
bias_e    = 0;

global(rer,xterms,rcrc,rie,scrc,rieest,sie);
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);
tseries   = vector (200,i,0);
gie       = vector (200,i,0);

sieinit(w) = {
  sie[1] = (rlnB+lnlnB)/2;
  sie[2] =  rlnB-lnlnB;
  sie[3] =  0;
  if (B<3, sie[1]=1; sie[2]=1);  /* renorminit will fill in sie[2,3] */
  recenter = 0;
  strmat = 3;
  for (s=4, strm, sie[s]=0);
  loopcount=0;
  recenterup;   /* recenter sexp(z) algorithm so that sexp(-1)=0, sexp(0)=1 */
}

init(initbase) = {
  local(s,y,z,lastabs, curabs, lastl);
  default(format, "g0.24");
  if (initbase==0, initbase=exp(1)); B=initbase;
  strm=12;
  /* approximation for base B, first derivative continuous */
  lnB = log(B);
  rlnB = 1/lnB;
  lnlnB = log(rlnB)*rlnB;
  L=0.5;
  curabs = 1;
  lastabs = 1;
  s=120; if (B<1.7, s=300); if (B<1.47, s=1800);
  while ((curabs<lastabs) | (s>0),
    lastabs = curabs;
    lastl=L;
    L=log(L)*rlnB;
    curabs = abs(lastl-L);
    s--;
  );
  Period = 2*Pi*I/(L*log(B)+log(log(B)));
  if (B<etaB,
    L2=0.5;
    curabs = 1;
    lastabs = 1;
    s=300;
    while ((curabs<lastabs) | (s>0),
      lastabs = curabs;
      lastl=L2;
      L2=exp(L2*lnB);
      curabs = abs(lastl-L2);
      s--;
    );
    Period2 = -2*Pi*I/(L2*log(B)+log(log(B)));
  );
  /* precision goal bits for loop(t) */
  z=1.0;
  precis=precision(z);
  y=1.0;while(z<>(z+y),y=y/2);
  goalbits = abs(log(sqrt(y)*50)/log(2))-10;  /* subtract 10 for the last iteration */
  print ("   base          " B);
  print ("   fixed point   " L);
  s=10^-((precis/2)-1.5);
  setsfuncprec(s); /* set sfuncprec; calculate value for scnt.  This is the number iterations required for the superfunction, isuperf routines   */
  LxlnB = L*lnB;
  rlogLxlnB = 1/log(LxlnB);
/*idelta=0.031*I;  */ /* works better for initialization renormup */
  idelta=fastdelta;
  imult = exp(-idelta*I*2*Pi);
  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);
  if ((B<1.47), xterminc=0); /* limited convergence range for isuperf(z), superf/riemaprx are not effected */
  centerat=0; if ((B>20),   centerat=-0.5);
  if ((B>5000), centerat=-0.65);
  strmmult=3/2;
  curprecision=0;
  sieinit;
  rtrm=3; for (s=1, rtrm, rie[s]=0);
  print ("   Pseudo Period " Period);
  default(format, "g0.32");
  if (B<etaB, calczerosf);
  initcheta; /* init base cheta stuff */
  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);
}

superf2(z) = {
/* complex valued superfunction for base B, generated from the fixed point L2 */
/* only valid for B<eta */
  local (y);
  if (B>etaB, print ("superf2 only valid for B<eta"),
    y=z+zerosf2+scnt2;
    y=L2-(L2*lnB)^y; /* LxlnB=L*log(B) */
    for (i=1,scnt2, y=log(y)*rlnB);
    return(y);
  );
}

isuperf2(z) = {
/* complex valued inverse superfunction for base B, generated from the fixed point L */
  local (y,sc,L2xlnB);
  y=z;
  for (i=1,scnt2, y=exp(y*lnB)); /* rlnB=1/log(B) */
  y=y-L2;
  L2xlnB = L2*lnB;
  sc=(L2xlnB)^(-scnt2-zerosf2);   /* LxlnB=log(L*lnB) */
  y=-y*sc;
  y=log(y)/log(L2xlnB);
  return(y);
}

calczerosf(z) = {
  local(s,t,lastt);
  zerosf2=1;
  s=1;
  t=1;
  lastt=1;
  while ((s<10) || abs(t)<abs(lastt),
    lastt=t;
    t = -(superf2(0)-1)/0.6;
    zerosf2 = zerosf2+t;
    s++;
  );
  zerosf1=real(isuperf(superf2(Period2/2)))+Period/2;
}

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);
  if (abs(imag(z))>1,
    if (imag(z)>0,
      print ("!!! WARNING, riemaprx(z) much better than sexp(z) for imag(z)>I !!!"),
      print ("!!! WARNING, splicesexp(z) much better than sexp(z) for imag(z)<-I !!!") );
  );
  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,w1,y,y1,z,z1,tot);
  if (imag(w)<0, print ("!!! WARNING, splicesexp(z) much better than riemaprx(z) for imag(z)<-I !!!"));
  if ((B>etaB) || gennewsexp,
    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);
  );
  if ((B<etaB) && (gennewsexp==0),
    /* generates the sexp from the attracting fixed point via */
    /* kneser mapping from the repelling fixed point */
    /* Period2  = Period - 2*I*imag(rie[1]); */
    /* if imag(w)=Period2/2, imag(riemaprxB(w))=0 */
    /* f(Period2/2 + w) = conj(f(conj(Period2/2 - w))) */
    /*                                      */
    crossover = Period - 2*I*imag(rie[1]); /* crossover = Period2; */
    w1 = conj(w)+crossover;
    z  = exp(w*2*Pi*I);
    z1 = exp(w1*2*Pi*I);
    tot = w+rie[1];
    y = z;
    y1 = z1;
    for (s=2,rtrmat,
      tot=tot+rie[s]*y+conj(rie[s]*y1);
      y=y*z;
      y1=y1*z1;
    );
    z = superf(tot);
    return(z);
  );
}

riemup(r) = {
  local(s,t,x,m,tot);

  if (r==0, r=rtrm);
  rtrm=r;
  rtrmat=r;
  rer       = vector (rtrm*2,i,0);
  xterms    = vector (rtrm*2,i,0);
  rcrc      = vector (rtrm*2,i,0);
  rie       = vector (rtrm,i,0);
  print (strmat " strm(s) out of " strm " sexp(z) generates " 2*rtrm " Riemann samples, scnt= " scnt);

  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;    /* 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 */
    rcrc[s]=exp(-2*Pi*I*x);
    rer[s] = isuperf(precision(sexp(x+idelta),precis))-x;    /* make sure we have full precision before 100's of isuperf log iterations */
  );

  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>=9) && (noboundchk==0),
      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 (renorm0,renorm0t,delt,renorm1,renorm2,sie1,sie2,sie1s,sie2s,denom,s);
  renorm0 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
  for (s=1,2,
    renorm0t = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
    delt = abs(renorm0t);
    if (delt<>0,
      sie1=sie[1];
      sie2=sie[2];

      sie[1]=sie[1]+delt;
      renorm1 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
      sie1s=(renorm0t-renorm1)/delt; /* slope of sie[1] with respect to renorm */
      sie[1]=sie1;

      sie[2]=sie[2]+delt;
      renorm2 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
      sie2s=(renorm0t-renorm2)/delt; /* slope of sie[2] with respect to renorm */
      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(renorm0);  */
      /*   a1*real(sie1s)+a2*real(sie2s) = real(renorm0);  */
      denom = real(sie2s)*imag(sie1s) - imag(sie2s)*real(sie1s);
      if (denom<>0,
        sie[1] = sie1 + (imag(renorm0t)*real(sie2s) - real(renorm0t)*imag(sie2s)) / denom;
        sie[2] = sie2 + (real(renorm0t)*imag(sie1s) - imag(renorm0t)*real(sie1s)) / denom;
      );
    );
  );
  return(renorm0);
}

renorminit(w) = {
  local (renorm0,delt,renorm1,renorm2,sie3,sie2,sie3s,sie2s,denom,s,mys);
/* allows smooth bases at idelta, essential for initializing bases near eta, without large values of recenter */
/* for larger bases, provides a smooth transition initialization at idelta, preserving sie[1] */
/* this allows convergence for bases up to 2 million */
  renorms = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
  delt = 5;
/*sie[1]=1;*/
  mys=10^(-precis+5);
  s=1;
  while ((delt>mys) && (s<100),
    renorm0 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
    delt = abs(renorm0);
    if (delt<>0,
      sie3=sie[3];
      sie2=sie[2];

      sie[3]=sie[3]+delt;
      renorm1 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
      sie3s=(renorm0-renorm1)/delt; /* slope of sie[3] with respect to renorm */
      sie[3]=sie3;

      sie[2]=sie[2]+delt;
      renorm2 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
      sie2s=(renorm0-renorm2)/delt; /* slope of sie[2] with respect to renorm */
      sie[2]=sie2;

      /*   solve for a1, a2, which are the delta values for sie[3], sie[2], respectively */
      /*   a1*imag(sie3s)+a2*imag(sie2s) = imag(renorm0);  */
      /*   a1*real(sie3s)+a2*real(sie2s) = real(renorm0);  */
      denom = real(sie2s)*imag(sie3s) - imag(sie2s)*real(sie3s);
      if (denom<>0,
        sie[3] = sie3 + (imag(renorm0)*real(sie2s) - real(renorm0)*imag(sie2s)) / denom;
        sie[2] = sie2 + (real(renorm0)*imag(sie3s) - imag(renorm0)*real(sie3s)) / denom;
      );
    );
    s++;
  );
  return(renorms);
}

/* renormup=2 for alternative renorminit routine */
recenterup(w) = {
  local (delta,newdelta,s,x,y,z);
  if ((abs(recenter)>0.5) && loopcount, recenter=0);  /* after initialization, set recenter=0 */
  /* renorm update */
  if (renormup,
    if ((renormup==2) && (loopcount>=1) && (centerat==0), sie[1]=1);
    if ((renormup==1) && (loopcount>=1), renorm=abs(renormsub), renorm=abs(renorminit);)
  );
  /* recenter update */
  if ((B>20) && (loopcount==0), x=slog(1,0.25), x=slog(1,0));
  recenter=precision(recenter+real(x),precis);
}

sexpup(st) = {
  local(s,t,x,z,tot,ideltai);
/* use the rie array, and regenerate the "sie" approximation */
  print (rtrmat " rtrm(s) out of " rtrm " riemaprx(z) generates " st " sexp samples");
  if (st==0, st=ceil(strm*(1/strmmult)));
  rieest    = vector (st,i,0);
  scrc      = vector (st,i,0);
  ideltai  = imag(idelta);
  for (t=1,st,
    x=-1/(st*2)+(t/st);
    scrc[t]=exp(Pi*I*x);
    x=scrc[t];
    if (slowmode || (imag(x)>=ideltai),
      rieest[t]=riemaprx(scrc[t]+centerat),
      rieest[t]=precision(sexp(scrc[t]+centerat),precis);  /* make sure sexp has full precision before 100's of isuperf log iterations, riemaprx naturally has precision */
    );
  );
  strm = floor(st*strmmult);
  sie       = vector (strm,i,0);
  strmat=strm;
  for (s=1,strm,
    tot=0;
    for (t=1,st,
      tot=tot+real(rieest[t]);
      rieest[t]=rieest[t]*conj(scrc[t]);
    );
    sie[s]=tot/st;
    if (noboundchk==0,
      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);
      );
    );
  );
}

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)+idelta;
    tot=tot+abs(sexp(z)-riemaprx(z));
  );
  tot=-log(tot/150)/log(2);
  default(format, "g0.10");
  print (tot " Riemann/sexp binary precision bits I=" idelta);
  print (loopcount "=loopcount " recenter " recenter/renorm " renorm);
  if (precis>67, default(format, "g0.64"), default(format, "g0.32"));
  return(tot);
}

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

setsfuncprec(newprec) = {
/* set sfuncprec; calculate value for scnt. */
/* This is the number iterations required for the superfunction, isuperf routines   */
  local(y,s);
  sfuncprec=newprec;
  y=0.5;
  s=0;
  while ((abs(y-L)>sfuncprec), y=log(y)*rlnB; s++; ); /* sfuncprec defaults to 1E-8 */
  scnt=s;
  if (B<etaB,
    y=0.5;
    s=0;
    while ((abs(y-L2)>sfuncprec), y=exp(y*lnB); s++; ); /* sfuncprec defaults to 1E-8 */
    scnt2=s;
  );
}

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

slog(z,est) = {
  local (x,y,s,slop,lastyz,curyz,lest,ly,jt);
/* inverse sexp(z), using "est" as initial estimate */
/* est is an optional parameter, important near the singularity */
/* the 0.01 radius used for the slope means won't converge within 0.01 of singularity */
  jt=0;
  if ((imag(z)==0) && (est==0),
    while (z>B,z=rlnB*log(z);jt++);
    if (z<0,z=exp(z*lnB);jt=-1);
  );
  lastyz=100;
  y=splicesexp(est);
  curyz=abs(y-z);
  lest=est+curyz*0.05;
  ly=splicesexp(lest);
  s=1;
  while ((curyz>sfuncprec) && ((curyz<lastyz) || (s<3)),
    est=precision(est,precis);
    y=precision(y,precis);
    slop=(y-ly)/(est-lest);
    lest=est;
    ly=y;
    est=est+(z-y)/slop;
    lastyz=curyz;
    y=splicesexp(est);
    curyz=abs(y-z);
    s++;
  );
  if (curyz>0.1, print (curyz " bad result, need better initial est, slog(z,est)"));
  return(est+jt);
}

/* split the imag(z) into three cases, for optimal precision in the complex plane */
splicesexp(z) = { if (imag(z)>=0.8, riemaprx(z), if (imag(z)<=-0.8, conj(riemaprx(conj(z))), sexp(z))); }

sexptaylor( w,r) = {
  local(rinv,s,t,x,y,z,tot,t_est,tcrc);
/* outputs tseries taylor series, the complex taylor series for sexp, */
/* centered at "w", radius "r", use stitchsexp for the complex plane */
/* if the taylor series for slog is desired, replace stitchsexp with slog */
  t_est    = vector (240,i,0);
  tcrc     = vector (240,i,0);
  if (r==0,r=1);
  rinv = 1/r;
  for(s=1, 240, x=-1/(120*2)+(s/120); tcrc[s]=exp(Pi*I*x); );

  /* uses the splicesexp to allow for centering anywhere in the complex plane */
  /* if the taylor series for slog is desired, replace stitchsexp with slog here */
  for (t=1,240, t_est[t] = splicesexp(w+r*tcrc[t]); );

  for (s=1,200,
    tot=0;
    for (t=1,240,
      tot=tot+t_est[t];
      t_est[t]=t_est[t]*conj(tcrc[t]);
    );
    tseries[s]=tot/240;
  );
  for (s=2,200, tseries[s]=tseries[s]*(rinv)^(s-1); );
  print ("sexp taylor series; first 60 terms of tseries[1..200] centered at " w);
  if ((imag(w)==0) && (real(w)>-2), tseries=real(tseries), default(format,"g0.22"));
  for (s=1,61,z=tseries[s];if (s>10, print1("a" s-1 "= "), print1("a" s-1 "=  "));if (real(z)<0, print(z), print(" " z)));
}

slogtaylor(w,r,est) = {
  local(rinv,s,t,x,y,z,tot,t_est,tcrc);
/* outputs tseries taylor series, the complex taylor series for sexp, */
/* centered at "w", radius "r", use stitchsexp for the complex plane */
/* if the taylor series for slog is desired, replace stitchsexp with slog */
/* est is an optional initial estimate for slog(w+r) */
  t_est    = vector (240,i,0);
  tcrc     = vector (240,i,0);
  if (r==0,r=1);
  rinv = 1/r;
  for(s=1, 240, x=-1/(120*2)+(s/120); tcrc[s]=exp(Pi*I*x); );

  /* uses the splicesexp to allow for centering anywhere in the complex plane */
  /* if the taylor series for slog is desired, replace stitchsexp with slog here */
  if (est==0, est=slog(w+r));
  for (t=1,240, est=slog(w+r*tcrc[t],est);t_est[t]=est);
  for (s=1,200,
    tot=0;
    for (t=1,240,
      tot=tot+t_est[t];
      t_est[t]=t_est[t]*conj(tcrc[t]);
    );
    tseries[s]=tot/240;
  );
  for (s=2,200, tseries[s]=tseries[s]*(rinv)^(s-1); );
  print ("slog taylor series; first 60 terms of tseries[1..200] centered at " w);
  if (imag(w)==0, tseries=real(tseries), default(format,"g0.22"));
  for (s=1,61,z=tseries[s];if (s>10, print1("a" s-1 "= "), print1("a" s-1 "=  "));if (real(z)<0, print(z), print(" " z)));
}

taprx(z) = {
/* sexp approximation from taylor series generated with the sexptaylor routine */
  local (i,s,y);
  y = z;
  s = tseries[1];
  for (i=2,200,
    s=s+tseries[i]*y;
    y=y*z;
  );
  return(s);
}

/* Cheta, sexp_eta stuff added to kneser.gp */
/* inicheta, cheta, sexpeta, genpoly, */
/* sexp_e because its the complex base change function, too much fun! */

initcheta(w) = {
  /* stuff added for base cheta, sexpeta, invcheta, invsexpeta, sexp_e */
  /* could be a little more agressive here, aiming for 75% precision, option for usematrix=0 */
  chterms   = 2*(floor(precis/2)-8)+1;
  chdelta   = (chterms-1)*2;
  fcheta = genpoly(chterms,chdelta);
  invchetr = imag(cheta(-chdelta+(chterms-1)*I/2));
  invprecis = 10000.*(10^-chterms);
  /* invprecis = 100*abs(chetaerr(I+0.5,chdelta)); */
  /* invprecis = 1E-47; */
  chetadlt=cheta(-chdelta);
  fsexpeta = genpoly(chterms,chdelta,1);
  sxpetadlt=sexpeta(chdelta);
  return(0);
}

cheta(w) = {
  local(i,x,y,e1);
  e1=exp(-1);
  x=w;
  i=chdelta;
  while (real(x)>+1/2,x--;i++);
  while (real(x)<-1/2,x++;i--);
  y = eval(fcheta);
  if (i>0, for (s=1,i, y=exp(y*e1)));
  if (i<0, for (s=1,-i,y=log(y)*exp(1)));
  return(y);
}

sexpeta(w) = {
  local(i,x,y,e,e1);
  e=exp(1);
  e1=exp(-1);
  x=w;
  i=-chdelta;
  while (real(x)>+1/2,x--;i++);
  while (real(x)<-1/2,x++;i--);
  y = eval(fsexpeta);
  if (i>0, for (s=1,i, y=exp(y*e1)));
  if (i<0, for (s=1,-i,y=log(y)*e));
  return(y);
}

genpoly(mterms,mdelta,forsexp) = {
  local(i,ht,xr,yr,initc,e,e1,fl,m25,s25,r25);
  xr  = vector (mterms,i,0);
  yr  = vector (mterms,i,0);
  if (mdelta==0,mdelta=chdelta);
  ht = (mterms-1)/2;
  e=exp(1);
  e1=1/e;
  if (forsexp==0,
    initc=2*e;
    for (s=1,mdelta-ht,initc=log(initc)*e);
    for (i=1,mterms, yr[mterms+1-i]=initc;initc=log(initc)*e);
  );
  if (forsexp,
    initc=1;
    for (s=1,mdelta-ht,initc=exp(initc*e1));
    for (i=1,mterms, yr[i]=initc;initc=exp(initc*e1));
  );
  if (usematrix==1,
    m25=matrix(ht,ht,i,j,i^(2*j));
    s25=matrix(ht,1,i,j,(yr[ht+1+i]+yr[ht+1-i])/2-yr[ht+1]);
    r25=matsolve(m25,s25);
    fl=0;
    for (s=1,ht,fl=fl+r25[s,1]*x^(s*2));
    m25=matrix(ht,ht,i,j,i^(2*j-1));
    s25=matrix(ht,1,i,j,(yr[ht+1+i]-yr[ht+1-i])/2);
    r25=matsolve(m25,s25);
    for (s=1,ht,fl=fl+r25[s,1]*x^(s*2-1));
    fl = fl + yr[ht+1];
  );
  if (usematrix==2,
    m25=matrix(mterms,mterms,i,j,if ((i==1+ht) && (j==1), 1, (i-1-ht)^(j-1)));
    s25=matrix(mterms,1,i,j,yr[i]);
    r25=matsolve(m25,s25);
    fl=0;
    for (s=1,mterms,fl=fl+r25[s,1]*x^(s-1));
  );
  if (usematrix==0,
    for(i=1,mterms,xr[i]=i-1-ht);
    fl = polinterpolate(xr,yr);
  );
  return(fl);
}

invcheta(z,est) = {
  local(i,t,x,mz,e,lastyz,curyz,lest,ly,y);
  e=exp(1);
  mz=z;
  t=0;
  if (est==0,
    x=0;
    mz=z;
    while ((abs(mz-chetadlt)>invchetr) && (t<1000),
      t++;
      mz=log(mz)*e;
    );
    if ((t==1000), print ("too close to e " z " " mz));
    est = 2/(chetadlt/e-1) - chdelta - 2/(mz/e-1);
  );
/* inverse sexp(z), using "est" as initial estimate */
/* est is an optional parameter, important near the singularity */
/* the 0.05 radius used for the slope means won't converge within 0.05 of singularity */
  lastyz=100;
  y=cheta(est);
  curyz=abs(y-mz);
  lest=est+curyz*0.05;
  ly=cheta(lest);
  i=40;
  while ((curyz>invprecis) && ((curyz<lastyz) || i),
    est=precision(est,precis);
    y=precision(y,precis);
    slop=(y-ly)/(est-lest);
    lest=est;
    ly=y;
    est=est+(mz-y)/slop;
    lastyz=curyz;
    y=cheta(est);
    curyz=abs(y-mz);
    i--;
  );
  return(est+t);
}

invsexpeta(z,est) = {
  local(i,t,x,mz,e,e1,lastyz,curyz,lest,ly,y);
  e=exp(1);
  e1=1/e;
  mz=z;
  t=0;
  if (est==0,
    x=0;
    mz=z;
    while ((abs(mz-sxpetadlt)>invchetr) && (t<1000),
      t++;
      mz=exp(mz*e1);
    );
    if ((t==1000), print ("too close to e " z " " mz));
    est = chdelta - 2/(1-sxpetadlt*e1) + 2/(1-mz*e1);
  );
/* inverse sexp(z), using "est" as initial estimate */
/* est is an optional parameter, important near the singularity */
/* the 0.05 radius used for the slope means won't converge within 0.05 of singularity */
  lastyz=100;
  y=sexpeta(est);
  curyz=abs(y-mz);
  lest=est+curyz*0.05;
  ly=sexpeta(lest);
  i=40;
  while ((curyz>invprecis) && ((curyz<lastyz) || i),

    est=precision(est,precis);
    y=precision(y,precis);
    slop=(y-ly)/(est-lest);
    lest=est;
    ly=y;
    est=est+(mz-y)/slop;
    lastyz=curyz;
    y=sexpeta(est);
    curyz=abs(y-mz);
    i--;
  );
  return(est-t);
}

/* base conversion sexp */
/* sexp_e(z)=return(log(log(log(cheta(z+bias_e)/exp(1)-1)))); */
sexp_e(z,w) = {
  local(s,x,y);
  if (bias_e==0,
    bias_e=invcheta((exp(exp(exp(1)))+1)*exp(1));  /* used for base conversion sexp_e */
    bias_e=bias_e-3;
  );
  /* while (real(z+bias_e+w)>4.584,w--); */ /* cheta(z) overflow protection */
  x=cheta(z+bias_e+w)/exp(1)-1;
  if ((w>0) && (imag(z)==0), for (s=1,w,x=log(x)));
  if ((w>0) && (imag(z)<>0),
    y=cheta(z+bias_e+w-1)/exp(1)-1;
    x=log(x);
/*  while (imag(y)>(imag(x)+Pi), x=x+2*Pi*I);  */
/*  while (imag(y)<(imag(x)-Pi), x=x-2*Pi*I);  */
    if (imag(y)>(imag(x)+Pi), x=x+2*Pi*I*floor((imag(y)-imag(x)+Pi)/(2*Pi)));
    if (imag(y)<(imag(x)-Pi), x=x-2*Pi*I*floor((imag(x)-imag(y)+Pi)/(2*Pi)));
    for (s=2,w,x=log(x));
  );
  return(x);
}

/* gfunc/gie/gtaylor/gaprx is a free extra taylor series routine for use by anything you want */
/* gfunc(z) = { cheta(invcheta(z)+0.5); } */
/* gfunc(z) = { sexp_e(z,3); } */
/* gfunc(z) = { invcheta(z) } */
/* gfunc(z) = { sexpeta(z) } */
/* gfunc(z) = { cheta(z) } */
gfunc(z)=cheta(z);

gtaylor( w,r,samples) = {
  local(rinv,s,t,x,y,z,tot,t_est,tcrc,halfsamples);
/* outputs gie taylor series, the complex taylor series for sexp, */
/* centered at "w", radius "r", use stitchsexp for the complex plane */
/* if the taylor series for slog is desired, replace stitchsexp with slog */
/* default halfsamples=240 */
  if (samples==0, samples=240);  /* no matter how many sample points, the default gie series size is 200 halfsamples */
  halfsamples=samples/2;
  t_est    = vector (samples,i,0);
  tcrc     = vector (samples,i,0);
  gcent    = w;
  if (r==0,r=1);
  rinv = 1/r;
  for(s=1, samples, x=-1/(samples)+(s/halfsamples); tcrc[s]=exp(Pi*I*x); );

  /* uses the gfunc to allow for centering anywhere in the complex plane */
  /* if the taylor series for slog is desired, replace stitchsexp with slog here */
  for (t=1,samples, t_est[t] = gfunc(w+r*tcrc[t]); );

  for (s=1,200,
    tot=0;
    for (t=1,samples,
      tot=tot+t_est[t];
      t_est[t]=t_est[t]*conj(tcrc[t]);
    );
    gie[s]=tot/(samples);
  );
  for (s=2,200, gie[s]=gie[s]*(rinv)^(s-1); );
  print ("sexp taylor series; first 60 halfsamples of gie[1..200] centered at " w);
  default(format,"g0.22");
  for (s=1,81,z=gie[s];if (s>10, print1("a" s-1 "= "), print1("a" s-1 "=  "));if (real(z)<0, print(z), print(" " z)));
}

gaprx(z) = {
/* inverse ffunc approximation from unit circle Taylor series.  Converges nicely for -1<imag(z)<1.  Real(z) over a very large range   */
  local (i,s,y);
  z = (z-gcent+recenter);
  y = z;
  s = gie[1];
  for (i=2,200,
    s=s+gie[i]*y;
    y=y*z;
  );
  return(s);
}

invgaprx(z,est) = {
  local (x,y,s,slop,lastyz,curyz,lest,ly);
/* inverse sexp(z), using "est" as initial estimate */
/* est is an optional parameter, important near the singularity */
/* the 0.01 radius used for the slope means won't converge within 0.01 of singularity */
  lastyz=100;
  y=gaprx(est);
  curyz=abs(y-z);
  lest=est+curyz*0.05;
  ly=gaprx(lest);
  s=1;
  while ((curyz>sfuncprec) && ((curyz<lastyz) || (s<3)),
    est=precision(est,precis);
    y=precision(y,precis);
    slop=(y-ly)/(est-lest);
    lest=est;
    ly=y;
    est=est+(z-y)/slop;
    lastyz=curyz;
    y=gaprx(est);
    curyz=abs(y-z);
    s++;
  );
  if (curyz>0.1, print (curyz " bad result, need better initial est, slog(z,est)"));
  return(est);
}

loop(t) = {
  local (rt,st,s,xmlt,y,precbits,lastloop,slowprecision);
  if (loopcount>0, t=t+loopcount);
  lastloop=0;
  if ((loopcount==0),
    if (t==1, lastloop=1); /* t=1, only one loop */
    loopcount++;
    /* this is the initialization loop, use a reasonable set of initial values */
    ideltaupd(12);
    setsfuncprec(1E-8);
    riemup(18); /* calculate a 18 term rie array from the initial sexp approximation */
    sexpup(12);/* calculate a 9 term sie array from the riemaprx */
    recenterup; /* recenter and renormzlize the sie array */
    lastprecision=0;
    curprecision = erroravg; /* calculate the current error term, to be used in the next iteration */
  );
  precbits=0;
  while ((lastloop==0),
    if ((curprecision<lastprecision) && (curprecision<goalbits),  /* if curprecision>goalbits, who cares? */
      loopcount=-1;
      print ("UNEXPECTED LOSS: curprecision<lastprecision. EXITING " curprecision);
      return(0);
/*    lastprecision=0; */
/*    curprecision=curprecision-8; */
    );
    loopcount++;
    if (loopcount==t, lastloop=1);
    precbits = curprecision+16;
    xmlt=loopradius; /* 0.5 is default radius in the loop */
    if (curprecision>goalbits, /* if-then-else */
      lastloop=1;                   /* IF: last iteration */
      xmlt=lastradius,              /* IF: use 0.65 for last iteration, and don't update precision */
      setsfuncprec(2^(-precbits));  /* ELSE: this is not the last iteration, update the precision, and the scnt count used for superf/isuperf */
    );
    /* st will become the new number of terms and sample points in the sexp taylor series array */
    y = precbits*xmlt + 1.5 + log(abs(sie[strmat]))/log(2);
    st = strmat + floor(y);
    if (st<12, st=12);
    ideltaupd(st); /* update idelta to match the value in st */
    /* determine the value to update for the number of points in the riemaprx taylor series */
    y = floor((log(2)/log(imult))*(precbits - 2.5 + log(abs(rie[rtrmat]))/log(2)));
    riemup(y);
    sexpup(st);
    recenterup; /* recenter and renormzlize the sie array */
    lastprecision=curprecision;
    curprecision = erroravg; /* calculate the current error term, to be used in the next iteration */
    if (debugloop,
      if (loopcount>=debugloop,
      /*ploth(t=-0.6,0.6,x=t+idelta*1.1;y=sexp(x)-riemaprx(x);[real(y),imag(y)]);*/
        ploth(t=-0.6,0.6,x=t+idelta*1.1;y=sexp(x)-riemaprx(x);z=isuperf(sexp(x))-isuperf(riemaprx(x));[real(y),imag(y),real(z),imag(z)]);
      /*ploth(t=-0.6,0.6,x=t+idelta*1.1;z=isuperf(sexp(x))-isuperf(riemaprx(x));[real(z),imag(z)]);*/
      );
    );
  );
  if (precis>67, default(format, "g0.64"));
  return(1);
}

base_e_2_10(w) = {
  init(exp(1));
  loop;
  dumparray;
  init(2);
  loop(0);
  dumparray;
  init(10);
  loop(0);
  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", loopcount "=loopcount " recenter " recenter/renorm " renorm);
  write ("kneser.txt", " ");
  write ("kneser.txt", "sie sexp Taylor series approximation terms");
  if (loopcount==-1, write ("kneser.txt", "UNEXPECTED PRECISION LOSS.  EXITED EARLY, curprecision<lastprecision " curprecision));
  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=t+idelta;z=sexp(z)-riemaprx(z);[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, z=t+0.16*I;z=sexp(z)-riemaprx(z);[real(z), imag(z)] );  */
}

help(w) = {
  print ("");
  print ("help menu for Kneser.gp program, 'morestuff' for second menu");
  print ("");
  print ("loop(n)      n iteration loops, each loop adds 7-8 bits precision");
  print ("init(base)   initializes program for base, default init; initializes base 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.012*I");
  print ("splicesexp(z) sexp/riemaprx(z) valid everywhere in the complex plane");
  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 ("cheta(z)     base eta functions: cheta(z) sexpeta(z) invcheta(z) invsexpeta(z)");
  print ("help         print this menu,  'morestuff'  for more special functions");
  print ("loop or loop(0), loop until optimal precision, ~104 bits, takes ~30 seconds");
  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, based on renormup");
  print ("new functions");
  print ("slog(z,est)      inverse sexp(z), using optional "est" as initial estimate");
  print ("sexptaylor(w,r)  generates tseries sexp taylor series, centered at w, radius r");
  print ("slogtaylor(w,r)  generates tseries slog taylor series, centered at w, radius r");
  print ("taprx(z)         evaluates tseries sexptaylor/slogtaylor series at z");
  print ("gfunc(z)         set to whatever you want, gtaylor(z,r) genseries gaprx(z)");
  print ("superf2(z)   only for B<eta, sexp(z) from attracting fixed point of base B");
  print ("             ");
  print ("default constants");
  print ("slowmode=0   if slowmode=1, the idelta value is updated, otherwise 0.,12*I");
  print ("renormup=1;  used in the recenterup routine, for experiments, try renormup=0");
  print ("centerat=0;  default center, other real values are used for larger bases");
  print ("gennewsexp=1 for B<etaB, generate the kneser mapping for newsexp");
  print ("gennewsexp=0 for B<etaB, generate the kneser mapping for sexp");
  print ("lastradius   radius of conversion for sexp(z) function, for last iteration");
  print ("\\p 67        default 67 digits initial precision, allows sexp to 32 digits");
  print ("\\p 134 \\p 28 higher precision to 134 digits.  very fast precision 28 digits");
  print ("");
}


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