mpf_sqrt.c

   1 /* LibTomFloat, multiple-precision floating-point library
   2  *
   3  * LibTomFloat is a library that provides multiple-precision
   4  * floating-point artihmetic as well as trigonometric functionality.
   5  *
   6  * This library requires the public domain LibTomMath to be installed.
   7  *
   8  * This library is free for all purposes without any express
   9  * gurantee it works
  10  *
  11  * Tom St Denis, tomstdenis@iahu.ca, http://float.libtomcrypt.org
  12  */
  13 #include <tomfloat.h>
  14
  15 /* compute using sqrt(x) = y_{n+1} = 0.5 * (y_{n} + x/y_{n}) */
  16 int  mpf_sqrt(mp_float *a, mp_float *b)
  17 {
  18    mp_float oldval, tmp, res;
  19    int      err, itts;
  20
  21    /* make sure it's positive */
  22    if (a->mantissa.sign == MP_NEG) {
  23       return MP_VAL;
  24    }
  25
  26    /* let's roll */
  27    if ((err = mpf_init_multi(b->radix, &oldval, &tmp, &res, NULL)) != MP_OKAY) {
  28       return err;
  29    }
  30
  31    /* get copy of a and make reasonable guestimte towards the sqrt */
  32    if ((err = mpf_copy(a, &res)) != MP_OKAY)                                     { goto __ERR; }
  33    mp_rshd(&(res.mantissa), res.mantissa.used >> 1);
  34    res.exp   = res.exp / 2;
  35    if ((err = mpf_normalize_to(&res, b->radix)) != MP_OKAY)                      { goto __ERR; }
  36
  37    /* number of iterations */
  38    itts = mpf_iterations(b);
  39
  40    while (itts--) {
  41         if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                          { goto __ERR; }
  42
  43         /* compute x/res */
  44         if ((err = mpf_div(a, &res, &tmp)) != MP_OKAY)                           { goto __ERR; }
  45         /* res + x/res */
  46         if ((err = mpf_add(&res, &tmp, &res)) != MP_OKAY)                        { goto __ERR; }
  47         /* 0.5 * (res + x/res) */
  48         if ((err = mpf_div_2(&res, &res)) != MP_OKAY)                            { goto __ERR; }
  49
  50         if (mpf_cmp(&res, &oldval) == MP_EQ) {
  51            break;
  52         }
  53    }
  54
  55    mpf_exch(&res, b);
  56 __ERR: mpf_clear_multi(&oldval, &tmp, &res, NULL);
  57    return err;
  58 }
  59