3 * IBM Accurate Mathematical Library
4 * Copyright (c) International Business Machines Corp., 2001
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.
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.
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.
20 /******************************************************************/
22 /* MODULE_NAME:mpatan.c */
24 /* FUNCTIONS:mpatan */
26 /* FILES NEEDED: mpa.h endian.h mpatan.h */
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), */
32 /******************************************************************/
36 void __mpsqrt(mp_no *, mp_no *, int);
38 void __mpatan(mp_no *x, mp_no *y, int p) {
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 */
60 __mp_dbl(x,&dx,p); dx=ABS(dx);
62 {if (dx>xm[m].d) break;}
64 mpone.e = mptwo.e = mptwoim1.e = 1;
65 mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE;
68 /* Reduce x m times */
70 if (m==0) __cpy(x,&mps,p);
73 __add(&mpone,&mpsm,&mpt1,p);
74 __mpsqrt(&mpt1,&mpt2,p);
75 __add(&mpt2,&mpt2,&mpt1,p);
76 __add(&mptwo,&mpsm,&mpt2,p);
77 __add(&mpt1,&mpt2,&mpt3,p);
78 __dvd(&mpsm,&mpt3,&mpt1,p);
81 __mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0];
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--) {
89 __dvd(&mpsm,&mptwoim1,&mpt1,p);
90 __mul(&mpsm,&mpt,&mpt2,p);
91 __sub(&mpt1,&mpt2,&mpt,p);
93 __mul(&mps,&mpt,&mpt1,p);
94 __sub(&mps,&mpt1,&mpt,p);
97 mptwoim1.d[1] = twom[m].d;
98 __mul(&mptwoim1,&mpt,y,p);