Merge remote-tracking branch 'upstream/master' into cmake2

.build_dependencies

@@ -8,7 +8,7 @@ make -j4 > /dev/null 2>&1
 make install
 cd ..
 
-wget http://www.mpfr.org/mpfr-current/mpfr-3.1.4.tar.bz2
+wget http://www.mpfr.org/mpfr-3.1.4/mpfr-3.1.4.tar.bz2
 tar -xf mpfr-3.1.4.tar.bz2
 cd mpfr-3.1.4
 ./configure --with-gmp=$HOME/deps --prefix=$HOME/deps --disable-static

.gitignore

@@ -1,5 +1,6 @@
 Makefile
 build/*
 libarb.so*
+libarb.a
 doc/build/*
 *.ppm

.travis.yml

@@ -6,7 +6,7 @@ addons:
       - texinfo
 
 os:
-  - osx
+  # - osx
   - linux
 
 compiler:

Makefile.in

@@ -13,8 +13,8 @@ QUIET_AR  = @echo '   ' AR  ' ' $@;
 AT=@
 
 BUILD_DIRS = fmpr arf mag arb arb_mat arb_poly arb_calc acb acb_mat acb_poly \
-   acb_calc acb_hypgeom acb_modular acb_dirichlet bernoulli hypgeom \
-   fmpz_extras bool_mat partitions \
+   acb_calc acb_hypgeom acb_modular dirichlet acb_dirichlet arb_hypgeom bernoulli hypgeom \
+   fmpz_extras bool_mat partitions dlog \
    $(EXTRA_BUILD_DIRS)
 
 TEMPLATE_DIRS = 

README.md

@@ -8,7 +8,7 @@ GNU Lesser General Public License (LGPL), version 2.1 or later.
 
 ![arb logo](http://fredrikj.net/blog/2015/01/arb-2-5-0-released/arbtext.png)
 
-Documentation: http://fredrikj.net/arb/
+Documentation: http://arblib.org
 
 Development updates: http://fredrikj.net/blog/
 
@@ -70,7 +70,7 @@ can rely on Arb's automatic error bound tracking to get an output
 that is guaranteed to be accurate -- no error analysis
 needs to be done by the user.
 
-For several other example programs, see: http://fredrikj.net/arb/examples.html
+For more example programs, see: http://arblib.org/examples.html
 
 ## General features
 
@@ -140,7 +140,7 @@ Arb depends on FLINT (http://flintlib.org/), either
 GMP (http://gmplib.org) or MPIR (http://mpir.org),
 and MPFR (http://mpfr.org). 
 
-See http://fredrikj.net/arb/setup.html for instructions
+See http://arblib.org/setup.html for instructions
 on building and installing Arb directly from the source code.
 Arb might also be available (or coming soon) as a package for
 your Linux distribution.

acb.h

@@ -48,12 +48,7 @@ acb_init(acb_t x)
     arb_init(acb_imagref(x));
 }
 
-ACB_INLINE void
-acb_clear(acb_t x)
-{
-    arb_clear(acb_realref(x));
-    arb_clear(acb_imagref(x));
-}
+void acb_clear(acb_t x);
 
 ACB_INLINE acb_ptr
 _acb_vec_init(slong n)
@@ -1014,6 +1009,14 @@ acb_printd(const acb_t z, slong digits)
     acb_fprintd(stdout, z, digits);
 }
 
+void acb_fprintn(FILE * fp, const acb_t z, long digits, ulong flags);
+
+ACB_INLINE void
+acb_printn(const acb_t x, slong digits, ulong flags)
+{
+    acb_fprintn(stdout, x, digits, flags);
+}
+
 void acb_randtest(acb_t z, flint_rand_t state, slong prec, slong mag_bits);
 
 void acb_randtest_special(acb_t z, flint_rand_t state, slong prec, slong mag_bits);
@@ -1132,6 +1135,11 @@ _acb_vec_get_unique_fmpz_vec(fmpz * res,  acb_srcptr vec, slong len)
 /* sort complex numbers in a nice-to-display order */
 void _acb_vec_sort_pretty(acb_ptr vec, slong len);
 
+/* roots of unity */
+void acb_nth_root(acb_t res, ulong order, slong prec);
+void _acb_vec_nth_roots_pe(acb_ptr z, slong p, slong e, slong len, slong step, slong prec);
+void _acb_vec_nth_roots(acb_ptr z, slong len, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

acb/clear.c

@@ -0,0 +1,37 @@
+/*
+    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.h"
+
+/* The clear methods are not defined as inlines since this inflates
+   the code size. We inline the content manually here to avoid function
+   call overhead. Question: would calling arb_clear twice here be worth
+   the call overhead, by improving branch prediction? */
+void
+acb_clear(acb_t x)
+{
+    ARF_DEMOTE(arb_midref(acb_realref(x)));
+
+    if (COEFF_IS_MPZ(ARF_EXP(arb_midref(acb_realref(x)))))
+        _fmpz_clear_mpz(ARF_EXP(arb_midref(acb_realref(x))));
+
+    if (COEFF_IS_MPZ(MAG_EXP(arb_radref(acb_realref(x)))))
+        _fmpz_clear_mpz(MAG_EXP(arb_radref(acb_realref(x))));
+
+    ARF_DEMOTE(arb_midref(acb_imagref(x)));
+
+    if (COEFF_IS_MPZ(ARF_EXP(arb_midref(acb_imagref(x)))))
+        _fmpz_clear_mpz(ARF_EXP(arb_midref(acb_imagref(x))));
+
+    if (COEFF_IS_MPZ(MAG_EXP(arb_radref(acb_imagref(x)))))
+        _fmpz_clear_mpz(MAG_EXP(arb_radref(acb_imagref(x))));
+}
+

acb/fprintn.c

@@ -0,0 +1,50 @@
+/*
+    Copyright (C) 2015 Fredrik Johansson
+    Copyright (C) 2015 Arb authors
+
+    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.h"
+
+void
+acb_fprintn(FILE * file, const acb_t z, long digits, ulong flags)
+{
+    if (arb_is_zero(acb_imagref(z)))
+    {
+        arb_fprintn(file, acb_realref(z), digits, flags);
+    }
+    else if (arb_is_zero(acb_realref(z)))
+    {
+        arb_fprintn(file, acb_imagref(z), digits, flags);
+        flint_fprintf(file, "*I");
+    }
+    else
+    {
+        arb_fprintn(file, acb_realref(z), digits, flags);
+
+        if ((arb_is_exact(acb_imagref(z)) || (flags & ARB_STR_NO_RADIUS))
+                && arf_sgn(arb_midref(acb_imagref(z))) < 0)
+        {
+            arb_t t;
+            arb_init(t);
+            arb_neg(t, acb_imagref(z));
+            flint_fprintf(file, " - ");
+            arb_fprintn(file, t, digits, flags);
+            arb_clear(t);
+        }
+        else
+        {
+            flint_fprintf(file, " + ");
+            arb_fprintn(file, acb_imagref(z), digits, flags);
+        }
+
+        flint_fprintf(file, "*I");
+    }
+}
+

acb/nth_root.c

@@ -0,0 +1,42 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+static void
+_acb_nth_root(acb_t res, ulong order, slong prec)
+{
+    fmpq_t t;
+    fmpq_init(t);
+    fmpq_set_si(t, 2, order);
+    arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+    fmpq_clear(t);
+}
+
+void
+acb_nth_root(acb_t res, ulong order, slong prec)
+{
+    switch (order)
+    {
+       case 1:
+           acb_one(res);
+           break;
+       case 2:
+           acb_set_si(res, -1);
+           break;
+       case 4:
+           acb_onei(res);
+           break;
+       default:
+           _acb_nth_root(res, order, prec);
+           break;
+    }
+}

acb/sinc.c

@@ -59,6 +59,8 @@ _acb_sinc_diffbound(acb_t res, const acb_t z, slong prec)
 
     arf_set(arb_midref(acb_realref(res)), arb_midref(acb_realref(z)));
     arf_set(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(z)));
+    mag_zero(arb_radref(acb_realref(res)));
+    mag_zero(arb_radref(acb_imagref(res)));
 
     _acb_sinc_direct(res, res, prec);
 

acb/vec_nth_root.c

@@ -0,0 +1,91 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+/* assume xs has size >= len * step */
+void
+_acb_vec_set_powers_step(acb_ptr xs, slong n, slong len, slong step, slong prec)
+{
+    slong i, j;
+    prec += n_clog(n, 2);
+
+    for (i = 0, j = 0; i < len; i++, j += step)
+    {
+        if (i == 0)
+            acb_one(xs + j);
+        else if (i == 1)
+            acb_nth_root(xs + j, n, prec);
+        else if (i % 2 == 0)
+            acb_mul(xs + j, xs + j / 2, xs + j / 2, prec);
+        else
+            acb_mul(xs + j, xs + j - step, xs + step, prec);
+    }
+}
+
+/* assume len = p^e and z has size >= len * step */
+void
+_acb_vec_nth_roots_pe(acb_ptr z, slong p, slong e, slong len, slong step, slong prec)
+{
+    if (e <= 1)
+    {
+        _acb_vec_set_powers_step(z, p, len, step, prec);
+    }
+    else
+    {
+        slong q, r;
+        _acb_vec_nth_roots_pe(z, p, e - 1, len / p, step * p, prec);
+        _acb_vec_set_powers_step(z, n_pow(p, e), p, step, prec);
+
+        for (q = p; q < len; q += p)
+            for (r = 1; r < p; r++)
+                acb_mul(z + (q + r) * step, z + q * step, z + r * step, prec);
+    }
+}
+
+void
+_acb_vec_nth_roots(acb_ptr z, slong len, slong prec)
+{
+   slong i, q;
+   n_factor_t fac;
+
+   acb_one(z + 0);
+
+   n_factor_init(&fac);
+   n_factor(&fac, len, 0);
+   q = 1;
+
+   for (i = 0; i < fac.num; i++)
+   {
+       slong p, e, pe, mp, mq;
+       p = fac.p[i];
+       e = fac.exp[i];
+       pe = n_pow(p, e);
+       mp = len / pe;
+       mq = len / q;
+
+       _acb_vec_nth_roots_pe(z, p, e, pe, mp, prec);
+
+       if (i > 0)
+       {
+           slong k, l;
+
+           for (k = mp; k < len; k += mp)
+           {
+               for (l = mq; l < len - k; l += mq)
+                   acb_mul(z + k + l, z + k, z + l, prec);
+               for (; l < len; l += mq)
+                   acb_mul(z + k + l - len, z + k, z + l, prec);
+           }
+       }
+       q *= pe;
+   }
+}

acb_dirichlet.h

@@ -1,6 +1,7 @@
 /*
     Copyright (C) 2015 Jonathan Bober
     Copyright (C) 2016 Fredrik Johansson
+    Copyright (C) 2016 Pascal Molin
 
     This file is part of Arb.
 
@@ -13,43 +14,98 @@
 #ifndef ACB_DIRICHLET_H
 #define ACB_DIRICHLET_H
 
+#ifdef ACB_DIRICHLET_INLINES_C
+#define ACB_DIRICHLET_INLINE
+#else
+#define ACB_DIRICHLET_INLINE static __inline__
+#endif
+
 #include "acb.h"
+#include "dirichlet.h"
 
 #ifdef __cplusplus
 extern "C" {
 #endif
 
-typedef struct
-{
-    ulong q;                /* modulus */
-    ulong q_even;           /* even part of modulus */
-    ulong q_odd;            /* odd part of modulus */
-    ulong phi_q;            /* phi(q) */
-    ulong phi_q_odd;        /* phi(q_odd) */
-    slong num;              /* number of odd prime factors */
-    ulong * primes;         /* odd primes p[k] */
-    ulong * exponents;      /* exponents e[k] */
-    ulong * generators;     /* generator for each odd prime p[k] */
-    ulong * PHI;            /* PHI(k) = phi(q_odd) / phi(p[k]^e[k])  */
-}
-acb_dirichlet_group_struct;
-
-typedef acb_dirichlet_group_struct acb_dirichlet_group_t[1];
-
 void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s,
     const signed char * chi, int mod, int reciprocal, slong prec);
 
 void acb_dirichlet_eta(acb_t res, const acb_t s, slong prec);
 
-void acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q);
+void acb_dirichlet_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec);
+void acb_dirichlet_pairing_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec);
 
-void acb_dirichlet_group_clear(acb_dirichlet_group_t G);
+/* precompute powers of a root of unity */
+typedef struct
+{
+    ulong order;
+    acb_t z;
+    slong size;
+    slong depth;
+    acb_ptr * Z;
+}
+acb_dirichlet_powers_struct;
 
-void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec);
+typedef acb_dirichlet_powers_struct acb_dirichlet_powers_t[1];
+
+void _acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong size, slong depth, slong prec);
+void acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec);
+void acb_dirichlet_powers_clear(acb_dirichlet_powers_t t);
+void acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec);
+
+void acb_dirichlet_chi(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong n, slong prec);
+void acb_dirichlet_chi_vec(acb_ptr v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv, slong prec);
+
+void acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec);
+void acb_dirichlet_qseries_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec);
+void acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec);
+void acb_dirichlet_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec);
+
+ulong acb_dirichlet_theta_length_d(ulong q, double x, slong prec);
+ulong acb_dirichlet_theta_length(ulong q, const arb_t x, slong prec);
+void mag_tail_kexpk2_arb(mag_t res, const arb_t a, ulong n);
+
+void _acb_dirichlet_theta_argument_at_arb(arb_t xt, ulong q, const arb_t t, slong prec);
+void acb_dirichlet_theta_arb(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t t, slong prec);
+void acb_dirichlet_ui_theta_arb(acb_t res, const dirichlet_group_t G, ulong a, const arb_t t, slong prec);
+
+void acb_dirichlet_gauss_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+void acb_dirichlet_gauss_sum_factor(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+void acb_dirichlet_gauss_sum_order2(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+void acb_dirichlet_gauss_sum_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+void acb_dirichlet_gauss_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+
+void acb_dirichlet_root_number_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+void acb_dirichlet_root_number(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+
+void acb_dirichlet_si_poly_evaluate(acb_t res, slong * v, slong len, const acb_t z, slong prec);
+
+void acb_dirichlet_jacobi_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec);
+ulong jacobi_one_prime(ulong p, ulong e, ulong pe, ulong cond);
+void acb_dirichlet_jacobi_sum_factor(acb_t res,  const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec);
+void acb_dirichlet_jacobi_sum_gauss(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec);
+void acb_dirichlet_jacobi_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec);
+void acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong b, slong prec);
+
+void acb_dirichlet_l_euler_product(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+
+void acb_dirichlet_l_hurwitz(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+void acb_dirichlet_l(acb_t res, const acb_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+
+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);
+
+/* utils */
+
+ACB_DIRICHLET_INLINE void
+acb_vec_printd(acb_srcptr vec, slong len, slong digits)
+{
+    slong i;
+    for (i = 0; i < len; i++)
+        acb_printd(vec + i, digits), flint_printf("\n");
+}
 
 #ifdef __cplusplus
 }
 #endif
 
 #endif
-

acb_dirichlet/arb_quadratic_powers.c

@@ -0,0 +1,38 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+/* TODO: BSGS can reduce to nv mul */
+void
+acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec)
+{
+    slong i;
+    arb_t dx, x2;
+    arb_init(dx);
+    arb_init(x2);
+    arb_set(dx, x);
+    arb_mul(x2, x, x, prec);
+    for (i = 0; i < nv; i++)
+    {
+        if (i == 0)
+            arb_one(v + i);
+        else if (i == 1)
+            arb_set_round(v + i, x, prec);
+        else
+        {
+            arb_mul(dx, dx, x2, prec);
+            arb_mul(v + i, v + i - 1, dx, prec);
+        }
+    }
+    arb_clear(dx);
+    arb_clear(x2);
+}

acb_dirichlet/chi.c

@@ -1,6 +1,7 @@
 /*
     Copyright (C) 2015 Jonathan Bober
     Copyright (C) 2016 Fredrik Johansson
+    Copyright (C) 2016 Pascal Molin
 
     This file is part of Arb.
 
@@ -12,116 +13,19 @@
 
 #include "acb_dirichlet.h"
 
-/* todo: modular arithmetic */
-
-static ulong
-chi_odd_exponent(const acb_dirichlet_group_t G, ulong m, ulong n)
-{
-    ulong x, k, pk, gk, logm, logn;
-
-    x = 0;
-
-    for (k = 0; k < G->num; k++)
-    {
-        pk = n_pow(G->primes[k], G->exponents[k]);
-        gk = G->generators[k] % pk;
-
-        logm = n_discrete_log_bsgs(m % pk, gk, pk);
-        logn = n_discrete_log_bsgs(n % pk, gk, pk);
-
-        x = (x + G->PHI[k] * logm * logn) % G->phi_q_odd;
-    }
-
-    return x;
-}
-
-static ulong
-chi_even_exponent(const acb_dirichlet_group_t G, ulong m, ulong n)
-{
-    ulong x;
-    ulong q_even = G->q_even;
-
-    if (q_even <= 2)
-        return 0;
-
-    x = 0;
-
-    if ((m % 4 == 3) && (n % 4 == 3))
-        x = q_even / 8;
-
-    if (q_even > 4)
-    {
-        ulong g2, logm, logn;
-
-        g2 = 5;
-
-        if (m % 4 == 3)
-        {
-            m = n_negmod(m, q_even);
-        }
-
-        if (n % 4 == 3)
-        {
-            n = n_negmod(n, q_even);
-        }
-
-        logm = n_discrete_log_bsgs(m % q_even, g2, q_even);
-        logn = n_discrete_log_bsgs(n % q_even, g2, q_even);
-        x += logm * logn;
-    }
-
-    return x % (q_even / 4);
-}
-
 void
-acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec)
+acb_dirichlet_chi(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong n, slong prec)
 {
-    fmpq_t t, u;
-
-    ulong odd_part, even_part;
-    ulong q_even = G->q_even;
-    ulong q_odd = G->q_odd;
-
-    if ((q_even > 1 && (n % 2 == 0)) || (q_odd > 1 && (n_gcd(q_odd, n) != 1)))
-    {
+    ulong expo;
+    expo = dirichlet_chi(G, chi, n);
+    if (expo == DIRICHLET_CHI_NULL)
         acb_zero(res);
-        return;
-    }
-
-    if (q_even > 2)
-    {
-        if (q_even == 4)
-        {
-            if (m % 4 == 3 && n % 4 == 3)
-                even_part = q_even / 2;  /* -1 */
-            else
-                even_part = 0;           /* 1 */
-        }
-        else
-        {
-            even_part = 4 * chi_even_exponent(G, m % q_even, n % q_even);
-        }
-    }
     else
     {
-        even_part = 0;
+        fmpq_t t;
+        fmpq_init(t);
+        fmpq_set_si(t, 2 * expo , G->expo);
+        arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+        fmpq_clear(t);
     }
-
-    if (q_odd > 1)
-        odd_part = chi_odd_exponent(G, m % q_odd, n % q_odd);
-    else
-        odd_part = 0;
-
-    fmpq_init(t);
-    fmpq_init(u);
-
-    fmpq_set_si(t, 2 * even_part, q_even);
-    fmpq_set_si(u, 2 * odd_part, G->phi_q_odd);
-    fmpq_add(t, t, u);
-
-    arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
-
-    fmpq_clear(t);
-    fmpq_clear(u);
 }
-

acb_dirichlet/chi_vec.c

@@ -0,0 +1,38 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_chi_vec(acb_ptr v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv, slong prec)
+{
+    slong k;
+    ulong * a, order;
+    acb_dirichlet_powers_t t;
+
+    a = flint_malloc(nv * sizeof(ulong));
+    order = dirichlet_order_char(G, chi);
+    dirichlet_chi_vec_order(a, G, chi, order, nv);
+
+    acb_dirichlet_powers_init(t, order, nv, prec);
+
+    acb_zero(v + 0);
+    for (k = 0; k < nv; k++)
+    {
+        if (a[k] != DIRICHLET_CHI_NULL)
+            acb_dirichlet_power(v + k, t, a[k], prec);
+        else
+            acb_zero(v + k);
+    }
+
+    acb_dirichlet_powers_clear(t);
+    flint_free(a);
+}

acb_dirichlet/eta.c

@@ -34,4 +34,3 @@ acb_dirichlet_eta(acb_t res, const acb_t s, slong prec)
         acb_clear(t);
     }
 }
-

acb_dirichlet/gauss_sum.c

@@ -0,0 +1,112 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+static void
+gauss_sum_non_primitive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, ulong cond, slong prec)
+{
+
+    slong k, mu = 1;
+    ulong NN0 = G->q / cond;
+
+    /* G(chi) = mu(N/N0)chi0(N/N0)G(chi0) */
+
+    if (NN0 % 2 == 0)
+    {
+        if (G->q % 4)
+            mu = -1;
+        else
+        {
+            acb_zero(res);
+            return;
+        }
+    }
+
+    for (k = G->neven; k < G->num; k++)
+    {
+        ulong p = G->P[k].p;
+
+        if (G->P[k].e > 1 && NN0 % (p*p) == 0)
+        {
+            acb_zero(res);
+            return;
+        }
+
+        if (NN0 % p == 0)
+            mu *= -1;
+    }
+
+    if (chi->n == 1)
+    {
+        acb_set_si(res, mu);
+    }
+    else
+    {
+
+        dirichlet_group_t G0;
+        dirichlet_char_t chi0;
+        acb_t z;
+
+        /* primitive char associated to chi */
+        dirichlet_subgroup_init(G0, G, cond);
+        dirichlet_char_init(chi0, G0);
+        dirichlet_char_lower(chi0, G0, chi, G);
+
+        acb_init(z);
+        acb_dirichlet_gauss_sum(z, G0, chi0, prec);
+
+        acb_dirichlet_chi(res, G0, chi0, NN0, prec);
+        acb_mul(res, res, z, prec);
+        acb_mul_si(res, res, mu, prec);
+
+        dirichlet_group_clear(G0);
+        dirichlet_char_clear(chi0);
+        acb_clear(z);
+    }
+}
+
+void
+acb_dirichlet_gauss_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    ulong cond = dirichlet_conductor_char(G, chi);
+    /* TODO: no need, factor also does it... */
+    if (cond != G->q)
+    {
+        gauss_sum_non_primitive(res, G, chi, cond, prec);
+    }
+    else if (dirichlet_char_is_real(G, chi))
+    {
+        acb_dirichlet_gauss_sum_order2(res, G, chi, prec);
+    }
+    else if (G->num  > 1 && G->num > G->neven)
+    {
+        acb_dirichlet_gauss_sum_factor(res, G, chi, prec);
+    }
+    else
+    {
+        /* must be non primitive */
+        if (acb_dirichlet_theta_length_d(G->q, 1, prec) > G->q)
+            acb_dirichlet_gauss_sum_naive(res, G, chi, prec);
+        else
+            acb_dirichlet_gauss_sum_theta(res, G, chi, prec);
+    }
+}
+
+void
+acb_dirichlet_gauss_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, slong prec)
+{
+    dirichlet_char_t chi;
+    dirichlet_char_init(chi, G);
+    dirichlet_char_log(chi, G, a);
+    acb_dirichlet_gauss_sum(res, G, chi, prec);
+    dirichlet_char_clear(chi);
+}

