mpf_ln.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 /*
  16
  17 Using the newton approximation y1 = y - (e^y - x)/e^y
  18
  19 Which converges quickly and we can reuse e^y so we only calc it once per loop
  20
  21 */
  22 int  mpf_ln(mp_float *a, mp_float *b)
  23 {
  24     mp_float  oldval, tmpey, tmpy, val;
  25     int       itts, err;
  26     long      k;
  27
  28     /* ensure positive */
  29     if (a->mantissa.sign == MP_NEG) {
  30        return MP_VAL;
  31     }
  32
  33     /* easy out for 0 */
  34     if (mpf_iszero(a) == MP_YES) {
  35        return mpf_const_d(b, 1);
  36     }
  37
  38     /* initialize temps */
  39     if ((err = mpf_init_multi(b->radix, &oldval, &tmpey, &tmpy, &val, NULL)) != MP_OKAY) {
  40        return err;
  41     }
  42
  43     /* initial guess */
  44     if ((err = mpf_const_e(&val)) != MP_OKAY)                                           { goto __ERR; }
  45     if ((err = mpf_sqr(&val, &val)) != MP_OKAY)                                         { goto __ERR; }
  46     if ((err = mpf_inv(&val, &tmpey)) != MP_OKAY)                                       { goto __ERR; }
  47     if ((err = mpf_copy(a, &tmpy)) != MP_OKAY)                                          { goto __ERR; }
  48     if ((err = mpf_normalize_to(&tmpy, b->radix)) != MP_OKAY)                           { goto __ERR; }
  49
  50     /* divide out e's */
  51     k = 0;
  52     while (mpf_cmp(&tmpy, &val) == MP_GT) {
  53         ++k;
  54         if ((err = mpf_mul(&tmpy, &tmpey, &tmpy)) != MP_OKAY)                           { goto __ERR; }
  55     }
  56     if ((err = mpf_const_d(&tmpy, k*2)) != MP_OKAY)                                     { goto __ERR; }
  57
  58     /* number of iterations */
  59     itts = mpf_iterations(b);
  60
  61     while (itts--) {
  62         if ((err = mpf_copy(&tmpy, &oldval)) != MP_OKAY)                                { goto __ERR; }
  63
  64         /* get e^y and save it */
  65         if ((err = mpf_exp(&tmpy, &tmpey)) != MP_OKAY)                                  { goto __ERR; }
  66
  67         /* now compute e^y - x */
  68         if ((err = mpf_sub(&tmpey, a, &val)) != MP_OKAY)                                { goto __ERR; }
  69
  70         /* now compute (e^y - x) / e^y */
  71         if ((err = mpf_div(&val, &tmpey, &val)) != MP_OKAY)                             { goto __ERR; }
  72
  73         /* y = y - (e^y - x)/e^y */
  74         if ((err = mpf_sub(&tmpy, &val, &tmpy)) != MP_OKAY)                             { goto __ERR; }
  75
  76         if (mpf_cmp(&tmpy, &oldval) == MP_EQ) {
  77            break;
  78         }
  79     }
  80     mpf_exch(&tmpy, b);
  81 __ERR: mpf_clear_multi(&oldval, &tmpey, &tmpy, &val, NULL);
  82    return err;
  83 }