From 5f8b33f52f38f11227612f45c3adb56e0cd3e96e Mon Sep 17 00:00:00 2001 From: zimmerma Date: Tue, 3 Apr 2012 15:40:49 +0000 Subject: added new function mpc_log10 git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1149 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- NEWS | 1 + doc/mpc.texi | 5 +- src/Makefile.am | 2 +- src/log10.c | 248 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/mpc.h | 1 + tests/Makefile.am | 8 +- tests/log.dat | 1 + tests/log10.dat | 172 +++++++++++++++++++++++++++++++++++++ tests/tlog10.c | 37 ++++++++ 9 files changed, 468 insertions(+), 7 deletions(-) create mode 100644 src/log10.c create mode 100644 tests/log10.dat create mode 100644 tests/tlog10.c diff --git a/NEWS b/NEWS index 980f19e..0b3887a 100644 --- a/NEWS +++ b/NEWS @@ -4,6 +4,7 @@ Changes in version 1.0: - 0^0, which returned (NaN,NaN) previously, now returns (1,+0). - removed compatibility with K&R compilers, untestable due to lack of such compilers + - New function mpc_log10. - Bug fixes: - div and norm now return a value indicating the effective rounding direction, as the other functions. diff --git a/doc/mpc.texi b/doc/mpc.texi index acb4576..a724877 100644 --- a/doc/mpc.texi +++ b/doc/mpc.texi @@ -959,11 +959,12 @@ rounded according to @var{rnd} with the precision of @var{rop}. @end deftypefun @deftypefun int mpc_log (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) -Set @var{rop} to the logarithm of @var{op}, +@deftypefunx int mpc_log10 (mpc_t @var{rop}, mpc_t @var{op}, mpc_rnd_t @var{rnd}) +Set @var{rop} to the natural and base-10 logarithm of @var{op} respectively, rounded according to @var{rnd} with the precision of @var{rop}. The principal branch is chosen, with the branch cut on the negative real axis, so that the imaginary part of the result lies in -@math{]-\pi , \pi]}. +@math{]-\pi , \pi]} and @math{]-\pi/log(10) , \pi/log(10)]} respectively. @end deftypefun diff --git a/src/Makefile.am b/src/Makefile.am index 0f0ddd6..5cea5c2 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -23,7 +23,7 @@ libmpc_la_SOURCES = mpc-impl.h abs.c acos.c acosh.c add.c add_fr.c \ add_si.c add_ui.c arg.c asin.c asinh.c atan.c atanh.c clear.c cmp.c \ cmp_si_si.c conj.c cos.c cosh.c div_2exp.c div.c div_fr.c div_ui.c exp.c \ fma.c fr_div.c fr_sub.c get_prec2.c get_prec.c get_version.c get_x.c \ - imag.c init2.c init3.c inp_str.c log.c mem.c mul_2exp.c mul.c \ + imag.c init2.c init3.c inp_str.c log.c log10.c mem.c mul_2exp.c mul.c \ mul_fr.c mul_i.c mul_si.c mul_ui.c neg.c norm.c out_str.c pow.c pow_fr.c \ pow_ld.c pow_d.c pow_si.c pow_ui.c pow_z.c proj.c real.c urandom.c set.c \ set_prec.c set_str.c set_x.c set_x_x.c sin.c sin_cos.c sinh.c sqr.c \ diff --git a/src/log10.c b/src/log10.c new file mode 100644 index 0000000..88dafd6 --- /dev/null +++ b/src/log10.c @@ -0,0 +1,248 @@ +/* mpc_log10 -- Take the base-10 logarithm of a complex number. + +Copyright (C) 2012 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#include "mpc-impl.h" + +/* Auxiliary functions which implement Ziv's strategy for special cases. + if flag = 0: compute only real part + if flag = 1: compute only imaginary + Exact cases should be dealt with separately. */ +static int +mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int flag, int nb) +{ + mp_prec_t prec = (MPFR_PREC_MIN > 4) ? MPFR_PREC_MIN : 4; + mpc_t tmp; + mpfr_t log10; + int ok = 0, ret; + + prec = mpfr_get_prec ((flag == 0) ? mpc_realref (rop) : mpc_imagref (rop)); + prec += 10; + mpc_init2 (tmp, prec); + mpfr_init2 (log10, prec); + while (ok == 0) + { + mpfr_set_ui (log10, 10, GMP_RNDN); /* exact since prec >= 4 */ + mpfr_log (log10, log10, GMP_RNDN); + /* In each case we have two roundings, thus the final value is + x * (1+u)^2 where x is the exact value, and |u| <= 2^(-prec-1). + Thus the error is always less than 3 ulps. */ + switch (nb) + { + case 0: /* imag <- atan2(y/x) */ + mpfr_atan2 (mpc_imagref (tmp), mpc_imagref (op), mpc_realref (op), + MPC_RND_IM (rnd)); + mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN); + ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDZ, + GMP_RNDZ, MPC_PREC_IM(rop) + + MPC_RND_IM (rnd) == GMP_RNDN); + if (ok) + ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp), + MPC_RND_IM (rnd)); + break; + case 1: /* real <- log(x) */ + mpfr_log (mpc_realref (tmp), mpc_realref (op), MPC_RND_RE (rnd)); + mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN); + ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDZ, + GMP_RNDZ, MPC_PREC_RE(rop) + + MPC_RND_RE (rnd) == GMP_RNDN); + if (ok) + ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp), + MPC_RND_RE (rnd)); + break; + case 2: /* imag <- pi */ + mpfr_const_pi (mpc_imagref (tmp), MPC_RND_IM (rnd)); + mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, GMP_RNDN); + ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, GMP_RNDZ, + GMP_RNDZ, MPC_PREC_IM(rop) + + MPC_RND_IM (rnd) == GMP_RNDN); + if (ok) + ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp), + MPC_RND_IM (rnd)); + break; + case 3: /* real <- log(y) */ + mpfr_log (mpc_realref (tmp), mpc_imagref (op), MPC_RND_RE (rnd)); + mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, GMP_RNDN); + ok = mpfr_can_round (mpc_realref (tmp), prec - 2, GMP_RNDZ, + GMP_RNDZ, MPC_PREC_RE(rop) + + MPC_RND_RE (rnd) == GMP_RNDN); + if (ok) + ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp), + MPC_RND_RE (rnd)); + break; + } + prec += prec / 2; + mpc_set_prec (tmp, prec); + mpfr_set_prec (log10, prec); + } + mpc_clear (tmp); + mpfr_clear (log10); + return ret; +} + +int +mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) +{ + int ok = 0, loops = 0, re_cmp, im_cmp, inex_re, inex_im, negative_zero; + mpfr_t w; + mpfr_prec_t prec; + mpfr_rnd_t rnd_im; + mpc_t ww; + mpc_rnd_t invrnd; + + /* special values: NaN and infinities: same as mpc_log */ + if (!mpc_fin_p (op)) /* real or imaginary parts are NaN or Inf */ + { + if (mpfr_nan_p (mpc_realref (op))) + { + if (mpfr_inf_p (mpc_imagref (op))) + /* (NaN, Inf) -> (+Inf, NaN) */ + mpfr_set_inf (mpc_realref (rop), +1); + else + /* (NaN, xxx) -> (NaN, NaN) */ + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); + inex_im = 0; /* Inf/NaN is exact */ + } + else if (mpfr_nan_p (mpc_imagref (op))) + { + if (mpfr_inf_p (mpc_realref (op))) + /* (Inf, NaN) -> (+Inf, NaN) */ + mpfr_set_inf (mpc_realref (rop), +1); + else + /* (xxx, NaN) -> (NaN, NaN) */ + mpfr_set_nan (mpc_realref (rop)); + mpfr_set_nan (mpc_imagref (rop)); + inex_im = 0; /* Inf/NaN is exact */ + } + else /* We have an infinity in at least one part. */ + { + /* (+Inf, y) -> (+Inf, 0) for finite positive-signed y */ + if (mpfr_inf_p (mpc_realref (op)) && mpfr_signbit (mpc_realref (op)) + == 0 && mpfr_number_p (mpc_imagref (op))) + inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), + mpc_realref (op), MPC_RND_IM (rnd)); + else + /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx)) + (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */ + inex_im = mpc_log10_aux (rop, op, rnd, 1, 0); + mpfr_set_inf (mpc_realref (rop), +1); + } + return MPC_INEX(0, inex_im); + } + + /* special cases: real and purely imaginary numbers */ + re_cmp = mpfr_cmp_ui (mpc_realref (op), 0); + im_cmp = mpfr_cmp_ui (mpc_imagref (op), 0); + if (im_cmp == 0) /* Im(op) = 0 */ + { + if (re_cmp == 0) /* Re(op) = 0 */ + { + if (mpfr_signbit (mpc_realref (op)) == 0) + inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), + mpc_realref (op), MPC_RND_IM (rnd)); + else + inex_im = mpc_log10_aux (rop, op, rnd, 1, 0); + mpfr_set_inf (mpc_realref (rop), -1); + inex_re = 0; /* -Inf is exact */ + } + else if (re_cmp > 0) + { + inex_re = mpfr_log10 (mpc_realref (rop), mpc_realref (op), + MPC_RND_RE (rnd)); + inex_im = mpfr_set (mpc_imagref (rop), mpc_imagref (op), + MPC_RND_IM (rnd)); + } + else /* log10(x + 0*i) for negative x */ + { /* op = x + 0*i; let w = -x = |x| */ + negative_zero = mpfr_signbit (mpc_imagref (op)); + if (negative_zero) + rnd_im = INV_RND (MPC_RND_IM (rnd)); + else + rnd_im = MPC_RND_IM (rnd); + ww->re[0] = *mpc_realref (op); + MPFR_CHANGE_SIGN (ww->re); + ww->im[0] = *mpc_imagref (op); + if (mpfr_cmp_ui (ww->re, 1) == 0) + inex_re = mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); + else + inex_re = mpc_log10_aux (rop, ww, rnd, 0, 1); + inex_im = mpc_log10_aux (rop, op, RNDC(0,rnd_im), 1, 2); + if (negative_zero) + { + mpc_conj (rop, rop, MPC_RNDNN); + inex_im = -inex_im; + } + } + return MPC_INEX(inex_re, inex_im); + } + else if (re_cmp == 0) + { + if (im_cmp > 0) + { + inex_re = mpc_log10_aux (rop, op, rnd, 0, 3); + inex_im = mpc_log10_aux (rop, op, rnd, 1, 2); + /* division by 2 does not change the ternary flag */ + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); + } + else + { + ww->re[0] = *mpc_realref (op); + ww->im[0] = *mpc_imagref (op); + MPFR_CHANGE_SIGN (ww->im); + inex_re = mpc_log10_aux (rop, ww, rnd, 0, 3); + invrnd = RNDC(0, INV_RND (MPC_RND_IM (rnd))); + inex_im = mpc_log10_aux (rop, op, invrnd, 1, 2); + /* division by 2 does not change the ternary flag */ + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); + mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); + inex_im = -inex_im; /* negate the ternary flag */ + } + return MPC_INEX(inex_re, inex_im); + } + + /* generic case: neither Re(op) nor Im(op) is NaN, Inf or zero */ + prec = MPC_PREC_RE(rop); + mpfr_init2 (w, prec); + mpc_init2 (ww, prec); + /* let op = x + iy; compute log(op)/log(10) */ + while (ok == 0) + { + loops ++; + prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; + mpfr_set_prec (w, prec); + mpc_set_prec (ww, prec); + + mpc_log (ww, op, MPC_RNDNN); + mpfr_set_ui (w, 10, GMP_RNDN); /* exact since prec >= 4 */ + mpfr_log (w, w, GMP_RNDN); + mpc_div_fr (ww, ww, w, MPC_RNDNN); + + ok = mpfr_can_round (mpc_realref (ww), prec - 2, GMP_RNDZ, GMP_RNDZ, + MPC_PREC_RE(rop) + MPC_RND_RE (rnd) == GMP_RNDN); + ok = ok && mpfr_can_round (mpc_imagref (ww), prec-2, GMP_RNDZ, GMP_RNDZ, + MPC_PREC_IM(rop) + MPC_RND_IM (rnd) == GMP_RNDN); + } + + inex_re = mpfr_set (mpc_realref(rop), mpc_realref (ww), MPC_RND_RE (rnd)); + inex_im = mpfr_set (mpc_imagref(rop), mpc_imagref (ww), MPC_RND_IM (rnd)); + mpfr_clear (w); + mpc_clear (ww); + return MPC_INEX(inex_re, inex_im); +} diff --git a/src/mpc.h b/src/mpc.h index ab87e7c..c98c3fb 100644 --- a/src/mpc.h +++ b/src/mpc.h @@ -188,6 +188,7 @@ __MPC_DECLSPEC int mpc_cmp (mpc_srcptr, mpc_srcptr); __MPC_DECLSPEC int mpc_cmp_si_si (mpc_srcptr, long int, long int); __MPC_DECLSPEC int mpc_exp (mpc_ptr, mpc_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_log (mpc_ptr, mpc_srcptr, mpc_rnd_t); +__MPC_DECLSPEC int mpc_log10 (mpc_ptr, mpc_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_sin (mpc_ptr, mpc_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_cos (mpc_ptr, mpc_srcptr, mpc_rnd_t); __MPC_DECLSPEC int mpc_sin_cos (mpc_ptr, mpc_ptr, mpc_srcptr, mpc_rnd_t, mpc_rnd_t); diff --git a/tests/Makefile.am b/tests/Makefile.am index d4f425e..dede832 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -29,10 +29,10 @@ LOADLIBES=$(DEFS) -I$(top_srcdir)/src -I$(top_builddir) $(CPPFLAGS) $(CFLAGS) -L check_PROGRAMS = tabs tacos tacosh tadd tadd_fr tadd_si tadd_ui targ tasin \ tasinh tatan tatanh tconj tcos tcosh tdiv tdiv_2exp tdiv_fr tdiv_ui texp tfma \ -tfr_div tfr_sub timag tio_str tlog tmul tmul_2exp tmul_fr tmul_i tmul_si \ -tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_ui tpow_z tprec \ -tproj treal treimref tset tsin tsin_cos tsinh tsqr tsqrt tstrtoc tsub tsub_fr \ -tsub_ui tswap ttan ttanh tui_div tui_ui_sub tget_version +tfr_div tfr_sub timag tio_str tlog tlog10 tmul tmul_2exp tmul_fr tmul_i \ +tmul_si tmul_ui tneg tnorm tpow tpow_ld tpow_d tpow_fr tpow_si tpow_ui tpow_z \ +tprec tproj treal treimref tset tsin tsin_cos tsinh tsqr tsqrt tstrtoc tsub \ +tsub_fr tsub_ui tswap ttan ttanh tui_div tui_ui_sub tget_version check_LTLIBRARIES=libmpc-tests.la libmpc_tests_la_SOURCES=mpc-tests.h random.c tgeneric.c read_data.c \ diff --git a/tests/log.dat b/tests/log.dat index 7222fb5..73818f3 100644 --- a/tests/log.dat +++ b/tests/log.dat @@ -95,6 +95,7 @@ # log(nan + i*inf) = +inf + i*nan 0 0 2 +inf 2 nan 2 nan 2 +inf N N +0 0 2 +inf 2 nan 2 nan 2 -inf N N # log(nan + i*nan) = nan + i*nan 0 0 2 nan 2 nan 2 nan 2 nan N N diff --git a/tests/log10.dat b/tests/log10.dat new file mode 100644 index 0000000..acd24dd --- /dev/null +++ b/tests/log10.dat @@ -0,0 +1,172 @@ +# Data test file for mpc_log10. +# +# Copyright (C) 2012 INRIA +# +# This file is part of GNU MPC. +# +# GNU MPC 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 3 of the License, or (at your +#o ption) any later version. +# +# GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . + +# See file sin.dat for the format description. + +# Special values, following ISO C99 standard, Annex G, +# more precisely Section 7.26 "Future library directions" + +# Rule [conj]: log10(conj(z)) = conj(log10(z)) + +# log10(nan + i*inf) = +inf + i*nan +0 0 2 +inf 2 nan 2 nan 2 +inf N N +0 0 2 +inf 2 nan 2 nan 2 -inf N N + +# log10(nan + i*nan) = nan + i*nan +0 0 2 nan 2 nan 2 nan 2 nan N N + +# log10(nan + i*y) = nan + i*nan for finite y +0 0 2 nan 2 nan 2 nan 2 1 N N +0 0 2 nan 2 nan 2 nan 2 +0 N N +0 0 2 nan 2 nan 2 nan 2 -0 N N +0 0 2 nan 2 nan 2 nan 2 -1 N N + +# log10(+/-inf + i*nan) = +inf + i*nan +0 0 2 +inf 2 nan 2 +inf 2 nan N N +0 0 2 +inf 2 nan 2 -inf 2 nan N N + +# log10(x + i*nan) = nan+i*nan for finite x +0 0 2 nan 2 nan 2 1 2 nan N N +0 0 2 nan 2 nan 2 +0 2 nan N N +0 0 2 nan 2 nan 2 -0 2 nan N N +0 0 2 nan 2 nan 2 -1 2 nan N N + +# log10(x + i*inf) = +inf + i*pi/2/log(10) for finite x +0 + 2 +inf 53 +0x15d47c4cb2fba1p-53 2 1 2 +inf N N +0 + 2 +inf 53 +0x15d47c4cb2fba1p-53 2 +0 2 +inf N N +0 + 2 +inf 53 +0x15d47c4cb2fba1p-53 2 -0 2 +inf N N +0 + 2 +inf 53 +0x15d47c4cb2fba1p-53 2 -1 2 +inf N N +# by [conj]: log10(x - i*inf) = +inf - i*pi/2/log(10) for finite x +0 - 2 +inf 53 -0x15d47c4cb2fba1p-53 2 1 2 -inf N N +0 - 2 +inf 53 -0x15d47c4cb2fba1p-53 2 +0 2 -inf N N +0 - 2 +inf 53 -0x15d47c4cb2fba1p-53 2 -0 2 -inf N N +0 - 2 +inf 53 -0x15d47c4cb2fba1p-53 2 -1 2 -inf N N + +# log10(-inf + i*y) = +inf + i*pi/log(10) for finite positive-signed y +0 + 2 +inf 53 0x15d47c4cb2fba1p-52 2 -inf 2 1 N N +0 + 2 +inf 53 0x15d47c4cb2fba1p-52 2 -inf 2 +0 N N +# by [conj]: log10(-inf+i*y) = +inf - i*pi/log(10) for finite negative-signed y +0 - 2 +inf 53 -0x15d47c4cb2fba1p-52 2 -inf 2 -1 N N +0 - 2 +inf 53 -0x15d47c4cb2fba1p-52 2 -inf 2 -0 N N + +# log10(+inf + i*y) = +inf + i*0 for finite positive-signed y +0 0 2 +inf 2 +0 2 +inf 2 1 N N +0 0 2 +inf 2 +0 2 +inf 2 +0 N N +# by [conj]: log10(+inf + i*y) = +inf - i*0 for finite negative-signed y +0 0 2 +inf 2 -0 2 +inf 2 -1 N N +0 0 2 +inf 2 -0 2 +inf 2 -0 N N + +# log10(-inf + i*inf) = +inf + i*(3*pi/4)/log(10) +0 + 2 +inf 53 0x105f5d39863cb9p-52 2 -inf 2 +inf N N +# by [conj]: log10(-inf - i*inf) = +inf - i*(3*pi/4)/log(10) +0 - 2 +inf 53 -0x105f5d39863cb9p-52 2 -inf 2 -inf N N + +# log10(+inf + i*inf) = +inf + i*(pi/4)/log(10) +0 + 2 +inf 53 0x15d47c4cb2fba1p-54 2 +inf 2 +inf N N +# by [conj]: log10(+inf - i*inf) = +inf - i*(pi/4)/log(10) +0 - 2 +inf 53 -0x15d47c4cb2fba1p-54 2 +inf 2 -inf N N + +# log10(-0 + i*0) = -inf + i*pi/log(10) +0 + 2 -inf 53 +0x15d47c4cb2fba1p-52 2 -0 2 +0 N N +0 - 2 -inf 53 -0x15d47c4cb2fba1p-52 2 -0 2 -0 N N +0 + 2 -inf 53 +0x15d47c4cb2fba1p-52 2 -0 2 +0 N U +0 + 2 -inf 53 -0xaea3e26597ddp-47 2 -0 2 -0 N U +0 - 2 -inf 53 +0xaea3e26597ddp-47 2 -0 2 +0 N D +0 - 2 -inf 53 -0x15d47c4cb2fba1p-52 2 -0 2 -0 N D + +# log10(+0 + i*0) = -inf + i*0 +0 0 2 -inf 2 +0 2 +0 2 +0 N N +0 0 2 -inf 2 -0 2 +0 2 -0 N N + +# log10(+1 +- i*0) = +0 +- i*0 +0 0 2 +0 2 +0 2 1 2 +0 N N +0 0 2 +0 2 -0 2 1 2 -0 N N + +# log10(10 +- i*0) = 1 +- i*0 +0 0 2 1 2 +0 4 10 2 +0 N N +0 0 2 1 2 -0 4 10 2 -0 N N + +# log10(100 +- i*0) = 1 +- i*0 +0 0 2 2 2 +0 5 100 2 +0 N N +0 0 2 2 2 -0 5 100 2 -0 N N + +# log10(-1 +- i*0) = +0 +- i*pi/log(10) +0 + 2 +0 53 +0x15d47c4cb2fba1p-52 2 -1 2 +0 N N +0 - 2 +0 53 -0x15d47c4cb2fba1p-52 2 -1 2 -0 N N + +# log10(x + i*y) with either x or y zero and the other non-zero +- 0 53 0x13afeb354b7d97p-52 2 0 5 0x11 2 0 N N +- + 53 0x13afeb354b7d97p-52 53 0x15d47c4cb2fba1p-53 2 0 5 0x11 N N +- + 53 0x1475c655fbc11p-48 53 0x15d47c4cb2fba1p-52 5 -0x13 2 +0 N N +- - 53 0x1475c655fbc11p-48 53 -0x15d47c4cb2fba1p-52 5 -0x13 2 -0 N N +- - 53 0x1475c655fbc11p-48 53 -0x15d47c4cb2fba1p-53 2 0 5 -0x13 N N + +- + 53 0x15c9a3209bf97fp-52 53 0x15d47c4cb2fba1p-52 5 -0x17 2 +0 Z U +- - 53 0x15c9a3209bf97fp-52 53 -0x15d47c4cb2fba1p-52 5 -0x17 2 -0 N D ++ - 53 0x2b93464137f3p-45 53 0xaea3e26597ddp-47 5 -0x17 2 +0 U Z +- - 53 0x15c9a3209bf97fp-52 53 -0x15d47c4cb2fba1p-52 5 -0x17 2 -0 D N +- - 53 0x1d1cda1a0c996dp-52 53 0xaea3e26597ddp-47 7 -0x42 2 +0 Z D ++ + 53 0xe8e6d0d064cb7p-51 53 -0xaea3e26597ddp-47 7 -0x42 2 -0 N Z ++ + 53 0xe8e6d0d064cb7p-51 53 0x15d47c4cb2fba1p-52 7 -0x42 2 +0 U N +- + 53 0x1d1cda1a0c996dp-52 53 -0xaea3e26597ddp-47 7 -0x42 2 -0 D U + +# huge values ++ 0 53 0x134413509f79ffp-44 2 0 2 0x1p1024 2 +0 U U ++ 0 53 0x134413509f79ffp-43 2 0 2 0x1p2048 2 +0 U U ++ 0 53 0x134413509f79ffp-42 2 0 2 0x1p4096 2 +0 U U ++ 0 53 0x134413509f79ffp-41 2 0 2 0x1p8192 2 +0 U U ++ 0 53 0x134413509f79ffp-40 2 0 2 0x1p16384 2 +0 U U ++ 0 53 0x134413509f79ffp-39 2 0 2 0x1p32768 2 +0 U U ++ 0 53 0x134413509f79ffp-38 2 0 2 0x1p65536 2 +0 U U ++ 0 53 0x134413509f79ffp-37 2 0 2 0x1p131072 2 +0 U U ++ 0 53 0x134413509f79ffp-36 2 0 2 0x1p262144 2 +0 U U ++ 0 53 0x134413509f79ffp-35 2 0 2 0x1p524288 2 +0 U U ++ 0 53 0x134413509f79ffp-34 2 0 2 0x1p1048576 2 +0 U U ++ 0 53 0x134413509f79ffp-33 2 0 2 0x1p2097152 2 +0 U U ++ 0 53 0x134413509f79ffp-32 2 0 2 0x1p4194304 2 +0 U U ++ 0 53 0x134413509f79ffp-31 2 0 2 0x1p8388608 2 +0 U U ++ 0 53 0x134413509f79ffp-30 2 0 2 0x1p16777216 2 +0 U U ++ 0 53 0x134413509f79ffp-29 2 0 2 0x1p33554432 2 +0 U U ++ 0 53 0x134413509f79ffp-28 2 0 2 0x1p67108864 2 +0 U U ++ 0 53 0x134413509f79ffp-27 2 0 2 0x1p134217728 2 +0 U U ++ 0 53 0x134413509f79ffp-26 2 0 2 0x1p268435456 2 +0 U U ++ 0 53 0x134413509f79ffp-25 2 0 2 0x1p536870912 2 +0 U U ++ + 53 0x13467bd3098defp-44 53 0x15d47c4cb2fba1p-54 2 0x1p1024 2 0x1p1024 U U ++ + 53 0x13454791d483f7p-43 53 0x15d47c4cb2fba1p-54 2 0x1p2048 2 0x1p2048 U U ++ + 53 0x1344ad7139fefbp-42 53 0x15d47c4cb2fba1p-54 2 0x1p4096 2 0x1p4096 U U ++ + 53 0x13446060ecbc7dp-41 53 0x15d47c4cb2fba1p-54 2 0x1p8192 2 0x1p8192 U U ++ + 53 0x9a21cec630d9fp-39 53 0x15d47c4cb2fba1p-54 2 0x1p16384 2 0x1p16384 U U ++ + 53 0x13442694b2ca9fp-39 53 0x15d47c4cb2fba1p-54 2 0x1p32768 2 0x1p32768 U U ++ + 53 0x13441cf2a9224fp-38 53 0x15d47c4cb2fba1p-54 2 0x1p65536 2 0x1p65536 U U ++ + 53 0x13441821a44e27p-37 2 0x3p-3 2 0x1p131072 2 0x1p131072 U U ++ + 53 0x134415b921e413p-36 2 0x3p-3 2 0x1p262144 2 0x1p262144 U U ++ + 53 0x13441484e0af09p-35 2 0x3p-3 2 0x1p524288 2 0x1p524288 U U ++ + 53 0x4d104fab00521p-32 2 0x3p-3 2 0x1p1048576 2 0x1p1048576 U U ++ + 53 0x9a209ced7e3a1p-32 2 0x3p-3 2 0x1p2097152 2 0x1p2097152 U U ++ + 53 0x1344137727a0a1p-32 2 0x3p-3 2 0x1p4194304 2 0x1p4194304 U U ++ + 53 0x13441363e38d5p-27 2 0x3p-3 2 0x1p8388608 2 0x1p8388608 U U ++ + 53 0x268826b483075p-27 2 0x3p-3 2 0x1p16777216 2 0x1p16777216 U U ++ + 53 0x4d104d55c1fb5p-27 2 0x3p-3 2 0x1p33554432 2 0x1p33554432 U U ++ + 53 0x9a209a983fe35p-27 2 0x3p-3 2 0x1p67108864 2 0x1p67108864 U U ++ + 53 0x4d104d474eecdp-25 2 0x3p-3 2 0x1p134217728 2 0x1p134217728 U U ++ + 53 0x9a209a89ccd4dp-25 2 0x3p-3 2 0x1p268435456 2 0x1p268435456 U U +# Due to intermediate overflow, the following result has wrong real part +# instead of the correct result. Since this may happen in other parts of the +# library as well, we do not consider it a bug for the moment. +# + + 53 0x13441350ec8a4dp-25 2 0x3p-3 2 0x1p536870912 2 0x1p536870912 U U diff --git a/tests/tlog10.c b/tests/tlog10.c new file mode 100644 index 0000000..5591b06 --- /dev/null +++ b/tests/tlog10.c @@ -0,0 +1,37 @@ +/* tlog10 -- test file for mpc_log10. + +Copyright (C) 2012 INRIA + +This file is part of GNU MPC. + +GNU MPC 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 3 of the License, or (at your +option) any later version. + +GNU MPC 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 this program. If not, see http://www.gnu.org/licenses/ . +*/ + +#include +#include "mpc-tests.h" + +int +main (void) +{ + DECL_FUNC (CC, f, mpc_log10); + + test_start (); + + data_check (f, "log10.dat"); + tgeneric (f, 2, 512, 7, 128); + + test_end (); + + return 0; +} -- cgit v1.2.1