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.


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

This is more of a question appropriate to the CS forums, but here is the short of it: it requires you to know how the compiler parses numbers. In particular, when the compiler sees something like -2147483648, it parses this as unary minus applied to 2147483648. But the latter (a positive number) is too large for the int data type, and so may be wrapped ...


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

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


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

No it's not OK. There are no negative literals per-se in C. "-12345" is interpreted as the combination of the positive literal 12345 with the unary minus operator. 0x80000000 is outside the range of int and will be interpreted as a literal of type unsigned int. On most compilers you could get around this by using a typecast #define INT_MIN ((int)(...


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

Floating-point exceptions may help you here. C support varies by implementation (compiler) but see GCC here: https://www.gnu.org/software/libc/manual/html_node/FP-Exceptions.html Python support is documented here: https://docs.python.org/2/library/fpectl.html I’ve only used these features a few times, and then only with the Intel compiler (https://software....


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

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

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

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.


2

Normally one does not try to detect loss of precision algorithmically, but rather analyzes and modifies algorithms to assess how they are affected by it. For instance, in your first example you would run a (forward) error analysis and figure out that the summation error is bounded by $3 \cdot 10^{10} \mathsf{u}$, where $\mathsf{u}$ is machine precision, or ...


2

The idea is that if you have $n$ bits $b = [b_{n-1}, \ldots, b_0]$, then you can represent $2^n$ different things with them, as all $b_k \in \{0, 1\}$. You can decide to go for the numbers $0, \ldots, 2^n - 1$, by saying that let the bit sequence $b$ represent the number $\sum_{k=0}^{n-1} 2^k b_k$. If you need negative numbers, then you could say that the ...


1

OP's comment on OP: yes,numbers are intergers For arbitrarily large integer calculations, MAPLE is an option worth trying. it has very efficient implementations for a number of such calculations.


1

Using something such as numpy.divide, you most likely will see no considerable difference from doing it in C++. Although casting might be something to consider. C++ may see some improvement in this area as Python will have to deal with this in real-time. With C++ this will not pose an issue. (Python will now cast any ints to floats for doing division as ...


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


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