Hi Nathalie (in English for everybody to be able to read it…)

Well, I use an integer random number generator and I cast the result to floating point afterwards… Anyway I guess all simple random number generators use integer numbers inside (integer divisions).

@Nathalie

That the lattice of representable numbers is finer near 0 than near 1 does not mean that a random generator for floating point numbers will selects numbers near 0 more often than numbers near 1. If it is a good generator, it will not fall into that trap.

@Oliver

The problem under consideration here is independent on the internal employment of integers. It's internally binary anyway.

@Olivier@Ulrich: This topic does not deal with random number generator, but with the use of real numbers in computing.

Well then… now that the question has changed I must say I don't understand it anymore ;-) I understand your point about floats being "more dense" near 0, but what exactly is the problem now ?

The problem is the accuracy of computational results. The systematic errors are not a white noise (with zero average) but encounter a bias depending on how close to zero the calculation take place which can dramatically change the computed results. It is especially a problem in non linear dynamics.

Nathalie,

it sounds as if you just discovered the phenomenon of numerical errors in numerical calculations. Is there any specific, non-commonplace, insight that you want to let us know?

Yes, the absolute errors scale with the size of the quantity of interest; however the relative error is constant over the entire floating point range. The floating point specifications (IEEE 754-2008, and ISO/IEC/IEEE 60559:2011) provide more detail than most humans would care to know. Goldbergs "What Every Computer Scientist Should Know About Floating-Point Arithmetic" [ACM Computing Surveys, vol 23, no 1, March 1991] is a good little read.

If your computations are extremely sensitive to finite precision errors, you can try to compute in extended precision: libquadmath from http://gcc.gnu.org implements __float128; other sources for extended/arbitrary precision packages can be found at http://crd-legacy.lbl.gov/~dhbailey/mpdist/ and https://gmplib.org/

@Peter, of course you can always use extended precision but it increases the computing time. The problem that I have is: I have something which diverges like 1/x: if I compute it near zero the accuracy is good but if I compute it near 1 i.e.1/(1-x) the numerical errors increase...do you have a solution to this problem?

@Natalie

As a rule, numerical integration over singularities is to be avoided. Often mathematical reformulations (e.g. combining contributions from both sides of the singularity) let singularities disappear.

Doing the singularity 'by brute force' is the last resort. If problems arise with the presently employed floating point numbers one has to go to longer numbers.

For some compilers (e.g. g++) one replaces (employing macros) type double by type long double and thus 64-bit real numbers by 80-byte real numbers. Execution time will increase by, say, 30% and hardly will become prohibitive. One gets a clear indication of the numerical error of the final result.

@Ulrich, I am not doing an integration over a singularity, How do you get a clear indication of the numerical error (not doing an integration, I deal with non linear dynamics)?

@Nathalie

The trivial, in most cases sufficient, test is to do the computation with 'double' and (in addition) with 'long double' numbers and find that the results differ only by empirically irrelevant values. In non-linear dynamical systems the exact trajetories are not empirically relevant and nevertheless one might want to compute them. Then one needs ever longer numbers to follow the exact trajectories a bit longer. An example of this type is worked out in the appended work.

Article Precision-dependent symmetry breaking in simulated motion of...

@Ulrich, thanks for your article, very helpful!

I am not quite sure where your problem is or why it is a problem in your case. But maybe my insights might still help.

The idea behind IEEE floating point numbers is that usually you will work with numbers in a certain range, i.e. they all have similar exponents. In these cases floating point numbers are always sufficient (since numerical methods are not exact on their own why bother with higher precision in a slightly wrong result?). Double precision floating point provides 16 decimal digits of accuracy which is sufficient in most applications. This also shows that around x=10e-17 the addition 1+x starts to break.

I once had the problem that I summed a lot of small numbers which eventually accumulated to a significantly larger number (at least they should have). Then I found out about the so-called Kahan summation algorithm (http://en.wikipedia.org/wiki/Kahan_summation_algorithm). This algorithm uses two floating point numbers, one with the significant part and one with the accumulated error. I wrote myself a class to use in C++. But, I only needed summation. Extending this approach to multiplication and division might be possible, but is not straight forward.

Extended precision will not help in your case as it will only represent up to 19 decimal digits. If you definitely need 30 digits of precision you will have to use the according number of bits. According to Wikipedia quad precision (i.e. 128 bit floating point numbers) will give you 34 decimal digits of precision. So, it is your choice between speed and accuracy. You cannot have both.

Here is another thought: If your range is clearly [-1,1] you can think about using a fixed point integer format. For the smallest representable number being 10^-30 you will need log(10^30)/log(2) ~ 99.6 bits. Add a sign to have positive and negative numbers and you will get 101 bit integers to represent fixed point numbers. In the end fixed point numbers might be faster than quad precision floating point operations. But I don't have any experiences on my own in this area.

It is hard to understand that among the standard functions for floating point numbers there is no shift (or step) function which takes a floating point argument x and an integer argument n and returns the |n|-th discrete neighbor of x in direction sign(n). In my own C++ tool set (named C+- , see www.ulrichmutze.de) there is such a function and a function that gives the distance between real numbers in just representable steps . It is very instructive to see that evaluation of various mathematically equivalent expressions gives results which differ by a few or by thousands of discrete steps.

Fortran has some functions in this area:

nearest (x, s) returns the nearest machine number different from X in the direction of the sign of s.

spacing(x) returns the spacing between x and its closest neighbor.

rrspacing(x) returns the reciprocal of the relative spacing of values between x and its neighbor.

There are also several IEEE functions supported, including IEEE_next_after(x,s) which returns the same result as nearest(x,s) if the processor uses IEEE floating point representation.