add fmprb_poly rising factorial; plus helpers overlaps and randtest

doc/source/fmpcb_poly.rst

@@ -93,8 +93,7 @@ Comparisons
 
     Returns nonzero iff *poly2* is contained in *poly1*.
 
-.. function:: int _fmpcb_poly_overlaps(const fmpcb_struct * poly1, long len1,
-        const fmpcb_struct * poly2, long len2)
+.. function:: int _fmpcb_poly_overlaps(const fmpcb_struct * poly1, long len1, const fmpcb_struct * poly2, long len2)
 
 .. function:: int fmpcb_poly_overlaps(const fmpcb_poly_t poly1, const fmpcb_poly_t poly2)
 

doc/source/fmprb_poly.rst

@@ -78,6 +78,14 @@ Input and output
     coefficient using *fmprb_printd*.
 
 
+Random generation
+-------------------------------------------------------------------------------
+
+.. function:: void fmprb_poly_randtest(fmprb_poly_t poly, flint_rand_t state, long len, long prec, long mag_bits)
+
+    Creates a random polynomial with length at most *len*.
+
+
 Comparisons
 -------------------------------------------------------------------------------
 
@@ -90,6 +98,12 @@ Comparisons
     Returns nonzero iff *A* and *B* are equal as polynomial balls, i.e. all
     coefficients have equal midpoint and radius.
 
+.. function:: int _fmprb_poly_overlaps(const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2)
+
+.. function:: int fmprb_poly_overlaps(const fmprb_poly_t poly1, const fmprb_poly_t poly2)
+
+    Returns nonzero iff *poly1* overlaps with *poly2*. The underscore
+    function requires that *len1* ist at least as large as *len2*.
 
 Arithmetic
 -------------------------------------------------------------------------------
@@ -341,4 +355,12 @@ Special functions
     Sets `f` to the series expansion of `\log(\Gamma(1-x))`, truncated to
     length `n`.
 
+.. function:: void _fmprb_poly_rfac_series_ui(fmprb_struct * res, const fmprb_struct * f, long flen, ulong r, long trunc, long prec)
+
+.. function:: void fmprb_poly_rfac_series_ui(fmprb_poly_t res, const fmprb_poly_t f, ulong r, long trunc, long prec)
+
+    Sets *res* to the rising factorial `(f) (f+1) (f+2) \cdots (f+r-1)`, truncated
+    to length *trunc*. The underscore method assumes that *flen*, *r* and *trunc*
+    are at least 1, and does not support aliasing. Uses binary splitting.
+
 

doc/source/index.rst

@@ -16,7 +16,7 @@ General information
 ::::::::::::::::::::
 
 .. toctree::
-   :maxdepth: 1
+   :maxdepth: 2
 
    overview.rst
    setup.rst
@@ -26,7 +26,7 @@ Module documentation
 ::::::::::::::::::::
 
 .. toctree::
-   :maxdepth: 1
+   :maxdepth: 2
 
    fmpr.rst
    fmprb.rst
@@ -41,7 +41,7 @@ Credits and references
 ::::::::::::::::::::::::
 
 .. toctree::
-   :maxdepth: 1
+   :maxdepth: 2
 
    credits.rst
 

fmprb_poly.h

