53
Something important to keep in mind is that Mathematica parses x / y as
Times[x, Power[y, -1]]
For actual floating point division, use Divide:
Divide[275., 6.]*1.03692775514337 // InputForm
(* 47.52585544407113 *)
which should agree with the C++ result.
38
Oleksandr is correct about the way evaluation works. a/b seems to be interpreted (parsed) directly as Times[a, Power[b,-1]], or more readably: $a\times b^{-1}$. Divide[a,b] is interpreted as is. Evaluation then proceeds from these forms, and the arithmetic is carried out differently for the two cases: either $a\times (1/b)$ or $a/b$.
Here are some ...
32
The backtick is a short-hand to mark the precision of your output. If it is not followed by any number, it denotes machine precision. You can denote arbitrary precision by including a number, as for example, 0.3`20.
By default, these are not displayed in StandardForm, which is why you see them only when copying, at which point it gets converted to InputForm....
30
The default value of
$NumberMarks
Automatic
means that ` should by default be used in arbitrary-precision but not machine-precision numbers. Arbitrary-precision numbers can contain an arbitrary number of digits e.g. :
Sqrt[3`21] == 1.73205080756887729353
Machine numbers contain the same number of digits and maintain no information on their precision, ...
28
What is wrong: a) you're using exact arithmetic. b) You keep iterating even if the point seems to be escaping.
Try this
ClearAll@prodOrb;
prodOrb[c_, maxIters_: 100, escapeRadius_: 1] :=
NestWhileList[#^2 + c &,
0.,
Abs[#] < escapeRadius &,
1,
maxIters
]
prodOrb[0. + 10. I]
prodOrb[0. + .1 I]
(if you don't need the entire list but ...
25
Normally Plot uses machine precision numbers; your $x^x$ expression is hitting the limit of the numbers that can be represented in machine precision right about $x>143$.
Note:
Solve[$MaxMachineNumber == x^x, x]
(* Out: {{x -> 143.016}} *)
You can increase the WorkingPrecision setting for Plot adequately, and the plot will be complete:
f[x_] = Sin[...
22
SetPrecision[] does this:
SetPrecision[0.1, ∞]
3602879701896397/36028797018963968
22
The result given by python is completely wrong, as are the results given by Mathematica for machine numbers. To get a correct result, you need to use extended precision numbers:
N[Sin[10^50], 10]
Precision @ %
-0.7896724934
10.
This result is correct to 10 digits. In order to get a correct result, Mathematica needs to use extra precision beyond the ...
21
There are two undocumented functions which control the tolerance of Equal and SameQ: Internal`$EqualTolerance and Internal`$SameQTolerance. I refer you to this answer of mine for further information.
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 ...
18
As the comments indicate, there is no completely hardware-based solution - but that doesn't mean you can't do some tweaking. The trick is always: stick with machine precision as long as you can, then switch to arbitrary precision only to refine your results.
Instead of making up an example (which is hard because Mathematica implements the above principle ...
18
Congratulations! You find one of subnormal positive double :)
Another example
f = Compile[{{t, _Real}}, 2.0^t];
f[-1074]
f[-1075]
5.*10^-324
0
MachineNumberQ@f[-1074]
True
This doesn't mean that CompiledFunction can work with arbitrary-precision numbers.
Update
Normally Mathematica prevents such numbers
2.^-1074.
% // MachineNumberQ
4.940656458413*...
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 ...
18
There are several important things about the way computer systems represent real numbers, which most of the time can be blithely ignored, just like the safety of bridges in the United States.
One important thing is that numbers are discrete. With regular machine precision (double precision), the mantissa has 53 bits, which provides a lot of resolution. ...
17
If you have an analytic formula for f[x_] := Erfc[x]*Exp[x^2] not using Erfc[x] you could do what you expect. However it is somewhat problematic to do in this form because Erfc[x] < $MinNumber for x == 27300.
$MinNumber
1.887662394852454*10^-323228468
N[Erfc[27280.], 20]
5.680044213569341*10^-323201264
Edit
A very good approximation of your ...
17
While it is a bad idea to compare floating point numbers, in this case I think something simpler is going on: you have an off-by-1 problem in your loop. See what your code does:
test[t_, dt_] := Module[{},
For[ti = dt, ti <= t, ti = ti + dt, Print[ti];];
Print[MemberQ[{0.01, 0.02}, ti]];
Return[0];
]
then, after
test[0.01, .001];
look at ti:
...
17
Just a bit of fun with @acl's code:
ArrayPlot[Table[
NestWhile[#^2 - (0. - 1 I) & , r + i I, Abs[#] < 2.0 &, 1, 10],
{r, -2, 2, 0.005},
{i, -2, 2, 0.005}]]
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 ...
16
Compilation of Total
As pointed out in a comment, it seems that the compensated summation form of Total can't be compiled. You can check this using CompilePrint - note the call to MainEvaluate.
Needs["CompiledFunctionTools`"]
CompilePrint@f2 (* from your question *)
Summation in Mathematica
There seem to be plenty of options for summing a list in ...
16
Without code and your actual results, this question cannot be answered. Here is one thing that might help: We have a compiler that can compile to C it can show you the code it creates. So why don't you try this?
a = 275.;
b = 6.;
c = 1.03692775514337
fC = Compile[{{a, _Real}, {b, _Real}, {c, _Real}},
a/b*c,
CompilationTarget -> "C"
]
Now we can ...
15
I've had to work with that kind of function (relying on cancellation of large terms) before, and the most practical workaround I could figure out to be able to evaluate the function numerically is to use its power expansion near the point of trouble (here, $+\infty$).
So, get a good look at the series expansion and find out how it works (or derive it on ...
15
I wonder whether I have understood your question correctly because I know you'll be aware of Clip
data =
Clip[#, {-$MaxMachineNumber, $MaxMachineNumber}] & /@ {0, Exp[1000.]}
(*
==> {0, 1.797693135*10^308}
*)
Precision /@ data
(*
==> {\[Infinity], MachinePrecision}
*)
data = RandomReal[10, {10, 2}]~Join~{{0, Exp[1000.]}};
Graphics[Point[...
answered Mar 14 '12 at 9:07
Sjoerd C. de Vries
61.9k12 gold badges168 silver badges305 bronze badges
14
For numerical evaluation, there is the rapidly-converging continued fraction (due to Jones and Thron):
$$\exp(x^2)\mathrm{erfc}(x)=\frac{2x}{\sqrt \pi}\cfrac{1}{2x^2+1-\cfrac{1\cdot2}{2x^2+5-\cfrac{3\cdot4}{2x^2+9-\cdots}}},\qquad x > 0$$
One can use the built-in function ContinuedFractionK[] with a suitable cut-off:
With[{x = N[30000], n = 10}, -2 x /(...
answered Apr 17 '12 at 13:34
14
You can define your own "precise equal" using Congruent (≡) (entered as Esc===Esc or \[Congruent]):
Congruent[x_, y_] := Equal @@ SetPrecision[{x, y}, Infinity]
Now
x ≡ y
False
for your example.
13
With a compiled version you get it so fast, that you can manipulate it in real time.
fc = Compile[{{in, _Complex, 0}, {c, _Complex, 0}},
Module[{iter = 0, max = 10, z = in},
While[iter++ < max,
If[Abs[z = z^2 + c] > 2.0,
Break[]
]
];
{Abs[z], iter}
], CompilationTarget -> "C", Parallelization -> True,
...
13
It seems I found my answer in OleksandrR's comment to this question. He says,
Bear in mind Equal applies an extra tolerance in Mathematica. The proper comparison is
Block[{Internal`$EqualTolerance = -Infinity}, 1 == 1 + $MachineEpsilon] (* False *)
Block[{Internal`$EqualTolerance = -Infinity}, 1 == 1 + $MachineEpsilon/2] (* True *)
In fact, the value ...
13
Numerics in Mathematica can be as precise as you like. However, precision comes at price; you pay for it in computation time and in additional coding effort.
In Mathematica there are several computational classes of non-complex numbers, which form a tree like this.
The computation you made was made with machine reals because you included 0.5 as a term. ...
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 ...
12
You would do well to understand the difference between tools that are intended for structural operations and those that are intended for mathematical operations. DeleteDuplicates is of the former, generally speaking. As such it is comparing the exact FullForm of the objects, or at least something close (caveat).
As b.gatessucks recommends in a comment you ...
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 ...
Only top voted, non community-wiki answers of a minimum length are eligible