acb_dirichlet/gauss_sum_factor.c

@@ -0,0 +1,70 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_gauss_sum_factor(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+
+    slong k;
+    acb_t tmp;
+
+    /* early check */
+    for (k = (G->neven == 2); k < G->num; k++)
+    {
+        /* if e > 1 and not primitive, 0 */
+        if (chi->log[k] % G->P[k].p == 0 && G->P[k].e > 1)
+        {
+            acb_zero(res);
+            return;
+        }
+    }
+
+    /* factor */
+    acb_one(res);
+    acb_init(tmp);
+
+    for (k = (G->neven == 2); k < G->num; k++)
+    {
+        ulong pe = G->P[k].pe.n;
+        dirichlet_group_t Gp;
+        dirichlet_char_t chip;
+
+        dirichlet_subgroup_init(Gp, G, pe);
+        dirichlet_char_init(chip, Gp);
+
+        chip->n = chi->n % pe;
+
+        if (k == 1 && G->neven == 2)
+        {
+            chip->log[0] = chi->log[0];
+            chip->log[1] = chi->log[1];
+        }
+        else
+            chip->log[0] = chi->log[k];
+
+        /* chi_pe(a, q/pe) * G_pe(a) */
+        acb_dirichlet_gauss_sum(tmp, Gp, chip, prec);
+        acb_mul(res, res, tmp, prec);
+
+        acb_dirichlet_chi(tmp, Gp, chip, (G->q / pe) % pe, prec);
+        acb_mul(res, res, tmp, prec);
+
+        dirichlet_char_clear(chip);
+        dirichlet_group_clear(Gp);
+    }
+
+    if (G->q_even == 2)
+        acb_neg(res, res);
+
+    acb_clear(tmp);
+}

acb_dirichlet/gauss_sum_naive.c

@@ -0,0 +1,31 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_gauss_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    acb_t z;
+    acb_ptr v;
+
+    v = _acb_vec_init(G->q);
+    acb_dirichlet_chi_vec(v, G, chi, G->q, prec);
+
+    acb_init(z);
+    acb_nth_root(z, G->q, prec);
+
+    _acb_poly_evaluate(res, v, G->q, z, prec);
+
+    acb_clear(z);
+    _acb_vec_clear(v, G->q);
+}

acb_dirichlet/gauss_sum_order2.c

@@ -0,0 +1,27 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_gauss_sum_order2(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    if (dirichlet_parity_char(G, chi))
+    {
+        arb_zero(acb_realref(res));
+        arb_sqrt_ui(acb_imagref(res), G->q, prec);
+    }
+    else
+    {
+        arb_zero(acb_imagref(res));
+        arb_sqrt_ui(acb_realref(res), G->q, prec);
+    }
+}

acb_dirichlet/gauss_sum_theta.c

@@ -0,0 +1,35 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_gauss_sum_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    ulong cond = dirichlet_conductor_char(G, chi);
+    if (cond < G->q || (G->q == 300 && (chi->n == 71 || chi->n == 131))
+            || (G->q == 600 && (chi->n == 11 || chi->n == 491)))
+    {
+        flint_printf("gauss_sum_theta: non available for non primitive character"
+                "or exceptional characters chi_300(71,.), chi_300(131,.), "
+                "chi_600(11,.) and chi_600(491,.)\n");
+        abort();
+    }
+    else
+    {
+        acb_t iq;
+        acb_init(iq);
+        acb_dirichlet_gauss_sum_order2(iq, G, chi, prec);
+        acb_dirichlet_root_number_theta(res, G, chi, prec);
+        acb_mul(res, res, iq, prec);
+        acb_clear(iq);
+    }
+}

acb_dirichlet/group_init.c

@@ -1,84 +0,0 @@
-/*
-    Copyright (C) 2015 Jonathan Bober
-    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"
-
-static ulong
-primitive_root_p_and_p2(ulong p)
-{
-    if (p == 40487)
-        return 10;
-
-#if FLINT_BITS == 64
-    if (p == UWORD(6692367337))
-        return 7;
-
-    if (p > UWORD(1000000000000))
-    {
-        printf("primitive root: p > 10^12 not implemented");
-        abort();
-    }
-#endif
-
-    return n_primitive_root_prime(p);
-}
-
-void
-acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q)
-{
-    slong k;
-    n_factor_t fac;
-
-    G->q = q;
-    G->q_odd = q;
-    G->q_even = 1;
-
-    while (G->q_odd % 2 == 0)
-    {
-        G->q_odd /= 2;
-        G->q_even *= 2;
-    }
-
-    n_factor_init(&fac);
-    n_factor(&fac, G->q_odd, 1);
-
-    G->num = fac.num;
-    G->primes = flint_malloc(G->num * sizeof(ulong));
-    G->exponents = flint_malloc(G->num * sizeof(ulong));
-    G->generators = flint_malloc(G->num * sizeof(ulong));
-    G->PHI = flint_malloc(G->num * sizeof(ulong));
-
-    for (k = 0; k < G->num; k++)
-    {
-        G->primes[k] = fac.p[k];
-        G->exponents[k] = fac.exp[k];
-    }
-
-    G->phi_q_odd = 1;
-    for (k = 0; k < G->num; k++)
-        G->phi_q_odd *= (G->primes[k] - 1) * n_pow(G->primes[k], G->exponents[k]-1);
-
-    if (G->q_even == 1)
-        G->phi_q = G->phi_q_odd;
-    else
-        G->phi_q = G->phi_q_odd * (G->q_even / 2);
-
-    for (k = 0; k < G->num; k++)
-    {
-        ulong phi;
-
-        G->generators[k] = primitive_root_p_and_p2(G->primes[k]);
-        phi = n_pow(G->primes[k], G->exponents[k] - 1) * (G->primes[k] - 1);
-        G->PHI[k] = G->phi_q_odd / phi;
-    }
-}
-

acb_dirichlet/jacobi_sum.c

@@ -0,0 +1,109 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+/* J_N(1,a) = sum on x = 1 mod some p | q */
+ulong
+jacobi_one_prime(ulong p, ulong e, ulong pe, ulong cond)
+{
+    if (e > 1 && cond % (p*p) == 0)
+    {
+        return 0;
+    }
+    else
+    {
+        slong r = (e > 1) ? pe / p : 1;
+        if (cond % p)
+            return r * (p - 2);
+        else
+            return -r;
+    }
+}
+
+static ulong
+jacobi_one(const dirichlet_group_t G, ulong cond)
+{
+    slong k, r = 1;
+
+    for (k = 0; k < G->num; k++)
+        r *= jacobi_one_prime(G->P[k].p, G->P[k].e,
+                G->P[k].pe.n, cond);
+    return r;
+}
+
+void
+acb_dirichlet_jacobi_sum(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)
+{
+    if (G->q_even > 1)
+    {
+        acb_zero(res);
+    }
+    else if (chi1->n == 1 || chi2->n == 1)
+    {
+        ulong cond = (chi1->n == 1) ? dirichlet_conductor_char(G, chi2) : dirichlet_conductor_char(G, chi1);
+        acb_set_si(res, jacobi_one(G, cond));
+    }
+    else if (nmod_mul(chi1->n, chi2->n, G->mod) == 1)
+    {
+        ulong n;
+        n = jacobi_one(G, dirichlet_conductor_char(G, chi1));
+        if (dirichlet_parity_char(G, chi1))
+            acb_set_si(res, -n);
+        else
+            acb_set_si(res, n);
+    }
+    else
+    {
+        if (G->q <= 150)
+            acb_dirichlet_jacobi_sum_naive(res, G, chi1, chi2, prec);
+        else if (G->num > 1)
+            acb_dirichlet_jacobi_sum_factor(res, G, chi1, chi2, prec);
+        else if (G->P[0].e > 1)
+            acb_dirichlet_jacobi_sum_naive(res, G, chi1, chi2, prec);
+        else
+            acb_dirichlet_jacobi_sum_gauss(res, G, chi1, chi2, prec);
+    }
+}
+
+void
+acb_dirichlet_jacobi_sum_ui(acb_t res, const dirichlet_group_t G, ulong a, ulong b, slong prec)
+{
+    if (G->q_even > 1)
+    {
+        acb_zero(res);
+    }
+    else if (a == 1 || b == 1)
+    {
+        ulong cond = (a == 1) ? dirichlet_conductor_ui(G, b) : dirichlet_conductor_ui(G, a);
+        acb_set_si(res, jacobi_one(G, cond));
+    }
+    else if (nmod_mul(a, b, G->mod) == 1)
+    {
+        ulong n;
+        n = jacobi_one(G, dirichlet_conductor_ui(G, a));
+        if (dirichlet_parity_ui(G, a))
+            acb_set_si(res, -n);
+        else
+            acb_set_si(res, n);
+    }
+    else
+    {
+        dirichlet_char_t chi1, chi2;
+        dirichlet_char_init(chi1, G);
+        dirichlet_char_init(chi2, G);
+        dirichlet_char_log(chi1, G, a);
+        dirichlet_char_log(chi2, G, b);
+        acb_dirichlet_jacobi_sum(res, G, chi1, chi2, prec);
+        dirichlet_char_clear(chi1);
+        dirichlet_char_clear(chi2);
+    }
+}

acb_dirichlet/jacobi_sum_factor.c

@@ -0,0 +1,76 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_jacobi_sum_factor(acb_t res,  const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)
+{
+    slong k;
+    acb_t tmp;
+    acb_init(tmp);
+    acb_one(res);
+
+    /* TODO: efficient subgroup */
+    for (k = 0; k < G->num; k++)
+    {
+        nmod_t pe;
+        ulong p, e, ap, bp;
+
+        p = G->P[k].p;
+        e = G->P[k].e;
+        pe = G->P[k].pe;
+        ap = chi1->n % pe.n;
+        bp = chi2->n % pe.n;
+
+        if (ap == 1 || bp == 1 || nmod_mul(ap, bp, pe) == 1)
+        {
+            slong r;
+            ulong cond;
+
+            cond = (ap == 1) ? dirichlet_conductor_char(G, chi2) : dirichlet_conductor_char(G, chi1);
+            r = jacobi_one_prime(p, e, pe.n, cond);
+
+            /* chi(a,-1) if ap * bp = 1 */
+            if (ap != 1 && bp != 1)
+                r *= n_jacobi_unsigned(ap, p);
+
+            acb_mul_si(res, res, r, prec);
+        }
+        else
+        {
+            dirichlet_group_t Gp;
+            dirichlet_char_t chi1p, chi2p;
+
+            dirichlet_group_init(Gp, pe.n);
+            dirichlet_char_init(chi1p, Gp);
+            dirichlet_char_init(chi2p, Gp);
+
+            chi1p->n = ap;
+            chi1p->log[0] = chi1->log[k];
+            chi2p->n = ap;
+            chi2p->log[0] = chi2->log[k];
+
+            /* TODO: work out gauss relations for e > 1 */
+            if (p <= 100 || e > 1)
+                acb_dirichlet_jacobi_sum_naive(tmp, Gp, chi1p, chi2p, prec);
+            else
+                acb_dirichlet_jacobi_sum_gauss(tmp, Gp, chi1p, chi2p, prec);
+
+            acb_mul(res, res, tmp, prec);
+
+            dirichlet_char_clear(chi1p);
+            dirichlet_char_clear(chi2p);
+            dirichlet_group_clear(Gp);
+        }
+    }
+    acb_clear(tmp);
+}

acb_dirichlet/jacobi_sum_gauss.c

@@ -0,0 +1,38 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+/* should use only for prime power modulus */
+void
+acb_dirichlet_jacobi_sum_gauss(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)
+{
+    /* J_q(a,b)G_q(ab) = G_q(a)G_q(b) */
+    acb_t tmp;
+    dirichlet_char_t chi12;
+
+    dirichlet_char_init(chi12, G);
+    dirichlet_char_mul(chi12, G, chi1, chi2);
+
+    acb_init(tmp);
+
+    acb_dirichlet_gauss_sum(res, G, chi1, prec);
+    if (chi2->n == chi1->n)
+        acb_set(tmp, res);
+    else
+        acb_dirichlet_gauss_sum(tmp, G, chi2, prec);
+    acb_mul(res, res, tmp, prec);
+    acb_dirichlet_gauss_sum(tmp, G, chi12, prec);
+    acb_div(res, res, tmp, prec);
+
+    dirichlet_char_clear(chi12);
+    acb_clear(tmp);
+}

acb_dirichlet/jacobi_sum_naive.c

@@ -0,0 +1,61 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_jacobi_sum_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi1, const dirichlet_char_t chi2, slong prec)
+{
+
+    ulong k1, k2, m1, m2, g, e, m;
+    ulong * v1, * v2;
+    slong *v;
+    nmod_t expo;
+    acb_t z;
+
+    v1 = flint_malloc(G->q * sizeof(ulong));
+    v2 = flint_malloc(G->q * sizeof(ulong));
+
+    dirichlet_vec_set_null(v1, G, G->q);
+    dirichlet_chi_vec_loop(v1, G, chi1, G->q);
+
+    dirichlet_vec_set_null(v2, G, G->q);
+    dirichlet_chi_vec_loop(v2, G, chi2, G->q);
+
+    nmod_init(&expo, G->expo);
+    m1 = dirichlet_order_char(G, chi1);
+    m2 = dirichlet_order_char(G, chi2);
+    g = m1 * m2 / n_gcd(m1, m2);
+    m = G->expo / g;
+
+    v = flint_malloc(g * sizeof(slong));
+
+    for (e = 0; e < g; e++)
+        v[e] = 0;
+
+    for (k1 = 2, k2 = G->q - 1; k2 > 1; k1++, k2--)
+    {
+        if (v1[k1] == DIRICHLET_CHI_NULL ||
+            v2[k2] == DIRICHLET_CHI_NULL)
+            continue;
+        e = nmod_add(v1[k1], v2[k2], expo) / m;
+        v[e]++;
+    }
+
+    acb_init(z);
+    acb_nth_root(z, g, prec);
+    acb_dirichlet_si_poly_evaluate(res, v, g, z, prec);
+
+    acb_clear(z);
+    flint_free(v);
+    flint_free(v2);
+    flint_free(v1);
+}

acb_dirichlet/l.c

@@ -1,5 +1,6 @@
 /*
     Copyright (C) 2016 Fredrik Johansson
+    Copyright (C) 2016 Pascal Molin
 
     This file is part of Arb.
 
@@ -13,39 +14,16 @@
 
 void
 acb_dirichlet_l(acb_t res, const acb_t s,
-    const acb_dirichlet_group_t G, ulong m, slong prec)
+    const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
 {
-    acb_t chi, t, u, a;
-    ulong k;
-
-    acb_init(chi);
-    acb_init(t);
-    acb_init(u);
-    acb_init(a);
-
-    acb_zero(t);
-
-    for (k = 1; k <= G->q; k++)
+    /* this cutoff is probably too conservative when q is large */
+    if (arf_cmp_d(arb_midref(acb_realref(s)), 8 + 0.5 * prec / log(prec)) >= 0)
     {
-        acb_dirichlet_chi(chi, G, m, k, prec);
-
-        if (!acb_is_zero(chi))
-        {
-            acb_set_ui(a, k);
-            acb_div_ui(a, a, G->q, prec);
-            acb_hurwitz_zeta(u, s, a, prec);
-            acb_addmul(t, chi, u, prec);
-        }
+        acb_dirichlet_l_euler_product(res, s, G, chi, prec);
+    }
+    else
+    {
+        acb_dirichlet_l_hurwitz(res, s, G, chi, prec);
     }
-
-    acb_set_ui(u, G->q);
-    acb_neg(a, s);
-    acb_pow(u, u, a, prec);
-    acb_mul(res, t, u, prec);
-
-    acb_clear(chi);
-    acb_clear(t);
-    acb_clear(u);
-    acb_clear(a);
 }
 

acb_dirichlet/l_euler_product.c

@@ -0,0 +1,128 @@
+/*
+    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"
+
+#define ONE_OVER_LOG2 1.4426950408889634
+
+void
+acb_dirichlet_l_euler_product(acb_t res, const acb_t s,
+    const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    arf_t left;
+    slong wp, powprec, left_s;
+    ulong p, p_limit;
+    double p_needed_approx, powmag, logp, errmag;
+    int is_real;
+    acb_t t, u, v, c;
+    mag_t err;
+
+    if (!acb_is_finite(s))
+    {
+        acb_indeterminate(res);
+        return;
+    }
+
+    arf_init(left);
+    arf_set_mag(left, arb_radref(acb_realref(s)));
+    arf_sub(left, arb_midref(acb_realref(s)), left, 2 * MAG_BITS, ARF_RND_FLOOR);
+
+    /* Require re(s) >= 2. */
+    if (arf_cmp_si(left, 2) < 0)
+    {
+        acb_indeterminate(res);
+        arf_clear(left);
+        return;
+    }
+
+    is_real = acb_is_real(s) && dirichlet_char_is_real(G, chi);
+
+    /* L(s) ~= 1. */
+    if (arf_cmp_si(left, prec) > 0)
+    {
+        acb_one(res);
+        mag_hurwitz_zeta_uiui(arb_radref(acb_realref(res)), prec, 2);
+        if (!is_real)
+            mag_set(arb_radref(acb_imagref(res)), arb_radref(acb_realref(res)));
+        acb_inv(res, res, prec);
+        arf_clear(left);
+        return;
+    }
+
+    left_s = arf_get_si(left, ARF_RND_FLOOR);
+
+    /* Heuristic. */
+    wp = prec + FLINT_BIT_COUNT(prec) + (prec / left_s) + 4;
+
+    /* Don't work too hard if a small s was passed as input. */
+    p_limit = 100 + prec * sqrt(prec);
+
+    /* Truncating at p ~= 2^(prec/s) gives an error of 2^-prec */
+    if (((double) prec) / left_s > 50.0)
+        p_needed_approx = pow(2.0, 50.0);
+    else
+        p_needed_approx = pow(2.0, ((double) prec) / left_s);
+
+    p_needed_approx = FLINT_MIN(p_limit, p_needed_approx);
+    /* todo: use this to determine whether to precompute chi */
+
+    acb_init(t);
+    acb_init(u);
+    acb_init(v);
+    acb_init(c);
+
+    acb_one(v);
+
+    for (p = 2; p < p_limit; p = n_nextprime(p, 1))
+    {
+        /* p^s */
+        logp = log(p);
+        powmag = left_s * logp * ONE_OVER_LOG2;
+
+        /* zeta(s,p) ~= 1/p^s + 1/((s-1) p^(s-1)) */
+        errmag = (log(left_s - 1.0) + (left_s - 1.0) * logp) * ONE_OVER_LOG2;
+        errmag = FLINT_MIN(powmag, errmag);
+
+        if (errmag > prec + 2)
+            break;
+
+        powprec = FLINT_MAX(wp - powmag, 8);
+
+        acb_dirichlet_chi(c, G, chi, p, powprec);
+
+        if (!acb_is_zero(c))
+        {
+            acb_set_ui(t, p);
+            acb_pow(t, t, s, powprec);
+            acb_set_round(u, v, powprec);
+            acb_div(t, u, t, powprec);
+            acb_mul(t, t, c, powprec);
+            acb_sub(v, v, t, wp);
+        }
+    }
+
+    mag_init(err);
+    mag_hurwitz_zeta_uiui(err, left_s, p);
+    if (is_real)
+        arb_add_error_mag(acb_realref(v), err);
+    else
+        acb_add_error_mag(v, err);
+    mag_clear(err);
+
+    acb_inv(res, v, prec);
+
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(v);
+    acb_clear(c);
+    arf_clear(left);
+}
+

acb_dirichlet/l_hurwitz.c

@@ -0,0 +1,80 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_l_hurwitz(acb_t res, const acb_t s,
+    const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    ulong order, chin, mult;
+    acb_t t, u, a;
+    acb_ptr z;
+    dirichlet_char_t cn;
+    int deflate;
+
+    /* remove pole in Hurwitz zeta at s = 1 */
+    deflate = 0;
+    if (acb_is_one(s))
+    {
+        if (dirichlet_char_is_principal(G, chi))
+        {
+            acb_indeterminate(res);
+            return;
+        }
+        deflate = 1;
+    }
+
+    dirichlet_char_init(cn, G);
+    acb_init(t);
+    acb_init(u);
+    acb_init(a);
+
+    dirichlet_char_one(cn, G);
+    acb_zero(t);
+
+    prec += n_clog(G->phi_q, 2);
+
+    order = dirichlet_order_char(G, chi);
+    mult = G->expo / order;
+    z = _acb_vec_init(order);
+    _acb_vec_nth_roots(z, order, prec);
+
+    do {
+        chin = dirichlet_pairing_char(G, chi, cn) / mult;
+
+        acb_set_ui(a, cn->n);
+        acb_div_ui(a, a, G->q, prec);
+
+        if (deflate == 0)
+            acb_hurwitz_zeta(u, s, a, prec);
+        else
+            _acb_poly_zeta_cpx_series(u, s, a, 1, 1, prec);
+
+        acb_addmul(t, z + chin, u, prec);
+
+    } while (dirichlet_char_next(cn, G) >= 0);
+
+    acb_set_ui(u, G->q);
+    acb_neg(a, s);
+    acb_pow(u, u, a, prec);
+    acb_mul(res, t, u, prec);
+
+    dirichlet_char_clear(cn);
+
+    _acb_vec_clear(z, order);
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(a);
+}
+

acb_dirichlet/l_jet.c

