From 7d04c8950426713bf42abea1d3f35a1ec7a63fcd Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 6 Nov 2014 14:29:56 +0100 Subject: [PATCH] implement mag_sub_lower properly and test+document it --- acb_modular/eta_sum.c | 2 - acb_modular/theta_sum.c | 2 - acb_poly/polylog_series.c | 3 - doc/source/mag.rst | 4 ++ hypgeom/bound.c | 36 ------------ mag.h | 2 + mag/sub_lower.c | 100 +++++++++++++++++++++++++++++++ mag/test/t-sub_lower.c | 120 ++++++++++++++++++++++++++++++++++++++ 8 files changed, 226 insertions(+), 43 deletions(-) create mode 100644 mag/sub_lower.c create mode 100644 mag/test/t-sub_lower.c diff --git a/acb_modular/eta_sum.c b/acb_modular/eta_sum.c index e3f487d4..99a7f925 100644 --- a/acb_modular/eta_sum.c +++ b/acb_modular/eta_sum.c @@ -45,8 +45,6 @@ acb_mul_approx(acb_t z, acb_t tmp1, acb_t tmp2, const acb_t x, const acb_t y, lo } } -void mag_sub_lower(mag_t z, const mag_t x, const mag_t y); - double mag_get_log2_d_approx(const mag_t x); #define PENTAGONAL(N) ((((N)+2)/2) * ((3*(N)+5)/2)/2) diff --git a/acb_modular/theta_sum.c b/acb_modular/theta_sum.c index 20caa024..217d1520 100644 --- a/acb_modular/theta_sum.c +++ b/acb_modular/theta_sum.c @@ -45,8 +45,6 @@ acb_mul_approx(acb_t z, acb_t tmp1, acb_t tmp2, const acb_t x, const acb_t y, lo } } -void mag_sub_lower(mag_t z, const mag_t x, const mag_t y); - double mag_get_log2_d_approx(const mag_t x) { diff --git a/acb_poly/polylog_series.c b/acb_poly/polylog_series.c index b82ccabc..a3431384 100644 --- a/acb_poly/polylog_series.c +++ b/acb_poly/polylog_series.c @@ -43,9 +43,6 @@ mag_log_ui(mag_t t, ulong n) } } -/* TODO: needs reimplementing */ -void mag_sub_lower(mag_t z, const mag_t x, const mag_t y); - long arb_get_si_lower(const arb_t x) { diff --git a/doc/source/mag.rst b/doc/source/mag.rst index 41c8faf6..cf8e6931 100644 --- a/doc/source/mag.rst +++ b/doc/source/mag.rst @@ -220,6 +220,10 @@ Arithmetic Sets *z* to a lower bound for `x + y`. +.. function:: void mag_sub_lower(mag_t z, const mag_t x, const mag_t y) + + Sets *z* to a lower bound for `\max(x-y, 0)`. + Fast, unsafe arithmetic ------------------------------------------------------------------------------- diff --git a/hypgeom/bound.c b/hypgeom/bound.c index dd4f60c0..5173917e 100644 --- a/hypgeom/bound.c +++ b/hypgeom/bound.c @@ -207,42 +207,6 @@ hypgeom_term_bound(mag_t Tn, const mag_t TK, long K, long A, long B, int r, cons mag_clear(num); } - -/* z = max(x-y, 0), rounding down [lower bound] */ -void -mag_sub_lower(mag_t z, const mag_t x, const mag_t y) -{ - if (mag_is_special(x) || mag_is_special(y)) - { - if (mag_is_zero(y)) - mag_set(z, x); - else if (mag_is_zero(x) || mag_is_inf(y)) - mag_zero(z); - else - mag_inf(z); - } - else - { - fmpr_t t, u; - - fmpr_init(t); - fmpr_init(u); - - mag_get_fmpr(t, x); - mag_get_fmpr(u, y); - - fmpr_sub(t, t, u, MAG_BITS, FMPR_RND_DOWN); - - if (fmpr_sgn(t) <= 0) - mag_zero(z); - else /* XXX: exact? */ - mag_set_fmpz_2exp_fmpz_lower(z, fmpr_manref(t), fmpr_expref(t)); - - fmpr_clear(t); - fmpr_clear(u); - } -} - long hypgeom_bound(mag_t error, int r, long A, long B, long K, const mag_t TK, const mag_t z, long tol_2exp) diff --git a/mag.h b/mag.h index b0928709..1124d380 100644 --- a/mag.h +++ b/mag.h @@ -331,6 +331,8 @@ mag_mul_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t y) } } +void mag_sub_lower(mag_t z, const mag_t x, const mag_t y); + /* Fast versions (no infs/nans, small exponents). Note that this applies to outputs too! */ diff --git a/mag/sub_lower.c b/mag/sub_lower.c new file mode 100644 index 00000000..607421df --- /dev/null +++ b/mag/sub_lower.c @@ -0,0 +1,100 @@ +/*============================================================================= + + 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 "arf.h" + +void +mag_sub_lower(mag_t z, const mag_t x, const mag_t y) +{ + if (mag_is_special(x) || mag_is_special(y)) + { + if (mag_is_zero(y)) + mag_set(z, x); + else if (mag_is_zero(x) || mag_is_inf(y)) + mag_zero(z); + else + mag_inf(z); + } + else + { + long shift, fix; + + shift = _fmpz_sub_small(MAG_EXPREF(x), MAG_EXPREF(y)); + + /* y > x */ + if (shift < 0) + { + mag_zero(z); + } + else if (shift == 0) + { + if (MAG_MAN(y) >= MAG_MAN(x)) + { + mag_zero(z); + } + else + { + MAG_MAN(z) = MAG_MAN(x) - MAG_MAN(y); + fix = MAG_BITS - FLINT_BIT_COUNT(MAG_MAN(z)); + MAG_MAN(z) <<= fix; + _fmpz_add_fast(MAG_EXPREF(z), MAG_EXPREF(x), -fix); + } + } + else + { + if (shift <= MAG_BITS) + { + mp_limb_t c = MAG_MAN(x) - (MAG_MAN(y) >> shift) - 1; + + /* too much cancellation -- compute precisely */ + if (c < (1UL << (MAG_BITS - 4))) + { + arf_t t, u; + arf_init(t); + arf_init(u); + arf_set_mag(t, x); + arf_set_mag(u, y); + arf_sub(t, t, u, MAG_BITS, ARF_RND_DOWN); + arf_get_mag_lower(z, t); + arf_clear(t); + arf_clear(u); + return; + } + + MAG_MAN(z) = c; + } + else + { + MAG_MAN(z) = MAG_MAN(x) - 1; + } + + fix = MAG_BITS - FLINT_BIT_COUNT(MAG_MAN(z)); + MAG_MAN(z) <<= fix; + _fmpz_add_fast(MAG_EXPREF(z), MAG_EXPREF(x), -fix); + } + } +} + diff --git a/mag/test/t-sub_lower.c b/mag/test/t-sub_lower.c new file mode 100644 index 00000000..147abc4c --- /dev/null +++ b/mag/test/t-sub_lower.c @@ -0,0 +1,120 @@ +/*============================================================================= + + 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("sub_lower...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmpr_t x, y, z, z2, w; + mag_t xb, yb, zb; + + fmpr_init(x); + fmpr_init(y); + fmpr_init(z); + fmpr_init(z2); + fmpr_init(w); + + mag_init(xb); + mag_init(yb); + mag_init(zb); + + mag_randtest_special(xb, state, 100); + mag_randtest_special(yb, state, 100); + + mag_get_fmpr(x, xb); + mag_get_fmpr(y, yb); + + fmpr_sub(z2, x, y, 100, FMPR_RND_DOWN); + if (fmpr_sgn(z2) < 0) + fmpr_zero(z2); + + fmpr_mul_ui(z, z2, 1023, MAG_BITS, FMPR_RND_DOWN); + fmpr_mul_2exp_si(z, z, -10); + + mag_sub_lower(zb, xb, yb); + mag_get_fmpr(w, zb); + + MAG_CHECK_BITS(xb) + MAG_CHECK_BITS(yb) + MAG_CHECK_BITS(zb) + + if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, 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("w = "); fmpr_print(w); printf("\n\n"); + abort(); + } + + if (n_randint(state, 2)) + { + mag_sub_lower(xb, xb, yb); + + if (!mag_equal(xb, zb)) + { + printf("FAIL (aliasing 1)\n\n"); + abort(); + } + } + else + { + mag_sub_lower(yb, xb, yb); + + if (!mag_equal(yb, zb)) + { + printf("FAIL (aliasing 2)\n\n"); + abort(); + } + } + + fmpr_clear(x); + fmpr_clear(y); + fmpr_clear(z); + fmpr_clear(z2); + fmpr_clear(w); + + mag_clear(xb); + mag_clear(yb); + mag_clear(zb); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +