# Why does this Quantity calculation lose precision?

(Sqrt[#] Sqrt[Quantity[2, "ElectronMass" ("SpeedOfLight")^2] + #] )/(
Quantity[1, "ElectronMass" ("SpeedOfLight")^2] + #) &[
Quantity[50, "Gigaelectronvolts"]]
(* 1.00000 *)
InputForm[%]
(* 0.999999999947777086554772593904225832936.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.

-
Just the way it displays, and the reason there's an option to change displayed precision. You can also wrap the whole thing in a SetAccuracy or SetPrecision to spew more digits. –  rasher May 1 at 8:06
I know I could do that, but what use is Mathematica's arbitrary-precision system if it gets the precision wrong? –  Escalona May 1 at 8:10
Look at the precision of ElectronMass - it's tracking precision like it should. –  rasher May 1 at 8:14
The precision of the result should be much higher than the precision of ElectronMass, though, because of the structure of the expression. I mean, in the real world, I do know the value to a much higher precision, so why doesn't Mathematica? –  Escalona May 1 at 8:21
@rasher I mean for the magnitude of the error. You say the precision of about 7 is correct, but according to my calculations it should be about 18 given that Mathematica uses the first-order approximation, and assuming the values and uncertainties for the physical constants are correct and up to date. Personally, and this is just my own opinion, I wouldn't trust either the QuantityUnits package or significance arithmetic for careful work where uncertainties need to be carefully considered. –  Oleksandr R. May 1 at 10:14

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.602176565000000036063791900824543565347.359612335017222*^-19 Quantity[1, "Joules"]
me -> 9.109382917074297695420817667.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.602176565000000036063791900824543565347.359612335017222*^-19,
me -> 9.109382917074297695420817667.0617630921243855*^-31,
c -> 299792458
};

(* padded to 30 places to prevent significance arithmetic from affecting results *)
nist = {
cov[eV, me] -> 0.999830 0.0000004030*^-31 0.00000003530*^-19,
var[eV] -> (0.00000003530*^-19)^2, eV -> 1.60217656530*^-19,
var[me] -> (0.0000004030*^-31)^2, me -> 9.1093829130*^-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.

-
N.B. strictly speaking it can be argued whether the first or second order expression is truly the "correct" one, because we do not have the coskewness tensor, and so both results rely on different approximations. Nonetheless, they are comparable to within an order of magnitude, and significantly different to the result at issue here. –  Oleksandr R. May 2 at 11:38

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} *)

-
Very nice answer. I knew of this paper, but I didn't expect that the definitions it presents are still the current ones (I thought that they had probably been improved since then). Incidentally, the compWithPrecision definitions for addition and square root seem to be tautological? Precision is exactly scaled Accuracy, so it raises the question of what is the actual definition for Accuracy here? For multiplication it is okay because the result is expressed in terms of the precision of the inputs. –  Oleksandr R. May 2 at 12:28
@OleksandrR. Thanks - I whipped this together awfully quickly when I saw your answer. I expanded it a bit in light of your comments. I think the two answers are clearly somewhat orthogonal. –  Mark McClure May 2 at 12:49
I definitely agree that the answers are complementary. However, I think you might have misunderstood my point slightly--Mathematica gets a result different to what I would really consider correct even under a first-order approximation and not considering correlation. I think I will have to look at this in more detail later on and write a third answer(!) pointing out where the differences arise. –  Oleksandr R. May 2 at 13:29

You can see the very same effect already without Quantity:

a = 1.07
(*
==> 1.000000
*)
b = 1.014
(*
==> 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.07, and there's no longer and indication that all those 1.07 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.

-
Well noticed! I think that has to be the final, correct answer. –  Oleksandr R. May 2 at 14:39
Using your suggestion and definitions from my answer--Sqrt@Apart[expr^2] /. mma /. x -> 50*^9` gives a precision of 16.866, i.e. a value and uncertainty of $0.999999999947777087(14)$. Precision has indeed been gained and the result certainly agrees well enough with that obtained by algebraic means. –  Oleksandr R. May 2 at 16:20