1 /*
2 Copyright (C) 2011 William Hart
3 Copyright (C) 2011 Sebastian Pancratz
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include <gmp.h>
14 #include "flint.h"
15 #include "fmpz.h"
16 #include "fmpz_vec.h"
17 #include "fmpz_poly.h"
18 #include "mpn_extras.h"
19
_fmpz_poly_xgcd_modular(fmpz_t r,fmpz * s,fmpz * t,const fmpz * poly1,slong len1,const fmpz * poly2,slong len2)20 void _fmpz_poly_xgcd_modular(fmpz_t r, fmpz * s, fmpz * t,
21 const fmpz * poly1, slong len1,
22 const fmpz * poly2, slong len2)
23 {
24 mp_ptr G, S, T, A, B, T1, T2;
25 fmpz_t prod;
26 int stabilised = 0, first;
27 mp_limb_t p;
28 flint_bitcnt_t s_bits = 0, t_bits = 0;
29
30 /* Compute resultant of input polys */
31 _fmpz_poly_resultant(r, poly1, len1, poly2, len2);
32
33 if (fmpz_is_zero(r))
34 return;
35
36 fmpz_init(prod);
37 fmpz_one(prod);
38
39 _fmpz_vec_zero(s, len2);
40 _fmpz_vec_zero(t, len1);
41
42 p = (UWORD(1) << (FLINT_BITS - 1));
43
44 G = _nmod_vec_init(4 * len1 + 5 * len2 - 2);
45 S = G + len2;
46 T = S + len2;
47 A = T + len1;
48 B = A + len1;
49 T1 = B + len2;
50 T2 = T1 + (len1 + len2 - 1);
51
52 _nmod_vec_zero(S, len2 + len1); /* S = T = 0 */
53
54 first = 1;
55
56 for (;;)
57 {
58 mp_limb_t R;
59 nmod_t mod;
60
61 /* Get next prime */
62 p = n_nextprime(p, 0);
63
64 /* Resultant mod p */
65 R = fmpz_fdiv_ui(r, p);
66
67 /* If p divides resultant or either leading coeff, discard p */
68 if ((fmpz_fdiv_ui(poly1 + len1 - 1, p) == WORD(0)) ||
69 (fmpz_fdiv_ui(poly2 + len2 - 1, p) == WORD(0)) || (R == 0))
70 continue;
71
72 nmod_init(&mod, p);
73
74 /* Reduce polynomials modulo p */
75 _fmpz_vec_get_nmod_vec(A, poly1, len1, mod);
76 _fmpz_vec_get_nmod_vec(B, poly2, len2, mod);
77
78 if (stabilised) /* CRT has stabilised, probably don't need more xgcds */
79 {
80 slong tlen;
81
82 /* Multiply out A*S + B*T to see if it is R mod p */
83 _fmpz_vec_get_nmod_vec(S, s, len2, mod);
84 _fmpz_vec_get_nmod_vec(T, t, len1, mod);
85
86 _nmod_poly_mul(T1, A, len1, S, len2, mod);
87 _nmod_poly_mul(T2, T, len1, B, len2, mod);
88 _nmod_vec_add(T1, T1, T2, len1 + len2 - 1, mod);
89 tlen = len1 + len2 - 1;
90 FMPZ_VEC_NORM(T1, tlen);
91
92 if (tlen == 1 && T1[0] == R) /* It is, so this prime is good */
93 fmpz_mul_ui(prod, prod, p);
94 else
95 stabilised = 0; /* It's not, keep going with xgcds */
96 }
97
98 if (!stabilised) /* Need to keep computing xgcds mod p */
99 {
100 mp_limb_t RGinv;
101
102 /* Compute xgcd mod p */
103 _nmod_poly_xgcd(G, S, T, A, len1, B, len2, mod);
104 RGinv = n_invmod(G[0], mod.n);
105 RGinv = n_mulmod2_preinv(RGinv, R, mod.n, mod.ninv);
106
107 /* Scale appropriately */
108 _nmod_vec_scalar_mul_nmod(S, S, len2, RGinv, mod);
109 _nmod_vec_scalar_mul_nmod(T, T, len1, RGinv, mod);
110
111 if (first) /* First time around set s and t to S and T */
112 {
113 _fmpz_vec_set_nmod_vec(s, S, len2, mod);
114 _fmpz_vec_set_nmod_vec(t, T, len1, mod);
115 fmpz_set_ui(prod, p);
116
117 stabilised = 1; /* Optimise the case where one prime is enough */
118 first = 0;
119 }
120 else /* Otherwise do CRT */
121 {
122 flint_bitcnt_t new_s_bits, new_t_bits;
123
124 _fmpz_poly_CRT_ui(s, s, len2, prod, S, len2, mod.n, mod.ninv, 1);
125 _fmpz_poly_CRT_ui(t, t, len1, prod, T, len1, mod.n, mod.ninv, 1);
126 fmpz_mul_ui(prod, prod, p);
127
128 /* Check to see if CRT has stabilised */
129 new_s_bits = FLINT_ABS(_fmpz_vec_max_bits(s, len2));
130 new_t_bits = FLINT_ABS(_fmpz_vec_max_bits(t, len1));
131
132 stabilised = (s_bits == new_s_bits && t_bits == new_t_bits);
133
134 s_bits = new_s_bits;
135 t_bits = new_t_bits;
136 }
137 }
138
139 if (stabilised)
140 {
141 slong bound1, bound2, bound;
142
143 bound1 = FLINT_BIT_COUNT(len2)
144 + FLINT_ABS(_fmpz_vec_max_bits(poly1, len1))
145 + FLINT_ABS(_fmpz_vec_max_bits(s, len2));
146 bound2 = FLINT_BIT_COUNT(len2)
147 + FLINT_ABS(_fmpz_vec_max_bits(poly2, len2))
148 + FLINT_ABS(_fmpz_vec_max_bits(t, len1));
149
150 bound = 4 + FLINT_MAX(fmpz_bits(r), FLINT_MAX(bound1, bound2));
151
152 if (fmpz_bits(prod) > bound)
153 break;
154 }
155 }
156
157 _nmod_vec_clear(G);
158 fmpz_clear(prod);
159 }
160
161 void
fmpz_poly_xgcd_modular(fmpz_t r,fmpz_poly_t s,fmpz_poly_t t,const fmpz_poly_t poly1,const fmpz_poly_t poly2)162 fmpz_poly_xgcd_modular(fmpz_t r, fmpz_poly_t s, fmpz_poly_t t,
163 const fmpz_poly_t poly1, const fmpz_poly_t poly2)
164 {
165 if (poly1->length < poly2->length)
166 {
167 fmpz_poly_xgcd_modular(r, t, s, poly2, poly1);
168 } else /* len1 >= len2 >= 0 */
169 {
170 const slong len1 = poly1->length;
171 const slong len2 = poly2->length;
172 fmpz *S, *T;
173 fmpz_poly_t temp1, temp2;
174
175 if (len1 == 0 || len2 == 0)
176 {
177 fmpz_zero(r);
178 }
179 else /* len1 >= len2 >= 1 */
180 {
181 if (s == poly1 || s == poly2)
182 {
183 fmpz_poly_init2(temp1, len2);
184 S = temp1->coeffs;
185 }
186 else
187 {
188 fmpz_poly_fit_length(s, len2);
189 S = s->coeffs;
190 }
191
192 if (t == poly1 || t == poly2)
193 {
194 fmpz_poly_init2(temp2, len1);
195 T = temp2->coeffs;
196 }
197 else
198 {
199 fmpz_poly_fit_length(t, len1);
200 T = t->coeffs;
201 }
202
203 _fmpz_poly_xgcd_modular(r, S, T, poly1->coeffs, len1,
204 poly2->coeffs, len2);
205
206 if (s == poly1 || s == poly2)
207 {
208 fmpz_poly_swap(s, temp1);
209 fmpz_poly_clear(temp1);
210 }
211
212 if (t == poly1 || t == poly2)
213 {
214 fmpz_poly_swap(t, temp2);
215 fmpz_poly_clear(temp2);
216 }
217
218 _fmpz_poly_set_length(s, len2);
219 _fmpz_poly_normalise(s);
220
221 _fmpz_poly_set_length(t, len1);
222 _fmpz_poly_normalise(t);
223 }
224 }
225 }
226
227