]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
aeb25823 | 3 | * written by International Business Machines Corp. |
2b778ceb | 4 | * Copyright (C) 2001-2021 Free Software Foundation, Inc. |
e4d82761 UD |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
50944bca | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License | |
5a82c748 | 17 | * along with this program; if not, see <https://www.gnu.org/licenses/>. |
e4d82761 UD |
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" | |
0e9be4db | 35 | #include <math.h> |
e4d82761 | 36 | |
31d3cc00 UD |
37 | #ifndef SECTION |
38 | # define SECTION | |
39 | #endif | |
40 | ||
e4d82761 | 41 | #include "mpatan.h" |
50944bca | 42 | |
31d3cc00 UD |
43 | void |
44 | SECTION | |
6295157a SP |
45 | __mpatan (mp_no *x, mp_no *y, int p) |
46 | { | |
6295157a | 47 | int i, m, n; |
e4d82761 | 48 | double dx; |
6295157a SP |
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 | }; | |
e4d82761 | 58 | |
6295157a | 59 | mp_no mps, mpsm, mpt, mpt1, mpt2, mpt3; |
50944bca | 60 | |
6295157a SP |
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); | |
0e9be4db | 69 | dx = fabs (dx); |
6295157a SP |
70 | for (m = 6; m > 0; m--) |
71 | { | |
72 | if (dx > __atan_xm[m].d) | |
73 | break; | |
74 | } | |
e4d82761 | 75 | } |
6295157a | 76 | mptwoim1.e = 1; |
c2d94018 | 77 | mptwoim1.d[0] = 1; |
50944bca | 78 | |
6295157a SP |
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 | { | |
107a5bf0 | 87 | __add (&__mpone, &mpsm, &mpt1, p); |
6295157a SP |
88 | __mpsqrt (&mpt1, &mpt2, p); |
89 | __add (&mpt2, &mpt2, &mpt1, p); | |
107a5bf0 | 90 | __add (&__mptwo, &mpsm, &mpt2, p); |
6295157a SP |
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]; | |
e4d82761 UD |
97 | } |
98 | ||
6295157a SP |
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 | { | |
c871eccd | 105 | mptwoim1.d[1] -= 2; |
6295157a SP |
106 | __dvd (&mpsm, &mptwoim1, &mpt1, p); |
107 | __mul (&mpsm, &mpt, &mpt2, p); | |
108 | __sub (&mpt1, &mpt2, &mpt, p); | |
e4d82761 | 109 | } |
6295157a SP |
110 | __mul (&mps, &mpt, &mpt1, p); |
111 | __sub (&mps, &mpt1, &mpt, p); | |
50944bca | 112 | |
6295157a SP |
113 | /* Compute Atan(x). */ |
114 | mptwoim1.d[1] = 1 << m; | |
115 | __mul (&mptwoim1, &mpt, y, p); | |
e4d82761 | 116 | } |