1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2021 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 <https://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
43void
44SECTION
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}
117