2
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ For starters, use 11/2 and 1335/4 instead of 5.5 and 333.75. $\endgroup$
    – J. M.'s eventual burnout ♦
    Jul 10, 2016 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, 2016 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.'s eventual burnout ♦
    Jul 10, 2016 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, 2016 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.'s eventual burnout ♦
    Jul 10, 2016 at 9:23

2 Answers 2

Reset to default
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 *)
Improve this answer
$\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

Improve this answer
$\endgroup$

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