Issue #7652: Clean up _mpd_qinvroot() and mark it LIBMPDEC_ONLY. Use the

algorithm from decimal.py for mpd_qsqrt().

1 files changed, 140 insertions, 135 deletions

diff --git a/Modules/_decimal/libmpdec/mpdecimal.c b/Modules/_decimal/libmpdec/mpdecimal.c
index b234400..c7d4cbd1 100644
--- a/Modules/_decimal/libmpdec/mpdecimal.c
+++ b/Modules/_decimal/libmpdec/mpdecimal.c
@@ -4298,7 +4298,7 @@ static const mpd_t _mpd_ln10 = {
 /*
  * Set 'result' to log(10).
  *   Ulp error: abs(result - log(10)) < ulp(log(10))
- *   Relative error : abs(result - log(10)) < 5 * 10**-prec * log(10)
+ *   Relative error: abs(result - log(10)) < 5 * 10**-prec * log(10)
  *
  * NOTE: The relative error is not derived from the ulp error, but
  * calculated separately using the fact that 23/10 < log(10) < 24/10.
@@ -7204,6 +7204,19 @@ nanresult:
     mpd_setspecial(r, MPD_POS, MPD_NAN);
 }
 
+/* LIBMPDEC_ONLY */
+/*
+ * Schedule the optimal precision increase for the Newton iteration.
+ *   v := input operand
+ *   z_0 := initial approximation
+ *   initprec := natural number such that abs(sqrt(v) - z_0) < 10**-initprec
+ *   maxprec := target precision
+ *
+ * For convenience the output klist contains the elements in reverse order:
+ *   klist := [k_n-1, ..., k_0], where
+ *     1) k_0 <= initprec and
+ *     2) abs(sqrt(v) - result) < 10**(-2*k_n-1 + 2) <= 10**-maxprec.
+ */
 static inline int
 invroot_schedule_prec(mpd_ssize_t klist[MPD_MAX_PREC_LOG2],
                       mpd_ssize_t maxprec, mpd_ssize_t initprec)
@@ -7224,29 +7237,27 @@ invroot_schedule_prec(mpd_ssize_t klist[MPD_MAX_PREC_LOG2],
 }
 
 /*
- * Initial approximation for the inverse square root.
- *
+ * Initial approximation for the inverse square root function.
  *   Input:
- *     v := 7 or 8 decimal digits with an implicit exponent of 10**-6,
- *          representing a number 1 <= x < 100.
- *
+ *     v := rational number, with 1 <= v < 100
+ *     vhat := floor(v * 10**6)
  *   Output:
- *     An approximation to 1/sqrt(v)
+ *     z := approximation to 1/sqrt(v), such that abs(z - 1/sqrt(v)) < 10**-3.
  */
 static inline void
-_invroot_init_approx(mpd_t *z, mpd_uint_t v)
+_invroot_init_approx(mpd_t *z, mpd_uint_t vhat)
 {
     mpd_uint_t lo = 1000;
     mpd_uint_t hi = 10000;
     mpd_uint_t a, sq;
 
-    assert(v >= lo*lo && v < (hi+1)*(hi+1));
+    assert(lo*lo <= vhat && vhat < (hi+1)*(hi+1));
 
     for(;;) {
         a = (lo + hi) / 2;
         sq = a * a;
-        if (v >= sq) {
-            if (v < sq + 2*a + 1) {
+        if (vhat >= sq) {
+            if (vhat < sq + 2*a + 1) {
                 break;
             }
             lo = a + 1;
@@ -7256,7 +7267,19 @@ _invroot_init_approx(mpd_t *z, mpd_uint_t v)
         }
     }
 
-    /* At this point a/1000 is an approximation to sqrt(v). */
+    /*
+     * After the binary search we have:
+     *  1) a**2 <= floor(v * 10**6) < (a + 1)**2
+     * This implies:
+     *  2) a**2 <= v * 10**6 < (a + 1)**2
+     *  3) a <= sqrt(v) * 10**3 < a + 1
+     * Since 10**3 <= a:
+     *  4) 0 <= 10**prec/a - 1/sqrt(v) < 10**-prec
+     * We have:
+     *  5) 10**3/a - 10**-3 < floor(10**9/a) * 10**-6 <= 10**3/a
+     * Merging 4) and 5):
+     *  6) abs(floor(10**9/a) * 10**-6 - 1/sqrt(v)) < 10**-3
+     */
     mpd_minalloc(z);
     mpd_clear_flags(z);
     z->data[0] = 1000000000UL / a;
@@ -7265,6 +7288,10 @@ _invroot_init_approx(mpd_t *z, mpd_uint_t v)
     mpd_setdigits(z);
 }
 
+/*
+ * Set 'result' to 1/sqrt(a).
+ *   Relative error: abs(result - 1/sqrt(a)) < 10**-prec * 1/sqrt(a)
+ */
 static void
 _mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
               uint32_t *status)
@@ -7282,7 +7309,7 @@ _mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
     mpd_ssize_t ideal_exp, shift;
     mpd_ssize_t adj, tz;
     mpd_ssize_t maxprec, fracdigits;
-    mpd_uint_t x, dummy;
+    mpd_uint_t vhat, dummy;
     int i, n;
 
 
@@ -7301,30 +7328,33 @@ _mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
         fracdigits = v->digits - 1;
         v->exp = -fracdigits;
         n = (v->digits > 7) ? 7 : (int)v->digits;
-        _mpd_get_msdigits(&dummy, &x, v, n);
+        /* Let vhat := floor(v * 10**(2*initprec)) */
+        _mpd_get_msdigits(&dummy, &vhat, v, n);
         if (n < 7) {
-            x *= mpd_pow10[7-n];
+            vhat *= mpd_pow10[7-n];
         }
     }
     else {
         fracdigits = v->digits - 2;
         v->exp = -fracdigits;
         n = (v->digits > 8) ? 8 : (int)v->digits;
-        _mpd_get_msdigits(&dummy, &x, v, n);
+        /* Let vhat := floor(v * 10**(2*initprec)) */
+        _mpd_get_msdigits(&dummy, &vhat, v, n);
         if (n < 8) {
-            x *= mpd_pow10[8-n];
+            vhat *= mpd_pow10[8-n];
         }
     }
     adj = (a->exp-v->exp) / 2;
 
     /* initial approximation */
-    _invroot_init_approx(z, x);
+    _invroot_init_approx(z, vhat);
 
     mpd_maxcontext(&maxcontext);
     mpd_maxcontext(&varcontext);
     varcontext.round = MPD_ROUND_TRUNC;
-    maxprec = ctx->prec + 2;
+    maxprec = ctx->prec + 1;
 
+    /* initprec == 3 */
     i = invroot_schedule_prec(klist, maxprec, 3);
     for (; i >= 0; i--) {
         varcontext.prec = 2*klist[i]+2;
@@ -7358,15 +7388,14 @@ _mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
     mpd_del(&t);
     if (v != &vtmp) mpd_del(v);
     *status |= (workstatus&MPD_Errors);
-    varcontext = *ctx;
-    varcontext.round = MPD_ROUND_HALF_EVEN;
-    mpd_qfinalize(result, &varcontext, status);
+    *status |= (MPD_Rounded|MPD_Inexact);
 }
 
 void
 mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
              uint32_t *status)
 {
+    mpd_context_t workctx;
 
     if (mpd_isspecial(a)) {
         if (mpd_qcheck_nan(result, a, ctx, status)) {
@@ -7391,75 +7420,29 @@ mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
         return;
     }
 
-    _mpd_qinvroot(result, a, ctx, status);
-}
-
-/*
- * Ensure correct rounding. Algorithm after Hull & Abrham, "Properly Rounded
- * Variable Precision Square Root", ACM Transactions on Mathematical Software,
- * Vol. 11, No. 3.
- */
-static void
-_mpd_fix_sqrt(mpd_t *result, const mpd_t *a, mpd_t *tmp,
-              const mpd_context_t *ctx, uint32_t *status)
-{
-    mpd_context_t maxctx;
-    MPD_NEW_CONST(u,0,0,1,1,1,5);
-
-    mpd_maxcontext(&maxctx);
-    u.exp = u.digits - ctx->prec + result->exp - 1;
-
-    _mpd_qsub(tmp, result, &u, &maxctx, status);
-    if (*status&MPD_Errors) goto nanresult;
-
-    _mpd_qmul(tmp, tmp, tmp, &maxctx, status);
-    if (*status&MPD_Errors) goto nanresult;
-
-    if (_mpd_cmp(tmp, a) == 1) {
-        u.exp += 1;
-        u.data[0] = 1;
-        _mpd_qsub(result, result, &u, &maxctx, status);
-    }
-    else {
-        _mpd_qadd(tmp, result, &u, &maxctx, status);
-        if (*status&MPD_Errors) goto nanresult;
-
-        _mpd_qmul(tmp, tmp, tmp, &maxctx, status);
-        if (*status&MPD_Errors) goto nanresult;
-
-        if (_mpd_cmp(tmp, a) == -1) {
-            u.exp += 1;
-            u.data[0] = 1;
-            _mpd_qadd(result, result, &u, &maxctx, status);
-        }
-    }
-
-    return;
-
-nanresult:
-    mpd_setspecial(result, MPD_POS, MPD_NAN);
+    workctx = *ctx;
+    workctx.prec += 2;
+    workctx.round = MPD_ROUND_HALF_EVEN;
+    _mpd_qinvroot(result, a, &workctx, status);
+    mpd_qfinalize(result, ctx, status);
 }
+/* END LIBMPDEC_ONLY */
 
+/* Algorithm from decimal.py */
 void
 mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
           uint32_t *status)
 {
-    uint32_t workstatus = 0;
-    mpd_context_t varcontext;
-    mpd_t *z = result;         /* current approximation */
-    MPD_NEW_STATIC(v,0,0,0,0); /* a, normalized to a number between 1 and 10 */
-    MPD_NEW_STATIC(vtmp,0,0,0,0);
-    MPD_NEW_STATIC(tmp,0,0,0,0);
-    mpd_ssize_t ideal_exp, shift;
-    mpd_ssize_t target_prec, fracdigits;
-    mpd_ssize_t a_exp, a_digits;
-    mpd_ssize_t adj, tz;
-    mpd_uint_t dummy, t;
+    mpd_context_t maxcontext;
+    MPD_NEW_STATIC(c,0,0,0,0);
+    MPD_NEW_STATIC(q,0,0,0,0);
+    MPD_NEW_STATIC(r,0,0,0,0);
+    MPD_NEW_CONST(two,0,0,1,1,1,2);
+    mpd_ssize_t prec, ideal_exp;
+    mpd_ssize_t l, shift;
     int exact = 0;
 
 
-    varcontext = *ctx;
-    varcontext.round = MPD_ROUND_HALF_EVEN;
     ideal_exp = (a->exp - (a->exp & 1)) / 2;
 
     if (mpd_isspecial(a)) {
@@ -7483,79 +7466,101 @@ mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
         return;
     }
 
-    if (!mpd_qcopy(&v, a, status)) {
-        mpd_seterror(result, MPD_Malloc_error, status);
-        goto finish;
-    }
+    mpd_maxcontext(&maxcontext);
+    prec = ctx->prec + 1;
 
-    a_exp = a->exp;
-    a_digits = a->digits;
+    if (!mpd_qcopy(&c, a, status)) {
+        goto malloc_error;
+    }
+    c.exp = 0;
 
-    /* normalize a to 1 <= v < 100 */
-    if ((v.digits+v.exp) & 1) {
-        fracdigits = v.digits - 1;
-        v.exp = -fracdigits;
-        _mpd_get_msdigits(&dummy, &t, &v, 3);
-        t = t < 100 ? t*10 : t;
-        t = t < 100 ? t*10 : t;
+    if (a->exp & 1) {
+        if (!mpd_qshiftl(&c, &c, 1, status)) {
+            goto malloc_error;
+        }
+        l = (a->digits >> 1) + 1;
     }
     else {
-        fracdigits = v.digits - 2;
-        v.exp = -fracdigits;
-        _mpd_get_msdigits(&dummy, &t, &v, 4);
-        t = t < 1000 ? t*10 : t;
-        t = t < 1000 ? t*10 : t;
-        t = t < 1000 ? t*10 : t;
+        l = (a->digits + 1) >> 1;
     }
-    adj = (a_exp-v.exp) / 2;
-
 
-    /* use excess digits */
-    target_prec = (a_digits > ctx->prec) ? a_digits : ctx->prec;
-    target_prec += 2;
-    varcontext.prec = target_prec + 3;
-
-    /* invroot is much faster for large numbers */
-    _mpd_qinvroot(&tmp, &v, &varcontext, &workstatus);
+    shift = prec - l;
+    if (shift >= 0) {
+        if (!mpd_qshiftl(&c, &c, 2*shift, status)) {
+            goto malloc_error;
+        }
+        exact = 1;
+    }
+    else {
+        exact = !mpd_qshiftr_inplace(&c, -2*shift);
+    }
 
-    varcontext.prec = target_prec;
-    _mpd_qdiv(NO_IDEAL_EXP, z, &one, &tmp, &varcontext, &workstatus);
+    ideal_exp -= shift;
 
+    /* find result = floor(sqrt(c)) using Newton's method */
+    if (!mpd_qshiftl(result, &one, prec, status)) {
+        goto malloc_error;
+    }
 
-    tz = mpd_trail_zeros(result);
-    if ((result->digits-tz)*2-1 <= v.digits) {
-        _mpd_qmul(&tmp, result, result, &varcontext, &workstatus);
-        if (workstatus&MPD_Errors) {
-            mpd_seterror(result, workstatus&MPD_Errors, status);
-            goto finish;
+    while (1) {
+        _mpd_qdivmod(&q, &r, &c, result, &maxcontext, &maxcontext.status);
+        if (mpd_isspecial(result) || mpd_isspecial(&q)) {
+            mpd_seterror(result, maxcontext.status&MPD_Errors, status);
+            goto out;
+        }
+        if (_mpd_cmp(result, &q) <= 0) {
+            break;
         }
-        exact = (_mpd_cmp(&tmp, &v) == 0);
+        _mpd_qadd_exact(result, result, &q, &maxcontext, &maxcontext.status);
+        if (mpd_isspecial(result)) {
+            mpd_seterror(result, maxcontext.status&MPD_Errors, status);
+            goto out;
+        }
+        _mpd_qdivmod(result, &r, result, &two, &maxcontext, &maxcontext.status);
     }
-    *status |= (workstatus&MPD_Errors);
 
-    if (!exact && !mpd_isspecial(result) && !mpd_iszero(result)) {
-        _mpd_fix_sqrt(result, &v, &tmp, &varcontext, status);
-        if (mpd_isspecial(result)) goto finish;
-        *status |= (MPD_Rounded|MPD_Inexact);
+    if (exact) {
+        _mpd_qmul_exact(&r, result, result, &maxcontext, &maxcontext.status);
+        if (mpd_isspecial(&r)) {
+            mpd_seterror(result, maxcontext.status&MPD_Errors, status);
+            goto out;
+        }
+        exact = (_mpd_cmp(&r, &c) == 0);
     }
 
-    result->exp += adj;
     if (exact) {
-        shift = ideal_exp - result->exp;
-        shift = (tz > shift) ? shift : tz;
-        if (shift > 0) {
+        if (shift >= 0) {
             mpd_qshiftr_inplace(result, shift);
-            result->exp += shift;
+        }
+        else {
+            if (!mpd_qshiftl(result, result, -shift, status)) {
+                goto malloc_error;
+            }
+        }
+        ideal_exp += shift;
+    }
+    else {
+        int lsd = mpd_lsd(result->data[0]);
+        if (lsd == 0 || lsd == 5) {
+            result->data[0] += 1;
         }
     }
 
+    result->exp = ideal_exp;
 
-finish:
-    mpd_del(&v);
-    mpd_del(&vtmp);
-    mpd_del(&tmp);
-    varcontext.prec = ctx->prec;
-    mpd_qfinalize(result, &varcontext, status);
+
+out:
+    mpd_del(&c);
+    mpd_del(&q);
+    mpd_del(&r);
+    maxcontext = *ctx;
+    maxcontext.round = MPD_ROUND_HALF_EVEN;
+    mpd_qfinalize(result, &maxcontext, status);
+    return;
+
+malloc_error:
+    mpd_seterror(result, MPD_Malloc_error, status);
+    goto out;
 }