improve _fmpz_poly_inv_series, especially for short input (includes i…

…nterface change)
fredrik-johansson · Apr 4, 2014 · 83d3bf0 · 83d3bf0
1 parent 10e5d44
commit 83d3bf0
Show file tree

Hide file tree

Showing 15 changed files with 506 additions and 91 deletions.
diff --git a/arith/number_of_partitions_vec.c b/arith/number_of_partitions_vec.c
@@ -51,7 +51,7 @@ arith_number_of_partitions_vec(fmpz * res, slong len)
     if (n + k < len) tmp[n + k] = WORD(-1);
     if (n + 3*k + 1 < len) tmp[n + 3*k + 1] = WORD(1);
 
-    _fmpz_poly_inv_series(res, tmp, len);
+    _fmpz_poly_inv_series(res, tmp, len, len);
 
     _fmpz_vec_clear(tmp, len);
 }
diff --git a/flint.h b/flint.h
@@ -319,6 +319,30 @@ void mpn_tdiv_q(mp_ptr qp,
       __tmp_root = __tmp_root->next; \
    }
 
+/* Newton iteration macros */
+#define FLINT_NEWTON_INIT(from, to) \
+    { \
+        slong __steps[FLINT_BITS], __i, __from, __to; \
+        __steps[__i = 0] = __to = (to); \
+        __from = (from); \
+        while (__to > __from) \
+            __steps[++__i] = (__to = (__to + 1) / 2); \
+
+#define FLINT_NEWTON_BASECASE(bc_to) { slong bc_to = __to;
+
+#define FLINT_NEWTON_END_BASECASE }
+
+#define FLINT_NEWTON_LOOP(step_from, step_to) \
+        { \
+            for (__i--; __i >= 0; __i--) \
+            { \
+                slong step_from = __steps[__i+1]; \
+                slong step_to = __steps[__i]; \
+
+#define FLINT_NEWTON_END_LOOP }}
+
+#define FLINT_NEWTON_END }
+
 int parse_fmt(int * floating, const char * fmt);
 
 size_t flint_printf(const char * str, ...); /* flint version of printf */

diff --git a/fmpz_poly.h b/fmpz_poly.h
@@ -47,7 +47,7 @@
  extern "C" {
 #endif
 
-#define FMPZ_POLY_INV_NEWTON_CUTOFF  32
+#define FMPZ_POLY_INV_NEWTON_CUTOFF 32
 
 /*  Type definitions *********************************************************/
 
@@ -717,21 +717,17 @@ _fmpz_poly_div_root(fmpz * Q, const fmpz * A, slong len, const fmpz_t c);
 
 /*  Power series division  ***************************************************/
 
-void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong n);
+void _fmpz_poly_inv_series_basecase(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n);
+
+void fmpz_poly_inv_series_basecase(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n);
+
+void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n);
 
 void fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n);
 
-static __inline__ void 
-_fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong n)
-{
-    _fmpz_poly_inv_series_newton(Qinv, Q, n);
-}
+void _fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n);
 
-static __inline__ void 
-fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n)
-{
-    fmpz_poly_inv_series_newton(Qinv, Q, n);
-}
+void fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n);
 
 void _fmpz_poly_div_series(fmpz * Q, const fmpz * A, const fmpz * B, slong n);
 

diff --git a/fmpz_poly/div_series.c b/fmpz_poly/div_series.c
@@ -41,7 +41,7 @@ _fmpz_poly_div_series(fmpz * Q, const fmpz * A, const fmpz * B, slong n)
     {
         fmpz * Binv = _fmpz_vec_init(n);
 
-        _fmpz_poly_inv_series(Binv, B, n);
+        _fmpz_poly_inv_series(Binv, B, n, n);
         _fmpz_poly_mullow(Q, A, n, Binv, n, n);
 
         _fmpz_vec_clear(Binv, n);

diff --git a/fmpz_poly/doc/fmpz_poly.txt b/fmpz_poly/doc/fmpz_poly.txt
@@ -1694,31 +1694,45 @@ int fmpz_poly_divides(fmpz_poly_t Q,
 
 *******************************************************************************
 
-void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong n)
+void _fmpz_poly_inv_series_basecase(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n)
 
-    Computes the first $n$ terms of the inverse power series of $Q$ using 
-    Newton iteration.
+    Computes the first $n$ terms of the inverse power series of
+    \code{(Q, lenQ)} using a recurrence.
 
-    Assumes that $n \geq 1$, that $Q$ has length at least $n$ and constant 
-    term~$\pm 1$.  Does not support aliasing.
+    Assumes that $n \geq 1$ and that $Q$ has constant term~$\pm 1$.
+    Does not support aliasing.
 
-id fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n)
+void fmpz_poly_inv_series_basecase(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n)
 
     Computes the first $n$ terms of the inverse power series of $Q$ 
-    using Newton iteration, assuming that $Q$ has constant term~$\pm 1$ 
+    using a recurrence, assuming that $Q$ has constant term~$\pm 1$ 
     and $n \geq 1$.
 
+void _fmpz_poly_inv_series_newton(fmpz * Qinv, const fmpz * Q, slong n)
+
+    Computes the first $n$ terms of the inverse power series of
+    \code{(Q, lenQ)} using Newton iteration.
+
+    Assumes that $n \geq 1$ and that $Q$ has constant term~$\pm 1$.
+    Does not support aliasing.
+
+void fmpz_poly_inv_series_newton(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong Qlen, slong n)
+
+    Computes the first $n$ terms of the inverse power series of $Q$ using
+    Newton iteration, assuming $Q$ has constant term~$\pm 1$ and $n \geq 1$.
+
 void _fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong n)
 
-    Computes the first $n$ terms of the inverse power series of $Q$.
+    Computes the first $n$ terms of the inverse power series of
+    \code{(Q, lenQ)}.
 
-    Assumes that $n \geq 1$, that $Q$ has length at least $n$ and constant 
-    term~$1$.  Does not support aliasing.
+    Assumes that $n \geq 1$ and that $Q$ has constant term~$\pm 1$.
+    Does not support aliasing.
 
 void fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n)
 
     Computes the first $n$ terms of the inverse power series of $Q$, 
-    assuming $Q$ has constant term~$1$ and $n \geq 1$.
+    assuming $Q$ has constant term~$\pm 1$ and $n \geq 1$.
 
 void _fmpz_poly_div_series(fmpz * Q, const fmpz * A, const fmpz * B)
 

diff --git a/fmpz_poly/eta_qexp.c b/fmpz_poly/eta_qexp.c
@@ -176,7 +176,7 @@ _fmpz_poly_eta_qexp(fmpz * f, slong e, slong len)
     {
         fmpz * t = _fmpz_vec_init(len);
         _fmpz_poly_eta_qexp(t, -e, len);
-        _fmpz_poly_inv_series(f, t, len);
+        _fmpz_poly_inv_series(f, t, len, len);
         _fmpz_vec_clear(t, len);
         return;
     }

diff --git a/fmpz_poly/inv_series.c b/fmpz_poly/inv_series.c
@@ -0,0 +1,67 @@
+/*=============================================================================
+
+    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) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpz_poly.h"
+
+void
+_fmpz_poly_inv_series(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n)
+{
+    if (Qlen <= 8 || n <= 24)
+        _fmpz_poly_inv_series_basecase(Qinv, Q, Qlen, n);
+    else
+        _fmpz_poly_inv_series_newton(Qinv, Q, Qlen, n);
+}
+
+void
+fmpz_poly_inv_series(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n)
+{
+    slong Qlen = Q->length;
+
+    Qlen = FLINT_MIN(Qlen, n);
+
+    if (Qlen == 0)
+    {
+        flint_printf("Exception (fmpz_poly_inv_series). Division by zero.\n");
+        abort();
+    }
+
+    if (Qinv != Q)
+    {
+        fmpz_poly_fit_length(Qinv, n);
+        _fmpz_poly_inv_series(Qinv->coeffs, Q->coeffs, Qlen, n);
+    }
+    else
+    {
+        fmpz_poly_t t;
+        fmpz_poly_init2(t, n);
+        _fmpz_poly_inv_series(t->coeffs, Q->coeffs, Qlen, n);
+        fmpz_poly_swap(Qinv, t);
+        fmpz_poly_clear(t);
+    }
+
+    _fmpz_poly_set_length(Qinv, n);
+    _fmpz_poly_normalise(Qinv);
+}
+
diff --git a/fmpz_poly/inv_series_basecase.c b/fmpz_poly/inv_series_basecase.c
@@ -0,0 +1,86 @@
+/*=============================================================================
+
+    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) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpz_poly.h"
+
+void
+_fmpz_poly_inv_series_basecase(fmpz * Qinv, const fmpz * Q, slong Qlen, slong n)
+{
+    Qlen = FLINT_MIN(Qlen, n);
+
+    fmpz_set(Qinv, Q);
+
+    if (Qlen == 1)
+    {
+        _fmpz_vec_zero(Qinv + 1, n - 1);
+    }
+    else
+    {
+        slong i, j;
+
+        for (i = 1; i < n; i++)
+        {
+            fmpz_mul(Qinv + i, Q + 1, Qinv + i - 1);
+
+            for (j = 2; j < FLINT_MIN(i + 1, Qlen); j++)
+                fmpz_addmul(Qinv + i, Q + j, Qinv + i - j);
+
+            if (fmpz_is_one(Qinv))
+                fmpz_neg(Qinv + i, Qinv + i);
+        }
+    }
+}
+
+void
+fmpz_poly_inv_series_basecase(fmpz_poly_t Qinv, const fmpz_poly_t Q, slong n)
+{
+    slong Qlen = Q->length;
+
+    Qlen = FLINT_MIN(Qlen, n);
+
+    if (Qlen == 0)
+    {
+        flint_printf("Exception (fmpz_poly_inv_series_basecase). Division by zero.\n");
+        abort();
+    }
+
+    if (Qinv != Q)
+    {
+        fmpz_poly_fit_length(Qinv, n);
+        _fmpz_poly_inv_series_basecase(Qinv->coeffs, Q->coeffs, Qlen, n);
+    }
+    else
+    {
+        fmpz_poly_t t;
+        fmpz_poly_init2(t, n);
+        _fmpz_poly_inv_series_basecase(t->coeffs, Q->coeffs, Qlen, n);
+        fmpz_poly_swap(Qinv, t);
+        fmpz_poly_clear(t);
+    }
+
+    _fmpz_poly_set_length(Qinv, n);
+    _fmpz_poly_normalise(Qinv);
+}
+