add arb_poly_log1p_series, acb_poly_log1p_series

acb_poly.h

@@ -568,9 +568,11 @@ void _acb_poly_sqrt_series(acb_ptr g, acb_srcptr h, slong hlen, slong len, slong
 void acb_poly_sqrt_series(acb_poly_t g, const acb_poly_t h, slong n, slong prec);
 
 void _acb_poly_log_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec);
-
 void acb_poly_log_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec);
 
+void _acb_poly_log1p_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec);
+void acb_poly_log1p_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec);
+
 void _acb_poly_atan_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec);
 
 void acb_poly_atan_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec);

acb_poly/log1p_series.c

@@ -0,0 +1,94 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_poly.h"
+
+void
+_acb_poly_log1p_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec)
+{
+    acb_t a;
+
+    flen = FLINT_MIN(flen, n);
+
+    acb_init(a);
+    acb_log1p(a, f, prec);
+
+    if (flen == 1)
+    {
+        _acb_vec_zero(res + 1, n - 1);
+    }
+    else if (n == 2)
+    {
+        acb_add_ui(res, f + 0, 1, prec);
+        acb_div(res + 1, f + 1, res + 0, prec);
+    }
+    else if (_acb_vec_is_zero(f + 1, flen - 2))  /* f = a + bx^d */
+    {
+        slong i, j, d = flen - 1;
+
+        acb_add_ui(res, f + 0, 1, prec);
+
+        for (i = 1, j = d; j < n; j += d, i++)
+        {
+            if (i == 1)
+                acb_div(res + j, f + d, res, prec);
+            else
+                acb_mul(res + j, res + j - d, res + d, prec);
+            _acb_vec_zero(res + j - d + 1, flen - 2);
+        }
+        _acb_vec_zero(res + j - d + 1, n - (j - d + 1));
+
+        for (i = 2, j = 2 * d; j < n; j += d, i++)
+            acb_div_si(res + j, res + j, i % 2 ? i : -i, prec);
+    }
+    else
+    {
+        acb_ptr f_diff, f_inv;
+        slong alloc;
+
+        alloc = n + flen;
+        f_inv = _acb_vec_init(alloc);
+        f_diff = f_inv + n;
+
+        acb_add_ui(f_diff, f, 1, prec);
+        _acb_vec_set(f_diff + 1, f + 1, flen - 1);
+        _acb_poly_inv_series(f_inv, f_diff, flen, n, prec);
+        _acb_poly_derivative(f_diff, f, flen, prec);
+        _acb_poly_mullow(res, f_inv, n - 1, f_diff, flen - 1, n - 1, prec);
+        _acb_poly_integral(res, res, n, prec);
+
+        _acb_vec_clear(f_inv, alloc);
+    }
+
+    acb_swap(res, a);
+    acb_clear(a);
+}
+
+void
+acb_poly_log1p_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec)
+{
+    slong flen = f->length;
+
+    if (flen == 0 || n == 0)
+    {
+        acb_poly_zero(res);
+        return;
+    }
+
+    if (flen == 1 /*&& !acb_contains_si(f->coeffs, -1)*/)
+        n = 1;
+
+    acb_poly_fit_length(res, n);
+    _acb_poly_log1p_series(res->coeffs, f->coeffs, flen, n, prec);
+    _acb_poly_set_length(res, n);
+    _acb_poly_normalise(res);
+}
+

acb_poly/test/t-log1p_series.c

@@ -0,0 +1,75 @@
+/*
+    Copyright (C) 2012 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("log1p_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* compare with log_series */
+    for (iter = 0; iter < 5000 * arb_test_multiplier(); iter++)
+    {
+        acb_poly_t a, b, c;
+        slong m, n, prec;
+
+        prec = 2 + n_randint(state, 200);
+        m = n_randint(state, 10);
+        n = n_randint(state, 10);
+
+        acb_poly_init(a);
+        acb_poly_init(b);
+        acb_poly_init(c);
+
+        acb_poly_randtest(a, state, 1 + n_randint(state, 10), prec, 5);
+        acb_poly_randtest(b, state, 1 + n_randint(state, 10), prec, 5);
+        acb_poly_randtest(c, state, 1 + n_randint(state, 10), prec, 5);
+
+        if (n_randint(state, 2))
+            acb_poly_log1p_series(b, a, n, prec);
+        else
+        {
+            acb_poly_set(b, a);
+            acb_poly_log1p_series(b, b, n, prec);
+        }
+
+        acb_poly_add_si(c, a, 1, prec);
+        acb_poly_log_series(c, c, m, prec);
+
+        acb_poly_truncate(b, FLINT_MIN(m, n));
+        acb_poly_truncate(c, FLINT_MIN(m, n));
+
+        if (!acb_poly_overlaps(b, c))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("m = %wd, n = %wd\n\n", m, n);
+            flint_printf("a = "); acb_poly_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = "); acb_poly_printd(b, 15); flint_printf("\n\n");
+            flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_poly_clear(a);
+        acb_poly_clear(b);
+        acb_poly_clear(c);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb_poly.h

@@ -544,9 +544,11 @@ void _arb_poly_sqrt_series(arb_ptr g,
 void arb_poly_sqrt_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec);
 
 void _arb_poly_log_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec);
-
 void arb_poly_log_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec);
 
+void _arb_poly_log1p_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec);
+void arb_poly_log1p_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec);
+
 void _arb_poly_atan_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec);
 
 void arb_poly_atan_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec);

arb_poly/log1p_series.c

