Annotation of OpenXM_contrib/gmp/mpfr/atan.c, Revision 1.1.1.1
1.1 ohara 1: /* mpfr_atan -- arc-tangent of a floating-point number
2:
3: Copyright 2001 Free Software Foundation.
4:
5: This file is part of the MPFR Library, and was contributed by Mathieu Dutour.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
8: it under the terms of the GNU Lesser General Public License as published by
9: the Free Software Foundation; either version 2.1 of the License, or (at your
10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Lesser General Public License
18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include <stdio.h>
23: #include <stdlib.h>
24: #include "gmp.h"
25: #include "gmp-impl.h"
26: #include "mpfr.h"
27: #include "mpfr-impl.h"
28:
29: #define CST 2.27 /* CST=1+ln(2.4)/ln(2) */
30: #define CST2 1.45 /* CST2=1/ln(2) */
31:
32: static int mpfr_atan_aux _PROTO((mpfr_ptr, mpz_srcptr, int, int));
33:
34: #undef B
35: #define A
36: #define A1 1
37: #define A2 2
38: #define C
39: #define C1 3
40: #define C2 2
41: #define NO_FACTORIAL
42: #define GENERIC mpfr_atan_aux
43: #include "generic.c"
44: #undef C
45: #undef C1
46: #undef C2
47: #undef A
48: #undef A1
49: #undef A2
50: #undef NO_FACTORIAL
51: #undef GENERIC
52:
53: int
54: mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode)
55: {
56: mpfr_t Pisur2;
57: mpfr_t xp;
58: mpfr_t arctgt;
59:
60: int comparaison, signe, suplement;
61:
62: mpfr_t t_arctan;
63: int i;
64: mpz_t ukz;
65: mpfr_t ukf;
66: mpfr_t sk,Ak;
67: mpz_t square;
68: mpfr_t tmp_arctan;
69: mpfr_t tmp, tmp2;
70: #ifdef DEBUG
71: mpfr_t tst;
72: #endif
73: int twopoweri;
74: int Prec;
75: int prec_x;
76: int prec_arctan;
77: int good = 0;
78: int realprec;
79: int estimated_delta;
80: /* calculation of the floor */
81: mp_exp_t exptol;
82:
83: int N0;
84: int logn;
85:
86: /* Trivial cases */
87: if (MPFR_IS_NAN(x))
88: {
89: MPFR_SET_NAN(arctangent);
90: MPFR_RET_NAN;
91: }
92:
93: if (MPFR_IS_INF(x))
94: {
95: MPFR_CLEAR_FLAGS(arctangent);
96: if (MPFR_SIGN(x) > 0) /* arctan(+inf) = Pi/2 */
97: mpfr_const_pi (arctangent, rnd_mode);
98: else /* arctan(-inf) = -Pi/2 */
99: {
100: if (rnd_mode == GMP_RNDU)
101: rnd_mode = GMP_RNDD;
102: else if (rnd_mode == GMP_RNDD)
103: rnd_mode = GMP_RNDU;
104: mpfr_const_pi (arctangent, rnd_mode);
105: }
106: MPFR_EXP(arctangent)--;
107: return 1; /* inexact */
108: }
109:
110: MPFR_CLEAR_FLAGS(arctangent);
111: if (MPFR_IS_ZERO(x))
112: {
113: mpfr_set_ui (arctangent, 0, GMP_RNDN);
114: return 0; /* exact result */
115: }
116:
117: signe = MPFR_SIGN(x);
118: prec_arctan = MPFR_PREC(arctangent);
119:
120: /* Set x_p=|x| */
121: mpfr_init2(xp, MPFR_PREC(x));
122: mpfr_set(xp, x, rnd_mode);
123: if (signe == -1)
124: MPFR_CHANGE_SIGN(xp);
125:
126: /* Other simple case arctang(-+1)=-+pi/4 */
127: comparaison=mpfr_cmp_ui(xp, 1);
128: if (comparaison == 0) {
129: mpfr_init2(Pisur2, prec_arctan);
130: mpfr_const_pi(Pisur2, rnd_mode);
131: mpfr_div_2ui(arctangent, Pisur2, 2, rnd_mode);
132: if (signe == -1)
133: MPFR_CHANGE_SIGN(arctangent);
134: mpfr_clear(Pisur2);
135: mpfr_clear(xp);
136: return 0; /* Result correct */
137: }
138: if (comparaison > 0)
139: suplement = 2;
140: else
141: suplement = 2-MPFR_EXP(xp);
142:
143: prec_x = _mpfr_ceil_log2 ((double) MPFR_PREC(x) / BITS_PER_MP_LIMB);
144: logn = _mpfr_ceil_log2 ((double) prec_x);
145: if (logn < 2) logn = 2;
146: realprec = prec_arctan + _mpfr_ceil_log2((double) prec_arctan) + 4;
147: mpz_init(ukz);
148: mpz_init(square);
149:
150:
151: while (!good){
152: N0 = _mpfr_ceil_log2((double) realprec + suplement + CST);
153: estimated_delta = 1 + suplement + _mpfr_ceil_log2((double) (3*N0-2));
154: Prec = realprec+estimated_delta;
155:
156: /* Initialisation */
157: mpfr_init2(sk,Prec);
158: mpfr_init2(ukf, Prec);
159: mpfr_init2(t_arctan, Prec);
160: mpfr_init2(tmp_arctan, Prec);
161: mpfr_init2(tmp, Prec);
162: mpfr_init2(tmp2, Prec);
163: mpfr_init2(Ak, Prec);
164: mpfr_init2(arctgt, Prec);
165:
166: #ifdef DEBUG
167: /* Tests */
168: mpfr_init2(tst, realprec);
169: #endif
170:
171: if (comparaison > 0)
172: {
173: mpfr_init2(Pisur2, Prec);
174: mpfr_const_pi(Pisur2, GMP_RNDN);
175: mpfr_div_2ui(Pisur2, Pisur2, 1, GMP_RNDN);
176: mpfr_ui_div(sk, 1, xp, GMP_RNDN);
177: }
178: else
179: mpfr_set(sk, xp, GMP_RNDN);
180:
181: /* Assignation */
182: mpfr_set_ui (tmp_arctan, 0, GMP_RNDN);
183: twopoweri = 1;
184: for(i = 0; i <= N0; i++){
185: mpfr_mul_2ui(tmp, sk, twopoweri, GMP_RNDN);
186: /* Calculation of trunc(tmp) --> mpz */
187: mpfr_trunc (ukf, tmp);
188: exptol = mpfr_get_z_exp (ukz, ukf);
189: if (exptol>0)
190: mpz_mul_2exp (ukz, ukz, exptol);
191: else
192: mpz_tdiv_q_2exp (ukz, ukz, (unsigned long int) (-exptol));
193:
194: /* Calculation of arctan(Ak) */
195: mpz_mul(square, ukz, ukz);
196: mpz_neg(square, square);
197: mpfr_atan_aux(t_arctan, square, 2*twopoweri, N0 - i);
198: mpfr_set_z(Ak, ukz, GMP_RNDN);
199: mpfr_div_2ui(Ak, Ak, twopoweri, GMP_RNDN);
200: mpfr_mul(t_arctan, t_arctan, Ak, GMP_RNDN);
201:
202: /* Addition and iteration */
203: mpfr_add(tmp_arctan, tmp_arctan, t_arctan, GMP_RNDN);
204: if (i<N0)
205: {
206: mpfr_sub(tmp, sk, Ak, GMP_RNDN);
207: mpfr_mul(tmp2, sk, Ak, GMP_RNDN);
208: mpfr_add_ui(tmp2, tmp2, 1, GMP_RNDN);
209: mpfr_div(sk, tmp, tmp2, GMP_RNDN);
210: twopoweri <<= 1;
211: }
212: }
213:
214: if (comparaison > 0)
215: {
216: mpfr_sub(arctgt, Pisur2, tmp_arctan, GMP_RNDN);
217: if (signe == -1)
218: MPFR_CHANGE_SIGN(arctgt);
219: }
220: else
221: {
222: mpfr_set(arctgt, tmp_arctan, GMP_RNDN);
223: if (signe == -1)
224: MPFR_CHANGE_SIGN(arctgt);
225: }
226:
227: #ifdef DEBUG
228: mpfr_set(tst, arctgt, rnd_mode);
229: #endif
230:
231: if (mpfr_can_round(arctgt, realprec, GMP_RNDN, rnd_mode, MPFR_PREC(arctangent)))
232: {
233: mpfr_set(arctangent, arctgt, rnd_mode);
234: good = 1;
235: realprec += 1;
236: }
237: else
238: {
239: realprec += _mpfr_ceil_log2 ((double) realprec);
240: }
241:
242: mpfr_clear(sk);
243: mpfr_clear(ukf);
244: mpfr_clear(t_arctan);
245: mpfr_clear(tmp_arctan);
246: mpfr_clear(tmp);
247: mpfr_clear(tmp2);
248: mpfr_clear(Ak);
249: mpfr_clear(arctgt);
250:
251: #ifdef DEBUG
252: mpfr_clear(tst);
253: #endif
254: if (comparaison > 0)
255: mpfr_clear(Pisur2);
256: }
257:
258: mpfr_clear(xp);
259: mpz_clear(ukz);
260: mpz_clear(square);
261:
262: return 1; /* inexact result */
263: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to atan.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr