How to find a Newton-like approximation for that function?

I want to find the complex fixpoint $t=b^t$ for real bases $b> \eta = \exp(\exp(-1))$.

added remark: I'm aware that there is a solution using branches of the Lambert-W-function, but I've no Maple/Mathematica and only a rough implementation in Pari/GP for real values. That motivated to try a solution via Newton/Raphson. And to understand and solve such an implementation (which has to deal with derivatives and complex values) is/was then my question here. See also my comment to Fabian's answer below. I seem to have solved it myself for the "principal branch" (see my own answer below) but it is still open for the general case of k'th branch. [end of remark]

What I have is a function depending on a parameter $\beta$ giving the auxiliary values
$$u = \frac{\beta}{ \sin(\beta) } *\exp( i * \beta)$$ $$t=\exp(u)$$ $$b= f(\beta) = \exp(u/t) = \exp(u * \exp(-u))$$

By this I can do an approximation given a base $B$ using binary search. I can find the bounds of an interval taking lower and upper-limit beta's $\beta_l = \epsilon$ and $\beta_u = \pi-\epsilon$ with small epsilons giving the lower and upper bases $b_l$ and $b_u$ respectively. Then comparing $b_m = f(\beta_m)$ where $\beta_m = (\beta_l + \beta_u)/2$ with my given base $B$ I can implement a binary search which approximates $b_m$ to $B$ arbitrarily well and having $\beta_m$ I can reconstruct u and the fixpoint t by the above formula.

However, that binary search needs surprisingly many iterations and I thought, possibly a Newton-like method for that approximation would be more efficient. But since I have complex values involved I do not even see the derivative and even less the formula how to involve that derivative in such an approximation-formula and how to apply this finally to actually do the iterations...

[update 4] moved my own findings into an own answer (as suggested in meta.***)

[update 1] (in this old plot I used the letter s instead of b)

• In your formulas $*$ is used for simple product or for convolution? – mpiktas Mar 3 '11 at 9:31
• @mpiktas: just multiplication – Gottfried Helms Mar 3 '11 at 12:18
• There is also a nice treatize by L. Euler on that question of finding a fixpoint in exponentiation. Unfortunately it seems that there is only the latin version and no translation yet. But it is nice to see the formulae in that old paper ... math.dartmouth.edu/~euler/docs/originals/E489.pdf – Gottfried Helms Mar 3 '11 at 13:30

3 Answers

It might be advantageous to reduce the problem to something which has been studied in the literature: the keyword is Lambert W function. The Lambert W function $W(z)$ is defined via $$z= W(z)e^{W(z)}.$$ On many systems there are already implementations of this function. In terms of the Lambert W function your $t$ is given via the relation $$t = - \frac{W(-\log b)}{\log b}.$$

We need to find a way to calculate $w=W(z)$ for $z<-1/e$. I follow an idea of this page. Let $z=-r$ (with $r>1/e$), $w=R e^{i \theta}$ ($R>0$). Writing $w e^{w} =z$ in polar coordinates, reduces to the set of equation $$r= R e^{R \cos \theta}$$ and $$\pi = \theta + R \sin \theta.$$ This set of equations can be solved for $R>0$ and $\theta \in [0,\pi]$.

The second equation can be rewritten as $R= (\pi - \theta)/\sin \theta$. Plugging this into the first equation yields $$r \sin \theta= (\pi - \theta) e^{(\pi - \theta) \cot \theta}.$$ This last equation can be solved using bisection (we know that the solution is between 0 and $\pi$).

• thanks! Just didn't think of the W-function although I knew of the Galidakis-site you mention. Anyway- what I was doing was (using Pari/GP, which has no W-implementation) searching that approximation and then I got aware, that I even would not know how to implement a Newton-approximation in this case. So that was the original motivation. (For the Lambert W I have only a real-version from some open-source collection; after your hint I asked Wolfram-alpha and got a solution using the "-1"-branch of the W - not usable with my Pari/GP-function) – Gottfried Helms Mar 1 '11 at 20:37
• But solving the equation $r \sin \theta= (\pi - \theta) e^{(\pi - \theta) \cot \theta}$ for $\theta$ with bisection should work...? – Fabian Mar 1 '11 at 20:40
• I think so. Well, meanwhile it has become late here, I'll come back to this tomorrow morning. – Gottfried Helms Mar 1 '11 at 20:58

(As suggested by a thread in meta.stackexchange I moved my own solution (which is however still partial) into a separate answer)

The following solves the problem of finding the Newton-iteration at least for the "principal" value.

Using $c(b) = ln(ln(b))$ extends not only the range for numerical evaluation but also simplifies the whole problem. First I remove the use of the named variable u by a function properly dependend on $\beta$ : $$u(\beta) = \frac{\beta}{ \sin(\beta) } *\exp( i * \beta) = \beta \cot(\beta) + i \beta$$

