High precision doubles in C++

[Sassy Rabbit]

What is the best/most efficient way to handle high precision, where by high, I mean a precision greater than 10^-18, in C++?

Help is much appreciated in advance

-Sassy

 

[CTho9305]

Is this a homework assignment? Did you look at the ranges that ordinary floats can express?

edit: In what range do you need 10^-18 precision? Near 0, near 1, or near 1 billion?

 

[Sassy Rabbit]

Originally posted by: CTho9305
Is this a homework assignment? Did you look at the ranges that ordinary floats can express?

edit: In what range do you need 10^-18 precision? Near 0, near 1, or near 1 billion?

Nope - not a homework assignment, so you won't be giving me an answer - it's for a simulation code I am working with. I am iterating the solution to a matrix equation and I want to ensure that the largest error between the solution at step n-1 and step n is on the order of 10^-18 before i stop the iteration. I need the ability to have numbers that are on the order of 1000-100 with at least 18 decimal places behind them. Long doubles are only getting me a solution with a precision of 10^-17.

 

[Sassy Rabbit]

Originally posted by: CTho9305
Is this a homework assignment? Did you look at the ranges that ordinary floats can express?

edit: In what range do you need 10^-18 precision? Near 0, near 1, or near 1 billion?

BTW - I do appreciate the help

 

[CTho9305]

Originally posted by: Sassy Rabbit

Originally posted by: CTho9305
Is this a homework assignment? Did you look at the ranges that ordinary floats can express?

edit: In what range do you need 10^-18 precision? Near 0, near 1, or near 1 billion?

Nope - not a homework assignment, so you won't be giving me an answer - it's for a simulation code I am working with. I am iterating the solution to a matrix equation and I want to ensure that the largest error between the solution at step n-1 and step n is on the order of 10^-18 before i stop the iteration. I need the ability to have numbers that are on the order of 1000-100 with at least 18 decimal places behind them. Long doubles are only getting me a solution with a precision of 10^-17.

Sounds like a problem. I don't know of any solutions off-hand.... and any solution will likely involve emulating the FP, which will be drastically slower than using natively-supported data types. Are you sure 64 significant bits really isn't enough precision? I'm having a hard time thinking of a case where an error of 1 part in 10^19 (10 sextillion?) would matter... that's more precision than the distance from earth to the sun in microns (1.5*10^17).

 

[Sassy Rabbit]

Originally posted by: CTho9305

Originally posted by: Sassy Rabbit

Originally posted by: CTho9305
Is this a homework assignment? Did you look at the ranges that ordinary floats can express?

edit: In what range do you need 10^-18 precision? Near 0, near 1, or near 1 billion?

Nope - not a homework assignment, so you won't be giving me an answer - it's for a simulation code I am working with. I am iterating the solution to a matrix equation and I want to ensure that the largest error between the solution at step n-1 and step n is on the order of 10^-18 before i stop the iteration. I need the ability to have numbers that are on the order of 1000-100 with at least 18 decimal places behind them. Long doubles are only getting me a solution with a precision of 10^-17.

Sounds like a problem. I don't know of any solutions off-hand.... and any solution will likely involve emulating the FP, which will be drastically slower than using natively-supported data types. Are you sure 64 significant bits really isn't enough precision? I'm having a hard time thinking of a case where an error of 1 part in 10^19 (10 sextillion?) would matter... that's more precision than the distance from earth to the sun in microns (1.5*10^17).

So I have this Fortan77 code that I inherited. Because dynamic memory allocation is a must (absent as far as I am aware in Fortan77 - the code itself has to be recompiled if the number of particles changes) and objects and the C++ vector class lend themselves really well to the particular application I am converting the code to C++ before I start to expand the functionality. In the original Fortran code, they have this level of precision and I was hoping to be able to maintain it - as I believe it has to do with the long term stability (energy conservation) of the N-body simulation.

 

[lousydood]

Have you tried using a Fortran77 to C compiler like f2c?

 

[Sassy Rabbit]

