diff --git a/acb.h b/acb.h index c0a3068d..65477298 100644 --- a/acb.h +++ b/acb.h @@ -602,6 +602,17 @@ void acb_dot(acb_t res, const acb_t initial, int subtract, void acb_approx_dot(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec); +void acb_dot_ui(acb_t res, const acb_t initial, int subtract, + acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec); +void acb_dot_si(acb_t res, const acb_t initial, int subtract, + acb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec); +void acb_dot_uiui(acb_t res, const acb_t initial, int subtract, + acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec); +void acb_dot_siui(acb_t res, const acb_t initial, int subtract, + acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec); +void acb_dot_fmpz(acb_t res, const acb_t initial, int subtract, + acb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec); + void acb_inv(acb_t z, const acb_t x, slong prec); void acb_div(acb_t z, const acb_t x, const acb_t y, slong prec); diff --git a/acb/dot_fmpz.c b/acb/dot_fmpz.c new file mode 100644 index 00000000..dd52cdfb --- /dev/null +++ b/acb/dot_fmpz.c @@ -0,0 +1,148 @@ +/* + Copyright (C) 2021 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 "acb.h" + +void +acb_dot_fmpz(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i, ssize, size, tmp_size; + mp_ptr ztmp; + fmpz v; + ulong av, al; + unsigned int bc; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + acb_zero(res); + else + { + acb_mul_fmpz(res, x, y, prec); + if (subtract) + acb_neg(res, res); + } + return; + } + else if (len <= 0) + { + acb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + tmp_size = 0; + for (i = 0; i < len; i++) + { + v = y[i * ystep]; + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + + if (v == 0) + { + ARF_XSIZE(arb_midref(t + i)) = 0; + ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO; + } + else if (!COEFF_IS_MPZ(v)) + { + av = FLINT_ABS(v); + count_leading_zeros(bc, av); + + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0); + } + else + { + __mpz_struct * z = COEFF_TO_PTR(v); + + ssize = z->_mp_size; + size = FLINT_ABS(ssize); + + av = z->_mp_d[size - 1]; + count_leading_zeros(bc, av); + + if (size == 1) + { + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, ssize < 0); + } + else if (size == 2) + { + al = z->_mp_d[0]; + + ARF_EXP(arb_midref(t + i)) = 2 * FLINT_BITS - bc; + + if (bc != 0) + { + av = (av << bc) | (al >> (FLINT_BITS - bc)); + al = al << bc; + } + + ARF_NOPTR_D(arb_midref(t + i))[0] = al; + ARF_NOPTR_D(arb_midref(t + i))[1] = av; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(2, ssize < 0); + } + else + { + if (bc != 0) + { + tmp_size += size; + /* use to flag tmp where we need tmp storage */ + MAG_MAN(arb_radref(t + i)) = bc; + } + + ARF_EXP(arb_midref(t + i)) = size * FLINT_BITS - bc; + ARF_PTR_D(arb_midref(t + i)) = z->_mp_d; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(size, ssize < 0); + } + } + } + + if (tmp_size != 0) + { + ztmp = TMP_ALLOC(sizeof(mp_limb_t) * tmp_size); + + for (i = 0; i < len; i++) + { + bc = MAG_MAN(arb_radref(t + i)); + + if (bc != 0) + { + size = ARF_SIZE(arb_midref(t + i)); + + mpn_lshift(ztmp, ARF_PTR_D(arb_midref(t + i)), size, bc); + ARF_PTR_D(arb_midref(t + i)) = ztmp; + ztmp += size; + } + + MAG_MAN(arb_radref(t + i)) = 0; + } + } + + TMP_END; + + arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec); + arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/acb/dot_si.c b/acb/dot_si.c new file mode 100644 index 00000000..ee24894b --- /dev/null +++ b/acb/dot_si.c @@ -0,0 +1,77 @@ +/* + Copyright (C) 2021 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 "acb.h" + +void +acb_dot_si(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + slong v; + ulong av; + unsigned int bc; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + acb_zero(res); + else + { + acb_mul_si(res, x, y[0], prec); + if (subtract) + acb_neg(res, res); + } + return; + } + else if (len <= 0) + { + acb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + v = y[i * ystep]; + + if (v == 0) + { + ARF_XSIZE(arb_midref(t + i)) = 0; + ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO; + } + else + { + av = FLINT_ABS(v); + count_leading_zeros(bc, av); + + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0); + } + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec); + arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/acb/dot_siui.c b/acb/dot_siui.c new file mode 100644 index 00000000..eb243525 --- /dev/null +++ b/acb/dot_siui.c @@ -0,0 +1,115 @@ +/* + Copyright (C) 2021 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 "acb.h" + +static void +arf_shallow_set_siui(arf_t res, ulong vhi, ulong vlo) +{ + int negative; + unsigned int bc; + + negative = ((slong) vhi) < 0; + + if (negative) + { + vhi = -vhi - (vlo != 0); + vlo = -vlo; + } + + if (vhi == 0) + { + if (vlo == 0) + { + ARF_XSIZE(res) = 0; + ARF_EXP(res) = ARF_EXP_ZERO; + } + else + { + count_leading_zeros(bc, vlo); + ARF_EXP(res) = FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative); + } + } + else if (vlo == 0) + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vhi << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative); + } + else + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + if (bc == 0) + ARF_NOPTR_D(res)[1] = vhi; + else + ARF_NOPTR_D(res)[1] = (vhi << bc) | (vlo >> (FLINT_BITS - bc)); + ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, negative); + } +} + +void +acb_dot_siui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + ulong vhi, vlo; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + acb_zero(res); + else + { + arf_t t; + arf_shallow_set_siui(t, y[1], y[0]); + arb_mul_arf(acb_realref(res), acb_realref(x), t, prec); + arb_mul_arf(acb_imagref(res), acb_imagref(x), t, prec); + if (subtract) + acb_neg(res, res); + } + return; + } + else if (len <= 0) + { + acb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + vlo = y[2 * i * ystep]; + vhi = y[2 * i * ystep + 1]; + + arf_shallow_set_siui(arb_midref(t + i), vhi, vlo); + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec); + arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/acb/dot_ui.c b/acb/dot_ui.c new file mode 100644 index 00000000..54d157d9 --- /dev/null +++ b/acb/dot_ui.c @@ -0,0 +1,75 @@ +/* + Copyright (C) 2021 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 "acb.h" + +void +acb_dot_ui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + ulong v; + unsigned int bc; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + acb_zero(res); + else + { + acb_mul_ui(res, x, y[0], prec); + if (subtract) + acb_neg(res, res); + } + return; + } + else if (len <= 0) + { + acb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + v = y[i * ystep]; + + if (v == 0) + { + ARF_XSIZE(arb_midref(t + i)) = 0; + ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO; + } + else + { + count_leading_zeros(bc, v); + + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = v << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, 0); + } + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec); + arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/acb/dot_uiui.c b/acb/dot_uiui.c new file mode 100644 index 00000000..c403ce70 --- /dev/null +++ b/acb/dot_uiui.c @@ -0,0 +1,106 @@ +/* + Copyright (C) 2021 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 "acb.h" + +static void +arf_shallow_set_uiui(arf_t res, ulong vhi, ulong vlo) +{ + unsigned int bc; + + if (vhi == 0) + { + if (vlo == 0) + { + ARF_XSIZE(res) = 0; + ARF_EXP(res) = ARF_EXP_ZERO; + } + else + { + count_leading_zeros(bc, vlo); + ARF_EXP(res) = FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0); + } + } + else if (vlo == 0) + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vhi << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0); + } + else + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + if (bc == 0) + ARF_NOPTR_D(res)[1] = vhi; + else + ARF_NOPTR_D(res)[1] = (vhi << bc) | (vlo >> (FLINT_BITS - bc)); + ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, 0); + } +} + +void +acb_dot_uiui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + ulong vhi, vlo; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + acb_zero(res); + else + { + arf_t t; + arf_shallow_set_uiui(t, y[1], y[0]); + arb_mul_arf(acb_realref(res), acb_realref(x), t, prec); + arb_mul_arf(acb_imagref(res), acb_imagref(x), t, prec); + if (subtract) + acb_neg(res, res); + } + return; + } + else if (len <= 0) + { + acb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + vlo = y[2 * i * ystep]; + vhi = y[2 * i * ystep + 1]; + + arf_shallow_set_uiui(arb_midref(t + i), vhi, vlo); + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(((arb_ptr) res) + 0, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 0, subtract, ((arb_srcptr) x) + 0, 2 * xstep, t, 1, len, prec); + arb_dot(((arb_ptr) res) + 1, (initial == NULL) ? NULL : ((arb_srcptr) initial) + 1, subtract, ((arb_srcptr) x) + 1, 2 * xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/acb/test/t-dot_fmpz.c b/acb/test/t-dot_fmpz.c new file mode 100644 index 00000000..3ec11248 --- /dev/null +++ b/acb/test/t-dot_fmpz.c @@ -0,0 +1,102 @@ +/* + Copyright (C) 2021 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 "acb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_siui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + acb_ptr x, y; + fmpz * w; + acb_t s1, s2, z; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _acb_vec_init(len); + y = _acb_vec_init(len); + w = _fmpz_vec_init(len); + acb_init(s1); + acb_init(s2); + acb_init(z); + + for (i = 0; i < len; i++) + { + acb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + fmpz_randtest(w + i, state, 1 + n_randint(state, 200)); + acb_set_fmpz(y + i, w + i); + } + + acb_randtest(s1, state, 200, 10); + acb_randtest(s2, state, 200, 10); + acb_randtest(z, state, 200, 10); + + acb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + acb_dot_fmpz(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + len - 1) : w, revy ? -1 : 1, + len, prec); + + if (!acb_overlaps(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i)); + flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(s1); + acb_clear(s2); + acb_clear(z); + _acb_vec_clear(x, len); + _acb_vec_clear(y, len); + _fmpz_vec_clear(w, len); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb/test/t-dot_si.c b/acb/test/t-dot_si.c new file mode 100644 index 00000000..7056d2e9 --- /dev/null +++ b/acb/test/t-dot_si.c @@ -0,0 +1,102 @@ +/* + Copyright (C) 2021 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 "acb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + acb_ptr x, y; + slong * w; + acb_t s1, s2, z; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _acb_vec_init(len); + y = _acb_vec_init(len); + w = flint_malloc(sizeof(ulong) * len); + acb_init(s1); + acb_init(s2); + acb_init(z); + + for (i = 0; i < len; i++) + { + acb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[i] = n_randtest(state); + acb_set_si(y + i, w[i]); + } + + acb_randtest(s1, state, 200, 10); + acb_randtest(s2, state, 200, 10); + acb_randtest(z, state, 200, 10); + + acb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + acb_dot_si(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + len - 1) : w, revy ? -1 : 1, + len, prec); + + if (!acb_equal(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i)); + flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(s1); + acb_clear(s2); + acb_clear(z); + _acb_vec_clear(x, len); + _acb_vec_clear(y, len); + flint_free(w); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb/test/t-dot_siui.c b/acb/test/t-dot_siui.c new file mode 100644 index 00000000..41d8f958 --- /dev/null +++ b/acb/test/t-dot_siui.c @@ -0,0 +1,107 @@ +/* + Copyright (C) 2021 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 "acb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_siui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + acb_ptr x, y; + ulong * w; + acb_t s1, s2, z; + fmpz_t c; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _acb_vec_init(len); + y = _acb_vec_init(len); + w = flint_malloc(2 * sizeof(ulong) * len); + acb_init(s1); + acb_init(s2); + acb_init(z); + fmpz_init(c); + + for (i = 0; i < len; i++) + { + acb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[2 * i] = n_randtest(state); + w[2 * i + 1] = n_randtest(state); + fmpz_set_signed_uiui(c, w[2 * i + 1], w[2 * i]); + acb_set_fmpz(y + i, c); + } + + acb_randtest(s1, state, 200, 10); + acb_randtest(s2, state, 200, 10); + acb_randtest(z, state, 200, 10); + + acb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + acb_dot_siui(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + 2 * len - 2) : w, revy ? -1 : 1, + len, prec); + + if (!acb_overlaps(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i)); + flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(s1); + acb_clear(s2); + acb_clear(z); + _acb_vec_clear(x, len); + _acb_vec_clear(y, len); + flint_free(w); + fmpz_clear(c); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb/test/t-dot_ui.c b/acb/test/t-dot_ui.c new file mode 100644 index 00000000..37d0d411 --- /dev/null +++ b/acb/test/t-dot_ui.c @@ -0,0 +1,102 @@ +/* + Copyright (C) 2021 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 "acb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + acb_ptr x, y; + ulong * w; + acb_t s1, s2, z; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _acb_vec_init(len); + y = _acb_vec_init(len); + w = flint_malloc(sizeof(ulong) * len); + acb_init(s1); + acb_init(s2); + acb_init(z); + + for (i = 0; i < len; i++) + { + acb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[i] = n_randtest(state); + acb_set_ui(y + i, w[i]); + } + + acb_randtest(s1, state, 200, 10); + acb_randtest(s2, state, 200, 10); + acb_randtest(z, state, 200, 10); + + acb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + acb_dot_ui(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + len - 1) : w, revy ? -1 : 1, + len, prec); + + if (!acb_equal(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i)); + flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(s1); + acb_clear(s2); + acb_clear(z); + _acb_vec_clear(x, len); + _acb_vec_clear(y, len); + flint_free(w); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/acb/test/t-dot_uiui.c b/acb/test/t-dot_uiui.c new file mode 100644 index 00000000..6b48ee4c --- /dev/null +++ b/acb/test/t-dot_uiui.c @@ -0,0 +1,107 @@ +/* + Copyright (C) 2021 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 "acb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_uiui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + acb_ptr x, y; + ulong * w; + acb_t s1, s2, z; + fmpz_t c; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _acb_vec_init(len); + y = _acb_vec_init(len); + w = flint_malloc(2 * sizeof(ulong) * len); + acb_init(s1); + acb_init(s2); + acb_init(z); + fmpz_init(c); + + for (i = 0; i < len; i++) + { + acb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[2 * i] = n_randtest(state); + w[2 * i + 1] = n_randtest(state); + fmpz_set_uiui(c, w[2 * i + 1], w[2 * i]); + acb_set_fmpz(y + i, c); + } + + acb_randtest(s1, state, 200, 10); + acb_randtest(s2, state, 200, 10); + acb_randtest(z, state, 200, 10); + + acb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + acb_dot_uiui(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + 2 * len - 2) : w, revy ? -1 : 1, + len, prec); + + if (!acb_overlaps(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); acb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", acb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); acb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(x + i)); + flint_printf("y[%wd] = ", i); acb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", acb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); acb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); acb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + acb_clear(s1); + acb_clear(s2); + acb_clear(z); + _acb_vec_clear(x, len); + _acb_vec_clear(y, len); + flint_free(w); + fmpz_clear(c); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb.h b/arb.h index 47969373..10d93c4a 100644 --- a/arb.h +++ b/arb.h @@ -435,6 +435,17 @@ void arb_dot(arb_t res, const arb_t initial, int subtract, void arb_approx_dot(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec); +void arb_dot_ui(arb_t res, const arb_t initial, int subtract, + arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec); +void arb_dot_si(arb_t res, const arb_t initial, int subtract, + arb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec); +void arb_dot_uiui(arb_t res, const arb_t initial, int subtract, + arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec); +void arb_dot_siui(arb_t res, const arb_t initial, int subtract, + arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec); +void arb_dot_fmpz(arb_t res, const arb_t initial, int subtract, + arb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec); + void arb_div(arb_t z, const arb_t x, const arb_t y, slong prec); void arb_div_arf(arb_t z, const arb_t x, const arf_t y, slong prec); void arb_div_si(arb_t z, const arb_t x, slong y, slong prec); diff --git a/arb/dot_fmpz.c b/arb/dot_fmpz.c new file mode 100644 index 00000000..56b1a5ff --- /dev/null +++ b/arb/dot_fmpz.c @@ -0,0 +1,147 @@ +/* + Copyright (C) 2021 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 "arb.h" + +void +arb_dot_fmpz(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i, ssize, size, tmp_size; + mp_ptr ztmp; + fmpz v; + ulong av, al; + unsigned int bc; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + arb_zero(res); + else + { + arb_mul_fmpz(res, x, y, prec); + if (subtract) + arb_neg(res, res); + } + return; + } + else if (len <= 0) + { + arb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + tmp_size = 0; + for (i = 0; i < len; i++) + { + v = y[i * ystep]; + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + + if (v == 0) + { + ARF_XSIZE(arb_midref(t + i)) = 0; + ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO; + } + else if (!COEFF_IS_MPZ(v)) + { + av = FLINT_ABS(v); + count_leading_zeros(bc, av); + + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0); + } + else + { + __mpz_struct * z = COEFF_TO_PTR(v); + + ssize = z->_mp_size; + size = FLINT_ABS(ssize); + + av = z->_mp_d[size - 1]; + count_leading_zeros(bc, av); + + if (size == 1) + { + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, ssize < 0); + } + else if (size == 2) + { + al = z->_mp_d[0]; + + ARF_EXP(arb_midref(t + i)) = 2 * FLINT_BITS - bc; + + if (bc != 0) + { + av = (av << bc) | (al >> (FLINT_BITS - bc)); + al = al << bc; + } + + ARF_NOPTR_D(arb_midref(t + i))[0] = al; + ARF_NOPTR_D(arb_midref(t + i))[1] = av; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(2, ssize < 0); + } + else + { + if (bc != 0) + { + tmp_size += size; + /* use to flag tmp where we need tmp storage */ + MAG_MAN(arb_radref(t + i)) = bc; + } + + ARF_EXP(arb_midref(t + i)) = size * FLINT_BITS - bc; + ARF_PTR_D(arb_midref(t + i)) = z->_mp_d; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(size, ssize < 0); + } + } + } + + if (tmp_size != 0) + { + ztmp = TMP_ALLOC(sizeof(mp_limb_t) * tmp_size); + + for (i = 0; i < len; i++) + { + bc = MAG_MAN(arb_radref(t + i)); + + if (bc != 0) + { + size = ARF_SIZE(arb_midref(t + i)); + + mpn_lshift(ztmp, ARF_PTR_D(arb_midref(t + i)), size, bc); + ARF_PTR_D(arb_midref(t + i)) = ztmp; + ztmp += size; + } + + MAG_MAN(arb_radref(t + i)) = 0; + } + } + + TMP_END; + + arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/arb/dot_si.c b/arb/dot_si.c new file mode 100644 index 00000000..286419ec --- /dev/null +++ b/arb/dot_si.c @@ -0,0 +1,76 @@ +/* + Copyright (C) 2021 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 "arb.h" + +void +arb_dot_si(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + slong v; + ulong av; + unsigned int bc; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + arb_zero(res); + else + { + arb_mul_si(res, x, y[0], prec); + if (subtract) + arb_neg(res, res); + } + return; + } + else if (len <= 0) + { + arb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + v = y[i * ystep]; + + if (v == 0) + { + ARF_XSIZE(arb_midref(t + i)) = 0; + ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO; + } + else + { + av = FLINT_ABS(v); + count_leading_zeros(bc, av); + + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = av << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, v < 0); + } + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/arb/dot_siui.c b/arb/dot_siui.c new file mode 100644 index 00000000..082b9646 --- /dev/null +++ b/arb/dot_siui.c @@ -0,0 +1,113 @@ +/* + Copyright (C) 2021 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 "arb.h" + +static void +arf_shallow_set_siui(arf_t res, ulong vhi, ulong vlo) +{ + int negative; + unsigned int bc; + + negative = ((slong) vhi) < 0; + + if (negative) + { + vhi = -vhi - (vlo != 0); + vlo = -vlo; + } + + if (vhi == 0) + { + if (vlo == 0) + { + ARF_XSIZE(res) = 0; + ARF_EXP(res) = ARF_EXP_ZERO; + } + else + { + count_leading_zeros(bc, vlo); + ARF_EXP(res) = FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative); + } + } + else if (vlo == 0) + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vhi << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, negative); + } + else + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + if (bc == 0) + ARF_NOPTR_D(res)[1] = vhi; + else + ARF_NOPTR_D(res)[1] = (vhi << bc) | (vlo >> (FLINT_BITS - bc)); + ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, negative); + } +} + +void +arb_dot_siui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + ulong vhi, vlo; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + arb_zero(res); + else + { + arf_t t; + arf_shallow_set_siui(t, y[1], y[0]); + arb_mul_arf(res, x, t, prec); + if (subtract) + arb_neg(res, res); + } + return; + } + else if (len <= 0) + { + arb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + vlo = y[2 * i * ystep]; + vhi = y[2 * i * ystep + 1]; + + arf_shallow_set_siui(arb_midref(t + i), vhi, vlo); + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/arb/dot_ui.c b/arb/dot_ui.c new file mode 100644 index 00000000..12909e37 --- /dev/null +++ b/arb/dot_ui.c @@ -0,0 +1,74 @@ +/* + Copyright (C) 2021 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 "arb.h" + +void +arb_dot_ui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + ulong v; + unsigned int bc; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + arb_zero(res); + else + { + arb_mul_ui(res, x, y[0], prec); + if (subtract) + arb_neg(res, res); + } + return; + } + else if (len <= 0) + { + arb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + v = y[i * ystep]; + + if (v == 0) + { + ARF_XSIZE(arb_midref(t + i)) = 0; + ARF_EXP(arb_midref(t + i)) = ARF_EXP_ZERO; + } + else + { + count_leading_zeros(bc, v); + + ARF_EXP(arb_midref(t + i)) = FLINT_BITS - bc; + ARF_NOPTR_D(arb_midref(t + i))[0] = v << bc; + ARF_XSIZE(arb_midref(t + i)) = ARF_MAKE_XSIZE(1, 0); + } + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/arb/dot_uiui.c b/arb/dot_uiui.c new file mode 100644 index 00000000..e42a1000 --- /dev/null +++ b/arb/dot_uiui.c @@ -0,0 +1,104 @@ +/* + Copyright (C) 2021 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 "arb.h" + +static void +arf_shallow_set_uiui(arf_t res, ulong vhi, ulong vlo) +{ + unsigned int bc; + + if (vhi == 0) + { + if (vlo == 0) + { + ARF_XSIZE(res) = 0; + ARF_EXP(res) = ARF_EXP_ZERO; + } + else + { + count_leading_zeros(bc, vlo); + ARF_EXP(res) = FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0); + } + } + else if (vlo == 0) + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vhi << bc; + ARF_XSIZE(res) = ARF_MAKE_XSIZE(1, 0); + } + else + { + count_leading_zeros(bc, vhi); + ARF_EXP(res) = 2 * FLINT_BITS - bc; + ARF_NOPTR_D(res)[0] = vlo << bc; + if (bc == 0) + ARF_NOPTR_D(res)[1] = vhi; + else + ARF_NOPTR_D(res)[1] = (vhi << bc) | (vlo >> (FLINT_BITS - bc)); + ARF_XSIZE(res) = ARF_MAKE_XSIZE(2, 0); + } +} + +void +arb_dot_uiui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) +{ + arb_ptr t; + slong i; + ulong vhi, vlo; + TMP_INIT; + + /* todo: fast fma and fmma (len=2) code */ + if (len <= 1) + { + if (initial == NULL) + { + if (len <= 0) + arb_zero(res); + else + { + arf_t t; + arf_shallow_set_uiui(t, y[1], y[0]); + arb_mul_arf(res, x, t, prec); + if (subtract) + arb_neg(res, res); + } + return; + } + else if (len <= 0) + { + arb_set_round(res, initial, prec); + return; + } + } + + TMP_START; + t = TMP_ALLOC(sizeof(arb_struct) * len); + + for (i = 0; i < len; i++) + { + vlo = y[2 * i * ystep]; + vhi = y[2 * i * ystep + 1]; + + arf_shallow_set_uiui(arb_midref(t + i), vhi, vlo); + + MAG_EXP(arb_radref(t + i)) = 0; + MAG_MAN(arb_radref(t + i)) = 0; + } + + arb_dot(res, initial, subtract, x, xstep, t, 1, len, prec); + + TMP_END; +} + diff --git a/arb/test/t-dot_fmpz.c b/arb/test/t-dot_fmpz.c new file mode 100644 index 00000000..1d1b1ea5 --- /dev/null +++ b/arb/test/t-dot_fmpz.c @@ -0,0 +1,102 @@ +/* + Copyright (C) 2021 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 "arb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_siui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + arb_ptr x, y; + fmpz * w; + arb_t s1, s2, z; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _arb_vec_init(len); + y = _arb_vec_init(len); + w = _fmpz_vec_init(len); + arb_init(s1); + arb_init(s2); + arb_init(z); + + for (i = 0; i < len; i++) + { + arb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + fmpz_randtest(w + i, state, 1 + n_randint(state, 200)); + arb_set_fmpz(y + i, w + i); + } + + arb_randtest(s1, state, 200, 10); + arb_randtest(s2, state, 200, 10); + arb_randtest(z, state, 200, 10); + + arb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + arb_dot_fmpz(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + len - 1) : w, revy ? -1 : 1, + len, prec); + + if (!arb_overlaps(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i)); + flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(s1); + arb_clear(s2); + arb_clear(z); + _arb_vec_clear(x, len); + _arb_vec_clear(y, len); + _fmpz_vec_clear(w, len); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb/test/t-dot_si.c b/arb/test/t-dot_si.c new file mode 100644 index 00000000..da91c4de --- /dev/null +++ b/arb/test/t-dot_si.c @@ -0,0 +1,102 @@ +/* + Copyright (C) 2021 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 "arb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + arb_ptr x, y; + slong * w; + arb_t s1, s2, z; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _arb_vec_init(len); + y = _arb_vec_init(len); + w = flint_malloc(sizeof(ulong) * len); + arb_init(s1); + arb_init(s2); + arb_init(z); + + for (i = 0; i < len; i++) + { + arb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[i] = n_randtest(state); + arb_set_si(y + i, w[i]); + } + + arb_randtest(s1, state, 200, 10); + arb_randtest(s2, state, 200, 10); + arb_randtest(z, state, 200, 10); + + arb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + arb_dot_si(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + len - 1) : w, revy ? -1 : 1, + len, prec); + + if (!arb_equal(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i)); + flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(s1); + arb_clear(s2); + arb_clear(z); + _arb_vec_clear(x, len); + _arb_vec_clear(y, len); + flint_free(w); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb/test/t-dot_siui.c b/arb/test/t-dot_siui.c new file mode 100644 index 00000000..3d7dc7af --- /dev/null +++ b/arb/test/t-dot_siui.c @@ -0,0 +1,107 @@ +/* + Copyright (C) 2021 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 "arb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_siui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + arb_ptr x, y; + ulong * w; + arb_t s1, s2, z; + fmpz_t c; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _arb_vec_init(len); + y = _arb_vec_init(len); + w = flint_malloc(2 * sizeof(ulong) * len); + arb_init(s1); + arb_init(s2); + arb_init(z); + fmpz_init(c); + + for (i = 0; i < len; i++) + { + arb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[2 * i] = n_randtest(state); + w[2 * i + 1] = n_randtest(state); + fmpz_set_signed_uiui(c, w[2 * i + 1], w[2 * i]); + arb_set_fmpz(y + i, c); + } + + arb_randtest(s1, state, 200, 10); + arb_randtest(s2, state, 200, 10); + arb_randtest(z, state, 200, 10); + + arb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + arb_dot_siui(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + 2 * len - 2) : w, revy ? -1 : 1, + len, prec); + + if (!arb_overlaps(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i)); + flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(s1); + arb_clear(s2); + arb_clear(z); + _arb_vec_clear(x, len); + _arb_vec_clear(y, len); + flint_free(w); + fmpz_clear(c); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb/test/t-dot_ui.c b/arb/test/t-dot_ui.c new file mode 100644 index 00000000..c2d87eab --- /dev/null +++ b/arb/test/t-dot_ui.c @@ -0,0 +1,102 @@ +/* + Copyright (C) 2021 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 "arb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + arb_ptr x, y; + ulong * w; + arb_t s1, s2, z; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _arb_vec_init(len); + y = _arb_vec_init(len); + w = flint_malloc(sizeof(ulong) * len); + arb_init(s1); + arb_init(s2); + arb_init(z); + + for (i = 0; i < len; i++) + { + arb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[i] = n_randtest(state); + arb_set_ui(y + i, w[i]); + } + + arb_randtest(s1, state, 200, 10); + arb_randtest(s2, state, 200, 10); + arb_randtest(z, state, 200, 10); + + arb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + arb_dot_ui(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + len - 1) : w, revy ? -1 : 1, + len, prec); + + if (!arb_equal(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i)); + flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(s1); + arb_clear(s2); + arb_clear(z); + _arb_vec_clear(x, len); + _arb_vec_clear(y, len); + flint_free(w); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb/test/t-dot_uiui.c b/arb/test/t-dot_uiui.c new file mode 100644 index 00000000..e7103165 --- /dev/null +++ b/arb/test/t-dot_uiui.c @@ -0,0 +1,107 @@ +/* + Copyright (C) 2021 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 "arb.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("dot_uiui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000 * arb_test_multiplier(); iter++) + { + arb_ptr x, y; + ulong * w; + arb_t s1, s2, z; + fmpz_t c; + slong i, len, prec; + int initial, subtract, revx, revy; + + len = n_randint(state, 5); + prec = 2 + n_randint(state, 200); + + initial = n_randint(state, 2); + subtract = n_randint(state, 2); + revx = n_randint(state, 2); + revy = n_randint(state, 2); + + x = _arb_vec_init(len); + y = _arb_vec_init(len); + w = flint_malloc(2 * sizeof(ulong) * len); + arb_init(s1); + arb_init(s2); + arb_init(z); + fmpz_init(c); + + for (i = 0; i < len; i++) + { + arb_randtest(x + i, state, 2 + n_randint(state, 200), 10); + w[2 * i] = n_randtest(state); + w[2 * i + 1] = n_randtest(state); + fmpz_set_uiui(c, w[2 * i + 1], w[2 * i]); + arb_set_fmpz(y + i, c); + } + + arb_randtest(s1, state, 200, 10); + arb_randtest(s2, state, 200, 10); + arb_randtest(z, state, 200, 10); + + arb_dot(s1, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (y + len - 1) : y, revy ? -1 : 1, + len, prec); + + arb_dot_uiui(s2, initial ? z : NULL, subtract, + revx ? (x + len - 1) : x, revx ? -1 : 1, + revy ? (w + 2 * len - 2) : w, revy ? -1 : 1, + len, prec); + + if (!arb_overlaps(s1, s2)) + { + flint_printf("FAIL\n\n"); + flint_printf("iter = %wd, len = %wd, prec = %wd\n\n", iter, len, prec); + + if (initial) + { + flint_printf("z = ", i); arb_printn(z, 100, ARB_STR_MORE); flint_printf(" (%wd)\n\n", arb_bits(z)); + } + + for (i = 0; i < len; i++) + { + flint_printf("x[%wd] = ", i); arb_printn(x + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(x + i)); + flint_printf("y[%wd] = ", i); arb_printn(y + i, 100, ARB_STR_MORE); flint_printf(" (%wd)\n", arb_bits(y + i)); + } + flint_printf("\n\n"); + flint_printf("s1 = "); arb_printn(s1, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_printf("s2 = "); arb_printn(s2, 100, ARB_STR_MORE); flint_printf("\n\n"); + flint_abort(); + } + + arb_clear(s1); + arb_clear(s2); + arb_clear(z); + _arb_vec_clear(x, len); + _arb_vec_clear(y, len); + flint_free(w); + fmpz_clear(c); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arf.h b/arf.h index bdc51d88..c9e62f32 100644 --- a/arf.h +++ b/arf.h @@ -105,7 +105,7 @@ arf_rnd_to_mpfr(arf_rnd_t rnd) #define ARF_XSIZE(x) ((x)->size) /* Construct size field value from size in limbs and sign bit. */ -#define ARF_MAKE_XSIZE(size, sgnbit) ((((mp_size_t) size) << 1) | sgnbit) +#define ARF_MAKE_XSIZE(size, sgnbit) ((((mp_size_t) size) << 1) | (sgnbit)) /* The limb size, and the sign bit. */ #define ARF_SIZE(x) (ARF_XSIZE(x) >> 1) diff --git a/doc/source/acb.rst b/doc/source/acb.rst index 7ceca2dc..84cb9d55 100644 --- a/doc/source/acb.rst +++ b/doc/source/acb.rst @@ -522,10 +522,8 @@ Dot product ------------------------------------------------------------------------------- .. function:: void acb_dot_precise(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec) - -.. function:: void acb_dot_simple(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec) - -.. function:: void acb_dot(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec) + void acb_dot_simple(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec) + void acb_dot(acb_t res, const acb_t s, int subtract, acb_srcptr x, slong xstep, acb_srcptr y, slong ystep, slong len, slong prec) Computes the dot product of the vectors *x* and *y*, setting *res* to `s + (-1)^{subtract} \sum_{i=0}^{len-1} x_i y_i`. @@ -564,6 +562,17 @@ Dot product The radii of the inputs are ignored (only the midpoints are read) and only the midpoint of the output is written. +.. function:: void acb_dot_ui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) + void acb_dot_si(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec) + void acb_dot_uiui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) + void acb_dot_siui(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) + void acb_dot_fmpz(acb_t res, const acb_t initial, int subtract, acb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec) + + Equivalent to :func:`acb_dot`, but with integers in the array *y*. + The *uiui* and *siui* versions take an array of double-limb integers + as input; the *siui* version assumes that these represent signed + integers in two's complement form. + Mathematical constants ------------------------------------------------------------------------------- diff --git a/doc/source/arb.rst b/doc/source/arb.rst index dc1fcc52..31ba1f90 100644 --- a/doc/source/arb.rst +++ b/doc/source/arb.rst @@ -866,10 +866,8 @@ Dot product ------------------------------------------------------------------------------- .. function:: void arb_dot_precise(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec) - -.. function:: void arb_dot_simple(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec) - -.. function:: void arb_dot(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec) + void arb_dot_simple(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec) + void arb_dot(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec) Computes the dot product of the vectors *x* and *y*, setting *res* to `s + (-1)^{subtract} \sum_{i=0}^{len-1} x_i y_i`. @@ -908,6 +906,18 @@ Dot product The radii of the inputs are ignored (only the midpoints are read) and only the midpoint of the output is written. +.. function:: void arb_dot_ui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) + void arb_dot_si(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const slong * y, slong ystep, slong len, slong prec) + void arb_dot_uiui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) + void arb_dot_siui(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const ulong * y, slong ystep, slong len, slong prec) + void arb_dot_fmpz(arb_t res, const arb_t initial, int subtract, arb_srcptr x, slong xstep, const fmpz * y, slong ystep, slong len, slong prec) + + Equivalent to :func:`arb_dot`, but with integers in the array *y*. + The *uiui* and *siui* versions take an array of double-limb integers + as input; the *siui* version assumes that these represent signed + integers in two's complement form. + + Powers and roots -------------------------------------------------------------------------------