From 270a65569bcdc553ee8982238f6912d2c16ec11d Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 18 Jun 2014 16:50:15 +0200 Subject: [PATCH] add mag_log1p --- mag.h | 3 ++ mag/log1p.c | 96 +++++++++++++++++++++++++++++++++++++++++++++ mag/test/t-log1p.c | 98 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 197 insertions(+) create mode 100644 mag/log1p.c create mode 100644 mag/test/t-log1p.c diff --git a/mag.h b/mag.h index 272cfac1..6abe3c66 100644 --- a/mag.h +++ b/mag.h @@ -648,6 +648,9 @@ static __inline__ void mag_set_d(mag_t z, double x) double mag_d_log_upper_bound(double x); double mag_d_log_lower_bound(double x); +/* TODO: document */ +void mag_log1p(mag_t z, const mag_t x); + /* TODO: document/test */ void mag_pow_ui(mag_t z, const mag_t x, ulong e); void mag_pow_ui_lower(mag_t z, const mag_t x, ulong e); diff --git a/mag/log1p.c b/mag/log1p.c new file mode 100644 index 00000000..04c9e90e --- /dev/null +++ b/mag/log1p.c @@ -0,0 +1,96 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "mag.h" +#include "double_extras.h" + +void +mag_log1p(mag_t z, const mag_t x) +{ + if (mag_is_special(x)) + { + if (mag_is_zero(x)) + mag_zero(z); + else + mag_inf(z); + } + else + { + fmpz exp = MAG_EXP(x); + + if (!COEFF_IS_MPZ(exp)) + { + /* Quick bound by x */ + if (exp < -10) + { + mag_set(z, x); + return; + } + else if (exp < 1000) + { + double t; + t = ldexp(MAG_MAN(x), exp - MAG_BITS); + t = (1.0 + t) * (1 + 1e-14); + t = mag_d_log_upper_bound(t); + mag_set_d(z, t); + return; + } + } + else if (fmpz_sgn(MAG_EXPREF(x)) < 0) + { + /* Quick bound by x */ + mag_set(z, x); + return; + } + + /* Now we must have x >= 2^1000 */ + /* Use log(2^(exp-1) * (2*v)) = exp*log(2) + log(2*v) */ + { + double t; + fmpz_t b; + mag_t u; + + mag_init(u); + fmpz_init(b); + + /* incrementing the mantissa gives an upper bound for x+1 */ + t = ldexp(MAG_MAN(x) + 1, 1 - MAG_BITS); + t = mag_d_log_upper_bound(t); + mag_set_d(u, t); + + /* log(2) < 744261118/2^30 */ + _fmpz_add_fast(b, MAG_EXPREF(x), -1); + fmpz_mul_ui(b, b, 744261118); + mag_set_fmpz(z, b); + _fmpz_add_fast(MAG_EXPREF(z), MAG_EXPREF(z), -30); + + mag_add(z, z, u); + + mag_clear(u); + fmpz_clear(b); + } + } +} + diff --git a/mag/test/t-log1p.c b/mag/test/t-log1p.c new file mode 100644 index 00000000..d31e3563 --- /dev/null +++ b/mag/test/t-log1p.c @@ -0,0 +1,98 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "mag.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("log1p...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmpr_t x, y, z, z2; + mag_t xb, yb; + + fmpr_init(x); + fmpr_init(y); + fmpr_init(z); + fmpr_init(z2); + + mag_init(xb); + mag_init(yb); + + mag_randtest_special(xb, state, 25); + mag_randtest_special(yb, state, 25); + + mag_log1p(yb, xb); + + mag_get_fmpr(x, xb); + mag_get_fmpr(y, yb); + + fmpr_log1p(z, x, MAG_BITS, FMPR_RND_UP); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + MAG_CHECK_BITS(xb) + MAG_CHECK_BITS(yb) + + if (!(fmpr_cmpabs(z, y) <= 0 && fmpr_cmpabs(y, z2) <= 0)) + { + printf("FAIL\n\n"); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("y = "); fmpr_print(y); printf("\n\n"); + printf("z = "); fmpr_print(z); printf("\n\n"); + printf("z2 = "); fmpr_print(z2); printf("\n\n"); + abort(); + } + + mag_log1p(xb, xb); + + if (!mag_equal(xb, yb)) + { + printf("FAIL (aliasing)\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(y); + fmpr_clear(z); + fmpr_clear(z2); + + mag_clear(xb); + mag_clear(yb); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +