\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) */
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 */

/* 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;

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);

init(initbase) = {
  local(s,x,y,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;
  sie[1] = (rlnB+lnlnB)/2;
  sie[2] =  rlnB-lnlnB;
  sie[3] =  0;
  sie[4] =  0;
  if (B<2, sie[1]=1; sie[2]=1);  /* renorminit will fill in sie[2,3] */
  recenter = 0;
  strmat = 4;
  while (sexp(0)>1,recenter=recenter-0.5);
  for (s=5, strm, sie[s]=0);
  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) */
  x=1.0;
  precis=precision(x);
  y=1.0;while(x<>(x+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=sfuncprec;
  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);  /* bernard, initial convergence fails for B>20000 */
  strmmult=3/2; /* if ((B>1000), strmmult=1); */ /* bernard, initial convergence fails for B>1000 using strmult, 35% performance for small bases */
  curprecision=0;
  loopcount=0;
  recenterup;   /* recenter sexp(z) algorithm so that sexp(-1)=0, sexp(0)=1 */
  rtrm=3; for (s=1, rtrm, rie[s]=0);
  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);
}

superf2(z) = {
/* complex valued superfunction for base B, generated from the fixed point L */
  local (L2xlnB,y);
  L2xlnB = L2*lnB;
  y=z+scnt;
  y=L2-((L2xlnB)^y); /* LxlnB=L*log(B) */
  for (i=1,scnt, y=log(y)*rlnB);
  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);
  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,y,z,tot);
  if (imag(w)<0, print ("!!! WARNING, splicesexp(z) much better than riemaprx(z) for imag(z)<-I !!!"));
  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);
}

riemaprxB(w) = {
  local(s,w1,y,y1,z,z1,tot);
  /* Period2/2  = Period/2 - I*imag(rie[1]); */
  /* this is also the crossover value */
  /* 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]);
  /* rossover = Period2/2;                   */ /* desperation, didn't work */
  /* rie[1]=real(rie[1])+Period/2-crossover; */ /* desperation, didn't work */
  w1 = conj(w-crossover)+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),
      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,delt,renorm1,renorm2,sie1,sie2,sie1s,sie2s,denom,s);
  renorm0 = (-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
  for (s=1,2,
    delt = abs(-shortsexp(0.5+idelta) + exp(lnB*shortsexp(-0.5+idelta)));
    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=(renorm0-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=(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[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(renorm0)*real(sie2s) - real(renorm0)*imag(sie2s)) / denom;
        sie[2] = sie2 + (real(renorm0)*imag(sie1s) - imag(renorm0)*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);
}

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 (abs(recenter)>0.5, recenter=0);  /* after initialization, set recenter=0 */
  /* renorm update */
  if (renormup && (loopcount==0), renorminit);
  if (renormup && (loopcount>=1), renorm=abs(renormsub));
  /* recenter update */
  x = slog(1,0);
  recenter=recenter+real(x);
}

sexpup(st) = {
  local(s,t,y,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 ((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 (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;
}

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

slog(z,est) = {
  local (x,y,s,tot,lastyz,curyz);
/* inverse sexp(z), using "est" as initial estimate */
/* could replace sexp with stichsexp if greater complex plane coverage is desired */
/* 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=sexp(est);
  curyz=abs(y-z);
  s=1;
  while ((curyz<lastyz) || (s<3),
    tot=sexp(est+0.01  );
    tot=sexp(est+0.01*I)*-I+tot;
    tot=sexp(est-0.01  )*-1+tot;
    tot=sexp(est-0.01*I)* I+tot;
    tot=100*tot/4;
    est=est+(z-y)/tot;
    lastyz=curyz;
    y=sexp(est);
    curyz=abs(y-z);
    s++;
  );
  if (curyz>0.1, print (curyz " bad result, need better initial est, slog(z,est)"));
  return(est);
}

/* 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))); }

gentaylor(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 ("complex sexp Taylor series centered at " w);
}

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

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=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; 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.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 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 ("gentaylor(w,r)   generates tseries sexp taylor series, centered at w, radius r");
  print ("taprx(z)         evaulates tseries taylor series from gentaylor at z");
  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 optional. centerat=-0.5");
  print ("             for base B>20.  after changing centerat, recenterup");
  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       set precision to 134 digits.  Then init(B) to reinitialize");
  print ("\\p 28        set precision to 28 digits for very fast iterative experiments");
  print ("");
}


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