prove conjecture; the limit of iterating is $\sqrt{z^2 - 2}$

Question

$$\lim_{n \to \infty} f_n(x)=x-\frac{1}{nx}\;\;\; g(x)=f_n^{on}(x)$$

The conjecture is for values of $|x|>\sqrt{2}$: $g(x) = \sqrt{z^2 - 2}$

This question comes from another matstack question/answer.

If $n=2^m$, then convergence is much quicker by starting with $f_n(x)=x-\frac{1}{nx}$ and then iterating the Taylor/Laurent series for $f(x) \mapsto f(f(x))$ m times. I start by generating the formal Taylor series after moving the fixed point from infinity to zero, $1/f(\frac{1}{x})\;$. Then we start iterating with $f_n$

$$f_n(x) = \frac{x}{1 - x^2/n} = x + \frac{x^3}{n} + \frac{x^5}{n^2} + \frac{x^7}{n^3} + \frac{x^9}{n^4} + \frac{x^{11}}{n^5} ...$$

Using this speedup with $m=\log_2(n)$, one can iterate the Taylor series for $f(x) \mapsto f(f(x))$ m times, rather then iterating $f^{on}$, but the two are of course identical. Its just that otherwise convergence is pretty slow, with n iterations to get accuracy to 1/n.

The formal Taylor series coefficients of $1/g(\frac{1}{x}) = \sqrt{z^2/(1-2z^2)}$ are: $x + x^3 + \frac{3 x^5 }{2 } + \frac{5 x^7 }{2 } + \frac{35 x^9 }{8 } + \frac{63 x^{11} }{8 } + \frac{231 x^{13} }{16 } + \frac{429 x^{15} }{16 } + \frac{6435 x^{17} }{128 } + \frac{12155 x^{19} }{128 } + \frac{46189 x^{21} }{256 } ...$

Empirically, for the $2^n$th iteration starting with $f_n$ in the limit above, the Taylor coefficients are accurate to approximately $O 2^{-n}$, so there is pretty good numerical computation evidence for the conjecture, but I have no idea how to prove it.

EDIT

