From 3a6ef78f789cb29a92e2ffb9dcf56642866556af Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sun, 8 Apr 2012 02:44:20 +0200 Subject: [PATCH] add arb_poly power series inversion --- arb_poly.h | 2 - arb_poly/inv_series.c | 123 +++++++++++++++++++++++++++++++++ arb_poly/test/t-inv_series.c | 127 +++++++++++++++++++++++++++++++++++ 3 files changed, 250 insertions(+), 2 deletions(-) create mode 100644 arb_poly/inv_series.c create mode 100644 arb_poly/test/t-inv_series.c diff --git a/arb_poly.h b/arb_poly.h index bf613ee2..0a24fc56 100644 --- a/arb_poly.h +++ b/arb_poly.h @@ -108,8 +108,6 @@ void arb_poly_neg(arb_poly_t z, const arb_poly_t x); void arb_poly_mul(arb_poly_t z, const arb_poly_t x, const arb_poly_t y); void arb_poly_mullow(arb_poly_t z, const arb_poly_t x, const arb_poly_t y, long n); -/* void arb_poly_inv_series(arb_poly_t z, const arb_poly_t x, long n); -*/ #endif diff --git a/arb_poly/inv_series.c b/arb_poly/inv_series.c new file mode 100644 index 00000000..4d929dc0 --- /dev/null +++ b/arb_poly/inv_series.c @@ -0,0 +1,123 @@ +/*============================================================================= + + This file is part of FLINT. + + FLINT 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. + + FLINT 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 FLINT; 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 "arb_poly.h" + +/* todo: fix arb_poly_add so this isn't needed */ +static void +_arb_poly_add_two(arb_poly_t t) +{ + if (fmpz_sgn(arb_poly_expref(t)) <= 0) + { + fmpz_t two; + fmpz_init(two); + fmpz_set_ui(two, 2UL); + fmpz_mul_2exp(two, two, -fmpz_get_si(arb_poly_expref(t))); + fmpz_add(t->coeffs, t->coeffs, two); + fmpz_clear(two); + } + else if (*arb_poly_expref(t) == 1) + { + fmpz_add_ui(t->coeffs, t->coeffs, 1); + } + else + { + fmpz_add_ui(arb_poly_radref(t), arb_poly_radref(t), 1); + } +} + +void +arb_poly_inv_series(arb_poly_t z, const arb_poly_t x, long n) +{ + long a[FLINT_BITS]; + long i, m; + arb_poly_t t, u; + arb_t one, r; + long xlen; + + xlen = x->length; + + if (xlen == 0) + { + printf("arb_poly_inv_series: division by zero\n"); + abort(); + } + + if (z == x) + { + arb_poly_t t; + arb_poly_init(t, arb_poly_prec(x)); + arb_poly_set(t, x); + arb_poly_inv_series(z, t, n); + arb_poly_clear(t); + return; + } + + a[i = 0] = n; + while (n >= 2) + a[++i] = (n = (n + 1) / 2); + + arb_poly_init(t, arb_poly_prec(z)); + arb_poly_init(u, arb_poly_prec(z)); + + /* first coefficient */ + arb_init(r, arb_poly_prec(z)); + arb_init(one, arb_poly_prec(z)); + arb_set_si(one, 1); + + fmpz_set(arb_midref(r), arb_poly_coeffs(x)); + fmpz_set(arb_radref(r), arb_poly_radref(x)); + fmpz_set(arb_expref(r), arb_poly_expref(x)); + + arb_div(r, one, r); + + _arb_poly_fit_length(z, 1); + fmpz_set(arb_poly_coeffs(z), arb_midref(r)); + fmpz_set(arb_poly_radref(z), arb_radref(r)); + fmpz_set(arb_poly_expref(z), arb_expref(r)); + z->length = 1; + + arb_clear(r); + arb_clear(one); + /* done first coefficient */ + + /* Newton iteration */ + for (i--; i >= 0; i--) + { + m = n; + n = a[i]; + + arb_poly_mullow(t, x, z, n); + arb_poly_neg(t, t); + + _arb_poly_add_two(t); + + arb_poly_mullow(u, t, z, n); + arb_poly_set(z, u); + } + + arb_poly_clear(t); + arb_poly_clear(u); +} diff --git a/arb_poly/test/t-inv_series.c b/arb_poly/test/t-inv_series.c new file mode 100644 index 00000000..00b2d265 --- /dev/null +++ b/arb_poly/test/t-inv_series.c @@ -0,0 +1,127 @@ +/*============================================================================= + + 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 "arb_poly.h" + +static int +is_invertible(arb_poly_t poly) +{ + if (poly->length == 0) + return 0; + + return !(fmpz_cmpabs(arb_poly_radref(poly), arb_poly_coeffs(poly)) >= 0); +} + +int main() +{ + long iter; + flint_rand_t state; + + printf("inv_series...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arb_poly_t a, b; + fmpq_poly_t x, y; + long n; + + arb_poly_init(a, 1 + n_randint(state, 200)); + arb_poly_init(b, 1 + n_randint(state, 200)); + + fmpq_poly_init(x); + fmpq_poly_init(y); + + do { + arb_poly_randtest(a, state, n_randint(state, 20), 9); + } while (!is_invertible(a)); + + arb_poly_randtest(b, state, n_randint(state, 20), 9); + + arb_poly_get_rand_fmpq_poly(x, state, a); + arb_poly_get_rand_fmpq_poly(y, state, b); + + n = 1 + n_randint(state, 20); + + arb_poly_inv_series(b, a, n); + fmpq_poly_inv_series(y, x, n); + + if (!arb_poly_contains_fmpq_poly(b, y)) + { + printf("FAIL: containment\n\n"); + printf("a = "); arb_poly_debug(a); printf("\n\n"); + printf("x = "); fmpq_poly_print(x); printf("\n\n"); + printf("b = "); arb_poly_debug(b); printf("\n\n"); + printf("y = "); fmpq_poly_print(y); printf("\n\n"); + abort(); + } + + arb_poly_clear(a); + arb_poly_clear(b); + + fmpq_poly_clear(x); + fmpq_poly_clear(y); + } + + /* aliasing */ + for (iter = 0; iter < 10000; iter++) + { + arb_poly_t a; + fmpq_poly_t x; + long n; + + arb_poly_init(a, 1 + n_randint(state, 200)); + fmpq_poly_init(x); + + do { + arb_poly_randtest(a, state, n_randint(state, 20), 9); + } while (!is_invertible(a)); + + arb_poly_get_rand_fmpq_poly(x, state, a); + + n = 1 + n_randint(state, 20); + + arb_poly_inv_series(a, a, n); + fmpq_poly_inv_series(x, x, n); + + if (!arb_poly_contains_fmpq_poly(a, x)) + { + printf("FAIL: aliasing\n\n"); + printf("a = "); arb_poly_debug(a); printf("\n\n"); + printf("x = "); fmpq_poly_print(x); printf("\n\n"); + abort(); + } + + arb_poly_clear(a); + fmpq_poly_clear(x); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}