File: hurwitz_precomp_eval.c

package info (click to toggle)
flint-arb 1%3A2.19.0-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 13,028 kB
  • sloc: ansic: 177,109; sh: 553; makefile: 288; python: 268
file content (77 lines) | stat: -rw-r--r-- 1,896 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
/*
    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"
#include "acb_poly.h"

void
acb_dirichlet_hurwitz_precomp_eval(acb_t res,
    const acb_dirichlet_hurwitz_precomp_t pre, ulong p, ulong q, slong prec)
{
    slong i;
    acb_t a, t;

    if (p > q)
    {
        flint_printf("hurwitz_precomp_eval: require p <= n\n");
        flint_abort();
    }

    if (pre->A == 0)
    {
        acb_init(a);
        acb_set_ui(a, p);
        acb_div_ui(a, a, q, prec);

        if (pre->deflate == 0)
            acb_hurwitz_zeta(res, &pre->s, a, prec);
        else
            _acb_poly_zeta_cpx_series(res, &pre->s, a, 1, 1, prec);

        acb_clear(a);
        return;
    }

    acb_init(a);
    acb_init(t);

    /* todo: avoid integer overflows below */
    if (p == q)
        i = pre->N - 1;
    else
        i = (pre->N * p) / q;
    acb_set_si(a, 2 * pre->N * p - q * (2 * i + 1));
    acb_div_ui(a, a, 2 * q * pre->N, prec);

    /* compute zeta(s,A+p/q) */
    _acb_poly_evaluate(res, pre->coeffs + i * pre->K, pre->K, a, prec);

    /* error bound */
    if (acb_is_real(&pre->s))
        arb_add_error_mag(acb_realref(res), &pre->err);
    else
        acb_add_error_mag(res, &pre->err);

    /* zeta(s,p/q) = (p/q)^-s + ... + (p/q+A-1)^-s zeta(s,A+p/q) */
    for (i = 0; i < pre->A; i++)
    {
        acb_set_ui(a, p);
        acb_div_ui(a, a, q, prec);
        acb_add_ui(a, a, i, prec);
        acb_neg(t, &pre->s);
        acb_pow(a, a, t, prec);
        acb_add(res, res, a, prec);
    }

    acb_clear(a);
    acb_clear(t);
}