sysdeps/ieee754/dbl-64/mpatan.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2014 Free Software Foundation, Inc.
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU Lesser General Public License as published by
   8  * the Free Software Foundation; either version 2.1 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU  Lesser General Public License
  17  * along with this program; if not, see <http://www.gnu.org/licenses/>.
  18  */
  19 /******************************************************************/
  20 /*                                                                */
  21 /* MODULE_NAME:mpatan.c                                           */
  22 /*                                                                */
  23 /* FUNCTIONS:mpatan                                               */
  24 /*                                                                */
  25 /* FILES NEEDED: mpa.h endian.h mpatan.h                          */
  26 /*               mpa.c                                            */
  27 /*                                                                */
  28 /* Multi-Precision Atan function subroutine, for precision p >= 4.*/
  29 /* The relative error of the result is bounded by 34.32*r**(1-p), */
  30 /* where r=2**24.                                                 */
  31 /******************************************************************/
  32
  33 #include "endian.h"
  34 #include "mpa.h"
  35
  36 #ifndef SECTION
  37 # define SECTION
  38 #endif
  39
  40 #include "mpatan.h"
  41
  42 void
  43 SECTION
  44 __mpatan (mp_no *x, mp_no *y, int p)
  45 {
  46   int i, m, n;
  47   double dx;
  48   mp_no mptwoim1 =
  49   {
  50     0,
  51     {
  52       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  53       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  54       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
  55     }
  56   };
  57
  58   mp_no mps, mpsm, mpt, mpt1, mpt2, mpt3;
  59
  60   /* Choose m and initiate mptwoim1.  */
  61   if (EX > 0)
  62     m = 7;
  63   else if (EX < 0)
  64     m = 0;
  65   else
  66     {
  67       __mp_dbl (x, &dx, p);
  68       dx = ABS (dx);
  69       for (m = 6; m > 0; m--)
  70         {
  71           if (dx > __atan_xm[m].d)
  72             break;
  73         }
  74     }
  75   mptwoim1.e = 1;
  76   mptwoim1.d[0] = 1;
  77
  78   /* Reduce x m times.  */
  79   __sqr (x, &mpsm, p);
  80   if (m == 0)
  81     __cpy (x, &mps, p);
  82   else
  83     {
  84       for (i = 0; i < m; i++)
  85         {
  86           __add (&mpone, &mpsm, &mpt1, p);
  87           __mpsqrt (&mpt1, &mpt2, p);
  88           __add (&mpt2, &mpt2, &mpt1, p);
  89           __add (&mptwo, &mpsm, &mpt2, p);
  90           __add (&mpt1, &mpt2, &mpt3, p);
  91           __dvd (&mpsm, &mpt3, &mpt1, p);
  92           __cpy (&mpt1, &mpsm, p);
  93         }
  94       __mpsqrt (&mpsm, &mps, p);
  95       mps.d[0] = X[0];
  96     }
  97
  98   /* Evaluate a truncated power series for Atan(s).  */
  99   n = __atan_np[p];
 100   mptwoim1.d[1] = __atan_twonm1[p].d;
 101   __dvd (&mpsm, &mptwoim1, &mpt, p);
 102   for (i = n - 1; i > 1; i--)
 103     {
 104       mptwoim1.d[1] -= 2;
 105       __dvd (&mpsm, &mptwoim1, &mpt1, p);
 106       __mul (&mpsm, &mpt, &mpt2, p);
 107       __sub (&mpt1, &mpt2, &mpt, p);
 108     }
 109   __mul (&mps, &mpt, &mpt1, p);
 110   __sub (&mps, &mpt1, &mpt, p);
 111
 112   /* Compute Atan(x).  */
 113   mptwoim1.d[1] = 1 << m;
 114   __mul (&mptwoim1, &mpt, y, p);
 115 }