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.


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 ...


9

There are aspects of modern computing systems that are inherently non-deterministic that can cause these kinds of differences. As long as the differences are very small in comparison with the required accuracy of your solutions, there probably isn't any reason to worry about this. An example of what can go wrong based on my own experience. Consider the ...


9

Their observation ''the form of the algorithm may needlessly introduce some numerical imprecision'' is correct. But their explanation ''This arises because, in eqn (3.14), a small term ($O(δt^2)$) is added to a difference of large terms ($O(δt^0)$), in order to generate the trajectory.'' is spurious. The true reason for the slight numerical instability of ...


6

There is a canonical answer to this question in GCC's wiki, which is presumably maintained, it is by far the most authoritative source of information for this type of question. This question, on the other hand, may eventually go out of date. This is all explained in more detail in the wiki, with examples. The following is essentially a quote from it to ...


6

Why not give GMPY2 a try? From the introduction: gmpy2 is a C-coded Python extension module that supports multiple-precision arithmetic. gmpy2 is the successor to the original gmpy module. The gmpy module only supported the GMP multiple-precision library. gmpy2 adds support for the MPFR (correctly rounded real floating-point arithmetic) and MPC (correctly ...


5

In double-precision floating-point arithmetic, the number $2^{-1022}$ is the smallest (normal) positive number that can be represented. Computing any smaller positive number results in a loss of precision, as it cannot be represented as a floating-point value. This either means the computed number will be denormal (partial loss of precision) or possibly it ...


5

I think, this one Krylov Subspace Methods in Finite Precision: A Unified Approach, Jens-Peter M. Zemke is also worth reading.


5

To apply Pedro's example to the equation $(3.14)$, assume your variables are stored with the following values: $$ \mathbf{r}(t) = 101 $$ $$ \mathbf{r}(t - \delta t) = 100 $$ $$ \delta t^2 \mathbf{a}(t) = 1.49 $$ From $(3.14)$ it should follow that $$ \mathbf{r}(t + \delta t) = 103.49 $$ but, since we can only use three digits, the result becomes ...


5

To first order approximation, icc does use -ffast-math by default. If you run icc --help, it will show you its options. One of the sections is Floating Point, which begins like this: Floating Point -------------- -fp-model <name> enable <name> floating point model variation [no-]except - enable/disable floating point ...


4

I think the main issue that makes this tricky is that trying to evaluate the power series for the hypergeometric function is numerically unstable due to catastrophic cancellation for one negative parameter and function value $\ll1$. Otherwise, having real parameters, and a real argument in $[0,1]$ should be fine. If you look at the DLMF entry for ...


4

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 about the log1p function read this other article from the same author.


3

The mentioned FFTW library might run in a non-deterministic mode. If you're using FFTW_MEASURE or FFTW_PATIENT mode, the programs checks at runtime, which parameter values work the quickest and then will use those parameters throughout the whole program. Because the run time will obviously fluctuate a bit, the parameters will be different and the result of ...


3

There are a number of symbolic algebra packages that will compute the exact matrix inverse; e.g. Sage, Mathematica, Maple, or Maxima. This becomes computationally expensive for large matrices, but there is a chance that a sparse $100 \times 100$ matrix would be within reach. I wouldn't rule it out without trying, since it's so easy to try.


3

You haven't told us the original source of this question, or precisely quoted it, so it's possible that something critical is missing from the statement of the problem. It's also quite possible that the solution is just wrong. However, a likely explanation of your confusion is in the idea of normalized binary floating point numbers with an implicit ...


3

Assume $m \ge n$. Consider the function $$ \begin{aligned}g(k) &= \frac{2^n-k-1}{2^n-k} \\ &= \frac{1-(k+1)2^{-n}}{1-k2^{-n}} \\ &= \left(1-(k+1)2^{-n} \right)\left(1 + k2^{-n} + k^2 2^{-2n} + O(2^{-3n})\right) \\ &= 1-2^{-n}-k 2^{-2n} + O(2^{-3n}) \end{aligned}$$ For large enough $n$, $g(1)$ and $g(2)$ differ by essentially $2^{-2n}$, and ...


3

Pedro already gives the important fact, namely cancellation. The point is that every number you compute with has an associated accuracy; for example, a single precision floating point number can only represent things up to approximately 8 digits of accuracy. If you have two numbers that are almost exactly the same but differ in the 7th digit, then the ...


2

Different vendors have different paths forward on this, and most of it is under NDA. If you really want to know, ask Intel, AMD, IBM, ARM, and maybe NVIDIA for briefings under NDA. Some of these vendors will have full support for vectorized integer instructions in the future so that integer will keep up with floating-point, and some of them may not. For ...


2

Let's suppose that by "fastest supercomputer in 2000", and as your metric of speed, you choose the High Performance LINPACK benchmark, as used in the TOP500 list. Depending on when "in 2000" means, you're either selecting as your metric of performance: 2.4 TFLOP/s, in June 2000 4.9 TFLOP/s, in November 2000 Setting aside all of the issues involved in using ...


2

The first thing to say is that such numbers exist. The "trivial" cases are where $n = F(n)$, namely $n=1,5$. After that the next solution is $n=12$, and many multiples of $12$ are also solutions. We will show below that infinitely many solutions exist. The second thing to note is that computing Fibonacci numbers $F(n)$ by the recurrence relation modulo $...


2

This may be more of a comment than an answer, but I can't comment. Yes on arbitrary precision. No on symbolic inverse, because at some point, you have to numerically evaluate it, and it is likely to be extremely numerically unstable. Even for 6 by 6 matrix, the symbolic inverse is starting to get pretty ungainly. Which leaves me the main item. Why do you "...


2

Parseval's theorem tells you that the dot product between two vectors equals the dot product of the Fourier transforms, possibly up to a constant. Consequently, you do not need to transform back the result of $a\ast b$ and $c \ast d$ and should be able to do the overall operation with just 4 FFTs.


2

Two's complement is always done with a specific bit width in mind, and you use modular arithmetic within that width. Basically, you should ignore that leading 1, it's the output carry and should be modulo'd away.


1

Alright I did some fixing and debugging. Quite a few mistakes I would say, I think you're a beginner in C/Programming ? Anyway find the code below... #include <math.h> #include <stdio.h> #define N 3 #define M 3 #define A 0.0 #define B 1.0 #define C 0.0 #define D 1.0 int main(int argc, char *argv[]){ double h_x = (double) 1/(N+1); double h_y = ...


1

While floating point rounding from async operations may be the issue, I suspect that it is something more banal. The use of uninitialized variable that is adding randomness to your otherwise deterministic code. It is a common issue that is often overlooked by developers because when you run in debug mode all variables are initialized to 0 on declaration. ...


1

While it's true that expression term evaluation order changes may very well occur due to multi-core/multi-thread processing scenarios, don't forget that there may be (even though it's a long shot) some kind of hardware design flaw at work. Remember the Pentium FDIV problem? (See https://en.wikipedia.org/wiki/Pentium_FDIV_bug ). Some time ago, I worked on pc-...


1

An $n^2$ algorithm can be reduced to something that behaves $O(n)$ if, for example, if $n$ is known to the compiler, and is a multiple of the vector size for the vector instructions (if any) supported by the processor. If the compiler can see all of this, it can unroll an inner loop and use vector instructions to do the work. This may reduce the overall ...


1

Firstly, I can't guarantee it, but I suspect that matlab will return a signed double rather than an unsigned, so you've probably only got 52 bits to play with. However to actually address the question, consider a clock that has stopped with the hour hand on 1 o'clock. It has a terrible precision, but if you look at it at 12.30, it's accurate to within half ...


1

You can use dec2bin(typecast((exp(1) + 1e12) - 1e12, 'uint64')) to obtain the binary representation of the double precision floating point. 100000000000101101111110000101010001011000101000101100000000000 100001001101101000110101001010010100010000000000101100000000000 100000000000101101111110000000000000000000000000000000000000000 ...


Only top voted, non community-wiki answers of a minimum length are eligible