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