# Cancellation in C++

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?

• What is the problem? Does the compiler complain or is the resultant valurnorna erong? – ghellquist Sep 22 '20 at 15:23
• @ghellquist It has something to do with cancellation. Somehow this code can lead to bad results if executet with double values. – Pepsilon7 Sep 22 '20 at 15:33
• @Pepsilon7 Can you give an example with a particular input value, the expected output, and the actual incorrect output? – Aaron Rotenberg Sep 22 '20 at 16:30
• std::cosh(x) is not computed as $\frac{e^x+e^{-x}}{2}$. For example, in the range $0\leq x\leq \ln(2)/2$ it is computed as $1+\frac{(e^x-1)^2}{2e^x}$, where the $e^x-1$ is not computed using std::exp and subtracting $1$, but using std::expm1. In the range $\ln(2)/2\leq x\leq 22$ is computed as $\frac{e^x+1/e^x}{2}$. From $22\leq x\leq \ln(\operatorname{maxdouble})$ it is computed as $e^x/2$. – plop Sep 22 '20 at 16:40
• @AaronRotenberg The problematic cases are when $\cosh(x)$ is large and when it is close to $1$. For example, $\log(\cosh(123456))$, when computed as std::log(std::cosh(123456)) overflows, while it is about $123455.3$, which is nowhere near to overflow. Similarly, $\log(\cosh(0.0000000000001234))$, when computed as std::log(std::cosh(0.0000000000001234)) it gives 0.0, while it is about $7.61378\times10^-{27}$, which is nowhere near to underflow. Those examples, when computed as in the comments above, give results close to the actual values. – plop Sep 22 '20 at 18:34

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]$$.