Merge 8.7

Undo Tcl-specific changes in bn_mp_sqrt.c, and re-enable the two disabled test-cases: This proves the observed crash with DIBIT_BIT=60 is caused by those Tcl-specific changes!

1 files changed, 4 insertions, 68 deletions

diff --git a/libtommath/bn_mp_sqrt.c b/libtommath/bn_mp_sqrt.c
index bbca158..ea93e74 100644
--- a/libtommath/bn_mp_sqrt.c
+++ b/libtommath/bn_mp_sqrt.c
@@ -12,20 +12,11 @@
  * SPDX-License-Identifier: Unlicense
  */
 
-#ifndef NO_FLOATING_POINT
-#include <math.h>
-#endif
-
 /* this function is less generic than mp_n_root, simpler and faster */
 int mp_sqrt(const mp_int *arg, mp_int *ret)
 {
    int res;
    mp_int t1, t2;
-   int i, j, k;
-#ifndef NO_FLOATING_POINT
-   volatile double d;
-   mp_digit dig;
-#endif
 
    /* must be positive */
    if (arg->sign == MP_NEG) {
@@ -33,14 +24,12 @@ int mp_sqrt(const mp_int *arg, mp_int *ret)
    }
 
    /* easy out */
-   if (mp_iszero(arg) == MP_YES) {
+   if (mp_iszero(arg)) {
       mp_zero(ret);
       return MP_OKAY;
    }
 
-   i = (arg->used / 2) - 1;
-   j = 2 * i;
-   if ((res = mp_init_size(&t1, i+2)) != MP_OKAY) {
+   if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) {
       return res;
    }
 
@@ -48,61 +37,8 @@ int mp_sqrt(const mp_int *arg, mp_int *ret)
       goto E2;
    }
 
-   for (k = 0; k < i; ++k) {
-      t1.dp[k] = (mp_digit) 0;
-   }
-
-#ifndef NO_FLOATING_POINT
-
-   /* Estimate the square root using the hardware floating point unit. */
-
-   d = 0.0;
-   for (k = arg->used-1; k >= j; --k) {
-      d = ldexp(d, DIGIT_BIT) + (double)(arg->dp[k]);
-   }
-
-   /*
-    * At this point, d is the nearest floating point number to the most
-    * significant 1 or 2 mp_digits of arg. Extract its square root.
-    */
-
-   d = sqrt(d);
-
-   /* dig is the most significant mp_digit of the square root */
-
-   dig = (mp_digit) ldexp(d, -DIGIT_BIT);
-
-   /*
-    * If the most significant digit is nonzero, find the next digit down
-    * by subtracting DIGIT_BIT times thie most significant digit.
-    * Subtract one from the result so that our initial estimate is always
-    * low.
-    */
-
-   if (dig) {
-      t1.used = i+2;
-      d -= ldexp((double) dig, DIGIT_BIT);
-      if (d >= 1.0) {
-         t1.dp[i+1] = dig;
-         t1.dp[i] = ((mp_digit) d) - 1;
-      } else {
-         t1.dp[i+1] = dig-1;
-         t1.dp[i] = MP_DIGIT_MAX;
-      }
-   } else {
-      t1.used = i+1;
-      t1.dp[i] = ((mp_digit) d) - 1;
-   }
-
-#else
-
-   /* Estimate the square root as having 1 in the most significant place. */
-
-   t1.used = i + 2;
-   t1.dp[i+1] = (mp_digit) 1;
-   t1.dp[i] = (mp_digit) 0;
-
-#endif
+   /* First approx. (not very bad for large arg) */
+   mp_rshd(&t1, t1.used/2);
 
    /* t1 > 0  */
    if ((res = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {