4
$\begingroup$

In my calculations I need some larger precision. But due to the fact that I iteratively refine the results to compensate for rounding errors accumulated in previous iteration, Mathematica's arbitraty precision arithmetic's roundoff error tracker kills all the results by converting them to effective zero — I can't even divide by such result without getting divide-by-zero error.

In my application a completely wrong intermediate result is much better (and is actually expected) than no result as I get with error tracking. I've tried using SetPrecision[] at some places in calculation, but in some cases I have to do it in many places, which makes the code unreadable.

So, is there any way to just get additional precision, but without tracking roundoff errors?

Improve this question
$\endgroup$
7
  • 3
    $\begingroup$ I can't picture a case in which your desired behavior would be desirable. Could you provide a minimal example of your specific application? $\endgroup$
    – MarcoB
    Commented Nov 2, 2016 at 16:24
  • $\begingroup$ So you just want to use a fixed and pre-determined number of digits in all calculations, and you want to avoid truncating the number of digits after each operation? And you need more than 15 digits? $\endgroup$
    – Szabolcs
    Commented Nov 2, 2016 at 18:51
  • $\begingroup$ @Szabolcs yes. Basically what e.g. x87's tbyte or IEEE 754's binary128 formats would provide, but potentially with even more precision and range. $\endgroup$
    – Ruslan
    Commented Nov 2, 2016 at 18:55
  • $\begingroup$ How important is performance? I am looking at the ComputerArithmetic package for the first time and it's quite interesting. It lets you define a number format and "arithmetic" with specific rules, then do calculations with it. I would expect it to be slow but I haven't tested how slow exactly. $\endgroup$
    – Szabolcs
    Commented Nov 2, 2016 at 18:58
  • 1
    $\begingroup$ @Szabolcs there is now a good answer. $\endgroup$
    – Ruslan
    Commented Nov 3, 2016 at 20:16

1 Answer 1

Reset to default
10
$\begingroup$

There are a few reasonable ways. I'll illustrate with an example of Newton iterations for square roots, take from this MathGroup post

r[x_, n_] := x - (x^2 - n)/(2*x)
x = 1.0`20;
two = 2.0`20;

First we run it with the usual arithmetic.

Table[x = r[x, two], {30}]

(* Out[680]= {1.4142135623730950488, 1.4142135623730950488, \
1.414213562373095049, 1.414213562373095049, 1.414213562373095049, \
1.41421356237309505, 1.41421356237309505, 1.41421356237309505, \
1.4142135623730950, 1.4142135623730950, 1.4142135623730950, \
1.4142135623730950, 1.414213562373095, 1.414213562373095, \
1.414213562373095, 1.41421356237310, 1.41421356237310, \
1.41421356237310, 1.4142135623731, 1.4142135623731, 1.4142135623731, \
1.4142135623731, 1.414213562373, 1.414213562373, 1.414213562373, \
1.41421356237, 1.41421356237, 1.41421356237, 1.4142135624, \
1.4142135624} *)

We can force fixed precision as below. I'll just show a few iterations.

x = 1.0`20; NumericalMath`FixedPrecisionEvaluate[
 Table[x = r[x, two], {10}], 20]

(* Out[681]= {1.5000000000000000000, 1.4166666666666666667, \
1.4142156862745098039, 1.4142135623746899106, 1.4142135623730950488, \
1.4142135623730950488, 1.4142135623730950488, 1.4142135623730950488, \
1.4142135623730950488, 1.4142135623730950488} *)

Alternatively, use SetPrecision explicitly to reset upward.

x = 1.0`20;
Table[x = SetPrecision[r[x, two], 20], {10}]

(* Out[683]= {1.5000000000000000000, 1.4166666666666666667, \
1.4142156862745098039, 1.4142135623746899106, 1.4142135623730950488, \
1.4142135623730950488, 1.4142135623730950488, 1.4142135623730950488, \
1.4142135623730950488, 1.4142135623730950488} *)

Last is to temporarily set min and max precisions to be equal. Block is good for this type of localized assignment.

x = 1.0`20;
Block[{$MinPrecision = 20, $MaxPrecision = 20}, 
 Table[x = r[x, two], {10}]]

(* Out[685]= {1.5000000000000000000, 1.4166666666666666667, \
1.4142156862745098039, 1.4142135623746899106, 1.4142135623730950488, \
1.4142135623730950488, 1.4142135623730950488, 1.4142135623730950488, \
1.4142135623730950488, 1.4142135623730950488} *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ This is great! A caveat with SetPrecision is that if the result has zero digits of precision, it'll give very precise zero instead of what you'd get with FixedPrecisionEvaluate: see e.g. SetPrecision[((2`2^2`1)^2`1)^2`1, 10]. $\endgroup$
    – Ruslan
    Commented Nov 3, 2016 at 20:13
  • $\begingroup$ Right. SetPrecision as used above is intentionally subverting the precision tracking system, so it should be done only in situations where one knows it is appropriate. $\endgroup$
    – Daniel Lichtblau
    Commented Nov 3, 2016 at 20:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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