15,843,937 members
4.00/5 (1 vote)
See more:
Hi,

GeoGrebra, C++ and C# all agree about what cosh(16.612) should be. But when I use a function I wrote myself, which uses an algorithm to calculate the hyperbolic cosine, I get a (very slightly) different result. I've tried setting Visual Studio's floating point handling to precise, strict and fast, but to no avail.

What could the proble be?
This is the one that everyone else agrees about:
`cosh(16.612000) = 8193509.0494933873414993`

This is my first approach:
`_cosh(16.612000) = 8193509.0494933854788542`

And my second:
`ecosh(16.612000) = 8193509.0494933798909187`

The last one uses cosh(x) = 0.5(e^x + e^(-x))
Here is my code:

```#include <cstdio>
#include <iostream>
#include <string>

long double e = 2.718281828459045090795598298427648842334747314453125L;

long double ecosh(long double x) {
//return 0.5 * (exp(x) + exp(-x));
return 0.5L * (pow(e, x) + pow(e, -x));
}

using namespace std;

using decimal = double;

constexpr bool b(int i) { return i > 0; }

static_assert(b(7), "hej");

auto f = [=](double x) {
return sqrt(x*x*x);
};

long double fac(long double x) {
if ((int)x > 1)
return x * fac(x - 1);

return 1;
}

//25/2 + 25/4 - 25/2/2*1.5 + 25/2/2/2*1.5
decimal _sqrt(decimal x, int iterations = 15) {
decimal y = x - x / 2 + x / 4;

for (int i = 1; i < 49; i++) {
y += x / pow(2, i) * 1.5 * pow(-1, i);
}
return y+0.0000000001;
}

long double _cosh2(double x) {
long double y = 0;
for (int i = 1; i < 70; i++) {
y += pow(x, (double)(2.0*i)) / fac(2.0*i);
}
return y+1;
}

long double _sinh2(double x) {
long double y = x;
for (int i = 1; i < 70; i++) {
y += pow(x, (double)(2.0*i+1)) / fac(2.0*i+1);
}
return y + 0;
}

int main(int argc, char** argv)
{
double r = 4.0;
printf("sqrt(%f) = %.25f\n\n", r, _sqrt(r));
printf("sqrt(%f) = %.25f\n\n", r, sqrt(r));
double d = 16.612;
printf("_cosh(%f) = %.16f\n cosh(%f) = %.16f\n", d, _cosh2(d), d, cosh(d));
printf("ecosh(%f) = %.16lf\n", d, ecosh(d));
printf("--------------------------------------\n");
printf("e() = %.75lf\ne   = %.75f\n", exp(1.0L), e);
printf("                       %.50f\n", _sinh2(6.125));
printf("                       %.50f\n", sinh(6.125));
getchar();
}```

What I have tried:

I've tried reconfiguring Visual Studio 2017.
Posted
Updated 24-Nov-18 6:10am
v2
nv3 1-Oct-18 17:01pm
How many digits of precision do you expect when computing in type double?
Dave Kreskowiak 1-Oct-18 21:54pm
This has nothing to do with any settings in Visual Studio. This is all up to your code.

No, I don't have any idea what it's going to be.
[no name] 1-Oct-18 23:14pm
The question is why do you get a different result to others? Apart from purely trivial errors this is a subtle question. Writing a math library is a very complex task and so far you have taken a very superficial look at it. You know nothing of how the other methods are calculating this internally or how close they themselves are to the "correct" result. Particularly you don't seem to have done a thorough analysis yourself. What is pow() up to for example? I would add that the result you get would seem perfectly adequate in confirming the algorithm works.
YvesDaoust 7-Dec-18 12:02pm
Why worry ? Your value is accurate to fifteen significant figures, the maximu that double precision provides.

## Solution 2

I tried your formula in my little test bed which uses VS17 and I got exactly the same results you did. I tried various values for the iteration count from 32 to 128 and all gave the same result. FWIW, above about 127 gives an NAN.

I think you are running into the numerical limits of 64-bit float pointing math. This is compounded by using two iterative calculations to compute the answer. In both calculations a constant value is used which can not be represented exactly (1 and the loop counter). This small inaccuracy in approximation is then magnified with each iteration of the loop. In the end, you are accurate to 15 digits and that is really good since 64-bit values with a 52-bit mantissa are said to have 15 to 17 significant digits of precision according to the wikipedia.

One thing to keep in mind is the long double is treated as a double by VS2017.

## Solution 4

Writing a good mathematical library is very difficult. The library must not only give accurate results, but must also preserve other attributes, such as monotonicity (i.e. if the mathematical function is monotonic, your approximation should also be monotonic), odd/even behaviour (i.e. f(x) == f(-x), or f(x) == -f(-x) for certain functions), etc.

For example:
1. The cosh() function is even, i.e. cosh(x) == cosh(-x). In order to guarantee this, many mathematical libraries will calculate y = x^2, cosh(x) == F(y).
2. The sinh() function is odd, i.e. sinh(x) == -sinh(-x). In order to guarantee this, many mathematical libraries will calculate y = x^2, sinh(x) = x * G(y).

Other issues, such as argument reduction, calculation of the function in a restricted range, and reconstructing the result are much too difficult to go into here. A good book on the subject is Elementary Functions - Algorithms and Implementation | Jean-Michel Muller | Springer[^]. You may also look at the source for a good-quality library at fdlibm[^].

deXo-fan 24-Nov-18 11:12am
I know about the even/odd stuff, but I didn't quite get what you meant by calculating y as x^2 and cosh(x) == F(y).
By F, do you mean an integrated f? The stem function?
Daniel Pfeffer 24-Nov-18 14:51pm
We have an even function, f(x). We wish to develop an approximation, F(x) to f.

When developing an approximation to an even function, you can sometimes get odd components with a very small coefficient. These are due to round off errors in the calculations, and they spoil the evenness of the approximation.

One way to avoid this is to replace 'x' with 'x^2', and calculate a different approximation to f, using the new parameter. This guarantees that the approximation is also even.

## Solution 5

Thank you very much, all of you. I must admit that I wasn't intending to write a math library, I was simply curious. But I'm still intrigued, and I will definitely look at all the libraries you guys have suggested, and that one book suggested by Daniel Pfeffer.

I wish you all the best.

## Solution 3

if you need higher precision you may use some math library like Boost.math and read the wikipedia article. It has also some interesting links.

Remark: it may be better to write
C++
`decimal y = x - x / 2. + x / 4.;`
to avoid rounding problems.

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Top Experts
Last 24hrsThis month
 Patrice T 110 CPallini 105 OriginalGriff 95 merano99 40 Pete O'Hanlon 35
 OriginalGriff 1,798 Richard Deeming 705 Dave Kreskowiak 561 Graeme_Grant 479 Patrice T 285

CodeProject, 20 Bay Street, 11th Floor Toronto, Ontario, Canada M5J 2N8 +1 (416) 849-8900