68

Control the Precision and Accuracy of Numerical Results This is an excellent question. Of course everyone could claim highest accuracy for their product. To deal with this situation there exist benchmarks to test for accuracy. One such benchmark is from NIST. This specific benchmark deals with the accuracy of statistical software for instance. The NIST ...


24

As it seems to depend on more than machine bits I'm curious what $MaxNumber various Mathematica installs have. If your setup is different please fill in system information and Log2 @ Log2 @ $MaxNumber // Round in the table below: $$\begin{array}{r|c|c|l|c} \text{OS} & \text{Bits} & \text{Version} & \text{\\\$MaxNumber} & \log_2\log_2\\ \...


24

I have figured out why you are getting the structure you are getting. The reason has to do with your initial choice of the angle, which you set at $0$ in the Nest[] statement. The actual image is generated by choosing the mean result of iterating the map for many initial values chosen uniformly at random in $[0,1]$. With $n = 50$ iterations and $m = 20$ ...


20

David Goldberg ("What every computer scientist should know about floating-point arithmetic", ACM Computing Surveys, Vol 23, No 1, March 1991, p 12, Th 4) gives pseudocode that is equivalent to log1p[x_Real] := With[{w = 1 + x}, If[w - 1 == 0, x, x * Log @ w/(w - 1)]] EDIT - Following Mark Adler's comments, I checked the binary representation of the results ...


20

If you have a value $x$ with an absolute uncertainty $dx$ the precision of $x$ is by definition: $$\text{Precision}(x) = - \log_{10}(dx/x)$$ That is why for $x=0$ the precision is always infinity. You cannot change this by changing $dx$. If you want to assign an absolute uncertainty to zero, you can set accuracy, which is defined as $$\text{Accuracy}(x) ...


18

