util/src/TclTk/tcl8.6.12/libtommath/bn_mp_sqrt.c

   1 #include "tommath_private.h"
   2 #ifdef BN_MP_SQRT_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 #ifndef NO_FLOATING_POINT
   7 #include <float.h>
   8 #include <math.h>
   9 #if (MP_DIGIT_BIT != 28) || (FLT_RADIX != 2) || (DBL_MANT_DIG != 53) || (DBL_MAX_EXP != 1024)
  10 #define NO_FLOATING_POINT
  11 #endif
  12 #endif
  13
  14 /* this function is less generic than mp_n_root, simpler and faster */
  15 mp_err mp_sqrt(const mp_int *arg, mp_int *ret)
  16 {
  17    mp_err err;
  18    mp_int t1, t2;
  19 #ifndef NO_FLOATING_POINT
  20    int i, j, k;
  21    volatile double d;
  22    mp_digit dig;
  23 #endif
  24
  25    /* must be positive */
  26    if (arg->sign == MP_NEG) {
  27       return MP_VAL;
  28    }
  29
  30    /* easy out */
  31    if (MP_IS_ZERO(arg)) {
  32       mp_zero(ret);
  33       return MP_OKAY;
  34    }
  35
  36 #ifndef NO_FLOATING_POINT
  37
  38    i = (arg->used / 2) - 1;
  39    j = 2 * i;
  40    if ((err = mp_init_size(&t1, i+2)) != MP_OKAY) {
  41       return err;
  42    }
  43
  44    if ((err = mp_init(&t2)) != MP_OKAY) {
  45       goto E2;
  46    }
  47
  48    for (k = 0; k < i; ++k) {
  49       t1.dp[k] = (mp_digit) 0;
  50    }
  51
  52    /* Estimate the square root using the hardware floating point unit. */
  53
  54    d = 0.0;
  55    for (k = arg->used-1; k >= j; --k) {
  56       d = ldexp(d, MP_DIGIT_BIT) + (double)(arg->dp[k]);
  57    }
  58
  59    /*
  60     * At this point, d is the nearest floating point number to the most
  61     * significant 1 or 2 mp_digits of arg. Extract its square root.
  62     */
  63
  64    d = sqrt(d);
  65
  66    /* dig is the most significant mp_digit of the square root */
  67
  68    dig = (mp_digit) ldexp(d, -MP_DIGIT_BIT);
  69
  70    /*
  71     * If the most significant digit is nonzero, find the next digit down
  72     * by subtracting MP_DIGIT_BIT times thie most significant digit.
  73     * Subtract one from the result so that our initial estimate is always
  74     * low.
  75     */
  76
  77    if (dig) {
  78       t1.used = i+2;
  79       d -= ldexp((double) dig, MP_DIGIT_BIT);
  80       if (d >= 1.0) {
  81          t1.dp[i+1] = dig;
  82          t1.dp[i] = ((mp_digit) d) - 1;
  83       } else {
  84          t1.dp[i+1] = dig-1;
  85          t1.dp[i] = MP_DIGIT_MAX;
  86       }
  87    } else {
  88       t1.used = i+1;
  89       t1.dp[i] = ((mp_digit) d) - 1;
  90    }
  91
  92 #else
  93
  94    if ((err = mp_init_copy(&t1, arg)) != MP_OKAY) {
  95       return err;
  96    }
  97
  98    if ((err = mp_init(&t2)) != MP_OKAY) {
  99       goto E2;
 100    }
 101
 102    /* First approx. (not very bad for large arg) */
 103    mp_rshd(&t1, t1.used/2);
 104
 105 #endif
 106
 107    /* t1 > 0  */
 108    if ((err = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
 109       goto E1;
 110    }
 111    if ((err = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
 112       goto E1;
 113    }
 114    if ((err = mp_div_2(&t1, &t1)) != MP_OKAY) {
 115       goto E1;
 116    }
 117    /* And now t1 > sqrt(arg) */
 118    do {
 119       if ((err = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
 120          goto E1;
 121       }
 122       if ((err = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
 123          goto E1;
 124       }
 125       if ((err = mp_div_2(&t1, &t1)) != MP_OKAY) {
 126          goto E1;
 127       }
 128       /* t1 >= sqrt(arg) >= t2 at this point */
 129    } while (mp_cmp_mag(&t1, &t2) == MP_GT);
 130
 131    mp_exch(&t1, ret);
 132
 133 E1:
 134    mp_clear(&t2);
 135 E2:
 136    mp_clear(&t1);
 137    return err;
 138 }
 139
 140 #endif