From 0998aa9821414be99a2d871aed696df271381da0 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 8 May 2014 14:48:36 +0200 Subject: [PATCH] more code and test code --- arb/get_mag_lbound.c | 112 +++++++++++++++++++++ arf.h | 140 ++++++++++++--------------- mag/set_arf.c => arf/get_mag.c | 18 ++-- doc/source/arb.rst | 7 ++ doc/source/arf.rst | 7 ++ doc/source/index.rst | 10 ++ doc/source/mag.rst | 122 +++++++++++++++++++++++ mag.h | 172 ++++++++++++++++++++------------- mag/add_2exp_fmpz.c | 8 +- mag/addmul.c | 4 +- mag/div.c | 58 +++++++++++ mag/mul.c | 2 +- mag/set_fmpr.c | 3 +- mag/test/t-add_2exp_fmpz.c | 100 +++++++++++++++++++ mag/test/t-addmul.c | 102 +++++++++++++++++++ mag/test/t-fast_add_2exp_si.c | 95 ++++++++++++++++++ mag/test/t-fast_addmul.c | 102 +++++++++++++++++++ mag/test/t-fast_mul.c | 95 ++++++++++++++++++ mag/test/t-mul.c | 95 ++++++++++++++++++ 19 files changed, 1092 insertions(+), 160 deletions(-) create mode 100644 arb/get_mag_lbound.c rename mag/set_arf.c => arf/get_mag.c (79%) create mode 100644 doc/source/arb.rst create mode 100644 doc/source/arf.rst create mode 100644 doc/source/mag.rst create mode 100644 mag/div.c create mode 100644 mag/test/t-add_2exp_fmpz.c create mode 100644 mag/test/t-addmul.c create mode 100644 mag/test/t-fast_add_2exp_si.c create mode 100644 mag/test/t-fast_addmul.c create mode 100644 mag/test/t-fast_mul.c create mode 100644 mag/test/t-mul.c diff --git a/arb/get_mag_lbound.c b/arb/get_mag_lbound.c new file mode 100644 index 00000000..58c281de --- /dev/null +++ b/arb/get_mag_lbound.c @@ -0,0 +1,112 @@ +/*============================================================================= + + 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" + +/* + +arb_get_mag_upper +arb_get_mag_lower + +arf_get_mag_upper +arf_get_mag_lower + + + +max(0, |mid|-rad) + +*/ + +static __inline__ void +_arb_get_mag_lower(mag_t z, const arf_t mid, const mag_t rad) +{ + if (mag_is_inf(rad) || arf_is_zero(mid)) + { + mag_zero(z); + } + else if (mag_is_zero(rad) || !arf_is_finite(mid)) + { + abort(); + /* arf_get_mag_lower(z, mid); */ + } + else + { + long shift, fix; + + shift = _fmpz_sub_small(MAG_EXPREF(mid), MAG_EXPREF(rad)); + + /* mid < rad */ + if (shift < 0) + { + mag_zero(z); + } + else if (shift <= 1) /* can be cancellation */ + { + arf_t t; + arf_init(t); + + arf_set_mag(t, rad); + + if (arf_sgn(mid) > 0) + { + arf_sub(t, mid, t, MAG_BITS, ARF_RND_DOWN); + } + else + { + arf_add(t, mid, t, MAG_BITS, ARF_RND_DOWN); + arf_neg(t, t); + } + + if (arf_sgn(t) <= 0) + mag_zero(z); + else + abort(); /* arf_get_mag_lower(z, t); */ + + arf_clear(t); + } + else + { + mp_limb_t m; + + ARF_GET_TOP_LIMB(m, mid); + + if (shift <= MAG_BITS) + m = m - (MAG_MAN(rad) >> shift) - 1; + else + m = m - 1; + + fix = !(m >> (MAG_BITS - 1)); + m <<= fix; + _fmpz_add_fast(MAG_EXPREF(z), MAG_EXPREF(mid), -fix); + } + } +} + +void +arb_get_mag_lower(mag_t z, const arb_t x) +{ + _arb_get_mag_lower(z, ARB_MIDREF(x), ARB_RADREF(x)); +} + diff --git a/arf.h b/arf.h index 5620ae53..d8c3ea3d 100644 --- a/arf.h +++ b/arf.h @@ -32,17 +32,12 @@ #include #include "flint.h" #include "fmpr.h" +#include "mag.h" #ifdef __cplusplus extern "C" { #endif -#define LIMB_ONE ((mp_limb_t) 1) -#define LIMB_ONES (-(mp_limb_t) 1) -#define LIMB_TOP (((mp_limb_t) 1) << (FLINT_BITS - 1)) - -#define MASK_LIMB(n, c) ((n) & (LIMB_ONES << (c))) - #define arf_rnd_t fmpr_rnd_t #define ARF_RND_DOWN FMPR_RND_DOWN #define ARF_RND_UP FMPR_RND_UP @@ -87,8 +82,8 @@ arf_rnd_to_mpfr(arf_rnd_t rnd) #define ARF_RESULT_INEXACT 1 /* Range where we can skip fmpz overflow checks for exponent manipulation. */ -#define ARF_MAX_LAGOM_EXP (COEFF_MAX / 4) -#define ARF_MIN_LAGOM_EXP (-ARF_MAX_LAGOM_EXP) +#define ARF_MAX_LAGOM_EXP MAG_MAX_LAGOM_EXP +#define ARF_MIN_LAGOM_EXP MAG_MIN_LAGOM_EXP /* Exponent values used to encode special values. */ #define ARF_EXP_ZERO 0 @@ -192,6 +187,15 @@ typedef const arf_struct * arf_srcptr; xptr = ARF_PTR_D(x); \ } while (0) +/* Assumes non-special! */ +#define ARF_GET_TOP_LIMB(lmb, x) \ + do { \ + mp_srcptr __xptr; \ + mp_size_t __xn; \ + ARF_GET_MPN_READONLY(__xptr, __xn, (x)); \ + (lmb) = __xptr[__xn - 1]; \ + } while (0) + /* Get mpn pointer xptr for writing *exactly* xn limbs to x. */ #define ARF_GET_MPN_WRITE(xptr, xn, x) \ do { \ @@ -742,74 +746,6 @@ void arf_randtest_not_zero(arf_t x, flint_rand_t state, long bits, long mag_bits void arf_randtest_special(arf_t x, flint_rand_t state, long bits, long mag_bits); -#define ADD2_FAST_MAX (COEFF_MAX / 4) -#define ADD2_FAST_MIN (-ADD2_FAST_MAX) - -static __inline__ void -_fmpz_set_fast(fmpz_t f, const fmpz_t g) -{ - if (!COEFF_IS_MPZ(*f) && !COEFF_IS_MPZ(*g)) - *f = *g; - else - fmpz_set(f, g); -} - -static __inline__ void -_fmpz_add_fast(fmpz_t z, const fmpz_t x, long c) -{ - fmpz ze, xe; - - ze = *z; - xe = *x; - - if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX)) - *z = xe + c; - else - fmpz_add_si(z, x, c); -} - -static __inline__ void -_fmpz_add2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c) -{ - fmpz ze, xe, ye; - - ze = *z; - xe = *x; - ye = *y; - - if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) && - (ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX)) - { - *z = xe + ye + c; - } - else - { - fmpz_add(z, x, y); - fmpz_add_si(z, z, c); - } -} - -static __inline__ void -_fmpz_sub2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c) -{ - fmpz ze, xe, ye; - - ze = *z; - xe = *x; - ye = *y; - - if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) && - (ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX)) - { - *z = xe - ye + c; - } - else - { - fmpz_sub(z, x, y); - fmpz_add_si(z, z, c); - } -} - #define MUL_MPFR_MIN_LIMBS 25 #define MUL_MPFR_MAX_LIMBS 10000 @@ -1143,6 +1079,58 @@ int arf_sqrt_fmpz(arf_t z, const fmpz_t x, long prec, arf_rnd_t rnd); int arf_rsqrt(arf_ptr z, arf_srcptr x, long prec, arf_rnd_t rnd); +/* Magnitude bounds */ + +void arf_get_mag(mag_t y, const arf_t x); + +static __inline__ void +arf_set_mag(arf_t y, const mag_t x) +{ + if (mag_is_zero(x)) + { + arf_zero(y); + } + else if (mag_is_inf(x)) + { + arf_pos_inf(y); + } + else + { + _fmpz_set_fast(ARF_EXPREF(y), MAG_EXPREF(x)); + ARF_DEMOTE(y); + ARF_XSIZE(y) = ARF_MAKE_XSIZE(1, 0); + ARF_NOPTR_D(y)[0] = MAG_MAN(x) << (FLINT_BITS - MAG_BITS); + } +} + +static __inline__ void +mag_init_set_arf(mag_t y, const arf_t x) +{ + mag_init(y); + arf_get_mag(y, x); +} + +static __inline__ void +mag_fast_init_set_arf(mag_t y, const arf_t x) +{ + if (ARF_IS_SPECIAL(x)) /* x == 0 */ + { + mag_fast_zero(y); + } + else + { + mp_srcptr xp; + mp_size_t xn; + + ARF_GET_MPN_READONLY(xp, xn, x); + + MAG_MAN(y) = (xp[xn - 1] >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE; + MAG_EXP(y) = ARF_EXP(x); + + MAG_FAST_ADJUST_ONE_TOO_LARGE(y); + } +} + #ifdef __cplusplus } #endif diff --git a/mag/set_arf.c b/arf/get_mag.c similarity index 79% rename from mag/set_arf.c rename to arf/get_mag.c index 87aef4da..e0022140 100644 --- a/mag/set_arf.c +++ b/arf/get_mag.c @@ -23,14 +23,11 @@ ******************************************************************************/ -#include "mag.h" +#include "arf.h" void -mag_set_arf(mag_t y, const arf_t x) +arf_get_mag(mag_t y, const arf_t x) { - mp_srcptr xp; - mp_size_t xn; - if (arf_is_zero(x)) { mag_zero(y); @@ -41,11 +38,14 @@ mag_set_arf(mag_t y, const arf_t x) } else { - ARF_GET_MPN_READONLY(xp, xn, x); + mp_limb_t t, u; - MAG_MAN(y) = (xp[xn - 1] >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE; - _fmpz_set_fast(MAG_EXPREF(y), ARF_EXPREF(x)); - MAG_ADJUST_ONE_TOO_LARGE(y); /* TODO: combine */ + ARF_GET_TOP_LIMB(t, x); + t = (t >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE; + u = t >> MAG_BITS; + t = (t >> u) + u; + _fmpz_add_fast(MAG_EXPREF(y), ARF_EXPREF(x), u); + MAG_MAN(y) = t; } } diff --git a/doc/source/arb.rst b/doc/source/arb.rst new file mode 100644 index 00000000..a25bfc08 --- /dev/null +++ b/doc/source/arb.rst @@ -0,0 +1,7 @@ +.. _arb: + +**arb.h** -- arbitrary-precision real floating-point balls +=============================================================================== + +TBD. + diff --git a/doc/source/arf.rst b/doc/source/arf.rst new file mode 100644 index 00000000..96c33dbc --- /dev/null +++ b/doc/source/arf.rst @@ -0,0 +1,7 @@ +.. _arf: + +**arf.h** -- arbitrary-precision real floating-point numbers +=============================================================================== + +TBD. + diff --git a/doc/source/index.rst b/doc/source/index.rst index 12e46aaf..849fc12e 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -47,6 +47,16 @@ Module documentation hypgeom.rst partitions.rst +Module documentation (Arb 2.x types) +:::::::::::::::::::::::::::::::::::: + +.. toctree:: + :maxdepth: 1 + + mag.rst + arf.rst + arb.rst + Credits and references :::::::::::::::::::::::: diff --git a/doc/source/mag.rst b/doc/source/mag.rst new file mode 100644 index 00000000..e5b9bec1 --- /dev/null +++ b/doc/source/mag.rst @@ -0,0 +1,122 @@ +.. _mag: + +**mag.h** -- fixed-precision unsigned floating-point numbers for bounds +=============================================================================== + +The :type:`mag_t` type is an unsigned floating-point type with a +fixed-precision mantissa (30 bits) and unlimited exponent range, suited for +representing and rigorously manipulating magnitude bounds efficiently. +Operations always produce a strict upper or lower bound, but for performance +reasons, no attempt is made to compute the best possible bound +(in general, a result may a few ulps larger/smaller than the optimal value). +The special values zero and positive infinity are supported (but not NaN). +Applications requiring more flexibility (such as correct rounding, or +higher precision) should use the :type:`arf_t` type instead. + +The exponent is represented as an :type:`fmpz`. Special "fast" methods are +provided for manipulating :type:`mag_t` objects under the assumption that +no infinities are present and that exponents stay within a certain range. + +Types, macros and constants +------------------------------------------------------------------------------- + +.. type:: mag_struct + + A :type:`mag_struct` holds a mantissa and an exponent. + Special values are encoded by the mantissa being set to zero. + +.. type:: mag_t + + A :type:`mag_t` is defined as an array of length one of type + :type:`mag_struct`, permitting a :type:`mag_t` to be passed by reference. + +Memory management +------------------------------------------------------------------------------- + +.. function:: void mag_init(mag_t x) + + Initializes the variable *x* for use. Its value is set to zero. + +.. function:: void mag_clear(mag_t x) + + Clears the variable *x*, freeing or recycling its allocated memory. + +.. function:: void mag_init_set(mag_t x, const mag_t y) + +.. function:: void mag_swap(mag_t x, mag_t y) + +.. function:: void mag_set(mag_t x, const mag_t y) + +Special values +------------------------------------------------------------------------------- + +.. function:: void mag_zero(mag_t x) + + Sets *x* to zero. + +.. function:: void mag_inf(mag_t x) + + Sets *x* to positive infinity. + +.. function:: int mag_is_special(const mag_t x) + + Returns nonzero iff *x* is zero or positive infinity. + +.. function:: int mag_is_zero(const mag_t x) + + Returns nonzero iff *x* is zero. + +.. function:: int mag_is_inf(const mag_t x) + + Returns nonzero iff *x* is positive infinity. + +Arithmetic +------------------------------------------------------------------------------- + +.. function:: void mag_mul(mag_t z, const mag_t x, const mag_t y) + + Sets `z` to an upper bound for `x + y`. + +.. function:: void mag_addmul(mag_t z, const mag_t x, const mag_t y) + + Sets `z` to an upper bound for `z + xy`. + +.. function:: void mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e) + + Sets `z` to an upper bound for `x + 2^e`. + +.. function:: void mag_div(mag_t z, const mag_t x, const mag_t y) + + Sets `z` to an upper bound for `x / y`. + +Fast versions +------------------------------------------------------------------------------- + +.. function:: void mag_fast_init_set(mag_t x, const mag_t y) + +.. function:: void mag_fast_zero(mag_t x) + +.. function:: int mag_fast_is_zero(const mag_t x) + +.. function:: void mag_fast_init_set_arf(mag_t y, const arf_t x) + +.. function:: void mag_fast_mul(mag_t z, const mag_t x, const mag_t y) + +.. function:: void mag_fast_addmul(mag_t z, const mag_t x, const mag_t y) + +.. function:: void mag_fast_add_2exp_si(mag_t z, const mag_t x, long e) + +Conversions +------------------------------------------------------------------------------- + +These functions are intended for debugging purposes: see the :doc:`arf ` +module for other conversion functions. + +.. function:: void mag_set_fmpr(mag_t y, const fmpr_t x) + + Sets *y* to an upper bound for *x*. + +.. function:: void mag_get_fmpr(fmpr_t y, const mag_t x) + + Sets *y* to exactly *x*. + diff --git a/mag.h b/mag.h index b49c3616..e66c7c06 100644 --- a/mag.h +++ b/mag.h @@ -28,37 +28,102 @@ #include #include "flint.h" -#include "arf.h" +#include "fmpz.h" +#include "fmpz_extras.h" #ifdef __cplusplus extern "C" { #endif -/* -The mag_t type is an unsigned floating-point type with a fixed-precision -mantissa (30 bits) and unlimited exponent range, suited for representing -magnitude bounds efficiently. +#define LIMB_ONE ((mp_limb_t) 1) +#define LIMB_ONES (-(mp_limb_t) 1) +#define LIMB_TOP (((mp_limb_t) 1) << (FLINT_BITS - 1)) +#define MASK_LIMB(n, c) ((n) & (LIMB_ONES << (c))) -Operations always produce a strict upper/lower bound, but for performance -reasons, no attempt is made to compute the best possible bound -(in general, a result may a few ulps larger/smaller than the optimal value). +#define MAG_MAX_LAGOM_EXP (COEFF_MAX / 4) +#define MAG_MIN_LAGOM_EXP (-MAG_MAX_LAGOM_EXP) -The special values zero and positive infinity are supported (but not NaN). +#define ADD2_FAST_MAX (COEFF_MAX / 4) +#define ADD2_FAST_MIN (-ADD2_FAST_MAX) -Applications requiring more flexibility (such as correct rounding, or -higher precision) should use the arf_t type instead. -*/ +/* TODO: rename these and move them to fmpz_extras */ -/* TODO: arf.h should depend on this, not the other way around */ +static __inline__ void +_fmpz_set_fast(fmpz_t f, const fmpz_t g) +{ + if (!COEFF_IS_MPZ(*f) && !COEFF_IS_MPZ(*g)) + *f = *g; + else + fmpz_set(f, g); +} + +static __inline__ void +_fmpz_add_fast(fmpz_t z, const fmpz_t x, long c) +{ + fmpz ze, xe; + + ze = *z; + xe = *x; + + if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX)) + *z = xe + c; + else + fmpz_add_si(z, x, c); +} + +static __inline__ void +_fmpz_add2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c) +{ + fmpz ze, xe, ye; + + ze = *z; + xe = *x; + ye = *y; + + if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) && + (ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX)) + { + *z = xe + ye + c; + } + else + { + fmpz_add(z, x, y); + fmpz_add_si(z, z, c); + } +} + +static __inline__ void +_fmpz_sub2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c) +{ + fmpz ze, xe, ye; + + ze = *z; + xe = *x; + ye = *y; + + if (!COEFF_IS_MPZ(ze) && (xe > ADD2_FAST_MIN && xe < ADD2_FAST_MAX) && + (ye > ADD2_FAST_MIN && ye < ADD2_FAST_MAX)) + { + *z = xe - ye + c; + } + else + { + fmpz_sub(z, x, y); + fmpz_add_si(z, z, c); + } +} + + +#define MAG_EXP_POS_INF (COEFF_MIN+1) + +/* Finite and with lagom big exponents. */ +#define MAG_IS_LAGOM(x) (MAG_EXP(x) >= MAG_MIN_LAGOM_EXP && \ + MAG_EXP(x) <= MAG_MAX_LAGOM_EXP) #define MAG_EXPREF(x) (&(x)->exp) #define MAG_EXP(x) ((x)->exp) #define MAG_MAN(x) ((x)->man) -/* Finite and with lagom big exponents. */ -#define MAG_IS_LAGOM(x) (MAG_EXP(x) >= ARF_MIN_LAGOM_EXP && \ - MAG_EXP(x) <= ARF_MAX_LAGOM_EXP) - #define MAG_BITS 30 #define MAG_ONE_HALF (1UL << (MAG_BITS - 1)) @@ -155,7 +220,7 @@ mag_swap(mag_t x, mag_t y) static __inline__ void mag_set(mag_t x, const mag_t y) { - _fmpz_set_fast(ARF_EXPREF(x), ARF_EXPREF(y)); + _fmpz_set_fast(MAG_EXPREF(x), MAG_EXPREF(y)); x->man = y->man; } @@ -166,13 +231,6 @@ mag_zero(mag_t x) MAG_MAN(x) = 0; } -static __inline__ void -mag_inf(mag_t x) -{ - fmpz_clear(MAG_EXPREF(x)); - MAG_EXP(x) = ARF_EXP_POS_INF; -} - static __inline__ int mag_is_special(const mag_t x) { @@ -185,6 +243,13 @@ mag_is_zero(const mag_t x) return (MAG_MAN(x) == 0) && (MAG_EXP(x) == 0); } +static __inline__ void +mag_inf(mag_t x) +{ + fmpz_clear(MAG_EXPREF(x)); + MAG_EXP(x) = MAG_EXP_POS_INF; +} + static __inline__ int mag_is_inf(const mag_t x) { @@ -193,23 +258,14 @@ mag_is_inf(const mag_t x) /* general versions */ -void mag_set_fmpr(mag_t x, const fmpr_t y); - -void mag_set_arf(mag_t y, const arf_t x); - -static __inline__ void -mag_init_set_arf(mag_t y, const arf_t x) -{ - mag_init(y); - mag_set_arf(y, x); -} - void mag_mul(mag_t z, const mag_t x, const mag_t y); void mag_addmul(mag_t z, const mag_t x, const mag_t y); void mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e); +void mag_div(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! */ @@ -233,27 +289,6 @@ mag_fast_is_zero(const mag_t x) return MAG_MAN(x) == 0; } -static __inline__ void -mag_fast_init_set_arf(mag_t y, const arf_t x) -{ - if (ARF_IS_SPECIAL(x)) /* x == 0 */ - { - mag_fast_zero(y); - } - else - { - mp_srcptr xp; - mp_size_t xn; - - ARF_GET_MPN_READONLY(xp, xn, x); - - MAG_MAN(y) = (xp[xn - 1] >> (FLINT_BITS - MAG_BITS)) + LIMB_ONE; - MAG_EXP(y) = ARF_EXP(x); - - MAG_FAST_ADJUST_ONE_TOO_LARGE(y); - } -} - static __inline__ void mag_fast_mul(mag_t z, const mag_t x, const mag_t y) { @@ -293,8 +328,7 @@ mag_fast_addmul(mag_t z, const mag_t x, const mag_t y) if (shift >= MAG_BITS) MAG_MAN(z)++; else - MAG_MAN(z) = MAG_MAN(z) + (MAG_FIXMUL(MAG_MAN(x), - MAG_MAN(y)) >> shift) + LIMB_ONE; + MAG_MAN(z) = MAG_MAN(z) + (MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) >> shift) + 1; } else { @@ -302,11 +336,11 @@ mag_fast_addmul(mag_t z, const mag_t x, const mag_t y) MAG_EXP(z) = e; if (shift >= MAG_BITS) - MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) - + (2 * LIMB_ONE); + MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) + 2; else - MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) - + (MAG_MAN(z) >> shift) + (2 * LIMB_ONE); + MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) + (MAG_MAN(z) >> shift) + 2; + + MAG_FAST_ADJUST_ONE_TOO_SMALL(z); } MAG_FAST_ADJUST_ONE_TOO_LARGE(z); @@ -327,7 +361,7 @@ mag_fast_add_2exp_si(mag_t z, const mag_t x, long e) long shift; shift = MAG_EXP(x) - e; - if (shift >= 0) + if (shift > 0) { MAG_EXP(z) = MAG_EXP(x); @@ -340,18 +374,21 @@ mag_fast_add_2exp_si(mag_t z, const mag_t x, long e) { shift = -shift; - MAG_EXP(z) = e; + MAG_EXP(z) = e + 1; if (shift >= MAG_BITS) - MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + LIMB_ONE; + MAG_MAN(z) = MAG_ONE_HALF + LIMB_ONE; else - MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + (MAG_MAN(x) >> shift); + MAG_MAN(z) = MAG_ONE_HALF + (MAG_MAN(x) >> (shift + 1)) + LIMB_ONE; } MAG_FAST_ADJUST_ONE_TOO_LARGE(z); } } +#include "fmpr.h" + +void mag_set_fmpr(mag_t x, const fmpr_t y); static __inline__ void mag_get_fmpr(fmpr_t x, const mag_t r) @@ -371,7 +408,6 @@ mag_get_fmpr(fmpr_t x, const mag_t r) } } - #ifdef __cplusplus } #endif diff --git a/mag/add_2exp_fmpz.c b/mag/add_2exp_fmpz.c index 7a635905..4934168b 100644 --- a/mag/add_2exp_fmpz.c +++ b/mag/add_2exp_fmpz.c @@ -46,7 +46,7 @@ mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e) shift = _fmpz_sub_small(MAG_EXPREF(x), e); - if (shift >= 0) + if (shift > 0) { _fmpz_set_fast(MAG_EXPREF(z), MAG_EXPREF(x)); @@ -59,12 +59,12 @@ mag_add_2exp_fmpz(mag_t z, const mag_t x, const fmpz_t e) { shift = -shift; - _fmpz_set_fast(MAG_EXPREF(z), e); + _fmpz_add_fast(MAG_EXPREF(z), e, 1); if (shift >= MAG_BITS) - MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + LIMB_ONE; + MAG_MAN(z) = MAG_ONE_HALF + LIMB_ONE; else - MAG_MAN(z) = (LIMB_ONE << MAG_BITS) + (MAG_MAN(x) >> shift); + MAG_MAN(z) = MAG_ONE_HALF + (MAG_MAN(x) >> (shift + 1)) + LIMB_ONE; } MAG_ADJUST_ONE_TOO_LARGE(z); diff --git a/mag/addmul.c b/mag/addmul.c index eeae7f37..146062b4 100644 --- a/mag/addmul.c +++ b/mag/addmul.c @@ -38,7 +38,7 @@ mag_addmul(mag_t z, const mag_t x, const mag_t y) } else if (mag_is_zero(x) || mag_is_zero(y)) { - mag_zero(z); + return; } else { @@ -70,6 +70,8 @@ mag_addmul(mag_t z, const mag_t x, const mag_t y) else MAG_MAN(z) = MAG_FIXMUL(MAG_MAN(x), MAG_MAN(y)) + (MAG_MAN(z) >> shift) + (2 * LIMB_ONE); + + MAG_ADJUST_ONE_TOO_SMALL(z); } MAG_ADJUST_ONE_TOO_LARGE(z); diff --git a/mag/div.c b/mag/div.c new file mode 100644 index 00000000..73695ee5 --- /dev/null +++ b/mag/div.c @@ -0,0 +1,58 @@ +/*============================================================================= + + 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" + +void +mag_div(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_is_inf(x)) + mag_inf(z); + else + mag_zero(z); + } + else + { + mp_limb_t q; + long fix; + +#if FLINT_BITS == 64 + q = (MAG_MAN(x) << MAG_BITS) / MAG_MAN(y) + 1; +#else + mp_limb_t hi, lo, r; + lo = MAG_MAN(x) << MAG_BITS; + hi = MAG_MAN(x) >> (FLINT_BITS - MAG_BITS); + udiv_qrnnd(q, r, hi, lo, MAG_MAN(y)); + q += 1; +#endif + + fix = q >> MAG_BITS; + q = (q >> fix) + fix; + _fmpz_sub2_fast(MAG_EXPREF(z), MAG_EXPREF(x), MAG_EXPREF(y), fix); + } +} + diff --git a/mag/mul.c b/mag/mul.c index e7af1575..060fbf29 100644 --- a/mag/mul.c +++ b/mag/mul.c @@ -30,7 +30,7 @@ mag_mul(mag_t z, const mag_t x, const mag_t y) { if (mag_is_special(x) || mag_is_special(y)) { - if (mag_is_inf(x) || mag_is_special(y)) + if (mag_is_inf(x) || mag_is_inf(y)) mag_inf(z); else mag_zero(z); diff --git a/mag/set_fmpr.c b/mag/set_fmpr.c index 7d3c67af..0564bbe7 100644 --- a/mag/set_fmpr.c +++ b/mag/set_fmpr.c @@ -24,6 +24,7 @@ ******************************************************************************/ #include "mag.h" +#include "arf.h" void mag_set_fmpr(mag_t x, const fmpr_t y) @@ -31,7 +32,7 @@ mag_set_fmpr(mag_t x, const fmpr_t y) arf_t t; arf_init(t); arf_set_fmpr(t, y); - mag_set_arf(x, t); + arf_get_mag(x, t); arf_clear(t); } diff --git a/mag/test/t-add_2exp_fmpz.c b/mag/test/t-add_2exp_fmpz.c new file mode 100644 index 00000000..b1f23d94 --- /dev/null +++ b/mag/test/t-add_2exp_fmpz.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" + +int main() +{ + long iter; + flint_rand_t state; + + printf("add_2exp_fmpz...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmpr_t x, y, z, z2, w; + mag_t xb, zb; + fmpz_t e; + + fmpr_init(x); + fmpr_init(y); + fmpr_init(z); + fmpr_init(z2); + fmpr_init(w); + + fmpz_init(e); + + mag_init(xb); + mag_init(zb); + + fmpr_randtest(x, state, 2 + n_randint(state, 200), 100); + fmpr_abs(x, x); + + fmpz_randtest(e, state, 100); + fmpr_one(y); + fmpr_mul_2exp_fmpz(y, y, e); + + fmpr_add(z, x, y, MAG_BITS + 10, FMPR_RND_DOWN); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + mag_set_fmpr(xb, x); + + mag_add_2exp_fmpz(zb, xb, e); + mag_get_fmpr(w, zb); + + MAG_CHECK_BITS(xb) + MAG_CHECK_BITS(zb) + + if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0)) + { + printf("FAIL\n\n"); + printf("x = "); fmpr_printd(x, 15); printf("\n\n"); + printf("y = "); fmpr_printd(y, 15); printf("\n\n"); + printf("z = "); fmpr_printd(z, 15); printf("\n\n"); + printf("w = "); fmpr_printd(w, 15); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(y); + fmpr_clear(z); + fmpr_clear(z2); + fmpr_clear(w); + + mag_clear(xb); + mag_clear(zb); + + fmpz_clear(e); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/mag/test/t-addmul.c b/mag/test/t-addmul.c new file mode 100644 index 00000000..2765c6fc --- /dev/null +++ b/mag/test/t-addmul.c @@ -0,0 +1,102 @@ +/*============================================================================= + + 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("addmul...."); + 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); + + fmpr_randtest(x, state, 2 + n_randint(state, 200), 100); + fmpr_randtest(y, state, 2 + n_randint(state, 200), 100); + fmpr_randtest(z, state, 2 + n_randint(state, 200), 100); + + fmpr_abs(x, x); + fmpr_abs(y, y); + fmpr_abs(z, z); + + mag_set_fmpr(xb, x); + mag_set_fmpr(yb, y); + mag_set_fmpr(zb, z); + + fmpr_addmul(z, x, y, MAG_BITS + 10, FMPR_RND_DOWN); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + mag_addmul(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_printd(x, 15); printf("\n\n"); + printf("y = "); fmpr_printd(y, 15); printf("\n\n"); + printf("z = "); fmpr_printd(z, 15); printf("\n\n"); + printf("w = "); fmpr_printd(w, 15); printf("\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; +} + diff --git a/mag/test/t-fast_add_2exp_si.c b/mag/test/t-fast_add_2exp_si.c new file mode 100644 index 00000000..d4beb660 --- /dev/null +++ b/mag/test/t-fast_add_2exp_si.c @@ -0,0 +1,95 @@ +/*============================================================================= + + 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("fast_add_2exp_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + fmpr_t x, y, z, z2, w; + mag_t xb, zb; + long e; + + fmpr_init(x); + fmpr_init(y); + fmpr_init(z); + fmpr_init(z2); + fmpr_init(w); + + mag_init(xb); + mag_init(zb); + + fmpr_randtest(x, state, 2 + n_randint(state, 200), 15); + fmpr_abs(x, x); + + e = n_randint(state, 10000) - n_randint(state, 10000); + fmpr_set_ui_2exp_si(y, 1, e); + + fmpr_add(z, x, y, FMPR_PREC_EXACT, FMPR_RND_DOWN); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + mag_set_fmpr(xb, x); + + mag_fast_add_2exp_si(zb, xb, e); + mag_get_fmpr(w, zb); + + MAG_CHECK_BITS(xb) + MAG_CHECK_BITS(zb) + + if (!(fmpr_cmpabs(z, w) <= 0 && fmpr_cmpabs(w, z2) <= 0)) + { + printf("FAIL\n\n"); + printf("x = "); fmpr_printd(x, 15); printf("\n\n"); + printf("y = "); fmpr_printd(y, 15); printf("\n\n"); + printf("z = "); fmpr_printd(z, 15); printf("\n\n"); + printf("w = "); fmpr_printd(w, 15); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(y); + fmpr_clear(z); + fmpr_clear(z2); + fmpr_clear(w); + + mag_clear(xb); + mag_clear(zb); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/mag/test/t-fast_addmul.c b/mag/test/t-fast_addmul.c new file mode 100644 index 00000000..4e373113 --- /dev/null +++ b/mag/test/t-fast_addmul.c @@ -0,0 +1,102 @@ +/*============================================================================= + + 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("fast_addmul...."); + 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); + + fmpr_randtest(x, state, 2 + n_randint(state, 200), 15); + fmpr_randtest(y, state, 2 + n_randint(state, 200), 15); + fmpr_randtest(z, state, 2 + n_randint(state, 200), 15); + + fmpr_abs(x, x); + fmpr_abs(y, y); + fmpr_abs(z, z); + + mag_set_fmpr(xb, x); + mag_set_fmpr(yb, y); + mag_set_fmpr(zb, z); + + fmpr_addmul(z, x, y, MAG_BITS + 10, FMPR_RND_DOWN); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + mag_fast_addmul(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_printd(x, 15); printf("\n\n"); + printf("y = "); fmpr_printd(y, 15); printf("\n\n"); + printf("z = "); fmpr_printd(z, 15); printf("\n\n"); + printf("w = "); fmpr_printd(w, 15); printf("\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; +} + diff --git a/mag/test/t-fast_mul.c b/mag/test/t-fast_mul.c new file mode 100644 index 00000000..f7fa457c --- /dev/null +++ b/mag/test/t-fast_mul.c @@ -0,0 +1,95 @@ +/*============================================================================= + + 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("fast_mul...."); + 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); + + fmpr_randtest(x, state, 2 + n_randint(state, 200), 15); + fmpr_randtest(y, state, 2 + n_randint(state, 200), 15); + + fmpr_mul(z, x, y, FMPR_PREC_EXACT, FMPR_RND_DOWN); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + mag_set_fmpr(xb, x); + mag_set_fmpr(yb, y); + mag_fast_mul(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_printd(x, 15); printf("\n\n"); + printf("y = "); fmpr_printd(y, 15); printf("\n\n"); + printf("z = "); fmpr_printd(z, 15); printf("\n\n"); + printf("w = "); fmpr_printd(w, 15); printf("\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; +} + diff --git a/mag/test/t-mul.c b/mag/test/t-mul.c new file mode 100644 index 00000000..7e75d567 --- /dev/null +++ b/mag/test/t-mul.c @@ -0,0 +1,95 @@ +/*============================================================================= + + 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("mul...."); + 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); + + fmpr_randtest(x, state, 2 + n_randint(state, 200), 100); + fmpr_randtest(y, state, 2 + n_randint(state, 200), 100); + + fmpr_mul(z, x, y, FMPR_PREC_EXACT, FMPR_RND_DOWN); + fmpr_mul_ui(z2, z, 1025, MAG_BITS, FMPR_RND_UP); + fmpr_mul_2exp_si(z2, z2, -10); + + mag_set_fmpr(xb, x); + mag_set_fmpr(yb, y); + mag_mul(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(); + } + + 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; +} +