improve rising factorial series for length 1, 2; fix a bug and add missing test file

fmpcb_poly/rising_ui_series.c

@@ -19,22 +19,12 @@
 =============================================================================*/
 /******************************************************************************
 
-    Copyright (C) 2012 Fredrik Johansson
+    Copyright (C) 2012, 2013 Fredrik Johansson
 
 ******************************************************************************/
 
 #include "fmpcb_poly.h"
-
+#include "gamma.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
 _fmpcb_poly_rising_ui_series_bsplit(fmpcb_ptr res,
@@ -55,8 +45,8 @@ _fmpcb_poly_rising_ui_series_bsplit(fmpcb_ptr res,
 
         long m = a + (b - a) / 2;
 
-        len1 = length(flen, m - a, trunc);
+        len1 = poly_pow_length(flen, m - a, trunc);
-        len2 = length(flen, b - m, trunc);
+        len2 = poly_pow_length(flen, b - m, trunc);
 
         L = _fmpcb_vec_init(len1 + len2);
         R = L + len1;
@@ -76,10 +66,22 @@ _fmpcb_poly_rising_ui_series(fmpcb_ptr res,
     fmpcb_srcptr f, long flen, ulong r,
         long trunc, long prec)
 {
-    _fmpcb_poly_rising_ui_series_bsplit(res, f, flen, 0, r, trunc, prec);
+    if (trunc == 1 || flen == 1)
+    {
+        gamma_rising_fmpcb_ui_bsplit(res, f, r, prec);
+        _fmpcb_vec_zero(res + 1, trunc - 1);
+    }
+    else if (trunc == 2)
+    {
+        gamma_rising2_fmpcb_ui(res, res + 1, f, r, prec);
+        fmpcb_mul(res + 1, res + 1, f + 1, prec);
+    }
+    else
+    {
+        _fmpcb_poly_rising_ui_series_bsplit(res, f, flen, 0, r, trunc, prec);
+    }
 }
 
-
 void
 fmpcb_poly_rising_ui_series(fmpcb_poly_t res, const fmpcb_poly_t f, ulong r, long trunc, long prec)
 {
@@ -97,20 +99,20 @@ fmpcb_poly_rising_ui_series(fmpcb_poly_t res, const fmpcb_poly_t f, ulong r, lon
         return;
     }
 
-    len = length(f->length, r, trunc);
+    len = poly_pow_length(f->length, r, trunc);
 
     if (f == res)
     {
         fmpcb_poly_t tmp;
         fmpcb_poly_init(tmp);
-        fmpcb_poly_rising_ui_series(tmp, f, r, trunc, prec);
+        fmpcb_poly_rising_ui_series(tmp, f, r, len, prec);
         fmpcb_poly_swap(tmp, res);
         fmpcb_poly_clear(tmp);
     }
     else
     {
         fmpcb_poly_fit_length(res, len);
-        _fmpcb_poly_rising_ui_series(res->coeffs, f->coeffs, f->length, r, trunc, prec);
+        _fmpcb_poly_rising_ui_series(res->coeffs, f->coeffs, f->length, r, len, prec);
         _fmpcb_poly_set_length(res, len);
         _fmpcb_poly_normalise(res);
     }

fmpcb_poly/test/t-rising_ui_series.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 "fmpcb_poly.h"
+
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("rising_ui_series....");
+    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;
+        fmpcb_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);
+
+        fmpcb_poly_init(f);
+        fmpcb_poly_init(g);
+        fmpcb_poly_init(h1);
+        fmpcb_poly_init(h2);
+        fmpcb_poly_init(h1h2);
+        fmpcb_poly_init(h3);
+
+        fmpcb_poly_randtest(f, state, 1 + n_randint(state, 20), bits, 4);
+        fmpcb_poly_set(g, f);
+
+        /* g = f + 1 */
+        if (g->length == 0)
+        {
+            fmpcb_poly_fit_length(g, 1);
+            fmpcb_set_ui(g->coeffs, a);
+            _fmpcb_poly_set_length(g, 1);
+            _fmpcb_poly_normalise(g);
+        }
+        else
+        {
+            fmpcb_add_ui(g->coeffs, g->coeffs, a, bits);
+            _fmpcb_poly_normalise(g);
+        }
+
+        fmpcb_poly_rising_ui_series(h1, f, a, trunc, bits);
+        fmpcb_poly_rising_ui_series(h2, g, b, trunc, bits);
+        fmpcb_poly_rising_ui_series(h3, f, a + b, trunc, bits);
+
+        fmpcb_poly_mullow(h1h2, h1, h2, trunc, bits);
+
+        if (!fmpcb_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 = "); fmpcb_poly_printd(f, 15); printf("\n\n");
+            printf("g = "); fmpcb_poly_printd(g, 15); printf("\n\n");
+            printf("h1 = "); fmpcb_poly_printd(h1, 15); printf("\n\n");
+            printf("h2 = "); fmpcb_poly_printd(h2, 15); printf("\n\n");
+            printf("h1h2 = "); fmpcb_poly_printd(h1h2, 15); printf("\n\n");
+            printf("h3 = "); fmpcb_poly_printd(h3, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmpcb_poly_rising_ui_series(f, f, a, trunc, bits);
+
+        if (!fmpcb_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 = "); fmpcb_poly_printd(f, 15); printf("\n\n");
+            printf("h1 = "); fmpcb_poly_printd(h1, 15); printf("\n\n");
+
+            abort();
+        }
+
+        fmpcb_poly_clear(f);
+        fmpcb_poly_clear(g);
+        fmpcb_poly_clear(h1);
+        fmpcb_poly_clear(h2);
+        fmpcb_poly_clear(h1h2);
+        fmpcb_poly_clear(h3);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

fmprb_poly/rising_ui_series.c

@@ -19,11 +19,12 @@
 =============================================================================*/
 /******************************************************************************
 
-    Copyright (C) 2012 Fredrik Johansson
+    Copyright (C) 2012, 2013 Fredrik Johansson
 
 ******************************************************************************/
 
 #include "fmprb_poly.h"
+#include "gamma.h"
 
 static void
 _fmprb_poly_rising_ui_series_bsplit(fmprb_ptr res,
@@ -65,7 +66,20 @@ _fmprb_poly_rising_ui_series(fmprb_ptr res,
     fmprb_srcptr f, long flen, ulong r,
         long trunc, long prec)
 {
-    _fmprb_poly_rising_ui_series_bsplit(res, f, flen, 0, r, trunc, prec);
+    if (trunc == 1 || flen == 1)
+    {
+        gamma_rising_fmprb_ui_bsplit(res, f, r, prec);
+        _fmprb_vec_zero(res + 1, trunc - 1);
+    }
+    else if (trunc == 2)
+    {
+        gamma_rising2_fmprb_ui(res, res + 1, f, r, prec);
+        fmprb_mul(res + 1, res + 1, f + 1, prec);
+    }
+    else
+    {
+        _fmprb_poly_rising_ui_series_bsplit(res, f, flen, 0, r, trunc, prec);
+    }
 }
 
 void
@@ -91,14 +105,14 @@ fmprb_poly_rising_ui_series(fmprb_poly_t res, const fmprb_poly_t f, ulong r, lon
     {
         fmprb_poly_t tmp;
         fmprb_poly_init(tmp);
-        fmprb_poly_rising_ui_series(tmp, f, r, trunc, prec);
+        fmprb_poly_rising_ui_series(tmp, f, r, len, prec);
         fmprb_poly_swap(tmp, res);
         fmprb_poly_clear(tmp);
     }
     else
     {
         fmprb_poly_fit_length(res, len);
-        _fmprb_poly_rising_ui_series(res->coeffs, f->coeffs, f->length, r, trunc, prec);
+        _fmprb_poly_rising_ui_series(res->coeffs, f->coeffs, f->length, r, len, prec);
         _fmprb_poly_set_length(res, len);
         _fmprb_poly_normalise(res);
     }