2 /*
3  * IBM Accurate Mathematical Library
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 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 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, write to the Free Software
18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19  */
20 /******************************************************************/
21 /*                                                                */
22 /* MODULE_NAME:mpatan.c                                           */
23 /*                                                                */
24 /* FUNCTIONS:mpatan                                               */
25 /*                                                                */
26 /* FILES NEEDED: mpa.h endian.h mpatan.h                          */
27 /*               mpa.c                                            */
28 /*                                                                */
29 /* Multi-Precision Atan function subroutine, for precision p >= 4.*/
30 /* The relative error of the result is bounded by 34.32*r**(1-p), */
31 /* where r=2**24.                                                 */
32 /******************************************************************/
34 #include "endian.h"
35 #include "mpa.h"
36 void mpsqrt(mp_no *, mp_no *, int);
38 void mpatan(mp_no *x, mp_no *y, int p) {
39 #include "mpatan.h"
41   int i,m,n;
42   double dx;
43   mp_no
44     mpone    = {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.0,
45                 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,
46                 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
47     mptwo    = {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.0,
48                 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,
49                 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
50     mptwoim1 = {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.0,
51                 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,
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}};
54   mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3;
56                       /* Choose m and initiate mpone, mptwo & mptwoim1 */
57     if      (EX>0) m=7;
58     else if (EX<0) m=0;
59     else {
60       __mp_dbl(x,&dx,p);  dx=ABS(dx);
61       for (m=6; m>0; m--)
62         {if (dx>xm[m].d) break;}
63     }
64     mpone.e    = mptwo.e    = mptwoim1.e = 1;
65     mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE;
66     mptwo.d[1] = TWO;
68                                  /* Reduce x m times */
69     mul(x,x,&mpsm,p);
70     if (m==0) cpy(x,&mps,p);
71     else {
72       for (i=0; i<m; i++) {
74         mpsqrt(&mpt1,&mpt2,p);
78         dvd(&mpsm,&mpt3,&mpt1,p);
79         cpy(&mpt1,&mpsm,p);
80       }
81       mpsqrt(&mpsm,&mps,p);    mps.d[0] = X[0];
82     }
84                     /* Evaluate a truncated power series for Atan(s) */
85     n=np[p];    mptwoim1.d[1] = twonm1[p].d;
86     dvd(&mpsm,&mptwoim1,&mpt,p);
87     for (i=n-1; i>1; i--) {
88       mptwoim1.d[1] -= TWO;
89       dvd(&mpsm,&mptwoim1,&mpt1,p);
90       mul(&mpsm,&mpt,&mpt2,p);
91       sub(&mpt1,&mpt2,&mpt,p);
92     }
93     mul(&mps,&mpt,&mpt1,p);
94     sub(&mps,&mpt1,&mpt,p);
96                           /* Compute Atan(x) */
97     mptwoim1.d[1] = twom[m].d;
98     mul(&mptwoim1,&mpt,y,p);
100   return;
101 }