mpf_exp.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 /* compute b = e^a using e^x == \sum_{n=0}^{\infty} {1 \over n!}x^n */
  17 int  mpf_exp(mp_float *a, mp_float *b)
  18 {
  19    mp_float  oldval, tmpx, tmpovern, tmp, res;
  20    int       err, itts;
  21    long      n;
  22
  23    /* initialize temps */
  24    if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, NULL)) != MP_OKAY) {
  25       return err;
  26    }
  27
  28    /* initlialize temps */
  29    /* all three start at one */
  30    if ((err = mpf_const_d(&res, 1)) != MP_OKAY)                                         { goto __ERR; }
  31    if ((err = mpf_const_d(&tmpovern, 1)) != MP_OKAY)                                    { goto __ERR; }
  32    if ((err = mpf_const_d(&tmpx, 1)) != MP_OKAY)                                        { goto __ERR; }
  33    n = 1;
  34
  35    /* get number of iterations */
  36    itts = mpf_iterations(b);
  37
  38    while (itts-- > 0) {
  39        if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                                  { goto __ERR; }
  40
  41        /* compute 1/n! as 1/(n-1)! * 1/n */
  42        if ((err = mpf_const_d(&tmp, n++)) != MP_OKAY)                                   { goto __ERR; }
  43        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
  44        if ((err = mpf_mul(&tmp, &tmpovern, &tmpovern)) != MP_OKAY)                      { goto __ERR; }
  45
  46        /* compute x^n as x^(n-1) * x */
  47        if ((err = mpf_mul(&tmpx, a, &tmpx)) != MP_OKAY)                                 { goto __ERR; }
  48
  49        /* multiply and sum them */
  50        if ((err = mpf_mul(&tmpovern, &tmpx, &tmp)) != MP_OKAY)                          { goto __ERR; }
  51        if ((err = mpf_add(&tmp, &res, &res)) != MP_OKAY)                                { goto __ERR; }
  52
  53        if (mpf_cmp(&oldval, &res) == MP_EQ) {
  54           break;
  55        }
  56    }
  57
  58    mpf_exch(&res, b);
  59 __ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, NULL);
  60    return err;
  61 }
  62