31
This is indeed called catastrophic cancellation. In fact, this particular case is very easy: rewrite the function using the equivalent, numerically stable expression $$ \frac{t}{1+\sqrt{1-t^2}}. $$
Since you probably need a reference, this is discussed in most numerical methods textbooks in relation to the formula for solving quadratic equations (that ...
26
Are these integers or floating point numbers? Assuming it's floating point, I would go with the first option. It's better to add the smaller numbers to each other, then add the bigger numbers later. With the second option, you'll end up adding a small number to a big number as i increases, which can lead to problems. Here's a good resource on floating point ...
25
The best solution that I know of is to program the symbolic expressions in Mathematica, Maple, or SymPy; all of the links go directly to the code generation documentation. All of the programs above can generate code in C or Fortran.
None of the programs above mentions accuracy in IEEE 754 arithmetic; in general, it would be difficult to anticipate all ...
25
We all know that
\begin{equation}
\exp(x) = \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + x + \frac12 x^2 + \dots
\end{equation}
implies that for $|x| \ll 1$, we have $\exp(x) \approx 1 + x$. This means that if we have to evaluate in floating point $\exp(x) -1$, for $|x| \ll 1$ catastrophic cancellation can occur.
This can be easily demonstrated in python:
>&...
24
animal_magic's answer is correct that you should add the numbers from smallest to largest, however I want to give an example to show why.
Assume we are working in a floating point format that gives us a staggering 3 digits of accuracy. Now we want to add ten numbers:
[1000, 1, 1, 1, 1, 1, 1, 1, 1, 1]
Of course the exact answer is 1009, but we can't get ...
20
Another way to do this is to first explicitly specify the precision you desire in the variable using the SELECTED_REAL_KIND intrinsic and then use this to define and initialize the variables. Something like:
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
REAL(dp) :: x
x = 1.0_dp
A nice advantage to doing it this way is that you can store the definition ...
18
Adding arbitrary floating point numbers will usually give some rounding error, and the rounding error will be proportional to the size of the result. If you calculate a single sum and start by adding the largest numbers first, the average result will be larger. So you would start adding with the smallest numbers.
But you get better result (and it runs ...
18
I think Higham's Accuracy and Stability of Numerical Algorithms addresses how one can analyze these types of problems. See Chapter 2, especially exercise 2.8.
In this answer I'd like to point out something that isn't really addressed in Higham's book (it doesn't seem to be very widely known, for that matter). If you are interested in proving properties of ...
17
If you're looking for a good bound on your rounding error, you don't necessarily need an aribtrary-precision library. You can use running error analysis instead.
I wasn't able to find a good online reference, but it's all described in Section 3.3 of Nick Higham's book "Accuracy and Stability of Numerical Algorithms". The idea is quite simple:
Re-factor ...
17
Double precision is fairly common on newer GPUs. For instance I own a NVIDIA GTX560 Ti (fairly low end when it comes to computing) that has no issue running ViennaCL in double precision. From here (section 4) it appears all NVIDIA cards from GTX4xx onward support double precision natively.
I would guess that the GROMACS information is simply outdated.
15
There is a straightforward solution with only two passes through the data:
First compute
$$K := \max_i\; a_i,$$
which tells you that, if there are $n$ terms, then
$$\sum_i e^{a_i} \le n e^K.$$
Since you presumably don't have $n$ anywhere near as large as even $10^{20}$, you should have no worry about overflowing in the computation of
$$\tau := \sum_i e^{...
15
For the general case, I'd use compensated summation (or Kahan summation). Unless the numbers are already sorted, sorting them will be much more expensive than adding them. Compensated summation is also more accurate than sorted summation or naive summation (see the previous link).
As for references, What every programmer should know about floating-point ...
15
This question is very closely related to (and possibly a duplicate of) Is variable scaling essential when solving some PDE problems numerically?.
There are still good practical reasons to nondimensionalize equations, if possible:
It reduces the number of independent parameters for parametric studies (which was one of the original reasons for ...
15
No, that is not guaranteed. If you are using a NETLIB BLAS without any optimizations, it it mostly true that the results are the same. But for any practical usage of BLAS and LAPACK one uses a highly optimized an parallel BLAS. The parallelization causes, even if it only works in parallel inside the vector registers of a CPU, that the order how the single ...
14
If you want an idea of just how far we are away from such a software package, please look at the 2001 LAPACK working note on computing Givens rotations reliably and efficiently. I would expect most non-specialists (and many specialists!) in numerical analysis to be surprised at just how much analysis went into solving such an ostensibly simple problem:
...
14
Here is R1, as computed in MATLAB:
1.0e+07 *
-7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830
-1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472
3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415
-5....
14
First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :)
Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...
12
You declared the variables as double precision, but you initialized them with single precision values.
You could have written:
X=1.0d0
Y=1.0d-1
Barron's answer below is another way of making a literal double precision, with the advantage that it allows you to change the precision of your variables at a later time.
12
The formula
$$\mathrm{logsum}(x,y)=\max(x,y)+\mathrm{log1p}(\exp(-\operatorname{abs}(x-y))$$
should be numerically stable. It generalizes to a numerically stable computation of
$$\log \sum_i e^{x_i} = \xi+ \log\sum_i e^{x_i-\xi},~~~\xi=\max_i x_i$$
In case the logsum is very close to zero and you want high relative accuracy, you can probably use
$$\...
12
Use the (IEEE standard) library function log1p, which should be present in all programming languages. The function log1p(x) returns $\log(1+x)$, and is implemented with particular attention to accuracy when $x$ is small. It is designed to solve exactly this kind of problem.
11
Code generation and precompilation of mathematical expressions is becoming more popular.
While symbolic packages like SymPy, Mathematica, and Maple may include code generation I'm not confident that any of them also think hard about numerics.
There are a couple of other projects one could look into that are interested both in symbolics and in numerics.
...
11
You can directly check for denormals using isnormal() from math.h (C99 or later, POSIX.1 or later). In Fortran, if the module ieee_arithmetic is available, you can use ieee_is_normal(). To be more precise about fuzzy equality, you have to consider the floating point representation of denormals and decide what you mean for the results to be good enough.
More ...
11
$p \log p$ won't suffer from precision loss anywhere in $[0,1]$, and won't suffer from exponent overflow near $0$ either. Thus, the fast, accurate way is
p ? p * -log(p) : 0
11
In the 1980's era Intel 80x86 architecture, there was a scalar floating point unit that had instructions like FSIN, FCOS, etc. for computing functions like sin and cos. These functions were implemented in microcode and might take 100's of CPU cycles to execute.
Later, Intel added Streaming SIMD Extensions (SSE) which gave the processor parallel floating ...
10
To keep precision while you add doubles together you need to use Kahan Summation, this is the software equivalent to having a carry register.
This is fine for most values, but if you are getting overflow then you are hitting the limits of IEEE 754 double-precision which would be about $e^{709.783}$. At this point you need a new representation. You can ...
10
It sounds like you want a way to evaluate how FPU-bound your code is, or how effectively you are using the FPU, rather than to count the number of flops according to same anachronistic definition of a "flop". In other words, you want a metric that reaches the same peak if every floating point unit is running at full capacity every cycle. Let's look at an ...
10
If you're looking for a good introduction, I would suggest David Goldberg's What Every Computer Scientist Should Know About Floating Point Arithmetic. It may be a bit too detailed, but it's available online for free.
If you've got a good library, I'd suggest Michael Overton's Numerical Computing with IEEE Floating Point Arithmetic, or the first few chapters ...
10
Some references on rounding error analysis of Krylov methods:
Gutknecht, Martin H., and Zdenvek Strakos. "Accuracy of two three-term and three two-term recurrences for Krylov space solvers." SIAM Journal on Matrix Analysis and Applications 22.1 (2000): 213-229.
Paige, Christopher C., and Zdenvek Strakos. "Residual and backward error bounds in minimum ...
10
The point they are making is that if you compute the difference
123457.1467 - 123456.659 = 0.4877
then you get an answer that seems to be quite different than the difference
123457.1 - 123456.7 = 0.4
that you get by subtracting the rounded numbers. The ~20% figure is because 0.4877 is more than 20% larger than 0.4.
This example shows that if you subtract ...
10
The numeric precision is not perfect. You get rounding errors during your computation. When working with floats, don't check if they are = 0, but check if their absolute distance to 0 is smaller than some epsilon.
Only top voted, non community-wiki answers of a minimum length are eligible
