20
In typical floating point implementations, the result of a single operation is produced as if the operation was performed with infinite precision, and then rounded to the nearest floating-point number.
Compare $a+b$ and $b+a$: The result of each operation performed with infinite precision is the same, therefore these identical infinite precision results ...
14
The usual trick to avoid this underflow is to compute with logarithms, using the identity $$ \log \prod_{i=1}^n p_i = \sum_{i=1}^n \log p_i. $$
That is, instead of using probabilities, you use their logarithms. Instead of multiplying them, you add them.
Another approach, which is not so common, is to normalize the product manually. Instead of keeping just ...
12
A LU decomposition of a $n \times n$ matrix can be computed in $O(M(n))$ time, where $M(n)$ is the time to multiply two $n \times n$ matrices. Therefore, you can find a solution to a system of $n$ linear equations in $n$ unknowns in $O(M(n))$ time. For instance, Strassen's algorithm achieves $M(n) = O(n^{2.8})$, which is faster than Gaussian elimination. ...
10
No, $O(n \lg q)$ running time is not achievable. It takes $\Omega(q)$ space even just to write out the answer, so any algorithm will necessarily have running time at least $\Omega(q)$.
However, you can find algorithms that are more efficient than the naive solution. The naive solution for evaluating a polynomial of degree $n$ at $q$ points takes $O(nq)$ ...
9
If you want to multiply two matrices $A$ and $B$ then observe that
$$\begin{pmatrix}I_n&A&\\&I_n&B\\&&I_n\end{pmatrix}^{-1}=
\begin{pmatrix}I_n&-A&AB\\&I_n&-B\\&&I_n\end{pmatrix}$$
which gives you $AB$ in the top-right block.
It follows that inversion is at least hard as multiplication.
EDIT: I had misread the ...
8
One of the problems with Newton's method is that it requires a division operation in each iteration, which is the slowest basic integer operation.
Newton's method for the reciprocal square root, however, doesn't. If $x$ is the number for which you want to find $\frac{1}{\sqrt x}$, iterate:
$$r_{i+1} = \frac{1}{2} r_i (3 - x r_i^2)$$
This is often ...
8
There actually is some research on improving the numerical stability of floating point expressions, the Herbie project. Herbie is a tool to automatically improve the accuracy of floating point expressions. It's not quite comprehensive, but it will find a lot of accuracy improving transformations automatically.
Cheers,
Alex Sanchez-Stern
8
The wikipedia article on Methods of computing square roots: base 2 presents a strikingly similar snippet of C-code [a], but the link to the source is dead. Let's try to do better.
The snippets from both wikipedia and the question are very similar to Martin Guy's widely circulated C implementation [b]. which contains a comment: From a book on programming ...
7
No. You could make a decent argument that every algorithm is in some sense a symbolic manipulation, though not in the sense normally used for symbolic computation - they are not explicitly algebraic, and more naturally find expression as logic and set theory (unsurprisingly). Taking the wider view of symbolic manipulation though, even numerical analysis is ...
6
There's nothing fundamentally hard about computing $\sin(10^{99})$. You simply compute $x = 10^{99} \bmod 2\pi$, then compute $\sin(x)$. (Why is this valid? It's because $\sin(x)=\sin(y)$ if $x\equiv y \pmod{2\pi}$.) It's not too hard to compute $x$ if you use a numerical representation that has enough digits of precision, and then to compute $\sin(x)$ ...
6
A quadratic program is an optimization problem where the goal is to minimize $y^T Q y + c^T y$ subject to $A y \leq b$. If $Q$ is positive definite, then this is a convex quadratic program and we can solve this problem in polynomial time using several methods, one being the ellipsoid method (originally due to Kozlov, Tarasov and Khachiyan [1]).
We can ...
6
Consider the case that x is small. (1 + x) has a rounding error; the result that you get is not (1 + x) but (1 + x') for some x' close to x. If x is very small, the relative difference between x' and x can be quite large. Trying to calculate log (1 + x) will calculate log (1 + x') which can have a large relative error. Instead the better formula calculates ...
6
It sounds extremely inefficient compared to floating point. We have a very good understanding of how to control the errors in floating-point calculations (e.g., adding small numbers before large ones, avoiding taking the difference of large numbers and so on) so the only benefit of the methods you're suggesting would appear to be that they offer increased ...
6
A "numerically stable" method calculates a result in a way that will not produce excessive rounding errors.
Given a small number x, let's say x = 0.00079, it is possible to calculate log (1 + x) with much higher precision than 1 + x. (Of course calculating log (1 + x) avoids adding 1 + x). In both cases the relative error will be small and about equal size,...
6
With some algebraic manipulation (as pointed out in @orlp's answer), we can deduce the following:
$$f(x) = x \tanh(\log(1+e^x)) \tag{1}$$
$$ = x\frac{(1+e^x)^2 - 1}{(1+e^x)^2 + 1} = x\frac{e^{2x} + 2e^x}{e^{2x} + 2e^x + 2}\tag{2}$$
$$ = x - \frac{2x}{(1 + e^x)^2 + 1} \tag{3}$$
Expression $(3)$ works great when $x$ is negative with very little loss of ...
5
Let me illustrate one problem which could happen, and one way to solve it. You want to distribute a given amount of cents $N$ into $k$ piles, in proportions $p_1,\ldots,p_k$, where $p_1,\ldots,p_k \geq 0$ and $p_1 + \cdots + p_k = 1$. The problem is that $Np_i$ need not be an integer. As a simple example, you might want to divide $20$ cents into $1/3:2/3$ ...
5
The answer will depend on the compiler. As @vonbrand wrote, "Given a good enough compiler, you might even get the very same object code." In particular, good compilers will do tail-call elimination. In some cases this can effectively transform the code into a for-loop. Your example looks like a good example of an instance where this could happen.
As @...
5
The Bailey–Borwein–Plouffe formula only works in hexadecimal. There might be other formulas for other bases, but I'm not aware of a decimal-based formula. If you want to obtain the $N$th decimal digit, you have to compute enough hexadecimal digits, there are no shortcuts.
5
You can use Newton's method or any of a number of other methods for finding approximations to roots of the polynomial $p(x) = x^2 -c$.
The rate of convergence for Newton's method will be quadratic, which means that the number of bits that are correct doubles in each iteration. This means $O(\lg n)$ iterations of Newton's method suffice.
Each iteration of ...
5
A paper Computing the Entropy of a Large Matrix by Thomas P. Wihler, Bänz Bessire, André Stefanov suggests approximating $x \lg x$ with a polynomial. Then you can use the trace of powers of the matrix to sum the results of applying that polynomial to each of the eigenvalues.
(The polynomial-via-power-and-trace thing works because density matrices have ...
5
"Numerical stability" is a much vaguer term than most people realise. We typically use it when referring to an approximation method, such as some kind of linear analysis, or numeric quadrature, or solving (possibly partial) differential equations. It refers to the property that given some appropriate assumptions (e.g. the inputs are reasonable), the ...
5
If you have multiple processors that can work in parallel, then you can calculate any power up to the power (2^k) in k steps. For example: To calculate $M^{15}$, you calculate:
Stage 1: Calculate $M^2$
Stage 2: Calculate $M^3 = M^2 * M$ and $M^4 = M^2 * M^2$
Stage 3: Calculate $M^7 = M^4 * M^3$ and $M^8 = M^4 * M^4$
Stage 4: Calculate $M^{15} = M^8 * M^...
5
Java uses IEEE 754 binary floating point representation, which dedicates 23 binary digits to the mantissa, that is normalized to begin with the first significant digit (omitted, to save space).
$0.00004_{10} = 0.00000000000000101001111100010110101100010001110001101101000111..._{2} = [1.]\color{red}{01001111100010110101100}010001110001101101000111..._{2} \...
5
The binary floating point format supported by computers is essentially similar to decimal scientific notation used by humans.
A floating-point number consists of a sign, mantissa (fixed width), and exponent (fixed width), like this:
+/- 1.0101010101 × 2^12345
sign ^mantissa^ ^exp^
Regular scientific notation has a similar format:
+/- 1.23456 × 10^...
5
There is what you want to achieve, and there is reality, and sometimes they are in conflict. First you check if your problem is a special case that can be solved quicker, for example a sparse matrix. Then you look for faster algorithms; LU decomposition will end up a bit faster. Then you investigate what Strassen can do for you (which is not very much; it ...
4
Taking the sine of large numbers is a numerically unstable operation.
Considering an argument like $10^{99}$, you can get a completely different value of the sine by adding, say $1$ to it. Think that this is a relative change of $10^{-99}$ !
Indeed,
$$|\sin(a+1)-\sin(a)|=|2\sin(\frac12)\cos(a+\frac12)|>0.95|\sin(a+\frac12)|,$$
so that you can find ...
4
This is basically an instance of the line segment intersection problem. One standard approach is to use a sweep line algorithm. For instance, the Bentley-Ottman algorithm would be a reasonable choice, and is not too difficult to implement.
At each iteration of the algorithm, we have some value of $x$, and we calculate what is the next largest value of $x$ ...
4
The partial convergents of the continued fraction of $x$ consists of all the best rational approximations of $x$; see Wikipedia, for example. A best rational approximation of $x$ is a rational number $p/q$ such that $\left|x-\frac{p}{q}\right| \leq \left|x-\frac{p'}{q'}\right|$ for all $q' \leq q$. Your $p/q$ is in particular a best rational approximation, ...
4
The relation between the Stern–Brocot tree and Farey sequences shows that if $0 < p/q < 1$ and $(p,q) = 1$ (i.e., $p/q$ is a reduced fraction) then $p/q$ is at the $q$th level of the tree. Since the running term of your algorithm is linear at the level in which you terminate, your algorithm takes time $O(q)$, where $p/q$ is the answer; but this is not ...
4
The phrase "machine epsilon" is misleading. In reality, it doesn't depend on the machine at all, but on the floating-point data types that the programmer chooses to use.
Only top voted, non community-wiki answers of a minimum length are eligible
