# 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?

• For starters, use 11/2 and 1335/4 instead of 5.5 and 333.75. Jul 10, 2016 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). Jul 10, 2016 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. Jul 10, 2016 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. Jul 10, 2016 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. Jul 10, 2016 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
(* -0.82739605995 *)
c/. values
(* 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