Return to sinh.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

Annotation of OpenXM_contrib/gmp/mpfr/sinh.c, Revision 1.1.1.1

1.1       ohara       1: /* mpfr_sinh -- hyperbolic sine
                      2:
                      3: Copyright 2001, 2002 Free Software Foundation, Inc.
                      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:
                     23: #include "gmp.h"
                     24: #include "gmp-impl.h"
                     25: #include "mpfr.h"
                     26: #include "mpfr-impl.h"
                     27:
                     28:  /* The computation of sinh is done by
                     29:
                     30:     sinh(x) = 1/2 [e^(x)-e^(-x)]
                     31:  */
                     32:
                     33: int
                     34: mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
                     35: {
                     36:     /****** Declarations ******/
                     37:     mpfr_t x;
                     38:     mp_prec_t Nxt = MPFR_PREC(xt);
                     39:     int flag_neg=0, inexact =0;
                     40:
                     41:     if (MPFR_IS_NAN(xt))
                     42:       {
                     43:         MPFR_SET_NAN(y);
                     44:         MPFR_RET_NAN;
                     45:       }
                     46:     MPFR_CLEAR_NAN(y);
                     47:
                     48:     if (MPFR_IS_INF(xt))
                     49:       {
                     50:         MPFR_SET_INF(y);
                     51:         MPFR_SET_SAME_SIGN(y, xt);
                     52:         MPFR_RET(0);
                     53:       }
                     54:
                     55:     MPFR_CLEAR_INF(y);
                     56:
                     57:     if (MPFR_IS_ZERO(xt))
                     58:       {
                     59:         MPFR_SET_ZERO(y);   /* sinh(0) = 0 */
                     60:         MPFR_SET_SAME_SIGN(y, xt);
                     61:         MPFR_RET(0);
                     62:       }
                     63:
                     64:     mpfr_init2 (x, Nxt);
                     65:     mpfr_set (x, xt, GMP_RNDN);
                     66:
                     67:     if(MPFR_SIGN(x)<0)
                     68:       {
                     69:         MPFR_CHANGE_SIGN(x);
                     70:         flag_neg=1;
                     71:       }
                     72:
                     73:     /* General case */
                     74:     {
                     75:     /* Declaration of the intermediary variable */
                     76:       mpfr_t t, te, ti;
                     77:       int d;
                     78:
                     79:       /* Declaration of the size variable */
                     80:       mp_prec_t Nx = Nxt;   /* Precision of input variable */
                     81:       mp_prec_t Ny = MPFR_PREC(y);   /* Precision of input variable */
                     82:
                     83:       mp_prec_t Nt;   /* Precision of the intermediary variable */
                     84:       long int err;  /* Precision of error */
                     85:
                     86:       /* compute the precision of intermediary variable */
                     87:       Nt = MAX(Nx, Ny);
                     88:       /* the optimal number of bits : see algorithms.ps */
                     89:       Nt = Nt + _mpfr_ceil_log2 (5) + _mpfr_ceil_log2 (Nt);
                     90:
                     91:       /* initialise of intermediary    variable */
                     92:       mpfr_init (t);
                     93:       mpfr_init (te);
                     94:       mpfr_init (ti);
                     95:
                     96:       /* First computation of sinh */
                     97:       do {
                     98:
                     99:        /* reactualisation of the precision */
                    100:
                    101:        mpfr_set_prec (t, Nt);
                    102:        mpfr_set_prec (te, Nt);
                    103:        mpfr_set_prec (ti, Nt);
                    104:
                    105:        /* compute sinh */
                    106:        mpfr_exp (te, x, GMP_RNDD);        /* exp(x) */
                    107:        mpfr_ui_div (ti, 1, te, GMP_RNDU); /* 1/exp(x) */
                    108:         mpfr_sub (t, te, ti, GMP_RNDN);    /* exp(x) - 1/exp(x) */
                    109:        mpfr_div_2ui (t, t, 1, GMP_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */
                    110:
                    111:         /* it may be that t is zero (in fact, it can only occur when te=1,
                    112:            and thus ti=1 too) */
                    113:
                    114:         if (MPFR_IS_ZERO(t))
                    115:           err = -1;
                    116:         else
                    117:           {
                    118:             /* calculation of the error */
                    119:             d = MPFR_EXP(te) - MPFR_EXP(t) + 2;
                    120:
                    121:             /* estimation of the error */
                    122:             /* err = Nt-(_mpfr_ceil_log2(1+pow(2,d)));*/
                    123:             err = Nt - (MAX(d,0) + 1);
                    124:           }
                    125:
                    126:        /* actualisation of the precision */
                    127:         Nt += 10;
                    128:
                    129:       } while ((err < 0) || !mpfr_can_round(t, err, GMP_RNDN, rnd_mode, Ny));
                    130:
                    131:       if (flag_neg == 1)
                    132:           MPFR_CHANGE_SIGN(t);
                    133:
                    134:       inexact = mpfr_set (y, t, rnd_mode);
                    135:       mpfr_clear (t);
                    136:       mpfr_clear (ti);
                    137:       mpfr_clear (te);
                    138:     }
                    139:     mpfr_clear (x);
                    140:
                    141:     return inexact;
                    142: }