# Making a calculation with high precision

I would like to make the following calculation:

1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.)] - 1


Mathematica 8 returns 0. The result is obviously not 0, but my calculation must be losing precision. How can I compute this value with sufficient precision?

• Just remove the . after 16 – rm -rf May 23 '13 at 14:36

For the sake of alternatives, here's another way. You can specify a precision with a backtick. This is called an "arbitrary-precision" number. Such numbers are treated differently than machine precision numbers.

Examples of 3 and 5 digits of precision:

16.3

16.0

16.5

16.000


Machine precision depends on the machine, but it's often nearly 16 digits.

$MachinePrecision  15.9546  The results of your calculation at different precisions are given below. Note they are more accurate than machine precision, even though the specified precisions are technically less. That's because Mathematica can use up to $MaxExtraPrecision extra digits of precision in its internal calculations.

1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.3)] - 1
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.4)] - 1
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.5)] - 1
1/Sqrt[1 - (150^2 10^(-4))/(9 10^16.15)] - 1

1.3*10^-17
1.25*10^-17
1.250*10^-17
1.2500000000000*10^-17


Which way is better, this one or using N and exact numbers as in the answers of @chuy and @m_goldberg? Well it depends. Normally N[exact] is sufficient, although in this case, it fails -- it gives 0. Arbitrary-precision calculations take longer, but Mathematica does try to keep track of the precision as it goes along. On the other hand, N[exact, prec] (@m_goldberg's answer) tracks the precision and, more importantly, it gives an answer with the target precision, if the internal precision limit $MaxExtraPrecision is not exceeded. Update: Alternative approach The difficulty in this calculation is a classic problem in numerical analysis, so please forgive me if this is well-known. The reason can be seen in the numerator of $${1 \over \sqrt{1-a}} -1 = {1 - \sqrt{1 - a} \over \sqrt{1 -a}} \,.$$ If$a$is very small, then we're subtracting two numbers close to one another, resulting in a loss pf precision. In fact,$\sqrt{1 - a}$may be equal to$1$when rounded to machine precision. A solution is to rationalize the numerator to get $${a \over \big(1 + \sqrt{1 - a} \big) \sqrt{1 -a}} = {a \over 1-a+\sqrt{1-a}}\,.$$ In this form, it evaluates with machine precision without a problem. a/(1 + Sqrt[1 - a] - a) /. a -> (150^2 10^(-4))/(9 10^16.)  1.25*10^-17  Even N has less trouble with this form. Here we limit the extra precision to force the difficulty earlier: Block[{$MaxExtraPrecision = 4},
N[a/(1 + Sqrt[1 - a] - a) /. a -> (150^2 10^(-4))/(9 10^16), 20]]

1.2500000000000000234*10^-17

Block[{$MaxExtraPrecision = 4}, N[1/Sqrt[1 - a] - 1 /. a -> (150^2 10^(-4))/(9 10^16), 20]]  N::meprec: Internal precision limit$MaxExtraPrecision = 4. reached while evaluating -1+200000000/Sqrt. >>

1.250000*10^-17


This same problem occurs in simple arithmetic:

1. - 1. + 10.^-16
1. + 10.^-16 - 1.

1.*10^-16
0.


In sum, another solution would be to choose a formula that avoids such numerical problems. (I don't know of a way to get Mathematica to figure this out automatically.) While it is rare that N cannot handle a computation, machine precision calculations are fast and to be desired if possible.

• Thank you very much for your answer! There is really usefull information in your answer! Thank's! – Thanos May 24 '13 at 6:07
• @Thanos Thanks. I added another solution that in some ways is better but more work for the user, and I hope I made it clear that I prefer m_goldberg's way of dealing with the given calculation – Michael E2 May 24 '13 at 13:14
• In particular, $MachineEpsilon is the smallest h such that 1. + h - 1. is nonzero. Smaller than that (as in Michael's example), and you get subtractive cancellation. – J. M.'s technical difficulties May 24 '13 at 13:15 Ditch the decimal point after the 16 (this makes it a machine precision number): 1/Sqrt[1 - (150^2 10^(-4))/(9 10^16)] - 1 yeilds: -1 + 200000000/Sqrt • Thank you very much for your answer! If I remove the . I don't get an actual number. I just need better precision. – Thanos May 23 '13 at 14:40 • What do you mean you don't get an actual number, @Thanos? This is a number, symbolically presented. I suppose you will say that Sqrt is not a number, too? In any event, look up N[], and note well what its second argument does. – J. M.'s technical difficulties May 23 '13 at 14:42 • @J.M.: he he he! It is a true number indeed! all I need is something in the format 1e-95 or something like that. – Thanos May 23 '13 at 14:43 • Hence my third sentence, @Thanos. Did you read the docs already? – J. M.'s technical difficulties May 23 '13 at 14:45 • @Thanos, then it seems you don't want the maximum precision (which your question explicitly asks for). – chuy May 23 '13 at 15:05 In V11.0 we try: N[1/Sqrt[1 - (150^2 10^(-4))/(9 10^16)] - 1]  we get an answer of 0. As Michael E2 points out we are getting subtractive cancellation causing the loss of significant digits. So we say $$a=\frac{22500}{10^4 \left(9 \cdot 10^{16}\right)}$$ In numerical work the Taylor Series is your friend, it is used often. We will try to get a better answer just with a series. We expand a around 0. Series[1/Sqrt[1 - a] - 1, {a, 0, 3}]//Normal  yields$\frac{5 a^3}{16}+\frac{3 a^2}{8}+\frac{a}{2}$This is numerically very stable as it has small numerators and no subtraction of nearly equal quantities. We Hornerize it by using HornerForm[a/2 + (3 a^2)/8 + (5 a^3)/16]  yields$a (1/2 + (3/8 + (5 a)/16) a) \$

N[a (1/2 + (3/8 + (5 a)/16) a) /. a -> (22500 10^(-4))/(9 10^16)]


yields $$1.25 \cdot 10^{-17}$$

You can get a high-precision approximation by using N with two arguments, the second argument specifying the precision you desire.

a = 1/Sqrt[1 - (150^2 10^(-4))/(9 10^16)] - 1;
Table[ N[a, p], {p, 17, 22}] // TableForm
`

1.2500000000000000*10^-17
1.25000000000000002*10^-17
1.250000000000000023*10^-17
1.2500000000000000234*10^-17
1.25000000000000002344*10^-17
1.250000000000000023438*10^-17

• Thank you very much for your answer! – Thanos May 24 '13 at 6:10