flint-sys 0.9.0

Bindings to the FLINT C library
Documentation 
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

/*
    Copyright (C) 2021 Fredrik Johansson

    This file is part of FLINT.

    FLINT 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 3 of the License, or
    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
*/

#include <math.h>
#include "ulong_extras.h"
#include "acb.h"
#include "arb_hypgeom.h"
#include "acb_hypgeom.h"

void
acb_hypgeom_rising_ui_rs(acb_t res, const acb_t x, ulong n, ulong m, slong prec)
{
    slong i, k, l, m0, climbs, climbs_max, wp;
    acb_ptr xpow;
    acb_t t, u;
    nn_ptr c;
    TMP_INIT;

    if (n <= 1)
    {
        if (n == 0)
            acb_one(res);
        else
            acb_set_round(res, x, prec);
        return;
    }

    TMP_START;

    if (m == 0 || m == -1)
    {
        if (n <= 6)
            m = 2;
        else if (n <= 16)
            m = 4;
        else if (n <= 40)
            m = 6;
        else
        {
            m0 = n_sqrt(n);
            m = 8 + 0.27 * pow(1.5 * FLINT_MAX(0, prec - 1024), 0.4);
            m = FLINT_MIN(m, m0);
            m = FLINT_MIN(m, 64);
        }
    }

    wp = ARF_PREC_ADD(prec, FLINT_BIT_COUNT(n));

    climbs_max = FLINT_BIT_COUNT(n - 1) * m;
    c = TMP_ALLOC(sizeof(ulong) * climbs_max * m);

    xpow = _acb_vec_init(m + 1);
    _acb_vec_set_powers(xpow, x, m + 1, wp);
    acb_init(t);
    acb_init(u);

    for (k = 0; k < n; k += m)
    {
        l = FLINT_MIN(m, n - k);
        climbs = FLINT_BIT_COUNT(k + l - 1) * l;
        climbs = (climbs + FLINT_BITS - 1) / FLINT_BITS;

        /* assumes l >= 2 */
        if (l == 1)
        {
            acb_add_ui(u, x, k, wp);
        }
        else
        {
            if (climbs == 1)
            {
                _arb_hypgeom_rising_coeffs_1(c, k, l);
                acb_dot_ui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
            }
            else if (climbs == 2)
            {
                _arb_hypgeom_rising_coeffs_2(c, k, l);
                acb_dot_uiui(u, xpow + l, 0, xpow, 1, c, 1, l, wp);
            }
            else
            {
                fmpz * f = (fmpz *) c;

                for (i = 0; i < l; i++)
                    fmpz_init(f + i);

                _arb_hypgeom_rising_coeffs_fmpz(f, k, l);

                acb_dot_fmpz(u, xpow + l, 0, xpow, 1, f, 1, l, wp);

                for (i = 0; i < l; i++)
                    fmpz_clear(f + i);
            }
        }

        if (k == 0)
            acb_swap(t, u);
        else
            acb_mul(t, t, u, wp);
    }

    acb_set_round(res, t, prec);

    acb_clear(t);
    acb_clear(u);
    _acb_vec_clear(xpow, m + 1);
    TMP_END;
}