add acb_poly series versions of l, hardy_theta and hardy_z

acb_dirichlet.h

@@ -21,6 +21,7 @@
 #endif
 
 #include "acb.h"
+#include "acb_poly.h"
 #include "dirichlet.h"
 
 #ifdef __cplusplus
@@ -131,6 +132,14 @@ void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const
 
 void acb_dirichlet_l_jet(acb_ptr res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec);
 
+void _acb_dirichlet_l_series(acb_ptr res, acb_srcptr s, slong slen,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    int deflate, slong len, slong prec);
+
+void acb_dirichlet_l_series(acb_poly_t res, const acb_poly_t s,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    int deflate, slong len, slong prec);
+
 void acb_dirichlet_hardy_theta(acb_ptr res, const acb_t t,
     const dirichlet_group_t G, const dirichlet_char_t chi,
     slong len, slong prec);
@@ -139,6 +148,10 @@ void acb_dirichlet_hardy_z(acb_ptr res, const acb_t t,
     const dirichlet_group_t G, const dirichlet_char_t chi,
     slong len, slong prec);
 
+void _acb_dirichlet_hardy_theta_series(acb_ptr res, acb_srcptr s, slong slen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec);
+void acb_dirichlet_hardy_theta_series(acb_poly_t res, const acb_poly_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec);
+void _acb_dirichlet_hardy_z_series(acb_ptr res, acb_srcptr s, slong slen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec);
+void acb_dirichlet_hardy_z_series(acb_poly_t res, const acb_poly_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec);
 
 /* utils */
 

acb_dirichlet/hardy_theta_series.c

@@ -0,0 +1,75 @@
+/*
+    Copyright (C) 2016 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_dirichlet.h"
+
+void
+_acb_dirichlet_hardy_theta_series(acb_ptr res, acb_srcptr s, slong slen,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    slong len, slong prec)
+{
+    slen = FLINT_MIN(slen, len);
+
+    if (len == 0)
+        return;
+
+    if (slen == 1)
+    {
+        acb_dirichlet_hardy_theta(res, s, G, chi, 1, prec);
+        _acb_vec_zero(res + 1, len - 1);
+    }
+    else
+    {
+        acb_ptr t, u;
+        t = _acb_vec_init(len);
+        u = _acb_vec_init(slen);
+
+        acb_dirichlet_hardy_theta(t, s, G, chi, len, prec);
+
+        /* compose with nonconstant part */
+        acb_zero(u);
+        _acb_vec_set(u + 1, s + 1, slen - 1);
+        _acb_poly_compose_series(res, t, len, u, slen, len, prec);
+
+        _acb_vec_clear(t, len);
+        _acb_vec_clear(u, slen);
+    }
+}
+
+void
+acb_dirichlet_hardy_theta_series(acb_poly_t res, const acb_poly_t s,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    slong len, slong prec)
+{
+    if (len == 0)
+    {
+        acb_poly_zero(res);
+        return;
+    }
+
+    acb_poly_fit_length(res, len);
+
+    if (s->length == 0)
+    {
+        acb_t t;
+        acb_init(t);
+        _acb_dirichlet_hardy_theta_series(res->coeffs, t, 1, G, chi, len, prec);
+        acb_clear(t);
+    }
+    else
+    {
+        _acb_dirichlet_hardy_theta_series(res->coeffs, s->coeffs, s->length, G, chi, len, prec);
+    }
+
+    _acb_poly_set_length(res, len);
+    _acb_poly_normalise(res);
+}
+

acb_dirichlet/hardy_z_series.c

@@ -0,0 +1,75 @@
+/*
+    Copyright (C) 2016 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_dirichlet.h"
+
+void
+_acb_dirichlet_hardy_z_series(acb_ptr res, acb_srcptr s, slong slen,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    slong len, slong prec)
+{
+    slen = FLINT_MIN(slen, len);
+
+    if (len == 0)
+        return;
+
+    if (slen == 1)
+    {
+        acb_dirichlet_hardy_z(res, s, G, chi, 1, prec);
+        _acb_vec_zero(res + 1, len - 1);
+    }
+    else
+    {
+        acb_ptr t, u;
+        t = _acb_vec_init(len);
+        u = _acb_vec_init(slen);
+
+        acb_dirichlet_hardy_z(t, s, G, chi, len, prec);
+
+        /* compose with nonconstant part */
+        acb_zero(u);
+        _acb_vec_set(u + 1, s + 1, slen - 1);
+        _acb_poly_compose_series(res, t, len, u, slen, len, prec);
+
+        _acb_vec_clear(t, len);
+        _acb_vec_clear(u, slen);
+    }
+}
+
+void
+acb_dirichlet_hardy_z_series(acb_poly_t res, const acb_poly_t s,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    slong len, slong prec)
+{
+    if (len == 0)
+    {
+        acb_poly_zero(res);
+        return;
+    }
+
+    acb_poly_fit_length(res, len);
+
+    if (s->length == 0)
+    {
+        acb_t t;
+        acb_init(t);
+        _acb_dirichlet_hardy_z_series(res->coeffs, t, 1, G, chi, len, prec);
+        acb_clear(t);
+    }
+    else
+    {
+        _acb_dirichlet_hardy_z_series(res->coeffs, s->coeffs, s->length, G, chi, len, prec);
+    }
+
+    _acb_poly_set_length(res, len);
+    _acb_poly_normalise(res);
+}
+

