2
$\begingroup$

This question already has an answer here:

I am trying to solve a simple expression:

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

$\endgroup$

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

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 1
    $\begingroup$ For starters, use 11/2 and 1335/4 instead of 5.5 and 333.75. $\endgroup$ – J. M. is away Jul 10 '16 at 8:57
  • $\begingroup$ @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). $\endgroup$ – syscreat Jul 10 '16 at 9:03
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M. is away Jul 10 '16 at 9:10
  • $\begingroup$ Thank you so much! So, as I understand, to let Mathematica correctly evaluate expressions, all variables should have rational and integer types. $\endgroup$ – syscreat Jul 10 '16 at 9:21
  • $\begingroup$ 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. $\endgroup$ – J. M. is away Jul 10 '16 at 9:23
3
$\begingroup$

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 *)
$\endgroup$
1
$\begingroup$

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

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.