Computing a sum with high precision

I am aware that a similar question has been asked, for example, here and here, but none of the proposed solutions solves my problem, so please wait before you mark this question as a duplicate.

Consider this simple problem where one has two bodies of masses mA, mB, their relative position BA, and the position CC of their center of mass:

mA = 398600.435436;
mB = 4902.800066;
CC = 0.12051741410138465477
BA = -0.00080817735147818490


One wants to know the positions, xB and xA, of the two bodies, which are given by

xA = (mA CC + mB CC - mB BA)/(mA + mB)
xB = (mB CC + mA (CC + BA))/(mA + mB)


If I evaluate xA and xB numerically in Mathematica 10.3.1.0, and then compute the relative distance between the two bodies, I should get BA. However, there is a small numerical difference:

In[47]:= xB - xA - BA
Out[47]= -2.168404344971009*10^-18


Do you know how to increase the numerical precision in order to decrease the difference above?

• I'm voting to close this question as off-topic because its really a question about numerical analysis, not Mathematica per se. – Daniel Lichtblau Jul 21 '18 at 14:56
• @Daniel Lichtblau I don't read it as a question about numerical analysis. Edward clearly understands that it's a precision issue. His question is specifically about increasing precision of a calculation in Mathematica. – John Doty Jul 21 '18 at 15:47
• @JohnDoty I guess, but is that a meaningful thing to do? There are cases where one does do that (which is actually something that will be in a talk I give next week). I don't see how raising precision applies here though-- the stated purpose makes me think one should use numerical analysis methods to determine when "small" means "treat as zero". – Daniel Lichtblau Jul 21 '18 at 16:22

Like this?

mA = 398600.4354360000000000000000000000000000000000000000000000000000000000;
mB = 4902.800066000000000000000000000000000000000000000000000000000000000000;
CC = 0.120517414101384654770000000000000000000000000000000000000000000000000;
BA = -0.00080817735147818490000000000000000000000000000000000000000000000000;
xA = (mA CC + mB CC - mB BA)/(mA + mB);
xB = (mB CC + mA (CC + BA))/(mA + mB);
xB - xA - BA


which instantly displays

0.*10^-68


but you almost certainly do not know that those four constants do have about 60 trailing zeros.

Note there are always several different ways of accomplishing anything in Mathematica. This is only one of the ways of telling Mathematica that you have lots of trailing zeros on a number, but it is perhaps the easiest way to really know and see what you have done and requires the least knowledge of the tens of thousands of Mathematica functions that you might not know.

Now if you do the following

Clear[mA, mB, CC, BA];
xA = (mA CC + mB CC - mB BA)/(mA + mB);
xB = (mB CC + mA (CC + BA))/(mA + mB);
FullSimplify[xB - xA - BA]


then that will instantly display

0


and that tells you that your final expression should be zero for infinitely precise values of your four quantities

• prec=100; mA = SetPrecision[398600.435436, prec]; etc. would be a less cumbersome way to increase the precision of the numbers. – Henrik Schumacher Jul 21 '18 at 14:42
• @Bill: That works, thank you. – Edward Jul 21 '18 at 14:42
• You could also use mA = 398600.435436100 – Carl Woll Jul 21 '18 at 18:36

You can use exact numbers:

r[x_] := Rationalize[x, 0]
mA = 398600.435436 // r;
mB = 4902.800066 // r;
CC = 0.12051741410138465477 // r;
BA = -0.00080817735147818490 // r;
xA = (mA CC + mB CC - mB BA)/(mA + mB);
xB = (mB CC + mA (CC + BA))/(mA + mB);
xB - xA - BA


Yielding:

0

• Or set precision: r[x_] := SetPrecision[x, 100], :) – Mariusz Iwaniuk Jul 21 '18 at 15:22
• @Mariusz Iwaniuk, SetPrecision[x,100] yields an approximate number of only 100 digits precision. Rationalize[x,0] yields an exact number. Neither, of course, reflects the precision of the input data, but the exact approach verifies the identity xB - xA - BA==0 exactly, while the approximate approach does so only approximately. – John Doty Jul 21 '18 at 15:36
• Ok.it was supposed to be: r[x_] := SetPrecision[x, Infinity].:P` – Mariusz Iwaniuk Jul 21 '18 at 16:19