@@ -107,10 +107,18 @@ int fmprb_poly_contains_fmpq_poly(const fmprb_poly_t poly1, const fmpq_poly_t po
 
 int fmprb_poly_equal(const fmprb_poly_t A, const fmprb_poly_t B);
 
+int _fmprb_poly_overlaps(const fmprb_struct * poly1, long len1, const fmprb_struct * poly2, long len2);
+
+int fmprb_poly_overlaps(const fmprb_poly_t poly1, const fmprb_poly_t poly2);
+
 /* IO */
 
 void fmprb_poly_printd(const fmprb_poly_t poly, long digits);
 
+/* Random generation */
+
+void fmprb_poly_randtest(fmprb_poly_t poly, flint_rand_t state, long len, long prec, long mag_bits);
+
 /* Arithmetic */
 
 void
@@ -314,5 +322,9 @@ void fmprb_poly_exp_series(fmprb_poly_t f, const fmprb_poly_t h, long n, long pr
 
 void fmprb_poly_log_gamma_series(fmprb_poly_t z, long n, long prec);
 
+void _fmprb_poly_rfac_series_ui(fmprb_struct * res, const fmprb_struct * f, long flen, ulong r, long trunc, long prec);
+
+void fmprb_poly_rfac_series_ui(fmprb_poly_t res, const fmprb_poly_t f, ulong r, long trunc, long prec);
+
 #endif
 

fmprb_poly/overlaps.c

@@ -0,0 +1,56 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+int
+_fmprb_poly_overlaps(const fmprb_struct * poly1, long len1,
+        const fmprb_struct * poly2, long len2)
+{
+    long i;
+
+    for (i = 0; i < len2; i++)
+        if (!fmprb_overlaps(poly1 + i, poly2 + i))
+            return 0;
+
+    for (i = len2; i < len1; i++)
+        if (!fmprb_contains_zero(poly1 + i))
+            return 0;
+
+    return 1;
+}
+
+int
+fmprb_poly_overlaps(const fmprb_poly_t poly1, const fmprb_poly_t poly2)
+{
+    long len1 = poly1->length;
+    long len2 = poly2->length;
+
+    if (len1 >= len2)
+        return _fmprb_poly_overlaps(poly1->coeffs, len1, poly2->coeffs, len2);
+    else
+        return _fmprb_poly_overlaps(poly2->coeffs, len2, poly1->coeffs, len1);
+}
+

fmprb_poly/randtest.c

@@ -0,0 +1,49 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+static void
+fmprb_randtest2(fmprb_t x, flint_rand_t state, long prec, long mag_bits)
+{
+    fmpr_randtest(fmprb_midref(x), state, prec, mag_bits);
+
+    fmpr_randtest_not_zero(fmprb_radref(x), state, FMPRB_RAD_PREC, 4);
+    fmpr_abs(fmprb_radref(x), fmprb_radref(x));
+    fmpz_add(fmpr_expref(fmprb_radref(x)), fmpr_expref(fmprb_radref(x)), fmpr_expref(fmprb_midref(x)));
+}
+
+void
+fmprb_poly_randtest(fmprb_poly_t poly, flint_rand_t state, long len, long prec, long mag_bits)
+{
+    long i;
+
+    fmprb_poly_fit_length(poly, len);
+    for (i = 0; i < len; i++)
+        fmprb_randtest2(poly->coeffs + i, state, prec, mag_bits);
+    _fmprb_poly_set_length(poly, len);
+    _fmprb_poly_normalise(poly);
+}
+

fmprb_poly/rfac_series_ui.c

@@ -0,0 +1,118 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+
+static __inline__ long length(long flen, ulong r, long trunc)
+{
+    mp_limb_t hi, lo;
+    umul_ppmm(hi, lo, flen - 1, r);
+    add_ssaaaa(hi, lo, hi, lo, 0, 1);
+    if (hi != 0 || lo > (mp_limb_t) LONG_MAX)
+        return trunc;
+    return FLINT_MIN(lo, trunc);
+}
+
+static void
+_fmprb_poly_rfac_series_ui_bsplit(fmprb_struct * res,
+    const fmprb_struct * f, long flen, ulong a, ulong b,
+        long trunc, long prec)
+{
+    flen = FLINT_MIN(flen, trunc);
+
+    if (b - a == 1)
+    {
+        fmprb_add_ui(res, f, a, prec);
+        _fmprb_vec_set(res + 1, f + 1, flen - 1);
+    }
+    else
+    {
+        fmprb_struct *L, *R;
+        long len1, len2;
+
+        long m = a + (b - a) / 2;
+
+        len1 = length(flen, m - a, trunc);
+        len2 = length(flen, b - m, trunc);
+
+        L = _fmprb_vec_init(len1 + len2);
+        R = L + len1;
+
+        _fmprb_poly_rfac_series_ui_bsplit(L, f, flen, a, m, trunc, prec);
+        _fmprb_poly_rfac_series_ui_bsplit(R, f, flen, m, b, trunc, prec);
+
+        _fmprb_poly_mullow(res, L, len1, R, len2,
+            FLINT_MIN(trunc, len1 + len2 - 1), prec);
+
+        _fmprb_vec_clear(L, len1 + len2);
+    }
+}
+
+void
+_fmprb_poly_rfac_series_ui(fmprb_struct * res,
+    const fmprb_struct * f, long flen, ulong r,
+        long trunc, long prec)
+{
+    _fmprb_poly_rfac_series_ui_bsplit(res, f, flen, 0, r, trunc, prec);
+}
+
+
+void
+fmprb_poly_rfac_series_ui(fmprb_poly_t res, const fmprb_poly_t f, ulong r, long trunc, long prec)
+{
+    long len;
+
+    if (f->length == 0 && r != 0)
+    {
+        fmprb_poly_zero(res);
+        return;
+    }
+
+    if (r == 0)
+    {
+        fmprb_poly_one(res);
+        return;
+    }
+
+    len = length(f->length, r, trunc);
+
+    if (f == res)
+    {
+        fmprb_poly_t tmp;
+        fmprb_poly_init(tmp);
+        fmprb_poly_rfac_series_ui(tmp, f, r, trunc, prec);
+        fmprb_poly_swap(tmp, res);
+        fmprb_poly_clear(tmp);
+    }
+    else
+    {
+        fmprb_poly_fit_length(res, len);
+        _fmprb_poly_rfac_series_ui(res->coeffs, f->coeffs, f->length, r, trunc, prec);
+        _fmprb_poly_set_length(res, len);
+        _fmprb_poly_normalise(res);
+    }
+}
+

fmprb_poly/test/t-rfac_series_ui.c

@@ -0,0 +1,128 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("rfac_series_ui....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* check rf(f, a) * rf(f + a, b) = rf(f, a + b) */
+    for (iter = 0; iter < 1000; iter++)
+    {
+        long bits, trunc;
+        ulong a, b;
+        fmprb_poly_t f, g, h1, h2, h1h2, h3;
+
+        bits = 2 + n_randint(state, 200);
+        trunc = 1 + n_randint(state, 20);
+        a = n_randint(state, 10);
+        b = n_randint(state, 10);
+
+        fmprb_poly_init(f);
+        fmprb_poly_init(g);
+        fmprb_poly_init(h1);
+        fmprb_poly_init(h2);
+        fmprb_poly_init(h1h2);
+        fmprb_poly_init(h3);
+
+        fmprb_poly_randtest(f, state, 1 + n_randint(state, 20), bits, 4);
+        fmprb_poly_set(g, f);
+
+        /* g = f + 1 */
+        if (g->length == 0)
+        {
+            fmprb_poly_fit_length(g, 1);
+            fmprb_set_ui(g->coeffs, a);
+            _fmprb_poly_set_length(g, 1);
+            _fmprb_poly_normalise(g);
+        }
+        else
+        {
+            fmprb_add_ui(g->coeffs, g->coeffs, a, bits);
+            _fmprb_poly_normalise(g);
+        }
+
+        fmprb_poly_rfac_series_ui(h1, f, a, trunc, bits);
+        fmprb_poly_rfac_series_ui(h2, g, b, trunc, bits);
+        fmprb_poly_rfac_series_ui(h3, f, a + b, trunc, bits);
+
+        fmprb_poly_mullow(h1h2, h1, h2, trunc, bits);
+
+        if (!fmprb_poly_overlaps(h1h2, h3))
+        {
+            printf("FAIL\n\n");
+            printf("bits = %ld\n", bits);
+            printf("trunc = %ld\n", trunc);
+            printf("a = %lu\n", a);
+            printf("b = %lu\n", a);
+
+            printf("f = "); fmprb_poly_printd(f, 15); printf("\n\n");
+            printf("g = "); fmprb_poly_printd(g, 15); printf("\n\n");
+            printf("h1 = "); fmprb_poly_printd(h1, 15); printf("\n\n");
+            printf("h2 = "); fmprb_poly_printd(h2, 15); printf("\n\n");
+            printf("h1h2 = "); fmprb_poly_printd(h1h2, 15); printf("\n\n");
+            printf("h3 = "); fmprb_poly_printd(h3, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmprb_poly_rfac_series_ui(f, f, a, trunc, bits);
+
+        if (!fmprb_poly_equal(f, h1))
+        {
+            printf("FAIL (aliasing)\n\n");
+
+            printf("bits = %ld\n", bits);
+            printf("trunc = %ld\n", trunc);
+            printf("a = %lu\n", a);
+            printf("b = %lu\n", a);
+
+            printf("f = "); fmprb_poly_printd(f, 15); printf("\n\n");
+            printf("h1 = "); fmprb_poly_printd(h1, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmprb_poly_clear(f);
+        fmprb_poly_clear(g);
+        fmprb_poly_clear(h1);
+        fmprb_poly_clear(h2);
+        fmprb_poly_clear(h1h2);
+        fmprb_poly_clear(h3);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}