2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001, 2011 Free Software Foundation
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.
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.
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/>.
19 /******************************************************************/
21 /* MODULE_NAME:mpatan.c */
23 /* FUNCTIONS:mpatan */
25 /* FILES NEEDED: mpa.h endian.h mpatan.h */
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), */
31 /******************************************************************/
42 void __mpsqrt(mp_no
*, mp_no
*, int);
46 __mpatan(mp_no
*x
, mp_no
*y
, int p
) {
51 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,
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}},
54 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,
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,0.0,
56 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}},
57 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,
58 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,
59 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
61 mp_no mps
,mpsm
,mpt
,mpt1
,mpt2
,mpt3
;
63 /* Choose m and initiate mpone, mptwo & mptwoim1 */
67 __mp_dbl(x
,&dx
,p
); dx
=ABS(dx
);
69 {if (dx
>__atan_xm
[m
].d
) break;}
71 mpone
.e
= mptwo
.e
= mptwoim1
.e
= 1;
72 mpone
.d
[0] = mpone
.d
[1] = mptwo
.d
[0] = mptwoim1
.d
[0] = ONE
;
75 /* Reduce x m times */
77 if (m
==0) __cpy(x
,&mps
,p
);
80 __add(&mpone
,&mpsm
,&mpt1
,p
);
81 __mpsqrt(&mpt1
,&mpt2
,p
);
82 __add(&mpt2
,&mpt2
,&mpt1
,p
);
83 __add(&mptwo
,&mpsm
,&mpt2
,p
);
84 __add(&mpt1
,&mpt2
,&mpt3
,p
);
85 __dvd(&mpsm
,&mpt3
,&mpt1
,p
);
88 __mpsqrt(&mpsm
,&mps
,p
); mps
.d
[0] = X
[0];
91 /* Evaluate a truncated power series for Atan(s) */
92 n
=__atan_np
[p
]; mptwoim1
.d
[1] = __atan_twonm1
[p
].d
;
93 __dvd(&mpsm
,&mptwoim1
,&mpt
,p
);
94 for (i
=n
-1; i
>1; i
--) {
96 __dvd(&mpsm
,&mptwoim1
,&mpt1
,p
);
97 __mul(&mpsm
,&mpt
,&mpt2
,p
);
98 __sub(&mpt1
,&mpt2
,&mpt
,p
);
100 __mul(&mps
,&mpt
,&mpt1
,p
);
101 __sub(&mps
,&mpt1
,&mpt
,p
);
103 /* Compute Atan(x) */
104 mptwoim1
.d
[1] = __atan_twom
[m
].d
;
105 __mul(&mptwoim1
,&mpt
,y
,p
);