I (https://math.stackexchange.com/users/39261/mick) will place a bounty for the following problem :

Basicly the inverse :

Suppose we are given $g(x)=\sqrt{x^2-2}$ and we are asked to find $f_n$ such that :

$$\lim_{n \to \infty} \;\;\; g(x)=f_n^{on}(x)$$

How do we solve such problems ??

EDIT

EDIT 2

As Sheldon's comment says $f_n(x)=\sqrt{x^2-\frac{2}{n}}$ is also a solution but I want to find the $f_n$ from the OP : $x - \frac{1}{nx}$ or another $f_n$ that is real-meromorphic on the entire complex plane.

EDIT 2

Answer 1

I'll solve a slightly more general problem with this post - that is, if $f_n(x)=x-\frac{1}{nx}$, what is the limit $$g(x)=\lim_{n\rightarrow\infty}f^{n}_n(x).$$ To do this, we first change the form of $f_n$ to a conjugation of other functions. If we define: $$f(x)=x-\frac{1}x$$ $$\alpha_n(x)=\sqrt{n}\cdot x$$ then $$f_n=\alpha_n^{-1}\circ f \circ \alpha_n$$ The advantage of this is that we can rewrite, where $g_n=f_n^n$ that $$g_n = \alpha_n^{-1}\circ f^n \circ \alpha$$ which is nice because it means we only have to study the iterates of $f$ to figure this out. However, this isn't too bad, since we only need tight bounds on $f$. In particular, consider the following equation: $$\Delta f^n=\frac{-1}{f^n}$$ where we take $\Delta_n f_n$ to be defined as $f^{n+1}-f^n$. I write it this way, because next, we can make a leap to a related differential equation; in particular, let's a function $I_x(n)$ which will approximate $f_n$. In particular, we want the following equations to hold: $$I_x(0)=x$$ $$I_x'(n)=\frac{-1}{I_x(n)}$$ This is a nice separable differential equation and has the solution (for positive $x$): $$I_x(n)=\sqrt{x^2-2n}.$$ The precise connection here is that $f^n(x)$ is basically the result of applying Euler's (forward) method to the same differential equation, with step size $1$. However, notice that we have the identity: $$\frac{I_{\sqrt{k}x}(kn)}{\sqrt{k}}=I_x(n)$$ which tells us that our differential equation scales in a certain way (see remarks at end if you want another way to see this) - but what's remarkable, is this means that $f^{kn}(\sqrt{k}x)$ uses Euler's method to approximate $I_{\sqrt{k}x}(n)$ with step size $1$ - but when we scale this result back by the above equation, we get that $\frac{f^{kn}(\sqrt{k}x)}{\sqrt{k}}$ is actually Euler's method, with step size $\frac{1}k$, approximating $I_x(n)$. As the differential equation is nicely behaved (i.e. Lipschitz continuous), we reach the result that, as we decrease the step size, we get $$\lim_{k\rightarrow\infty}\frac{f^{kn}(\sqrt{k}x)}{\sqrt{k}}=I_{x}(n)$$ setting $n=1$ yields: $$\lim_{k\rightarrow\infty}\frac{f^{k}(\sqrt{k}x)}{\sqrt{k}}=I_{x}(1)$$ and substituting $$g(x)=\sqrt{x^2-2}$$ as desired.

---

To really be clear, it's wise to note that when we're talking about Euler's method, what we're really doing is taking the points given from the scaling operation $\frac{1}{\sqrt{k}}f^{kn}(\sqrt{k}x)$ - which are points $kn$ which are $(n,f)$ pairs: $$\left(0,x\right),\left(\frac{1}{k},\frac{1}{\sqrt{k}}f(\sqrt{k}x)\right),\left(\frac{2}k,\frac{1}{\sqrt{k}}f^2(\sqrt{k}x)\right),\ldots$$ and saying, "Gee, that looks like Euler's method". Though, given the identity on $I_x(n)$ it is clear that this must work, it's good to try an example to be sure. Notice that the second point above ought to be equal to $x+\frac{1}kI_x'(0)$, as would be Euler's method. This ends up expanding as: $$x+\frac{1}kI_x'(0)=\frac{1}{\sqrt{k}}f(\sqrt{k}x)$$ $$x-\frac{1}{kx}=\frac{1}{\sqrt{k}}(\sqrt{k}x-\frac{1}{\sqrt{k}x})$$ $$x-\frac{1}{kx}=1-\frac{1}{kx}$$ Since all pairs of points have this same basic form, we can prove algebraically (and inductively) that $\frac{f^{kn}(\sqrt{k}x)}{\sqrt{n}}$ is actually Euler's method approximating $I_x(n)$ with step size $\frac{1}k$ which does suffice to prove the claim. This warped my mind for an hour, but I think I'm okay now.

---

A particular generalization of this that we can make is that if we let $$f_n(x)=x+\frac{\alpha}{2nx}$$ and, correspondingly $$f(x)=x+\frac{\alpha}{2x}$$ then we have $$\Delta f^n=\frac{-\alpha}{f^n}$$ and correspondingly $$I_x'(n)=\frac{\alpha}{2I_x(n)}$$ which has the solution $$I_x(n)=\sqrt{x^2+\alpha n}$$ and by the same argument, we prove that $$\lim_{n\rightarrow\infty}f_n^n(x)=\sqrt{x^2+\alpha}$$

Answer 2

Pick an $x > \sqrt{2}$ and $n > \dfrac{1}{x^2-2}$ and define a sequence by $x_0 = x$, $x_{k+1} = f_n(x_{k})$ for $k \ge 0$.

Rearrange $x_{k+1} = f_n(x_k) = x_k - \dfrac{1}{nx_k}$, to get $x_k(x_k-x_{k+1}) = \dfrac{1}{n}$.

Clearly, $x_0 \ge \sqrt{x^2}$. Now, suppose that $x_{m-1} \ge \sqrt{x^2-\dfrac{2(m-1)}{n}}$ for some $1 \le m \le n$.

Then, $x_{m-1} \ge \sqrt{x^2-2} > \dfrac{1}{\sqrt{n}}$, and so, $x_m = x_{m-1} - \dfrac{1}{nx_{m-1}} > 0$. Furthermore, we have:

$x^2-x_m^2 = \displaystyle\sum_{k = 0}^{m-1}(x_k^2-x_{k+1}^2) = \sum_{k = 0}^{m-1}\left[(2x_k-(x_k-x_{k+1}))(x_k-x_{k+1})\right]$ $= \displaystyle\sum_{k = 0}^{m-1}\left[2x_k(x_k-x_{k+1}) - (x_k-x_{k+1})^2\right] = \sum_{k = 0}^{m-1}\left[\dfrac{2}{n} - \dfrac{1}{n^2x_k^2}\right] = \dfrac{2m}{n} - \dfrac{1}{n^2}\sum_{k = 0}^{m-1}\dfrac{1}{x_k^2} \le \dfrac{2m}{n}$.

Hence, $x_m^2 \ge x^2 - \dfrac{2m}{n}$. Then, since $x_m > 0$, we have $x_m \ge \sqrt{x^2-\dfrac{2m}{n}}$.

So by induction, we have $x_n \ge \sqrt{x^2-2}$. (Note: all of that was necessary to ensure that $x_n > 0$ , otherwise, $x^2-x_n^2 \le 2$ doesn't imply that $x_n \ge \sqrt{x^2-2}$.)

Then, since $x^2-x_n^2 = 2 - \dfrac{1}{n^2}\displaystyle\sum_{k = 0}^{n-1}\dfrac{1}{x_k^2} \ge 2 - \dfrac{1}{nx_n^2} \ge 2 - \dfrac{1}{n(x^2-2)}$, we have $x_n^2 \le x^2-2 + \dfrac{1}{n(x^2-2)}$, i.e. $x_n \le \sqrt{x^2-2 + \dfrac{1}{n(x^2-2)}}$.

Therefore, $\sqrt{x^2-2} \le f_n^n(x) \le \sqrt{x^2-2 + \dfrac{1}{n(x^2-2)}}$ for all $n > \dfrac{1}{x^2-2}$.

Taking the limit as $n \to \infty$ and squeezing yields $g(x) = \displaystyle\lim_{n \to \infty}f_n^n(x) = \sqrt{x^2-2}$, as desired.

Note that the proof that $g(x) = \sqrt{x^2-2}$ for $x < -\sqrt{2}$ is similar.

Answer 3

@Mick comments: As Sheldon's comment says $f_n(x)=\sqrt{x^2-\frac{2}{n}}$ is also a solution but I want to find the $f_n$ from the OP : $x - \frac{1}{nx}$ ...

The answer to Mick's comments might be as follows, which would be another solution to the Op's question. Here is a formal equation for the Abel function for $g(z)$

$$f(z)=z-\frac{1}{zn}$$ Let $\alpha_f(z)$ be the formal Abel function solution such that $$\alpha_f(f(z)) = \alpha_f(z)+1\;\;\;g(z) = f(z)^{on}$$ Therefore $$ \alpha_g(z) = \frac{\alpha_f(z)}{n};\;\;\; \alpha_g(f(z))=\alpha_g(z)+\frac{1}{n};\;\;\; \alpha_g(g(z))=\alpha_g(z)+1$$

To calculate $\alpha_f(z)$, we use the method of Ecalle, which works for parabolic fixed points. First, we move the fixed point of $f_n$ from infinity to zero. $$ f_n(z) = z-\frac{1}{nz} \mapsto \frac{1}{f_n(1/z)} $$ $$ \frac{1}{f_n(1/z)} = \frac{z}{1 - z^2/n} = z + \frac{z^3}{n} + \frac{z^5}{n^2} + \frac{z^7}{n^3} + \frac{z^9}{n^4} + \frac{z^{11}}{n^5} ...$$

Then using the attatched pari-gp program, one gets the assymptotic equation for the Abel function, $\alpha_f(z)$, and as noted, $\alpha_g(z)=\alpha_f(z)/n$. The first few terms of the formal $\alpha_g(z)$ series are printed below. To get arbitrarily accurate convergence, one would iterate $f^{-1}(z)$ a few times, $\alpha_g(z)=\alpha_g(f^{o-m}(z))+\frac{m}{n}$

$$\alpha_f(z) = \frac{-n}{2z^2} + \frac{\ln(z)}{2} + \frac{-z^2}{8n} + \frac{5z^4}{96n^2} + \frac{-7z^6}{288n^3} + \frac{-1z^8}{1280n^4} + \frac{671z^{10}}{28800n^5} + \frac{-9607z^{12}}{483840n^6} ... $$ $$\alpha_g(z) = \frac{-1}{2z^2} + \frac{\ln(z)}{2n} + \frac{-z^2}{8n^2} + \frac{5z^4}{96n^3} + \frac{-7z^6}{288n^4} + \frac{-1z^8}{1280n^5} + \frac{671z^{10}}{28800n^6} + \frac{-9607z^{12}}{483840n^7} ... $$ $$ \alpha_g(g(z)) = \alpha_g(z)+1$$

For arbitrary fractional iterations of $g(z)$, for finite values of n, we make due with numerical methods. Here, since we moved the fixed point to zero, $g(z) \mapsto \frac{1}{g(1/z)}$ $$g^{oz} = \alpha^{-1}(z)$$ $$g(z)=\alpha_g^{-1}(\alpha_g(z)+1)$$

Conveniently, as n gets arbitrarily large, the Abel function converges, which leads directly to the solution.

$$\alpha(z)=\frac{-1}{2z^2};\;\;\; \alpha^{-1}(z) = -\sqrt{\frac{-1}{2z}}$$

And then with a little algebra $$g(z)=\alpha_g^{-1}(\alpha_g(z)+1);\;=\sqrt{\frac{z^2}{1-2z^2}}$$ Finally, moving the fixed point back to infinity yields $$g(z) \mapsto \frac{1}{g(1/z)} = \sqrt{z^2-2}$$

Here is the pari-gp program I used to calculate the formal Abelseries of $\alpha_g(z)$, which is based on Ecalle's method. See Will Jagy's math overflow post

\ps 24
/* kn is just n, but I use n already so ... */
kx=0; forstep (n=0,12,1,kx=kx+x^(2*n+1)/kn^n);
abelseries(fz,n) = {
  local(i,z,ns,rem,logfzx);
  kabel=0;
  klog=0;
  kfz=fz; /* save fz */
  logfzx=Ser(log(fz/x));
  negterms=1;
  while (polcoeff(fz,negterms+1)==0,negterms++);
  print("terms with negative coeffients= "negterms);
  for (i=-negterms,n,
    if (i==0, klog=acoeff, kabel=kabel+acoeff*x^i);
    rem = Ser(subst(kabel,x,fz) - kabel + klog*logfzx - 1);
    z=polcoeff(rem,i+negterms);
    z=subst(z,acoeff,x);
    ns=-polcoeff(z,0)/polcoeff(z,1);
    kabel=subst(kabel,acoeff,ns);
    klog=subst(klog,acoeff,ns);
  );
  return([kabel,klog]);
}
kz = abelseries(kx,20);
kg = kz/kn;

Answer 4

Mick wrote: "but I want to find the $f_n$ from the OP : $x - \frac{1}{nx}$". This answer gives an analytic function for $g_n(z) = f_n^{on}(z)$, where $$f_n(z)= z-\frac{1}{nz}$$ This analytic function converges everywhere in the complex plane except for at the real axis between $[-\sqrt{2},\sqrt{2}]$. As part of this answer, we assume some details, that for purposes of brevity, will not be proven, but the rest of the post assumes it. Conjecture: all of the zeros and poles of $g_n(x)$ are at the real axis between $[-\sqrt{2},\sqrt{2}]$.

So, the first thing we do is map to the entire complex plane to a unit circle, where the boundary of the unit circle has all of the poles of $g_n(z)$. Assuming the conjecture, that tells us the boundary of the unit circle must be mapped to the the real axis, between $[-\sqrt{2},\sqrt{2}]$. The function $\kappa$ maps the unit circle to $[-\sqrt{2},\sqrt{2}]$,

$$\kappa(z) = \frac{z + 1/z}{\sqrt{2}}$$

Here is a picture of how $\kappa(z)$ maps the complex plane to a unit circle, first we have the unit circle, with circles of various radii shown.

And here, we have the corresponding $\kappa(z)=\frac{z+1/z}{\sqrt{2}}$ for the circles above. The boundary of the unit circle is mapped to the real axis, between $[-\sqrt{2},\sqrt{2}]$.

So then the analytic funtion we are interested in is $$m_n(z) = \frac{1}{g_n(\kappa(z))}$$

Once we have $m_n(z)$, we can generate $g_n(z)$ via the following formula $$g_n(z) = \frac{1}{m_n(\kappa^{-1}(z))}$$

Since the limit as n gets arbitrarily large is proven in other answers, then we have $$ \lim_{n \to \infty} \frac{1}{m_n(\kappa^{-1}(z))} = \sqrt{z^2-2} $$ $$ \lim_{n \to \infty} m_n(z) = \frac{1}{\sqrt{\kappa(z)^2-2}} $$

Interestingly, with a little bit of algebra, we get the formal power series for the limiting value of $m_n(z)$, as n gets arbitrarily large. $$ \frac{1}{\sqrt{\kappa(z)^2-2}} = \frac{1}{\sqrt{\frac{(z+1/z)^2}{2}-2}} = \frac{-z+1/z}{\sqrt{2}} $$ $$ \lim_{n \to \infty} m_n(z) = \sqrt{2}\cdot(z+z^3+z^5+z^7+z^9+...) $$

Notice all of the terms for in the series for the limit are 1, and that this matches the series calculated below, as would be expected. Also, the limiting series converges and is analytic for $\sqrt{z^2-2}$ except at $\pm \sqrt{2}$. At this point,the human and computer Algebra gets a little onerous.. The algorithm starts by using the mapping for the fixed point from infinity to zero, from previous answers. Here: $$f_n(z) \mapsto \frac{1}{r_n(1/z)}$$. $$r_n(z) = \frac{1}{f_n(1/z)} = \frac{1}{\sqrt{1/z-z/n} }= z + \frac{z^3}{n} + \frac{z^5}{n^2} + \frac{z^7}{n^3} + \frac{z^9}{n^4} + \frac{z^{11}}{n^5}...$$ $$f_n^{on}(z) = \frac{1}{r_n^{on}(1/z)}$$ Finally, we substitute in $\kappa(z)=\frac{z+1/z}{\sqrt{2}}$, but we are working with $1/z$ so we substitute $z \mapsto \frac{1}{\kappa(z)}$. Here, the algebra begins to get complicated, due to the apparently unavoidable occurrence of a $\sqrt{2}$ term. Then we generate the $r_n^{o n}(\frac{1}{\kappa(z)})$ function. $$\frac{1}{\kappa(z)} = \sqrt{2}\cdot(z - z^3 + z^5 - z^7 + z^9 - z^{11} ...)$$ $$m_n(z) = r_n^{o n}\left(\frac{1}{\kappa(z)}\right) $$

The next step generates the pattern for all values of n, which I did by using pari-gp's polynomial interpolation, for each term. The first few terms of the analytic mapping function, $m_n(z)$ then are as follows: $$ \frac{m_n(z)}{\sqrt{2}} = z + z^3 + \left( 1-\frac{2}{n} \right) z^5 + \left( 1-\frac{6}{n}+\frac{4}{n^2} \right) z^7 + $$ $$\left( 1-\frac{38}{3n}+\frac{22}{n^2}-\frac{28}{3n^3} \right)z^9 +$$ $$\left( 1-\frac{22}{n}+\frac{230}{3n^2}-\frac{76}{n^3}+\frac{64}{3n^4} \right) z^{11} +$$ $$\left( 1-\frac{172}{5n}+\frac{598}{3n^2}-\frac{376}{n^3}+\frac{752}{3n^4}-\frac{208}{5n^5} \right) z^{13} +$$ $$\left( 1-\frac{748}{15n}+\frac{2202}{5n^2}-\frac{3976}{3n^3}+\frac{1640}{n^4}-\frac{3824}{5n^5}+\frac{288}{5n^6} \right) z^{15} +$$ $$\left( 1-\frac{2404}{35n}+\frac{4302}{5n^2}-\frac{57592}{15n^3}+\frac{7454}{n^4}-\frac{96344}{15n^5}+\frac{10408}{5n^6}-\frac{2272}{35n^7} \right) z^{17} +$$ $$\left( 1-\frac{636}{7n}+\frac{10810}{7n^2}-\frac{430088}{45n^3}+\frac{405658}{15n^4}-\frac{331864}{9n^5}+\frac{67480}{3n^6}-\frac{1575584}{315n^7}+\frac{15488}{35n^8} \right) z^{19} ...$$

Then the assertion is this series, if calculated to an infinite number of terms converges inside a unit circle, and with the addition of $\kappa^{-1}$ mapping shown below, then this is the function for $f_n^{on}(z)$ for all values of n. I don't have a closed form, but I have a pari-gp program that will compute the series, for up to 100 terms or so, and have verified that the series works and converges as expected.

$$f(z) = z+\frac{1}{nz};\;\;\;\; \kappa(z)= \frac{z+1/z}{\sqrt{2}};\;\;\;\;f_n^{on}(z) = \frac{1}{m_n(\kappa^{-1}(z))}$$ $$ \lim_{n \to \infty} m_n(z) = \sqrt{2}\cdot(z+z^3+z^5+z^7+z^9+...) = \frac{\sqrt{2}}{-z+1/z} $$

And also, notice the pattern for the Taylor series coefficients, that would support the limit equation as n gets arbitrarily large, since each term converges to 1 as n gets arbitrarily large. This is expected and required if the conjecture that iterating $f^{on}$ converges to $\sqrt{z^2-2}$ everywhere in the complex plane, except at the real axis between $[-\sqrt{2},\sqrt{2}]$, which is mapped to the boundary of the unit circle. Each Taylor series coefficient is a Laurent series with a leading term of $1-\frac{a_1}{n}-\frac{a_2}{n^2}...$, and all of the Laurent terms eventually decay to zero. Convergence with a finite number of terms depends on the radius, since the series itself is pretty poorly behaved. But with 100 terms, and a radius of 0.9, one can expect accuracy of a few parts in 10,000. Also, the Taylor series terms themselves take a long time to converge, as the iteration count gets arbitrarily large. The chaotic region at the boundary of the unit circle is pretty cool too, and I posted a couple of pictures below.

This is a cool image of the z^99th Taylor series coefficient, $a_{99}z^{99}$. This is a logarithmic plot of $a_{99}$ for n between n10^0 to n=10^6. The term has a 48 term Laurent series, and converges to +1, like all terms, as n gets larger than 5000. But notice how the term oscillates between roughly -1 and +1, before slowly converging to +1 as the iteration count gets larger than 5000.

Here are two more images, showing what $f^{on}(z)$ looks like for a very modest value of n=8, where $f(z)=z+\frac{1}{8z}$. Red is real and green is imaginary. Here we graph $f^{o 8}(\kappa(z))$, at a radius of 0.9, where the posted series is visually indistinguishable with 100 terms. Superimposed on the graph is the limit as n gets arbitrarily large, $\frac{-z+1/z}{\sqrt{2}}$. The second graph is at the unit circle, which is mapped to $[-\sqrt{2},\sqrt{2}]$. This second graph shows the effect of the ~256 chaotic zeros and poles of $f^{o8}(z)$, for even this very modest value of n=8. Superimposed, you can see $\frac{-z+1/z}{\sqrt{2}}$ following the imag(z) path between $\pm i\sqrt{2}$. At the unit circle, $f_n^{on}$ no longer converges as n gets arbitrarily large, but at all radii less than the unit circle, the conjecture is that $f_n^{on}$ does converge as n gets arbitrarily large.

Here is the pari-gp program.

\ps 45
/* \ps 101 works too, just slower */
gon(z,i)={
  local(y,n);
  y=z*1.0;
  if (i==0,i=kt);
  for (n=1,kt,y=y-1/(y*kt));
  return(y);
}
sqrtz2m2(z) = {
  z*sqrt(1-2/z^2);
}
kx=0;
forstep (n=0,50,1,kx=kx+x^(2*n+1)/kn^n);
kx=Ser(kx);
invfn(z) = {
  local(zx,zy);
  zx = z/sqrt(2); /* bernard, the scaling changed ... */
  zy = zx * (1 + sqrt(1-zx^-2));
  zy = 1/zy;
  return(zy);
}
evalzd(z,lzd) = {
  local(y);
  if (lzd==0, lzd=zd);
  y = sqrt(0.5)/subst(lzd,x,invfn(z));
  return(y);
}

genzd(nt,returnzd) = {
  local(n,ln2,oddterms,zf);
  kx=0;
  forstep (n=0,50,1,kx=kx+x^(2*n+1)/kn^n);
  kx=Ser(kx);
  /* this zopt/x = (2/x) maps  */
  zoptsq2 = x; forstep(n=3,99,2,zoptsq2=zoptsq2+x^n);
  zoptsq2=Ser(zoptsq2);
  ln2 = floor(log(nt)/log(2) + 0.5);
  if ((2^ln2)==nt, oddterms=0, oddterms=1);
  if (oddterms==0,
    kt = nt;
    zf=kx;
    for(n=1,ln2, zf=subst(zf,x,zf));
    zf=subst(zf,kn,2^(ln2+1));
  ,
    kt = nt;
    zf=kx;
    for(n=2,nt, zf=subst(zf,x,kx));
    /* the correct equation would be nt, but 2*nt moves sqrt(2)->1 */
    /* here, zf(x^2) =>                                            */
    zf=subst(zf,kn,2*nt);
    /* since we are calculating the ratio, divide by x             */
  );
  /* maps all of the singularities of our function to the unit circle */
  zf = subst(zf,x,Ser(2/(x+1/x)))/2;
  zd = Pol(zf);
  zf = Pol(zf-zoptsq2);
  if (returnzd<>0, return(zf), return(zd));
}

/* this works best with \ps 45 or something like that */
formalzdest(i,returnzd) = {
  local(n,z,lgt,xr,yr,sr);
  if ((i % 2)==0, i=i-1);
  lgt=(i-1)/2;
  xr= vector(lgt); yr=vector(lgt); sr=vector(lgt);
  for (kt=1,lgt, sr[kt]=genzd(kt,returnzd); print(kt););
  if (returnzd<>0, zdest=0;, zdest=x+x^3;);
  print("generating polynomials... Latex version printing below");
  forstep (n=5,i,2,
    for (kt=1,lgt,
      zd=sr[kt];
      xr[kt]=kt;
      yr[kt]=polcoeff(zd,n)*kt^((n-3)/2);
    );
    z=polinterpolate(xr,yr);
    z=subst(z,x,nt);
    z=x^n*z/nt^((n-3)/2);
    zdest=zdest+z;
  );
  for(n=0,i,
    printzn(zdest,n);
  );
  return(zdest);
}

prtpoly(wtaylor,t) = {
  local(s,z);
  if (t==0,t=20);
  print1("{func"kt "=");
  z=polcoeff(wtaylor,0);
  if (z<>0,
    if (real(z)<0, print("");print("      "z), print("");print("        " z));
  );
  for (s=1,t,
    z=polcoeff(wtaylor,s);
    if (z<>0,
      print("");
      if (s>9, print1("+x^" s), print1("+x^ " s));
      print1("* ");
      if (real(z)<0, print1(z), print1(" " z));
    );
  );
  print(" }");
}
wrtpoly(wtaylor,t,name) = {
  local(s,z);
  if (t==0,t=20);
  if (name==0,
    write("joe.txt","{func"kt "=");
  ,
    write("joe.txt","{"name"=");
  );
  z=polcoeff(wtaylor,0);
  if (z<>0, write("joe.txt",z); );
  for (s=1,t,
    z=polcoeff(wtaylor,s);
    if (z<>0, write("joe.txt","+x^" s "* "z); );
  );
  write ("joe.txt"," }");
}

printzn(z,n) = {
local(yz,y);
  yz=subst(polcoeff(z,n),nt,1/nt);
  if (yz<>0,
    yz=subst(yz,nt,x);
    print1("+ z^"n);
    if (yz<>1,
      print1(" \\left( ");
      for (n=0,20,
        y=polcoeff(yz,n);
        if (y<>0,
          if ((n==0) && (y==1),
            print1(1);
          ,
            if (numerator(y)>0, print1("+");, y=-y;print1("-"));
            if (denominator(y)<>1,
              print1("\\frac{"numerator(y)"}{"denominator(y)"n");
            ,
              print1("\\frac{"numerator(y)"}{n");
            );
            if (n<>1,print1("^"n"}"), print1("}"););
          );
        );
      );
      print("\\right)");
    );
  );
}

help(w) = {
print ("genzd(n);       /* gen zd series for n,  speedup if n=2^k   */");
print ("formalzdest(i); /* generate formal zdest series for i terms */");
print ("evalzd(z,zd);   /* if lzd=0, evaluates zd for z             */");
print ("kt;             /* default iteration count, used by gon     */");
print ("gon(z,i);       /* exact evaluates f^on, i==0, kt iteations */");
print ("sqrtz2m2(z);    /* sqrt(z^2-2) with correct sqrt picked     */");
print ("kx;             /* formal series 1/f(1/x) for x-1/(kn*x)    */");
print ("prtpoly(s);     /* prints a polynomial or series            */");
print ("invfn(z);       /* inverse of kappa, (x+1/x)/sqrt(2)        */");
print ("printzn(z,n)    /* prints nth coeffient of series using nt  */");
print ("zoptsq2         /* optimal limit for zd as kt to infinity   */");
print ("");
print ("example");
print ("genzd(64);");
print ("gon(2,64);");
print ("evalzd(2,zd)");
print ("evalzd(2,zd) is the same as sqrt(0.5)/subst(zd,x,invfn(2))");
print ("formalzdest(13);");  /* this generates the series I posted */
}
help();