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
117
118
119
120
121
122
123
124
125
126
127
128
129

/*
    Copyright (C) 2021 Albin Ahlbäck
    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 "fmpz.h"
#include "fmpz_vec.h"
#include "fmpz_poly.h"

/* Up to this order the coefficients fit inside small fmpz */
#define SMALL_BOUND (FLINT_BITS == 64 ? 20 : 12)

static void
_fmpz_vec_powers(fmpz * res, ulong n, slong len)
{
    slong ix;

    if (len != 0)
        fmpz_set_ui(res + 0, n == 0);

    for (ix = 1; ix < len; ix += 2)
    {
        fmpz_set_ui(res + ix, ix);
        fmpz_pow_ui(res + ix, res + ix, n);
    }

    for (ix = 2; ix < len; ix += 2)
        fmpz_mul_2exp(res + ix, res + ix / 2, n);
}

/* Only valid for when 1 <= len <= n / 2 */
static void
_fmpz_vec_binomials(fmpz * res, ulong n, slong len)
{
    slong ix, jx;

    fmpz_one(res + 0);
    fmpz_set_ui(res + 1, n);

    for (ix = 2, jx = 1; ix < len; ix++, jx++)
    {
        fmpz_mul_ui(res + ix, res + jx, n - jx);
        fmpz_divexact_si(res + ix, res + ix, ix);
    }
}

static void
_fmpz_poly_eulerian_polynomial_series(fmpz * res, ulong n)
{
    slong ix, m;
    fmpz * tmp;

    m = (n + 1) / 2;
    tmp = _fmpz_vec_init(2 * m + 1);
    _fmpz_vec_binomials(tmp, n + 1, m);
    for (ix = 1; ix < m; ix += 2)
        fmpz_neg(tmp + ix, tmp + ix);
    _fmpz_vec_powers(tmp + m, n, m + 1);
    _fmpz_poly_mullow(res, tmp, m, tmp + m + 1, m, m);
    _fmpz_vec_clear(tmp, 2 * m + 1);
}

static void
_fmpz_poly_eulerian_polynomial_rec(fmpz * res, ulong n)
{
    slong jx;
    ulong ix;

    fmpz_one(res);
    for (ix = 1; ix <= FLINT_MIN(n / 2, SMALL_BOUND / 2); ix++)
        _fmpz_demote(res + ix);

    for (ix = 3; ix <= n; ix++)
    {
        if (ix <= SMALL_BOUND)
        {
            if (ix & 1)
                res[ix / 2] = (ix + 1) * res[ix / 2 - 1];
            for (jx = ix / 2 - 1; jx >= 1; jx--)
                res[jx] = (ix - jx) * res[jx - 1] + (jx + 1) * res[jx];
        }
        else
        {
            if (ix & 1)
                fmpz_mul_ui(res + ix / 2, res + ix / 2 - 1, ix + 1);
            for (jx = ix / 2 - 1; jx >= 1; jx--)
            {
                fmpz_mul_ui(res + jx, res + jx, jx + 1);
                fmpz_addmul_ui(res + jx, res + (jx - 1), ix - jx);
            }
        }
    }
}

static void
_fmpz_poly_eulerian_polynomial(fmpz * res, ulong n)
{
    ulong ix;

    if (n < 32)
        _fmpz_poly_eulerian_polynomial_rec(res, n);
    else
        _fmpz_poly_eulerian_polynomial_series(res, n);

    for (ix = 0; ix < n / 2; ix++)
        fmpz_set(res + (n - ix - 1), res + ix);
}

void
fmpz_poly_eulerian_polynomial(fmpz_poly_t poly, ulong n)
{
    /* A_0, A_1 = 1 */
    if (n <= 1)
    {
        fmpz_poly_one(poly);
        return;
    }

    fmpz_poly_fit_length(poly, n);
    _fmpz_poly_eulerian_polynomial(poly->coeffs, n);
    _fmpz_poly_set_length(poly, n);
}