From 4d66d4968a30b94070f5e667c03c18864c3b1dfa Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Mon, 5 May 2014 10:52:51 +0200 Subject: [PATCH] more multiplication and addmul/submul methods --- arf.h | 85 +++++++++++++++++++ arf/addmul.c | 164 +++++++++++++++++++++++++++++++++++ arf/mul_rnd_down.c | 55 ++++++++++++ arf/submul.c | 168 ++++++++++++++++++++++++++++++++++++ arf/test/t-addmul.c | 179 +++++++++++++++++++++++++++++++++++++++ arf/test/t-addmul_fmpz.c | 135 +++++++++++++++++++++++++++++ arf/test/t-addmul_si.c | 134 +++++++++++++++++++++++++++++ arf/test/t-addmul_ui.c | 133 +++++++++++++++++++++++++++++ arf/test/t-mul_fmpz.c | 124 +++++++++++++++++++++++++++ arf/test/t-mul_si.c | 123 +++++++++++++++++++++++++++ arf/test/t-mul_ui.c | 123 +++++++++++++++++++++++++++ arf/test/t-submul.c | 179 +++++++++++++++++++++++++++++++++++++++ arf/test/t-submul_fmpz.c | 135 +++++++++++++++++++++++++++++ arf/test/t-submul_si.c | 134 +++++++++++++++++++++++++++++ arf/test/t-submul_ui.c | 133 +++++++++++++++++++++++++++++ 15 files changed, 2004 insertions(+) create mode 100644 arf/addmul.c create mode 100644 arf/submul.c create mode 100644 arf/test/t-addmul.c create mode 100644 arf/test/t-addmul_fmpz.c create mode 100644 arf/test/t-addmul_si.c create mode 100644 arf/test/t-addmul_ui.c create mode 100644 arf/test/t-mul_fmpz.c create mode 100644 arf/test/t-mul_si.c create mode 100644 arf/test/t-mul_ui.c create mode 100644 arf/test/t-submul.c create mode 100644 arf/test/t-submul_fmpz.c create mode 100644 arf/test/t-submul_si.c create mode 100644 arf/test/t-submul_ui.c diff --git a/arf.h b/arf.h index 505aa06c..0054018f 100644 --- a/arf.h +++ b/arf.h @@ -915,6 +915,33 @@ int arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec); ? arf_mul_rnd_down(z, x, y, prec) \ : arf_mul_rnd_any(z, x, y, prec, rnd)) +static __inline__ int +arf_mul_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd) +{ + arf_t t; + arf_init_set_ui(t, y); /* no need to free */ + return arf_mul(z, x, t, prec, rnd); +} + +static __inline__ int +arf_mul_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd) +{ + arf_t t; + arf_init_set_si(t, y); /* no need to free */ + return arf_mul(z, x, t, prec, rnd); +} + +int arf_mul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd); + +static __inline__ int +arf_mul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd) +{ + if (!COEFF_IS_MPZ(*y)) + return arf_mul_si(z, x, *y, prec, rnd); + else + return arf_mul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd); +} + #define ARF_ADD_STACK_ALLOC 40 #define ARF_ADD_TLS_ALLOC 1000 @@ -967,6 +994,64 @@ int arf_sub_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd); int arf_sub_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd); int arf_sub_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd); +int arf_addmul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd); + +static __inline__ int +arf_addmul_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd) +{ + arf_t t; + arf_init_set_ui(t, y); /* no need to free */ + return arf_addmul(z, x, t, prec, rnd); +} + +static __inline__ int +arf_addmul_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd) +{ + arf_t t; + arf_init_set_si(t, y); /* no need to free */ + return arf_addmul(z, x, t, prec, rnd); +} + +int arf_addmul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd); + +static __inline__ int +arf_addmul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd) +{ + if (!COEFF_IS_MPZ(*y)) + return arf_addmul_si(z, x, *y, prec, rnd); + else + return arf_addmul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd); +} + +int arf_submul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd); + +static __inline__ int +arf_submul_ui(arf_ptr z, arf_srcptr x, ulong y, long prec, arf_rnd_t rnd) +{ + arf_t t; + arf_init_set_ui(t, y); /* no need to free */ + return arf_submul(z, x, t, prec, rnd); +} + +static __inline__ int +arf_submul_si(arf_ptr z, arf_srcptr x, long y, long prec, arf_rnd_t rnd) +{ + arf_t t; + arf_init_set_si(t, y); /* no need to free */ + return arf_submul(z, x, t, prec, rnd); +} + +int arf_submul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd); + +static __inline__ int +arf_submul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, long prec, arf_rnd_t rnd) +{ + if (!COEFF_IS_MPZ(*y)) + return arf_submul_si(z, x, *y, prec, rnd); + else + return arf_submul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd); +} + int arf_div(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd); static __inline__ int diff --git a/arf/addmul.c b/arf/addmul.c new file mode 100644 index 00000000..ab116b58 --- /dev/null +++ b/arf/addmul.c @@ -0,0 +1,164 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_addmul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd) +{ + mp_size_t xn, yn, zn, tn, alloc; + mp_srcptr xptr, yptr, zptr; + mp_ptr tptr, tptr2; + fmpz_t texp; + long shift; + int tsgnbit, inexact; + ARF_MUL_TMP_DECL + + if (arf_is_special(x) || arf_is_special(y) || arf_is_special(z)) + { + if (arf_is_zero(z)) + { + return arf_mul(z, x, y, prec, rnd); + } + else if (arf_is_finite(x) && arf_is_finite(y)) + { + return arf_set_round(z, z, prec, rnd); + } + else + { + /* todo: speed up */ + arf_t t; + arf_init(t); + arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + inexact = arf_add(z, z, t, prec, rnd); + arf_clear(t); + return inexact; + } + } + + tsgnbit = ARF_SGNBIT(x) ^ ARF_SGNBIT(y); + ARF_GET_MPN_READONLY(xptr, xn, x); + ARF_GET_MPN_READONLY(yptr, yn, y); + ARF_GET_MPN_READONLY(zptr, zn, z); + + fmpz_init(texp); + + _fmpz_add2_fast(texp, ARF_EXPREF(x), ARF_EXPREF(y), 0); + shift = _fmpz_sub_small(ARF_EXPREF(z), texp); + + alloc = tn = xn + yn; + ARF_MUL_TMP_ALLOC(tptr2, alloc) + tptr = tptr2; + + ARF_MPN_MUL(tptr, xptr, xn, yptr, yn); + + tn -= (tptr[0] == 0); + tptr += (tptr[0] == 0); + + if (shift >= 0) + inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z), + tptr, tn, tsgnbit, shift, prec, rnd); + else + inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp, + zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd); + + ARF_MUL_TMP_FREE(tptr2, alloc) + fmpz_clear(texp); + + return inexact; +} + +int +arf_addmul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd) +{ + mp_size_t xn, yn, zn, tn, alloc; + mp_srcptr xptr, yptr, zptr; + mp_ptr tptr, tptr2; + fmpz_t texp, yexp; + long shift; + int tsgnbit, ysgnbit, inexact; + ARF_MUL_TMP_DECL + + yn = FLINT_ABS(y->_mp_size); + + if (arf_is_special(x) || yn == 0 || arf_is_special(z)) + { + if (arf_is_zero(z)) + { + return arf_mul_mpz(z, x, y, prec, rnd); + } + else if (arf_is_finite(x)) + { + return arf_set_round(z, z, prec, rnd); + } + else + { + /* todo: speed up */ + arf_t t; + arf_init(t); + arf_mul_mpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + inexact = arf_add(z, z, t, prec, rnd); + arf_clear(t); + return inexact; + } + } + + ARF_GET_MPN_READONLY(xptr, xn, x); + + yptr = y->_mp_d; + ysgnbit = (y->_mp_size < 0); + *yexp = yn * FLINT_BITS; + + ARF_GET_MPN_READONLY(zptr, zn, z); + + fmpz_init(texp); + + tsgnbit = ARF_SGNBIT(x) ^ ysgnbit; + + alloc = tn = xn + yn; + ARF_MUL_TMP_ALLOC(tptr2, alloc) + tptr = tptr2; + + ARF_MPN_MUL(tptr, xptr, xn, yptr, yn); + + shift = (tptr[tn - 1] == 0) * FLINT_BITS; + tn -= (tptr[tn - 1] == 0); + + _fmpz_add2_fast(texp, ARF_EXPREF(x), yexp, -shift); + shift = _fmpz_sub_small(ARF_EXPREF(z), texp); + + if (shift >= 0) + inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z), + tptr, tn, tsgnbit, shift, prec, rnd); + else + inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp, + zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd); + + ARF_MUL_TMP_FREE(tptr2, alloc) + fmpz_clear(texp); + + return inexact; +} + diff --git a/arf/mul_rnd_down.c b/arf/mul_rnd_down.c index 0fe7b8b4..42aa1944 100644 --- a/arf/mul_rnd_down.c +++ b/arf/mul_rnd_down.c @@ -232,3 +232,58 @@ arf_mul_rnd_down(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec) } } +int +arf_mul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd) +{ + mp_size_t xn, yn; + long fix, shift; + int sgnbit, inexact; + + yn = FLINT_ABS(y->_mp_size); + xn = ARF_XSIZE(x); + xn >>= 1; + sgnbit = ARF_SGNBIT(x) ^ (y->_mp_size < 0); + + /* Either operand is a special value. */ + if (xn == 0 || yn == 0) + { + if (arf_is_finite(x)) + { + arf_zero(z); + } + else + { + arf_t t; + arf_init_set_si(t, mpz_sgn(y)); + arf_mul(z, x, t, prec, rnd); + arf_clear(t); + } + return 0; + } + else + { + mp_size_t zn, alloc; + mp_srcptr xptr, yptr; + mp_ptr tmp; + ARF_MUL_TMP_DECL + + ARF_GET_MPN_READONLY(xptr, xn, x); + yptr = y->_mp_d; + + alloc = zn = xn + yn; + ARF_MUL_TMP_ALLOC(tmp, alloc) + + ARF_MPN_MUL(tmp, xptr, xn, yptr, yn); + + shift = yn * FLINT_BITS - (tmp[zn - 1] == 0) * FLINT_BITS; + zn -= (tmp[zn - 1] == 0); + + inexact = _arf_set_round_mpn(z, &fix, tmp, zn, sgnbit, prec, rnd); + + _fmpz_add_fast(ARF_EXPREF(z), ARF_EXPREF(x), fix + shift); + ARF_MUL_TMP_FREE(tmp, alloc) + + return inexact; + } +} + diff --git a/arf/submul.c b/arf/submul.c new file mode 100644 index 00000000..b91fbf03 --- /dev/null +++ b/arf/submul.c @@ -0,0 +1,168 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_submul(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd) +{ + mp_size_t xn, yn, zn, tn, alloc; + mp_srcptr xptr, yptr, zptr; + mp_ptr tptr, tptr2; + fmpz_t texp; + long shift; + int tsgnbit, inexact; + ARF_MUL_TMP_DECL + + if (arf_is_special(x) || arf_is_special(y) || arf_is_special(z)) + { + if (arf_is_zero(z)) + { + /* TODO: make more efficient */ + arf_mul(z, x, y, ARF_PREC_EXACT, rnd); + return arf_neg_round(z, z, prec, rnd); + } + else if (arf_is_finite(x) && arf_is_finite(y)) + { + return arf_set_round(z, z, prec, rnd); + } + else + { + /* todo: speed up */ + arf_t t; + arf_init(t); + arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + inexact = arf_sub(z, z, t, prec, rnd); + arf_clear(t); + return inexact; + } + } + + tsgnbit = ARF_SGNBIT(x) ^ ARF_SGNBIT(y) ^ 1; + ARF_GET_MPN_READONLY(xptr, xn, x); + ARF_GET_MPN_READONLY(yptr, yn, y); + ARF_GET_MPN_READONLY(zptr, zn, z); + + fmpz_init(texp); + + _fmpz_add2_fast(texp, ARF_EXPREF(x), ARF_EXPREF(y), 0); + shift = _fmpz_sub_small(ARF_EXPREF(z), texp); + + alloc = tn = xn + yn; + ARF_MUL_TMP_ALLOC(tptr2, alloc) + tptr = tptr2; + + ARF_MPN_MUL(tptr, xptr, xn, yptr, yn); + + tn -= (tptr[0] == 0); + tptr += (tptr[0] == 0); + + if (shift >= 0) + inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z), + tptr, tn, tsgnbit, shift, prec, rnd); + else + inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp, + zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd); + + ARF_MUL_TMP_FREE(tptr2, alloc) + fmpz_clear(texp); + + return inexact; +} + +int +arf_submul_mpz(arf_ptr z, arf_srcptr x, const mpz_t y, long prec, arf_rnd_t rnd) +{ + mp_size_t xn, yn, zn, tn, alloc; + mp_srcptr xptr, yptr, zptr; + mp_ptr tptr, tptr2; + fmpz_t texp, yexp; + long shift; + int tsgnbit, ysgnbit, inexact; + ARF_MUL_TMP_DECL + + yn = FLINT_ABS(y->_mp_size); + + if (arf_is_special(x) || yn == 0 || arf_is_special(z)) + { + if (arf_is_zero(z)) + { + /* TODO: make more efficient */ + arf_mul_mpz(z, x, y, ARF_PREC_EXACT, rnd); + return arf_neg_round(z, z, prec, rnd); + } + else if (arf_is_finite(x)) + { + return arf_set_round(z, z, prec, rnd); + } + else + { + /* todo: speed up */ + arf_t t; + arf_init(t); + arf_mul_mpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + inexact = arf_sub(z, z, t, prec, rnd); + arf_clear(t); + return inexact; + } + } + + ARF_GET_MPN_READONLY(xptr, xn, x); + + yptr = y->_mp_d; + ysgnbit = (y->_mp_size > 0); + *yexp = yn * FLINT_BITS; + + ARF_GET_MPN_READONLY(zptr, zn, z); + + fmpz_init(texp); + + tsgnbit = ARF_SGNBIT(x) ^ ysgnbit; + + alloc = tn = xn + yn; + ARF_MUL_TMP_ALLOC(tptr2, alloc) + tptr = tptr2; + + ARF_MPN_MUL(tptr, xptr, xn, yptr, yn); + + shift = (tptr[tn - 1] == 0) * FLINT_BITS; + tn -= (tptr[tn - 1] == 0); + + _fmpz_add2_fast(texp, ARF_EXPREF(x), yexp, -shift); + shift = _fmpz_sub_small(ARF_EXPREF(z), texp); + + if (shift >= 0) + inexact = _arf_add_mpn(z, zptr, zn, ARF_SGNBIT(z), ARF_EXPREF(z), + tptr, tn, tsgnbit, shift, prec, rnd); + else + inexact = _arf_add_mpn(z, tptr, tn, tsgnbit, texp, + zptr, zn, ARF_SGNBIT(z), -shift, prec, rnd); + + ARF_MUL_TMP_FREE(tptr2, alloc) + fmpz_clear(texp); + + return inexact; +} + diff --git a/arf/test/t-addmul.c b/arf/test/t-addmul.c new file mode 100644 index 00000000..c76393dc --- /dev/null +++ b/arf/test/t-addmul.c @@ -0,0 +1,179 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_addmul_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_add(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("addmul...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, y, z, v; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(y); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + arf_randtest_special(y, state, 2000, 100); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(y)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 5)) + { + case 0: + r1 = arf_addmul(z, x, y, prec, rnd); + r2 = arf_addmul_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); arf_print(y); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + case 1: + r1 = arf_addmul(z, x, x, prec, rnd); + r2 = arf_addmul_naive(v, x, x, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL (aliasing 1)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + case 2: + r2 = arf_addmul_naive(v, v, v, prec, rnd); + r1 = arf_addmul(z, z, z, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing 2)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + case 3: + r2 = arf_addmul_naive(v, v, y, prec, rnd); + r1 = arf_addmul(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing 3)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = "); arf_print(y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_addmul_naive(v, x, v, prec, rnd); + r1 = arf_addmul(z, x, z, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL (aliasing 4)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(y); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-addmul_fmpz.c b/arf/test/t-addmul_fmpz.c new file mode 100644 index 00000000..0e9fa6cb --- /dev/null +++ b/arf/test/t-addmul_fmpz.c @@ -0,0 +1,135 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_addmul_fmpz_naive(arf_t z, const arf_t x, const fmpz_t y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul_fmpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_add(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("addmul_fmpz...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, z, v; + fmpz_t y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + fmpz_init(y); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + fmpz_randtest(y, state, 2000); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_addmul_fmpz(z, x, y, prec, rnd); + r2 = arf_addmul_fmpz_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); fmpz_print(y); printf("\n\n"); + printf("z = "); arf_debug(z); printf("\n\n"); + printf("v = "); arf_debug(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_addmul_fmpz_naive(v, v, y, prec, rnd); + r1 = arf_addmul_fmpz(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = "); fmpz_print(y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + fmpz_clear(y); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-addmul_si.c b/arf/test/t-addmul_si.c new file mode 100644 index 00000000..efed5570 --- /dev/null +++ b/arf/test/t-addmul_si.c @@ -0,0 +1,134 @@ +/*============================================================================= + + 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 "arf.h" +#include "long_extras.h" + +int +arf_addmul_si_naive(arf_t z, const arf_t x, long y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul_si(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_add(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("addmul_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, z, v; + long y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + y = z_randtest(state); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_addmul_si(z, x, y, prec, rnd); + r2 = arf_addmul_si_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = %ld", y); printf("\n\n"); + printf("z = "); arf_debug(z); printf("\n\n"); + printf("v = "); arf_debug(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_addmul_si_naive(v, v, y, prec, rnd); + r1 = arf_addmul_si(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = %ld", y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-addmul_ui.c b/arf/test/t-addmul_ui.c new file mode 100644 index 00000000..d8f1c8ae --- /dev/null +++ b/arf/test/t-addmul_ui.c @@ -0,0 +1,133 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_addmul_ui_naive(arf_t z, const arf_t x, ulong y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul_ui(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_add(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("addmul_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, z, v; + ulong y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + y = n_randtest(state); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_addmul_ui(z, x, y, prec, rnd); + r2 = arf_addmul_ui_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = %lu", y); printf("\n\n"); + printf("z = "); arf_debug(z); printf("\n\n"); + printf("v = "); arf_debug(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_addmul_ui_naive(v, v, y, prec, rnd); + r1 = arf_addmul_ui(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = %lu", y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-mul_fmpz.c b/arf/test/t-mul_fmpz.c new file mode 100644 index 00000000..ffd5fed8 --- /dev/null +++ b/arf/test/t-mul_fmpz.c @@ -0,0 +1,124 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_mul_fmpz_naive(arf_t z, const arf_t x, const fmpz_t y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int r; + arf_init(t); + arf_set_fmpz(t, y); + r = arf_mul(z, x, t, prec, rnd); + arf_clear(t); + return r; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("mul_fmpz...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arf_t x, z, v; + fmpz_t y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + fmpz_init(y); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 10); + fmpz_randtest(y, state, 2000); + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10)) + prec = ARF_PREC_EXACT; + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_mul_fmpz(z, x, y, prec, rnd); + r2 = arf_mul_fmpz_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); fmpz_print(y); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_mul_fmpz_naive(v, x, y, prec, rnd); + r1 = arf_mul_fmpz(x, x, y, prec, rnd); + if (!arf_equal(x, v) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); fmpz_print(y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + fmpz_clear(y); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-mul_si.c b/arf/test/t-mul_si.c new file mode 100644 index 00000000..0fb5b732 --- /dev/null +++ b/arf/test/t-mul_si.c @@ -0,0 +1,123 @@ +/*============================================================================= + + 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 "arf.h" +#include "long_extras.h" + +int +arf_mul_si_naive(arf_t z, const arf_t x, long y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int r; + arf_init(t); + arf_set_si(t, y); + r = arf_mul(z, x, t, prec, rnd); + arf_clear(t); + return r; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("mul_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arf_t x, z, v; + long y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 10); + y = z_randtest(state); + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0) + prec = ARF_PREC_EXACT; + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_mul_si(z, x, y, prec, rnd); + r2 = arf_mul_si_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); printf("%ld", y); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_mul_si_naive(v, x, y, prec, rnd); + r1 = arf_mul_si(x, x, y, prec, rnd); + if (!arf_equal(x, v) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); printf("%ld", y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-mul_ui.c b/arf/test/t-mul_ui.c new file mode 100644 index 00000000..302dad3f --- /dev/null +++ b/arf/test/t-mul_ui.c @@ -0,0 +1,123 @@ +/*============================================================================= + + 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 "arf.h" +#include "ulong_extras.h" + +int +arf_mul_ui_naive(arf_t z, const arf_t x, ulong y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int r; + arf_init(t); + arf_set_ui(t, y); + r = arf_mul(z, x, t, prec, rnd); + arf_clear(t); + return r; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("mul_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arf_t x, z, v; + ulong y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 10); + y = n_randtest(state); + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0) + prec = ARF_PREC_EXACT; + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_mul_ui(z, x, y, prec, rnd); + r2 = arf_mul_ui_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); printf("%ld", y); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_mul_ui_naive(v, x, y, prec, rnd); + r1 = arf_mul_ui(x, x, y, prec, rnd); + if (!arf_equal(x, v) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); printf("%ld", y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-submul.c b/arf/test/t-submul.c new file mode 100644 index 00000000..3eb605bc --- /dev/null +++ b/arf/test/t-submul.c @@ -0,0 +1,179 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_submul_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_sub(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("submul...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, y, z, v; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(y); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + arf_randtest_special(y, state, 2000, 100); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(y)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 5)) + { + case 0: + r1 = arf_submul(z, x, y, prec, rnd); + r2 = arf_submul_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); arf_print(y); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + case 1: + r1 = arf_submul(z, x, x, prec, rnd); + r2 = arf_submul_naive(v, x, x, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL (aliasing 1)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + case 2: + r2 = arf_submul_naive(v, v, v, prec, rnd); + r1 = arf_submul(z, z, z, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing 2)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + case 3: + r2 = arf_submul_naive(v, v, y, prec, rnd); + r1 = arf_submul(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing 3)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = "); arf_print(y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_submul_naive(v, x, v, prec, rnd); + r1 = arf_submul(z, x, z, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL (aliasing 4)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(y); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-submul_fmpz.c b/arf/test/t-submul_fmpz.c new file mode 100644 index 00000000..fc3c43f7 --- /dev/null +++ b/arf/test/t-submul_fmpz.c @@ -0,0 +1,135 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_submul_fmpz_naive(arf_t z, const arf_t x, const fmpz_t y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul_fmpz(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_sub(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("submul_fmpz...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, z, v; + fmpz_t y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + fmpz_init(y); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + fmpz_randtest(y, state, 2000); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_submul_fmpz(z, x, y, prec, rnd); + r2 = arf_submul_fmpz_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); fmpz_print(y); printf("\n\n"); + printf("z = "); arf_debug(z); printf("\n\n"); + printf("v = "); arf_debug(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_submul_fmpz_naive(v, v, y, prec, rnd); + r1 = arf_submul_fmpz(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = "); fmpz_print(y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + fmpz_clear(y); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-submul_si.c b/arf/test/t-submul_si.c new file mode 100644 index 00000000..c6e6855b --- /dev/null +++ b/arf/test/t-submul_si.c @@ -0,0 +1,134 @@ +/*============================================================================= + + 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 "arf.h" +#include "long_extras.h" + +int +arf_submul_si_naive(arf_t z, const arf_t x, long y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul_si(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_sub(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("submul_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, z, v; + long y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + y = z_randtest(state); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_submul_si(z, x, y, prec, rnd); + r2 = arf_submul_si_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = %ld", y); printf("\n\n"); + printf("z = "); arf_debug(z); printf("\n\n"); + printf("v = "); arf_debug(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_submul_si_naive(v, v, y, prec, rnd); + r1 = arf_submul_si(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = %ld", y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-submul_ui.c b/arf/test/t-submul_ui.c new file mode 100644 index 00000000..a2afc196 --- /dev/null +++ b/arf/test/t-submul_ui.c @@ -0,0 +1,133 @@ +/*============================================================================= + + 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 "arf.h" + +int +arf_submul_ui_naive(arf_t z, const arf_t x, ulong y, long prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + + arf_init(t); + arf_mul_ui(t, x, y, ARF_PREC_EXACT, ARF_RND_DOWN); + + inexact = arf_sub(z, z, t, prec, rnd); + + arf_clear(t); + + return inexact; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("submul_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000; iter++) + { + arf_t x, z, v; + ulong y; + long prec, r1, r2; + arf_rnd_t rnd; + + arf_init(x); + arf_init(z); + arf_init(v); + + for (iter2 = 0; iter2 < 100; iter2++) + { + arf_randtest_special(x, state, 2000, 100); + y = n_randtest(state); + arf_randtest_special(z, state, 2000, 100); + arf_set(v, z); + + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(z)) < 10) + { + prec = ARF_PREC_EXACT; + } + + switch (n_randint(state, 4)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + default: rnd = ARF_RND_CEIL; break; + } + + switch (n_randint(state, 2)) + { + case 0: + r1 = arf_submul_ui(z, x, y, prec, rnd); + r2 = arf_submul_ui_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + printf("FAIL!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = %lu", y); printf("\n\n"); + printf("z = "); arf_debug(z); printf("\n\n"); + printf("v = "); arf_debug(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_submul_ui_naive(v, v, y, prec, rnd); + r1 = arf_submul_ui(z, z, y, prec, rnd); + if (!arf_equal(v, z) || r1 != r2) + { + printf("FAIL (aliasing)!\n"); + printf("prec = %ld, rnd = %d\n\n", prec, rnd); + printf("y = %lu", y); printf("\n\n"); + printf("v = "); arf_print(v); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}