29
votes
Accepted
Is half precision supported by modern architecture?
Intel support for IEEE float16 storage format
Intel supports IEEE half as a storage type in processors since Ivy Bridge (2013). Storage type means you can get a memory/cache capacity/bandwidth ...
- 2,136
28
votes
Accepted
Numerically stable way of computing angles between vectors
(I have tested this approach before, and I remember it worked correctly, but I haven't tested it specifically for this question.)
As far as I can tell, both $\|\mathbf{v}_1\times \mathbf{v}_2\|$ and $...
- 11.5k
26
votes
Accepted
How can I avoid catastrophic cancellation?
Take
\begin{align}
1-\sqrt{ 1-x^2} &= (1-\sqrt{ 1-x^2})\frac{1+\sqrt{ 1-x^2}}{1+\sqrt{ 1-x^2}}\\
&= \frac{x^2}{1+\sqrt{ 1-x^2}}
\end{align}
So
\begin{align}
y = x\sqrt{\frac{1}{2+2\sqrt{1-x^2}}...
- 725
15
votes
Numerically stable way of computing angles between vectors
The efficient answer to this question is, not too surprisingly, in another note by Velvel Kahan:
$$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\...
- 3,165
13
votes
Accepted
Computing $\frac{x - y}{x - z}$ when $x,y,z$ are close to each other
If your inputs are $x,y,z$, this computation is not unstable, but ill-conditioned. That's worse, because it means that a small change in your input (such as a previous approximation as a floating-...
- 12.6k
12
votes
Accepted
Evaluating $\log(\exp(x)+1)$ for negative $x$
Use the (IEEE standard) library function log1p, which should be present in all programming languages. The function log1p(x) ...
- 12.6k
10
votes
Numerical derivative and finite difference coefficients: any update of the Fornberg method?
Overview
Good question. There is a paper entitled "Improving the accuracy of the matrix differentiation method
for arbitrary collocation points" by R. Baltensperger. It's no big deal in my opinion, ...
- 3,237
10
votes
Accepted
Analytical convergent sequence and numerical divergent sequence
Jean-Michel Muller, et. al., "Handbook of Floating-Point Arithmetic 2nd ed.", Birkhäuser 2018, gives the following example due to Muller, specifically constructed to deliver incorrect results with ...
- 1,895
10
votes
Accepted
Accuracy loss in single-precision Euclidean norm computation
The squares are harmless (as long as they don't overflow/underflow), because the relative perturbation they introduce is of the order of the machine precision $u\approx 10^{-8}$. Your troubles here ...
- 12.6k
9
votes
Accepted
Is it possible to proof a-b+b = a for all double floating-point numbers?
You can sometimes prove such results (or get counterexamples) using an SMT solver such as Z3 that supports floating point arithmetic. Here is a proof of a version of your theorem that says $|((x+y)-y)-...
- 11.5k
8
votes
Accepted
Why can ill-conditioned linear systems be solved precisely?
Added after my initial answer:
It appears to me now that the author of the referenced paper is giving condition numbers (apparently 2-norm condition numbers but possibly infinity-norm condition ...
- 19k
8
votes
Is half precision supported by modern architecture?
In my opinion, not very uniformly. Low precision arithmetic seems to have gained some traction in machine learning, but there's varying definitions for what people mean by low precision. There's the ...
- 5,076
8
votes
Is half precision supported by modern architecture?
The accepted answer provides an overview. I'll add a few more details about support in NVIDIA processors. The support I'm describing here is 16 bit, IEEE 754 compliant, floating point arithmetic ...
- 181
8
votes
Unexpected result when summing sorted (and unsorted) positive floating point numbers
Very interesting problem! I might have a partial answer.
To start, I replicated a simple C++ demo that can reproduce the effect
...
- 996
6
votes
What good are hard-sphere event-driven molecular dynamics simulations in the face of chaos?
Yes, it is possible to show that the statistical behavior of the approximate system will reach that of the "exact" system. (This is true even though hard-sphere dynamics do not accurately describe ...
- 3,513
6
votes
Analytical convergent sequence and numerical divergent sequence
In my opinion, an example could be the calculation of the minimal solution of a three term recurrence relation (TTRR) $$y_{n+1} +a_{n}y_{n}+b_{n}y_{n-1} = 0, \quad n=1,2,3,\ldots$$.
For example, the ...
- 6,199
6
votes
What algorithm for solving a set of stiff ODEs would be easiest to port to high precision floating point arithmetic?
I prefer to always take a benchmark-informed approach to this. Look at the SciML Benchmarks:
https://benchmarks.sciml.ai/html/StiffODE/VanDerPol.html
https://benchmarks.sciml.ai/html/StiffODE/ROBER....
- 12.5k
5
votes
Accepted
How to deal with big numbers in intermediate calculations?
I think there's a simple way to do this. You have a rational function of identical cosh/sinh terms, where every expression is a homogeneous polynomial in cosh/sinh, and the only problem is that these ...
- 11.5k
5
votes
Numerical stability in the product of many matrices
Orthogonal matrices are about as well-conditioned as you can get, but numerical errors still occur. One common error is loss of orthogonality. A fix for this could be to re-orthogonalize your columns ...
- 3,101
5
votes
Numerically stable and fast sum of last K elements in sequence
Very interesting question!
LAPACK-inspired adaptive strategy
This reminds me of a bug that was found in a LAPACK routine (rank-revealing QR) related to 'downdating' norms: essentially, you are given ...
- 12.6k
5
votes
Accepted
Comparison of integrals with a function:
There is no need for numerical computation here.
First, $T(q)$ is a well-known function, the logarithmic integral. Repeated integration by parts gives an asymptotic expansion
$$\mathrm{Li}(q) = \...
- 166
4
votes
Evaluating $\log(\exp(x)+1)$ for negative $x$
If you want to use your own implementation to evaluate the function, take a look on this article. Otherwise, just go with log1p as suggested by @Frederico Poloni in his answer.
For more information ...
- 226
4
votes
Numerical derivative and finite difference coefficients: any update of the Fornberg method?
Assuming you are trying to differentiate a numerical implementation of a continuous function, there are a large number of methods:
1) Automatic differentiation. The most accurate and general method. ...
- 2,165
4
votes
Accepted
log-sum-exp trick for signed/complex numbers
After thinking about this some more, I can answer this one myself!
I don't think the complex plane makes the log-sum-exp trick appreciably different, at least in Cartesian coordinates. In particular,...
- 1,433
4
votes
Comparison of integrals with a function:
You write $S(q)$ and $T(q)$ as integrals, but it is easier to think of them as solutions of ODEs:
$$
S'(q) = \sin^2\left(\frac{π\Gamma(q)}{2q}\right)
$$
with initial conditions
$$
S(2) = 0,
$$
and ...
- 56.8k
4
votes
Is half precision supported by modern architecture?
You can find out if your hardware supports half-precision via:
...
- 2,165
4
votes
Accepted
What are the Exact Rules for Significant Figures, Precision, and Uncertainty?
The rules of significant figures are rule-of-thumb way to communicate errors and should only be seen as a primite first step to talk about uncertainties and measurement errors.
You gave the excellent ...
- 3,065
4
votes
Accepted
Robust unit test for reciprocal approximation
Compilers that support turning division into a multiplication with the reciprocal with some optimization flag settings (e.g. -ffast-math or ...
- 1,895
4
votes
Float equality tolerance for single and half precision
Those proposed tolerances look fine, but in my (opinionated) view this is really a problem with no satisfying solution, as the comments also argue.
For most algorithms, the error bounds one gets look ...
- 12.6k
3
votes
Half precision in Fortran
The Fortran standard doesn't specify what precisions and ranges of floating point numbers need be supported by a given compiler. All it says is that at least two different kinds of real numbers need ...
- 626
Only top scored, non community-wiki answers of a minimum length are eligible