Getting more precision in math.h

boyfarrell · Nov 28, 2005

Hello everybody,

I need to calculate hyperbolic sin of large numbers, e.g. sinh(1000). However I start to get infinites after sinh(710). How can I do this?

Here is some code that calculates sinh in sinh(700) region for float, double and long double numbers. Why does long double return infinities at the same point as double does?

Code:

#import <math.h>
#import <stdio.h>

int main(int argc, char *argv[])
{
	float flo;
	printf("\nfloat precision\n");
	
	int i;
	for(i=0; i < 20; i++)
	{
		flo = sinh(700.0+i);
		printf("sinh(%f) = %g\n ",(700.0+i), flo);
	}
	
	double doub;
	printf("\ndouble precision\n");
	for(i=0; i < 20; i++)
	{
	doub = sinh(700.0+i);
	printf("sinh(%lf) = %g\n ",(700.0+i), doub);
	}
	
	long double ldoub;
	long double x;
	printf("\nlong double precision\n");

		for(x=700; x < 720; x++)
		{
			ldoub = sinhl(x);
			printf("sinhl(%Lg) = %Lg\n", x , ldoub);
		}
		

    return 0;
}

Here is the output:

Code:

[Session started at 2005-11-25 17:26:05 +0000.]

float precision
sinh(700.000000) = inf
sinh(701.000000) = inf
sinh(702.000000) = inf
sinh(703.000000) = inf
sinh(704.000000) = inf
sinh(705.000000) = inf
sinh(706.000000) = inf
sinh(707.000000) = inf
sinh(708.000000) = inf
sinh(709.000000) = inf
sinh(710.000000) = inf
sinh(711.000000) = inf
sinh(712.000000) = inf
sinh(713.000000) = inf
sinh(714.000000) = inf
sinh(715.000000) = inf
sinh(716.000000) = inf
sinh(717.000000) = inf
sinh(718.000000) = inf
sinh(719.000000) = inf

double precision
sinh(700.000000) = 5.07116e+303
sinh(701.000000) = 1.37848e+304
sinh(702.000000) = 3.74711e+304
sinh(703.000000) = 1.01857e+305
sinh(704.000000) = 2.76876e+305
sinh(705.000000) = 7.52627e+305
sinh(706.000000) = 2.04585e+306
sinh(707.000000) = 5.5612e+306
sinh(708.000000) = 1.51169e+307
sinh(709.000000) = 4.1092e+307
sinh(710.000000) = 1.117e+308
sinh(711.000000) = inf
sinh(712.000000) = inf
sinh(713.000000) = inf
sinh(714.000000) = inf
sinh(715.000000) = inf
sinh(716.000000) = inf
sinh(717.000000) = inf
sinh(718.000000) = inf
sinh(719.000000) = inf

long double precision
sinhl(700) = 5.07116e+303
sinhl(701) = 1.37848e+304
sinhl(702) = 3.74711e+304
sinhl(703) = 1.01857e+305
sinhl(704) = 2.76876e+305
sinhl(705) = 7.52627e+305
sinhl(706) = 2.04585e+306
sinhl(707) = 5.5612e+306
sinhl(708) = 1.51169e+307
sinhl(709) = 4.1092e+307
sinhl(710) = 1.117e+308
sinhl(711) = inf
sinhl(712) = inf
sinhl(713) = inf
sinhl(714) = inf
sinhl(715) = inf
sinhl(716) = inf
sinhl(717) = inf
sinhl(718) = inf
sinhl(719) = inf

sinh has exited with status 0.

Any ideas? Daniel.

Viro · Nov 28, 2005

Just a quick question. You realize that sin takes radians instead of degrees as an argument?

boyfarrell · Nov 28, 2005

Yeah. These aren't angles. I'm using sinh because it makes my equation tidy, I'm grouping all my exp together.

Any idea on how do I get bigger floating point numbers?

Daniel.

ElDiabloConCaca · Nov 28, 2005

To answer one question, standards dictate:

sizeof(long double) >= sizeof(double)

More often than not, though, their sizes are the same. Under some Linux distros, a long double is 12 bytes, on others, it's 8 (for reference, usually the sizeof(double) == 8 on most machines we're familiar with). Still others, like AIX, it can be 16 bytes.

Chances are, on your system, sizeof(long double) == sizeof(double). You can run a quick check with that code to see for sure. In that case, using a long double won't get you any more precision than a regular double, I believe. That's probably why you're seeing long double top out at the same spot double tops out at.

boyfarrell · Nov 28, 2005

Things are getting strange ... here is the output:

Code:

sizeof(double) = 8
sizeof(long double) = 16

Is there a mistake, I'm calling sinhl from the math.h library rather than sinh,

Code:

Can you spot a mistake in the code?[CODE]	double doub;
	double a;
	printf("\ndouble precision\n");
	for(a=700.0; a < 720.0; a++)
	{
	doub = sinh(a);
	printf("sinh(%lf) = %g\n ", a, doub);
	}
	
	long double ldoub;
	long double x;
	printf("\nlong double precision\n");

		for(x=700.0; x < 720.0; x++)
		{
			ldoub = sinhl(x);
			printf("sinhl(%Lg) = %Lg\n", x , ldoub);
		}
		
		printf("\nsizeof(double) = %d\n", sizeof(double) );
		printf("sizeof(long double) = %d\n", sizeof(long double) );

ElDiabloConCaca · Nov 28, 2005

Hmmm... strange. You can indeed support a long double that is longer than a double. I can't spot the error in the code, but I thought I was on to something when I found the second prototype listed here:

Some places definte the prototype of the sinhl() function as follows:

long double sinhl(long double x);

And others define it like this:

long double sinhl(double x);

So it's not clear to me whether the argument should be a double or a long double. It only makes sense that you would be passing a long double as well, but you never know -- can you try again with your code, but instead pass a double instead of a long double into the sinhl function? I am suspecting this won't have any effect at all, though.

I just ran your code on my machine, and this is what I get when x is defined as a long double:

Code:

double precision
sinh(700.000000) = 5.07116e+303
 sinh(701.000000) = 1.37848e+304
 sinh(702.000000) = 3.74711e+304
 sinh(703.000000) = 1.01857e+305
 sinh(704.000000) = 2.76876e+305
 sinh(705.000000) = 7.52627e+305
 sinh(706.000000) = 2.04585e+306
 sinh(707.000000) = 5.5612e+306
 sinh(708.000000) = 1.51169e+307
 sinh(709.000000) = 4.1092e+307
 sinh(710.000000) = 1.117e+308
 sinh(711.000000) = inf
 sinh(712.000000) = inf
 sinh(713.000000) = inf
 sinh(714.000000) = inf
 sinh(715.000000) = inf
 sinh(716.000000) = inf
 sinh(717.000000) = inf
 sinh(718.000000) = inf
 sinh(719.000000) = inf
 
long double precision
sinhl(700) = 0
sinhl(701) = 0
sinhl(702) = 0
sinhl(703) = 0
sinhl(704) = 0
sinhl(705) = 0
sinhl(706) = 0
sinhl(707) = 0
sinhl(708) = 0
sinhl(709) = 0
sinhl(710) = 0
sinhl(711) = 0
sinhl(712) = 0
sinhl(713) = 0
sinhl(714) = 0
sinhl(715) = 0
sinhl(716) = 0
sinhl(717) = 0
sinhl(718) = 0
sinhl(719) = 0

sizeof(double) = 8
sizeof(long double) = 16

...and this when I define x as a double:

Code:

double precision
sinh(700.000000) = 5.07116e+303
 sinh(701.000000) = 1.37848e+304
 sinh(702.000000) = 3.74711e+304
 sinh(703.000000) = 1.01857e+305
 sinh(704.000000) = 2.76876e+305
 sinh(705.000000) = 7.52627e+305
 sinh(706.000000) = 2.04585e+306
 sinh(707.000000) = 5.5612e+306
 sinh(708.000000) = 1.51169e+307
 sinh(709.000000) = 4.1092e+307
 sinh(710.000000) = 1.117e+308
 sinh(711.000000) = inf
 sinh(712.000000) = inf
 sinh(713.000000) = inf
 sinh(714.000000) = inf
 sinh(715.000000) = inf
 sinh(716.000000) = inf
 sinh(717.000000) = inf
 sinh(718.000000) = inf
 sinh(719.000000) = inf
 
long double precision
sinhl(700) = 5.07116e+303
sinhl(701) = 1.37848e+304
sinhl(702) = 3.74711e+304
sinhl(703) = 1.01857e+305
sinhl(704) = 2.76876e+305
sinhl(705) = 7.52627e+305
sinhl(706) = 2.04585e+306
sinhl(707) = 5.5612e+306
sinhl(708) = 1.51169e+307
sinhl(709) = 4.1092e+307
sinhl(710) = 1.117e+308
sinhl(711) = inf
sinhl(712) = inf
sinhl(713) = inf
sinhl(714) = inf
sinhl(715) = inf
sinhl(716) = inf
sinhl(717) = inf
sinhl(718) = inf
sinhl(719) = inf

sizeof(double) = 8
sizeof(long double) = 16

That doesn't get us any closer to the reason, but it goes to show that my machine apparently doesn't like sinhl taking a long double as an argument... or something. This is indeed strange.

Perhaps it has something to do with the %Lg format?

boyfarrell · Nov 28, 2005

Yeah, no change.

long double sinhl(long double) == long double sinhl(double)

Does C have a 64 bit math library?

Daniel.

Mikuro · Nov 28, 2005

I could be completely wrong, but doesn't %Lf mean long float, and isn't a "long float" the same as a regular double? I'm not sure there's any way to display a 16-byte double in a format string.

ElDiabloConCaca · Nov 28, 2005

%lf means long float... the capital 'L' forces it to a long double. %Le is actually preferred (apparently %g has been deprecated or something -- they recommend using %e).

boyfarrell · Nov 30, 2005

I have run the above code on a 64-bit Linux Workstation and it computes the number correctly! I guess there is some compiler option checked that ignores long doubles. Does anybody have an idea on where I can go from here (apart from going to buy a quad!)

Daniel.

Mikuro · Nov 30, 2005

I don't know if there's any way to get the standard operations to work, but you could use third-party math libraries. I've never done any work like this in C, but I've been doing it in REALbasic for years using Bob Delaney's Extended Plugin, which is based on NTL. It supports integers and floats with arbitrary precision limited only by available memory.

GMP is another possiblility.

Getting more precision in math.h

boyfarrell

Registered

Viro

Registered

boyfarrell

Registered

ElDiabloConCaca

U.S.D.A. Prime

boyfarrell

Registered

ElDiabloConCaca

U.S.D.A. Prime

boyfarrell

Registered

Mikuro

Crotchety UI Nitpicker

ElDiabloConCaca

U.S.D.A. Prime

boyfarrell

Registered

Mikuro

Crotchety UI Nitpicker