@@ -0,0 +1,119 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include "acb_poly.h"
+
+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)
+{
+    ulong order, chin, mult, phi;
+    acb_t a;
+    acb_ptr t, u, z;
+    dirichlet_char_t cn;
+    int deflate_hurwitz;
+
+    if (len <= 0)
+        return;
+
+    if (G->q == 1)
+    {
+        acb_init(a);
+        acb_one(a);
+        _acb_poly_zeta_cpx_series(res, s, a, deflate, len, prec);
+        acb_clear(a);
+        return;
+    }
+
+    if (len == 1 && !(deflate && dirichlet_char_is_principal(G, chi)))
+    {
+        acb_dirichlet_l(res, s, G, chi, prec);
+        return;
+    }
+
+    if (dirichlet_char_is_principal(G, chi))
+        deflate_hurwitz = deflate;
+    else
+        deflate_hurwitz = acb_is_one(s);
+
+    dirichlet_char_init(cn, G);
+    t = _acb_vec_init(len);
+    u = _acb_vec_init(len + 2);
+    acb_init(a);
+
+    dirichlet_char_one(cn, G);
+
+    prec += n_clog(G->phi_q, 2);
+
+    order = dirichlet_order_char(G, chi);
+    mult = G->expo / order;
+    z = _acb_vec_init(order);
+    _acb_vec_nth_roots(z, order, prec);
+
+    phi = 0;
+    do
+    {
+        chin = dirichlet_pairing_char(G, chi, cn) / mult;
+        acb_set_ui(a, cn->n);
+        acb_div_ui(a, a, G->q, prec);
+        _acb_poly_zeta_cpx_series(u, s, a, deflate_hurwitz, len, prec);
+        _acb_vec_scalar_addmul(t, u, len, z + chin, prec);
+        phi++;
+    }
+    while (dirichlet_char_next(cn, G) >= 0);
+
+    if (dirichlet_char_is_principal(G, chi) && deflate)
+    {
+        /* res = t * q^(-(s+x)) + [phi(q) * (q^(-(s+x)) - q^-1) / ((s+x)-1)] */
+        if (acb_is_one(s))
+        {
+            acb_set_ui(a, G->q);
+            _acb_poly_acb_invpow_cpx(u, a, s, len + 1, prec);
+            _acb_poly_mullow(res, t, len, u, len, len, prec);
+            acb_set_ui(u, phi);
+            _acb_vec_scalar_addmul(res, u + 1, len, u, prec);
+        }
+        else
+        {
+            acb_sub_ui(u, s, 1, prec);
+            acb_one(u + 1);
+
+            acb_set_ui(a, G->q);
+            _acb_poly_acb_invpow_cpx(u + 2, a, s, len, prec);
+            _acb_poly_mullow(res, t, len, u + 2, len, len, prec);
+
+            acb_inv(a, a, prec);
+            acb_sub(u + 2, u + 2, a, prec);
+            _acb_poly_div_series(t, u + 2, len, u, 2, len, prec);
+
+            acb_set_ui(u, phi);
+            _acb_vec_scalar_addmul(res, t, len, u, prec);
+        }
+    }
+    else
+    {
+        /* res = t * q^(-(s+x)) */
+        acb_set_ui(a, G->q);
+        _acb_poly_acb_invpow_cpx(u, a, s, len, prec);
+        _acb_poly_mullow(res, t, len, u, len, len, prec);
+    }
+
+    dirichlet_char_clear(cn);
+
+    _acb_vec_clear(z, order);
+    _acb_vec_clear(t, len);
+    _acb_vec_clear(u, len + 2);
+    acb_clear(a);
+}
+

acb_dirichlet/pairing.c

@@ -0,0 +1,29 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_pairing(acb_t res, const dirichlet_group_t G, ulong m, ulong n, slong prec)
+{
+    ulong expo;
+    expo = dirichlet_pairing(G, m, n);
+    if (expo == DIRICHLET_CHI_NULL)
+        acb_zero(res);
+    else
+    {
+        fmpq_t t;
+        fmpq_init(t);
+        fmpq_set_si(t, 2 * expo, G->expo);
+        arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+        fmpq_clear(t);
+    }
+}

acb_dirichlet/pairing_conrey.c

@@ -0,0 +1,29 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_pairing_char(acb_t res, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b, slong prec)
+{
+    ulong expo;
+    expo = dirichlet_pairing_char(G, a, b);
+    if (expo == DIRICHLET_CHI_NULL)
+        acb_zero(res);
+    else
+    {
+        fmpq_t t;
+        fmpq_init(t);
+        fmpq_set_si(t, 2 * expo, G->expo);
+        arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+        fmpq_clear(t);
+    }
+}

acb_dirichlet/power.c

@@ -0,0 +1,36 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec)
+{
+    if (!t->depth)
+    {
+        acb_pow_ui(z, t->z, n, prec);
+    }
+    else
+    {
+        slong k;
+        ulong r;
+        r = n % t->size;
+        n = n / t->size;
+        acb_set(z, t->Z[0] + r);
+        for (k = 1; k < t->depth && n; k++)
+        {
+            r = n % t->size;
+            n = n / t->size;
+            acb_mul(z, z, t->Z[k] + r, prec);
+        }
+    }
+}

acb_dirichlet/powers_clear.c

@@ -0,0 +1,24 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_powers_clear(acb_dirichlet_powers_t t)
+{
+    slong k;
+
+    for (k = 0; k < t->depth; k++)
+        _acb_vec_clear(t->Z[k], t->size);
+
+    flint_free(t->Z);
+    acb_clear(t->z);
+}

acb_dirichlet/powers_init.c

@@ -0,0 +1,51 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include "acb_poly.h"
+
+void
+_acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong size, slong depth, slong prec)
+{
+    slong k;
+    t->order = order;
+    acb_init(t->z);
+    acb_nth_root(t->z, order, prec);
+    t->size = size;
+    t->depth = depth;
+    if (depth)
+    {
+        acb_struct * z;
+        z = t->z + 0;
+        t->Z = flint_malloc(depth * sizeof(acb_ptr));
+        for (k = 0; k < depth; k++)
+        {
+            t->Z[k] = _acb_vec_init(size);
+            _acb_vec_set_powers(t->Z[k], z, size, prec);
+            z = t->Z[k] + 1;
+        }
+    }
+    else
+    {
+        t->Z = NULL;
+    }
+}
+
+void
+acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec)
+{
+    slong size, depth;
+
+    depth = (num > order) ? 1 : n_flog(order, num);
+    size = n_root(order, depth) + 1;
+
+    _acb_dirichlet_powers_init(t, order, size, depth, prec);
+}

acb_dirichlet/profile/p-gauss.c

@@ -0,0 +1,128 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 <string.h>
+#include "acb_dirichlet.h"
+#include "profiler.h"
+
+#define LOG 0
+#define CSV 1
+#define JSON 2
+
+typedef void (*do_f) (acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
+
+int main(int argc, char *argv[])
+{
+    int out;
+    ulong n, nref, maxq = 5000;
+    ulong * rand;
+    flint_rand_t state;
+    slong r, nr;
+
+    int l, nf = 5;
+    do_f func[5] = { acb_dirichlet_gauss_sum_naive,
+        acb_dirichlet_gauss_sum_factor,
+        acb_dirichlet_gauss_sum_theta,
+        acb_dirichlet_gauss_sum,
+        acb_dirichlet_gauss_sum };
+    char * name[5] = { "naive", "factor", "theta", "default", "precomp" };
+
+    int i, ni = 5;
+    ulong qmin[5] =  {  3,  50, 500, 1000, 10000 };
+    ulong qmax[5] =  { 50, 100, 599, 1050, 10030 };
+
+    int j, nj = 3;
+    slong prec[3] = { 64, 128, 1024 };
+
+    if (argc < 2)
+        out = LOG;
+    else if (!strcmp(argv[1], "json"))
+        out = JSON;
+    else if (!strcmp(argv[1], "csv"))
+        out = CSV;
+    else if (!strcmp(argv[1], "log"))
+        out = LOG;
+    else
+    {
+        printf("usage: %s [log|csv|json]\n", argv[0]);
+        abort();
+    }
+
+    if (out == CSV)
+        flint_printf("# %-12s, %7s, %7s, %7s, %7s\n","name", "prec", "qmin", "qmax", "time");
+
+    for (j = 0; j < nj; j++)
+    {
+
+    for (i = 0; i < ni; i++)
+    {
+
+        if (out == LOG)
+            flint_printf("G_q(a) at prec %wu for all %wu <= q <= %wu....\n", prec[j], qmin[i], qmax[i]);
+
+        for (l = 0; l < nf; l++)
+        {
+            ulong q;
+
+            if (qmin[i] > 2000 && l == 0)
+                continue;
+
+            if (out == LOG)
+                flint_printf("%-14s ...  ", name[l]);
+            else if (out == CSV)
+                flint_printf("%-12s, %7d, %7d, %7d,   ", name[l], prec[j], qmin[i], qmax[i]);
+            else if (out == JSON)
+                flint_printf("{ \"name\": \"%s\", \"prec\": %d, \"qmin\": %d, \"qmax\": %d, \"time\": ",
+                        name[l], prec[j], qmin[i], qmax[i]);
+
+            TIMEIT_ONCE_START
+
+                for (q = qmin[i]; q <= qmax[i]; q++)
+                {
+                    acb_dirichlet_group_t G;
+                    acb_dirichlet_char_t chi;
+                    acb_t res;
+
+                    if (q % 4 == 2)
+                        continue;
+
+                    acb_dirichlet_group_init(G, q);
+                    if (l == 4)
+                        acb_dirichlet_group_dlog_precompute(G, 1);
+                    acb_dirichlet_char_init(chi, G);
+                    acb_dirichlet_char_first_primitive(chi, G);
+                    acb_init(res);
+
+                        do {
+                            func[l](res, G, chi, prec[j]);
+                        } while (acb_dirichlet_char_next_primitive(chi, G) >= 0);
+
+                    acb_clear(res);
+                    acb_dirichlet_char_clear(chi);
+                    if (l == 4)
+                        acb_dirichlet_group_dlog_clear(G);
+                    acb_dirichlet_group_clear(G);
+                }
+
+            TIMEIT_ONCE_STOP
+
+                if (out == JSON)
+                    flint_printf("}\n");
+                else
+                    flint_printf("\n");
+        }
+
+    }
+
+    }
+    flint_cleanup();
+    return EXIT_SUCCESS;
+}

acb_dirichlet/profile/p-jacobi.c

@@ -0,0 +1,124 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 <string.h>
+#include "acb_dirichlet.h"
+#include "profiler.h"
+
+#define LOG 0
+#define CSV 1
+#define JSON 2
+
+typedef void (*do_f) (acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec);
+
+int main(int argc, char *argv[])
+{
+    int out;
+    ulong n, nref, maxq = 5000;
+    ulong * rand;
+    flint_rand_t state;
+    slong r, nr;
+
+    int l, nf = 4;
+    do_f func[4] = { acb_dirichlet_jacobi_sum_naive,
+        acb_dirichlet_jacobi_sum_factor,
+        acb_dirichlet_jacobi_sum_gauss,
+        acb_dirichlet_jacobi_sum };
+    char * name[4] = { "naive", "factor", "gauss", "default" };
+
+    int i, ni = 4;
+    ulong qmin[4] =  {  3,  51, 501, 10001 };
+    ulong qmax[4] =  { 30, 100, 550, 10050 };
+
+    int j, nj = 3;
+    slong prec[3] = { 64, 128, 1024 };
+
+    if (argc < 2)
+        out = LOG;
+    else if (!strcmp(argv[1], "json"))
+        out = JSON;
+    else if (!strcmp(argv[1], "csv"))
+        out = CSV;
+    else if (!strcmp(argv[1], "log"))
+        out = LOG;
+    else
+    {
+        printf("usage: %s [log|csv|json]\n", argv[0]);
+        abort();
+    }
+
+    if (out == CSV)
+        flint_printf("# %-12s, %7s, %7s, %7s, %7s\n","name", "prec", "qmin", "qmax", "time");
+
+    for (j = 0; j < nj; j++)
+    {
+
+    for (i = 0; i < ni; i++)
+    {
+
+        if (out == LOG)
+            flint_printf("J_q(a,b) at prec %wu for all %wu <= q <= %wu....\n", prec[j], qmin[i], qmax[i]);
+
+        for (l = 0; l < nf; l++)
+        {
+            ulong q;
+
+            if (out == LOG)
+                flint_printf("%-14s ...  ", name[l]);
+            else if (out == CSV)
+                flint_printf("%-12s, %7d, %7d, %7d,   ", name[l], prec[j], qmin[i], qmax[i]);
+            else if (out == JSON)
+                flint_printf("{ \"name\": \"%s\", \"prec\": %d, \"qmin\": %d, \"qmax\": %d, \"time\": ",
+                        name[l], prec[j], qmin[i], qmax[i]);
+
+            TIMEIT_ONCE_START
+
+                for (q = qmin[i]; q <= qmax[i]; q+=2)
+                {
+                    acb_dirichlet_group_t G;
+                    acb_dirichlet_char_t chi1, chi2;
+                    acb_t res;
+
+                    acb_dirichlet_group_init(G, q);
+                    acb_dirichlet_char_init(chi1, G);
+                    acb_dirichlet_char_init(chi2, G);
+                    acb_init(res);
+
+                    do {
+                        slong c = 0;
+
+                        acb_dirichlet_char_set(chi2, G, chi1);
+                        do {
+                            func[l](res, G, chi1, chi2, prec[j]);
+                        } while (acb_dirichlet_char_next(chi2, G) >= 0 && c++ < 30);
+
+                    } while (acb_dirichlet_char_next(chi1, G) >= 0);
+
+                    acb_clear(res);
+                    acb_dirichlet_char_clear(chi1);
+                    acb_dirichlet_char_clear(chi2);
+                    acb_dirichlet_group_clear(G);
+                }
+
+            TIMEIT_ONCE_STOP
+
+                if (out == JSON)
+                    flint_printf("}\n");
+                else
+                    flint_printf("\n");
+        }
+
+    }
+
+    }
+    flint_cleanup();
+    return EXIT_SUCCESS;
+}

acb_dirichlet/profile/p-thetanull.c

@@ -0,0 +1,205 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 <string.h>
+#include "acb_dirichlet.h"
+
+static int usage(char *argv[])
+{
+    printf("usage: %s [--progress <step>] [--all] [--odd] [--prec <bits>] [--onlymod] {qmin qmax | --value q a}\n", argv[0]);
+    return 1;
+}
+
+static void
+value(ulong q, ulong a, slong prec, slong digits)
+{
+    acb_dirichlet_group_t G;
+    acb_dirichlet_char_t chi;
+    acb_t res;
+    arb_t one;
+
+    acb_dirichlet_group_init(G, q);
+    acb_dirichlet_char_init(chi, G);
+    acb_dirichlet_char(chi, G, a);
+    acb_init(res);
+    arb_init(one);
+
+    arb_one(one);
+    acb_dirichlet_theta_arb(res, G, chi, one, prec);
+
+    acb_printd(res, digits);
+    flint_printf("\n");
+
+    arb_clear(one);
+    acb_clear(res);
+    acb_dirichlet_group_clear(G);
+    acb_dirichlet_char_clear(chi);
+}
+
+static void
+check_q(ulong q, int odd, slong prec, slong digits, int onlymod)
+{
+    slong s;
+    ulong k, len;
+    acb_dirichlet_group_t G;
+    acb_dirichlet_conrey_t x;
+    acb_dirichlet_group_init(G, q);
+    acb_dirichlet_conrey_init(x, G);
+    acb_ptr theta, z;
+    arb_t z1;
+    arb_ptr t;
+
+    theta = _acb_vec_init(G->phi_q);
+    z =  _acb_vec_init(G->phi_q);
+
+    arb_init(z1);
+    arb_const_pi(z1, prec);
+    arb_div_ui(z1, z1, q, prec);
+    arb_neg(z1, z1);
+    arb_exp(z1, z1, prec);
+
+    /* len = acb_dirichlet_theta_length_d(q, 1., prec); */
+    t = _arb_vec_init(q);
+    acb_dirichlet_arb_quadratic_powers(t, q, z1, prec);
+
+    for (s = 0; s <= odd; s++)
+    {
+        k = 0;
+        acb_dirichlet_conrey_one(x, G);
+        do {
+            acb_set_arb(z + k, t + x->n);
+            k++;
+        } while (acb_dirichlet_conrey_next(x, G) >= 0);
+
+        acb_dirichlet_dft_conrey(theta, z, G, prec);
+
+        /* check zeros */
+        acb_dirichlet_conrey_one(x, G);
+        for (k = 0; k < G->phi_q; k++)
+        {
+            if (acb_contains_zero(theta + k)
+                    && acb_dirichlet_conrey_conductor(G, x) == q
+                    && acb_dirichlet_conrey_parity(G, x) == s)
+            {
+                if (onlymod)
+                {
+                    flint_printf("%wu,%wu\n",q,x->n);
+                }
+                else
+                {
+                    flint_printf("\ntheta null q = %wu, n = %wu\n",q, x->n);
+                    acb_printd(theta + k, digits);
+                    flint_printf("\n");
+                }
+            }
+            acb_dirichlet_conrey_next(x, G);
+        }
+
+        if (odd)
+        {
+            /* change t for odd characters */
+            for (k = 0; k < q; k++)
+                arb_mul_ui(t + k, t + k, k, prec);
+        }
+
+    }
+
+    _arb_vec_clear(t, q);
+    _acb_vec_clear(theta, G->phi_q);
+    _acb_vec_clear(z, G->phi_q);
+    arb_clear(z1);
+    acb_dirichlet_conrey_clear(x);
+    acb_dirichlet_group_clear(G);
+}
+
+
+int main(int argc, char *argv[])
+{
+    int i, all = 0, odd = 0, onlymod = 0, eval = 0;
+    slong prec = 50, digits = 10;
+    slong step = 0;
+    ulong q, qmin, qmax;
+    n_primes_t iter;
+
+    if (argc < 3)
+        return usage(argv);
+
+    for (i = 1; i < argc - 2; i++)
+    {
+        if (!strcmp(argv[i],"--eval"))
+            eval = 1;
+        if (!strcmp(argv[i],"--all"))
+            all = 1;
+        if (!strcmp(argv[i],"--onlymod"))
+            onlymod = 1;
+        else if (!strcmp(argv[i],"--odd"))
+            odd = 1;
+        else if (!strcmp(argv[i],"--progress"))
+        {
+            i++;
+            step = atol(argv[i]);
+        }
+        else if (!strcmp(argv[i],"--prec"))
+        {
+            i++;
+            prec = atol(argv[i]);
+            digits = floor(prec * 0.3);
+        }
+    }
+
+    if (argc < i + 2)
+        return usage(argv);
+
+    qmin = atol(argv[i]);
+    qmax = atol(argv[i+1]);
+
+    if (eval)
+    {
+        value(qmin, qmax, prec, digits);
+    }
+    else
+    {
+        if (all)
+        {
+            for (q = qmin; q <= qmax; q++)
+            {
+                if (q % 4 == 2)
+                    continue;
+                check_q(q, odd, prec, digits, onlymod);
+                if (step && q % step == 0)
+                    flint_printf("[%wu]",q);
+            }
+        }
+        else
+        {
+            ulong p;
+            slong it = 0;
+            /* look for vanishing theta values for prime power moduli */
+            n_primes_init(iter);
+
+            while ((p = n_primes_next(iter)) < qmax)
+            {
+
+                for (q = p; q < qmin; q*= p);
+                for (; q < qmax; q *= p)
+                    check_q(q, odd, prec, digits, onlymod);
+
+                if (step && (it++ % step == 0))
+                    flint_printf("[%wu]",p);
+            }
+
+            n_primes_clear(iter);
+        }
+    }
+
+    flint_cleanup();
+    return EXIT_SUCCESS;
+}

acb_dirichlet/profile/p-vec.c

@@ -0,0 +1,200 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 <string.h>
+#include "acb_dirichlet.h"
+#include "profiler.h"
+
+#define LOG 0
+#define CSV 1
+#define JSON 2
+
+typedef void (*do_f) (ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv);
+
+static void
+do_empty(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
+{
+    return;
+}
+
+static void
+do_dlog_primeloop(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
+{
+    slong k, l;
+
+    for (k = 0; k < nv; k++)
+        v[k] = 0;
+
+    for (l = G->neven; l < G->num; l++)
+    {
+        acb_dirichlet_prime_group_struct P = G->P[l];
+        dlog_vec_loop_add(v, nv, P.g, chi->expo[l], P.pe, P.phi, chi->order);
+    }
+    acb_dirichlet_ui_vec_set_null(v, G, nv);
+}
+
+static void
+do_eratos(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
+{
+	slong k, p, pmax;
+	n_primes_t iter;
+
+	n_primes_init(iter);
+
+	pmax = (nv < G->q) ? nv : G->q;
+	v[1] = 0;
+
+	while ((p = n_primes_next(iter)) < pmax)
+	{
+		if (G->q % p == 0)
+		{
+			for (k = p; k < nv; k += p)
+				v[k] = ACB_DIRICHLET_CHI_NULL;
+		}
+		else
+		{
+			long chip;
+			chip = acb_dirichlet_ui_chi(G, chi, p);
+
+			for (k = p; k < nv; k += p)
+				if (v[k] != -1)
+					 v[k] = nmod_add(v[k], chip, chi->order);
+		}
+	}
+
+	n_primes_clear(iter);
+}
+
+int main(int argc, char *argv[])
+{
+    int out;
+    ulong n, nref, maxq = 5000;
+    ulong * rand;
+    flint_rand_t state;
+    slong r, nr;
+
+    int l, nf = 9;
+    do_f func[9] = { do_empty, acb_dirichlet_ui_chi_vec_loop, do_dlog_primeloop,
+        acb_dirichlet_ui_chi_vec_primeloop, do_eratos,
+        acb_dirichlet_ui_chi_vec,
+        acb_dirichlet_ui_chi_vec,
+        acb_dirichlet_ui_chi_vec,
+        acb_dirichlet_ui_chi_vec };
+    char * name[9] = { "char only", "big loop", "prime loops",
+        "prime dlog_vec", "manual eratos", "default",
+        "precomp 1", "precomp 20", "precomp 100" };
+
+    int i, ni = 5;
+    ulong qmin[5] =  {   2, 1000, 3000, 10000, 100000 };
+    ulong qmax[5] =  { 500, 2000, 5000, 12000, 100500 };
+
+    int j, nj = 3;
+    slong nv[3] = { 50, 300, 2000 };
+
+    nr = 20;
+
+    flint_randinit(state);
+
+    rand = flint_malloc(nr * sizeof(ulong));
+
+    for (r = 0; r < nr; r++)
+        rand[r] = n_randprime(state, 42, 0);
+
+    if (argc < 2)
+        out = LOG;
+    else if (!strcmp(argv[1], "json"))
+        out = JSON;
+    else if (!strcmp(argv[1], "csv"))
+        out = CSV;
+    else if (!strcmp(argv[1], "log"))
+        out = LOG;
+    else
+    {
+        printf("usage: %s [log|csv|json]\n", argv[0]);
+        abort();
+    }
+
+    if (out == CSV)
+        flint_printf("# %-12s, %7s, %7s, %7s, %7s\n","name", "num", "qmin", "qmax", "time");
+
+    for (j = 0; j < nj; j++)
+    {
+
+        ulong * v;
+        v = flint_malloc(nv[j] * sizeof(ulong));
+
+        for (i = 0; i < ni; i++)
+        {
+
+            if (out == LOG)
+                flint_printf("%wu * ui_chi(rand, 1..%wu) for all %wu <= q <= %wu....\n", nr, nv[j], qmin[i], qmax[i]);
+
+            for (l = 0; l < nf; l++)
+            {
+                ulong q;
+
+                /* eratos too slow */
+                if (l == 4 && i > 2)
+                    continue;
+
+                if (out == LOG)
+                    flint_printf("%-14s ...  ", name[l]);
+                else if (out == CSV)
+                    flint_printf("%-12s, %7d, %7d, %7d,   ", name[l], nv[j], qmin[i], qmax[i]);
+                else if (out == JSON)
+                    flint_printf("{ \"name\": \"%s\", \"num\": %d, \"qmin\": %d, \"qmax\": %d, \"time\": ",
+                            name[l], nv[j], qmin[i], qmax[i]);
+
+                TIMEIT_ONCE_START
+
+                for (q = qmin[i]; q <= qmax[i]; q++)
+                {
+                    acb_dirichlet_group_t G;
+                    acb_dirichlet_char_t chi;
+
+                    acb_dirichlet_group_init(G, q);
+                    acb_dirichlet_char_init(chi, G);
+
+                    if (l >= 6)
+                        acb_dirichlet_group_dlog_precompute(G, (l == 6) ? 1 : (l==7) ? 20 : 100);
+
+                    for (r = 0; r < nr; r++)
+                    {
+                        acb_dirichlet_char(chi, G, rand[r] % q);
+                        func[l](v, G, chi, nv[j]);
+                    }
+
+                    if (l >= 6)
+                        acb_dirichlet_group_dlog_clear(G);
+
+                    acb_dirichlet_char_clear(chi);
+                    acb_dirichlet_group_clear(G);
+                }
+
+                TIMEIT_ONCE_STOP
+
+                    if (out == JSON)
+                        flint_printf("}\n");
+                    else
+                        flint_printf("\n");
+            }
+
+        }
+        flint_free(v);
+
+    }
+
+    flint_free(rand);
+    flint_randclear(state);
+
+    flint_cleanup();
+    return EXIT_SUCCESS;
+}

acb_dirichlet/qseries_arb.c

@@ -0,0 +1,42 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+/* sum a_k x^(k^2), assume a[0] = 0 */
+void
+acb_dirichlet_qseries_arb(acb_t res, acb_srcptr a, const arb_t x, slong len, slong prec)
+{
+    slong k;
+    arb_t xk2, dx, x2;
+
+    arb_init(xk2);
+    arb_init(dx);
+    arb_init(x2);
+
+    arb_set(dx, x);
+    arb_set(xk2, dx);
+    arb_mul(x2, dx, dx, prec);
+
+    acb_mul_arb(res, a + 1, xk2, prec);
+
+    /* TODO: reduce prec */
+    for (k = 2; k < len; k++)
+    {
+        arb_mul(dx, dx, x2, prec);
+        arb_mul(xk2, xk2, dx, prec);
+        acb_addmul_arb(res, a + k, xk2, prec);
+    }
+
+    arb_clear(xk2);
+    arb_clear(x2);
+    arb_clear(dx);
+}

acb_dirichlet/qseries_arb_powers.c

@@ -0,0 +1,100 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_qseries_arb_powers_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
+{
+    slong k;
+    arb_t xk2, dx, x2;
+    acb_t zk;
+    arb_init(xk2);
+    arb_init(dx);
+    arb_init(x2);
+    acb_init(zk);
+
+    arb_set(dx, x);
+    arb_set(xk2, dx);
+    arb_mul(x2, dx, dx, prec);
+
+    acb_set_arb(res, xk2);
+
+    /* TODO: reduce prec */
+    for (k = 2; k < len; k++)
+    {
+        arb_mul(dx, dx, x2, prec);
+        arb_mul(xk2, xk2, dx, prec);
+        if (a[k] != DIRICHLET_CHI_NULL)
+        {
+            acb_dirichlet_power(zk, z, a[k], prec);
+            if (parity)
+                acb_mul_si(zk, zk, k, prec);
+            acb_addmul_arb(res, zk, xk2, prec);
+        }
+    }
+
+    arb_clear(xk2);
+    arb_clear(x2);
+    arb_clear(dx);
+    acb_clear(zk);
+}
+
+/* small order, multiply by chi at the end */
+void
+acb_dirichlet_qseries_arb_powers_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
+{
+    slong k;
+    ulong order = z->order;
+    arb_t xk2, kxk2, dx, x2;
+    acb_ptr t;
+    arb_init(xk2);
+    arb_init(dx);
+    arb_init(x2);
+    arb_init(kxk2);
+
+    arb_set(dx, x);
+    arb_set(xk2, x);
+    arb_mul(x2, x, x, prec);
+
+    t = _acb_vec_init(order);
+    _acb_vec_zero(t, order);
+
+    arb_set(acb_realref(t + 0), xk2);
+
+    /* TODO: reduce precision at each step */
+    for (k = 2; k < len; k++)
+    {
+        arb_mul(dx, dx, x2, prec);
+        arb_mul(xk2, xk2, dx, prec);
+        if (a[k] != DIRICHLET_CHI_NULL)
+        {
+           if (parity)
+           {
+               arb_mul_si(kxk2, xk2, k, prec);
+               arb_add(acb_realref(t + a[k]), acb_realref(t + a[k]), kxk2, prec);
+           }
+           else
+           {
+               arb_add(acb_realref(t + a[k]), acb_realref(t + a[k]), xk2, prec);
+           }
+        }
+    }
+
+    /* now Hörner */
+    _acb_poly_evaluate(res, t, order, z->z, prec);
+
+    _acb_vec_clear(t, order);
+    arb_clear(xk2);
+    arb_clear(x2);
+    arb_clear(dx);
+}

acb_dirichlet/root_number.c

@@ -0,0 +1,52 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_root_number_theta(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    arb_t x;
+    acb_t eps;
+
+    arb_init(x);
+    arb_one(x);
+    acb_dirichlet_theta_arb(res, G, chi, x, prec);
+    acb_init(eps);
+    acb_conj(eps, res);
+    acb_div(res, res, eps, prec);
+
+    arb_clear(x);
+    acb_clear(eps);
+}
+
+void
+acb_dirichlet_root_number(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    if (dirichlet_conductor_char(G, chi) < G->q)
+    {
+        flint_printf("root number: need primitive character\n");
+        abort();
+    }
+    else if (G->num > 1)
+    {
+        acb_t iq;
+        acb_init(iq);
+        acb_dirichlet_gauss_sum_order2(iq, G, chi, prec);
+        acb_dirichlet_gauss_sum(res, G, chi, prec);
+        acb_div(res, res, iq, prec);
+        acb_clear(iq);
+    }
+    else
+    {
+        acb_dirichlet_root_number_theta(res, G, chi, prec);
+    }
+}

acb_dirichlet/si_poly_evaluate.c

@@ -0,0 +1,62 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+/* bsgs evaluation */
+void
+acb_dirichlet_si_poly_evaluate(acb_t res, slong * v, slong len, const acb_t z, slong prec)
+{
+    slong k, r, m;
+    acb_t sq;
+    acb_ptr zk;
+
+    if (len < 3)
+    {
+        if (len == 0)
+        {
+            acb_zero(res);
+        }
+        else if (len == 1)
+        {
+            acb_set_si(res, v[0]);
+        }
+        else if (len == 2)
+        {
+            acb_mul_si(res, z, v[1], prec);
+            acb_add_si(res, res, v[0], prec);
+        }
+        return;
+    }
+
+    m = n_sqrt(len) + 1;
+
+    zk = _acb_vec_init(m + 1);
+    _acb_vec_set_powers(zk, z, m + 1, prec);
+
+    acb_init(sq);
+    acb_zero(res);
+
+    k = len - 1;
+    r = k % m;
+    for (; k >= 0; r = m - 1)
+    {
+        acb_zero(sq);
+        for (; r >= 0; r--, k--)
+            acb_addmul_si(sq, zk + r, v[k], prec);
+        acb_mul(res, res, zk + m, prec);
+        acb_add(res, res, sq, prec);
+    }
+
+    _acb_vec_clear(zk, m + 1);
+    acb_clear(sq);
+}

acb_dirichlet/test/t-chi.c

@@ -23,13 +23,15 @@ int main()
     for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
     {
         acb_t zn1, zn2, zn1n2, zn1zn2;
-        acb_dirichlet_group_t G;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
         ulong q, m, n1, n2, iter2;
         int res;
 
         q = 1 + n_randint(state, 1000);
 
-        acb_dirichlet_group_init(G, q);
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
         acb_init(zn1);
         acb_init(zn2);
         acb_init(zn1n2);
@@ -42,12 +44,31 @@ int main()
                 m = 1 + n_randint(state, q);
             } while (n_gcd(q, m) != 1);
 
+            dirichlet_char_log(chi, G, m);
+
             n1 = n_randint(state, 1000);
             n2 = n_randint(state, 1000);
 
-            acb_dirichlet_chi(zn1, G, m, n1, 53);
-            acb_dirichlet_chi(zn2, G, m, n2, 53);
-            acb_dirichlet_chi(zn1n2, G, m, n1 * n2, 53);
+            acb_dirichlet_chi(zn1, G, chi, n1, 53);
+            acb_dirichlet_pairing(zn2, G, m, n1, 53);
+            if (!acb_overlaps(zn1, zn2))
+            {
+                flint_printf("FAIL: overlap\n\n");
+                flint_printf("q = %wu\n\n", q);
+                flint_printf("m = %wu\n\n", m);
+                flint_printf("n = %wu\n\n", n1);
+                flint_printf("char = "); acb_printd(zn1, 15); flint_printf("\n\n");
+                flint_printf("pairing = "); acb_printd(zn2, 15); flint_printf("\n\n");
+                dirichlet_char_print(G, chi);
+                dirichlet_char_log(chi, G, m);
+                flint_printf("log(m) = "); dirichlet_char_print(G, chi);
+                dirichlet_char_log(chi, G, n1);
+                flint_printf("log(n1) = "); dirichlet_char_print(G, chi);
+                abort();
+            }
+
+            acb_dirichlet_pairing(zn2, G, m, n2, 53);
+            acb_dirichlet_pairing(zn1n2, G, m, n1 * n2, 53);
             acb_mul(zn1zn2, zn1, zn2, 53);
 
             if (!acb_overlaps(zn1n2, zn1zn2))
@@ -74,7 +95,7 @@ int main()
             {
                 if (n_gcd(q, m) == 1)
                 {
-                    acb_dirichlet_chi(zn2, G, m, n1, 53);
+                    acb_dirichlet_pairing(zn2, G, m, n1, 53);
                     acb_add(zn1, zn1, zn2, 53);
                 }
             }
@@ -96,7 +117,8 @@ int main()
             }
         }
 
-        acb_dirichlet_group_clear(G);
+        dirichlet_group_clear(G);
+        dirichlet_char_clear(chi);
         acb_clear(zn1);
         acb_clear(zn2);
         acb_clear(zn1n2);
@@ -108,4 +130,3 @@ int main()
     flint_printf("PASS\n");
     return EXIT_SUCCESS;
 }
-

acb_dirichlet/test/t-gauss.c

@@ -0,0 +1,81 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 prec = 128;
+    ulong q;
+
+    flint_printf("gauss....");
+    fflush(stdout);
+
+    /* check Gauss sums */
+
+    for (q = 3; q < 250; q ++)
+    {
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+
+        acb_t s1, s2, s3, s4;
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
+
+        acb_init(s1);
+        acb_init(s2);
+        acb_init(s3);
+        acb_init(s4);
+        dirichlet_char_one(chi, G);
+
+        while (1) {
+
+            acb_dirichlet_gauss_sum_naive(s1, G, chi, prec);
+            acb_dirichlet_gauss_sum(s2, G, chi, prec);
+            acb_dirichlet_gauss_sum_factor(s3, G, chi, prec);
+            if (dirichlet_conductor_char(G, chi) == G->q)
+                acb_dirichlet_gauss_sum_theta(s4, G, chi, prec);
+            else
+                acb_set(s4, s1);
+
+            if (!acb_overlaps(s1, s2)
+                    || !acb_overlaps(s1, s3)
+                    || !acb_overlaps(s1, s4))
+            {
+                flint_printf("FAIL: G(chi_%wu(%wu))\n\n", q, chi->n);
+                flint_printf("\nnaive ");
+                acb_printd(s1, 25);
+                flint_printf("\ndefault ");
+                acb_printd(s2, 25);
+                flint_printf("\nfactor ");
+                acb_printd(s3, 25);
+                flint_printf("\ntheta ");
+                acb_printd(s4, 25);
+                abort();
+            }
+
+            if (dirichlet_char_next(chi, G) < 0)
+                break;
+        }
+        acb_clear(s1);
+        acb_clear(s2);
+        acb_clear(s3);
+        acb_clear(s4);
+
+        dirichlet_group_clear(G);
+        dirichlet_char_clear(chi);
+    }
+
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_dirichlet/test/t-jacobi.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 prec = 128;
+    ulong q;
+
+    flint_printf("jacobi....");
+    fflush(stdout);
+
+    /* check Jacobi sums */
+
+    for (q = 29 * 29; q > 1; q = q%2 ? 3*q+1 : q/2)
+    {
+        slong m1, m2;
+        dirichlet_group_t G;
+        dirichlet_char_t chi1, chi2;
+
+        acb_t s1, s2;
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi1, G);
+        dirichlet_char_init(chi2, G);
+
+        acb_init(s1);
+        acb_init(s2);
+
+        dirichlet_char_one(chi1, G);
+
+        for (m1 = 0; m1 < 50; m1++)
+        {
+
+            dirichlet_char_one(chi2, G);
+
+            for (m2 = 0; m2 < 50; m2++)
+            {
+
+                acb_dirichlet_jacobi_sum_naive(s1, G, chi1, chi2, prec);
+                acb_dirichlet_jacobi_sum(s2, G, chi1, chi2, prec);
+
+                if (!acb_overlaps(s1, s2))
+                {
+                    flint_printf("FAIL: J_%wu(%wu,%wu)",
+                            q, chi1->n, chi2->n);
+                    flint_printf("\nnaive ");
+                    acb_printd(s1, 25);
+                    flint_printf("\ndefault ");
+                    acb_printd(s2, 25);
+                    flint_printf("\n");
+                    flint_printf("cond = %wu, %wu, %wu\n",
+                            dirichlet_conductor_char(G, chi1),
+                            dirichlet_conductor_char(G, chi2),
+                            dirichlet_conductor_ui(G, nmod_mul(chi1->n, chi2->n, G->mod))
+                            );
+                    abort();
+                }
+                if (dirichlet_char_next(chi2, G) < 0)
+                    break;
+
+            }
+
+            if (dirichlet_char_next(chi1, G) < 0)
+                break;
+
+        }
+
+        acb_clear(s1);
+        acb_clear(s2);
+        dirichlet_group_clear(G);
+        dirichlet_char_clear(chi1);
+        dirichlet_char_clear(chi2);
+    }
+
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_dirichlet/test/t-l.c

