Annotation of OpenXM_contrib/gmp/mpfr/tanh.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_tanh -- hyperbolic tangent
! 2:
! 3: Copyright 2001, 2002 Free Software Foundation.
! 4:
! 5: This file is part of the MPFR Library.
! 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 "gmp.h"
! 23: #include "gmp-impl.h"
! 24: #include "mpfr.h"
! 25: #include "mpfr-impl.h"
! 26:
! 27: /* The computation of cosh is done by
! 28:
! 29: tanh= [e^(x)^2-1]/+[e^(x)^2+1]
! 30: */
! 31:
! 32: int
! 33: #if __STDC__
! 34: mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mp_rnd_t rnd_mode)
! 35: #else
! 36: mpfr_tanh (y, xt, rnd_mode)
! 37: mpfr_ptr y;
! 38: mpfr_srcptr xt;
! 39: mp_rnd_t rnd_mode;
! 40: #endif
! 41: {
! 42:
! 43: /****** Declaration ******/
! 44: mpfr_t x;
! 45: mp_prec_t Nxt = MPFR_PREC(xt);
! 46: int flag_neg=0, inexact=0;
! 47:
! 48: /* Special value checking */
! 49:
! 50: if (MPFR_IS_NAN(xt))
! 51: {
! 52: MPFR_SET_NAN(y);
! 53: MPFR_RET_NAN;
! 54: }
! 55: MPFR_CLEAR_NAN(y);
! 56:
! 57: if (MPFR_IS_INF(xt))
! 58: {
! 59: if (MPFR_SIGN(xt) > 0)
! 60: return mpfr_set_si(y,1,rnd_mode); /* tanh(inf) = 1 */
! 61: else
! 62: return mpfr_set_si(y,-1,rnd_mode); /* tanh(-inf) = -1 */
! 63: }
! 64: MPFR_CLEAR_INF(y);
! 65:
! 66: /* tanh(0) = 0 */
! 67: if (MPFR_IS_ZERO(xt))
! 68: {
! 69: MPFR_SET_ZERO(y);
! 70: MPFR_SET_SAME_SIGN(y,xt);
! 71: MPFR_RET(0);
! 72: }
! 73:
! 74: mpfr_init2(x,Nxt);
! 75: mpfr_set(x,xt,GMP_RNDN);
! 76:
! 77: if (MPFR_SIGN(x) < 0)
! 78: {
! 79: MPFR_CHANGE_SIGN(x);
! 80: flag_neg=1;
! 81: }
! 82:
! 83: /* General case */
! 84: {
! 85: /* Declaration of the intermediary variable */
! 86: mpfr_t t, te, ta,tb;
! 87: int d;
! 88:
! 89: /* Declaration of the size variable */
! 90: mp_prec_t Nx = Nxt; /* Precision of input variable */
! 91: mp_prec_t Ny = MPFR_PREC(y); /* Precision of input variable */
! 92:
! 93: mp_prec_t Nt; /* Precision of the intermediary variable */
! 94: long int err; /* Precision of error */
! 95:
! 96: /* compute the precision of intermediary variable */
! 97: Nt=MAX(Nx,Ny);
! 98: /* the optimal number of bits : see algorithms.ps */
! 99: Nt = Nt+_mpfr_ceil_log2(9)+_mpfr_ceil_log2(Nt);
! 100:
! 101: /* initialise of intermediary variable */
! 102: mpfr_init(t);
! 103: mpfr_init(te);
! 104: mpfr_init(ta);
! 105: mpfr_init(tb);
! 106:
! 107:
! 108: /* First computation of cosh */
! 109: do {
! 110:
! 111: /* reactualisation of the precision */
! 112: mpfr_set_prec(t,Nt);
! 113: mpfr_set_prec(te,Nt);
! 114: mpfr_set_prec(ta,Nt);
! 115: mpfr_set_prec(tb,Nt);
! 116:
! 117: /* compute tanh */
! 118: mpfr_mul_2ui(te,x,1,GMP_RNDN); /* 2x */
! 119: mpfr_exp(te,te,GMP_RNDN); /* exp(2x) */
! 120: mpfr_add_ui(ta,te,1,GMP_RNDD); /* exp(2x) + 1*/
! 121: mpfr_sub_ui(tb,te,1,GMP_RNDU); /* exp(2x) - 1*/
! 122: mpfr_div(t,tb,ta,GMP_RNDN); /* (exp(2x)-1)/(exp(2x)+1)*/
! 123:
! 124:
! 125: /* calculation of the error*/
! 126: d = MPFR_EXP(te)-MPFR_EXP(t);
! 127:
! 128: /* estimation of the error */
! 129: /*err = Nt-(_mpfr_ceil_log2(7+pow(2,d+1)));*/
! 130: err = Nt-(MAX(d+1,3)+1);
! 131:
! 132: /* actualisation of the precision */
! 133: Nt += 10;
! 134:
! 135: } while ((err <0)||!mpfr_can_round(t,err,GMP_RNDN,rnd_mode,Ny));
! 136:
! 137: if (flag_neg==1)
! 138: MPFR_CHANGE_SIGN(t);
! 139:
! 140: inexact = mpfr_set(y,t,rnd_mode);
! 141: mpfr_clear(t);
! 142: mpfr_clear(te);
! 143: mpfr_clear(ta);
! 144: mpfr_clear(tb);
! 145: }
! 146: mpfr_clear(x);
! 147: return inexact;
! 148:
! 149: }
! 150:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>