From 109afaab0f444b9059432ade48f48df1bef39173 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sat, 18 Feb 2017 15:56:00 +0100 Subject: [PATCH] add arf_sosq (sum of squares) --- arf.h | 2 + arf/sosq.c | 82 +++++++++++++++++++++++ arf/test/t-sosq.c | 162 +++++++++++++++++++++++++++++++++++++++++++++ doc/source/arf.rst | 5 ++ 4 files changed, 251 insertions(+) create mode 100644 arf/sosq.c create mode 100644 arf/test/t-sosq.c diff --git a/arf.h b/arf.h index 50ca49a9..f9202d26 100644 --- a/arf.h +++ b/arf.h @@ -1063,6 +1063,8 @@ arf_submul_fmpz(arf_ptr z, arf_srcptr x, const fmpz_t y, slong prec, arf_rnd_t r return arf_submul_mpz(z, x, COEFF_TO_PTR(*y), prec, rnd); } +int arf_sosq(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd); + int arf_div(arf_ptr z, arf_srcptr x, arf_srcptr y, slong prec, arf_rnd_t rnd); ARF_INLINE int diff --git a/arf/sosq.c b/arf/sosq.c new file mode 100644 index 00000000..21dbec1f --- /dev/null +++ b/arf/sosq.c @@ -0,0 +1,82 @@ +/* + Copyright (C) 2017 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arf.h" + +int +arf_sosq(arf_t res, const arf_t a, const arf_t b, slong prec, arf_rnd_t rnd) +{ + if (arf_is_special(a) || arf_is_special(b)) + { + if (arf_is_zero(a)) + return arf_mul(res, b, b, prec, rnd); + + if (arf_is_zero(b)) + return arf_mul(res, a, a, prec, rnd); + + if (arf_is_nan(a) || arf_is_nan(b)) + arf_nan(res); + else + arf_pos_inf(res); + + return 0; + } + else + { + mp_srcptr ap, bp; + int inexact; + mp_ptr tmp, aap, bbp; + mp_size_t an, bn, aan, bbn, alloc; + slong shift; + fmpz_t texp, uexp; + ARF_MUL_TMP_DECL + + ARF_GET_MPN_READONLY(ap, an, a); + ARF_GET_MPN_READONLY(bp, bn, b); + + fmpz_init(texp); + fmpz_init(uexp); + + _fmpz_add2_fast(texp, ARF_EXPREF(a), ARF_EXPREF(a), 0); + _fmpz_add2_fast(uexp, ARF_EXPREF(b), ARF_EXPREF(b), 0); + shift = _fmpz_sub_small(texp, uexp); + + aan = 2 * an; + bbn = 2 * bn; + alloc = aan + bbn; + + ARF_MUL_TMP_ALLOC(tmp, alloc) + aap = tmp; + bbp = tmp + aan; + + ARF_MPN_MUL(aap, ap, an, ap, an) + aan -= (aap[0] == 0); + aap += (aap[0] == 0); + + ARF_MPN_MUL(bbp, bp, bn, bp, bn) + bbn -= (bbp[0] == 0); + bbp += (bbp[0] == 0); + + if (shift >= 0) + inexact = _arf_add_mpn(res, aap, aan, 0, texp, + bbp, bbn, 0, shift, prec, rnd); + else + inexact = _arf_add_mpn(res, bbp, bbn, 0, uexp, + aap, aan, 0, -shift, prec, rnd); + + ARF_MUL_TMP_FREE(tmp, alloc) + fmpz_clear(texp); + fmpz_clear(uexp); + + return inexact; + } +} + diff --git a/arf/test/t-sosq.c b/arf/test/t-sosq.c new file mode 100644 index 00000000..8191382b --- /dev/null +++ b/arf/test/t-sosq.c @@ -0,0 +1,162 @@ +/* + Copyright (C) 2012 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arf.h" + +int +arf_sosq_naive(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd) +{ + arf_t t; + int inexact; + arf_init(t); + arf_mul(t, x, x, ARF_PREC_EXACT, ARF_RND_DOWN); + arf_mul(z, y, y, ARF_PREC_EXACT, ARF_RND_DOWN); + inexact = arf_add(z, z, t, prec, rnd); + arf_clear(t); + return inexact; +} + +int main() +{ + slong iter, iter2; + flint_rand_t state; + + flint_printf("sosq...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) + { + arf_t x, y, z, v; + slong 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_randtest_special(v, 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, 5)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + case 3: rnd = ARF_RND_CEIL; break; + default: rnd = ARF_RND_NEAR; break; + } + + switch (n_randint(state, 5)) + { + case 0: + r1 = arf_sosq(z, x, y, prec, rnd); + r2 = arf_sosq_naive(v, x, y, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + flint_printf("FAIL!\n"); + flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); + flint_printf("x = "); arf_print(x); flint_printf("\n\n"); + flint_printf("y = "); arf_print(y); flint_printf("\n\n"); + flint_printf("z = "); arf_print(z); flint_printf("\n\n"); + flint_printf("v = "); arf_print(v); flint_printf("\n\n"); + flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); + abort(); + } + break; + + case 1: + r1 = arf_sosq(z, x, x, prec, rnd); + r2 = arf_sosq_naive(v, x, x, prec, rnd); + if (!arf_equal(z, v) || r1 != r2) + { + flint_printf("FAIL (aliasing 1)!\n"); + flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); + flint_printf("x = "); arf_print(x); flint_printf("\n\n"); + flint_printf("z = "); arf_print(z); flint_printf("\n\n"); + flint_printf("v = "); arf_print(v); flint_printf("\n\n"); + flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); + abort(); + } + break; + + case 2: + r2 = arf_sosq_naive(v, x, y, prec, rnd); + r1 = arf_sosq(x, x, y, prec, rnd); + if (!arf_equal(v, x) || r1 != r2) + { + flint_printf("FAIL (aliasing 2)!\n"); + flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); + flint_printf("v = "); arf_print(v); flint_printf("\n\n"); + flint_printf("x = "); arf_print(x); flint_printf("\n\n"); + flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); + abort(); + } + break; + + case 3: + r2 = arf_sosq_naive(v, x, y, prec, rnd); + r1 = arf_sosq(y, x, y, prec, rnd); + if (!arf_equal(v, y) || r1 != r2) + { + flint_printf("FAIL (aliasing 3)!\n"); + flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); + flint_printf("x = "); arf_print(x); flint_printf("\n\n"); + flint_printf("y = "); arf_print(y); flint_printf("\n\n"); + flint_printf("v = "); arf_print(v); flint_printf("\n\n"); + flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); + abort(); + } + break; + + default: + r2 = arf_sosq_naive(v, x, x, prec, rnd); + r1 = arf_sosq(x, x, x, prec, rnd); + if (!arf_equal(v, x) || r1 != r2) + { + flint_printf("FAIL (aliasing 4)!\n"); + flint_printf("prec = %wd, rnd = %d\n\n", prec, rnd); + flint_printf("x = "); arf_print(x); flint_printf("\n\n"); + flint_printf("v = "); arf_print(v); flint_printf("\n\n"); + flint_printf("z = "); arf_print(z); flint_printf("\n\n"); + flint_printf("r1 = %wd, r2 = %wd\n", r1, r2); + abort(); + } + break; + } + } + + arf_clear(x); + arf_clear(y); + arf_clear(z); + arf_clear(v); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/doc/source/arf.rst b/doc/source/arf.rst index 9529ebee..3c0312b0 100644 --- a/doc/source/arf.rst +++ b/doc/source/arf.rst @@ -596,6 +596,11 @@ Addition and multiplication Sets `z = z - x \times y`, rounded to *prec* bits in the direction specified by *rnd*, returning nonzero iff the operation is inexact. +.. function:: int arf_sosq(arf_t z, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd) + + Sets `z = x^2 + y^2`, rounded to *prec* bits in the direction specified by *rnd*, + returning nonzero iff the operation is inexact. + Summation -------------------------------------------------------------------------------