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