@@ -0,0 +1,145 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+
+#define nq 5
+#define nx 3
+
+int main()
+{
+
+    slong i, j;
+    ulong q[nq] = {3, 5, 61, 91, 800};
+    ulong m[nq] = {2, 4, 11, 2, 3};
+    slong prec = 150;
+
+    acb_ptr x;
+    /* cannot test at s = 1 with hurwitz */
+    const char * x_r[nx] = { "1", "0.5", "0.5" };
+    const char * x_i[nx] = { "1", "0", "6" };
+
+    acb_t ref, res;
+    /*
+    default(realprecision, 54)
+    X = [ 1 + I, 1/2, 1/2 + 6 * I ]
+    C = [Mod(2,3),Mod(4,5),Mod(11,61),Mod(2,91),Mod(3,800)]
+    v = concat([ [lfun(c,x) | x<-X] | c<-C])
+    apply(z->printf("\"%s\",\n",real(z)),v)
+    apply(z->printf("\"%s\",\n",imag(z)),v)
+    */
+    const char * ref_r[nq * nx] = {
+        "0.655527984002548033786648216345221087939439503905627469",
+        "0.480867557696828626181220063235589877776829730832371863",
+        "1.56831301727577320813799211138797101541772722814204757",
+        "0.521271244517346991221550773660594765406476858135844321",
+        "0.231750947504015755883383661760877226427886696409005898",
+        "0.275543455389521803395512886745330595086898302178508437",
+        "0.598221809458540554839300433683735304093606595684903281",
+        "0.489264190003695740292779374874163221990017067040417393",
+        "0.573331076412428980263984182365365715292560207445592018",
+        "0.510279695870740409778738767334484809708615155539404548",
+        "0.635626509594367380604827545000418331455019188562281349",
+        "0.129304857274642475564179442785425797926079767522671163",
+        "1.18088858810025653590356481638012816019876881487868657",
+        "2.17175778983760437737667738885959799183430688287297767",
+        "3.41568550810774629867945639900431994221065497147578087"
+    };
+    const char * ref_i[nq * nx] = {
+        "0.220206044893215842652155131408935133596486560067476474",
+        "0",
+        "-0.969458654385732175077973304161399773236957587792986099",
+        "0.354614573267731219242838516784605303288232150760467423",
+        "0",
+        "-0.995392028773643947872231871832838543767598301887892796",
+        "1.04370497561090171487193145841005574472705644411957863",
+        "-0.108902811943905225853677097712717212629591264759957602",
+        "-0.232114369998608907337769019848201890558327186146689311",
+        "-0.133300066189980774635445078240315148544665020358019145",
+        "0.0119464572932630291870372694406253796888930803905106876",
+        "-0.567660589679294457801153713636532209809112025502518666",
+        "-0.654079942571300523223917775358845549990877148918886474",
+        "0.970337207245832214408380510793679653538607483205616894",
+        "-1.43652482351673593824956935036654893593947145947637807"
+    };
+
+    flint_printf("l....");
+    fflush(stdout);
+
+    x = _acb_vec_init(nx);
+
+    for (j = 0; j < nx; j++)
+    {
+        if (arb_set_str(acb_realref(x + j), x_r[j], prec) ||
+            arb_set_str(acb_imagref(x + j), x_i[j], prec) )
+        {
+            flint_printf("error while setting x[%ld] <- %s+I*%s\n",
+                    j, x_r[j], x_i[j]);
+            abort();
+        }
+
+    }
+
+    acb_init(ref);
+    acb_init(res);
+
+    for (i = 0; i < nq; i++)
+    {
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+
+        dirichlet_group_init(G, q[i]);
+        dirichlet_char_init(chi, G);
+
+        dirichlet_char_log(chi, G, m[i]);
+
+        for (j = 0; j < nx; j++)
+        {
+
+            if (arb_set_str(acb_realref(ref), ref_r[i * nx + j], prec - 10) ||
+                    arb_set_str(acb_imagref(ref), ref_i[i * nx + j], prec - 10) )
+            {
+                flint_printf("error while setting ref <- %s+I*%s\n",
+                        ref_r[i * nx + j], ref_i[i * nx + j]);
+                abort();
+            }
+
+            acb_dirichlet_l_hurwitz(res, x + j, G, chi, prec + 10);
+
+            if (!acb_contains(ref, res))
+            {
+                    flint_printf("FAIL:\n\n");
+                    flint_printf("q = %wu\n", q[i]);
+                    flint_printf("m = %wu\n", m[i]);
+                    flint_printf("x = ");
+                    acb_printd(x, 54);
+                    flint_printf("\nref = ");
+                    acb_printd(ref, 54);
+                    flint_printf("\nl(chi,x) = ");
+                    acb_printd(res, 54);
+                    flint_printf("\n\n");
+                    abort();
+            }
+
+        }
+        dirichlet_char_clear(chi);
+        dirichlet_group_clear(G);
+
+    }
+    acb_clear(ref);
+    acb_clear(res);
+    _acb_vec_clear(x, nx);
+
+
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_dirichlet/test/t-l_euler_product.c

@@ -0,0 +1,85 @@
+/*
+    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"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("l_euler_product....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 500 * arb_test_multiplier(); iter++)
+    {
+        acb_t s, t, u;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+        ulong q, k;
+        slong prec;
+
+        acb_init(s);
+        acb_init(t);
+        acb_init(u);
+
+        q = 1 + n_randint(state, 50);
+        prec = 2 + n_randint(state, 500);
+        k = n_randint(state, n_euler_phi(q));
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
+        dirichlet_char_index(chi, G, k);
+
+        if (n_randint(state, 2))
+        {
+            acb_set_ui(s, 2 + n_randint(state, prec + 50));
+        }
+        else
+        {
+            acb_randtest(s, state, 2 + n_randint(state, 200), 2);
+            arb_abs(acb_realref(s), acb_realref(s));
+            if (n_randint(state, 2))
+                acb_add_ui(s, s, n_randtest(state), prec);
+        }
+
+        if (n_randint(state, 2))
+            acb_dirichlet_l_hurwitz(t, s, G, chi, prec);
+        else
+            acb_dirichlet_l_euler_product(t, s, G, chi, prec * 1.5);
+
+        acb_dirichlet_l_euler_product(u, s, G, chi, prec);
+
+        if (!acb_overlaps(t, u))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("iter = %ld  q = %lu  k = %lu  prec = %ld\n\n", iter, q, k, prec);
+            flint_printf("s = "); acb_printn(s, 100, 0); flint_printf("\n\n");
+            flint_printf("t = "); acb_printn(t, 100, 0); flint_printf("\n\n");
+            flint_printf("u = "); acb_printn(u, 100, 0); flint_printf("\n\n");
+            abort();
+        }
+
+        dirichlet_char_clear(chi);
+        dirichlet_group_clear(G);
+        acb_clear(s);
+        acb_clear(t);
+        acb_clear(u);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_dirichlet/test/t-l_jet.c

@@ -0,0 +1,201 @@
+/*
+    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"
+
+/* Laurent expansions at s = 1 of first 10 principal L-functions */
+
+/* with mpmath:
+chis = [[1],[0,1],[0,1,1],[0,1,0,1],[0,1,1,1,1],[0,1,0,0,0,1],[0,1,1,1,1,1,1],
+ [0,1,0,1,0,1,0,1],[0,1,1,0,1,1,0,1,1],[0,1,0,1,0,0,0,1,0,1]]
+
+mp.dps = 40
+for chi in chis:
+    phi = chi.count(1); q = len(chi)
+    L = lambda s: dirichlet(s, chi) - phi/((s-1)*q)
+    c0 = taylor(L, 1, 0, method="quad")
+    c1 = taylor(L, 1, 5, singular=True)[1:]
+    for c in c0 + c1:
+        print nstr(c, 20) + ",",
+    print
+*/
+
+#define TESTQ 10
+#define TESTLEN 6
+
+static const double laurent_data[TESTQ][TESTLEN] = {
+    {0.57721566490153286061, 0.072815845483676724861, -0.0048451815964361592423,
+        -0.00034230573671722431103, 0.000096890419394470835728, -6.6110318108421891813e-6},
+    {0.63518142273073908501, 0.11634237461305384831, -0.018765738937942729408,
+        0.00061334298434914532242, 0.00042338142025747308027, -0.00010545096447379519004},
+    {0.7510145394903918042, 0.058764477744540050414, -0.019011359100973296683,
+        0.0056382252365739175151, -0.0009550480622176659462, 0.000021808301216554848718},
+    {0.63518142273073908501, 0.11634237461305384831, -0.018765738937942729408,
+        0.00061334298434914532242, 0.00042338142025747308027, -0.00010545096447379519004},
+    {0.78366011440804636341, -0.014977808062405260803, 0.0090104707969118845102,
+        0.003603799084856807634, -0.0029351216034181476022, 0.00093077685173004747355},
+    {0.60655632993184433857, 0.2095885418562151802, -0.060844893711330538429,
+        0.0068080382961291386117, 0.0022236616427578346453, -0.0013581825996235430782},
+    {0.77274344835207292411, -0.047596894381510269689, 0.035406039531261788462,
+        -0.0054159870134630085898, -0.0019749752308692423114, 0.0014492998471928196325},
+    {0.63518142273073908501, 0.11634237461305384831, -0.018765738937942729408,
+        0.00061334298434914532242, 0.00042338142025747308027, -0.00010545096447379519004},
+    {0.7510145394903918042, 0.058764477744540050414, -0.019011359100973296683,
+        0.0056382252365739175151, -0.0009550480622176659462, 0.000021808301216554848718},
+    {0.66908892942800130547, 0.16801639259476784034, -0.072611999814034642781,
+        0.024624650443138705595, -0.004951850872731033514, -0.00020178815459414925709}
+};
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("l_jet....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* test Laurent series at s = 1 */
+    {
+        acb_t s, t;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+        acb_ptr vec;
+        ulong q;
+        slong i;
+
+        acb_init(s);
+        acb_init(t);
+        vec = _acb_vec_init(TESTLEN);
+
+        acb_one(s);
+        for (q = 1; q <= TESTQ; q++)
+        {
+            dirichlet_group_init(G, q);
+            dirichlet_char_init(chi, G);
+
+            acb_dirichlet_l_jet(vec, s, G, chi, 1, TESTLEN, 100);
+
+            for (i = 0; i < TESTLEN; i++)
+            {
+                acb_set_d(t, laurent_data[q - 1][i]);
+                mag_set_d(arb_radref(acb_realref(t)),
+                    fabs(laurent_data[q - 1][i]) * 1e-14);
+
+                if (!acb_overlaps(vec + i, t))
+                {
+                    flint_printf("FAIL: Laurent series\n\n");
+                    flint_printf("q = %wu  i = %wd\n\n", q, i);
+                    flint_printf("r1 = "); acb_printn(vec + i, 50, 0); flint_printf("\n\n");
+                    flint_printf("r2 = "); acb_printn(t, 50, 0); flint_printf("\n\n");
+                    abort();
+                }
+            }
+
+            dirichlet_char_clear(chi);
+            dirichlet_group_clear(G);
+        }
+
+        acb_clear(s);
+        acb_clear(t);
+        _acb_vec_clear(vec, TESTLEN);
+    }
+
+    /* test self-consistency */
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_t s;
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+        acb_ptr vec1, vec2;
+        slong len1, len2;
+        slong prec1, prec2;
+        int deflate1, deflate2;
+        ulong q, k;
+        slong i;
+
+        len1 = n_randint(state, 5);
+        len2 = n_randint(state, 5);
+        prec1 = 2 + n_randint(state, 100);
+        prec2 = 2 + n_randint(state, 100);
+        deflate1 = n_randint(state, 2);
+        deflate2 = n_randint(state, 2);
+        q = 1 + n_randint(state, 20);
+        k = n_randint(state, n_euler_phi(q));
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
+        dirichlet_char_index(chi, G, k);
+
+        acb_init(s);
+        vec1 = _acb_vec_init(len1);
+        vec2 = _acb_vec_init(len2);
+
+        if (n_randint(state, 4) == 0)
+            acb_one(s);
+        else
+            acb_randtest(s, state, 2 + n_randint(state, 200), 2);
+
+        acb_dirichlet_l_jet(vec1, s, G, chi, deflate1, len1, prec1);
+        acb_dirichlet_l_jet(vec2, s, G, chi, deflate2, len2, prec2);
+
+        if (deflate1 != deflate2 && dirichlet_char_is_principal(G, chi))
+        {
+            /* add or subtract phi(q)/((s+x-1)q) */
+            acb_t t, u;
+            acb_init(t);
+            acb_init(u);
+
+            acb_set_ui(t, n_euler_phi(q));
+            acb_div_ui(t, t, q, prec1);
+            acb_sub_ui(u, s, 1, prec1);
+
+            for (i = 0; i < len1; i++)
+            {
+                acb_div(t, t, u, prec1);
+                if (deflate1)
+                    acb_add(vec1 + i, vec1 + i, t, prec1);
+                else
+                    acb_sub(vec1 + i, vec1 + i, t, prec1);
+                acb_neg(t, t);
+            }
+
+            acb_clear(t);
+            acb_clear(u);
+        }
+
+        for (i = 0; i < FLINT_MIN(len1, len2); i++)
+        {
+            if (!acb_overlaps(vec1 + i, vec2 + i))
+            {
+                flint_printf("FAIL: overlap\n\n");
+                flint_printf("iter = %wd  q = %wu  k = %wu  i = %wd\n\n", iter, q, k, i);
+                flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n");
+                flint_printf("r1 = "); acb_printn(vec1 + i, 50, 0); flint_printf("\n\n");
+                flint_printf("r2 = "); acb_printn(vec2 + i, 50, 0); flint_printf("\n\n");
+                abort();
+            }
+        }
+
+        dirichlet_char_clear(chi);
+        dirichlet_group_clear(G);
+        acb_clear(s);
+        _acb_vec_clear(vec1, len1);
+        _acb_vec_clear(vec2, len2);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_dirichlet/test/t-thetanull.c

@@ -0,0 +1,119 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 prec = 64;
+    ulong q;
+
+    flint_printf("thetanull....");
+    fflush(stdout);
+
+    /* check the only theta functions
+     *   theta(chi) = sum chi(k)* k^odd * exp(-Pi * k^2 / q)
+     * vanishing at 1 correspond to two specific
+     * characters of moduli 300 and 600 + conjugates
+     */
+
+    for (q = 3; q < 1000; q ++)
+    {
+        dirichlet_group_t G;
+        dirichlet_char_t chi;
+        ulong * v, nv, k;
+
+        acb_t sum;
+        acb_ptr z;
+
+        arb_t eq;
+        arb_ptr t, kt, tt;
+
+        if (q % 4 == 2)
+            /* no primitive character mod q */
+            continue;
+
+        dirichlet_group_init(G, q);
+        dirichlet_char_init(chi, G);
+
+        z = _acb_vec_init(G->expo);
+        _acb_vec_nth_roots(z, G->expo, prec);
+
+        nv = acb_dirichlet_theta_length_d(q, 1, prec);
+        v = flint_malloc(nv * sizeof(ulong));
+
+        arb_init(eq);
+        arb_const_pi(eq, prec);
+        arb_div_ui(eq, eq, q, prec);
+        arb_neg(eq, eq);
+        arb_exp(eq, eq, prec);
+
+        t = _arb_vec_init(nv);
+        acb_dirichlet_arb_quadratic_powers(t, nv, eq, prec);
+
+        kt = _arb_vec_init(nv);
+        for (k = 1; k < nv; k++)
+            arb_mul_ui(kt + k, t + k, k, prec);
+
+        /* theta function on primitive characters */
+        acb_init(sum);
+        dirichlet_char_first_primitive(chi, G);
+
+        do {
+
+            acb_zero(sum);
+            dirichlet_chi_vec(v, G, chi, nv);
+
+            tt = dirichlet_parity_char(G, chi) ? kt : t;
+
+            for (k = 1; k < nv; k++)
+                if (v[k] != DIRICHLET_CHI_NULL)
+                    acb_addmul_arb(sum, z + v[k], tt + k, prec);
+
+            if ((q == 300 && (chi->n == 71 || chi->n == 131))
+                    || (q == 600 && (chi->n == 11 || chi->n == 491)))
+            {
+                if (!acb_contains_zero(sum))
+                {
+                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n);
+                acb_printd(sum, 10);
+                flint_printf("\n");
+                dirichlet_char_print(G, chi);
+                flint_printf("\n");
+                abort();
+                }
+            }
+            else if (acb_contains_zero(sum))
+            {
+                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n);
+                acb_printd(sum, 10);
+                flint_printf("\n");
+                dirichlet_char_print(G, chi);
+                flint_printf("\n");
+                abort();
+            }
+
+        } while (dirichlet_char_next_primitive(chi, G) >= 0);
+
+        _acb_vec_clear(z, G->expo);
+        _arb_vec_clear(kt, nv);
+        _arb_vec_clear(t, nv);
+        acb_clear(sum);
+        arb_clear(eq);
+        flint_free(v);
+        dirichlet_group_clear(G);
+        dirichlet_char_clear(chi);
+    }
+
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_dirichlet/theta_arb.c