Update: yode points out in a comment below that there are log1p() and expm1() functions in Mathematica! However they are hidden. They are simply: Internal`Log1p Internal`Expm1 They operate only on numeric inputs. (I don't know at what version these showed up.) Compare this plot to the one further below using the Log function. The errors really are all ...


17

Actually this is not a duplicate. The prior question is about underflows that require massive bignums to represent at machine precision, and that much is present here as well. So what @J.M notes is certainly a part of the issue. But it is not entirely a matter of cancellation error and a need to represent a very small number. The problem is a combination ...


14

The OP linked to an answer of mine for interpolating over general point sets; for constructing a single interpolating function, a slight modification of my procedure is needed. (In particular, you don't need centripetal or chord-length parametrization in this case.) Let's start with some data: data = {{0, 0}, {1/10, 3/10}, {1/2, 3/5}, {1, -1/5}, {2, 3}, {3,...


14

The last column seems in error. Here's a workaround for the sample problem, although it does not fix NDSolve`FiniteDifferenceDerivative: dx2 = dx1; dx2[[All, -1]] = -Reverse@dx1[[All, 1]] (* {-0.50000000000000000000, 0.25989153247414500869, -0.29289321881345247560, 0.36161567304292239214, -0.50000000000000000000, 0.80995720221088751026, -1....


13

Here are my thoughts: Q1 Machine numbers: For machine numbers, what you describe is correct. I would just add that you can use InputForm or FullForm to see all the digits if desired: N[Pi] % //InputForm 3.14159 3.141592653589793 Extended precision numbers: For extended precision numbers (i.e., numbers whose precision is a number, and not ...


13

Sqrt is using exact methods in an effort to pull out "small" squares. This is going to take time. A direct approach, as already noted, would do the square root numerically. For purposes of Floor extraction, it suffices to use as precision half the digit size. floorSqrt[n_Integer] := Module[{prec = RealExponent[1.001*n]/2}, Floor[Sqrt[SetPrecision[n, ...


12

The determinant of the matrix A in this case is about 10^282 The determinant isn't very useful, but the condition number is: You can use SingularValuesList to get the largest and the smallest singular value. If the ratio between the two is too large, the matrix is ill-conditioned. Solving an ill-conditioned linear system will still give "exact" results (...


12

As others have mentioned, the wrong result is given by N[expr] and the errors are due to cancellation. Let's discuss a bit why N[expr, 3] is able to give a good result. Mathematica can do computations with inexact ( = floating point) numbers in two ways: Using the computers native floating point arithmetic, which is very fast, but has no precision ...


12

Here is my implementation of quadric spline arbitrary precision interpolation which is defined exactly as in Interpolation with options Method->"Spline", InterpolationOrder->2. Theoretical background Quadric spline interpolation for n datapoints is defined as a Piecewise function which consists of n-2 parabolas which are spliced together in the ...


12

Precision is the principal representation of numerical error Except for numbers that are equal to zero, error in arbitrary-precision numbers is stored internally as its precision. For numbers equal to zero, the accuracy is stored, because the precision turns out to be undefined (even if Precision[zero] is defined). One way to view zeros are as a form of ...


11

For the sake of alternatives, here's another way. You can specify a precision with a backtick. This is called an "arbitrary-precision" number. Such numbers are treated differently than machine precision numbers. Examples of 3 and 5 digits of precision: 16.`3 16.0 16.`5 16.000 Machine precision depends on the machine, but it's often nearly 16 ...


11

This question looks as a duplicate of these questions: How to create internally optimized expression for computing with high WorkingPrecision? How to work with Experimental`NumericalFunction? An internally optimized version of the original function can be created as follows: n = 500; f = Experimental`CreateNumericalFunction[{a, b}, Unevaluated[Nest[# ...


11

This is not an answer. But I don't believe we should close this question as "easily found in the documentation". Numerics in Mathematica is an extremely complicated and mostly undocumented subject, where several mathematical concepts run up against each other in subtle and non-trivial ways. I have been thinking for some time that we ought to address this ...


11

From SetPrecision: See also Precision and Accuracy.


11

Not directly addressing the question about GMP library, but you can get a Sqrt much faster starting with an extended precision float. n = 10^1000000 - 3^2095903; (g = Floor[Sqrt[n]]) // AbsoluteTiming // First 6.22126 (h = Floor[Sqrt[N[n, Log[10, n]]]]) // AbsoluteTiming // First 0.086576 g == h True using this you maybe want a validity check, (...


10

The problem with using N to generate arbitrary precision numbers is that is it will generate more digits of accuracy than asked for: N[Sqrt[3], 10] // InputForm (* 1.73205080756887729352744634150587236694`10. *) N[Sqrt[3], 28] // InputForm (* 1.73205080756887729352744634150587236694`28. *) Probably the square root algorithm converges quadratically, so ...


10

Consider any numerical integration method $I^*(f,a,b)$ that approximates the exact integral $I$ of a function $f$ over an interval $[a,b]$. It will be implemented by a computation represented by, say, I[f[x], {x, a, b}]. Conceiving an NIntegrate[] command in this way breaks the error into two independent components. The first is the error of the method $I^* -...


10

There are a few reasonable ways. I'll illustrate with an example of Newton iterations for square roots, take from this MathGroup post r[x_, n_] := x - (x^2 - n)/(2*x) x = 1.0`20; two = 2.0`20; First we run it with the usual arithmetic. Table[x = r[x, two], {30}] (* Out[680]= {1.4142135623730950488, 1.4142135623730950488, \ 1.414213562373095049, 1....


9

This looks like a duplicate of my question here Global precision setting . Since I formulated it in a slightly different way, coming from Maple it might however be worth quoting the answer given by Mr. Wizard here. It worked for me in Mathematica 9 if there are no arbitary precision numbers in the notebook (defined as e.g. 2.`16): $PreRead = (# /. ...


9

I think that is because 1. and 1.*^-15 are machine-precision numbers, and Mathematica does NOT do precision-tracking on machine-precision calculations. I suggust using the arbitrary-precision numbers instead. Through[{Precision, Accuracy}@#] & /@ {1., 1.*^-15} {{MachinePrecision, 15.9546}, {MachinePrecision, 30.9546}} Now we specify a precision to ...


9

The exact equality comparison returns unevaluated. Root[1 + 6 #1 - 12 #1^2 - 32 #1^3 + 16 #1^4 + 32 #1^5 &, 5] == root[1] (* Out[900]= Root[1 + 6 #1 - 12 #1^2 - 32 #1^3 + 16 #1^4 + 32 #1^5 &, 5] == Cos[(2 \[Pi])/11] *) The numerical values agree to all significant digits. N[Root[1 + 6 #1 - 12 #1^2 - 32 #1^3 + 16 #1^4 + 32 #1^5 &, 5], 20] (...


9

Incidentally, here is a higher resolution picture of the image produced by heropup's answer, plotted over both positive and negative values of the two circle map parameters: Blue represents negative, yellow represents positive. You can view a far larger 6001x6001 pixel image (about 29MB) here. Here is the code used to generate the upper right quarter of ...


8

Regarding your last question: in the docs for FindRoot it says that FindRoot continues until either of the goals specified by AccuracyGoal or PrecisionGoal is achieved. The same thing is mentioned in the docs for NMinimize. On the other hand, the docs for NDSolve say AccuracyGoal effectively specifies the absolute local error allowed at each step ...


8

Problem The problem with Log[1. + 1.*^-15] not yielding 1.*^-15 is not due to Log, but to MachinePrecision inputs, which I think the OP implied in the question statement: 1 + 1.*^-15 % - 1 (* 1. 1.11022*10^-15 *) So Log[1 + 1.*^-15] does return the right answer, 1.11022*10^-15, for the actual input. Solution Here is a simple way to get log1p-type ...


8

LogLogPlot[{Internal`Log1p[x], Log[1 + x]}, {x, 1*^-17, 1*^-14}] ps:Of course,maybe you need Internal`Expm1,too.


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