acb_dirichlet/l_jet.c

@@ -27,6 +27,7 @@ acb_dirichlet_l_jet(acb_ptr res, const acb_t s,
     if (len <= 0)
         return;
 
+    /* todo: reflection formula */
     if (G->q == 1)
     {
         acb_init(a);

acb_dirichlet/l_series.c

@@ -0,0 +1,75 @@
+/*
+    Copyright (C) 2016 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_dirichlet.h"
+
+void
+_acb_dirichlet_l_series(acb_ptr res, acb_srcptr s, slong slen,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    int deflate, slong len, slong prec)
+{
+    slen = FLINT_MIN(slen, len);
+
+    if (len == 0)
+        return;
+
+    if (slen == 1 && !deflate)
+    {
+        acb_dirichlet_l(res, s, G, chi, prec);
+        _acb_vec_zero(res + 1, len - 1);
+    }
+    else
+    {
+        acb_ptr t, u;
+        t = _acb_vec_init(len);
+        u = _acb_vec_init(slen);
+
+        acb_dirichlet_l_jet(t, s, G, chi, deflate, len, prec);
+
+        /* compose with nonconstant part */
+        acb_zero(u);
+        _acb_vec_set(u + 1, s + 1, slen - 1);
+        _acb_poly_compose_series(res, t, len, u, slen, len, prec);
+
+        _acb_vec_clear(t, len);
+        _acb_vec_clear(u, slen);
+    }
+}
+
+void
+acb_dirichlet_l_series(acb_poly_t res, const acb_poly_t s,
+    const dirichlet_group_t G, const dirichlet_char_t chi,
+    int deflate, slong len, slong prec)
+{
+    if (len == 0)
+    {
+        acb_poly_zero(res);
+        return;
+    }
+
+    acb_poly_fit_length(res, len);
+
+    if (s->length == 0)
+    {
+        acb_t t;
+        acb_init(t);
+        _acb_dirichlet_l_series(res->coeffs, t, 1, G, chi, deflate, len, prec);
+        acb_clear(t);
+    }
+    else
+    {
+        _acb_dirichlet_l_series(res->coeffs, s->coeffs, s->length, G, chi, deflate, len, prec);
+    }
+
+    _acb_poly_set_length(res, len);
+    _acb_poly_normalise(res);
+}
+

