From 3bc3c9e895afd860d2690c990882edf7237c1563 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 25 Feb 2016 02:28:19 +0100 Subject: [PATCH] remove some obsolete fmprb methods --- doc/source/fmprb.rst | 31 ----- fmprb.h | 285 ------------------------------------------- fmprb/agm.c | 81 ------------ fmprb/root.c | 102 ---------------- fmprb/test/t-agm.c | 102 ---------------- fmprb/test/t-trim.c | 91 -------------- fmprb/trim.c | 68 ----------- 7 files changed, 760 deletions(-) delete mode 100644 fmprb/agm.c delete mode 100644 fmprb/root.c delete mode 100644 fmprb/test/t-agm.c delete mode 100644 fmprb/test/t-trim.c delete mode 100644 fmprb/trim.c diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index fa2b1ac3..4ebb796c 100644 --- a/doc/source/fmprb.rst +++ b/doc/source/fmprb.rst @@ -225,18 +225,6 @@ Radius and interval operations Sets *x* to a ball containing the interval `[a, b]`. We require that `a \le b`. -.. function:: slong fmprb_rel_error_bits(const fmprb_t x) - - Returns the effective relative error of *x* measured in bits, defined as - the difference between the position of the top bit in the radius - and the top bit in the midpoint, plus one. - The result is clamped between plus/minus *FMPR_PREC_EXACT*. - -.. function:: slong fmprb_rel_accuracy_bits(const fmprb_t x) - - Returns the effective relative accuracy of *x* measured in bits, - equal to the negative of the return value from *fmprb_rel_error_bits*. - .. function:: slong fmprb_bits(const fmprb_t x) Returns the number of bits needed to represent the absolute value @@ -244,14 +232,6 @@ Radius and interval operations sufficient to represent *x* exactly. Returns 0 if the midpoint of *x* is a special value. -.. function:: void fmprb_trim(fmprb_t y, const fmprb_t x) - - Sets *y* to a trimmed copy of *x*: rounds *x* to a number of bits - equal to the accuracy of *x* (as indicated by its radius), - plus a few guard bits. The resulting ball is guaranteed to - contain *x*, but is more economical if *x* has - less than full accuracy. - .. function:: int fmprb_get_unique_fmpz(fmpz_t z, const fmprb_t x) If *x* contains a unique integer, sets *z* to that value and returns @@ -491,14 +471,3 @@ Powers and roots Sets *z* to the reciprocal square root of *x*, rounded to *prec* bits. At high precision, this is faster than computing a square root. -.. function:: void fmprb_root(fmprb_t z, const fmprb_t x, ulong k, slong prec) - - Sets *z* to the *k*-th root of *x*, rounded to *prec* bits. - As currently implemented, this function is only fast for small - fixed *k*. For large *k* it is better to use :func:`fmprb_pow_fmpq` - or :func:`fmprb_pow`. - -.. function:: void fmprb_agm(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec) - - Sets *z* to the arithmetic-geometric mean of *x* and *y*. - diff --git a/fmprb.h b/fmprb.h index 9716a15f..d2190bf2 100644 --- a/fmprb.h +++ b/fmprb.h @@ -75,27 +75,6 @@ fmprb_clear(fmprb_t x) fmpr_clear(fmprb_radref(x)); } -static __inline__ fmprb_ptr -_fmprb_vec_init(slong n) -{ - slong i; - fmprb_ptr v = (fmprb_ptr) flint_malloc(sizeof(fmprb_struct) * n); - - for (i = 0; i < n; i++) - fmprb_init(v + i); - - return v; -} - -static __inline__ void -_fmprb_vec_clear(fmprb_ptr v, slong n) -{ - slong i; - for (i = 0; i < n; i++) - fmprb_clear(v + i); - flint_free(v); -} - static __inline__ int fmprb_is_exact(const fmprb_t x) { @@ -165,8 +144,6 @@ fmprb_set(fmprb_t x, const fmprb_t y) void fmprb_set_round(fmprb_t z, const fmprb_t x, slong prec); -void fmprb_trim(fmprb_t y, const fmprb_t x); - static __inline__ void fmprb_swap(fmprb_t x, fmprb_t y) { @@ -324,10 +301,6 @@ void fmprb_submul_ui(fmprb_t z, const fmprb_t x, ulong y, slong prec); void fmprb_submul_si(fmprb_t z, const fmprb_t x, slong y, slong prec); void fmprb_submul_fmpz(fmprb_t z, const fmprb_t x, const fmpz_t y, slong prec); -void fmprb_root(fmprb_t z, const fmprb_t x, ulong k, slong prec); - -void fmprb_agm(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec); - static __inline__ void fmprb_print(const fmprb_t x) { @@ -491,38 +464,6 @@ void fmprb_set_interval_fmpr(fmprb_t x, const fmpr_t a, const fmpr_t b, slong pr void fmprb_union(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec); -static __inline__ slong -fmprb_rel_error_bits(const fmprb_t x) -{ - fmpz_t midmag, radmag; - slong result; - - if (fmpr_is_zero(fmprb_radref(x))) - return -FMPR_PREC_EXACT; - if (fmpr_is_special(fmprb_midref(x)) || fmpr_is_special(fmprb_radref(x))) - return FMPR_PREC_EXACT; - - fmpz_init(midmag); - fmpz_init(radmag); - - fmpr_abs_bound_lt_2exp_fmpz(midmag, fmprb_midref(x)); - fmpr_abs_bound_lt_2exp_fmpz(radmag, fmprb_radref(x)); - fmpz_add_ui(radmag, radmag, 1); - - result = _fmpz_sub_small(radmag, midmag); - - fmpz_clear(midmag); - fmpz_clear(radmag); - - return result; -} - -static __inline__ slong -fmprb_rel_accuracy_bits(const fmprb_t x) -{ - return -fmprb_rel_error_bits(x); -} - static __inline__ slong fmprb_bits(const fmprb_t x) { @@ -542,232 +483,6 @@ void fmprb_randtest_special(fmprb_t x, flint_rand_t state, slong prec, slong mag void fmprb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const fmprb_t x, slong bits); -#define DEF_CACHED_CONSTANT(name, comp_func) \ - TLS_PREFIX slong name ## _cached_prec = 0; \ - TLS_PREFIX fmprb_t name ## _cached_value; \ - void name ## _cleanup(void) \ - { \ - fmprb_clear(name ## _cached_value); \ - name ## _cached_prec = 0; \ - } \ - void name(fmprb_t x, slong prec) \ - { \ - if (name ## _cached_prec < prec) \ - { \ - if (name ## _cached_prec == 0) \ - { \ - fmprb_init(name ## _cached_value); \ - flint_register_cleanup_function(name ## _cleanup); \ - } \ - comp_func(name ## _cached_value, prec); \ - name ## _cached_prec = prec; \ - } \ - fmprb_set_round(x, name ## _cached_value, prec); \ - } - -/* vector functions */ - -static __inline__ void -_fmprb_vec_zero(fmprb_ptr A, slong n) -{ - slong i; - for (i = 0; i < n; i++) - fmprb_zero(A + i); -} - -static __inline__ int -_fmprb_vec_is_zero(fmprb_srcptr vec, slong len) -{ - slong i; - for (i = 0; i < len; i++) - if (!fmprb_is_zero(vec + i)) - return 0; - return 1; -} - -static __inline__ void -_fmprb_vec_set(fmprb_ptr res, fmprb_srcptr vec, slong len) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_set(res + i, vec + i); -} - -static __inline__ void -_fmprb_vec_set_round(fmprb_ptr res, fmprb_srcptr vec, slong len, slong prec) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_set_round(res + i, vec + i, prec); -} - -static __inline__ void -_fmprb_vec_swap(fmprb_ptr res, fmprb_ptr vec, slong len) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_swap(res + i, vec + i); -} - -static __inline__ void -_fmprb_vec_neg(fmprb_ptr B, fmprb_srcptr A, slong n) -{ - slong i; - for (i = 0; i < n; i++) - fmprb_neg(B + i, A + i); -} - -static __inline__ void -_fmprb_vec_sub(fmprb_ptr C, fmprb_srcptr A, - fmprb_srcptr B, slong n, slong prec) -{ - slong i; - for (i = 0; i < n; i++) - fmprb_sub(C + i, A + i, B + i, prec); -} - -static __inline__ void -_fmprb_vec_add(fmprb_ptr C, fmprb_srcptr A, - fmprb_srcptr B, slong n, slong prec) -{ - slong i; - for (i = 0; i < n; i++) - fmprb_add(C + i, A + i, B + i, prec); -} - -static __inline__ void -_fmprb_vec_scalar_mul(fmprb_ptr res, fmprb_srcptr vec, - slong len, const fmprb_t c, slong prec) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_mul(res + i, vec + i, c, prec); -} - -static __inline__ void -_fmprb_vec_scalar_div(fmprb_ptr res, fmprb_srcptr vec, - slong len, const fmprb_t c, slong prec) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_div(res + i, vec + i, c, prec); -} - -static __inline__ void -_fmprb_vec_scalar_mul_fmpz(fmprb_ptr res, fmprb_srcptr vec, - slong len, const fmpz_t c, slong prec) -{ - slong i; - fmpr_t t; - fmpr_init(t); - fmpr_set_fmpz(t, c); - for (i = 0; i < len; i++) - fmprb_mul_fmpr(res + i, vec + i, t, prec); - fmpr_clear(t); -} - -static __inline__ void -_fmprb_vec_scalar_mul_2exp_si(fmprb_ptr res, fmprb_srcptr src, slong len, slong c) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_mul_2exp_si(res + i, src + i, c); -} - -static __inline__ void -_fmprb_vec_scalar_addmul(fmprb_ptr res, fmprb_srcptr vec, - slong len, const fmprb_t c, slong prec) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_addmul(res + i, vec + i, c, prec); -} - -static __inline__ void -_fmprb_vec_get_abs_ubound_fmpr(fmpr_t bound, fmprb_srcptr vec, - slong len, slong prec) -{ - fmpr_t t; - slong i; - - if (len < 1) - { - fmpr_zero(bound); - } - else - { - fmprb_get_abs_ubound_fmpr(bound, vec, prec); - fmpr_init(t); - for (i = 1; i < len; i++) - { - fmprb_get_abs_ubound_fmpr(t, vec + i, prec); - if (fmpr_cmp(t, bound) > 0) - fmpr_set(bound, t); - } - fmpr_clear(t); - } -} - -static __inline__ slong -_fmprb_vec_bits(fmprb_srcptr x, slong len) -{ - slong i, b, c; - - b = 0; - for (i = 0; i < len; i++) - { - c = fmprb_bits(x + i); - b = FLINT_MAX(b, c); - } - - return b; -} - -static __inline__ void -_fmprb_vec_set_powers(fmprb_ptr xs, const fmprb_t x, slong len, slong prec) -{ - slong i; - - for (i = 0; i < len; i++) - { - if (i == 0) - fmprb_one(xs + i); - else if (i == 1) - fmprb_set_round(xs + i, x, prec); - else if (i % 2 == 0) - fmprb_mul(xs + i, xs + i / 2, xs + i / 2, prec); - else - fmprb_mul(xs + i, xs + i - 1, x, prec); - } -} - -static __inline__ void -_fmprb_vec_add_error_fmpr_vec(fmprb_ptr res, fmpr_srcptr err, slong len) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_add_error_fmpr(res + i, err + i); -} - -static __inline__ void -_fmprb_vec_indeterminate(fmprb_ptr vec, slong len) -{ - slong i; - for (i = 0; i < len; i++) - { - fmpr_nan(fmprb_midref(vec + i)); - fmpr_pos_inf(fmprb_radref(vec + i)); - } -} - -static __inline__ void -_fmprb_vec_trim(fmprb_ptr res, fmprb_srcptr vec, slong len) -{ - slong i; - for (i = 0; i < len; i++) - fmprb_trim(res + i, vec + i); -} - #ifdef __cplusplus } #endif diff --git a/fmprb/agm.c b/fmprb/agm.c deleted file mode 100644 index 79625142..00000000 --- a/fmprb/agm.c +++ /dev/null @@ -1,81 +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" - -void -fmprb_agm(fmprb_t z, const fmprb_t x, const fmprb_t y, slong prec) -{ - fmprb_t t, u, v, w; - - if (fmprb_contains_negative(x) || fmprb_contains_negative(y)) - { - fmprb_indeterminate(z); - return; - } - - if (fmprb_is_zero(x) || fmprb_is_zero(y)) - { - fmprb_zero(z); - return; - } - - fmprb_init(t); - fmprb_init(u); - fmprb_init(v); - fmprb_init(w); - - fmprb_set(t, x); - fmprb_set(u, y); - - while (!fmprb_overlaps(t, u) && - !fmprb_contains_nonpositive(t) && - !fmprb_contains_nonpositive(u)) - { - fmprb_add(v, t, u, prec); - fmprb_mul_2exp_si(v, v, -1); - - fmprb_mul(w, t, u, prec); - fmprb_sqrt(w, w, prec); - - fmprb_swap(v, t); - fmprb_swap(w, u); - } - - if (!fmprb_is_finite(t) || !fmprb_is_finite(u)) - { - fmprb_indeterminate(z); - } - else - { - fmprb_union(z, t, u, prec); - } - - fmprb_clear(t); - fmprb_clear(u); - fmprb_clear(v); - fmprb_clear(w); -} - diff --git a/fmprb/root.c b/fmprb/root.c deleted file mode 100644 index 86b11f62..00000000 --- a/fmprb/root.c +++ /dev/null @@ -1,102 +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) 2012, 2013 Fredrik Johansson - -******************************************************************************/ - -#include "fmprb.h" - -void -fmprb_root(fmprb_t z, const fmprb_t x, ulong k, slong prec) -{ - slong r; - - if (k == 1) - { - fmprb_set_round(z, x, prec); - return; - } - - if (k == 2) - { - fmprb_sqrt(z, x, prec); - return; - } - - if (k == -2) - { - fmprb_rsqrt(z, x, prec); - return; - } - - if (fmprb_contains_negative(x)) - { - fmprb_indeterminate(z); - return; - } - - if (fmprb_is_exact(x)) - { - r = fmpr_root(fmprb_midref(z), fmprb_midref(x), k, prec, FMPR_RND_DOWN); - fmpr_set_error_result(fmprb_radref(z), fmprb_midref(z), r); - } - else - { - fmpr_t t, u, err; - - fmpr_init(t); - fmpr_init(u); - fmpr_init(err); - - /* lower point */ - fmpr_sub(err, fmprb_midref(x), fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_DOWN); - - if (fmpr_is_zero(err)) - { - fmpr_root(err, fmprb_midref(x), k, FMPRB_RAD_PREC, FMPR_RND_UP); - r = fmpr_root(fmprb_midref(z), fmprb_midref(x), k, prec, FMPR_RND_DOWN); - fmpr_add_error_result(fmprb_radref(z), err, fmprb_midref(z), r, - FMPRB_RAD_PREC, FMPR_RND_UP); - } - else - { - /* derivative x^(1/k) / (x k) at lower point */ - fmpr_root(t, err, k, FMPRB_RAD_PREC, FMPR_RND_UP); - fmpr_mul_ui(u, err, k, FMPRB_RAD_PREC, FMPR_RND_DOWN); - fmpr_divappr_abs_ubound(t, t, u, FMPRB_RAD_PREC); - - /* multiplied by distance */ - fmpr_mul(err, t, fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_UP); - - r = fmpr_root(fmprb_midref(z), fmprb_midref(x), k, prec, FMPR_RND_DOWN); - fmpr_add_error_result(fmprb_radref(z), err, fmprb_midref(z), r, - FMPRB_RAD_PREC, FMPR_RND_UP); - } - - fmpr_clear(t); - fmpr_clear(u); - fmpr_clear(err); - } - - fmprb_adjust(z); -} - diff --git a/fmprb/test/t-agm.c b/fmprb/test/t-agm.c deleted file mode 100644 index 9008ab9d..00000000 --- a/fmprb/test/t-agm.c +++ /dev/null @@ -1,102 +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) 2012 Fredrik Johansson - -******************************************************************************/ - -#include "fmprb.h" - -int main() -{ - slong iter; - flint_rand_t state; - - flint_printf("agm...."); - fflush(stdout); - - flint_randinit(state); - - for (iter = 0; iter < 100000; iter++) - { - fmprb_t a, b, c; - fmpq_t q, r; - mpfr_t t, u; - slong prec = 2 + n_randint(state, 200); - - fmprb_init(a); - fmprb_init(b); - fmprb_init(c); - fmpq_init(q); - fmpq_init(r); - mpfr_init2(t, prec + 100); - mpfr_init2(u, prec + 100); - - fmprb_randtest(a, state, 1 + n_randint(state, 200), 3); - fmprb_randtest(b, state, 1 + n_randint(state, 200), 3); - fmprb_randtest(c, state, 1 + n_randint(state, 200), 3); - - fmprb_agm(c, a, b, prec); - - if (fmprb_equal(a, b)) - { - if (!fmprb_contains(c, a)) - { - flint_printf("FAIL: containment (identity)\n\n"); - flint_printf("a = "); fmprb_print(a); flint_printf("\n\n"); - flint_printf("b = "); fmprb_print(b); flint_printf("\n\n"); - flint_printf("c = "); fmprb_print(c); flint_printf("\n\n"); - abort(); - } - } - else - { - fmprb_get_rand_fmpq(q, state, a, 1 + n_randint(state, 200)); - fmprb_get_rand_fmpq(r, state, b, 1 + n_randint(state, 200)); - fmpq_get_mpfr(t, q, MPFR_RNDN); - fmpq_get_mpfr(u, r, MPFR_RNDN); - mpfr_agm(t, t, u, MPFR_RNDN); - - if (!fmprb_contains_mpfr(c, t)) - { - flint_printf("FAIL: containment\n\n"); - flint_printf("a = "); fmprb_print(a); flint_printf("\n\n"); - flint_printf("b = "); fmprb_print(b); flint_printf("\n\n"); - flint_printf("c = "); fmprb_print(c); flint_printf("\n\n"); - abort(); - } - } - - fmprb_clear(a); - fmprb_clear(b); - fmprb_clear(c); - fmpq_clear(q); - fmpq_clear(r); - mpfr_clear(t); - mpfr_clear(u); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} - diff --git a/fmprb/test/t-trim.c b/fmprb/test/t-trim.c deleted file mode 100644 index 6cc33f25..00000000 --- a/fmprb/test/t-trim.c +++ /dev/null @@ -1,91 +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() -{ - slong iter; - flint_rand_t state; - - flint_printf("trim...."); - fflush(stdout); - flint_randinit(state); - - for (iter = 0; iter < 100000; iter++) - { - fmprb_t x, y; - slong acc1, acc2; - int accuracy_ok; - - fmprb_init(x); - fmprb_init(y); - - fmprb_randtest_special(x, state, 1000, 100); - - fmprb_trim(y, x); - - if (!fmprb_contains(y, x)) - { - flint_printf("FAIL (containment):\n\n"); - flint_printf("x = "); fmprb_print(x); flint_printf("\n\n"); - flint_printf("y = "); fmprb_print(y); flint_printf("\n\n"); - abort(); - } - - acc1 = fmprb_rel_accuracy_bits(x); - acc2 = fmprb_rel_accuracy_bits(y); - - accuracy_ok = (acc1 < 0 && acc2 <= acc1) || (acc2 >= acc1 - 1); - - if (!accuracy_ok) - { - flint_printf("FAIL (accuracy):\n\n"); - flint_printf("x: %wd, y = %wd\n\n", acc1, acc2); - flint_printf("x = "); fmprb_print(x); flint_printf("\n\n"); - flint_printf("y = "); fmprb_print(y); flint_printf("\n\n"); - abort(); - } - - fmprb_trim(x, x); - - if (!fmprb_equal(y, x)) - { - flint_printf("FAIL (aliasing):\n\n"); - flint_printf("x = "); fmprb_print(x); flint_printf("\n\n"); - flint_printf("y = "); fmprb_print(y); flint_printf("\n\n"); - abort(); - } - - fmprb_clear(x); - fmprb_clear(y); - } - - flint_randclear(state); - flint_cleanup(); - flint_printf("PASS\n"); - return EXIT_SUCCESS; -} - diff --git a/fmprb/trim.c b/fmprb/trim.c deleted file mode 100644 index 256dad0d..00000000 --- a/fmprb/trim.c +++ /dev/null @@ -1,68 +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" - -#define TRIM_PADDING 16 - -void -fmprb_trim(fmprb_t y, const fmprb_t x) -{ - if (fmpr_is_zero(fmprb_radref(x)) || fmpr_is_special(fmprb_midref(x))) - { - fmprb_set(y, x); - } - else if (fmpr_is_special(fmprb_radref(x))) - { - /* midpoint must be finite, so set to 0 +/- inf */ - fmprb_zero_pm_inf(y); - } - else - { - slong bits, accuracy; - - bits = fmprb_bits(x); - accuracy = fmprb_rel_accuracy_bits(x); - - if (accuracy < -TRIM_PADDING) - { - /* set to 0 +/- rad */ - if (fmpr_sgn(fmprb_midref(x)) < 0) - fmpr_sub(fmprb_radref(y), fmprb_radref(x), fmprb_midref(x), FMPRB_RAD_PREC, FMPR_RND_UP); - else - fmpr_add(fmprb_radref(y), fmprb_radref(x), fmprb_midref(x), FMPRB_RAD_PREC, FMPR_RND_UP); - fmpr_zero(fmprb_midref(y)); - } - else if (accuracy < bits - TRIM_PADDING) - { - fmprb_set_round(y, x, FLINT_MAX(0, accuracy) + TRIM_PADDING); - } - else - { - fmprb_set(y, x); - } - } -} -