mpf_cos.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 cos x == \sum_{n=0}^{\infty} ((-1)^n/(2n)!) * x^2n */
  16 int  mpf_cos(mp_float *a, mp_float *b)
  17 {
  18    mp_float oldval, tmpovern, tmp, tmpx, res, sqr;
  19    int      oddeven, err, itts;
  20    long     n;
  21    /* initialize temps */
  22    if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
  23       return err;
  24    }
  25
  26    /* initlialize temps */
  27    /* three start at one, sqr is the square of a */
  28    if ((err = mpf_const_d(&res, 1)) != MP_OKAY)                                         { goto __ERR; }
  29    if ((err = mpf_const_d(&tmpovern, 1)) != MP_OKAY)                                    { goto __ERR; }
  30    if ((err = mpf_const_d(&tmpx, 1)) != MP_OKAY)                                        { goto __ERR; }
  31    if ((err = mpf_sqr(a, &sqr)) != MP_OKAY)                                             { goto __ERR; }
  32
  33    /* this is the denom counter.  Goes up by two per pass */
  34    n       = 0;
  35
  36    /* we alternate between adding and subtracting */
  37    oddeven = 1;
  38
  39    /* get number of iterations */
  40    itts = mpf_iterations(b);
  41
  42    while (itts-- > 0) {
  43        if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                                  { goto __ERR; }
  44        /* compute 1/(2n)! from 1/(2(n-1))! by multiplying by (1/n)(1/(n+1)) */
  45        if ((err = mpf_const_d(&tmp, ++n)) != MP_OKAY)                                   { goto __ERR; }
  46        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
  47        if ((err = mpf_mul(&tmpovern, &tmp, &tmpovern)) != MP_OKAY)                      { goto __ERR; }
  48        /* we do this twice */
  49        if ((err = mpf_const_d(&tmp, ++n)) != MP_OKAY)                                   { goto __ERR; }
  50        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
  51        if ((err = mpf_mul(&tmpovern, &tmp, &tmpovern)) != MP_OKAY)                      { goto __ERR; }
  52
  53        /* now multiply a into tmpx twice */
  54        if ((err = mpf_mul(&tmpx, &sqr, &tmpx)) != MP_OKAY)                              { goto __ERR; }
  55
  56        /* now multiply the two */
  57        if ((err = mpf_mul(&tmpx, &tmpovern, &tmp)) != MP_OKAY)                          { goto __ERR; }
  58
  59        /* now depending on if this is even or odd we add/sub */
  60        oddeven ^= 1;
  61        if (oddeven  == 1) {
  62           if ((err = mpf_add(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
  63        } else {
  64           if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
  65        }
  66
  67        if (mpf_cmp(&res, &oldval) == MP_EQ) {
  68           break;
  69        }
  70    }
  71    mpf_exch(&res, b);
  72 __ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
  73    return err;
  74 }