Mathematics arbitrary precision evaluation

Question

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. — J. M.'s missing motivation, Jul 10 '16 at 08:57
@J.M. (1335/4 - a^2)b^6 + a^2(11a^2b^2 - 121b^4 - 2) + 11/2b^8 + a/( 2.0*b) returns the same result, zero (with and without using N function). — syscreat, Jul 10 '16 at 09: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.'s missing motivation, Jul 10 '16 at 09: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 09:21
Well, yes, either use exact numbers, or if you really must have arbitrary-precision inexact numbers, use a backtick, e.g. 5.5`128 or 2`128. — J. M.'s missing motivation, Jul 10 '16 at 09:23

score 3 · Accepted Answer · answered Jul 10 '16 at 10:54

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 *)

score 1 · Answer 2 · answered Jul 10 '16 at 11:41

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

Mathematics arbitrary precision evaluation

2 Answers2