Then I consider the replacement of the binary search using $v(\beta) = \ln(u(\beta))-u(\beta)$ by the Newton/Raphson-iteration-scheme to approximate to a zero of $$v(\beta,b) = v(\beta) - c(b)$$

Btw., here v has a an interesting powerseries: its coefficients are rational scalings of bernoulli-numbers, resp. zeta at negative arguments: $$v(x)=-1 + \sum_{k=2}^{\infty} (2 i)^k*(k+1)/k! \zeta (1-k) x^k$$

That value c(b) to which I want to approximate defines the constant deviation from v(x) so I include that in the constant term: $$v(x,b)=-1 - c(b) + \sum_{k=2}^{\infty} (2 i)^k*(k+1)/k! \zeta (1-k) x^k$$

The derivative by termwise differentiation is $$dv(x) = v'(x) = \sum_{k=2}^{\infty} (2 i)^k*(k+1)*k/k! \zeta (1-k) x^{k-1}$$ and is the same as of $v(x,b)$

Now the Newton-iteration for finding the root of $v(\beta,b)$ is $$\beta_0 = \frac{\pi}2 \text{ and } \beta_{n+1} = \beta_n - \frac{v(\beta_n,b)}{v'(\beta_n)}$$ and with some n'th iteration I get $\beta_n$ which provides approximations to the parameters $$t \approx t_n = \exp(u(\beta_n))$$ $$B \approx b_n = \exp(\exp(v(\beta_n,B)))$$ where also $$B = t^{\frac1t} \approx t_n \ ^{\frac1{t_n}}$$ and B is the given real base for which I search the complex fixpoint t.

Numerically this iteration converges well with few (say 64) terms; I need now only 4 to 8 iterations to find a well approximated t for base b=4 where I needed hundreds of iterations using the binary search to get accuracy of, say, fourty digits.

(I've not yet tried to extend this to other complex solutions in t, so this is still open)

[update2] The selection of another branch (and meaning of another fixpoint) is as simple as possible: just adjust the search-interval for the iterative function to the appropriate interval of $k*\pi+\epsilon \ldots (k+1)*\pi-\epsilon$ Using the Pari/GP-builtin binary search with the new functions, I get the first seven fixpoints as follows

$$\begin{array} {rrr} k & log(t) & t \ \ (=Fixpoint) & error: \ \ \ 4 - t^{\frac1t} \\ 0 & 0.0886671+1.51223*I & 0.0639598+1.09084*I & 2.04128E-202 \\ 1 & 1.20116+4.44867*I & -0.866450-3.20904*I & 8.36925E-201 \\ 2 & 1.73066+7.63096*I & 1.24841+5.50457*I & 1.49014E-200 \\ 3 & 2.07153+10.8062*I & -1.49429-7.79501*I & 7.14449E-201 \\ 4 & 2.32409+13.9723*I & 1.67648+10.0789*I & 1.16353E-200 \\ 5 & 2.52508+17.1324*I & -1.82146-12.3584*I & 1.01043E-199 \\ 6 & 2.69214+20.2884*I & 1.94197+14.6350*I & 1.85757E-200 \end{array}$$

[update] Hmm, the needed number of terms m in the powerseries of v(x) to get, for instance, accuracy to 40 decimal digits can be fairly well estimated by the following formula for bases b where 2 <= b <= 2000 .
$$m \approx 62.448 * \ln(\ln(b)) + 100 \text{ (formula estimated using Excel) }$$

For a fixed base b one can determine the accuracy to number of digits depending on number of terms n of the power series v(ß). Here is a plot for b=2, b=4, b=16, b=256

For instance, for base B=4 I get using 150 terms in the power series of v(x) the complex fixpoint t to 40 digits using 6 iterations as

t =  0.06395981030134317249921085931098535528761 +
1.090843376539957576373598843800166909991*I
4 - t^(1/t) = 1.67294184828 E-46


You are interested in evaluating the function $\exp(-W_k(-\log\,b))$ for various $k$; we can restrict to nonnegative $k$ here due to the conjugate relationship between positive and negative branches. A truncation of formula 86 from this article can serve as a starting point for approximating $W_k(z)$, which can be further polished with Newton-Raphson or Halley:

$$W_k(z)\approx \log\,z+2\pi ik-\log(\log\,z+2\pi ik)-\log\left(1-\frac{\log\,\log\,z}{\log\,z}\right)\left(1+\frac1{\log\,z-\log\,\log\,z}\right)$$

Your course of action, then is to take $z=-\log\,b$, pick your $k$, evaluate the approximation given above, polish with an iterative method, and finally negate and exponentiate.