Cancellation in C++

Question

I am trying to figure out what the problem with the following expression in C++ is:

y=std::log(std::cosh(x));

My first intention was that there might occure a Cancellation due to the cosh(x) part, because it is definde as $\frac{e^x+e^{-x}}{2}$ and the computation of $e^x$ with double x results in Cancellation.

Am I on the right track? Or is there something different that causes Cancellation?

Answer 1

Let me summarize what I wrote in the comments. It is not a complete answer, since the intervals on which to apply each formula still need to be investigated.

It is enough to assume $x\geq0$, since $\cosh(x)$ is even.

The type of cancellation that occurs when evaluating, in finite precision floating point (FP2), the expression

$$\log(\cosh(x))$$

are:

1. When $x$ is large, in which case $\cosh(x)$ makes it even larger but $\log$ would bring the value back down. FP2 is more sparse for larger values. So, one should prevent $\cosh(x)$ making the value large.
2. When $x$ is small the $\cosh(x)$ is close to $1$. This is cool on its own. FP2 is densest near $1$, but then $\log$ becomes close to $0$. In this case it is better to approximate the function $\log(1+x)$ and the function $\cosh(x)-1$ and compose those.

So, for $x$ large one can approximate $\cosh(x)$ by $\frac{e^x}{2}$. Composing with $\log(x)$ one gets $x-\log(2)$.

For $x$ small we can write $\cosh(x)=1+\frac{(e^x-1)^2}{2e^x}$ and compute $\log(\cosh(x))$ by composing $\log(1+x)$ and $\frac{(e^x-1)^2}{2e^x}$. The latter would be computed by approximating $e^x-1$ directly and not by evaluating $e^x$ and subtracting $1$. A C++ implementation could compute $\log(1+x)$ using std::log1p(x) and $e^x-1$ using std::expm1(x).

Finally, one needs to investigate what would be good value $x_1,x_2$ such that one would use

1. the last computation on the interval $[0,x_1]$,
2. the computation std::log(std::cosh(x)) on $(x_1,x_2]$
3. and x-std::log(2) on the $[x_2,\infty]$.