Why does this Quantity calculation lose precision?

Question

(Sqrt[#] Sqrt[Quantity[2, "ElectronMass" ("SpeedOfLight")^2] + #] )/(
  Quantity[1, "ElectronMass" ("SpeedOfLight")^2] + #) &[
 Quantity[50, "Gigaelectronvolts"]]
(* 1.00000 *)
InputForm[%]
(* 0.99999999994777708655477259390422583293`6.882488164549814 *)

The precision is only 6 digits, which is anomalously low. I'm expecting about 0.999999999948, so that InputForm magnitude is not just a floating-point error. I'm aware that the physical constants involved, including the ones on which the units depend, are not themselves known to very high precision. But it seems to me that the structure of the expression should allow Mathematica to keep a great deal more precision than it shows me.

Answer 1

This is not an answer, but more of a comment to help motivate the question. Apparently the phrasing of the question, as it currently is, is not convincing for many people. However, this is not really related to Quantity. Perhaps it can be an "answer" in the sense that it provides an alternative to relying on the built-in methods.

Definitions

The values Mathematica has for the electron mass and electronvolt are

eV -> 1.60217656500000003606379190082454356534`7.359612335017222*^-19 Quantity[1, "Joules"]
me -> 9.10938291707429769542081766`7.0617630921243855*^-31 Quantity[1, "Kilograms"]

which corresponds in more usual terms to $\mathrm{eV}=1.602176565(44)\times10^{-19}\mathrm{J}$ and $m_\mathrm{e}=9.109382917(87)\times10^{-31}\mathrm{kg}$. These values are more or less correct, although I'm not sure where they come from since the error on the electron mass is about a fifth as large as that for the most recent NIST recommended value. (The electronvolt value seems good, though.)

For reference, NIST can tell us that the current best values are $\mathrm{eV}=1.602176565(35)\times10^{-19}\mathrm{J}$ and $m_\mathrm{e}=9.10938291(40)\times10^{-31}\mathrm{kg}$, and the coefficient of correlation between these is $r=0.9998$. The latter is critically important because it tells us that these uncertainties actually come from the same source, so if we treat the errors as independent, we will be double-counting. Significance arithmetic has no way to deal with correlated uncertainties, so in my view that is already a good enough reason not to use it for this type of calculation.

Here are the values in code:

mma = {
  eV -> 1.60217656500000003606379190082454356534`7.359612335017222*^-19,
  me -> 9.10938291707429769542081766`7.0617630921243855*^-31,
  c -> 299792458
 };

(* padded to 30 places to prevent significance arithmetic from affecting results *)
nist = {
  cov[eV, me] -> 0.9998`30 0.00000040`30*^-31 0.000000035`30*^-19,
  var[eV] -> (0.000000035`30*^-19)^2, eV -> 1.602176565`30*^-19,
  var[me] -> (0.00000040`30*^-31)^2, me -> 9.10938291`30*^-31,
  c -> 299792458
};

The expression is this:

expr = (Sqrt[x eV]*Sqrt[2*c^2*me + x eV])/(c^2*me + x eV);

The expectation value

The mean (expected) value of the expression is almost the same whichever set of values we take for the constants:

mmaresult = expr /. mma /. x -> 50*^9
 (* -> 0.999999999947777086555... *)

expr /. nist /. x -> 50*^9
 (* -> 0.999999999947777086636... *)

We see that the result matches to 18 places. To pre-empt some of the below, this is a symptom of the fact that the expression is not as sensitive to the uncertainty in the constants as Mathematica claims.

For the below, we will use the NIST values.

The uncertainties

Using this package, the "correct" value of the uncertainty is:

Sqrt@PropagateCovariance[expr, {eV, me}, "ExpansionOrder" -> 2] /. nist /. x -> 50*^9
 (* -> 5.9845*10^-17 *)

The value can therefore be stated as $0.999999999947777087(60)$.

Given the approximations it makes, the result significance arithmetic should give is:

Sqrt@PropagateCovariance[expr, {eV, me}, "ExpansionOrder" -> 1, 
  "InitialCovarianceMatrix" -> DiagonalMatrix[{var[eV], var[me]}]
 ] /. nist /. x -> 50*^9
 (* -> 5.1225*10^-18 *)

It corresponds to a value of $0.9999999999477770866(51)$, so the uncertainty is somewhat underestimated. (Actually, the effect of double-counting is not really significant in the first-order expansion.)

The value significance arithmetic actually gives is:

10^-Precision[mmaresult]
 (* -> 8.7382*10^-8 *)

In other words, Mathematica would have us believe that the proper result is $1.000000000 (87)$. But, as we have seen, this is not correct, and the uncertainty is overestimated by many orders of magnitude.

Explanation?

This is really just a guess, but significance arithmetic was never intended to be used for this kind of quantitative work. Its main function is to ensure that Mathematica knows when and where it loses precision, so that it can increase the working precision sufficiently to avoid introducing any spurious values into the calculation. For this reason, it may be that significance arithmetic is actually very conservative in certain cases, such as where a first-order approximation is known not to be sufficient. As long as there are no incorrect digits in the result, significance arithmetic can be said to have done its job. And this is certainly the case here, even if the result is not the one we would have wanted.

Answer 2

Mathematica arrives at the particular precision that it does using significance arithmetic and associated propagation of precision. This answer describes how to double check that. Oleksandr's answer makes a good case for the assertion that significance arithmetic is insufficient for this problem as it fails to account for the correlation in error between quantities involved in the computation. I agree. In general, though, precision is meant to be conservative so we expect underestimates, which is good. One can always call SetPrecision, if justified. Since significance computations are essentially first order, this happens in other contexts - for example, in iterative dynamics in the neighborhood of a point a super-attractive fixed point where $f'(x_0)=0$.

At any rate, expressed as a computation involving dimensionless quantities you're trying to understand the following.

normedGEV = Quantity[1, "Gigaelectronvolts"]/
  Quantity[1, "Kilograms"*"Meters"^2/"Seconds"^2];
normedEMSOL = Quantity[1, "ElectronMass" ("SpeedOfLight")^2]/
  Quantity[1, "Kilograms"*"Meters"^2/"Seconds"^2];
result = Sqrt[50 normedGEV] Sqrt[2 normedEMSOL + 50 normedGEV]/
  (normedEMSOL + 50 normedGEV)
Precision[result]

(* Out: 
    1.000000
    7.05858
*)

Note that the inputs normedGEV and normedEMSOL have associated precisions. We need to start with these precisions to understand the propagation of precision throughout the computation. Whether you think these precisions can be improved based on other sources is a separate question.

Precision /@ {normedGEV, normedEMSOL}
(* Out: {7.35961, 7.06176} *)

Now, the propagation of Precision throughout a computation is described in great detail in the paper "Precise numerical computation" by Mark Sofroniou and Giulia Spaletta. Here are the science direct link and a preprint version. The following definitions are based on that paper.

Attributes[compWithPrecision] = {HoldFirst};
scale[x_] := Log[10, Abs[x]];
compWithPrecision[x_ + y_] := {x + y, Precision[x + y], 
  scale[x + y] + Accuracy[x + y]};
compWithPrecision[x_*y_] := {x*y, Precision[x*y], 
  -Log[10, 10^(-Precision[x]) + 10^(-Precision[y])]};
compWithPrecision[Sqrt[x_]] := {Sqrt[x], Precision[Sqrt[x]], 
  Accuracy[Sqrt[x]] + scale[Sqrt[x]]};

Certainly, this could be expanded. Note, for example, that Precision of product is easily expressible in terms of the Precision of the inputs. Precision of a sum is defined in terms of the Accuracy of its inputs but then Accuracy is defined via

accuracy[x_ + y_] := -Log[10, 10^(-Accuracy[x]) + 10^(-Accuracy[y])]

We can go back and forth easily via

precision[x] == scale[x] + accuracy[x]

We can now simply perform the computations based these definitions, checking that Precision agrees with the stated definitions at each step.

step1 = compWithPrecision[2 normedEMSOL + 50 normedGEV]
(* Out: {8.011047*10^-9, 7.3596, 7.3596} *)

step2 = compWithPrecision[Sqrt[First[step1]]]
(* Out: {0.00008950445, 7.66063, 7.66063} *)

step3 = compWithPrecision[Sqrt[50 normedGEV]]
(* Out: {0.00008950354, 7.66064, 7.66064} *)

step4 = compWithPrecision[First[step2]*First[step3]]
(* Out: {8.010965*10^-9, 7.35961, 7.35961} *)

step5 = compWithPrecision[normedEMSOL + 50 normedGEV]
(* Out: {8.010965*10^-9, 7.35961, 7.35961} *) 

final = compWithPrecision[First[step4]/First[step5]]
(* Out: {1.000000, 7.05858, 7.05858} *)

Answer 3

You can see the very same effect already without Quantity:

a = 1.0`7
(*
==> 1.000000
*)
b = 1.0`14
(*
==> 1.0000000000000
*)
a/Sqrt[a]-Sqrt[a]+b
(*
1.000000
*)

You might argue that there should be more precision here because a cancels out completely. But the point is that as soon as a is evaluated, all Mathematica has is a value of 1.0`7, and there's no longer and indication that all those 1.0`7 come from the same source. Therefore you cannot gain precision from cancelling out, since at the time that cancelling out happens, the algorithm cannot know that the two values are not independent. If you do a simplification before inserting values, you get the desired precision:

Block[{a, b}, Simplify[a/Sqrt[a] - Sqrt[a] + b]]
(*
==> 1.0000000000000
*)

Now in your case, the very same happens: You evaluate the Electron mass twice, but after those evaluations, Mathematica's evaluation algorithm no longer knows that those two values both describe approximations of the same exact value, the electron mass; it has to assume that you've possibly got two independent masses whose approximate value and uncertainty happen to be the same (but whose true value may actually differ). Therefore you cannot gain any precision from cancelling out.