Return to acos.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

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: }