Annotation of OpenXM_contrib/gmp/mpfr/asin.c, Revision 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