sysdeps/ieee754/dbl-64/mpatan.c

   1 /*
   2  * IBM Accurate Mathematical Library
   3  * written by International Business Machines Corp.
   4  * Copyright (C) 2001-2016 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 #include <math.h>
  36
  37 #ifndef SECTION
  38 # define SECTION
  39 #endif
  40
  41 #include "mpatan.h"
  42
  43 void
  44 SECTION
  45 __mpatan (mp_no *x, mp_no *y, int p)
  46 {
  47   int i, m, n;
  48   double dx;
  49   mp_no mptwoim1 =
  50   {
  51     0,
  52     {
  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, 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       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
  56     }
  57   };
  58
  59   mp_no mps, mpsm, mpt, mpt1, mpt2, mpt3;
  60
  61   /* Choose m and initiate mptwoim1.  */
  62   if (EX > 0)
  63     m = 7;
  64   else if (EX < 0)
  65     m = 0;
  66   else
  67     {
  68       __mp_dbl (x, &dx, p);
  69       dx = fabs (dx);
  70       for (m = 6; m > 0; m--)
  71         {
  72           if (dx > __atan_xm[m].d)
  73             break;
  74         }
  75     }
  76   mptwoim1.e = 1;
  77   mptwoim1.d[0] = 1;
  78
  79   /* Reduce x m times.  */
  80   __sqr (x, &mpsm, p);
  81   if (m == 0)
  82     __cpy (x, &mps, p);
  83   else
  84     {
  85       for (i = 0; i < m; i++)
  86         {
  87           __add (&__mpone, &mpsm, &mpt1, p);
  88           __mpsqrt (&mpt1, &mpt2, p);
  89           __add (&mpt2, &mpt2, &mpt1, p);
  90           __add (&__mptwo, &mpsm, &mpt2, p);
  91           __add (&mpt1, &mpt2, &mpt3, p);
  92           __dvd (&mpsm, &mpt3, &mpt1, p);
  93           __cpy (&mpt1, &mpsm, p);
  94         }
  95       __mpsqrt (&mpsm, &mps, p);
  96       mps.d[0] = X[0];
  97     }
  98
  99   /* Evaluate a truncated power series for Atan(s).  */
 100   n = __atan_np[p];
 101   mptwoim1.d[1] = __atan_twonm1[p].d;
 102   __dvd (&mpsm, &mptwoim1, &mpt, p);
 103   for (i = n - 1; i > 1; i--)
 104     {
 105       mptwoim1.d[1] -= 2;
 106       __dvd (&mpsm, &mptwoim1, &mpt1, p);
 107       __mul (&mpsm, &mpt, &mpt2, p);
 108       __sub (&mpt1, &mpt2, &mpt, p);
 109     }
 110   __mul (&mps, &mpt, &mpt1, p);
 111   __sub (&mps, &mpt1, &mpt, p);
 112
 113   /* Compute Atan(x).  */
 114   mptwoim1.d[1] = 1 << m;
 115   __mul (&mptwoim1, &mpt, y, p);
 116 }