Mfrhtyj

I'm really not too sure how to go on about this:

Problem. Julia (or take MatLab etc. (I guess.)) computes

In [1]: exp(pi*sqrt(67)/6) - sqrt(5280)

Out[1]: 6.121204876308184e-8

Only by looking closely: How many digits are correct? (What could you do to compute $16$ decimals correctly?)

So I know that subtraction is somewhat dangerous numerically - though how do I solve this problem in a precise manner?

Using a multi-precision package like the CAS Magma you can compute a value where perhaps the last 10 of 200 digits are polluted by truncation and rounding errors. The loop using different target digit numbers

for N in [0..5] do 

  RR := RealField(40+32*N); 

  Exp(Pi(RR)*(RR!67)^(1/2)/6) - (RR!5280)^(1/2); 

end for;

gives the result

0.000000061212043348401464497838533606

0.00000006121204334840146449783853360999610568722086907308823071985234

0.0000000612120433484014644978385336099961056872208690730882307198520281414760909167862633336012797951

0.00000006121204334840146449783853360999610568722086907308823071985202814147609091678626333360127979547576192009412710141319214385241

0.00000006121204334840146449783853360999610568722086907308823071985202814147609091678626333360127979547576192009412710141319214385238281932550279412177187787689095831

0.0000000612120433484014644978385336099961056872208690730882307198520281414760909167862633336012797954757619200941271014131921438523828193255027941217718778768909589395497833016775104278900532217412

This shows empirically that all but the last 5 digits are correct. Compared to your floating point result

0.000000061212048763081839

we see that 7 digits are correct. Using the modified formula by Andrei
$$
frac{e^{2pisqrt{67}/6}-5280}{e^{pisqrt{67}/6}+sqrt{5280}}
$$
gives the floating point result

0.000000061212049284920693 

which is further away from the true result, it does not cure the catastrophic cancellation.

If you have two numbers $a,b$ that coincide in $d$ leading digits, then for any smooth function $f$ you will also get coincidence in about $d$ leading digits of $f(a)$ and $f(b)$. This means that the difference $f(a)-f(b)$ has $15-d$ correct digits. There is no way to force more precision, as that is not present in the original floating point values.

Note that this is different from the situation where you want to compute $f(a+h)-f(a)$ where $h$ is explicitly known. Then you can apply Taylor formulas or just a mean value argument like $f(a+h)-f(a)approx f(a+h/2)cdot h$. This will have almost as many correct digits as $h$ has for $|h|<10^{-5}cdot |a|$.

To assess the accuracy of your computation analytically you need to understand the errors in each piece. The subtraction is a problem because the numbers are represented by a fixed number of bits of mantissa. This represents a fractional error in the number of about $2^{-n}$ where $n$ is the number of bits. If we know $n$ for the computer in use we know what is left after the subtraction. The math functions should be accurate to (just about) the same number of bits.

We can compute $sqrt{5280} approx 72.6636$ so the size of the difference is about $frac {6.212cdot 10^{-8}}{72.6636}approx 8.5cdot 10^{-10}$ of the numbers themselves. The base $2$ log of this is about $-30$ so you are losing $30$ bits to the subtraction. If you use $64$ bit floats according to the IEEE specification, that leaves $23$ bits of accuracy. As $2^{-23} approx 10^{-7}$ you have about seven decimal digits of accuracy.