diff --git a/arf.h b/arf.h index bbba6050..3a4468cd 100644 --- a/arf.h +++ b/arf.h @@ -777,6 +777,46 @@ _fmpz_add2_fast(fmpz_t z, const fmpz_t x, const fmpz_t y, long c) r3 += r2 < t1; \ } while (0) +#define ARF_MPN_MUL(_z, _x, _xn, _y, _yn) \ + if ((_xn) == (_yn)) \ + { \ + if ((_xn) == 1) \ + { \ + umul_ppmm((_z)[1], (_z)[0], (_x)[0], (_y)[0]); \ + } \ + else if ((_xn) == 2) \ + { \ + mp_limb_t __x1, __x0, __y1, __y0; \ + __x0 = (_x)[0]; \ + __x1 = (_x)[1]; \ + __y0 = (_y)[0]; \ + __y1 = (_y)[1]; \ + nn_mul_2x2((_z)[3], (_z)[2], (_z)[1], (_z)[0], __x1, __x0, __y1, __y0); \ + } \ + else if ((_x) == (_y)) \ + { \ + mpn_sqr((_z), (_x), (_xn)); \ + } \ + else \ + { \ + mpn_mul_n((_z), (_x), (_y), (_xn)); \ + } \ + } \ + else if ((_xn) > (_yn)) \ + { \ + if ((_yn) == 1) \ + (_z)[(_xn) + (_yn) - 1] = mpn_mul_1((_z), (_x), (_xn), (_y)[0]); \ + else \ + mpn_mul((_z), (_x), (_xn), (_y), (_yn)); \ + } \ + else \ + { \ + if ((_xn) == 1) \ + (_z)[(_xn) + (_yn) - 1] = mpn_mul_1((_z), (_y), (_yn), (_x)[0]); \ + else \ + mpn_mul((_z), (_y), (_yn), (_x), (_xn)); \ + } + #define ARF_MUL_STACK_ALLOC 40 #define ARF_MUL_TLS_ALLOC 1000 @@ -872,6 +912,8 @@ int _arf_add_mpn(arf_t z, mp_srcptr xp, mp_size_t xn, int xsgnbit, int arf_add(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd); +int arf_sub(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd); + #ifdef __cplusplus } #endif diff --git a/arf/sub.c b/arf/sub.c new file mode 100644 index 00000000..609dc9c3 --- /dev/null +++ b/arf/sub.c @@ -0,0 +1,88 @@ +/*============================================================================= + + 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_sub_special(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd) +{ + if (arf_is_zero(x)) + { + if (arf_is_zero(y)) + { + arf_zero(z); + return 0; + } + else + return arf_neg_round(z, y, prec, rnd); + } + else if (arf_is_zero(y)) + { + return arf_set_round(z, x, prec, rnd); + } + else if (arf_is_nan(x) || arf_is_nan(y) + || (arf_is_pos_inf(x) && arf_is_pos_inf(y)) + || (arf_is_neg_inf(x) && arf_is_neg_inf(y))) + { + arf_nan(z); + return 0; + } + else if (arf_is_special(x)) + { + arf_set(z, x); + return 0; + } + else + { + arf_neg(z, y); + return 0; + } +} + +int +arf_sub(arf_ptr z, arf_srcptr x, arf_srcptr y, long prec, arf_rnd_t rnd) +{ + mp_size_t xn, yn; + mp_srcptr xptr, yptr; + long shift; + + if (arf_is_special(x) || arf_is_special(y)) + { + return arf_sub_special(z, x, y, prec, rnd); + } + + shift = _fmpz_sub_small(ARF_EXPREF(x), ARF_EXPREF(y)); + + ARF_GET_MPN_READONLY(xptr, xn, x); + ARF_GET_MPN_READONLY(yptr, yn, y); + + if (shift >= 0) + return _arf_add_mpn(z, xptr, xn, ARF_SGNBIT(x), ARF_EXPREF(x), + yptr, yn, ARF_SGNBIT(y) ^ 1, shift, prec, rnd); + else + return _arf_add_mpn(z, yptr, yn, ARF_SGNBIT(y) ^ 1, ARF_EXPREF(y), + xptr, xn, ARF_SGNBIT(x), -shift, prec, rnd); +} + diff --git a/arf/test/t-sub.c b/arf/test/t-sub.c new file mode 100644 index 00000000..e42bf375 --- /dev/null +++ b/arf/test/t-sub.c @@ -0,0 +1,182 @@ +/*============================================================================= + + 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_sub_naive(arf_t z, const arf_t x, const arf_t y, long prec, arf_rnd_t rnd) +{ + fmpr_t a, b; + long r; + + fmpr_init(a); + fmpr_init(b); + + arf_get_fmpr(a, x); + arf_get_fmpr(b, y); + + r = fmpr_sub(a, a, b, prec, rnd); + + arf_set_fmpr(z, a); + + fmpr_clear(a); + fmpr_clear(b); + + return (r == FMPR_RESULT_EXACT) ? 0 : 1; +} + +int main() +{ + long iter, iter2; + flint_rand_t state; + + printf("sub...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; 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); + prec = 2 + n_randint(state, 2000); + + if (n_randint(state, 10) == 0 && + fmpz_bits(ARF_EXPREF(x)) < 10 && + fmpz_bits(ARF_EXPREF(y)) < 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_sub(z, x, y, prec, rnd); + r2 = arf_sub_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_sub(z, x, x, prec, rnd); + r2 = arf_sub_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_sub_naive(v, x, x, prec, rnd); + r1 = arf_sub(x, x, x, prec, rnd); + if (!arf_equal(v, x) || r1 != r2) + { + printf("FAIL (aliasing 2)!\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 3: + r2 = arf_sub_naive(v, x, y, prec, rnd); + r1 = arf_sub(x, x, y, prec, rnd); + if (!arf_equal(x, v) || r1 != r2) + { + printf("FAIL (aliasing 3)!\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("v = "); arf_print(v); printf("\n\n"); + printf("r1 = %ld, r2 = %ld\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_sub_naive(v, y, x, prec, rnd); + r1 = arf_sub(x, y, x, prec, rnd); + if (!arf_equal(x, 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("y = "); arf_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(y); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}