1 /*
2 Copyright (C) 2011 William Hart
3 Copyright (C) 2012 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 <https://www.gnu.org/licenses/>.
11 */
12
13 #include <stdlib.h>
14 #include "fmpz_vec.h"
15 #include "fmpz_mod_poly.h"
16
_fmpz_mod_poly_xgcd_euclidean(fmpz * G,fmpz * S,fmpz * T,const fmpz * A,slong lenA,const fmpz * B,slong lenB,const fmpz_t invB,const fmpz_t p)17 slong _fmpz_mod_poly_xgcd_euclidean(fmpz *G, fmpz *S, fmpz *T,
18 const fmpz *A, slong lenA,
19 const fmpz *B, slong lenB,
20 const fmpz_t invB, const fmpz_t p)
21 {
22 _fmpz_vec_zero(G, lenB);
23 _fmpz_vec_zero(S, lenB - 1);
24 _fmpz_vec_zero(T, lenA - 1);
25
26 if (lenB == 1)
27 {
28 fmpz_set(G + 0, B + 0);
29 fmpz_one(T + 0);
30 return 1;
31 }
32 else
33 {
34 fmpz *Q, *R;
35 slong lenQ, lenR;
36 TMP_INIT;
37 TMP_START;
38
39 FMPZ_VEC_TMP_INIT(Q, 2*lenA);
40 R = Q + lenA;
41
42 _fmpz_mod_poly_divrem(Q, R, A, lenA, B, lenB, invB, p);
43 lenR = lenB - 1;
44 FMPZ_VEC_NORM(R, lenR);
45
46 if (lenR == 0)
47 {
48 _fmpz_vec_set(G, B, lenB);
49 fmpz_one(T + 0);
50
51 FMPZ_VEC_TMP_CLEAR(Q, 2*lenA);
52
53 TMP_END;
54 return lenB;
55 }
56 else
57 {
58 fmpz_t inv;
59 fmpz *D, *U, *V1, *V3, *W;
60 slong lenD, lenU, lenV1, lenV3, lenW;
61
62 fmpz_init(inv);
63 FMPZ_VEC_TMP_INIT(W, FLINT_MAX(5*lenB, lenA + lenB));
64 D = W + lenB;
65 U = D + lenB;
66 V1 = U + lenB;
67 V3 = V1 + lenB;
68
69 lenU = 0;
70 _fmpz_vec_set(D, B, lenB);
71 lenD = lenB;
72 fmpz_one(V1 + 0);
73 lenV1 = 1;
74 lenV3 = 0;
75 FMPZ_VEC_SWAP(V3, lenV3, R, lenR);
76
77 do {
78 fmpz_invmod(inv, V3 + (lenV3 - 1), p);
79 _fmpz_mod_poly_divrem_basecase(Q, D, D, lenD, V3, lenV3, inv, p);
80 lenQ = lenD - lenV3 + 1;
81 lenD = lenV3 - 1;
82 FMPZ_VEC_NORM(D, lenD);
83
84 if (lenV1 >= lenQ)
85 _fmpz_mod_poly_mul(W, V1, lenV1, Q, lenQ, p);
86 else
87 _fmpz_mod_poly_mul(W, Q, lenQ, V1, lenV1, p);
88 lenW = lenQ + lenV1 - 1;
89
90 _fmpz_mod_poly_sub(U, U, lenU, W, lenW, p);
91 lenU = FLINT_MAX(lenU, lenW);
92 FMPZ_VEC_NORM(U, lenU);
93
94 FMPZ_VEC_SWAP(U, lenU, V1, lenV1);
95 FMPZ_VEC_SWAP(D, lenD, V3, lenV3);
96
97 } while (lenV3 != 0);
98
99 _fmpz_vec_set(G, D, lenD);
100 _fmpz_vec_set(S, U, lenU);
101 {
102 lenQ = lenA + lenU - 1;
103
104 _fmpz_mod_poly_mul(Q, A, lenA, S, lenU, p);
105 _fmpz_mod_poly_neg(Q, Q, lenQ, p);
106 _fmpz_mod_poly_add(Q, G, lenD, Q, lenQ, p);
107
108 _fmpz_mod_poly_divrem(T, W, Q, lenQ, B, lenB, invB, p);
109 }
110
111 FMPZ_VEC_TMP_CLEAR(W, FLINT_MAX(5*lenB, lenA + lenB));
112 FMPZ_VEC_TMP_CLEAR(Q, 2*lenA);
113 fmpz_clear(inv);
114
115 TMP_END;
116 return lenD;
117 }
118 }
119 }
120
fmpz_mod_poly_xgcd_euclidean(fmpz_mod_poly_t G,fmpz_mod_poly_t S,fmpz_mod_poly_t T,const fmpz_mod_poly_t A,const fmpz_mod_poly_t B,const fmpz_mod_ctx_t ctx)121 void fmpz_mod_poly_xgcd_euclidean(fmpz_mod_poly_t G, fmpz_mod_poly_t S,
122 fmpz_mod_poly_t T, const fmpz_mod_poly_t A,
123 const fmpz_mod_poly_t B, const fmpz_mod_ctx_t ctx)
124 {
125 if (A->length < B->length)
126 {
127 fmpz_mod_poly_xgcd_euclidean(G, T, S, B, A, ctx);
128 }
129 else /* lenA >= lenB >= 0 */
130 {
131 const slong lenA = A->length, lenB = B->length;
132 const fmpz * p = fmpz_mod_ctx_modulus(ctx);
133 fmpz_t inv;
134
135 fmpz_init(inv);
136 if (lenA == 0) /* lenA = lenB = 0 */
137 {
138 fmpz_mod_poly_zero(G, ctx);
139 fmpz_mod_poly_zero(S, ctx);
140 fmpz_mod_poly_zero(T, ctx);
141 }
142 else if (lenB == 0) /* lenA > lenB = 0 */
143 {
144 fmpz_invmod(inv, fmpz_mod_poly_lead(A, ctx), p);
145 fmpz_mod_poly_scalar_mul_fmpz(G, A, inv, ctx);
146 fmpz_mod_poly_zero(T, ctx);
147 fmpz_mod_poly_set_fmpz(S, inv, ctx);
148 }
149 else /* lenA >= lenB >= 1 */
150 {
151 fmpz *g, *s, *t;
152 slong lenG;
153
154 if (G == A || G == B)
155 {
156 g = _fmpz_vec_init(FLINT_MIN(lenA, lenB));
157 }
158 else
159 {
160 fmpz_mod_poly_fit_length(G, FLINT_MIN(lenA, lenB), ctx);
161 g = G->coeffs;
162 }
163 if (S == A || S == B)
164 {
165 s = _fmpz_vec_init(lenB);
166 }
167 else
168 {
169 fmpz_mod_poly_fit_length(S, lenB, ctx);
170 s = S->coeffs;
171 }
172 if (T == A || T == B)
173 {
174 t = _fmpz_vec_init(lenA);
175 }
176 else
177 {
178 fmpz_mod_poly_fit_length(T, lenA, ctx);
179 t = T->coeffs;
180 }
181
182 fmpz_invmod(inv, fmpz_mod_poly_lead(B, ctx), p);
183 lenG = _fmpz_mod_poly_xgcd_euclidean(g, s, t,
184 A->coeffs, lenA, B->coeffs, lenB, inv, p);
185
186 if (G == A || G == B)
187 {
188 _fmpz_vec_clear(G->coeffs, G->alloc);
189 G->coeffs = g;
190 G->alloc = FLINT_MIN(lenA, lenB);
191 }
192 if (S == A || S == B)
193 {
194 _fmpz_vec_clear(S->coeffs, S->alloc);
195 S->coeffs = s;
196 S->alloc = lenB;
197 }
198 if (T == A || T == B)
199 {
200 _fmpz_vec_clear(T->coeffs, T->alloc);
201 T->coeffs = t;
202 T->alloc = lenA;
203 }
204
205 _fmpz_mod_poly_set_length(G, lenG);
206 _fmpz_mod_poly_set_length(S, FLINT_MAX(lenB - lenG, 1));
207 _fmpz_mod_poly_set_length(T, FLINT_MAX(lenA - lenG, 1));
208 _fmpz_mod_poly_normalise(S);
209 _fmpz_mod_poly_normalise(T);
210
211 if (!fmpz_is_one(fmpz_mod_poly_lead(G, ctx)))
212 {
213 fmpz_invmod(inv, fmpz_mod_poly_lead(G, ctx), p);
214 fmpz_mod_poly_scalar_mul_fmpz(G, G, inv, ctx);
215 fmpz_mod_poly_scalar_mul_fmpz(S, S, inv, ctx);
216 fmpz_mod_poly_scalar_mul_fmpz(T, T, inv, ctx);
217 }
218 }
219 fmpz_clear(inv);
220 }
221 }
222
223