mpf_atan.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 /* y = y - (tan(y) - x)/(tan(y+0.1)-tan(x)) */
  16
  17 int  mpf_atan(mp_float *a, mp_float *b)
  18 {
  19    mp_float oldval, tmp, tmpx, res, sqr;
  20    int      oddeven, ires, err, itts;
  21    long     n;
  22
  23    /* ensure -1 <= a <= 1 */
  24    if ((err = mpf_cmp_d(a, -1, &ires)) != MP_OKAY) {
  25       return err;
  26    }
  27    if (ires == MP_LT) {
  28       return MP_VAL;
  29    }
  30
  31    if ((err = mpf_cmp_d(a, 1, &ires)) != MP_OKAY) {
  32       return err;
  33    }
  34    if (ires == MP_GT) {
  35       return MP_VAL;
  36    }
  37
  38    /* easy out if a == 0 */
  39    if (mpf_iszero(a) == MP_YES) {
  40       return mpf_const_d(b, 1);
  41    }
  42
  43    /* now a != 0 */
  44
  45    /* initialize temps */
  46    if ((err = mpf_init_multi(b->radix, &oldval, &tmpx, &tmp, &res, &sqr, NULL)) != MP_OKAY) {
  47       return err;
  48    }
  49
  50    /* initlialize temps */
  51    /* res = 0 */
  52    /* tmpx = 1/a */
  53    if ((err = mpf_inv(a, &tmpx)) != MP_OKAY)                                            { goto __ERR; }
  54
  55    /* sqr = a^2 */
  56    if ((err = mpf_sqr(a, &sqr)) != MP_OKAY)                                             { goto __ERR; }
  57
  58    /* this is the denom counter.  Goes up by two per pass */
  59    n       = 1;
  60
  61    /* we alternate between adding and subtracting */
  62    oddeven = 0;
  63
  64    /* get number of iterations */
  65    itts = mpf_iterations(b);
  66
  67    while (itts-- > 0) {
  68        if ((err = mpf_copy(&res, &oldval)) != MP_OKAY)                                  { goto __ERR; }
  69
  70        /* compute 1/(2n-1) */
  71        if ((err = mpf_const_d(&tmp, (2*n++ - 1))) != MP_OKAY)                           { goto __ERR; }
  72        if ((err = mpf_inv(&tmp, &tmp)) != MP_OKAY)                                      { goto __ERR; }
  73
  74        /* now multiply a into tmpx twice */
  75        if ((err = mpf_mul(&tmpx, &sqr, &tmpx)) != MP_OKAY)                              { goto __ERR; }
  76
  77        /* now multiply the two */
  78        if ((err = mpf_mul(&tmpx, &tmp, &tmp)) != MP_OKAY)                               { goto __ERR; }
  79
  80        /* now depending on if this is even or odd we add/sub */
  81        oddeven ^= 1;
  82        if (oddeven  == 1) {
  83           if ((err = mpf_add(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
  84        } else {
  85           if ((err = mpf_sub(&res, &tmp, &res)) != MP_OKAY)                             { goto __ERR; }
  86        }
  87
  88        if (mpf_cmp(&oldval, &res) == MP_EQ) {
  89           break;
  90        }
  91    }
  92    mpf_exch(&res, b);
  93 __ERR: mpf_clear_multi(&oldval, &tmpx, &tmp, &res, &sqr, NULL);
  94    return err;
  95
  96 }