How can we calculate $\ln(x) - 1$ at $x \approx e$ using standard double precision arithmetic?

Question

This question is inspired by the fact that decent mathematical libraries in programming languages compute sine of double precision approximation of $\pi$ correctly to the last bit (example is in Julia):

> sin(Float64(π))
1.2246467991473532e-16

As a fun fact, it can be used to extract $\pi$ to quadruple precision using double precision sine function.

I was wondering, can we easily implement function $\ln(x) - 1$ that would yield correct values for $x \approx e$? By "easily" I mean without explicitly specifying $e$ to quadruple precision, preferably by just cleverly using available functions from standard double precision math libraries.

Straightforwardly calculating $\ln(x) - 1$ of double precision approximation of $e$ yields $0.0$:

> log(Float64(ℯ)) - 1
0.0

while the actual value is

> log(BigFloat(Float64(ℯ))) - 1
-5.318237706605891370519835446436883266985903310522696816605415323809030383416923e-17

---

Clarification I would like to clarify that the idea is not to explicitly use calculations beyond double machine precision. The trick is not to have a series expansion of the function near $x \approx e$ and expand it in terms of $x - e$, the hard part is computing $x - e$ within the double precision arithmetic if $x$ is already approximately $e$ to double precision.

For comparison, if $x = 3.141592653589793$ (actually, its binary approximation, that's important) then we can compute $x - \pi$ as $x - \pi \approx - \sin x$ and standard math library function will give us the correct result, all $\sim 17$ digits. I wondered, can we get the same with $e$ without explicitly going beyond double precision, computing $\ln(x) - 1$ is an equivalent task. Answers of emacs drives me nuts and Somos do explictly use beyond double precision calculation of $x - e$.

Answer 1

What you can try is expanding the Mercator series of $\ln$. Let $x\approx e$, say $x = e+\Delta x$ with $|\Delta x|\ll 1$:

$$\begin{align} \ln x - 1 = \ln\left(1+\frac{\Delta x}e\right) &= \sum_{k=1}^\infty (-1)^{n+1}\frac{\Delta x^n}{ne^n} \\ &= \frac{\Delta x}{e} - \frac{\Delta x^2}{2e^2} +\frac{\Delta x^3}{3e^3} \mp\cdots \end{align}$$ $\def\artanh{\operatorname{artanh}}$

You might get faster convergence using $\artanh$:

$$\begin{align} \ln x - 1 &= 2\artanh\frac{\Delta x}{2e+\Delta x} \\ &= 2\xi + \frac23\xi^3+\frac25\xi^5+\frac27\xi^7+\cdots \\ &= 2\xi \left(1 + \xi^2\left(\frac13 + \xi^2\left(\frac15+\xi^2\left(\frac17 + \cdots \right)\right)\right)\right)\tag 2\\ \end{align}$$ with $\xi=\Delta x/(2e+\Delta x)$. As indicated by the last line, evaluation is performed using the Horner scheme. Start with the highest power $n$ of $\xi$ such that $$n > 1+\frac{53\ln 2}{1+\ln 2 - \ln \Delta x}$$

For example, if $\Delta x = 10^{-8}$, then $n> 2.82$. The next-greater odd number is $n=3$ so that 2 terms of $(2)$ are enough to achieve IEEE-754 double precision.

Answer 2

You can try to do a Padé approximant. For example, using Mathematica:

PadeApproximant[Log[x] - 1, {x, E, 2}] // Simplify

gives $$f(x)\approx\frac{3(x^2-e^2)}{x^2+4ex+e^2}$$

This is very fast to compute (only simple operations) and reasonably close to the real value but will be off by some amount, it all depends on which precision you need. You can add more terms if needed (I used only quadratic approximant).

Answer 3

The problem may be easier than you think. Use the PARI/GP code

? bp=51; N=2^bp; e=exp(1); x=round(e*N)/N; print(y=x/e-1);
-5.3182377066058912291017433304155628559E-17

which seems to match your expected result.

Use the fact that if $\,x\approx e,\,$ then $\,y:=x/e\!-\!1\approx 0\,$ which implies that $$\log(x)\!-\!1 = \log(x/e) = \log(1\!+\!y) = y \!+\! O(y^2) \approx y = x/e\!-\!1.$$

NOTE the fact that if the bit precision $51$ is replaced by some other integer close to it, then the result will be different although it will still be of the same order of magnitude.

NOTE the Wikipedia article double precision states

Significant precision: 53 bits (52 explictly stored)

The $51$ bits which the PARI/GP computation uses is less than $53$ bits and therefore it does not exceed double precision.

NOTE the Mathematica code

Sin[Pi // N] // InputForm
(* 1.2246467991473532*^-16 *)

matches the Julia result and I can get the same result with Python using numpy, however, so far I can't get Mathematica or Python to produce a similar result for $\,\log(e)-1.$

Answer 4

Since $$\log(x)-1=\sum_{n=1}^\infty (-1)^{n+1}\frac {(x-e)^n}{n\,e^n}$$ the best Padé approximants will be $[p+1,p]$.

For example $$\color{blue}{P_2=t \,\frac {30+21t+t^2 }{30+36t+9t^2 } }\qquad \text{with} \qquad t=\frac {x-e}e$$ whose error is $\frac {t^6}{600}$.

For $x \in \left\{\frac{95 e}{100},\frac{105 e}{100}\right\}$, the maximum absolute error is $3.0\times 10^{-11}$.

For $x \in \left\{\frac{99 e}{100},\frac{101 e}{100}\right\}$, the maximum absolute error is $1.8\times 10^{-15}$.

The next one seems to be very good

$$\color{blue}{P_3=\frac t {12} \,\,\frac{420+510 t+140 t^2+3 t^3 } {35+60 t+30 t^2+4 t^3}}\qquad \text{with} \qquad t=\frac {x-e}e$$ whose error is $\frac {t^8}{9800}$

For $x \in \left\{\frac{95 e}{100},\frac{105 e}{100}\right\}$, the maximum absolute error is $4.8\times 10^{-15}$.

For $x \in \left\{\frac{99 e}{100},\frac{101 e}{100}\right\}$, the maximum absolute error is $6.4\times 10^{-17}$.

Answer 5

I see no other way but to explicitly use in the algorithm the difference between $e$ and its Float64 trunction, $de = \mathrm{Float64}(e) - e \approx -1.4456468917292502 \cdot 10^{-16}$. Then we just use the identity

$$ \ln(x) - 1 = \operatorname{ln1p}((x - e) / e) $$

where $\operatorname{ln1p} = \ln(1+x)$ is a standard library function. In Julia:

function log1m(x)
  de = -1.4456468917292502e-16
  dx = x - ℯ
  return log1p((dx+de) / ℯ)
end

> log1m(ℯ)
-5.318237706605892e-17

This should compute $\ln(x)-1$ over the entire domain to full double precision.