diff --git a/acb.h b/acb.h index 2924df3e..e7beca38 100644 --- a/acb.h +++ b/acb.h @@ -460,21 +460,47 @@ acb_mul_2exp_si(acb_t z, const acb_t x, long e) static __inline__ void acb_addmul(acb_t z, const acb_t x, const acb_t y, long prec) { - acb_t t; - acb_init(t); - acb_mul(t, x, y, prec); - acb_add(z, z, t, prec); - acb_clear(t); + if (arb_is_zero(acb_imagref(y))) + { + arb_addmul(acb_imagref(z), acb_imagref(x), acb_realref(y), prec); + arb_addmul(acb_realref(z), acb_realref(x), acb_realref(y), prec); + } + else if (arb_is_zero(acb_imagref(x))) + { + arb_addmul(acb_imagref(z), acb_imagref(y), acb_realref(x), prec); + arb_addmul(acb_realref(z), acb_realref(y), acb_realref(x), prec); + } + else + { + acb_t t; + acb_init(t); + acb_mul(t, x, y, prec); + acb_add(z, z, t, prec); + acb_clear(t); + } } static __inline__ void acb_submul(acb_t z, const acb_t x, const acb_t y, long prec) { - acb_t t; - acb_init(t); - acb_mul(t, x, y, prec); - acb_sub(z, z, t, prec); - acb_clear(t); + if (arb_is_zero(acb_imagref(y))) + { + arb_submul(acb_imagref(z), acb_imagref(x), acb_realref(y), prec); + arb_submul(acb_realref(z), acb_realref(x), acb_realref(y), prec); + } + else if (arb_is_zero(acb_imagref(x))) + { + arb_submul(acb_imagref(z), acb_imagref(y), acb_realref(x), prec); + arb_submul(acb_realref(z), acb_realref(y), acb_realref(x), prec); + } + else + { + acb_t t; + acb_init(t); + acb_mul(t, x, y, prec); + acb_sub(z, z, t, prec); + acb_clear(t); + } } static __inline__ void @@ -681,35 +707,17 @@ _acb_vec_sub(acb_ptr res, acb_srcptr vec1, acb_srcptr vec2, long len, long prec) static __inline__ void _acb_vec_scalar_submul(acb_ptr res, acb_srcptr vec, long len, const acb_t c, long prec) { - if (len > 0) - { - long i; - acb_t t; - acb_init(t); - for (i = 0; i < len; i++) - { - acb_mul(t, vec + i, c, prec); - acb_sub(res + i, res + i, t, prec); - } - acb_clear(t); - } + long i; + for (i = 0; i < len; i++) + acb_submul(res + i, vec + i, c, prec); } static __inline__ void _acb_vec_scalar_addmul(acb_ptr res, acb_srcptr vec, long len, const acb_t c, long prec) { - if (len > 0) - { - long i; - acb_t t; - acb_init(t); - for (i = 0; i < len; i++) - { - acb_mul(t, vec + i, c, prec); - acb_add(res + i, res + i, t, prec); - } - acb_clear(t); - } + long i; + for (i = 0; i < len; i++) + acb_addmul(res + i, vec + i, c, prec); } static __inline__ void diff --git a/arb_poly.h b/arb_poly.h index 1a91b494..0669d84f 100644 --- a/arb_poly.h +++ b/arb_poly.h @@ -27,6 +27,7 @@ #define ARB_POLY_H #include "arb.h" +#include "acb.h" #include "fmpz_poly.h" #include "fmpq_poly.h" @@ -526,29 +527,28 @@ void arb_poly_compose_series_brent_kung(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, long n, long prec); + +void _arb_poly_evaluate_acb_horner(acb_t res, arb_srcptr f, long len, const acb_t x, long prec); +void arb_poly_evaluate_acb_horner(acb_t res, const arb_poly_t f, const acb_t a, long prec); + +void _arb_poly_evaluate_acb_rectangular(acb_t y, arb_srcptr poly, long len, const acb_t x, long prec); +void arb_poly_evaluate_acb_rectangular(acb_t res, const arb_poly_t f, const acb_t a, long prec); + +void _arb_poly_evaluate_acb(acb_t res, arb_srcptr f, long len, const acb_t x, long prec); +void arb_poly_evaluate_acb(acb_t res, const arb_poly_t f, const acb_t a, long prec); + +void _arb_poly_evaluate2_acb_horner(acb_t y, acb_t z, arb_srcptr f, long len, const acb_t x, long prec); +void arb_poly_evaluate2_acb_horner(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, long prec); + +void _arb_poly_evaluate2_acb_rectangular(acb_t y, acb_t z, arb_srcptr f, long len, const acb_t x, long prec); +void arb_poly_evaluate2_acb_rectangular(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, long prec); + +void _arb_poly_evaluate2_acb(acb_t y, acb_t z, arb_srcptr f, long len, const acb_t x, long prec); +void arb_poly_evaluate2_acb(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, long prec); + /* TBD: -void _arb_poly_evaluate_fmpcb_horner(fmpcb_t res, arb_srcptr f, long len, const fmpcb_t x, long prec); -void arb_poly_evaluate_fmpcb_horner(fmpcb_t res, const arb_poly_t f, const fmpcb_t a, long prec); - -void _arb_poly_evaluate_fmpcb_rectangular(fmpcb_t y, arb_srcptr poly, long len, const fmpcb_t x, long prec); -void arb_poly_evaluate_fmpcb_rectangular(fmpcb_t res, const arb_poly_t f, const fmpcb_t a, long prec); - -void _arb_poly_evaluate_fmpcb(fmpcb_t res, arb_srcptr f, long len, const fmpcb_t x, long prec); -void arb_poly_evaluate_fmpcb(fmpcb_t res, const arb_poly_t f, const fmpcb_t a, long prec); - -void _arb_poly_evaluate2_fmpcb_horner(fmpcb_t y, fmpcb_t z, arb_srcptr f, long len, const fmpcb_t x, long prec); -void arb_poly_evaluate2_fmpcb_horner(fmpcb_t y, fmpcb_t z, const arb_poly_t f, const fmpcb_t x, long prec); - -void _arb_poly_evaluate2_fmpcb_rectangular(fmpcb_t y, fmpcb_t z, arb_srcptr f, long len, const fmpcb_t x, long prec); -void arb_poly_evaluate2_fmpcb_rectangular(fmpcb_t y, fmpcb_t z, const arb_poly_t f, const fmpcb_t x, long prec); - -void _arb_poly_evaluate2_fmpcb(fmpcb_t y, fmpcb_t z, arb_srcptr f, long len, const fmpcb_t x, long prec); -void arb_poly_evaluate2_fmpcb(fmpcb_t y, fmpcb_t z, const arb_poly_t f, const fmpcb_t x, long prec); - - - void _arb_poly_gamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long prec); void arb_poly_gamma_series(arb_poly_t res, const arb_poly_t f, long n, long prec); diff --git a/arb_poly/evaluate2_acb.c b/arb_poly/evaluate2_acb.c new file mode 100644 index 00000000..b0c126d2 --- /dev/null +++ b/arb_poly/evaluate2_acb.c @@ -0,0 +1,39 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +void +_arb_poly_evaluate2_acb(acb_t y, acb_t z, arb_srcptr f, long len, const acb_t x, long prec) +{ + _arb_poly_evaluate2_acb_rectangular(y, z, f, len, x, prec); +} + +void +arb_poly_evaluate2_acb(acb_t r, acb_t s, const arb_poly_t f, const acb_t a, long prec) +{ + _arb_poly_evaluate2_acb(r, s, f->coeffs, f->length, a, prec); +} + diff --git a/arb_poly/evaluate2_acb_horner.c b/arb_poly/evaluate2_acb_horner.c new file mode 100644 index 00000000..91a070f3 --- /dev/null +++ b/arb_poly/evaluate2_acb_horner.c @@ -0,0 +1,87 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +void +_arb_poly_evaluate2_acb_horner(acb_t y, acb_t z, + arb_srcptr poly, long len, const acb_t x, long prec) +{ + if (len == 0) + { + acb_zero(y); + acb_zero(z); + } + else if (len == 1) + { + acb_set_round_arb(y, poly + 0, prec); + acb_zero(z); + } + else if (acb_is_zero(x)) + { + acb_set_round_arb(y, poly + 0, prec); + acb_set_round_arb(z, poly + 1, prec); + } + else if (len == 2) + { + acb_mul_arb(y, x, poly + 1, prec); + acb_add_arb(y, y, poly + 0, prec); + acb_set_round_arb(z, poly + 1, prec); + } + else + { + acb_t t, u, v; + long i; + + acb_init(t); + acb_init(u); + acb_init(v); + + acb_set_round_arb(u, poly + len - 1, prec); + acb_zero(v); + + for (i = len - 2; i >= 0; i--) + { + acb_mul(t, v, x, prec); + acb_add(v, u, t, prec); + acb_mul(t, u, x, prec); + acb_add_arb(u, t, poly + i, prec); + } + + acb_swap(y, u); + acb_swap(z, v); + + acb_clear(t); + acb_clear(u); + acb_clear(v); + } +} + +void +arb_poly_evaluate2_acb_horner(acb_t r, acb_t s, const arb_poly_t f, const acb_t a, long prec) +{ + _arb_poly_evaluate2_acb_horner(r, s, f->coeffs, f->length, a, prec); +} + diff --git a/arb_poly/evaluate2_acb_rectangular.c b/arb_poly/evaluate2_acb_rectangular.c new file mode 100644 index 00000000..fda3f7f2 --- /dev/null +++ b/arb_poly/evaluate2_acb_rectangular.c @@ -0,0 +1,120 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +void +_arb_poly_evaluate2_acb_rectangular(acb_t y, acb_t z, + arb_srcptr poly, long len, const acb_t x, long prec) +{ + long i, j, m, r; + acb_ptr xs; + acb_t s, t; + arb_t c; + + if (len < 3) + { + if (len == 0) + { + acb_zero(y); + acb_zero(z); + } + else if (len == 1) + { + acb_set_round_arb(y, poly + 0, prec); + acb_zero(z); + } + else if (len == 2) + { + acb_mul_arb(y, x, poly + 1, prec); + acb_add_arb(y, y, poly + 0, prec); + acb_set_round_arb(z, poly + 1, prec); + } + return; + } + + m = n_sqrt(len) + 1; + m *= 1; + + r = (len + m - 1) / m; + + xs = _acb_vec_init(m + 1); + acb_init(s); + acb_init(t); + arb_init(c); + + _acb_vec_set_powers(xs, x, m + 1, prec); + + acb_set_arb(y, poly + (r - 1) * m); + for (j = 1; (r - 1) * m + j < len; j++) + acb_addmul_arb(y, xs + j, poly + (r - 1) * m + j, prec); + + for (i = r - 2; i >= 0; i--) + { + acb_set_arb(s, poly + i * m); + for (j = 1; j < m; j++) + acb_addmul_arb(s, xs + j, poly + i * m + j, prec); + + acb_mul(y, y, xs + m, prec); + acb_add(y, y, s, prec); + } + + len -= 1; + r = (len + m - 1) / m; + arb_mul_ui(acb_realref(z), poly + (r - 1) * m + 1, (r - 1) * m + 1, FMPR_PREC_EXACT); + arb_zero(acb_imagref(z)); + for (j = 1; (r - 1) * m + j < len; j++) + { + arb_mul_ui(c, poly + (r - 1) * m + j + 1, (r - 1) * m + j + 1, FMPR_PREC_EXACT); + acb_addmul_arb(z, xs + j, c, prec); + } + + for (i = r - 2; i >= 0; i--) + { + arb_mul_ui(acb_realref(s), poly + i * m + 1, i * m + 1, FMPR_PREC_EXACT); + arb_zero(acb_imagref(s)); + + for (j = 1; j < m; j++) + { + arb_mul_ui(c, poly + i * m + j + 1, i * m + j + 1, FMPR_PREC_EXACT); + acb_addmul_arb(s, xs + j, c, prec); + } + + acb_mul(z, z, xs + m, prec); + acb_add(z, z, s, prec); + } + + _acb_vec_clear(xs, m + 1); + acb_clear(s); + acb_clear(t); + arb_clear(c); +} + +void +arb_poly_evaluate2_acb_rectangular(acb_t r, acb_t s, const arb_poly_t f, const acb_t a, long prec) +{ + _arb_poly_evaluate2_acb_rectangular(r, s, f->coeffs, f->length, a, prec); +} + diff --git a/arb_poly/evaluate_acb.c b/arb_poly/evaluate_acb.c new file mode 100644 index 00000000..2ffa2361 --- /dev/null +++ b/arb_poly/evaluate_acb.c @@ -0,0 +1,40 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +void +_arb_poly_evaluate_acb(acb_t res, arb_srcptr f, long len, + const acb_t x, long prec) +{ + _arb_poly_evaluate_acb_rectangular(res, f, len, x, prec); +} + +void +arb_poly_evaluate_acb(acb_t res, const arb_poly_t f, const acb_t a, long prec) +{ + _arb_poly_evaluate_acb(res, f->coeffs, f->length, a, prec); +} + diff --git a/arb_poly/evaluate_acb_horner.c b/arb_poly/evaluate_acb_horner.c new file mode 100644 index 00000000..576f0f67 --- /dev/null +++ b/arb_poly/evaluate_acb_horner.c @@ -0,0 +1,73 @@ +/*============================================================================= + + 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) 2010 Sebastian Pancratz + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +void +_arb_poly_evaluate_acb_horner(acb_t y, arb_srcptr f, long len, + const acb_t x, long prec) +{ + if (len == 0) + { + acb_zero(y); + } + else if (len == 1 || acb_is_zero(x)) + { + acb_set_round_arb(y, f, prec); + } + else if (len == 2) + { + acb_mul_arb(y, x, f + 1, prec); + acb_add_arb(y, y, f + 0, prec); + } + else + { + long i = len - 1; + acb_t t, u; + + acb_init(t); + acb_init(u); + acb_set_arb(u, f + i); + + for (i = len - 2; i >= 0; i--) + { + acb_mul(t, u, x, prec); + acb_add_arb(u, t, f + i, prec); + } + + acb_swap(y, u); + + acb_clear(t); + acb_clear(u); + } +} + +void +arb_poly_evaluate_acb_horner(acb_t res, const arb_poly_t f, const acb_t a, long prec) +{ + _arb_poly_evaluate_acb_horner(res, f->coeffs, f->length, a, prec); +} + diff --git a/arb_poly/evaluate_acb_rectangular.c b/arb_poly/evaluate_acb_rectangular.c new file mode 100644 index 00000000..8d6500df --- /dev/null +++ b/arb_poly/evaluate_acb_rectangular.c @@ -0,0 +1,89 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +void +_arb_poly_evaluate_acb_rectangular(acb_t y, arb_srcptr poly, + long len, const acb_t x, long prec) +{ + long i, j, m, r; + acb_ptr xs; + acb_t s, t, c; + + if (len < 3) + { + if (len == 0) + { + acb_zero(y); + } + else if (len == 1) + { + acb_set_round_arb(y, poly + 0, prec); + } + else if (len == 2) + { + acb_mul_arb(y, x, poly + 1, prec); + acb_add_arb(y, y, poly + 0, prec); + } + return; + } + + m = n_sqrt(len) + 1; + r = (len + m - 1) / m; + + xs = _acb_vec_init(m + 1); + acb_init(s); + acb_init(t); + acb_init(c); + + _acb_vec_set_powers(xs, x, m + 1, prec); + + acb_set_arb(y, poly + (r - 1) * m); + for (j = 1; (r - 1) * m + j < len; j++) + acb_addmul_arb(y, xs + j, poly + (r - 1) * m + j, prec); + + for (i = r - 2; i >= 0; i--) + { + acb_set_arb(s, poly + i * m); + for (j = 1; j < m; j++) + acb_addmul_arb(s, xs + j, poly + i * m + j, prec); + + acb_mul(y, y, xs + m, prec); + acb_add(y, y, s, prec); + } + + _acb_vec_clear(xs, m + 1); + acb_clear(s); + acb_clear(t); + acb_clear(c); +} + +void +arb_poly_evaluate_acb_rectangular(acb_t res, const arb_poly_t f, const acb_t a, long prec) +{ + _arb_poly_evaluate_acb_rectangular(res, f->coeffs, f->length, a, prec); +} + diff --git a/arb_poly/test/t-evaluate2_acb_rectangular.c b/arb_poly/test/t-evaluate2_acb_rectangular.c new file mode 100644 index 00000000..580a1bfc --- /dev/null +++ b/arb_poly/test/t-evaluate2_acb_rectangular.c @@ -0,0 +1,87 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("evaluate2_acb_rectangular...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arb_poly_t f, g; + acb_t x, y1, z1, y2, z2; + + acb_init(x); + acb_init(y1); + acb_init(z1); + acb_init(y2); + acb_init(z2); + arb_poly_init(f); + arb_poly_init(g); + + acb_randtest(x, state, 2 + n_randint(state, 1000), 5); + arb_poly_randtest(f, state, 2 + n_randint(state, 100), 2 + n_randint(state, 1000), 5); + arb_poly_derivative(g, f, 2 + n_randint(state, 1000)); + + arb_poly_evaluate2_acb_rectangular(y1, z1, f, x, 2 + n_randint(state, 1000)); + + arb_poly_evaluate_acb_horner(y2, f, x, 2 + n_randint(state, 1000)); + arb_poly_evaluate_acb_horner(z2, g, x, 2 + n_randint(state, 1000)); + + if (!acb_overlaps(y1, y2) || !acb_overlaps(z1, z2)) + { + printf("FAIL\n\n"); + printf("f = "); arb_poly_printd(f, 15); printf("\n\n"); + printf("g = "); arb_poly_printd(g, 15); printf("\n\n"); + printf("x = "); acb_printd(x, 15); printf("\n\n"); + printf("y1 = "); acb_printd(y1, 15); printf("\n\n"); + printf("z1 = "); acb_printd(z1, 15); printf("\n\n"); + printf("y2 = "); acb_printd(y2, 15); printf("\n\n"); + printf("z2 = "); acb_printd(z2, 15); printf("\n\n"); + abort(); + } + + arb_poly_clear(f); + arb_poly_clear(g); + acb_clear(x); + acb_clear(y1); + acb_clear(z1); + acb_clear(y2); + acb_clear(z2); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb_poly/test/t-evaluate_acb_rectangular.c b/arb_poly/test/t-evaluate_acb_rectangular.c new file mode 100644 index 00000000..0d265b9f --- /dev/null +++ b/arb_poly/test/t-evaluate_acb_rectangular.c @@ -0,0 +1,75 @@ +/*============================================================================= + + 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, 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("evaluate_acb_rectangular...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arb_poly_t f; + acb_t x, y1, y2; + + acb_init(x); + acb_init(y1); + acb_init(y2); + arb_poly_init(f); + + acb_randtest(x, state, 2 + n_randint(state, 1000), 5); + arb_poly_randtest(f, state, 2 + n_randint(state, 100), 2 + n_randint(state, 1000), 5); + + arb_poly_evaluate_acb_rectangular(y1, f, x, 2 + n_randint(state, 1000)); + arb_poly_evaluate_acb_horner(y2, f, x, 2 + n_randint(state, 1000)); + + if (!acb_overlaps(y1, y2)) + { + printf("FAIL\n\n"); + printf("f = "); arb_poly_printd(f, 15); printf("\n\n"); + printf("x = "); acb_printd(x, 15); printf("\n\n"); + printf("y1 = "); acb_printd(y1, 15); printf("\n\n"); + printf("y2 = "); acb_printd(y2, 15); printf("\n\n"); + abort(); + } + + arb_poly_clear(f); + acb_clear(x); + acb_clear(y1); + acb_clear(y2); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +