mpf_invsqrt.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 /* using newtons method we have 1/sqrt(x) = Y_{n+1} = y_n * ((3 - xy^2_n)/2) */
  16 int  mpf_invsqrt(mp_float *a, mp_float *b)
  17 {
  18    mp_float oldval, tmp1, tmp2, const_3;
  19    int err, itts;
  20
  21    /* ensure a is not zero or negative */
  22    if ((mp_iszero(&(a->mantissa)) == MP_YES) || (a->mantissa.sign == MP_NEG)) {
  23       return MP_VAL;
  24    }
  25
  26    /* get number of iterations */
  27    itts = mpf_iterations(b);
  28
  29    /* init temps */
  30    if ((err = mpf_init_multi(b->radix, &oldval, &tmp1, &tmp2, &const_3, NULL)) != MP_OKAY) {
  31       return err;
  32    }
  33
  34    /* const_3 */
  35    if ((err = mpf_const_d(&const_3, 3)) != MP_OKAY)                               { goto __ERR; }
  36
  37    /* tmp1 == reasonable guess at sqrt */
  38    if ((err = mpf_copy(a, &tmp1)) != MP_OKAY)                                     { goto __ERR; }
  39    mp_rshd(&(tmp1.mantissa), tmp1.mantissa.used>>1);
  40    if ((err = mpf_normalize_to(&tmp1, b->radix)) != MP_OKAY)                      { goto __ERR; }
  41
  42    while (itts-- > 0) {
  43        /* grap copy of tmp1 for early out */
  44        if ((err = mpf_copy(&tmp1, &oldval)) != MP_OKAY)                           { goto __ERR; }
  45
  46        /* first tmp2 = y^2 == tmp1^2 */
  47        if ((err = mpf_sqr(&tmp1, &tmp2)) != MP_OKAY)                              { goto __ERR; }
  48        /* multiply by x, tmp1 * a */
  49        if ((err = mpf_mul(&tmp2, a, &tmp2)) != MP_OKAY)                           { goto __ERR; }
  50        /* 3 - xy^2_n == 3 - tmp1 */
  51        if ((err = mpf_sub(&const_3, &tmp2, &tmp2)) != MP_OKAY)                    { goto __ERR; }
  52        /* halve it */
  53        if ((err = mpf_div_2(&tmp2, &tmp2)) != MP_OKAY)                            { goto __ERR; }
  54        /* multiply by y_n and feedback */
  55        if ((err = mpf_mul(&tmp1, &tmp2, &tmp1)) != MP_OKAY)                       { goto __ERR; }
  56
  57        /* early out if stable */
  58        if (mpf_cmp(&oldval, &tmp1) == MP_EQ) {
  59           break;
  60        }
  61    }
  62
  63    mpf_exch(&tmp1, b);
  64 __ERR: mpf_clear_multi(&oldval, &tmp1, &tmp2, &const_3, NULL);
  65    return err;
  66 }
  67