From 6c81cba454ca8bf0ec2af7a84de8e576acc04c23 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Tue, 21 Jul 2015 18:05:30 +0200 Subject: [PATCH] port the div_2expm1_ui code completely to the arb type --- fmprb/div_2expm1.c => arb/div_2expm1_ui.c | 61 +++++++------- arb/fmprb_wrappers.c | 37 --------- arb/test/t-div_2expm1_ui.c | 21 ++++- doc/source/fmprb.rst | 4 - fmprb.h | 3 - fmprb/test/t-div_2expm1_ui.c | 96 ----------------------- 6 files changed, 50 insertions(+), 172 deletions(-) rename fmprb/div_2expm1.c => arb/div_2expm1_ui.c (61%) delete mode 100644 arb/fmprb_wrappers.c delete mode 100644 fmprb/test/t-div_2expm1_ui.c diff --git a/fmprb/div_2expm1.c b/arb/div_2expm1_ui.c similarity index 61% rename from fmprb/div_2expm1.c rename to arb/div_2expm1_ui.c index cd49b16f..d9b59c73 100644 --- a/fmprb/div_2expm1.c +++ b/arb/div_2expm1_ui.c @@ -19,61 +19,64 @@ =============================================================================*/ /****************************************************************************** - Copyright (C) 2012 Fredrik Johansson + Copyright (C) 2015 Fredrik Johansson ******************************************************************************/ -#include "fmprb.h" +#include "arb.h" - - -/* TODO: fix/optimize for n > prec */ void -fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec) +arb_div_2expm1_ui(arb_t y, const arb_t x, ulong n, long prec) { if (n < FLINT_BITS) { - fmprb_div_ui(y, x, (1UL << n) - 1, prec); + arb_div_ui(y, x, (1UL << n) - 1, prec); } - else if (n < 1024 + prec / 32 || n > LONG_MAX / 2) + else if (n < 1024 + prec / 32 || n > LONG_MAX / 4) { - fmprb_t t; - fmprb_init(t); - fmprb_one(t); - fmpz_set_ui(fmpr_expref(fmprb_midref(t)), n); - fmprb_sub_ui(t, t, 1, prec); - fmprb_div(y, x, t, prec); - fmprb_clear(t); + arb_t t; + fmpz_t e; + + arb_init(t); + fmpz_init_set_ui(e, n); + + arb_one(t); + arb_mul_2exp_fmpz(t, t, e); + arb_sub_ui(t, t, 1, prec); + arb_div(y, x, t, prec); + + arb_clear(t); + fmpz_clear(e); } else { - fmprb_t s, t; + arb_t s, t; long i, b; - fmprb_init(s); - fmprb_init(t); + arb_init(s); + arb_init(t); /* x / (2^n - 1) = sum_{k>=1} x * 2^(-k*n)*/ - - fmprb_mul_2exp_si(s, x, -n); - fmprb_set(t, s); + arb_mul_2exp_si(s, x, -n); + arb_set(t, s); b = 1; for (i = 2; i <= prec / n + 1; i++) { - fmprb_mul_2exp_si(t, t, -n); - fmprb_add(s, s, t, prec); + arb_mul_2exp_si(t, t, -n); + arb_add(s, s, t, prec); b = i; } /* error bound: sum_{k>b} x * 2^(-k*n) <= x * 2^(-b*n - (n-1)) */ - fmprb_mul_2exp_si(t, x, -b*n - (n-1)); - fmprb_abs(t, t); - fmprb_add_error(s, t); + arb_mul_2exp_si(t, x, -b*n - (n-1)); + arb_abs(t, t); + arb_add_error(s, t); - fmprb_set(y, s); + arb_set(y, s); - fmprb_clear(s); - fmprb_clear(t); + arb_clear(s); + arb_clear(t); } } + diff --git a/arb/fmprb_wrappers.c b/arb/fmprb_wrappers.c deleted file mode 100644 index 4645347e..00000000 --- a/arb/fmprb_wrappers.c +++ /dev/null @@ -1,37 +0,0 @@ -/*============================================================================= - - 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 "arb.h" - -void arb_div_2expm1_ui(arb_t z, const arb_t x, ulong n, long prec) -{ - fmprb_t t; - fmprb_init(t); - arb_get_fmprb(t, x); - fmprb_div_2expm1_ui(t, t, n, prec); - arb_set_fmprb(z, t); - fmprb_clear(t); -} - diff --git a/arb/test/t-div_2expm1_ui.c b/arb/test/t-div_2expm1_ui.c index 0c51d061..ec545cfa 100644 --- a/arb/test/t-div_2expm1_ui.c +++ b/arb/test/t-div_2expm1_ui.c @@ -49,12 +49,26 @@ int main() prec = 2 + n_randint(state, 10000); arb_randtest(a, state, 1 + n_randint(state, 10000), 10); - n = 1 + (n_randtest(state) % (10 * prec)); + + if (n_randint(state, 2)) + n = 1 + (n_randtest(state) % (10 * prec)); + else + n = n_randtest(state); arb_div_2expm1_ui(b, a, n, prec); arb_one(c); - arb_mul_2exp_si(c, c, n); + if (n >= (1UL << (FLINT_BITS-1))) + { + arb_mul_2exp_si(c, c, (1UL << (FLINT_BITS-2))); + arb_mul_2exp_si(c, c, (1UL << (FLINT_BITS-2))); + arb_mul_2exp_si(c, c, n - (1UL << (FLINT_BITS-1))); + } + else + { + arb_mul_2exp_si(c, c, n); + } + arb_sub_ui(c, c, 1, prec); arb_div(c, a, c, prec); @@ -64,13 +78,14 @@ int main() if (!arb_overlaps(b, c)) { printf("FAIL: containment\n\n"); + printf("n = %lu\n", n); printf("a = "); arb_print(a); printf("\n\n"); printf("b = "); arb_print(b); printf("\n\n"); printf("c = "); arb_print(c); printf("\n\n"); abort(); } - if ((acc2 < FLINT_MIN(prec, acc1) - 10) && + if (n > 0 && (acc2 < FLINT_MIN(prec, acc1) - 10) && !(acc1 == -ARF_PREC_EXACT && acc2 == -ARF_PREC_EXACT)) { printf("FAIL: poor accuracy\n\n"); diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index 9e8cb49d..19e4bf0d 100644 --- a/doc/source/fmprb.rst +++ b/doc/source/fmprb.rst @@ -437,10 +437,6 @@ Arithmetic where the triangle inequality has been applied to the numerator and the reverse triangle inequality has been applied to the denominator. -.. function:: void fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec) - - Sets `y = x / (2^n - 1)`, rounded to *prec* bits. - .. function:: void fmprb_addmul(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec) .. function:: void fmprb_addmul_ui(fmprb_t z, const fmprb_t x, ulong y, long prec) diff --git a/fmprb.h b/fmprb.h index 1436dbb5..7797e713 100644 --- a/fmprb.h +++ b/fmprb.h @@ -292,9 +292,6 @@ fmprb_inv(fmprb_t y, const fmprb_t x, long prec) fmprb_ui_div(y, 1, x, prec); } -void fmprb_div_2expm1_ui(fmprb_t y, const fmprb_t x, ulong n, long prec); - - void fmprb_mul_fmpr_naive(fmprb_t z, const fmprb_t x, const fmpr_t y, long prec); void fmprb_mul_main_naive(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec); void fmprb_mul_naive(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec); diff --git a/fmprb/test/t-div_2expm1_ui.c b/fmprb/test/t-div_2expm1_ui.c deleted file mode 100644 index dbab3ed3..00000000 --- a/fmprb/test/t-div_2expm1_ui.c +++ /dev/null @@ -1,96 +0,0 @@ -/*============================================================================= - - 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) 2013 Fredrik Johansson - -******************************************************************************/ - -#include "fmprb.h" - -int main() -{ - long iter; - flint_rand_t state; - - printf("div_2expm1_ui...."); - fflush(stdout); - - flint_randinit(state); - - for (iter = 0; iter < 100000; iter++) - { - fmprb_t a, b, c; - ulong n; - long prec, acc1, acc2; - fmpz_t t; - - fmprb_init(a); - fmprb_init(b); - fmprb_init(c); - fmpz_init(t); - - prec = 2 + n_randint(state, 10000); - fmprb_randtest(a, state, 1 + n_randint(state, 10000), 10); - n = 1 + (n_randtest(state) % (10 * prec)); - - fmprb_div_2expm1_ui(b, a, n, prec); - - fmprb_one(c); - fmpz_set_ui(fmpr_expref(fmprb_midref(c)), n); - fmprb_sub_ui(c, c, 1, prec); - fmprb_div(c, a, c, prec); - - acc1 = fmprb_rel_accuracy_bits(a); - acc2 = fmprb_rel_accuracy_bits(b); - - if (!fmprb_overlaps(b, c)) - { - printf("FAIL: containment\n\n"); - printf("a = "); fmprb_print(a); printf("\n\n"); - printf("b = "); fmprb_print(b); printf("\n\n"); - printf("c = "); fmprb_print(c); printf("\n\n"); - abort(); - } - - if ((acc2 < FLINT_MIN(prec, acc1) - 10) && - !(acc1 == -FMPR_PREC_EXACT && acc2 == -FMPR_PREC_EXACT)) - { - printf("FAIL: poor accuracy\n\n"); - printf("prec=%ld, acc1=%ld, acc2=%ld\n\n", prec, acc1, acc2); - printf("n = %lu\n\n", n); - printf("a = "); fmprb_print(a); printf("\n\n"); - printf("b = "); fmprb_print(b); printf("\n\n"); - printf("c = "); fmprb_print(c); printf("\n\n"); - abort(); - } - - fmprb_clear(a); - fmprb_clear(b); - fmprb_clear(c); - fmpz_clear(t); - } - - flint_randclear(state); - flint_cleanup(); - printf("PASS\n"); - return EXIT_SUCCESS; -} -