diff options
| -rw-r--r-- | Makefile.am | 2 | ||||
| -rw-r--r-- | mpfr.h | 2 | ||||
| -rw-r--r-- | mpfr.texi | 8 | ||||
| -rw-r--r-- | sinh_cosh.c | 144 | ||||
| -rw-r--r-- | tests/Makefile.am | 2 | ||||
| -rw-r--r-- | tests/reuse.c | 1 | ||||
| -rw-r--r-- | tests/tsinh_cosh.c | 129 |
7 files changed, 286 insertions, 2 deletions
diff --git a/Makefile.am b/Makefile.am index a01c44723..af023fcac 100644 --- a/Makefile.am +++ b/Makefile.am @@ -7,7 +7,7 @@ include_HEADERS = mpfr.h mpf2mpfr.h lib_LTLIBRARIES = libmpfr.la -libmpfr_la_SOURCES = mpfr.h mpf2mpfr.h mpfr-gmp.h mpfr-impl.h mpfr-longlong.h exceptions.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c modf.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c factorial.c cosh.c sinh.c tanh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c fms.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c signbit.c copysign.c setsign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c free_cache.c si_op.c cmp_ld.c set_ui_2exp.c set_si_2exp.c set_uj.c set_sj.c get_sj.c get_uj.c get_z.c iszero.c cache.c sqr.c int_ceil_log2.c isqrt.c strtofr.c pow_z.c logging.c mulders.c get_f.c round_p.c erfc.c atan2.c subnormal.c const_catalan.c root.c gen_inverse.h sec.c csc.c cot.c eint.c sech.c csch.c coth.c round_near_x.c constant.c abort_prec_max.c stack_interface.c lngamma.c zeta_ui.c set_d64.c get_d64.c jn.c yn.c remquo.c get_patches.c add_d.c sub_d.c d_sub.c mul_d.c div_d.c d_div.c +libmpfr_la_SOURCES = mpfr.h mpf2mpfr.h mpfr-gmp.h mpfr-impl.h mpfr-longlong.h exceptions.c extract.c uceil_exp2.c uceil_log2.c ufloor_log2.c add.c add1.c add_ui.c agm.c clear.c cmp.c cmp_abs.c cmp_si.c cmp_ui.c comparisons.c div_2exp.c div_2si.c div_2ui.c div.c div_ui.c dump.c eq.c exp10.c exp2.c exp3.c exp.c frac.c get_d.c get_exp.c get_str.c init.c inp_str.c isinteger.c isinf.c isnan.c isnum.c const_log2.c log.c modf.c mul_2exp.c mul_2si.c mul_2ui.c mul.c mul_ui.c neg.c next.c out_str.c const_pi.c pow.c pow_si.c pow_ui.c print_raw.c print_rnd_mode.c random2.c random.c reldiff.c round_prec.c set.c setmax.c setmin.c set_d.c set_dfl_prec.c set_exp.c set_rnd.c set_f.c set_prc_raw.c set_prec.c set_q.c set_si.c set_str.c set_str_raw.c set_ui.c set_z.c sqrt.c sqrt_ui.c sub.c sub1.c sub_ui.c rint.c ui_div.c ui_sub.c urandomb.c get_z_exp.c swap.c factorial.c cosh.c sinh.c tanh.c sinh_cosh.c acosh.c asinh.c atanh.c atan.c cmp2.c exp_2.c asin.c const_euler.c cos.c sin.c tan.c fma.c fms.c hypot.c log1p.c expm1.c log2.c log10.c ui_pow.c ui_pow_ui.c minmax.c dim.c signbit.c copysign.c setsign.c gmp_op.c init2.c acos.c sin_cos.c set_nan.c set_inf.c powerof2.c gamma.c set_ld.c get_ld.c cbrt.c volatile.c fits_s.h fits_sshort.c fits_sint.c fits_slong.c fits_u.h fits_ushort.c fits_uint.c fits_ulong.c fits_uintmax.c fits_intmax.c get_si.c get_ui.c zeta.c cmp_d.c erf.c inits.c inits2.c clears.c sgn.c check.c sub1sp.c version.c mpn_exp.c mpfr-gmp.c mp_clz_tab.c sum.c add1sp.c free_cache.c si_op.c cmp_ld.c set_ui_2exp.c set_si_2exp.c set_uj.c set_sj.c get_sj.c get_uj.c get_z.c iszero.c cache.c sqr.c int_ceil_log2.c isqrt.c strtofr.c pow_z.c logging.c mulders.c get_f.c round_p.c erfc.c atan2.c subnormal.c const_catalan.c root.c gen_inverse.h sec.c csc.c cot.c eint.c sech.c csch.c coth.c round_near_x.c constant.c abort_prec_max.c stack_interface.c lngamma.c zeta_ui.c set_d64.c get_d64.c jn.c yn.c remquo.c get_patches.c add_d.c sub_d.c d_sub.c mul_d.c div_d.c d_div.c libmpfr_la_LIBADD = @LIBOBJS@ @@ -505,6 +505,8 @@ __MPFR_DECLSPEC int mpfr_asinh _MPFR_PROTO((mpfr_ptr,mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_cosh _MPFR_PROTO((mpfr_ptr,mpfr_srcptr, mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_sinh _MPFR_PROTO((mpfr_ptr,mpfr_srcptr, mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_tanh _MPFR_PROTO((mpfr_ptr,mpfr_srcptr, mpfr_rnd_t)); +__MPFR_DECLSPEC int mpfr_sinh_cosh _MPFR_PROTO ((mpfr_ptr, mpfr_ptr, + mpfr_srcptr, mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_sech _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t)); __MPFR_DECLSPEC int mpfr_csch _MPFR_PROTO ((mpfr_ptr, mpfr_srcptr,mpfr_rnd_t)); @@ -1497,6 +1497,14 @@ Set @var{rop} to the hyperbolic cosine, sine or tangent of @var{op}, rounded in the direction @var{rnd}. @end deftypefun +@deftypefun int mpfr_sinh_cosh (mpfr_t @var{sop}, mpfr_t @var{cop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) +Set simultaneously @var{sop} to the hyperbolic sine of @var{op} and + @var{cop} to the hyperbolic cosine of @var{op}, +rounded in the direction @var{rnd} with the corresponding precision of +@var{sop} and @var{cop} which must be different variables. +Return 0 iff both results are exact. +@end deftypefun + @deftypefun int mpfr_sech (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) @deftypefunx int mpfr_csch (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) @deftypefunx int mpfr_coth (mpfr_t @var{rop}, mpfr_t @var{op}, mp_rnd_t @var{rnd}) diff --git a/sinh_cosh.c b/sinh_cosh.c new file mode 100644 index 000000000..d9ebf31a1 --- /dev/null +++ b/sinh_cosh.c @@ -0,0 +1,144 @@ +/* mpfr_sinh_cosh -- hyperbolic sine and cosine + +Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007 Free Software Foundation, Inc. +Contributed by the Arenaire and Cacao projects, INRIA. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, +MA 02110-1301, USA. */ + +#define MPFR_NEED_LONGLONG_H +#include "mpfr-impl.h" + + /* The computations are done by + cosh(x) = 1/2 [e^(x)+e^(-x)] + sinh(x) = 1/2 [e^(x)-e^(-x)] + Adapted from mpfr_sinh.c */ + +int +mpfr_sinh_cosh (mpfr_ptr sh, mpfr_ptr ch, mpfr_srcptr xt, mp_rnd_t rnd_mode) +{ + mpfr_t x; + int inexact, inexact_sh, inexact_ch; + + MPFR_ASSERTN (sh != ch); + + MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", xt, xt, rnd_mode), + ("sh[%#R]=%R ch[%#R]=%R inexact=%d", sh, sh, ch, ch, inexact)); + + if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt))) + { + if (MPFR_IS_NAN (xt)) + { + MPFR_SET_NAN (ch); + MPFR_SET_NAN (sh); + MPFR_RET_NAN; + } + else if (MPFR_IS_INF (xt)) + { + MPFR_SET_INF (sh); + MPFR_SET_SAME_SIGN (sh, xt); + MPFR_SET_INF (ch); + MPFR_SET_POS (ch); + MPFR_RET (0); + } + else /* xt is zero */ + { + MPFR_ASSERTD (MPFR_IS_ZERO (xt)); + MPFR_SET_ZERO (sh); /* sinh(0) = 0 */ + MPFR_SET_SAME_SIGN (sh, xt); + return mpfr_set_ui (ch, 1, rnd_mode); /* cosh(0) = 1 */ + } + } + + MPFR_TMP_INIT_ABS (x, xt); + + { + mpfr_t s, c, ti; + mp_exp_t d; + mp_prec_t N; /* Precision of the intermediary variables */ + long int err; /* Precision of error */ + MPFR_ZIV_DECL (loop); + MPFR_SAVE_EXPO_DECL (expo); + MPFR_GROUP_DECL (group); + + MPFR_SAVE_EXPO_MARK (expo); + + /* compute the precision of intermediary variable */ + N = MPFR_PREC (ch); + N = MAX (N, MPFR_PREC (sh)); + N = MAX (N, MPFR_PREC (x)); + /* the optimal number of bits : see algorithms.ps */ + N = N + MPFR_INT_CEIL_LOG2 (N) + 4; + + /* initialise of intermediary variables */ + MPFR_GROUP_INIT_3 (group, N, s, c, ti); + + /* First computation of sinh_cosh */ + MPFR_ZIV_INIT (loop, N); + for (;;) { + /* compute sinh_cosh */ + mpfr_clear_flags (); + mpfr_exp (s, x, GMP_RNDD); /* exp(x) */ + /* exp(x) can overflow! */ + /* BUG/TODO/FIXME: exp can overflow but sinh or cosh may be + representable! */ + if (MPFR_UNLIKELY (mpfr_overflow_p ())) { + inexact_ch = mpfr_overflow (ch, rnd_mode, MPFR_SIGN_POS); + inexact_sh = mpfr_overflow (sh, rnd_mode, MPFR_SIGN (xt)); + MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW); + break; + } + d = MPFR_GET_EXP (s); + mpfr_ui_div (ti, 1, s, GMP_RNDU); /* 1/exp(x) */ + mpfr_add (c, s, ti, GMP_RNDU); /* exp(x) + 1/exp(x) */ + mpfr_sub (s, s, ti, GMP_RNDN); /* exp(x) - 1/exp(x) */ + mpfr_div_2ui (c, c, 1, GMP_RNDN); /* 1/2(exp(x) + 1/exp(x)) */ + mpfr_div_2ui (s, s, 1, GMP_RNDN); /* 1/2(exp(x) - 1/exp(x)) */ + + /* it may be that s is zero (in fact, it can only occur when exp(x)=1, + and thus ti=1 too) */ + if (MPFR_IS_ZERO (s)) + err = N; /* double the precision */ + else + { + /* calculation of the error */ + d = d - MPFR_GET_EXP (s) + 2; + /* error estimate: err = N-(__gmpfr_ceil_log2(1+pow(2,d)));*/ + err = N - (MAX (d, 0) + 1); + if (MPFR_LIKELY (MPFR_CAN_ROUND (s, err, MPFR_PREC (sh), rnd_mode) \ + && (MPFR_CAN_ROUND (c, err, MPFR_PREC (ch), rnd_mode)))) + { + inexact_sh = mpfr_set4 (sh, s, rnd_mode, MPFR_SIGN (xt)); + inexact_ch = mpfr_set (ch, c, rnd_mode); + break; + } + } + /* actualisation of the precision */ + N += err; + MPFR_ZIV_NEXT (loop, N); + MPFR_GROUP_REPREC_3 (group, N, s, c, ti); + } + MPFR_ZIV_FREE (loop); + MPFR_GROUP_CLEAR (group); + MPFR_SAVE_EXPO_FREE (expo); + } + /* now, let's raise the flags if needed */ + inexact = mpfr_check_range (sh, inexact_sh, rnd_mode); + inexact = mpfr_check_range (ch, inexact_ch, rnd_mode) || inexact; + + return inexact; +} diff --git a/tests/Makefile.am b/tests/Makefile.am index 9a7a1d56a..40f92653d 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,6 +1,6 @@ AUTOMAKE_OPTIONS = 1.6 gnu $(top_builddir)/ansi2knr -check_PROGRAMS = tversion tinternals tinits tisqrt tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tadd_d tsub_d td_sub tmul_d tdiv_d td_div tgmpop tsi_op tmul_2exp tfma tfms tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac tmodf texp texp2 texp10 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal tlngamma tlgamma tzeta_ui tget_ld_2exp tget_set_d64 tj0 tj1 tjn ty0 ty1 tyn tremquo tl2b +check_PROGRAMS = tversion tinternals tinits tisqrt tsgn tcheck tisnan texceptions tset_exp tset tabs tset_d tset_f tset_q tset_si tset_str tset_z tset_ld tset_sj tswap tcopysign tcmp tcmp2 tcmpabs tcmp_d tcmp_ld tcomparisons teq tadd tsub tmul tdiv tsub1sp tadd1sp tadd_ui tsub_ui tcmp_ui tdiv_ui tmul_ui tsqrt_ui tui_div tui_sub tadd_d tsub_d td_sub tmul_d tdiv_d td_div tgmpop tsi_op tmul_2exp tfma tfms tsum tdim tminmax tnext tfits tget_d tget_d_2exp tget_z tget_str tget_sj tout_str tinp_str toutimpl tcan_round tround_prec tsqrt tconst_log2 tconst_pi tconst_euler trandom ttrunc trint tfrac tmodf texp texp2 texp10 texpm1 tlog tlog2 tlog10 tlog1p tpow tui_pow tpow3 tcosh tsinh ttanh tsinh_cosh tacosh tasinh tatanh thyperbolic tasin tacos tcos tatan tsin ttan tsin_cos tagm thypot tfactorial tgamma terf tcbrt tzeta mpf_compat mpfr_compat reuse tsqr tstrtofr tpow_z tget_f tconst_catalan troot tsec tcsc tcot teint tcoth tcsch tsech tstckintc tsubnormal tlngamma tlgamma tzeta_ui tget_ld_2exp tget_set_d64 tj0 tj1 tjn ty0 ty1 tyn tremquo tl2b EXTRA_DIST = tgeneric.c tgeneric_ui.c mpf_compat.h inp_str.data diff --git a/tests/reuse.c b/tests/reuse.c index f6fe715a0..df6de4716 100644 --- a/tests/reuse.c +++ b/tests/reuse.c @@ -604,6 +604,7 @@ main (void) test3a (mpfr_modf, "mpfr_modf", p, rnd); test3a (mpfr_sin_cos, "mpfr_sin_cos", p, rnd); + test3a (mpfr_sinh_cosh, "mpfr_sinh_cosh", p, rnd); } tests_end_mpfr (); diff --git a/tests/tsinh_cosh.c b/tests/tsinh_cosh.c new file mode 100644 index 000000000..a7a747039 --- /dev/null +++ b/tests/tsinh_cosh.c @@ -0,0 +1,129 @@ +/* Test file for mpfr_sinh_cosh. + +Copyright 2007 Free Software Foundation, Inc. +Contributed by the Arenaire and Cacao projects, INRIA. + +This file is part of the MPFR Library. + +The MPFR Library is free software; you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 2.1 of the License, or (at your +option) any later version. + +The MPFR Library is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with the MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, +MA 02110-1301, USA. */ + +#include <stdio.h> +#include <stdlib.h> + +#include "mpfr-test.h" + +static void +failed (mpfr_t x, mpfr_t esh, mpfr_t gsh, mpfr_t ech, mpfr_t gch) +{ + fprintf (stderr, "error : mpfr_sinh_cosh (x) x = "); + mpfr_out_str (stderr, 10, 0, x, GMP_RNDD); + fprintf (stderr, "\nsinh(x) expected "); + mpfr_out_str (stderr, 10, 0, esh, GMP_RNDD); + fprintf (stderr, "\n got "); + mpfr_out_str (stderr, 10, 0, gsh, GMP_RNDD); + fprintf (stderr, "\ncosh(x) expected "); + mpfr_out_str (stderr, 10, 0, ech, GMP_RNDD); + fprintf (stderr, "\n got "); + mpfr_out_str (stderr, 10, 0, gch, GMP_RNDD); + fprintf (stderr, "\n"); + + mpfr_clears (x, esh, gsh, ech, gch, (void*)0); + exit (-1); +} + +/* check against sinh, cosh */ +static void +check (mpfr_t x, mp_rnd_t rnd) +{ + mpfr_t s, c, sx, cx; + int isc, is, ic; + + mpfr_inits2 (MPFR_PREC(x), s, c, sx, cx, (void *)0); + + isc = mpfr_sinh_cosh (sx, cx, x, rnd); + is = mpfr_sinh (s, x, rnd); + ic = mpfr_cosh (c, x, rnd); + + if (!mpfr_equal_p (s, sx) || !mpfr_equal_p (c, cx)) + failed (x, s, sx, c, cx); + MPFR_ASSERTN (isc = is || ic); + + mpfr_clears (s, c, sx, cx, (void *)0); +} + +static void +check_nans (void) +{ + mpfr_t x, sh, ch; + + mpfr_init2 (x, 123); + mpfr_init2 (sh, 123); + mpfr_init2 (ch, 123); + + /* nan */ + mpfr_set_nan (x); + mpfr_sinh_cosh (sh, ch, x, GMP_RNDN); + MPFR_ASSERTN (mpfr_nan_p (sh)); + MPFR_ASSERTN (mpfr_nan_p (ch)); + + /* +inf */ + mpfr_set_inf (x, 1); + mpfr_sinh_cosh (sh, ch, x, GMP_RNDN); + MPFR_ASSERTN (mpfr_inf_p (sh)); + MPFR_ASSERTN (mpfr_sgn (sh) > 0); + MPFR_ASSERTN (mpfr_inf_p (ch)); + MPFR_ASSERTN (mpfr_sgn (ch) > 0); + + /* -inf */ + mpfr_set_inf (x, -1); + mpfr_sinh_cosh (sh, ch, x, GMP_RNDN); + MPFR_ASSERTN (mpfr_inf_p (sh)); + MPFR_ASSERTN (mpfr_sgn (sh) < 0); + MPFR_ASSERTN (mpfr_inf_p (ch)); + MPFR_ASSERTN (mpfr_sgn (ch) > 0); + + mpfr_clear (x); + mpfr_clear (sh); + mpfr_clear (ch); +} + +int +main (int argc, char *argv[]) +{ + int i; + mpfr_t x; + tests_start_mpfr (); + + check_nans (); + + /* check against values given by sinh(x), cosh(x) */ + mpfr_init2 (x, 53); + mpfr_set_str (x, "FEDCBA987654321p-48", 16, GMP_RNDN); + for (i = 0; i < 10; ++i) + { + /* x = i - x / 2 : boggle sign and bits */ + mpfr_ui_sub (x, i, x, GMP_RNDD); + mpfr_div_2ui (x, x, 2, GMP_RNDD); + + check (x, GMP_RNDN); + check (x, GMP_RNDU); + check (x, GMP_RNDD); + } + mpfr_clear (x); + + tests_end_mpfr (); + return 0; +} |