@@ -0,0 +1,104 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include "acb_poly.h"
+
+void
+_acb_dirichlet_theta_arb_smallorder(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
+{
+    ulong order;
+    ulong * a;
+    int parity;
+    acb_dirichlet_powers_t z;
+
+    parity = dirichlet_parity_char(G, chi);
+    order = dirichlet_order_char(G, chi);
+    a = flint_malloc(len * sizeof(ulong));
+    dirichlet_chi_vec_order(a, G, chi, order, len);
+
+    _acb_dirichlet_powers_init(z, order, 0, 0, prec);
+    acb_dirichlet_qseries_arb_powers_smallorder(res, xt, parity, a, z, len, prec);
+    acb_dirichlet_powers_clear(z);
+
+    flint_free(a);
+}
+
+void
+_acb_dirichlet_theta_arb_series(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
+{
+    acb_ptr a;
+    a = _acb_vec_init(len);
+    acb_dirichlet_chi_vec(a, G, chi, len, prec);
+    if (dirichlet_parity_char(G, chi))
+    {
+        slong k;
+        for (k = 2; k < len; k++)
+            acb_mul_si(a + k, a + k, k, prec);
+    }
+    acb_dirichlet_qseries_arb(res, a, xt, len, prec);
+    _acb_vec_clear(a, len);
+}
+
+void
+_acb_dirichlet_theta_arb_naive(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t xt, slong len, slong prec)
+{
+    ulong order, * a;
+    acb_dirichlet_powers_t z;
+    int parity;
+
+    parity = dirichlet_parity_char(G, chi);
+    order = dirichlet_order_char(G, chi);
+    a = flint_malloc(len * sizeof(ulong));
+    dirichlet_chi_vec_order(a, G, chi, order, len);
+
+    acb_dirichlet_powers_init(z, order, len, prec);
+
+    acb_dirichlet_qseries_arb_powers_naive(res, xt, parity, a, z, len, prec);
+
+    acb_dirichlet_powers_clear(z);
+    flint_free(a);
+}
+
+void
+acb_dirichlet_theta_arb(acb_t res, const dirichlet_group_t G, const dirichlet_char_t chi, const arb_t t, slong prec)
+{
+    slong len;
+    ulong order;
+    arb_t xt;
+    mag_t e;
+
+    len = acb_dirichlet_theta_length(G->q, t, prec);
+
+    arb_init(xt);
+    _acb_dirichlet_theta_argument_at_arb(xt, G->q, t, prec);
+
+    mag_init(e);
+    mag_tail_kexpk2_arb(e, xt, len);
+
+    arb_neg(xt, xt);
+    arb_exp(xt, xt, prec);
+
+    /* TODO: tune this limit */
+    order = dirichlet_order_char(G, chi);
+    if (order < 30)
+        _acb_dirichlet_theta_arb_smallorder(res, G, chi, xt, len, prec);
+    else
+        _acb_dirichlet_theta_arb_naive(res, G, chi, xt, len, prec);
+
+    arb_add_error_mag(acb_realref(res), e);
+    arb_add_error_mag(acb_imagref(res), e);
+
+    mag_clear(e);
+
+    acb_mul_2exp_si(res, res, 1);
+    arb_clear(xt);
+}

acb_dirichlet/theta_length.c

@@ -0,0 +1,80 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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"
+#include <math.h>
+#define PI   3.14159265358
+#define LOG2 0.69314718055
+
+ulong
+acb_dirichlet_theta_length_d(ulong q, double t, slong prec)
+{
+    double a, la;
+    a = PI / (double)q * t * t;
+    la = (a < .3) ? -log(2*a*(1-a)) : .8;
+    la = ((double)prec * LOG2 + la) / a;
+    return ceil(sqrt(la)+.5);
+}
+
+ulong
+acb_dirichlet_theta_length(ulong q, const arb_t t, slong prec)
+{
+    double dt;
+    ulong len;
+    arf_t at;
+    arf_init(at);
+    arb_get_lbound_arf(at, t, 53);
+    dt = arf_get_d(at, ARF_RND_DOWN);
+    len = acb_dirichlet_theta_length_d(q, dt, prec);
+    arf_clear(at);
+    return len;
+}
+
+/* bound for sum_{k>n} k*exp(-a k^2) */
+void
+mag_tail_kexpk2_arb(mag_t res, const arb_t a, ulong n)
+{
+    mag_t m;
+    mag_init(m);
+    arb_get_mag_lower(m, a);
+    /* a < 1/4 */
+    if (mag_cmp_2exp_si(m, -2) <= 0)
+    {
+        mag_t c;
+        mag_init(c);
+        mag_mul_ui_lower(res, m, n*n-n+1);  /* todo: possible overflow */
+        mag_expinv(res, res);
+        /* c = 2a(1+2a) */
+        mag_mul_2exp_si(m, m, 1);
+        mag_one(c);
+        mag_add_lower(c, m, c);
+        mag_mul_lower(c, m, c);
+        mag_div(res, res, c);
+        mag_clear(c);
+    }
+    else
+    {
+        mag_mul_ui_lower(res, m, n*n-n-1);  /* todo: possible overflow */
+        mag_expinv(res, res);
+        mag_mul_ui(res, res, 2);
+    }
+    mag_clear(m);
+}
+
+/* a(t) = Pi / q * t^2 */
+void
+_acb_dirichlet_theta_argument_at_arb(arb_t xt, ulong q, const arb_t t, slong prec)
+{
+    arb_const_pi(xt, prec);
+    arb_div_ui(xt, xt, q, prec);
+    arb_mul(xt, xt, t, prec);
+    arb_mul(xt, xt, t, prec);
+}

acb_dirichlet/ui_theta_arb.c

@@ -0,0 +1,25 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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_ui_theta_arb(acb_t res, const dirichlet_group_t G, ulong a, const arb_t t, slong prec)
+{
+    dirichlet_char_t chi;
+
+    dirichlet_char_init(chi, G);
+    dirichlet_char_log(chi, G, a);
+
+    acb_dirichlet_theta_arb(res, G, chi, t, prec);
+
+    dirichlet_char_clear(chi);
+}

acb_hypgeom.h

@@ -105,6 +105,7 @@ int acb_hypgeom_u_use_asymp(const acb_t z, slong prec);
 void acb_hypgeom_m_asymp(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
 void acb_hypgeom_m_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
 void acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
+void acb_hypgeom_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
 
 void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec);
 void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec);

acb_hypgeom/gamma_lower.c

@@ -14,7 +14,7 @@
 void
 acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized, slong prec)
 {
-    acb_t s1, nz, t;
+    acb_t s1, nz, t, u;
 
     if (!acb_is_finite(s) || !acb_is_finite(z))
     {
@@ -25,6 +25,7 @@ acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized
     acb_init(s1);
     acb_init(nz);
     acb_init(t);
+    acb_init(u);
 
     acb_add_ui(s1, s, 1, prec);
     acb_neg(nz, z);
@@ -32,17 +33,17 @@ acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized
     if (regularized == 0)
     {
         /* \gamma(s, z) = s^-1 z^s 1F1(s, 1+s, -z) */
-        acb_hypgeom_m(res, s, s1, nz, 0, prec);
+        acb_hypgeom_m(u, s, s1, nz, 0, prec);
         acb_pow(t, z, s, prec);
-        acb_mul(res, res, t, prec);
-        acb_div(res, res, s, prec);
+        acb_mul(u, u, t, prec);
+        acb_div(res, u, s, prec);
     }
     else if (regularized == 1)
     {
         /* P(s, z) = z^s \gamma^{*}(s, z) */
-        acb_hypgeom_m(res, s, s1, nz, 1, prec);
+        acb_hypgeom_m(u, s, s1, nz, 1, prec);
         acb_pow(t, z, s, prec);
-        acb_mul(res, res, t, prec);
+        acb_mul(res, u, t, prec);
     }
     else if (regularized == 2)
     {
@@ -53,4 +54,6 @@ acb_hypgeom_gamma_lower(acb_t res, const acb_t s, const acb_t z, int regularized
     acb_clear(s1);
     acb_clear(nz);
     acb_clear(t);
+    acb_clear(u);
 }
+

acb_hypgeom/m.c

@@ -311,3 +311,9 @@ acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regula
     acb_set_round(res, res, prec);
 }
 
+void
+acb_hypgeom_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec)
+{
+    acb_hypgeom_m(res, a, b, z, regularized, prec);
+}
+

acb_poly.h

@@ -69,6 +69,8 @@ ACB_POLY_INLINE slong acb_poly_degree(const acb_poly_t poly)
     return poly->length - 1;
 }
 
