My cosh versus their cosh

Question

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 1-Oct-18 10:38am

deXo-fan

Updated 24-Nov-18 5:10am

KarstenK

v2

Add a Solution

Comments

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.

4 solutions

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.

Posted 24-Nov-18 5:11am

deXo-fan

Add a Solution

Add your solution here

Treat my content as plain text, not as HTML

Preview 0

…

Existing Members

Sign in to your account

...or Join us

Download, Vote, Comment, Publish.

Your Email
Password
Forgot your password?

Your Email
This email is in use. Do you need your password?
Optional Password

I have read and agree to the Terms of Service and Privacy Policy
Please subscribe me to the CodeProject newsletters

When answering a question please:

Read the question carefully.
Understand that English isn't everyone's first language so be lenient of bad spelling and grammar.
If a question is poorly phrased then either ask for clarification, ignore it, or edit the question and fix the problem. Insults are not welcome.
Don't tell someone to read the manual. Chances are they have and don't get it. Provide an answer or move on to the next question.

Let's work to help developers, not make them feel stupid.

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

Rick York · Accepted Answer · 2018-10-01T20:32:00

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.

Daniel Pfeffer · Accepted Answer · 2018-10-03T04:10:00

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[^].

KarstenK · Accepted Answer · 2018-10-02T00:13:00

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.

You also should read this article about quadruple precision calculations because it has some example code.

Remark: it may be better to write

C++

decimal y = x - x / 2. + x / 4.;

to avoid rounding problems.

My cosh versus their cosh

4 solutions

Solution 2

Solution 4

Solution 5

Solution 3

Add your solution here

Preview 0