1 #ifndef MPZ_POLY_H_
2 #define MPZ_POLY_H_
3
4 // IWYU pragma: no_include "double_poly.h"
5 // (only the fwd-decl is needed)
6 #include <stdio.h>
7 #include <stdint.h>
8 #include <gmp.h>
9 #include "macros.h"
10
11 #define xxxMPZ_POLY_TIMINGS
12 // for timings of roots mod p (beware, this is not thread-safe)
13
14 #ifndef DOUBLE_POLY_H_
15 typedef struct double_poly_s * double_poly_ptr;
16 typedef const struct double_poly_s * double_poly_srcptr;
17 #endif
18
19 #ifdef __cplusplus
20 #include <string>
21 #include <istream> // std::istream // IWYU pragma: keep
22 #include <ostream> // std::ostream // IWYU pragma: keep
23 extern "C" {
24 #endif
25
26 /* maximum degree we can reconstruct using mpz_poly_mul_tc_interpolate */
27 #define MAX_TC_DEGREE 19
28
29
30 /* Note, deg = -1 means P=0; otherwise, one should have coeff[deg] != 0.
31 Warning: a polynomial of degree d needs d+1 allocation. */
32
33 struct mpz_poly_s {
34 int alloc;
35 int deg;
36 mpz_t *coeff;
37 };
38 #ifndef DOUBLE_POLY_H_
39 /* double_poly.h forward-declares these. Don't do it twice */
40 typedef struct mpz_poly_s * mpz_poly_ptr;
41 typedef const struct mpz_poly_s * mpz_poly_srcptr;
42 #endif
43 typedef struct mpz_poly_s mpz_poly[1];
44
45 /* -------------------------------------------------------------------------- */
46
47 /* [LI] This should also be renamed to mpz_polymodF. */
48 /* And may be moved outside mpz_poly.[ch]? */
49
50 /* Let F(x) be a (non-monic) polynomial of degree d:
51 F(x) = f_d x^d + f_{d-1} x^{d-1} + .... + f_1 x + f_0
52 The following type represents a polynomial modulo F(x):
53 If P is such an element, it means: P = P.p / f_d^P.v */
54 typedef struct {
55 mpz_poly p;
56 int v;
57 } polymodF_struct_t;
58
59 typedef polymodF_struct_t polymodF_t[1];
60 /* -------------------------------------------------------------------------- */
61
62 /* Special structure to represent coeffs of a polynomial in base p^k. */
63 typedef struct {
64 int deg;
65 char **coeff;
66 } poly_base_struct_t;
67
68 typedef poly_base_struct_t poly_base_t[1];
69
70 /* Note on parallelism.
71 *
72 * Some functions in this API can optionally use openmp. We want to make
73 * sure that the code that calls these openmp-enabled functions does so
74 * **willingly**. And we don't want to pollute the "normal" interface.
75 *
76 * The chosen solution is as follows.
77 *
78 * - the functions declared here, and use as plain (e.g.) mpz_poly_mul
79 * resolve to something that does **NOT** use openmp.
80 *
81 * - to use openmp, include mpz_poly_parallel.hpp instead, and call the
82 * member functions of an mpz_poly_parallel_info object. E.g.
83 *
84 * mpz_poly_parallel_info inf;
85 * // (maybe add some configuration code for the inf object, if the
86 * // need for that ever appears)
87 * inf.mpz_poly_mul(....)
88 *
89 * More detail on can be found in mpz_poly_parallel.hpp and mpz_poly.cpp
90 */
91
92 /* Management of the structure, set and print coefficients. */
93 void mpz_poly_init(mpz_poly_ptr, int d);
94 void mpz_poly_realloc (mpz_poly_ptr f, int nc);
95 void mpz_poly_set(mpz_poly_ptr g, mpz_poly_srcptr f);
96 void mpz_poly_swap (mpz_poly_ptr f, mpz_poly_ptr g);
97 void mpz_poly_clear(mpz_poly_ptr f);
mpz_poly_degree(mpz_poly_srcptr f)98 static inline int mpz_poly_degree(mpz_poly_srcptr f) { return f->deg; }
99
100 void mpz_poly_cleandeg(mpz_poly_ptr f, int deg);
101 void mpz_poly_setcoeffs(mpz_poly_ptr f, mpz_t * coeffs, int d);
102 void mpz_poly_set_zero(mpz_poly_ptr f);
103 void mpz_poly_set_xi(mpz_poly_ptr f, int i);
104 void mpz_poly_set_mpz(mpz_poly_ptr f, mpz_srcptr z);
105 void mpz_poly_set_double_poly(mpz_poly_ptr g, double_poly_srcptr f);
106
107 void mpz_poly_init_set_ab (mpz_poly_ptr rel, int64_t a, uint64_t b);
108 void mpz_poly_init_set_mpz_ab (mpz_poly_ptr rel, mpz_srcptr a, mpz_srcptr b);
109
110 void mpz_poly_setcoeff(mpz_poly_ptr f, int i, mpz_srcptr z);
111 void mpz_poly_setcoeff_si(mpz_poly_ptr f, int i, long z);
112 void mpz_poly_setcoeff_ui(mpz_poly_ptr f, int i, unsigned long z);
113 void mpz_poly_setcoeff_int64(mpz_poly_ptr f, int i, int64_t z);
114 void mpz_poly_setcoeff_uint64(mpz_poly_ptr f, int i, uint64_t z);
115 void mpz_poly_setcoeff_double(mpz_poly_ptr f, int i, double z);
116 void mpz_poly_getcoeff(mpz_ptr res, int i, mpz_poly_srcptr f);
117
118 /* functions for Joux-Lercier and generalized Joux-Lercier */
119 int mpz_poly_setcoeffs_counter(mpz_poly_ptr f, int* max_abs_coeffs, unsigned long *next_counter, int deg, unsigned long counter, unsigned int bound);
120 void mpz_poly_setcoeffs_counter_print_error_code(int error_code);
121 unsigned long mpz_poly_getcounter(mpz_poly_ptr f, unsigned int bound);
122 unsigned long mpz_poly_cardinality(int deg, unsigned int bound);
123
124 /* return the leading coefficient of f */
mpz_poly_lc(mpz_poly_srcptr f)125 static inline mpz_srcptr mpz_poly_lc (mpz_poly_srcptr f) {
126 ASSERT(f->deg >= 0);
127 return f->coeff[f->deg];
128 }
129
130 /* Print functions */
131 int mpz_poly_asprintf(char ** res, mpz_poly_srcptr f);
132 /* Print coefficients of f.
133 * endl = 1 if "\n" at the end of fprintf. */
134 void mpz_poly_fprintf_endl (FILE *fp, mpz_poly_srcptr f, int endl);
135 void mpz_poly_fprintf(FILE *fp, mpz_poly_srcptr f);
136 void mpz_poly_fprintf_coeffs (FILE *fp, mpz_poly_srcptr f, const char sep);
137 void mpz_poly_fscanf_coeffs (FILE *fp, mpz_poly_ptr f, const char sep);
138 void mpz_poly_fprintf_cado_format (FILE *fp, mpz_poly_srcptr f,
139 const char letter, const char *pre);
140 void mpz_poly_print_raw(mpz_poly_srcptr f);
141 #ifdef MPZ_POLY_TIMINGS
142 void print_timings_pow_mod_f_mod_p();
143 #endif
144 /* Tests and comparison functions */
145 int mpz_poly_cmp (mpz_poly_srcptr, mpz_poly_srcptr);
146 int mpz_poly_normalized_p (mpz_poly_srcptr f);
147 int mpz_poly_is_monic (mpz_poly_srcptr f);
148
149 /* Polynomial arithmetic */
150 void mpz_poly_to_monic(mpz_poly_ptr g, mpz_poly_srcptr f);
151 void mpz_poly_neg(mpz_poly_ptr f, mpz_poly_srcptr g);
152 void mpz_poly_add(mpz_poly_ptr f, mpz_poly_srcptr g, mpz_poly_srcptr h);
153 void mpz_poly_sub(mpz_poly_ptr f, mpz_poly_srcptr g, mpz_poly_srcptr h);
154 void mpz_poly_add_ui(mpz_poly_ptr g, mpz_poly_srcptr f, unsigned long a);
155 void mpz_poly_sub_ui(mpz_poly_ptr g, mpz_poly_srcptr f, unsigned long a);
156 void mpz_poly_add_mpz(mpz_poly_ptr f, mpz_poly_srcptr g, mpz_srcptr a);
157 void mpz_poly_sub_mpz(mpz_poly_ptr f, mpz_poly_srcptr g, mpz_srcptr a);
158 void mpz_poly_sub_mod_mpz(mpz_poly_ptr f, mpz_poly_srcptr g, mpz_poly_srcptr h, mpz_srcptr m);
159 void mpz_poly_mul(mpz_poly_ptr f, mpz_poly_srcptr g, mpz_poly_srcptr h);
160 void mpz_poly_mul_mpz(mpz_poly_ptr Q, mpz_poly_srcptr P, mpz_srcptr a);
161 void mpz_poly_divexact_mpz(mpz_poly_ptr Q, mpz_poly_srcptr P, mpz_srcptr a);
162 int mpz_poly_divisible_mpz (mpz_poly_srcptr P, mpz_srcptr a);
163 void mpz_poly_translation (mpz_poly_ptr, mpz_poly_srcptr, mpz_srcptr);
164 void mpz_poly_rotation (mpz_poly_ptr, mpz_poly_srcptr, mpz_poly_srcptr, mpz_srcptr, int);
165 void mpz_poly_addmul_si (mpz_poly_ptr, mpz_poly_srcptr, long);
166 void mpz_poly_mul_si (mpz_poly_ptr, mpz_poly_srcptr, long);
167 void mpz_poly_divexact_ui (mpz_poly_ptr, mpz_poly_srcptr, unsigned long);
168 void mpz_poly_rotation_int64 (mpz_poly_ptr, mpz_poly_srcptr, mpz_poly_srcptr, const int64_t, int);
169 void mpz_poly_makemonic_mod_mpz (mpz_poly_ptr Q, mpz_poly_srcptr P, mpz_srcptr m);
170 void barrett_precompute_inverse (mpz_ptr invm, mpz_srcptr m);
171 int mpz_poly_mod_f_mod_mpz(mpz_poly_ptr R, mpz_poly_srcptr f, mpz_srcptr m, mpz_srcptr invf, mpz_srcptr invm);
172 int mpz_poly_mod_mpz(mpz_poly_ptr R, mpz_poly_srcptr A, mpz_srcptr m, mpz_srcptr invm);
173 int mpz_poly_mod_mpz_lazy (mpz_poly_ptr R, mpz_poly_srcptr A, mpz_srcptr m);
174 void mpz_poly_mul_mod_f_mod_mpz(mpz_poly_ptr Q, mpz_poly_srcptr P1, mpz_poly_srcptr P2, mpz_poly_srcptr f, mpz_srcptr m, mpz_srcptr invf, mpz_srcptr invm);
175 void mpz_poly_mul_mod_f (mpz_poly_ptr Q, mpz_poly_srcptr P1, mpz_poly_srcptr P2, mpz_poly_srcptr f);
176 void mpz_poly_reduce_frac_mod_f_mod_mpz (mpz_poly_ptr num, mpz_poly_ptr denom, mpz_poly_srcptr F, mpz_srcptr m);
177 int mpz_poly_div_qr (mpz_poly_ptr q, mpz_poly_ptr r, mpz_poly_srcptr f, mpz_poly_srcptr g, mpz_srcptr p);
178 int mpz_poly_div_r (mpz_poly_ptr h, mpz_poly_srcptr f, mpz_srcptr p);
179 int mpz_poly_div_qr_z (mpz_poly_ptr q, mpz_poly_ptr r, mpz_poly_srcptr f, mpz_poly_srcptr g);
180 int mpz_poly_div_r_z (mpz_poly_ptr r, mpz_poly_srcptr f, mpz_poly_srcptr g);
181 int mpz_poly_divexact (mpz_poly_ptr q, mpz_poly_srcptr h, mpz_poly_srcptr f, mpz_srcptr p);
182 void mpz_poly_div_2_mod_mpz(mpz_poly_ptr f, mpz_poly_srcptr g, mpz_srcptr m);
183 void mpz_poly_div_xi(mpz_poly_ptr g, mpz_poly_srcptr f, int i);
184 void mpz_poly_mul_xi(mpz_poly_ptr g, mpz_poly_srcptr f, int i);
185 void mpz_poly_mul_xplusa(mpz_poly_ptr g, mpz_poly_srcptr f, mpz_srcptr a);
186
187
188 void mpz_poly_eval(mpz_ptr res, mpz_poly_srcptr f, mpz_srcptr x);
189 void mpz_poly_eval_ui (mpz_ptr res, mpz_poly_srcptr f, unsigned long x);
190 void mpz_poly_eval_diff_ui (mpz_ptr res, mpz_poly_srcptr f, unsigned long x);
191 void mpz_poly_eval_diff (mpz_ptr res, mpz_poly_srcptr f, mpz_srcptr x);
192 void mpz_poly_eval_poly(mpz_poly_ptr res, mpz_poly_srcptr f, mpz_poly_srcptr x);
193 void mpz_poly_eval_diff_poly (mpz_poly_ptr res, mpz_poly_srcptr f, mpz_poly_srcptr x);
194 void mpz_poly_eval_mod_mpz(mpz_t res, mpz_poly_srcptr f, mpz_srcptr x, mpz_srcptr m);
195 int mpz_poly_is_root(mpz_poly_srcptr poly, mpz_srcptr root, mpz_srcptr modulus);
196 void mpz_poly_eval_several_mod_mpz(mpz_ptr * res, mpz_poly_srcptr * f, int k, mpz_srcptr x, mpz_srcptr m);
197 void polymodF_mul(polymodF_t Q, const polymodF_t P1, const polymodF_t P2, mpz_poly_srcptr F);
198 void mpz_poly_sqr_mod_f_mod_mpz(mpz_poly_ptr Q, mpz_poly_srcptr P, mpz_poly_srcptr f, mpz_srcptr m, mpz_srcptr invf, mpz_srcptr invm);
199 void mpz_poly_pow_ui(mpz_poly_ptr B, mpz_poly_srcptr A, unsigned long n);
200 void mpz_poly_pow_ui_mod_f(mpz_poly_ptr B, mpz_poly_srcptr A, unsigned long n, mpz_poly_srcptr f);
201 void mpz_poly_pow_mod_f_mod_ui(mpz_poly_ptr Q, mpz_poly_srcptr P, mpz_poly_srcptr f, mpz_srcptr a, unsigned long p);
202 void mpz_poly_pow_mod_f_mod_mpz(mpz_poly_ptr Q, mpz_poly_srcptr P, mpz_poly_srcptr f, mpz_srcptr a, mpz_srcptr p);
203 void mpz_poly_pow_ui_mod_f_mod_mpz (mpz_poly_ptr Q, mpz_poly_srcptr P, mpz_poly_srcptr f, unsigned long a, mpz_srcptr p);
204 void mpz_poly_derivative(mpz_poly_ptr df, mpz_poly_srcptr f);
205 mpz_poly* mpz_poly_base_modp_init (mpz_poly_srcptr P0, int p, unsigned long *K, int l);
206 void mpz_poly_base_modp_clear (mpz_poly *P, int l);
207 void mpz_poly_base_modp_lift(mpz_poly_ptr a, mpz_poly *P, int k, mpz_srcptr pk);
208 size_t mpz_poly_sizeinbase (mpz_poly_srcptr f, int base);
209 size_t mpz_poly_size (mpz_poly_srcptr f);
210 void mpz_poly_infinity_norm(mpz_ptr in, mpz_poly_srcptr f);
211 size_t mpz_poly_totalsize (mpz_poly_srcptr f);
212 void mpz_poly_gcd_mpz (mpz_poly_ptr h, mpz_poly_srcptr f, mpz_poly_srcptr g, mpz_srcptr p);
213 // compute f = GCD(f,g) mod N. If this fails, put the factor in the last
214 // given argument.
215 int mpz_poly_pseudogcd_mpz(mpz_poly_ptr , mpz_poly_ptr , mpz_srcptr , mpz_ptr);
216 void mpz_poly_xgcd_mpz(mpz_poly_ptr gcd, mpz_poly_srcptr f, mpz_poly_srcptr g, mpz_poly_ptr u, mpz_poly_ptr v, mpz_srcptr p);
217 void mpz_poly_homography (mpz_poly_ptr Fij, mpz_poly_srcptr F, int64_t H[4]);
218 void mpz_poly_homogeneous_eval_siui (mpz_ptr v, mpz_poly_srcptr f, const int64_t i, const uint64_t j);
219 void mpz_poly_content (mpz_ptr c, mpz_poly_srcptr F);
220 void mpz_poly_resultant(mpz_ptr res, mpz_poly_srcptr p, mpz_poly_srcptr q);
221 void mpz_poly_discriminant(mpz_ptr res, mpz_poly_srcptr f);
222 int mpz_poly_squarefree_p(mpz_poly_srcptr f);
223 int mpz_poly_is_irreducible_z(mpz_poly_srcptr f);
224
225 int mpz_poly_number_of_real_roots(mpz_poly_srcptr f);
226
227 struct mpz_poly_with_m_s {
228 mpz_poly f;
229 int m;
230 };
231 typedef struct mpz_poly_with_m_s mpz_poly_with_m[1];
232 typedef struct mpz_poly_with_m_s * mpz_poly_with_m_ptr;
233 typedef const struct mpz_poly_with_m_s * mpz_poly_with_m_srcptr;
234
235 struct mpz_poly_factor_list_s {
236 mpz_poly_with_m * factors;
237 int alloc;
238 int size;
239 };
240 typedef struct mpz_poly_factor_list_s mpz_poly_factor_list[1];
241 typedef struct mpz_poly_factor_list_s * mpz_poly_factor_list_ptr;
242 typedef const struct mpz_poly_factor_list_s * mpz_poly_factor_list_srcptr;
243
244 void mpz_poly_factor_list_init(mpz_poly_factor_list_ptr l);
245 void mpz_poly_factor_list_clear(mpz_poly_factor_list_ptr l);
246 void mpz_poly_factor_list_flush(mpz_poly_factor_list_ptr l);
247 void mpz_poly_factor_list_push(mpz_poly_factor_list_ptr l, mpz_poly_srcptr f, int m);
248 void mpz_poly_factor_list_fprintf(FILE* fp, mpz_poly_factor_list_srcptr l);
249 int mpz_poly_factor_sqf(mpz_poly_factor_list_ptr lf, mpz_poly_srcptr f, mpz_srcptr p);
250 int mpz_poly_factor_ddf(mpz_poly_factor_list_ptr lf, mpz_poly_srcptr f0, mpz_srcptr p);
251 int mpz_poly_factor_edf(mpz_poly_factor_list_ptr lf, mpz_poly_srcptr f, int k, mpz_srcptr p, gmp_randstate_t rstate);
252
253 /* output is sorted by degree and lexicographically */
254 int mpz_poly_factor(mpz_poly_factor_list lf, mpz_poly_srcptr f, mpz_srcptr p, gmp_randstate_t rstate);
255 int mpz_poly_is_irreducible(mpz_poly_srcptr f, mpz_srcptr p);
256
257 /* lift from a factor list mod ell to a factor list mod ell2.
258 * ell does not need to be prime, provided all factors considered are
259 * unitary.
260 *
261 * ell and ell2 must be powers of the same prime, with ell2 <= ell^2
262 * (NOTE that this is not checked)
263 */
264 int mpz_poly_factor_list_lift(mpz_poly_factor_list_ptr fac, mpz_poly_srcptr f, mpz_srcptr ell, mpz_srcptr ell2);
265
266 /* This computes the ell-adic lifts of the factors of f, assuming
267 * we have no multiplicities, using Newton lifting.
268 * This requires that f be monic
269 *
270 * I'm terribly lazy, so at the moment this is working only for prec==2.
271 * Extending to arbitrary p is an easy exercise.
272 *
273 * The output is sorted based on the order of the factors mod p (that is,
274 * factors are the lifts of the factors returned by mpz_poly_factor mod
275 * p, in the same order).
276 */
277 int mpz_poly_factor_and_lift_padically(mpz_poly_factor_list_ptr fac, mpz_poly_srcptr f, mpz_srcptr ell, int prec, gmp_randstate_t rstate);
278
279 #ifdef __cplusplus
280 }
281 #endif
282
283 #ifdef __cplusplus
284 /* This is sort of a generic way to write a c++ equivalent to the C type.
285 * The first-class citizen in the cado-nfs code is (still) the C type, so
286 * we're definitely bound to have a few infelicities here:
287 * - the type name can't be the same because of the size-1 array trick
288 * in C.
289 * - the C type is embedded as a member x for the same reason.
290 * - most operations on the C type should go through the member x
291 * (however, the conversions we have to _ptr and _srcptr can ease
292 * things a bit).
293 */
294 struct cxx_mpz_poly {
295 mpz_poly x;
cxx_mpz_polycxx_mpz_poly296 cxx_mpz_poly() { mpz_poly_init(x, -1); }
cxx_mpz_polycxx_mpz_poly297 cxx_mpz_poly(int deg) { mpz_poly_init(x, deg); }
degreecxx_mpz_poly298 inline int degree() const { return x->deg; } /* handy */
cxx_mpz_polycxx_mpz_poly299 cxx_mpz_poly(mpz_poly_srcptr f) { mpz_poly_init(x, -1); mpz_poly_set(x, f); }
~cxx_mpz_polycxx_mpz_poly300 ~cxx_mpz_poly() { mpz_poly_clear(x); }
cxx_mpz_polycxx_mpz_poly301 cxx_mpz_poly(cxx_mpz_poly const & o) {
302 mpz_poly_init(x, -1);
303 mpz_poly_set(x, o.x);
304 }
305 cxx_mpz_poly & operator=(cxx_mpz_poly const & o) {
306 mpz_poly_set(x, o.x);
307 return *this;
308 }
309 #if __cplusplus >= 201103L
cxx_mpz_polycxx_mpz_poly310 cxx_mpz_poly(cxx_mpz_poly && o) {
311 mpz_poly_init(x, -1);
312 mpz_poly_swap(x, o.x);
313 }
314 cxx_mpz_poly& operator=(cxx_mpz_poly && o) {
315 mpz_poly_swap(x, o.x);
316 return *this;
317 }
318 #endif
mpz_poly_ptrcxx_mpz_poly319 operator mpz_poly_ptr() { return x; }
mpz_poly_srcptrcxx_mpz_poly320 operator mpz_poly_srcptr() const { return x; }
321 mpz_poly_ptr operator->() { return x; }
322 mpz_poly_srcptr operator->() const { return x; }
323 std::string print_poly(std::string const& var) const;
324 };
325
326 std::ostream& operator<<(std::ostream& o, cxx_mpz_poly const & f);
327 std::istream& operator>>(std::istream& in, cxx_mpz_poly & f);
328
329 #if GNUC_VERSION_ATLEAST(4,3,0)
330 extern void mpz_poly_init(cxx_mpz_poly & pl, int) __attribute__((error("mpz_poly_init must not be called on a mpz_poly reference -- it is the caller's business (via a ctor)")));
331 extern void mpz_poly_clear(cxx_mpz_poly & pl) __attribute__((error("mpz_poly_clear must not be called on a mpz_poly reference -- it is the caller's business (via a dtor)")));
332 #endif
333
334 #endif
335
336 #endif /* MPZ_POLY_H_ */
337