@@ -0,0 +1,94 @@
+/*
+    Copyright (C) 2017 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_poly.h"
+
+void
+_arb_poly_log1p_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)
+{
+    arb_t a;
+
+    flen = FLINT_MIN(flen, n);
+
+    arb_init(a);
+    arb_log1p(a, f, prec);
+
+    if (flen == 1)
+    {
+        _arb_vec_zero(res + 1, n - 1);
+    }
+    else if (n == 2)
+    {
+        arb_add_ui(res, f + 0, 1, prec);
+        arb_div(res + 1, f + 1, res + 0, prec);
+    }
+    else if (_arb_vec_is_zero(f + 1, flen - 2))  /* f = a + bx^d */
+    {
+        slong i, j, d = flen - 1;
+
+        arb_add_ui(res, f + 0, 1, prec);
+
+        for (i = 1, j = d; j < n; j += d, i++)
+        {
+            if (i == 1)
+                arb_div(res + j, f + d, res, prec);
+            else
+                arb_mul(res + j, res + j - d, res + d, prec);
+            _arb_vec_zero(res + j - d + 1, flen - 2);
+        }
+        _arb_vec_zero(res + j - d + 1, n - (j - d + 1));
+
+        for (i = 2, j = 2 * d; j < n; j += d, i++)
+            arb_div_si(res + j, res + j, i % 2 ? i : -i, prec);
+    }
+    else
+    {
+        arb_ptr f_diff, f_inv;
+        slong alloc;
+
+        alloc = n + flen;
+        f_inv = _arb_vec_init(alloc);
+        f_diff = f_inv + n;
+
+        arb_add_ui(f_diff, f, 1, prec);
+        _arb_vec_set(f_diff + 1, f + 1, flen - 1);
+        _arb_poly_inv_series(f_inv, f_diff, flen, n, prec);
+        _arb_poly_derivative(f_diff, f, flen, prec);
+        _arb_poly_mullow(res, f_inv, n - 1, f_diff, flen - 1, n - 1, prec);
+        _arb_poly_integral(res, res, n, prec);
+
+        _arb_vec_clear(f_inv, alloc);
+    }
+
+    arb_swap(res, a);
+    arb_clear(a);
+}
+
+void
+arb_poly_log1p_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)
+{
+    slong flen = f->length;
+
+    if (flen == 0 || n == 0)
+    {
+        arb_poly_zero(res);
+        return;
+    }
+
+    if (flen == 1 /*&& !arb_contains_si(f->coeffs, -1)*/)
+        n = 1;
+
+    arb_poly_fit_length(res, n);
+    _arb_poly_log1p_series(res->coeffs, f->coeffs, flen, n, prec);
+    _arb_poly_set_length(res, n);
+    _arb_poly_normalise(res);
+}
+

arb_poly/test/t-log1p_series.c

@@ -0,0 +1,75 @@
+/*
+    Copyright (C) 2012 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "arb_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("log1p_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* compare with log_series */
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+    {
+        arb_poly_t a, b, c;
+        slong m, n, prec;
+
+        prec = 2 + n_randint(state, 200);
+        m = n_randint(state, 10);
+        n = n_randint(state, 10);
+
+        arb_poly_init(a);
+        arb_poly_init(b);
+        arb_poly_init(c);
+
+        arb_poly_randtest(a, state, 1 + n_randint(state, 10), prec, 5);
+        arb_poly_randtest(b, state, 1 + n_randint(state, 10), prec, 5);
+        arb_poly_randtest(c, state, 1 + n_randint(state, 10), prec, 5);
+
+        if (n_randint(state, 2))
+            arb_poly_log1p_series(b, a, n, prec);
+        else
+        {
+            arb_poly_set(b, a);
+            arb_poly_log1p_series(b, b, n, prec);
+        }
+
+        arb_poly_add_si(c, a, 1, prec);
+        arb_poly_log_series(c, c, m, prec);
+
+        arb_poly_truncate(b, FLINT_MIN(m, n));
+        arb_poly_truncate(c, FLINT_MIN(m, n));
+
+        if (!arb_poly_overlaps(b, c))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("m = %wd, n = %wd\n\n", m, n);
+            flint_printf("a = "); arb_poly_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = "); arb_poly_printd(b, 15); flint_printf("\n\n");
+            flint_printf("c = "); arb_poly_printd(c, 15); flint_printf("\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        arb_poly_clear(b);
+        arb_poly_clear(c);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/acb_poly.rst

@@ -746,6 +746,12 @@ Elementary functions
     The underscore method supports aliasing of the input and output
     arrays. It requires that *flen* and *n* are greater than zero.
 
+.. function:: void _acb_poly_log1p_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec)
+
+.. function:: void acb_poly_log1p_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec)
+
+    Computes the power series `\log(1+f)`, with better accuracy when the constant term of *f* is small.
+
 .. function:: void _acb_poly_atan_series(acb_ptr res, acb_srcptr f, slong flen, slong n, slong prec)
 
 .. function:: void acb_poly_atan_series(acb_poly_t res, const acb_poly_t f, slong n, slong prec)

doc/source/arb_poly.rst

@@ -791,6 +791,12 @@ Powers and elementary functions
     The underscore method supports aliasing of the input and output
     arrays. It requires that *flen* and *n* are greater than zero.
 
+.. function:: void _arb_poly_log1p_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)
+
+.. function:: void arb_poly_log1p_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)
+
+    Computes the power series `\log(1+f)`, with better accuracy when the constant term of *f* is small.
+
 .. function:: void _arb_poly_atan_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)
 
 .. function:: void arb_poly_atan_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)