1 /*
2 Copyright (C) 2011 Sebastian Pancratz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <stdlib.h>
13 #include "flint.h"
14 #include "fmpz.h"
15 #include "fmpz_vec.h"
16 #include "fmpz_poly.h"
17 #include "fmpq_poly.h"
18
_fmpq_poly_xgcd(fmpz * G,fmpz_t denG,fmpz * S,fmpz_t denS,fmpz * T,fmpz_t denT,const fmpz * A,const fmpz_t denA,slong lenA,const fmpz * B,const fmpz_t denB,slong lenB)19 void _fmpq_poly_xgcd(fmpz *G, fmpz_t denG,
20 fmpz *S, fmpz_t denS, fmpz *T, fmpz_t denT,
21 const fmpz *A, const fmpz_t denA, slong lenA,
22 const fmpz *B, const fmpz_t denB, slong lenB)
23 {
24 int alloc = 0;
25 fmpz *primA, *primB, *C, *D;
26 fmpz_t cA, cB;
27 slong lenG, lenC, lenD;
28
29 fmpz_init(cA);
30 fmpz_init(cB);
31
32 _fmpz_vec_content(cA, A, lenA);
33 _fmpz_vec_content(cB, B, lenB);
34
35 if (fmpz_is_one(cA))
36 {
37 if (fmpz_is_one(cB))
38 {
39 primA = (fmpz *) A;
40 primB = (fmpz *) B;
41 }
42 else
43 {
44 alloc |= 1;
45 primA = (fmpz *) A;
46 primB = _fmpz_vec_init(lenB);
47 _fmpz_vec_scalar_divexact_fmpz(primB, B, lenB, cB);
48 }
49 }
50 else
51 {
52 if (fmpz_is_one(cA))
53 {
54 alloc |= 2;
55 primA = _fmpz_vec_init(lenA);
56 primB = (fmpz *) B;
57 _fmpz_vec_scalar_divexact_fmpz(primA, A, lenA, cA);
58 }
59 else
60 {
61 alloc |= 3;
62 primA = _fmpz_vec_init(lenA + lenB);
63 primB = primA + lenA;
64 _fmpz_vec_scalar_divexact_fmpz(primA, A, lenA, cA);
65 _fmpz_vec_scalar_divexact_fmpz(primB, B, lenB, cB);
66 }
67 }
68
69 _fmpz_poly_gcd(G, primA, lenA, primB, lenB);
70
71 for (lenG = lenB - 1; !G[lenG]; lenG--) ;
72 lenG++;
73
74 if (lenG > 1)
75 {
76 alloc |= 4;
77 lenC = lenA - lenG + 1;
78 lenD = lenB - lenG + 1;
79 C = _fmpz_vec_init(lenC + lenD);
80 D = C + lenC;
81 _fmpz_poly_div(C, primA, lenA, G, lenG, 0);
82 _fmpz_poly_div(D, primB, lenB, G, lenG, 0);
83 }
84 else
85 {
86 lenC = lenA;
87 lenD = lenB;
88 C = primA;
89 D = primB;
90 }
91
92 _fmpz_poly_xgcd(denG, S, T, C, lenC, D, lenD);
93
94 if (!fmpz_is_one(denA))
95 _fmpz_vec_scalar_mul_fmpz(S, S, lenD, denA);
96 fmpz_mul(cA, cA, denG);
97 fmpz_mul(denS, cA, G + (lenG - 1));
98
99 if (!fmpz_is_one(denB))
100 _fmpz_vec_scalar_mul_fmpz(T, T, lenC, denB);
101 fmpz_mul(cB, cB, denG);
102 fmpz_mul(denT, cB, G + (lenG - 1));
103
104 _fmpz_vec_zero(S + lenD, lenB - lenD);
105 _fmpz_vec_zero(T + lenC, lenA - lenC);
106
107 _fmpq_poly_canonicalise(S, denS, lenD);
108 _fmpq_poly_canonicalise(T, denT, lenC);
109
110 fmpz_set(denG, G + (lenG - 1));
111
112 if ((alloc & 3) == 1)
113 _fmpz_vec_clear(primB, lenB);
114 else if ((alloc & 3) == 2)
115 _fmpz_vec_clear(primA, lenA);
116 else if ((alloc & 3) == 3)
117 _fmpz_vec_clear(primA, lenA + lenB);
118
119 if ((alloc & 4))
120 _fmpz_vec_clear(C, lenC + lenD);
121
122 fmpz_clear(cA);
123 fmpz_clear(cB);
124 }
125
fmpq_poly_xgcd(fmpq_poly_t G,fmpq_poly_t S,fmpq_poly_t T,const fmpq_poly_t A,const fmpq_poly_t B)126 void fmpq_poly_xgcd(fmpq_poly_t G, fmpq_poly_t S, fmpq_poly_t T,
127 const fmpq_poly_t A, const fmpq_poly_t B)
128 {
129 if (G == S || G == T || S == T)
130 {
131 flint_printf("Exception (fmpq_poly_xgcd). Output arguments aliased.\n");
132 flint_abort();
133 }
134
135 if (A->length < B->length)
136 {
137 fmpq_poly_xgcd(G, T, S, B, A);
138 } else
139 {
140 slong lenA = A->length, lenB = B->length, lenG = lenB;
141
142 if (lenA == 0) /* lenA = lenB = 0 */
143 {
144 fmpq_poly_zero(G);
145 fmpq_poly_zero(S);
146 fmpq_poly_zero(T);
147 }
148 else if (lenB == 0) /* lenA > lenB = 0 */
149 {
150 fmpq_poly_make_monic(G, A);
151 fmpq_poly_zero(T);
152 fmpq_poly_fit_length(S, 1);
153 _fmpq_poly_set_length(S, 1);
154 if (fmpz_sgn(A->coeffs + (lenA - 1)) > 0)
155 {
156 fmpz_set(S->coeffs, A->den);
157 fmpz_set(S->den, A->coeffs + (lenA - 1));
158 }
159 else
160 {
161 fmpz_neg(S->coeffs, A->den);
162 fmpz_neg(S->den, A->coeffs + (lenA - 1));
163 }
164 fmpq_poly_canonicalise(S);
165 }
166 else if (lenB == 1) /* lenA >= lenB = 1 */
167 {
168 fmpq_poly_set_ui(G, 1);
169 fmpq_poly_zero(S);
170 fmpq_poly_fit_length(T, 1);
171 _fmpq_poly_set_length(T, 1);
172 if (fmpz_sgn(B->coeffs) > 0)
173 {
174 fmpz_set(T->coeffs, B->den);
175 fmpz_set(T->den, B->coeffs);
176 }
177 else
178 {
179 fmpz_neg(T->coeffs, B->den);
180 fmpz_neg(T->den, B->coeffs);
181 }
182 }
183 else /* lenA >= lenB >= 2 */
184 {
185 /* Aliasing */
186 if (G == A || G == B)
187 {
188 fmpq_poly_t tG;
189
190 fmpq_poly_init2(tG, lenG);
191 fmpq_poly_xgcd(tG, S, T, A, B);
192 fmpq_poly_swap(tG, G);
193 fmpq_poly_clear(tG);
194 }
195 else if (S == A || S == B)
196 {
197 fmpq_poly_t tS;
198
199 fmpq_poly_init2(tS, lenB);
200 fmpq_poly_xgcd(G, tS, T, A, B);
201 fmpq_poly_swap(tS, S);
202 fmpq_poly_clear(tS);
203 }
204 else if (T == A || T == B)
205 {
206 fmpq_poly_t tT;
207
208 fmpq_poly_init2(tT, lenA);
209 fmpq_poly_xgcd(G, S, tT, A, B);
210 fmpq_poly_swap(tT, T);
211 fmpq_poly_clear(tT);
212 }
213 else /* no aliasing */
214 {
215
216 fmpq_poly_fit_length(G, lenG);
217 fmpq_poly_fit_length(S, lenB);
218 fmpq_poly_fit_length(T, lenA);
219
220 _fmpq_poly_xgcd(G->coeffs, G->den, S->coeffs, S->den, T->coeffs, T->den,
221 A->coeffs, A->den, lenA, B->coeffs, B->den, lenB);
222
223 _fmpq_poly_set_length(G, lenG);
224 _fmpq_poly_set_length(S, lenB);
225 _fmpq_poly_set_length(T, lenA);
226
227 _fmpq_poly_normalise(G);
228 _fmpq_poly_normalise(S);
229 _fmpq_poly_normalise(T);
230 }
231 }
232 }
233 }
234
235