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?

$\endgroup$
  • 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 Nov 2 '16 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 Nov 2 '16 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 Nov 2 '16 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 Nov 2 '16 at 18:58
  • 1
    $\begingroup$ @Szabolcs there is now a good answer. $\endgroup$ – Ruslan Nov 3 '16 at 20:16
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} *)
$\endgroup$
  • $\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 Nov 3 '16 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 Nov 3 '16 at 20:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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