mp_sqrt.c

   1 #include "tommath_private.h"
   2 #ifdef MP_SQRT_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* this function is less generic than mp_n_root, simpler and faster */
   7 mp_err mp_sqrt(const mp_int *arg, mp_int *ret)
   8 {
   9    mp_err err;
  10    mp_int t1, t2;
  11
  12    /* must be positive */
  13    if (mp_isneg(arg)) {
  14       return MP_VAL;
  15    }
  16
  17    /* easy out */
  18    if (mp_iszero(arg)) {
  19       mp_zero(ret);
  20       return MP_OKAY;
  21    }
  22
  23    if ((err = mp_init_copy(&t1, arg)) != MP_OKAY) {
  24       return err;
  25    }
  26
  27    if ((err = mp_init(&t2)) != MP_OKAY) {
  28       goto LBL_ERR2;
  29    }
  30
  31    /* First approx. (not very bad for large arg) */
  32    mp_rshd(&t1, t1.used/2);
  33
  34    /* t1 > 0  */
  35    if ((err = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
  36       goto LBL_ERR1;
  37    }
  38    if ((err = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
  39       goto LBL_ERR1;
  40    }
  41    if ((err = mp_div_2(&t1, &t1)) != MP_OKAY) {
  42       goto LBL_ERR1;
  43    }
  44    /* And now t1 > sqrt(arg) */
  45    do {
  46       if ((err = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
  47          goto LBL_ERR1;
  48       }
  49       if ((err = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
  50          goto LBL_ERR1;
  51       }
  52       if ((err = mp_div_2(&t1, &t1)) != MP_OKAY) {
  53          goto LBL_ERR1;
  54       }
  55       /* t1 >= sqrt(arg) >= t2 at this point */
  56    } while (mp_cmp_mag(&t1, &t2) == MP_GT);
  57
  58    mp_exch(&t1, ret);
  59
  60 LBL_ERR1:
  61    mp_clear(&t2);
  62 LBL_ERR2:
  63    mp_clear(&t1);
  64    return err;
  65 }
  66
  67 #endif