1 /*
2 This file is part of FLINT.
3
4 FLINT is free software: you can redistribute it and/or modify it under
5 the terms of the GNU Lesser General Public License (LGPL) as published
6 by the Free Software Foundation; either version 2.1 of the License, or
7 (at your option) any later version. See <https://www.gnu.org/licenses/>.
8 */
9
10 /******************************************************************************
11
12 Authored 2016 by Daniel S. Roche; US Government work in the public domain.
13
14 ******************************************************************************/
15
16 #include "flint.h"
17 #include "fmpz_mod_poly.h"
18 #include "fmpz.h"
19 #include "fmpz_vec.h"
20
_fmpz_mod_poly_minpoly_bm(fmpz * poly,const fmpz * seq,slong len,const fmpz_t p)21 slong _fmpz_mod_poly_minpoly_bm(fmpz* poly,
22 const fmpz* seq, slong len, const fmpz_t p)
23 {
24 fmpz *buf, *curpoly, *prevpoly;
25 slong curlen, prevlen;
26 fmpz_t disc;
27 slong i, m;
28
29 buf = _fmpz_vec_init(len + 1);
30 curpoly = poly;
31 _fmpz_vec_zero(curpoly, len + 1);
32 prevpoly = buf;
33 fmpz_init(disc);
34
35 fmpz_one(curpoly + 0);
36 curlen = 1;
37 fmpz_one(prevpoly + 0);
38 prevlen = 1;
39 m = -1; /* m = last switching point */
40
41 for (i = 0; i < len; ++i)
42 {
43 /* compute next discrepancy */
44 _fmpz_vec_dot(disc, curpoly, seq + (i - curlen + 1), curlen);
45 fmpz_mod(disc, disc, p);
46
47 if (fmpz_is_zero(disc));
48 else if (i - m <= curlen - prevlen)
49 {
50 /* quick update; no switch, curlen doesn't change. */
51 slong pos = (curlen - prevlen) - (i - m);
52 _fmpz_vec_scalar_addmul_fmpz(curpoly + pos,
53 prevpoly, prevlen, disc);
54 }
55 else
56 {
57 /* switching update */
58 slong pos = (i - m) - (curlen - prevlen);
59 _fmpz_vec_scalar_mul_fmpz(prevpoly, prevpoly, prevlen, disc);
60 _fmpz_poly_add(prevpoly + pos,
61 prevpoly + pos, FLINT_MAX(0, prevlen - pos), curpoly, curlen);
62 prevlen = curlen + pos;
63
64 fmpz_sub(disc, p, disc);
65 fmpz_invmod(disc, disc, p);
66 _fmpz_mod_poly_scalar_mul_fmpz(curpoly, curpoly, curlen, disc, p);
67
68 FMPZ_VEC_SWAP(curpoly, curlen, prevpoly, prevlen);
69 m = i;
70 }
71 }
72
73 /* make curpoly monic and copy to poly if necessary */
74 fmpz_invmod(disc, curpoly + (curlen - 1), p);
75 _fmpz_mod_poly_scalar_mul_fmpz(poly, curpoly, curlen, disc, p);
76
77 _fmpz_vec_clear(buf, len + 1);
78 fmpz_clear(disc);
79
80 return curlen;
81 }
82