acb_dirichlet/test/t-hardy_theta_series.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2013 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_dirichlet.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("hardy_theta_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
+    {
+        slong m, n1, n2, bits1, bits2, bits3;
+        acb_poly_t S, A, B, C;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+        ulong q;
+
+        bits1 = 2 + n_randint(state, 200);
+        bits2 = 2 + n_randint(state, 200);
+        bits3 = 2 + n_randint(state, 200);
+
+        m = 1 + n_randint(state, 8);
+        n1 = 1 + n_randint(state, 8);
+        n2 = 1 + n_randint(state, 8);
+
+        do {
+            q = 1 + n_randint(state, 15);
+        } while (q % 4 == 2);
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
+        do {
+            dirichlet_char_index(chi, G, n_randint(state, G->phi_q));
+        } while (!dirichlet_char_is_primitive(G, chi));
+
+        acb_poly_init(S);
+        acb_poly_init(A);
+        acb_poly_init(B);
+        acb_poly_init(C);
+
+        acb_poly_randtest(S, state, m, bits1, 3);
+
+        acb_dirichlet_hardy_theta_series(A, S, G, chi, n1, bits2);
+        acb_poly_set(B, S);  /* aliasing */
+        acb_dirichlet_hardy_theta_series(B, B, G, chi, n2, bits3);
+
+        acb_poly_set(C, A);
+        acb_poly_truncate(C, FLINT_MIN(n1, n2));
+        acb_poly_truncate(B, FLINT_MIN(n1, n2));
+
+        if (!acb_poly_overlaps(B, C))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("S = "); acb_poly_printd(S, 15); flint_printf("\n\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");
+            abort();
+        }
+
+        acb_poly_clear(S);
+        acb_poly_clear(A);
+        acb_poly_clear(B);
+        acb_poly_clear(C);
+
+        dirichlet_char_clear(chi);
+        dirichlet_group_clear(G);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_dirichlet/test/t-hardy_z_series.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2013 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_dirichlet.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("hardy_z_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
+    {
+        slong m, n1, n2, bits1, bits2, bits3;
+        acb_poly_t S, A, B, C;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+        ulong q;
+
+        bits1 = 2 + n_randint(state, 200);
+        bits2 = 2 + n_randint(state, 200);
+        bits3 = 2 + n_randint(state, 200);
+
+        m = 1 + n_randint(state, 8);
+        n1 = 1 + n_randint(state, 8);
+        n2 = 1 + n_randint(state, 8);
+
+        do {
+            q = 1 + n_randint(state, 15);
+        } while (q % 4 == 2);
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
+        do {
+            dirichlet_char_index(chi, G, n_randint(state, G->phi_q));
+        } while (!dirichlet_char_is_primitive(G, chi));
+
+        acb_poly_init(S);
+        acb_poly_init(A);
+        acb_poly_init(B);
+        acb_poly_init(C);
+
+        acb_poly_randtest(S, state, m, bits1, 3);
+
+        acb_dirichlet_hardy_z_series(A, S, G, chi, n1, bits2);
+        acb_poly_set(B, S);  /* aliasing */
+        acb_dirichlet_hardy_z_series(B, B, G, chi, n2, bits3);
+
+        acb_poly_set(C, A);
+        acb_poly_truncate(C, FLINT_MIN(n1, n2));
+        acb_poly_truncate(B, FLINT_MIN(n1, n2));
+
+        if (!acb_poly_overlaps(B, C))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("S = "); acb_poly_printd(S, 15); flint_printf("\n\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");
+            abort();
+        }
+
+        acb_poly_clear(S);
+        acb_poly_clear(A);
+        acb_poly_clear(B);
+        acb_poly_clear(C);
+
+        dirichlet_char_clear(chi);
+        dirichlet_group_clear(G);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_dirichlet/test/t-l_series.c

@@ -0,0 +1,83 @@
+/*
+    Copyright (C) 2013 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_dirichlet.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("l_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100 * arb_test_multiplier(); iter++)
+    {
+        slong m, n1, n2, bits1, bits2, bits3;
+        int deflate;
+        acb_poly_t S, A, B, C;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+
+        bits1 = 2 + n_randint(state, 200);
+        bits2 = 2 + n_randint(state, 200);
+        bits3 = 2 + n_randint(state, 200);
+
+        m = 1 + n_randint(state, 8);
+        n1 = 1 + n_randint(state, 8);
+        n2 = 1 + n_randint(state, 8);
+
+        dirichlet_group_init(G, 1 + n_randint(state, 12));
+        dirichlet_char_init(chi, G);
+        dirichlet_char_index(chi, G, n_randint(state, G->phi_q));
+
+        acb_poly_init(S);
+        acb_poly_init(A);
+        acb_poly_init(B);
+        acb_poly_init(C);
+
+        deflate = n_randint(state, 2);
+        acb_poly_randtest(S, state, m, bits1, 3);
+
+        acb_dirichlet_l_series(A, S, G, chi, deflate, n1, bits2);
+        acb_poly_set(B, S);  /* aliasing */
+        acb_dirichlet_l_series(B, B, G, chi, deflate, n2, bits3);
+
+        acb_poly_set(C, A);
+        acb_poly_truncate(C, FLINT_MIN(n1, n2));
+        acb_poly_truncate(B, FLINT_MIN(n1, n2));
+
+        if (!acb_poly_overlaps(B, C))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("S = "); acb_poly_printd(S, 15); flint_printf("\n\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");
+            abort();
+        }
+
+        acb_poly_clear(S);
+        acb_poly_clear(A);
+        acb_poly_clear(B);
+        acb_poly_clear(C);
+
+        dirichlet_char_clear(chi);
+        dirichlet_group_clear(G);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb_dirichlet.rst

@@ -419,6 +419,13 @@ L-functions
     gives the regular part of the Laurent expansion.
     When *chi* is non-principal, *deflate* has no effect.
 
+.. function:: void _acb_dirichlet_l_series(acb_ptr res, acb_srcptr s, slong slen, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)
+
+.. function:: void acb_dirichlet_l_series(acb_poly_t res, const acb_poly_t s, const dirichlet_group_t G, const dirichlet_char_t chi, int deflate, slong len, slong prec)
+
+    Sets *res* to the power series `L(s,\chi)` where *s* is a given power series, truncating the result to length *len*.
+    See :func:`acb_dirichlet_l_jet` for the meaning of the *deflate* flag.
+
 Hardy Z-functions
 -------------------------------------------------------------------------------
 
@@ -442,3 +449,15 @@ Hardy Z-functions
     `Z(t) = e^{i \theta(t)} L(1/2+it)`, which is real-valued for real *t*.
     The first *len* terms in the Taylor expansion are written to the output.
 
+.. function:: void _acb_dirichlet_hardy_theta_series(acb_ptr res, acb_srcptr t, slong tlen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)
+
+.. function:: void acb_dirichlet_hardy_theta_series(acb_poly_t res, const acb_poly_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)
+
+    Sets *res* to the power series `\theta(t)` where *t* is a given power series, truncating the result to length *len*.
+
+.. function:: void _acb_dirichlet_hardy_z_series(acb_ptr res, acb_srcptr t, slong tlen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)
+
+.. function:: void acb_dirichlet_hardy_z_series(acb_poly_t res, const acb_poly_t t, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec)
+
+    Sets *res* to the power series `Z(t)` where *t* is a given power series, truncating the result to length *len*.
+