Annotation of OpenXM_contrib/gmp/mpfr/acos.c, Revision 1.1.1.1
1.1 ohara 1: /* mpfr_acos -- arc-cosinus 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: int
30: mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode)
31: {
32: mpfr_t xp;
33: mpfr_t arcc;
34:
35: int signe, suplement;
36:
37: mpfr_t tmp;
38: int Prec;
39: int prec_acos;
40: int good = 0;
41: int realprec;
42: int compared;
43: int inexact = 0;
44:
45: /* Trivial cases */
46: if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
47: {
48: MPFR_SET_NAN(acos);
49: MPFR_RET_NAN;
50: }
51:
52: /* Set x_p=|x| */
53: signe = MPFR_SIGN(x);
54: mpfr_init2 (xp, MPFR_PREC(x));
55: mpfr_abs (xp, x, rnd_mode);
56:
57: compared = mpfr_cmp_ui (xp, 1);
58:
59: if (compared > 0) /* acos(x) = NaN for x > 1 */
60: {
61: mpfr_clear (xp);
62: MPFR_SET_NAN(acos);
63: MPFR_RET_NAN;
64: }
65:
66: if (compared == 0)
67: {
68: mpfr_clear (xp);
69: if (signe > 0) /* acos(+1) = 0 */
70: return mpfr_set_ui (acos, 0, rnd_mode);
71: else /* acos(-1) = Pi */
72: {
73: mpfr_const_pi (acos, rnd_mode);
74: return 1; /* inexact */
75: }
76: }
77:
78: if (MPFR_IS_ZERO(x)) /* acos(0)=Pi/2 */
79: {
80: mpfr_clear (xp);
81: mpfr_const_pi (acos, rnd_mode);
82: MPFR_EXP(acos)--;
83: return 1; /* inexact */
84: }
85:
86: prec_acos = MPFR_PREC(acos);
87: mpfr_ui_sub (xp, 1, xp, GMP_RNDD);
88:
89: if (signe > 0)
90: suplement = 2 - 2 * MPFR_EXP(xp);
91: else
92: suplement = 2 - MPFR_EXP(xp);
93:
94: realprec = prec_acos + 10;
95:
96: while (!good)
97: {
98: Prec = realprec+suplement;
99:
100: /* Initialisation */
101: mpfr_init2 (tmp, Prec);
102: mpfr_init2 (arcc, Prec);
103:
104: mpfr_mul (tmp, x, x, GMP_RNDN);
105: mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN);
106: mpfr_sqrt (tmp, tmp, GMP_RNDN);
107: mpfr_div (tmp, x, tmp, GMP_RNDN);
108: mpfr_atan (arcc, tmp, GMP_RNDN);
109: mpfr_const_pi (tmp, GMP_RNDN);
110: mpfr_div_2ui (tmp, tmp, 1, GMP_RNDN);
111: mpfr_sub (arcc, tmp, arcc, GMP_RNDN);
112:
113: if (mpfr_can_round (arcc, realprec, GMP_RNDN, rnd_mode, MPFR_PREC(acos)))
114: {
115: inexact = mpfr_set (acos, arcc, rnd_mode);
116: good = 1;
117: }
118: else
119: realprec += _mpfr_ceil_log2 ((double) realprec);
120:
121: mpfr_clear (tmp);
122: mpfr_clear (arcc);
123: }
124:
125: mpfr_clear (xp);
126: return inexact;
127: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>