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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
/*
Copyright (C) 2025, Vincent Neiger, Éric Schost
Copyright (C) 2025, Mael Hostettler
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 "nmod.h"
#include "nmod_poly.h"
#include "nmod_vec.h"
void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly,
nn_srcptr v, const nmod_geometric_progression_t G, slong len, nmod_t mod)
{
FLINT_ASSERT(len <= G->len);
FLINT_ASSERT((G->function & UWORD(2)) == UWORD(2));
if (len == 1)
{
poly[0] = v[0];
return;
}
/** step1: Newton interpolation
* [Bostan - Schost, J.Complexity 2005, Section 5.1]
* -> The coefficients of the interpolant, in the Newton basis associated
* to the geometric progression 1, q, q**2, q**3, etc., are obtained as
* c_0 / q_0, ..., c_{len-1} / q_{len-1}, where c_0, ..., c_{len-1} are
* the first coefficients of the product
* (sum_{i=0}^{len-1} v[i]/u_i x**i) (sum_{i=0}^{len-1} (-1)**i q_i/u_i x**i)
* where v[i] is the element at index `i` in the input values `v`,
* and u_i = prod_{1 <= k <= i} (q**k - 1),
* and q_i = q**(i*(i-1)/2)
* -> With the precomputed data, these are the `len` coefficients of
* f * G->int_f mod x**len
* where f = sum_{i=0}^{len-1} v[i] * G->int_s1[i] x**i
*/
/* val = valuation of output poly in Newton basis */
slong val = 0;
for (; val < len; val++)
if (v[val] != 0)
break;
if (val == len)
{
_nmod_vec_zero(poly, len);
return;
}
slong f_len, h_len;
nn_ptr f = _nmod_vec_init(len - val);
nn_ptr h = _nmod_vec_init(len);
/* actual length of f */
for (f_len = len; f_len > val; f_len--)
if (v[f_len-1] != 0)
break;
f_len = f_len - val;
/* f = sum_{i=val}^{len-1} v[i] * G->int_s1[i] x**{i-val} */
/* == sum_{i=0}^{f_len-1} v[i+val] * G->int_s1[i+val] x**i */
for (slong i = 0; i < f_len; i++)
f[i] = nmod_mul(v[i+val], G->int_s1[i+val], mod);
/* h = (x**val * f) * G->int_f mod x**len */
/* == x**val (f * G->int_f mod x**(len-val)) */
/* note: len - val is <= G->int_f->length, since G->int_f has */
/* length G->len >= len (all its coefficients are nonzero) */
_nmod_poly_mullow(h+val, G->int_f->coeffs, len - val, f, f_len, len - val, mod);
/* for Newton interpolation, here we should compute h[i] = h[i]/q_i */
/* yet this "/q_i" will simplify with another operation just below, */
/* so we just leave h as it is */
/** step2: Newton basis -> monomial basis
* [Bostan - Schost, J.Complexity 2005, Section 5.2]
* -> Convert h[i]/q_i to monomial basis, through the transposed
* x**len-truncated multiplication of two polynomials
* sum_{i=0}^{len-1} uu_i x**i
* and sum_{i=0}^{len-1} (-1)**i * (h[i]/q_i)*q_i/uu_i x**i)
* where q_i = q**(i*(i-1)/2) as above,
* and uu_i = prod_{1 <= k <= i} q**(k-1) / (1 - q**k)
* == (-1)**i * q_i / u_i, for u_i as above
* (in the paper this is prod q**k / (1 - q**k), is this a typo?)
* This gives `len` coefficients that must then be scaled
* by (-1)**i * uu_i / q_i == 1 / u_i
* -> Transposing the truncated product of poly1,poly2 of degree < len,
* mullow_t(res, poly1, poly2, len)
* simply amounts to a (non-tranposed) mullow and reversals:
* mullow(res, poly1, rev(poly2, len), len)
* res = rev(res, len)
* -> With the precomputed data, the two polynomials have coefficients
* uu_i == G->int_f[i] and (-1)**i * h[i]/uu_i == h[i] * G->int_s2[i]
* meaning that we want to compute
* mullow_t(res, G->int_f, F, len)
* where F = sum_{i=0}^{len-1} h[i] * G->int_s2[i] x**i,
* and then scale by 1/u_i == G->int_s1[i]
*/
/* valuation of h is at least val, see if it is higher */
for (; val < len; val++)
if (h[val] != 0)
break;
/* actual length of h */
for (h_len = len; h_len > val; h_len--)
if (h[h_len-1] != 0)
break;
/* compute reversed and scaled f */
for (slong i = 0; i < h_len-val; i++)
f[i] = nmod_mul(h[h_len-1-i], G->int_s2[h_len-1-i], mod);
/* transposed short product */
_nmod_poly_mullow(h+len-h_len, G->int_f->coeffs, h_len, f, h_len-val, h_len, mod);
/* final scaling */
_nmod_vec_zero(poly+h_len, len-h_len);
for (slong i = 0; i < h_len; i++)
poly[i] = nmod_mul(h[len-1-i], G->int_s1[i], mod);
_nmod_vec_clear(f);
_nmod_vec_clear(h);
}
void nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly,
nn_srcptr v, const nmod_geometric_progression_t G, slong len)
{
FLINT_ASSERT((G->function & UWORD(2)) == UWORD(2));
nmod_poly_fit_length(poly, len);
_nmod_poly_set_length(poly, len);
_nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly->coeffs, v, G, len, G->mod);
_nmod_poly_normalise(poly);
}
void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly,
ulong r, nn_srcptr ys, slong len)
{
nmod_geometric_progression_t G;
_nmod_geometric_progression_init_function(G, r, len, poly->mod, UWORD(2));
nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, ys, G, len);
nmod_geometric_progression_clear(G);
}