Excuse my lack of vocabulary for I have no formal training in this field, which is also why I ask this question - it may be trivial or it may be impossible.

I want to evaluate an expression in the following form down to a binary float with "correct" rounding and I'm worried that double-rounding will corrupt my result:

$$\sqrt[c]{\frac{a}{b}} \qquad a \in \mathbb N \quad b,c \in \mathbb N_{> 0}$$

By binary float I mostly mean a plain old double, but since I use MPFR it could have any finite precision.

By correct rounding I mean as if the complete expression is evaluated to infinite precision and then rounded correctly, without any chance of "double rounding" corrupting the result.

I use MPFR and GMP to work with these expressions using arbitrary sized integers and floats, and so far I divide the work into two steps, rounding to a larger float in between:

$$t = a/b \qquad t:\mathrm{e.g.\ binary128}$$ $$y = \sqrt[c]{t} \qquad y: \mathrm{e.g.\ binary64}$$

Using MPFR I can calculate both $a/b$ and $\sqrt[c] t$ correctly rounded to any arbitrary precision, but I don't know how to do these two steps in a way that guarantees correct rounding for the complete expression!

Here is a code example of what I'm currently doing. Since it's written in C it is a bit clumsy (no overloads or default arguments) but I hope you can see the steps even if you're not fluent in this language or the libraries I use:

extern unsigned long int a, b, c;     // b and c != 0
MPFR_DECL_INIT( y, DBL_MANT_DIG    ); // Fits a double exactly
MPFR_DECL_INIT( t, DBL_MANT_DIG*2  ); // Twice the precision I want
mpfr_set_ui   ( t, a,    MPFR_RNDN ); // Exact conversion
mpfr_div_ui   ( t, t, b, MPFR_RNDN ); // Correctly rounded
mpfr_rootn_ui ( y, t, c, MPFR_RNDN ); // N'th root correctly rounded

I'm expecting one of these to hold true:

  • Twice the precision for the intermediate is not enough
  • Twice the precision for the intermediate is more than needed - I only need N extra bits
  • I need to do it differently but it's still possible in a finite time and space
  • It's not possible to guarantee correct rounding because of reason
  • 1
    $\begingroup$ The reason double rounding sometimes gives different results is that it is deciding whether $(a/b)^{1/c}$ and $(a/b)^{1/c}(1+\epsilon_1)(1+\epsilon_2)^{1/c}$ are nearest to the same of different floating point number $x$, and $x+\frac12\epsilon(x)$ is the decision boundary. But no matter what, it is conceivable that there might exist $a,b,c$ with $(a/b)^{1/c}<x+\frac12\epsilon(x)<(a/b)^{1/c}(1+\epsilon_1)(1+\epsilon_2)^{1/c}$ causing them to be rounded to $x$ and $x+\epsilon(x)$ respectively, it can only be ruled out in very special cases with very careful analysis. Usually nothing can be done $\endgroup$
    – Kirill
    Aug 18 '19 at 0:11
  • $\begingroup$ Because the roundoff errors $\epsilon_{1,2}$ are smaller than $\epsilon^2$ when using twice the precision, it is also very unlikely to happen. $\endgroup$
    – Kirill
    Aug 18 '19 at 0:13
  • 2
    $\begingroup$ The first statement will hold most of the time. For the few single operations $+,\times,/,√$ for which the required bit lengths in intermediate operations have been proven that are needed to avoid double-rounding issues, one needs at least $2n+i$ bits in intermediate computation when arguments and final result use $n$ bits, where $i \in {0,1,2}$ Canonical reference: Samuel A. Figueroa, "When is double rounding innocuous?", ACM SIGNUM Newsletter , Vol. 30, No. 3, July 1995, pp. 21-26. It stands to reason that compound operations, such as the case at hand, will require additional bits. $\endgroup$
    – njuffa
    Aug 18 '19 at 20:11

Your Answer

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

Browse other questions tagged or ask your own question.