Originally posted by: lousydood
Have you tried using a Fortran77 to C compiler like f2c?

I'm not exactly sure what this would accomplish - just compiling the source so it was compatible with C++?

 

[CTho9305]

Ah, that's unfortunate. I can imagine that it would be nice if you could verify your results by using the same precision as the original code. Have you tried comp.lang.fortran to see what scientific-computing gurus think? They're probably your best bet (maybe not that particular newsgroup...but some place with HPC-type people).

 

[aCynic2]

You could try BCD. Not exactly the best solution, but you're trying to work outside of C/C++'s native ability.

There is one idea I use for mapping space (I mean astronomical space). It uses the concept of base number systems, but you can't express a decimal natively with it. It's meant to overcome the granularity of floating point types in C/C++.

If you recall, a base 10 number (an example) is expressed naturally as: 1024. But, you can also think of it as:

1*10^3 + 0*10^2 + 2*10^1 + 4*10^0 = 1024.

Decimals are expressed as (example): .014

0*10^-1 + 1*10^-2 + 4*10^-3


Now, I came up with a base c number system for mapping a galaxy, where c is the approximate speed of light in a vacuum or 3,000,000,000. It's conversion works like so:

n3*c^4+n2*c^2+n1*c^1+n0*c^0

The fractionals for this is: d1*c^1+d2*c^-2...

This solution allows me to map our entire galaxy down to the meter in the data space of 3 long ints for each axis.

 

[lousydood]

Well, the f2c compiler comes with libraries which should give the Fortran-like math support.

Otherwise I think you're stuck with interval arithmetic or some other way to track error in computation. No fun. This is assuming that the Fortran code even deals with this stuff correctly.

Libraries like the GNU Multi Precision may offer the kind of numbers you want, but they won't be fast.

 

[Sassy Rabbit]

Originally posted by: CTho9305
Ah, that's unfortunate. I can imagine that it would be nice if you could verify your results by using the same precision as the original code. Have you tried comp.lang.fortran to see what scientific-computing gurus think? They're probably your best bet (maybe not that particular newsgroup...but some place with HPC-type people).

Actually - what you said earlier about the numbers being very large - close to a billion or very small did help - i recast the problem a little bit to get the numbers ranging from -5 to 5 and that bought me the precision I needed. Now if I can get that fortran code to compile so I can compare it

Thanks again.

 

[Sassy Rabbit]

Originally posted by: aCynic2
You could try BCD. Not exactly the best solution, but you're trying to work outside of C/C++'s native ability.

There is one idea I use for mapping space (I mean astronomical space). It uses the concept of base number systems, but you can't express a decimal natively with it. It's meant to overcome the granularity of floating point types in C/C++.

If you recall, a base 10 number (an example) is expressed naturally as: 1024. But, you can also think of it as:

1*10^3 + 0*10^2 + 2*10^1 + 4*10^0 = 1024.

Decimals are expressed as (example): .014

0*10^-1 + 1*10^-2 + 4*10^-3


Now, I came up with a base c number system for mapping a galaxy, where c is the approximate speed of light in a vacuum or 3,000,000,000. It's conversion works like so:

n3*c^4+n2*c^2+n1*c^1+n0*c^0

The fractionals for this is: d1*c^1+d2*c^-2...

This solution allows me to map our entire galaxy down to the meter in the data space of 3 long ints for each axis.

This is quite a cool idea - and while i did get my problem working I am keeping this in mind for possible use in the future.

 

[Sassy Rabbit]

Originally posted by: lousydood
Well, the f2c compiler comes with libraries which should give the Fortran-like math support.

Otherwise I think you're stuck with interval arithmetic or some other way to track error in computation. No fun. This is assuming that the Fortran code even deals with this stuff correctly.

Libraries like the GNU Multi Precision may offer the kind of numbers you want, but they won't be fast.

OK the really funny part about this, not that you would personally know, is that the Fortran code won't currently compile - bleh. Anyway - appreciate the help but I got the problem fixed.