Annotation of OpenXM_contrib/gmp/mpfr/asin.c, Revision 1.1.1.1
1.1 ohara 1: /* mpfr_asin -- arc-sinus 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_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
31: {
32: mpfr_t xp;
33: mpfr_t arcs;
34:
35: int signe, suplement;
36:
37: mpfr_t tmp;
38: int Prec;
39: int prec_asin;
40: int good = 0;
41: int realprec;
42: int estimated_delta;
43: int compared;
44:
45: /* Trivial cases */
46: if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
47: {
48: MPFR_SET_NAN(asin);
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_set (xp, x, rnd_mode);
56: if (signe == -1)
57: MPFR_CHANGE_SIGN(xp);
58:
59: compared = mpfr_cmp_ui (xp, 1);
60:
61: if (compared > 0) /* asin(x) = NaN for |x| > 1 */
62: {
63: MPFR_SET_NAN(asin);
64: mpfr_clear (xp);
65: MPFR_RET_NAN;
66: }
67:
68: if (compared == 0) /* x = 1 or x = -1 */
69: {
70: if (signe > 0) /* asin(+1) = Pi/2 */
71: mpfr_const_pi (asin, rnd_mode);
72: else /* asin(-1) = -Pi/2 */
73: {
74: if (rnd_mode == GMP_RNDU)
75: rnd_mode = GMP_RNDD;
76: else if (rnd_mode == GMP_RNDD)
77: rnd_mode = GMP_RNDU;
78: mpfr_const_pi (asin, rnd_mode);
79: mpfr_neg (asin, asin, rnd_mode);
80: }
81: MPFR_EXP(asin)--;
82: mpfr_clear (xp);
83: return 1; /* inexact */
84: }
85:
86: if (MPFR_IS_ZERO(x)) /* x = 0 */
87: {
88: mpfr_set_ui (asin, 0, GMP_RNDN);
89: mpfr_clear(xp);
90: return 0; /* exact result */
91: }
92:
93: prec_asin = MPFR_PREC(asin);
94: mpfr_ui_sub (xp, 1, xp, GMP_RNDD);
95:
96: suplement = 2 - MPFR_EXP(xp);
97: #ifdef DEBUG
98: printf("suplement=%d\n", suplement);
99: #endif
100: realprec = prec_asin + 10;
101:
102: while (!good)
103: {
104: estimated_delta = 1 + suplement;
105: Prec = realprec+estimated_delta;
106:
107: /* Initialisation */
108: mpfr_init2 (tmp, Prec);
109: mpfr_init2 (arcs, Prec);
110:
111: #ifdef DEBUG
112: printf("Prec=%d\n", Prec);
113: printf(" x=");
114: mpfr_out_str (stdout, 2, 0, x, GMP_RNDN);
115: printf ("\n");
116: #endif
117: mpfr_mul (tmp, x, x, GMP_RNDN);
118: #ifdef DEBUG
119: printf(" x^2=");
120: mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
121: printf ("\n");
122: #endif
123: mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN);
124: #ifdef DEBUG
125: printf(" 1-x^2=");
126: mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
127: printf ("\n");
128: printf("10: 1-x^2=");
129: mpfr_out_str (stdout, 10, 0, tmp, GMP_RNDN);
130: printf ("\n");
131: #endif
132: mpfr_sqrt (tmp, tmp, GMP_RNDN);
133: #ifdef DEBUG
134: printf(" sqrt(1-x^2)=");
135: mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
136: printf ("\n");
137: printf("10: sqrt(1-x^2)=");
138: mpfr_out_str (stdout, 10, 0, tmp, GMP_RNDN);
139: printf ("\n");
140: #endif
141: mpfr_div (tmp, x, tmp, GMP_RNDN);
142: #ifdef DEBUG
143: printf("x/sqrt(1-x^2)=");
144: mpfr_out_str (stdout, 2, 0, tmp, GMP_RNDN);
145: printf ("\n");
146: #endif
147: mpfr_atan (arcs, tmp, GMP_RNDN);
148: #ifdef DEBUG
149: printf("atan(x/..x^2)=");
150: mpfr_out_str (stdout, 2, 0, arcs, GMP_RNDN);
151: printf ("\n");
152: #endif
153: if (mpfr_can_round (arcs, realprec, GMP_RNDN, rnd_mode, MPFR_PREC(asin)))
154: {
155: mpfr_set (asin, arcs, rnd_mode);
156: #ifdef DEBUG
157: printf("asin =");
158: mpfr_out_str (stdout, 2, prec_asin, asin, GMP_RNDN);
159: printf ("\n");
160: #endif
161: good = 1;
162: }
163: else
164: {
165: realprec += _mpfr_ceil_log2 ((double) realprec);
166: #ifdef DEBUG
167: printf("RETRY\n");
168: #endif
169: }
170: mpfr_clear (tmp);
171: mpfr_clear (arcs);
172: }
173:
174: mpfr_clear (xp);
175:
176: return 1; /* inexact result */
177: }
178:
179:
180:
181:
182:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to asin.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr