mpf_sin.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 sin x == \sum_{n=0}^{\infty} ((-1)^n/(2n+1)!) * x^(2n+1) */
  16 int  mpf_sin(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
  22    /* initialize temps */
  23    if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
  24       return err;
  25    }
  26
  27    /* initlialize temps */
  28    /* this is the 1/n! which starts at 1 */
  29    if ((err = mpf_const_d(&tmpovern, 1)) != MP_OKAY)                                    { goto __ERR; }
  30
  31    /* the square of the input, used to save multiplications later */
  32    if ((err = mpf_sqr(a, &sqr)) != MP_OKAY)                                             { goto __ERR; }
  33
  34    /* tmpx starts at the input, so we copy and normalize */
  35    if ((err = mpf_copy(a, &tmpx)) != MP_OKAY)                                           { goto __ERR; }
  36    tmpx.radix = b->radix;
  37    if ((err = mpf_normalize(&tmpx)) != MP_OKAY)                                         { goto __ERR; }
  38
  39    /* the result starts off at a as we skip a term in the series */
  40    if ((err = mpf_copy(&tmpx, &res)) != MP_OKAY)                                        { goto __ERR; }
  41
  42    /* this is the denom counter.  Goes up by two per pass */
  43    n       = 1;
  44
  45    /* we alternate between adding and subtracting */
  46    oddeven = 1;
  47
  48    /* get number of iterations */
  49    itts = mpf_iterations(b);
  50
  51    while (itts-- > 0) {
  52        if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                                  { goto __ERR; }
  53
  54        /* compute 1/(2n)! from 1/(2(n-1))! by multiplying by (1/n)(1/(n+1)) */
  55        if ((err = mpf_const_d(&tmp, ++n)) != MP_OKAY)                                   { goto __ERR; }
  56        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
  57        if ((err = mpf_mul(&tmpovern, &tmp, &tmpovern)) != MP_OKAY)                      { goto __ERR; }
  58        /* we do this twice */
  59        if ((err = mpf_const_d(&tmp, ++n)) != MP_OKAY)                                   { goto __ERR; }
  60        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
  61        if ((err = mpf_mul(&tmpovern, &tmp, &tmpovern)) != MP_OKAY)                      { goto __ERR; }
  62
  63        /* now multiply sqr into tmpx */
  64        if ((err = mpf_mul(&tmpx, &sqr, &tmpx)) != MP_OKAY)                              { goto __ERR; }
  65
  66        /* now multiply the two */
  67        if ((err = mpf_mul(&tmpx, &tmpovern, &tmp)) != MP_OKAY)                          { goto __ERR; }
  68
  69        /* now depending on if this is even or odd we add/sub */
  70        oddeven ^= 1;
  71        if (oddeven  == 1) {
  72           if ((err = mpf_add(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
  73        } else {
  74           if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
  75        }
  76
  77        if (mpf_cmp(&res, &oldval) == MP_EQ) {
  78           break;
  79        }
  80    }
  81    mpf_exch(&res, b);
  82 __ERR: mpf_clear_multi(&oldval, &tmpx, &tmpovern, &tmp, &res, &sqr, NULL);
  83    return err;
  84 }