+slong acb_poly_valuation(const acb_poly_t poly);
+
 ACB_POLY_INLINE int
 acb_poly_is_zero(const acb_poly_t z)
 {
@@ -173,6 +175,10 @@ void acb_poly_set(acb_poly_t dest, const acb_poly_t src);
 
 void acb_poly_set_round(acb_poly_t dest, const acb_poly_t src, slong prec);
 
+void acb_poly_set_trunc(acb_poly_t res, const acb_poly_t poly, slong n);
+
+void acb_poly_set_trunc_round(acb_poly_t res, const acb_poly_t poly, slong n, slong prec);
+
 void acb_poly_set_arb_poly(acb_poly_t poly, const arb_poly_t re);
 
 void acb_poly_set2_arb_poly(acb_poly_t poly, const arb_poly_t re, const arb_poly_t im);
@@ -232,6 +238,12 @@ void _acb_poly_sub(acb_ptr res, acb_srcptr poly1, slong len1,
 void acb_poly_sub(acb_poly_t res, const acb_poly_t poly1,
               const acb_poly_t poly2, slong prec);
 
+void acb_poly_add_series(acb_poly_t res, const acb_poly_t poly1,
+              const acb_poly_t poly2, slong len, slong prec);
+
+void acb_poly_sub_series(acb_poly_t res, const acb_poly_t poly1,
+              const acb_poly_t poly2, slong len, slong prec);
+
 ACB_POLY_INLINE void
 acb_poly_neg(acb_poly_t res, const acb_poly_t poly)
 {

acb_poly/add_series.c

@@ -0,0 +1,32 @@
+/*
+    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_poly.h"
+
+void
+acb_poly_add_series(acb_poly_t res, const acb_poly_t poly1,
+              const acb_poly_t poly2, slong len, slong prec)
+{
+    slong len1, len2;
+
+    len1 = poly1->length;
+    len2 = poly2->length;
+
+    len1 = FLINT_MIN(len1, len);
+    len2 = FLINT_MIN(len2, len);
+    len = FLINT_MAX(len1, len2);
+
+    acb_poly_fit_length(res, len);
+    _acb_poly_add(res->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, prec);
+    _acb_poly_set_length(res, len);
+    _acb_poly_normalise(res);
+}
+

acb_poly/set_trunc.c

@@ -0,0 +1,35 @@
+/*
+    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_poly.h"
+
+
+void
+acb_poly_set_trunc(acb_poly_t res, const acb_poly_t poly, slong n)
+{
+    if (poly == res)
+    {
+        acb_poly_truncate(res, n);
+    }
+    else
+    {
+        slong rlen;
+
+        rlen = FLINT_MIN(n, poly->length);
+        while (rlen > 0 && acb_is_zero(poly->coeffs + rlen - 1))
+            rlen--;
+
+        acb_poly_fit_length(res, rlen);
+        _acb_vec_set(res->coeffs, poly->coeffs, rlen);
+        _acb_poly_set_length(res, rlen);
+    }
+}
+

acb_poly/set_trunc_round.c

@@ -0,0 +1,36 @@
+/*
+    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_poly.h"
+
+
+void
+acb_poly_set_trunc_round(acb_poly_t res, const acb_poly_t poly, slong n, slong prec)
+{
+    if (poly == res)
+    {
+        acb_poly_truncate(res, n);
+        _acb_vec_set_round(res->coeffs, res->coeffs, res->length, prec);
+    }
+    else
+    {
+        slong rlen;
+
+        rlen = FLINT_MIN(n, poly->length);
+        while (rlen > 0 && acb_is_zero(poly->coeffs + rlen - 1))
+            rlen--;
+
+        acb_poly_fit_length(res, rlen);
+        _acb_vec_set_round(res->coeffs, poly->coeffs, rlen, prec);
+        _acb_poly_set_length(res, rlen);
+    }
+}
+

acb_poly/sub_series.c

@@ -0,0 +1,32 @@
+/*
+    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_poly.h"
+
+void
+acb_poly_sub_series(acb_poly_t res, const acb_poly_t poly1,
+              const acb_poly_t poly2, slong len, slong prec)
+{
+    slong len1, len2;
+
+    len1 = poly1->length;
+    len2 = poly2->length;
+
+    len1 = FLINT_MIN(len1, len);
+    len2 = FLINT_MIN(len2, len);
+    len = FLINT_MAX(len1, len2);
+
+    acb_poly_fit_length(res, len);
+    _acb_poly_sub(res->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, prec);
+    _acb_poly_set_length(res, len);
+    _acb_poly_normalise(res);
+}
+

acb_poly/test/t-add_series.c

@@ -0,0 +1,83 @@
+/*
+    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("add_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_poly_t a, b, c, d;
+        slong len, prec;
+
+        acb_poly_init(a);
+        acb_poly_init(b);
+        acb_poly_init(c);
+        acb_poly_init(d);
+
+        acb_poly_randtest(a, state, 1 + n_randint(state, 10), 100, 10);
+        acb_poly_randtest(b, state, 1 + n_randint(state, 10), 100, 10);
+        acb_poly_randtest(c, state, 1 + n_randint(state, 10), 100, 10);
+        acb_poly_randtest(d, state, 1 + n_randint(state, 10), 100, 10);
+        prec = 2 + n_randint(state, 100);
+        len = n_randint(state, 10);
+
+        acb_poly_add_series(c, a, b, len, prec);
+        acb_poly_add(d, a, b, prec);
+        acb_poly_truncate(d, len);
+
+        if (!acb_poly_equal(c, d))
+        {
+            flint_printf("FAIL\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");
+            flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
+            flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_poly_set(d, a);
+        acb_poly_add_series(d, d, b, len, prec);
+        if (!acb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        acb_poly_set(d, b);
+        acb_poly_add_series(d, a, d, len, prec);
+        if (!acb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        acb_poly_clear(a);
+        acb_poly_clear(b);
+        acb_poly_clear(c);
+        acb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_poly/test/t-set_trunc_round.c

@@ -0,0 +1,79 @@
+/*
+    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_poly.h"
+
+int
+main(void)
+{
+    int iter;
+    flint_rand_t state;
+
+    flint_printf("set_trunc_round....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_poly_t a, b, c, d, e;
+        slong n, prec;
+
+        acb_poly_init(a);
+        acb_poly_init(b);
+        acb_poly_init(c);
+        acb_poly_init(d);
+        acb_poly_init(e);
+
+        acb_poly_randtest(a, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        acb_poly_randtest(b, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        acb_poly_randtest(c, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        acb_poly_randtest(d, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        acb_poly_randtest(e, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+
+        n = n_randint(state, 10);
+        prec = 2 + n_randint(state, 200);
+
+        acb_poly_set_trunc(b, a, n);
+        acb_poly_set_round(b, b, prec);
+
+        acb_poly_set_round(c, a, prec);
+        acb_poly_set_trunc(c, c, n);
+
+        acb_poly_set_trunc_round(d, a, n, prec);
+
+        acb_poly_set(e, a);
+        acb_poly_set_trunc_round(e, e, n, prec);
+
+        if (!acb_poly_equal(b, c) || !acb_poly_equal(c, d) || !acb_poly_equal(d, e))
+        {
+            flint_printf("FAIL\n\n");
+            acb_poly_printd(a, 50), flint_printf("\n\n");
+            acb_poly_printd(b, 50), flint_printf("\n\n");
+            acb_poly_printd(c, 50), flint_printf("\n\n");
+            acb_poly_printd(d, 50), flint_printf("\n\n");
+            acb_poly_printd(e, 50), flint_printf("\n\n");
+            abort();
+        }
+
+        acb_poly_clear(a);
+        acb_poly_clear(b);
+        acb_poly_clear(c);
+        acb_poly_clear(d);
+        acb_poly_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return 0;
+}
+

acb_poly/test/t-sub_series.c

@@ -0,0 +1,83 @@
+/*
+    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("sub_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_poly_t a, b, c, d;
+        slong len, prec;
+
+        acb_poly_init(a);
+        acb_poly_init(b);
+        acb_poly_init(c);
+        acb_poly_init(d);
+
+        acb_poly_randtest(a, state, 1 + n_randint(state, 10), 100, 10);
+        acb_poly_randtest(b, state, 1 + n_randint(state, 10), 100, 10);
+        acb_poly_randtest(c, state, 1 + n_randint(state, 10), 100, 10);
+        acb_poly_randtest(d, state, 1 + n_randint(state, 10), 100, 10);
+        prec = 2 + n_randint(state, 100);
+        len = n_randint(state, 10);
+
+        acb_poly_sub_series(c, a, b, len, prec);
+        acb_poly_sub(d, a, b, prec);
+        acb_poly_truncate(d, len);
+
+        if (!acb_poly_equal(c, d))
+        {
+            flint_printf("FAIL\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");
+            flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
+            flint_printf("c = "); acb_poly_printd(c, 15); flint_printf("\n\n");
+            abort();
+        }
+
+        acb_poly_set(d, a);
+        acb_poly_sub_series(d, d, b, len, prec);
+        if (!acb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        acb_poly_set(d, b);
+        acb_poly_sub_series(d, a, d, len, prec);
+        if (!acb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        acb_poly_clear(a);
+        acb_poly_clear(b);
+        acb_poly_clear(c);
+        acb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_poly/valuation.c

@@ -0,0 +1,25 @@
+/*
+    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_poly.h"
+
+slong
+acb_poly_valuation(const acb_poly_t poly)
+{
+    slong i, len = poly->length;
+
+    for (i = 0; i < len; i++)
+        if (!acb_is_zero(poly->coeffs + i))
+            return i;
+
+    return -1;
+}
+

arb.h

@@ -54,12 +54,7 @@ arb_init(arb_t x)
     mag_init(arb_radref(x));
 }
 
-ARB_INLINE void
-arb_clear(arb_t x)
-{
-    arf_clear(arb_midref(x));
-    mag_clear(arb_radref(x));
-}
+void arb_clear(arb_t x);
 
 ARB_INLINE arb_ptr
 _arb_vec_init(slong n)

arb/clear.c

@@ -0,0 +1,28 @@
+/*
+    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 "arb.h"
+
+/* The clear methods are not defined as inlines since this inflates
+   the code size. We inline the content manually here to avoid function
+   call overhead. */
+void
+arb_clear(arb_t x)
+{
+    ARF_DEMOTE(arb_midref(x));
+
+    if (COEFF_IS_MPZ(ARF_EXP(arb_midref(x))))
+        _fmpz_clear_mpz(ARF_EXP(arb_midref(x)));
+
+    if (COEFF_IS_MPZ(MAG_EXP(arb_radref(x))))
+        _fmpz_clear_mpz(MAG_EXP(arb_radref(x)));
+}
+

arb_hypgeom.h

@@ -0,0 +1,103 @@
+/*
+    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/>.
+*/
+
+#ifndef ARB_HYPGEOM_H
+#define ARB_HYPGEOM_H
+
+#include "arb.h"
+#include "arb_poly.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void arb_hypgeom_pfq(arb_t res, arb_srcptr a, slong p, arb_srcptr b, slong q,
+    const arb_t z, int regularized, slong prec);
+
+void arb_hypgeom_0f1(arb_t res, const arb_t a, const arb_t z, int regularized, slong prec);
+void arb_hypgeom_m(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec);
+void arb_hypgeom_1f1(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec);
+void arb_hypgeom_u(arb_t res, const arb_t a, const arb_t b, const arb_t z, slong prec);
+void arb_hypgeom_2f1(arb_t res, const arb_t a, const arb_t b, const arb_t c, const arb_t z, int regularized, slong prec);
+
+void arb_hypgeom_erf(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_erf_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_erf_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_erfc(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_erfc_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_erfc_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_erfi(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_erfi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_erfi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_fresnel(arb_t res1, arb_t res2, const arb_t z, int normalized, slong prec);
+void _arb_hypgeom_fresnel_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, int normalized, slong len, slong prec);
+void arb_hypgeom_fresnel_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, int normalized, slong len, slong prec);
+
+void arb_hypgeom_ei(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_ei_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_ei_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_si(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_si_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_si_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_ci(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_ci_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_ci_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_shi(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_shi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_shi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_chi(arb_t res, const arb_t z, slong prec);
+void _arb_hypgeom_chi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec);
+void arb_hypgeom_chi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec);
+
+void arb_hypgeom_li(arb_t res, const arb_t z, int offset, slong prec);
+void _arb_hypgeom_li_series(arb_ptr g, arb_srcptr h, slong hlen, int offset, slong len, slong prec);
+void arb_hypgeom_li_series(arb_poly_t g, const arb_poly_t h, int offset, slong len, slong prec);
+
+void arb_hypgeom_bessel_j(arb_t res, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_y(arb_t res, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_jy(arb_t res1, arb_t res2, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_i(arb_t res, const arb_t nu, const arb_t z, slong prec);
+void arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec);
+
+void arb_hypgeom_airy(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const arb_t z, slong prec);
+void arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec);
+void arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime,
+    arb_poly_t bi, arb_poly_t bi_prime, const arb_poly_t z, slong len, slong prec);
+void _arb_hypgeom_airy_series(arb_ptr ai, arb_ptr ai_prime,
+    arb_ptr bi, arb_ptr bi_prime, arb_srcptr z, slong zlen, slong len, slong prec);
+
+void arb_hypgeom_expint(arb_t res, const arb_t s, const arb_t z, slong prec);
+
+void arb_hypgeom_gamma_lower(arb_t res, const arb_t s, const arb_t z, int regularized, slong prec);
+void _arb_hypgeom_gamma_lower_series(arb_ptr g, const arb_t s, arb_srcptr h, slong hlen, int regularized, slong n, slong prec);
+void arb_hypgeom_gamma_lower_series(arb_poly_t g, const arb_t s, const arb_poly_t h, int regularized, slong n, slong prec);
+
+void arb_hypgeom_gamma_upper(arb_t res, const arb_t s, const arb_t z, int regularized, slong prec);
+void _arb_hypgeom_gamma_upper_series(arb_ptr g, const arb_t s, arb_srcptr h, slong hlen, int regularized, slong n, slong prec);
+void arb_hypgeom_gamma_upper_series(arb_poly_t g, const arb_t s, const arb_poly_t h, int regularized, slong n, slong prec);
+
+void arb_hypgeom_beta_lower(arb_t res, const arb_t a, const arb_t c, const arb_t z, int regularized, slong prec);
+void arb_hypgeom_beta_lower_series(arb_poly_t res, const arb_t a, const arb_t b, const arb_poly_t z, int regularized, slong len, slong prec);
+void _arb_hypgeom_beta_lower_series(arb_ptr res, const arb_t a, const arb_t b, arb_srcptr z, slong zlen, int regularized, slong len, slong prec);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
+

arb_hypgeom/airy.c

@@ -0,0 +1,46 @@
+/*
+    Copyright (C) 2015 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_hypgeom.h"
+#include "acb_hypgeom.h"
+
+void
+arb_hypgeom_airy(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const arb_t z, slong prec)
+{
+    acb_struct tmp[5];
+
+    if (ai != NULL) acb_init(tmp);
+    if (aip != NULL) acb_init(tmp + 1);
+    if (bi != NULL) acb_init(tmp + 2);
+    if (bip != NULL) acb_init(tmp + 3);
+
+    acb_init(tmp + 4);
+    acb_set_arb(tmp + 4, z);
+
+    acb_hypgeom_airy(ai ? tmp : NULL,
+                     aip ? tmp + 1 : NULL,
+                     bi ? tmp + 2 : NULL,
+                     bip ? tmp + 3 : NULL,
+                     tmp + 4, prec);
+
+    if (ai != NULL) arb_set(ai, acb_realref(tmp));
+    if (aip != NULL) arb_set(aip, acb_realref(tmp + 1));
+    if (bi != NULL) arb_set(bi, acb_realref(tmp + 2));
+    if (bip != NULL) arb_set(bip, acb_realref(tmp + 3));
+
+    if (ai != NULL) acb_clear(tmp);
+    if (aip != NULL) acb_clear(tmp + 1);
+    if (bi != NULL) acb_clear(tmp + 2);
+    if (bip != NULL) acb_clear(tmp + 3);
+
+    acb_clear(tmp + 4);
+}
+

arb_hypgeom/airy_jet.c

@@ -0,0 +1,50 @@
+/*
+    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 "arb_hypgeom.h"
+
+static void
+airy_recurrence(arb_ptr ai, const arb_t z, slong len, slong prec)
+{
+    slong k;
+
+    if (len >= 3)
+    {
+        arb_mul(ai + 2, ai, z, prec);
+        arb_mul_2exp_si(ai + 2, ai + 2, -1);
+    }
+
+    for (k = 3; k < len; k++)
+    {
+        arb_mul(ai + k, ai + k - 2, z, prec);
+        arb_add(ai + k, ai + k, ai + k - 3, prec);
+        arb_div_ui(ai + k, ai + k, k * (k - 1), prec);
+    }
+}
+
+void
+arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec)
+{
+    if (len <= 0)
+        return;
+
+    if (len == 1)
+    {
+        arb_hypgeom_airy(ai, NULL, bi, NULL, z, prec);
+        return;
+    }
+
+    arb_hypgeom_airy(ai, ai ? (ai + 1) : NULL, bi, bi ? (bi + 1) : NULL, z, prec);
+
+    if (ai != NULL) airy_recurrence(ai, z, len, prec);
+    if (bi != NULL) airy_recurrence(bi, z, len, prec);
+}
+

arb_hypgeom/airy_series.c

@@ -0,0 +1,127 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_airy_series(arb_ptr ai, arb_ptr ai_prime,
+    arb_ptr bi, arb_ptr bi_prime, arb_srcptr z, slong zlen, slong len, slong prec)
+{
+    arb_ptr t, u, v;
+    slong tlen = len + ((ai_prime != NULL) || (bi_prime != NULL));
+
+    zlen = FLINT_MIN(zlen, len);
+
+    if (zlen <= 0)
+        return;
+
+    if (zlen == 1)
+    {
+        arb_hypgeom_airy(ai, ai_prime, bi, bi_prime, z, prec);
+        return;
+    }
+
+    t = _arb_vec_init(tlen);
+    u = _arb_vec_init(tlen);
+    v = _arb_vec_init(len);
+
+    arb_hypgeom_airy_jet((ai || ai_prime) ? t : NULL,
+                         (bi || bi_prime) ? u : NULL, z, tlen, prec);
+
+    /* compose with nonconstant part */
+    arb_zero(v);
+    _arb_vec_set(v + 1, z + 1, zlen - 1);
+
+    if (ai != NULL) _arb_poly_compose_series(ai, t, len, v, zlen, len, prec);
+    if (bi != NULL) _arb_poly_compose_series(bi, u, len, v, zlen, len, prec);
+
+    /* todo: use chain rule to avoid compositions for derivatives? */
+    if (ai_prime != NULL)
+    {
+        _arb_poly_derivative(t, t, len + 1, prec);
+        _arb_poly_compose_series(ai_prime, t, len, v, zlen, len, prec);
+    }
+
+    if (bi_prime != NULL)
+    {
+        _arb_poly_derivative(u, u, len + 1, prec);
+        _arb_poly_compose_series(bi_prime, u, len, v, zlen, len, prec);
+    }
+
+    _arb_vec_clear(t, tlen);
+    _arb_vec_clear(u, tlen);
+    _arb_vec_clear(v, len);
+}
+
+void
+arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime,
+    arb_poly_t bi, arb_poly_t bi_prime, const arb_poly_t z, slong len, slong prec)
+{
+    if (len == 0)
+    {
+        if (ai != NULL) arb_poly_zero(ai);
+        if (ai_prime != NULL) arb_poly_zero(ai_prime);
+        if (bi != NULL) arb_poly_zero(bi);
+        if (bi_prime != NULL) arb_poly_zero(bi_prime);
+        return;
+    }
+
+    if (z->length <= 1)
+        len = 1;
+
+    if (ai != NULL) arb_poly_fit_length(ai, len);
+    if (ai_prime != NULL) arb_poly_fit_length(ai_prime, len);
+    if (bi != NULL) arb_poly_fit_length(bi, len);
+    if (bi_prime != NULL) arb_poly_fit_length(bi_prime, len);
+
+    if (z->length == 0)
+    {
+        arb_t t;
+        arb_init(t);
+        _arb_hypgeom_airy_series(
+            ai ? ai->coeffs : NULL, ai_prime ? ai_prime->coeffs : NULL,
+            bi ? bi->coeffs : NULL, bi_prime ? bi_prime->coeffs : NULL,
+            t, 1, len, prec);
+        arb_clear(t);
+    }
+    else
+    {
+        _arb_hypgeom_airy_series(
+            ai ? ai->coeffs : NULL, ai_prime ? ai_prime->coeffs : NULL,
+            bi ? bi->coeffs : NULL, bi_prime ? bi_prime->coeffs : NULL,
+            z->coeffs, z->length, len, prec);
+    }
+
+    if (ai != NULL)
+    {
+        _arb_poly_set_length(ai, len);
+        _arb_poly_normalise(ai);
+    }
+
+    if (ai_prime != NULL)
+    {
+        _arb_poly_set_length(ai_prime, len);
+        _arb_poly_normalise(ai_prime);
+    }
+
+    if (bi != NULL)
+    {
+        _arb_poly_set_length(bi, len);
+        _arb_poly_normalise(bi);
+    }
+
+    if (bi_prime != NULL)
+    {
+        _arb_poly_set_length(bi_prime, len);
+        _arb_poly_normalise(bi_prime);
+    }
+}
+

arb_hypgeom/beta_lower_series.c

@@ -0,0 +1,114 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_beta_lower_series(arb_ptr res,
+    const arb_t a, const arb_t b, arb_srcptr z, slong zlen, int regularized, 
+    slong len, slong prec)
+{
+    arb_ptr t, u, v;
+    arb_t c, d, e;
+
+    zlen = FLINT_MIN(zlen, len);
+
+    if (zlen == 1)
+    {
+        arb_hypgeom_beta_lower(res, a, b, z, regularized, prec);
+        _arb_vec_zero(res + 1, len - 1);
+        return;
+    }
+
+    t = _arb_vec_init(len);
+    u = _arb_vec_init(len);
+    v = _arb_vec_init(zlen - 1);
+    arb_init(c);
+    arb_init(d);
+    arb_init(e);
+
+    arb_hypgeom_beta_lower(d, a, b, z, regularized, prec);
+
+    if (regularized)
+    {
+        /* todo: except in special cases, we already compute a bunch of
+           gamma functions in beta_lower, so we could avoid recomputing them */
+        arb_add(e, a, b, prec);
+        arb_gamma(e, e, prec);
+        arb_rgamma(c, a, prec);
+        arb_mul(e, e, c, prec);
+        arb_rgamma(c, b, prec);
+        arb_mul(e, e, c, prec);
+    }
+
+    /* u = (1-z)^(b-1) */
+    _arb_vec_neg(t, z, zlen);
+    arb_add_ui(t, t, 1, prec);
+    arb_sub_ui(c, b, 1, prec);
+    _arb_poly_pow_arb_series(u, t, FLINT_MIN(zlen, len - 1), c, len - 1, prec);
+
+    /* t = z^(a-1) */
+    arb_sub_ui(c, a, 1, prec);
+    _arb_poly_pow_arb_series(t, z, FLINT_MIN(zlen, len - 1), c, len - 1, prec);
+
+    /* v = z' */
+    _arb_poly_derivative(v, z, zlen, prec);
+
+    _arb_poly_mullow(res, t, len - 1, u, len - 1, len - 1, prec);
+    _arb_poly_mullow(t, res, len - 1, v, zlen - 1, len - 1, prec);
+
+    _arb_poly_integral(res, t, len, prec);
+
+    if (regularized)
+    {
+        _arb_vec_scalar_mul(res, res, len, e, prec);
+    }
+
+    arb_set(res, d);
+
+    _arb_vec_clear(t, len);
+    _arb_vec_clear(u, len);
+    _arb_vec_clear(v, zlen - 1);
+    arb_clear(c);
+    arb_clear(d);
+    arb_clear(e);
+}
+
+void arb_hypgeom_beta_lower_series(arb_poly_t res,
+    const arb_t a, const arb_t b, const arb_poly_t z, int regularized,
+    slong len, slong prec)
+{
+    if (len == 0)
+    {
+        arb_poly_zero(res);
+        return;
+    }
+        
+    arb_poly_fit_length(res, len);
+
+    if (z->length == 0)
+    {
+        arb_t t;
+        arb_init(t);
+        _arb_hypgeom_beta_lower_series(res->coeffs, a, b, t, 1,
+            regularized, len, prec);
+        arb_clear(t);
+    }
+    else
+    {
+        _arb_hypgeom_beta_lower_series(res->coeffs, a, b, z->coeffs,
+            z->length, regularized, len, prec);
+    }
+
+    _arb_poly_set_length(res, len);
+    _arb_poly_normalise(res);
+}
+

arb_hypgeom/chi_series.c

@@ -0,0 +1,80 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_chi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    arb_t c;
+
+    if (!arb_is_positive(h))
+    {
+        _arb_vec_indeterminate(g, len);
+        return;
+    }
+
+    arb_init(c);
+    arb_hypgeom_chi(c, h, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else
+    {
+        arb_ptr t, u, v;
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(len);
+        v = _arb_vec_init(len);
+
+        /* Chi(h(x)) = integral([h'(x) / h(x)] cosh(h(x)) */
+        _arb_poly_cosh_series(t, h, hlen, len - 1, prec);
+        _arb_poly_derivative(u, h, hlen, prec);
+        _arb_poly_mullow(v, t, len - 1, u, len - 1, len - 1, prec);
+        _arb_poly_div_series(u, v, len - 1, h, hlen, len - 1, prec);
+        _arb_poly_integral(g, u, len, prec);
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, len);
+        _arb_vec_clear(v, len);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_chi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    if (hlen == 0)
+    {
+        arb_poly_inv_series(g, h, len, prec);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_chi_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/ci_series.c

@@ -0,0 +1,80 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_ci_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    arb_t c;
+
+    if (!arb_is_positive(h))
+    {
+        _arb_vec_indeterminate(g, len);
+        return;
+    }
+
+    arb_init(c);
+    arb_hypgeom_ci(c, h, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else
+    {
+        arb_ptr t, u, v;
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(len);
+        v = _arb_vec_init(len);
+
+        /* Ci(h(x)) = integral([h'(x) / h(x)] cos(h(x)) */
+        _arb_poly_cos_series(t, h, hlen, len - 1, prec);
+        _arb_poly_derivative(u, h, hlen, prec);
+        _arb_poly_mullow(v, t, len - 1, u, len - 1, len - 1, prec);
+        _arb_poly_div_series(u, v, len - 1, h, hlen, len - 1, prec);
+        _arb_poly_integral(g, u, len, prec);
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, len);
+        _arb_vec_clear(v, len);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_ci_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    if (hlen == 0)
+    {
+        arb_poly_inv_series(g, h, len, prec);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_ci_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/ei_series.c

@@ -0,0 +1,80 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_ei_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    arb_t c;
+
+    if (arb_contains_zero(h))
+    {
+        _arb_vec_indeterminate(g, len);
+        return;
+    }
+
+    arb_init(c);
+    arb_hypgeom_ei(c, h, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else
+    {
+        arb_ptr t, u, v;
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(len);
+        v = _arb_vec_init(len);
+
+        /* Ei(h(x)) = integral([h'(x) / h(x)] exp(h(x)) */
+        _arb_poly_exp_series(t, h, hlen, len - 1, prec);
+        _arb_poly_derivative(u, h, hlen, prec);
+        _arb_poly_mullow(v, t, len - 1, u, len - 1, len - 1, prec);
+        _arb_poly_div_series(u, v, len - 1, h, hlen, len - 1, prec);
+        _arb_poly_integral(g, u, len, prec);
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, len);
+        _arb_vec_clear(v, len);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_ei_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    if (hlen == 0)
+    {
+        arb_poly_inv_series(g, h, len, prec);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_ei_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/erf_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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_erf_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    arb_t c;
+    arb_init(c);
+
+    arb_hypgeom_erf(c, h, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else
+    {
+        arb_ptr t, u;
+        slong ulen;
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(len);
+
+        /* erf(h(x)) = integral(h'(x) exp(-h(x)^2)) * (2/sqrt(pi)) */
+        ulen = FLINT_MIN(len, 2 * hlen - 1);
+
+        _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec);
+        _arb_vec_neg(u, u, ulen);
+        _arb_poly_exp_series(u, u, ulen, len, prec);
+        _arb_poly_derivative(t, h, hlen, prec);
+        _arb_poly_mullow(g, u, len, t, hlen - 1, len, prec);
+        _arb_poly_integral(g, g, len, prec);
+
+        arb_const_sqrt_pi(t, prec);
+        arb_inv(t, t, prec);
+        arb_mul_2exp_si(t, t, 1);
+        _arb_vec_scalar_mul(g, g, len, t, prec);
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, len);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_erf_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (hlen == 0 || len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_erf_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/erfc_series.c

@@ -0,0 +1,82 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_erfc_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    arb_t c;
+    arb_init(c);
+
+    arb_hypgeom_erfc(c, h, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else
+    {
+        arb_ptr t, u;
+        slong ulen;
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(len);
+
+        /* erfc(h(x)) = integral(-h'(x) exp(-h(x)^2)) * (2/sqrt(pi)) */
+        ulen = FLINT_MIN(len, 2 * hlen - 1);
+
+        _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec);
+        _arb_vec_neg(u, u, ulen);
+        _arb_poly_exp_series(u, u, ulen, len, prec);
+        _arb_poly_derivative(t, h, hlen, prec);
+        _arb_poly_mullow(g, u, len, t, hlen - 1, len, prec);
+        _arb_poly_integral(g, g, len, prec);
+
+        arb_const_sqrt_pi(t, prec);
+        arb_inv(t, t, prec);
+        arb_mul_2exp_si(t, t, 1);
+        _arb_vec_scalar_mul(g, g, len, t, prec);
+        _arb_vec_neg(g, g, len);
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, len);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_erfc_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    if (hlen == 0)
+    {
+        arb_poly_one(g);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_erfc_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/erfi_series.c

@@ -0,0 +1,74 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_erfi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    arb_t c;
+    arb_init(c);
+
+    arb_hypgeom_erfi(c, h, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else
+    {
+        arb_ptr t, u;
+        slong ulen;
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(len);
+
+        /* erfi(h(x)) = integral(h'(x) exp(h(x)^2)) * (2/sqrt(pi)) */
+        ulen = FLINT_MIN(len, 2 * hlen - 1);
+
+        _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec);
+        _arb_poly_exp_series(u, u, ulen, len, prec);
+        _arb_poly_derivative(t, h, hlen, prec);
+        _arb_poly_mullow(g, u, len, t, hlen - 1, len, prec);
+        _arb_poly_integral(g, g, len, prec);
+
+        arb_const_sqrt_pi(t, prec);
+        arb_inv(t, t, prec);
+        arb_mul_2exp_si(t, t, 1);
+        _arb_vec_scalar_mul(g, g, len, t, prec);
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, len);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_erfi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (hlen == 0 || len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_erfi_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/fresnel_series.c

@@ -0,0 +1,112 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_fresnel_series(arb_ptr s, arb_ptr c,
+    arb_srcptr h, slong hlen, int normalized, slong len, slong prec)
+{
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        arb_hypgeom_fresnel(s, c, h, normalized, prec);
+
+        if (s != NULL) _arb_vec_zero(s + 1, len - 1);
+        if (c != NULL) _arb_vec_zero(c + 1, len - 1);
+    }
+    else
+    {
+        arb_t s0, c0;
+        arb_ptr t, u, v;
+        slong ulen;
+
+        arb_init(s0);
+        arb_init(c0);
+
+        arb_hypgeom_fresnel((s != NULL) ? s0 : NULL,
+                            (c != NULL) ? c0 : NULL, h, normalized, prec);
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(len);
+        v = _arb_vec_init(len);
+
+        /* normalized: */
+        /* C(h(x)) = integral(h'(x) cos(-(pi/2) h(x)^2)) */
+        /* S(h(x)) = -integral(h'(x) sin(-(pi/2) h(x)^2)) */
+        ulen = FLINT_MIN(len, 2 * hlen - 1);
+
+        _arb_poly_mullow(u, h, hlen, h, hlen, ulen, prec);
+        _arb_vec_neg(u, u, ulen);
+
+        if (normalized)
+        {
+            _arb_vec_scalar_mul_2exp_si(u, u, ulen, -1);
+            _arb_poly_sin_cos_pi_series(u, v, u, ulen, len, prec);
+        }
+        else
+        {
+            _arb_poly_sin_cos_series(u, v, u, ulen, len, prec);
+        }
+
+        _arb_poly_derivative(t, h, hlen, prec);
+
+        if (s != NULL)
+        {
+            _arb_poly_mullow(s, u, len, t, hlen - 1, len, prec);
+            _arb_poly_integral(s, s, len, prec);
+            _arb_vec_neg(s, s, len);
+            arb_swap(s, s0);
+        }
+
+        if (c != NULL)
+        {
+            _arb_poly_mullow(c, v, len, t, hlen - 1, len, prec);
+            _arb_poly_integral(c, c, len, prec);
+            arb_swap(c, c0);
+        }
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, len);
+        _arb_vec_clear(v, len);
+
+        arb_clear(s0);
+        arb_clear(c0);
+    }
+}
+
+void
+arb_hypgeom_fresnel_series(arb_poly_t s, arb_poly_t c,
+    const arb_poly_t h, int normalized, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (hlen == 0 || len == 0)
+    {
+        if (s != NULL) arb_poly_zero(s);
+        if (c != NULL) arb_poly_zero(c);
+        return;
+    }
+
+    if (s != NULL) arb_poly_fit_length(s, len);
+    if (c != NULL) arb_poly_fit_length(c, len);
+
+    _arb_hypgeom_fresnel_series((s != NULL) ? s->coeffs : NULL,
+                                (c != NULL) ? c->coeffs : NULL,
+                        h->coeffs, hlen, normalized, len, prec);
+
+    if (s != NULL) _arb_poly_set_length(s, len);
+    if (c != NULL) _arb_poly_set_length(c, len);
+    if (s != NULL) _arb_poly_normalise(s);
+    if (c != NULL) _arb_poly_normalise(c);
+}
+

arb_hypgeom/gamma_lower_series.c

@@ -0,0 +1,128 @@
+/*
+    Copyright (C) 2016 Arb authors
+
+    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"
+#include "arb_hypgeom.h"
+
+void
+_arb_hypgeom_gamma_lower_series(arb_ptr g, const arb_t s, arb_srcptr h, slong hlen, int regularized, slong n, slong prec)
+{
+    arb_t c;
+
+    hlen = FLINT_MIN(hlen, n);
+
+    if (regularized == 2 && arb_is_int(s) && arb_is_nonpositive(s))
+    {
+        arb_t ns;
+        arb_init(ns);
+        arb_neg(ns, s);
+        if (g == h)
+        {
+            arb_ptr t = _arb_vec_init(hlen);
+            _arb_vec_set(t, h, hlen);
+            _arb_poly_pow_arb_series(g, t, hlen, ns, n, prec);
+            _arb_vec_clear(t, hlen);
+        }
+        else
+        {
+            _arb_poly_pow_arb_series(g, h, hlen, ns, n, prec);
+        }
+        arb_clear(ns);
+        return;
+    }
+
+    arb_init(c);
+    arb_hypgeom_gamma_lower(c, s, h, regularized, prec);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, n - 1);
+    }
+    else
+    {
+        arb_ptr t, u, v;
+        arb_ptr w = NULL;
+
+        t = _arb_vec_init(n);
+        u = _arb_vec_init(n);
+        v = _arb_vec_init(n);
+
+        if (regularized == 2)
+        {
+            w = _arb_vec_init(n);
+            arb_neg(t, s);
+            _arb_poly_pow_arb_series(w, h, hlen, t, n, prec);
+        }
+
+        /* gamma(s, h(x)) = integral(h'(x) h(x)^(s-1) exp(-h(x)) */
+        arb_sub_ui(u, s, 1, prec);
+        _arb_poly_pow_arb_series(t, h, hlen, u, n, prec);
+        _arb_poly_derivative(u, h, hlen, prec);
+        _arb_poly_mullow(v, t, n, u, hlen - 1, n, prec);
+        _arb_vec_neg(t, h, hlen);
+        _arb_poly_exp_series(t, t, hlen, n, prec);
+        _arb_poly_mullow(g, v, n, t, n, n, prec);
+        _arb_poly_integral(g, g, n, prec);
+
+        if (regularized == 1)
+        {
+            arb_rgamma(t, s, prec);
+            _arb_vec_scalar_mul(g, g, n, t, prec);
+        }
+        else if (regularized == 2)
+        {
+            arb_rgamma(t, s, prec);
+            _arb_vec_scalar_mul(g, g, n, t, prec);
+            _arb_vec_set(u, g, n);
+            _arb_poly_mullow(g, u, n, w, n, n, prec);
+            _arb_vec_clear(w, n);
+        }
+
+        _arb_vec_clear(t, n);
+        _arb_vec_clear(u, n);
+        _arb_vec_clear(v, n);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+
+void
+arb_hypgeom_gamma_lower_series(arb_poly_t g, const arb_t s,
+        const arb_poly_t h, int regularized, slong n, slong prec)
+{
+    slong hlen = h->length;
+
+    if (n == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    arb_poly_fit_length(g, n);
+
+    if (hlen == 0)
+    {
+        arb_t t;
+        arb_init(t);
+        _arb_hypgeom_gamma_lower_series(g->coeffs, s, t, 1, regularized, n, prec);
+        arb_clear(t);
+    }
+    else
+    {
+        _arb_hypgeom_gamma_lower_series(
+                g->coeffs, s, h->coeffs, hlen, regularized, n, prec);
+    }
+
+    _arb_poly_set_length(g, n);
+    _arb_poly_normalise(g);
+}

arb_hypgeom/gamma_upper_series.c

@@ -0,0 +1,106 @@
+/*
+    Copyright (C) 2015 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"
+#include "arb_hypgeom.h"
+
+void
+_arb_hypgeom_gamma_upper_series(arb_ptr g, const arb_t s, arb_srcptr h, slong hlen, int regularized, slong n, slong prec)
+{
+    arb_t c;
+    arb_init(c);
+
+    arb_hypgeom_gamma_upper(c, s, h, regularized, prec);
+
+    hlen = FLINT_MIN(hlen, n);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, n - 1);
+    }
+    else
+    {
+        arb_ptr t, u, v;
+        arb_ptr w = NULL;
+
+        t = _arb_vec_init(n);
+        u = _arb_vec_init(n);
+        v = _arb_vec_init(n);
+
+        if (regularized == 2)
+        {
+            w = _arb_vec_init(n);
+            arb_neg(t, s);
+            _arb_poly_pow_arb_series(w, h, hlen, t, n, prec);
+        }
+
+        /* Gamma(s, h(x)) = -integral(h'(x) h(x)^(s-1) exp(-h(x)) */
+        arb_sub_ui(u, s, 1, prec);
+        _arb_poly_pow_arb_series(t, h, hlen, u, n, prec);
+        _arb_poly_derivative(u, h, hlen, prec);
+        _arb_poly_mullow(v, t, n, u, hlen - 1, n, prec);
+        _arb_vec_neg(t, h, hlen);
+        _arb_poly_exp_series(t, t, hlen, n, prec);
+        _arb_poly_mullow(g, v, n, t, n, n, prec);
+        _arb_poly_integral(g, g, n, prec);
+        _arb_vec_neg(g, g, n);
+
+        if (regularized == 1)
+        {
+            arb_rgamma(t, s, prec);
+            _arb_vec_scalar_mul(g, g, n, t, prec);
+        }
+        else if (regularized == 2)
+        {
+            _arb_vec_set(u, g, n);
+            _arb_poly_mullow(g, u, n, w, n, n, prec);
+            _arb_vec_clear(w, n);
+        }
+
+        _arb_vec_clear(t, n);
+        _arb_vec_clear(u, n);
+        _arb_vec_clear(v, n);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_gamma_upper_series(arb_poly_t g, const arb_t s,
+        const arb_poly_t h, int regularized, slong n, slong prec)
+{
+    slong hlen = h->length;
+
+    if (n == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    arb_poly_fit_length(g, n);
+
+    if (hlen == 0)
+    {
+        arb_t t;
+        arb_init(t);
+        _arb_hypgeom_gamma_upper_series(g->coeffs, s, t, 1, regularized, n, prec);
+        arb_clear(t);
+    }
+    else
+    {
+        _arb_hypgeom_gamma_upper_series(
+                g->coeffs, s, h->coeffs, hlen, regularized, n, prec);
+    }
+
+    _arb_poly_set_length(g, n);
+    _arb_poly_normalise(g);
+}

arb_hypgeom/li_series.c

@@ -0,0 +1,82 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_li_series(arb_ptr g, arb_srcptr h, slong hlen, int offset, slong len, slong prec)
+{
+    arb_t c;
+
+    if (!arb_is_positive(h) || arb_contains_si(h, 1))
+    {
+        _arb_vec_indeterminate(g, len);
+        return;
+    }
+
+    arb_init(c);
+    arb_hypgeom_li(c, h, offset, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else if (len == 2)
+    {
+        arb_log(g, h, prec);
+        arb_div(g + 1, h + 1, g, prec);
+    }
+    else
+    {
+        arb_ptr t, u;
+
+        t = _arb_vec_init(len);
+        u = _arb_vec_init(hlen);
+
+        /* li(h(x)) = integral(h'(x) / log(h(x))) */
+        _arb_poly_log_series(t, h, hlen, len - 1, prec);
+        _arb_poly_derivative(u, h, hlen, prec);
+        _arb_poly_div_series(g, u, hlen - 1, t, len - 1, len - 1, prec);
+        _arb_poly_integral(g, g, len, prec);
+
+        _arb_vec_clear(t, len);
+        _arb_vec_clear(u, hlen);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_li_series(arb_poly_t g, const arb_poly_t h, int offset, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    if (hlen == 0)
+    {
+        arb_poly_inv_series(g, h, len, prec);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_li_series(g->coeffs, h->coeffs, hlen, offset, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/shi_series.c

@@ -0,0 +1,61 @@
+/*
+    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 "arb_hypgeom.h"
+#include "acb_hypgeom.h"
+
+/* todo: use a sinch function? */
+void
+_arb_hypgeom_shi_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        arb_hypgeom_shi(g, h, prec);
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else
+    {
+        acb_ptr t;
+        slong i;
+
+        t = _acb_vec_init(len);
+
+        for (i = 0; i < hlen; i++)
+            arb_set(acb_realref(t + i), h + i);
+
+        _acb_hypgeom_shi_series(t, t, hlen, len, prec);
+
+        for (i = 0; i < len; i++)
+            arb_swap(g + i, acb_realref(t + i));
+
+        _acb_vec_clear(t, len);
+    }
+}
+
+void
+arb_hypgeom_shi_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (hlen == 0 || len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_shi_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/si_series.c

@@ -0,0 +1,70 @@
+/*
+    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 "arb_hypgeom.h"
+
+void
+_arb_hypgeom_si_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)
+{
+    arb_t c;
+
+    arb_init(c);
+    arb_hypgeom_si(c, h, prec);
+
+    hlen = FLINT_MIN(hlen, len);
+
+    if (hlen == 1)
+    {
+        _arb_vec_zero(g + 1, len - 1);
+    }
+    else if (len == 2)
+    {
+        arb_sinc(g, h, prec);
+        arb_mul(g + 1, g, h + 1, prec);
+    }
+    else
+    {
+        arb_ptr t, u;
+
+        t = _arb_vec_init(len - 1);
+        u = _arb_vec_init(hlen - 1);
+
+        /* Si(h(x)) = integral(h'(x) sinc(h(x))) */
+        _arb_poly_sinc_series(t, h, hlen, len - 1, prec);
+        _arb_poly_derivative(u, h, hlen, prec);
+        _arb_poly_mullow(g, t, len - 1, u, hlen - 1, len - 1, prec);
+        _arb_poly_integral(g, g, len, prec);
+
+        _arb_vec_clear(t, len - 1);
+        _arb_vec_clear(u, hlen - 1);
+    }
+
+    arb_swap(g, c);
+    arb_clear(c);
+}
+
+void
+arb_hypgeom_si_series(arb_poly_t g, const arb_poly_t h, slong len, slong prec)
+{
+    slong hlen = h->length;
+
+    if (hlen == 0 || len == 0)
+    {
+        arb_poly_zero(g);
+        return;
+    }
+
+    arb_poly_fit_length(g, len);
+    _arb_hypgeom_si_series(g->coeffs, h->coeffs, hlen, len, prec);
+    _arb_poly_set_length(g, len);
+    _arb_poly_normalise(g);
+}
+

arb_hypgeom/wrappers.c

@@ -0,0 +1,462 @@
+/*
+    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 "arb_hypgeom.h"
+#include "acb_hypgeom.h"
+
+void
+arb_hypgeom_erf(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_erf(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_erfc(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_erfc(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_erfi(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_erfi(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_fresnel(arb_t res1, arb_t res2, const arb_t z, int normalized, slong prec)
+{
+    if (!arb_is_finite(z))
+    {
+        if (res1 != NULL) arb_indeterminate(res1);
+        if (res2 != NULL) arb_indeterminate(res2);
+    }
+    else
+    {
+        acb_t t, u;
+        acb_init(t);
+        acb_init(u);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_fresnel(res1 ? t : NULL, res2 ? u : NULL, t, normalized, prec);
+        if (res1 != NULL) arb_swap(res1, acb_realref(t));
+        if (res2 != NULL) arb_swap(res2, acb_realref(u));
+        acb_clear(t);
+        acb_clear(u);
+    }
+}
+
+void
+arb_hypgeom_ei(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_ei(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_si(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_si(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_ci(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z) || !arb_is_positive(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_ci(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_shi(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_shi(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_chi(arb_t res, const arb_t z, slong prec)
+{
+    if (!arb_is_finite(z) || !arb_is_positive(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_chi(t, t, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_li(arb_t res, const arb_t z, int offset, slong prec)
+{
+    if (!arb_is_finite(z) || !arb_is_nonnegative(z))
+    {
+        arb_indeterminate(res);
+    }
+    else
+    {
+        acb_t t;
+        acb_init(t);
+        arb_set(acb_realref(t), z);
+        acb_hypgeom_li(t, t, offset, prec);
+        arb_swap(res, acb_realref(t));
+        acb_clear(t);
+    }
+}
+
+void
+arb_hypgeom_0f1(arb_t res, const arb_t a, const arb_t z, int regularized, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), a);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_0f1(t, t, u, regularized, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_m(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec)
+{
+    acb_t t, u, v;
+    acb_init(t);
+    acb_init(u);
+    acb_init(v);
+    arb_set(acb_realref(t), a);
+    arb_set(acb_realref(u), b);
+    arb_set(acb_realref(v), z);
+    acb_hypgeom_m(t, t, u, v, regularized, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(v);
+}
+
+void
+arb_hypgeom_1f1(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec)
+{
+    arb_hypgeom_m(res, a, b, z, regularized, prec);
+}
+
+void
+arb_hypgeom_u(arb_t res, const arb_t a, const arb_t b, const arb_t z, slong prec)
+{
+    acb_t t, u, v;
+    acb_init(t);
+    acb_init(u);
+    acb_init(v);
+    arb_set(acb_realref(t), a);
+    arb_set(acb_realref(u), b);
+    arb_set(acb_realref(v), z);
+    acb_hypgeom_u(t, t, u, v, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(v);
+}
+
+void
+arb_hypgeom_2f1(arb_t res, const arb_t a, const arb_t b, const arb_t c, const arb_t z, int regularized, slong prec)
+{
+    acb_t t, u, v, w;
+    acb_init(t);
+    acb_init(u);
+    acb_init(v);
+    acb_init(w);
+    arb_set(acb_realref(t), a);
+    arb_set(acb_realref(u), b);
+    arb_set(acb_realref(v), c);
+    arb_set(acb_realref(w), z);
+    acb_hypgeom_2f1(t, t, u, v, w, regularized, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(v);
+    acb_clear(w);
+}
+
+void
+arb_hypgeom_pfq(arb_t res, arb_srcptr a, slong p, arb_srcptr b, slong q, const arb_t z, int regularized, slong prec)
+{
+    acb_ptr t;
+    slong i;
+    t = _acb_vec_init(p + q + 1);
+    for (i = 0; i < p; i++)
+        arb_set(acb_realref(t + i), a + i);
+    for (i = 0; i < q; i++)
+        arb_set(acb_realref(t + p + i), b + i);
+    arb_set(acb_realref(t + p + q), z);
+    acb_hypgeom_pfq(t, t, p, t + p, q, t + p + q, regularized, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    _acb_vec_clear(t, p + q + 1);
+}
+
+void
+arb_hypgeom_bessel_j(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_j(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_bessel_y(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_y(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_bessel_jy(arb_t res1, arb_t res2, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_jy(t, u, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res1, acb_realref(t));
+    else
+        arb_indeterminate(res1);
+    if (acb_is_finite(u) && acb_is_real(u))
+        arb_swap(res2, acb_realref(u));
+    else
+        arb_indeterminate(res2);
+    acb_clear(t);
+    acb_clear(u);
+}
+void
+arb_hypgeom_bessel_i(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_i(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), nu);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_bessel_k(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_expint(arb_t res, const arb_t s, const arb_t z, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), s);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_expint(t, t, u, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_gamma_lower(arb_t res, const arb_t s, const arb_t z, int regularized, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), s);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_gamma_lower(t, t, u, regularized, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_gamma_upper(arb_t res, const arb_t s, const arb_t z, int regularized, slong prec)
+{
+    acb_t t, u;
+    acb_init(t);
+    acb_init(u);
+    arb_set(acb_realref(t), s);
+    arb_set(acb_realref(u), z);
+    acb_hypgeom_gamma_upper(t, t, u, regularized, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+}
+
+void
+arb_hypgeom_beta_lower(arb_t res, const arb_t a, const arb_t b, const arb_t z, int regularized, slong prec)
+{
+    acb_t t, u, v;
+    acb_init(t);
+    acb_init(u);
+    acb_init(v);
+    arb_set(acb_realref(t), a);
+    arb_set(acb_realref(u), b);
+    arb_set(acb_realref(v), z);
+    acb_hypgeom_beta_lower(t, t, u, v, regularized, prec);
+    if (acb_is_finite(t) && acb_is_real(t))
+        arb_swap(res, acb_realref(t));
+    else
+        arb_indeterminate(res);
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(v);
+}
+

arb_poly.h

@@ -65,6 +65,10 @@ void arb_poly_set(arb_poly_t poly, const arb_poly_t src);
 
 void arb_poly_set_round(arb_poly_t poly, const arb_poly_t src, slong prec);
 
+void arb_poly_set_trunc(arb_poly_t res, const arb_poly_t poly, slong n);
+
+void arb_poly_set_trunc_round(arb_poly_t res, const arb_poly_t poly, slong n, slong prec);
+
 /* Basic manipulation */
 
 ARB_POLY_INLINE slong arb_poly_length(const arb_poly_t poly)
@@ -77,6 +81,8 @@ ARB_POLY_INLINE slong arb_poly_degree(const arb_poly_t poly)
     return poly->length - 1;
 }
 
+slong arb_poly_valuation(const arb_poly_t poly);
+
 ARB_POLY_INLINE int
 arb_poly_is_zero(const arb_poly_t z)
 {
@@ -210,6 +216,12 @@ void _arb_poly_sub(arb_ptr res, arb_srcptr poly1, slong len1,
 void arb_poly_sub(arb_poly_t res, const arb_poly_t poly1,
               const arb_poly_t poly2, slong prec);
 
+void arb_poly_add_series(arb_poly_t res, const arb_poly_t poly1,
+              const arb_poly_t poly2, slong len, slong prec);
+
+void arb_poly_sub_series(arb_poly_t res, const arb_poly_t poly1,
+              const arb_poly_t poly2, slong len, slong prec);
+
 ARB_POLY_INLINE void
 arb_poly_neg(arb_poly_t res, const arb_poly_t poly)
 {

arb_poly/add_series.c

@@ -0,0 +1,32 @@
+/*
+    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 "arb_poly.h"
+
+void
+arb_poly_add_series(arb_poly_t res, const arb_poly_t poly1,
+              const arb_poly_t poly2, slong len, slong prec)
+{
+    slong len1, len2;
+
+    len1 = poly1->length;
+    len2 = poly2->length;
+
+    len1 = FLINT_MIN(len1, len);
+    len2 = FLINT_MIN(len2, len);
+    len = FLINT_MAX(len1, len2);
+
+    arb_poly_fit_length(res, len);
+    _arb_poly_add(res->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, prec);
+    _arb_poly_set_length(res, len);
+    _arb_poly_normalise(res);
+}
+

arb_poly/set_trunc.c

@@ -0,0 +1,35 @@
+/*
+    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 "arb_poly.h"
+
+
+void
+arb_poly_set_trunc(arb_poly_t res, const arb_poly_t poly, slong n)
+{
+    if (poly == res)
+    {
+        arb_poly_truncate(res, n);
+    }
+    else
+    {
+        slong rlen;
+
+        rlen = FLINT_MIN(n, poly->length);
+        while (rlen > 0 && arb_is_zero(poly->coeffs + rlen - 1))
+            rlen--;
+
+        arb_poly_fit_length(res, rlen);
+        _arb_vec_set(res->coeffs, poly->coeffs, rlen);
+        _arb_poly_set_length(res, rlen);
+    }
+}
+

arb_poly/set_trunc_round.c

@@ -0,0 +1,36 @@
+/*
+    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 "arb_poly.h"
+
+
+void
+arb_poly_set_trunc_round(arb_poly_t res, const arb_poly_t poly, slong n, slong prec)
+{
+    if (poly == res)
+    {
+        arb_poly_truncate(res, n);
+        _arb_vec_set_round(res->coeffs, res->coeffs, res->length, prec);
+    }
+    else
+    {
+        slong rlen;
+
+        rlen = FLINT_MIN(n, poly->length);
+        while (rlen > 0 && arb_is_zero(poly->coeffs + rlen - 1))
+            rlen--;
+
+        arb_poly_fit_length(res, rlen);
+        _arb_vec_set_round(res->coeffs, poly->coeffs, rlen, prec);
+        _arb_poly_set_length(res, rlen);
+    }
+}
+

arb_poly/sub_series.c

@@ -0,0 +1,32 @@
+/*
+    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 "arb_poly.h"
+
+void
+arb_poly_sub_series(arb_poly_t res, const arb_poly_t poly1,
+              const arb_poly_t poly2, slong len, slong prec)
+{
+    slong len1, len2;
+
+    len1 = poly1->length;
+    len2 = poly2->length;
+
+    len1 = FLINT_MIN(len1, len);
+    len2 = FLINT_MIN(len2, len);
+    len = FLINT_MAX(len1, len2);
+
+    arb_poly_fit_length(res, len);
+    _arb_poly_sub(res->coeffs, poly1->coeffs, len1, poly2->coeffs, len2, prec);
+    _arb_poly_set_length(res, len);
+    _arb_poly_normalise(res);
+}
+

arb_poly/test/t-add_series.c

@@ -0,0 +1,83 @@
+/*
+    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("add_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        arb_poly_t a, b, c, d;
+        slong len, prec;
+
+        arb_poly_init(a);
+        arb_poly_init(b);
+        arb_poly_init(c);
+        arb_poly_init(d);
+
+        arb_poly_randtest(a, state, 1 + n_randint(state, 10), 100, 10);
+        arb_poly_randtest(b, state, 1 + n_randint(state, 10), 100, 10);
+        arb_poly_randtest(c, state, 1 + n_randint(state, 10), 100, 10);
+        arb_poly_randtest(d, state, 1 + n_randint(state, 10), 100, 10);
+        prec = 2 + n_randint(state, 100);
+        len = n_randint(state, 10);
+
+        arb_poly_add_series(c, a, b, len, prec);
+        arb_poly_add(d, a, b, prec);
+        arb_poly_truncate(d, len);
+
+        if (!arb_poly_equal(c, d))
+        {
+            flint_printf("FAIL\n\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");
+            flint_printf("c = "); arb_poly_printd(c, 15); flint_printf("\n\n");
+            abort();
+        }
+
+        arb_poly_set(d, a);
+        arb_poly_add_series(d, d, b, len, prec);
+        if (!arb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        arb_poly_set(d, b);
+        arb_poly_add_series(d, a, d, len, prec);
+        if (!arb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        arb_poly_clear(b);
+        arb_poly_clear(c);
+        arb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arb_poly/test/t-set_trunc_round.c

@@ -0,0 +1,79 @@
+/*
+    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 "arb_poly.h"
+
+int
+main(void)
+{
+    int iter;
+    flint_rand_t state;
+
+    flint_printf("set_trunc_round....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        arb_poly_t a, b, c, d, e;
+        slong n, prec;
+
+        arb_poly_init(a);
+        arb_poly_init(b);
+        arb_poly_init(c);
+        arb_poly_init(d);
+        arb_poly_init(e);
+
+        arb_poly_randtest(a, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        arb_poly_randtest(b, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        arb_poly_randtest(c, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        arb_poly_randtest(d, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+        arb_poly_randtest(e, state, n_randint(state, 10), 2 + n_randint(state, 200), 10);
+
+        n = n_randint(state, 10);
+        prec = 2 + n_randint(state, 200);
+
+        arb_poly_set_trunc(b, a, n);
+        arb_poly_set_round(b, b, prec);
+
+        arb_poly_set_round(c, a, prec);
+        arb_poly_set_trunc(c, c, n);
+
+        arb_poly_set_trunc_round(d, a, n, prec);
+
+        arb_poly_set(e, a);
+        arb_poly_set_trunc_round(e, e, n, prec);
+
+        if (!arb_poly_equal(b, c) || !arb_poly_equal(c, d) || !arb_poly_equal(d, e))
+        {
+            flint_printf("FAIL\n\n");
+            arb_poly_printd(a, 50), flint_printf("\n\n");
+            arb_poly_printd(b, 50), flint_printf("\n\n");
+            arb_poly_printd(c, 50), flint_printf("\n\n");
+            arb_poly_printd(d, 50), flint_printf("\n\n");
+            arb_poly_printd(e, 50), flint_printf("\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        arb_poly_clear(b);
+        arb_poly_clear(c);
+        arb_poly_clear(d);
+        arb_poly_clear(e);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return 0;
+}
+

arb_poly/test/t-sub_series.c

@@ -0,0 +1,83 @@
+/*
+    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("sub_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        arb_poly_t a, b, c, d;
+        slong len, prec;
+
+        arb_poly_init(a);
+        arb_poly_init(b);
+        arb_poly_init(c);
+        arb_poly_init(d);
+
+        arb_poly_randtest(a, state, 1 + n_randint(state, 10), 100, 10);
+        arb_poly_randtest(b, state, 1 + n_randint(state, 10), 100, 10);
+        arb_poly_randtest(c, state, 1 + n_randint(state, 10), 100, 10);
+        arb_poly_randtest(d, state, 1 + n_randint(state, 10), 100, 10);
+        prec = 2 + n_randint(state, 100);
+        len = n_randint(state, 10);
+
+        arb_poly_sub_series(c, a, b, len, prec);
+        arb_poly_sub(d, a, b, prec);
+        arb_poly_truncate(d, len);
+
+        if (!arb_poly_equal(c, d))
+        {
+            flint_printf("FAIL\n\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");
+            flint_printf("c = "); arb_poly_printd(c, 15); flint_printf("\n\n");
+            abort();
+        }
+
+        arb_poly_set(d, a);
+        arb_poly_sub_series(d, d, b, len, prec);
+        if (!arb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 1)\n\n");
+            abort();
+        }
+
+        arb_poly_set(d, b);
+        arb_poly_sub_series(d, a, d, len, prec);
+        if (!arb_poly_equal(d, c))
+        {
+            flint_printf("FAIL (aliasing 2)\n\n");
+            abort();
+        }
+
+        arb_poly_clear(a);
+        arb_poly_clear(b);
+        arb_poly_clear(c);
+        arb_poly_clear(d);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

arb_poly/valuation.c

@@ -0,0 +1,25 @@
+/*
+    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 "arb_poly.h"
+
+slong
+arb_poly_valuation(const arb_poly_t poly)
+{
+    slong i, len = poly->length;
+
+    for (i = 0; i < len; i++)
+        if (!arb_is_zero(poly->coeffs + i))
+            return i;
+
+    return -1;
+}
+

arf.h

@@ -228,12 +228,7 @@ arf_init(arf_t x)
     ARF_XSIZE(x) = 0;
 }
 
-ARF_INLINE void
-arf_clear(arf_t x)
-{
-    fmpz_clear(ARF_EXPREF(x));
-    ARF_DEMOTE(x);
-}
+void arf_clear(arf_t x);
 
 ARF_INLINE void
 arf_zero(arf_t x)
@@ -348,6 +343,12 @@ int arf_cmp(const arf_t x, const arf_t y);
 
 int arf_cmpabs(const arf_t x, const arf_t y);
 
+int arf_cmp_si(const arf_t x, slong y);
+
+int arf_cmp_ui(const arf_t x, ulong y);
+
+int arf_cmp_d(const arf_t x, double y);
+
 ARF_INLINE void
 arf_swap(arf_t y, arf_t x)
 {

arf/clear.c

@@ -0,0 +1,23 @@
+/*
+    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 "arf.h"
+
+void
+arf_clear(arf_t x)
+{
+    ARF_DEMOTE(x);
+
+    /* fmpz_clear(), but avoids a redundant store */
+    if (COEFF_IS_MPZ(ARF_EXP(x)))
+        _fmpz_clear_mpz(ARF_EXP(x));
+}
+

arf/cmp.c

@@ -59,3 +59,25 @@ arf_cmp(const arf_t x, const arf_t y)
     return 0;
 }
 
+int arf_cmp_si(const arf_t x, slong y)
+{
+    arf_t t;
+    arf_init_set_si(t, y); /* no need to free */
+    return arf_cmp(x, t);
+}
+
+int arf_cmp_ui(const arf_t x, ulong y)
+{
+    arf_t t;
+    arf_init_set_ui(t, y); /* no need to free */
+    return arf_cmp(x, t);
+}
+
+int arf_cmp_d(const arf_t x, double y)
+{
+    arf_t t;
+    arf_init(t); /* no need to free */
+    arf_set_d(t, y);
+    return arf_cmp(x, t);
+}
+

dirichlet.h

@@ -0,0 +1,161 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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/>.
+*/
+
+#ifndef DIRICHLET_H
+#define DIRICHLET_H
+
+#ifdef DIRICHLET_INLINES_C
+#define DIRICHLET_INLINE
+#else
+#define DIRICHLET_INLINE static __inline__
+#endif
+
+#include "dlog.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#define MAX_FACTORS 15
+
+/* should this dlog pointer be in the prime or the global group? */
+typedef struct
+{
+    ulong p;    /* underlying prime */
+    int e;      /* exponent */
+    nmod_t pe;  /* modulus */
+    nmod_t phi;  /* phi(p^e) */
+    ulong g;    /* conrey generator */
+    dlog_precomp_struct * dlog;  /* precomputed data for discrete log mod p^e */
+}
+dirichlet_prime_group_struct;
+
+typedef struct
+{
+    ulong q;                /* modulus */
+    ulong q_even;           /* even part of modulus */
+    nmod_t mod;             /* modulus with precomputed inverse */
+    ulong rad_q;            /* radical = product of odd primes */
+    ulong phi_q;            /* phi(q) = group size */
+    slong neven;            /* number of even components (in 0,1,2)*/
+    slong num;              /* number of prime components (even + odd) */
+    ulong expo;             /* exponent = largest order in G */
+    dirichlet_prime_group_struct * P;
+    ulong * generators;     /* generators lifted mod q */
+    ulong * PHI;            /* PHI(k) = expo / phi(k)        */
+}
+dirichlet_group_struct;
+
+typedef dirichlet_group_struct dirichlet_group_t[1];
+
+DIRICHLET_INLINE ulong
+dirichlet_group_size(const dirichlet_group_t G)
+{
+    return G->phi_q;
+}
+
+void dirichlet_group_init(dirichlet_group_t G, ulong q);
+void dirichlet_subgroup_init(dirichlet_group_t H, const dirichlet_group_t G, ulong h);
+void dirichlet_group_clear(dirichlet_group_t G);
+void dirichlet_group_dlog_precompute(dirichlet_group_t G, ulong num);
+void dirichlet_group_dlog_clear(dirichlet_group_t G);
+
+/* properties of elements without log */
+
+ulong dirichlet_number_primitive(const dirichlet_group_t G);
+ulong dirichlet_conductor_ui(const dirichlet_group_t G, ulong a);
+int dirichlet_parity_ui(const dirichlet_group_t G, ulong a);
+ulong dirichlet_order_ui(const dirichlet_group_t G, ulong a);
+
+/* characters, keep both number and log */
+typedef struct
+{
+    ulong n;           /* number */
+    ulong * log;       /* s.t. prod generators[k]^log[k] = number */
+}
+dirichlet_char_struct;
+
+typedef dirichlet_char_struct dirichlet_char_t[1];
+
+void dirichlet_char_init(dirichlet_char_t x, const dirichlet_group_t G);
+void dirichlet_char_clear(dirichlet_char_t x);
+void dirichlet_char_print(const dirichlet_group_t G, const dirichlet_char_t x);
+
+DIRICHLET_INLINE void
+dirichlet_char_set(dirichlet_char_t x, const dirichlet_group_t G, const dirichlet_char_t y)
+{
+    slong k;
+    x->n = y->n;
+    for (k = 0; k < G->num; k++)
+        x->log[k] = y->log[k];
+}
+
+DIRICHLET_INLINE int
+dirichlet_char_eq(const dirichlet_char_t x, const dirichlet_char_t y)
+{
+    return (x->n == y->n);
+}
+
+int dirichlet_char_eq_deep(const dirichlet_group_t G, const dirichlet_char_t x, const dirichlet_char_t y);
+int dirichlet_parity_char(const dirichlet_group_t G, const dirichlet_char_t x);
+ulong dirichlet_conductor_char(const dirichlet_group_t G, const dirichlet_char_t x);
+ulong dirichlet_order_char(const dirichlet_group_t G, const dirichlet_char_t x);
+
+void dirichlet_char_log(dirichlet_char_t x, const dirichlet_group_t G, ulong m);
+ulong dirichlet_char_exp(dirichlet_char_t x, const dirichlet_group_t G);
+
+void dirichlet_char_index(dirichlet_char_t x, const dirichlet_group_t G, ulong j);
+ulong dirichlet_index_char(const dirichlet_group_t G, const dirichlet_char_t x);
+
+void dirichlet_char_one(dirichlet_char_t x, const dirichlet_group_t G);
+void dirichlet_char_first_primitive(dirichlet_char_t x, const dirichlet_group_t G);
+
+int dirichlet_char_next(dirichlet_char_t x, const dirichlet_group_t G);
+int dirichlet_char_next_primitive(dirichlet_char_t x, const dirichlet_group_t G);
+
+void dirichlet_char_mul(dirichlet_char_t c, const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b);
+void dirichlet_char_pow(dirichlet_char_t c, const dirichlet_group_t G, const dirichlet_char_t a, ulong n);
+
+void dirichlet_char_lower(dirichlet_char_t y, const dirichlet_group_t H, const dirichlet_char_t x, const dirichlet_group_t G);
+void dirichlet_char_lift(dirichlet_char_t y, const dirichlet_group_t G, const dirichlet_char_t x, const dirichlet_group_t H);
+
+#define DIRICHLET_CHI_NULL UWORD_MAX
+
+ulong dirichlet_pairing(const dirichlet_group_t G, ulong m, ulong n);
+ulong dirichlet_pairing_char(const dirichlet_group_t G, const dirichlet_char_t a, const dirichlet_char_t b);
+
+DIRICHLET_INLINE int
+dirichlet_char_is_principal(const dirichlet_group_t G, const dirichlet_char_t chi)
+{
+    return (chi->n == 1);
+}
+
+DIRICHLET_INLINE int
+dirichlet_char_is_real(const dirichlet_group_t G, const dirichlet_char_t chi)
+{
+    return G->q <= 4 || (nmod_mul(chi->n, chi->n, G->mod) == 1);
+}
+
+ulong dirichlet_chi(const dirichlet_group_t G, const dirichlet_char_t chi, ulong n);
+
+void dirichlet_vec_set_null(ulong *v, const dirichlet_group_t G, slong nv);
+void dirichlet_chi_vec_loop(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv);
+void dirichlet_chi_vec_primeloop(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv);
+void dirichlet_chi_vec(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv);
+void dirichlet_chi_vec_loop_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv);
+void dirichlet_chi_vec_primeloop_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv);
+void dirichlet_chi_vec_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

dirichlet/chi.c

@@ -0,0 +1,33 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 "dirichlet.h"
+
+ulong
+dirichlet_chi(const dirichlet_group_t G, const dirichlet_char_t chi, ulong n)
+{
+    if (n_gcd(G->q, n) > 1)
+    {
+        return DIRICHLET_CHI_NULL;
+    }
+    else
+    {
+        ulong v;
+        dirichlet_char_t x;
+        dirichlet_char_init(x, G);
+        dirichlet_char_log(x, G, n);
+
+        v = dirichlet_pairing_char(G, chi, x);
+
+        dirichlet_char_clear(x);
+        return v;
+    }
+}

dirichlet/chi_vec.c

@@ -0,0 +1,30 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    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 "dirichlet.h"
+
+void
+dirichlet_chi_vec(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, slong nv)
+{
+    if (2 * nv > G->phi_q)
+        dirichlet_chi_vec_loop(v, G, chi, nv);
+    else
+        dirichlet_chi_vec_primeloop(v, G, chi, nv);
+}
+
+void
+dirichlet_chi_vec_order(ulong *v, const dirichlet_group_t G, const dirichlet_char_t chi, ulong order, slong nv)
+{
+    if (2 * nv > G->phi_q)
+        dirichlet_chi_vec_loop_order(v, G, chi, order, nv);
+    else
+        dirichlet_chi_vec_primeloop_order(v, G, chi, order, nv);
+}