# Mathematics arbitrary precision evaluation [duplicate]

I am trying to solve a simple expression:

where a = 77617, and b = 33096.

Wolfram|Alpha returns a correct result, using following form:

a = 77617, b = 33096,
c = (333.75 - a^2) * b^6 + a^2 * (11 * a^2 * b^2 - 121 * b^4 - 2) + 5.5 * b^8 + a/(2.0 * b)


-54767/66192

which is approximately -0.827396.

Mathematica, when I evaluate the following code:

ClearAll[a, b, c]
a = N[77617, 128];
b = N[33096, 128];
N[(333.75 - a^2)*b^6 + a^2*(11*a^2*b^2 - 121*b^4 - 2) + 5.5*b^8 + a/(
2.0*b), 128]


returns 0.

Could you please to point how to get a correct result using Mathematica?

## marked as duplicate by Mr.Wizard♦Jul 10 '16 at 12:04

• For starters, use 11/2 and 1335/4 instead of 5.5 and 333.75. – J. M. will be back soon Jul 10 '16 at 8:57
• @J.M. (1335/4 - a^2)*b^6 + a^2*(11*a^2*b^2 - 121*b^4 - 2) + 11/2*b^8 + a/( 2.0*b) returns the same result, zero (with and without using N function). – syscreat Jul 10 '16 at 9:03
• I did not notice the 2.0; replace that with just 2. What you must know here is that numbers with decimal points and no precision indicators are machine precision numbers, which was why your original expression was not being evaluated to arbitrary precision. – J. M. will be back soon Jul 10 '16 at 9:10
• Thank you so much! So, as I understand, to let Mathematica correctly evaluate expressions, all variables should have rational and integer types. – syscreat Jul 10 '16 at 9:21
• Well, yes, either use exact numbers, or if you really must have arbitrary-precision inexact numbers, use a backtick, e.g. 5.5128 or 2128. – J. M. will be back soon Jul 10 '16 at 9:23

In working with approximate numbers, it is important to ensure that no terms you introduce (e.g. 5.5 or 333.75) are lower precision than you need to work to.

Wherever you have exact coefficients, it is good practice to specify these as rational to avoid this problem.

The particular problem you have here is that the terms in b^6 and b^8 are almost equal in magnitude but opposite in sign, requiring high precision to evaluate the result without catastrophic loss of precision.

We can investigate the behaviour with the following code

values[n_] := N[{a -> 77617, b -> 33096}, n];
c = (333 + 3/4 - a^2)*b^6 + a^2*(11*a^2*b^2 - 121*b^4 - 2) +  11/2*b^8 + a/(2*b);


For example

c /. values[∞]
(* -(54767/66192) *)
c /. values[50]
(* -0.82739605995 *)
c/. values[30]
(* 0.*10^8 *)


I would do it this way.

c = With[{a = 77617, b = 33096},
Map[
Rationalize,
Hold[(333.75 - a^2)*b^6 + a^2*(11*a^2*b^2 - 121*b^4 - 2) + 5.5*b^8 + a/(2.0*b)],
{-1}] // ReleaseHold]


-(54767/66192)

I think this is likely to be pretty close to the way Wolfram|Alpha does it.

If you really want the result to 128 decimal places, you can now evaluate

 N[c, 128]


-0.8273960599468213681411650954798162919990331157843848199178148416727\ 0969301426154218032390621223108532753202803964225284022238337