1
$\begingroup$

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?

$\endgroup$
13
  • $\begingroup$ What is the problem? Does the compiler complain or is the resultant valurnorna erong? $\endgroup$ – ghellquist Sep 22 '20 at 15:23
  • $\begingroup$ @ghellquist It has something to do with cancellation. Somehow this code can lead to bad results if executet with double values. $\endgroup$ – Pepsilon7 Sep 22 '20 at 15:33
  • 1
    $\begingroup$ @Pepsilon7 Can you give an example with a particular input value, the expected output, and the actual incorrect output? $\endgroup$ – Aaron Rotenberg Sep 22 '20 at 16:30
  • $\begingroup$ 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$. $\endgroup$ – plop Sep 22 '20 at 16:40
  • 2
    $\begingroup$ @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. $\endgroup$ – plop Sep 22 '20 at 18:34
2
$\begingroup$

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.