1 /*
2 Copyright (C) 2020 Daniel Schultz
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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "n_poly.h"
13
14
_n_fq_poly_gcd_euclidean_inplace_(mp_limb_t * A,slong Alen,mp_limb_t * B,slong Blen,const fq_nmod_ctx_t ctx,mp_limb_t * tmp)15 slong _n_fq_poly_gcd_euclidean_inplace_(
16 mp_limb_t * A, slong Alen,
17 mp_limb_t * B, slong Blen,
18 const fq_nmod_ctx_t ctx,
19 mp_limb_t * tmp)
20 {
21 slong d = fq_nmod_ctx_degree(ctx);
22 nmod_t mod = fq_nmod_ctx_mod(ctx);
23 slong i;
24 mp_limb_t * u = tmp;
25 mp_limb_t * q0 = u + d;
26 mp_limb_t * q1 = q0 + d;
27 mp_limb_t * t = q1 + d;
28
29 again:
30
31 if (Alen < 2 || Blen < 2)
32 {
33 if (Alen < 1)
34 {
35 if (Blen < 1)
36 return 0;
37
38 _n_fq_inv(u, B + d*(Blen - 1), ctx, t);
39 for (i = 0; i + 1 < Blen; i++)
40 _n_fq_mul(B + d*i, B + d*i, u, ctx, t);
41 _n_fq_one(B + d*(Blen - 1), d);
42 return -Blen - 1;
43 }
44
45 if (Blen < 1)
46 {
47 _n_fq_inv(u, A + d*(Alen - 1), ctx, t);
48 for (i = 0; i + 1 < Alen; i++)
49 _n_fq_mul(A + d*i, A + d*i, u, ctx, t);
50 _n_fq_one(A + d*(Alen - 1), d);
51 return Alen;
52 }
53
54 if (Blen < 2)
55 {
56 _n_fq_one(B + d*0, d);
57 return -1 - 1;
58 }
59
60 FLINT_ASSERT(Alen < 2);
61
62 _n_fq_one(A + d*0, d);
63 return 1;
64 }
65
66 if (Alen > Blen)
67 {
68 /* Q = A[Alen-1]/B[Blen-1] x + */
69 _n_fq_inv(u, B + d*(Blen - 1), ctx, t);
70 _n_fq_mul(q1, A + d*(Alen - 1), u, ctx, t);
71 _n_fq_mul(q0, q1, B + d*(Blen - 2), ctx, t);
72 _n_fq_sub(q0, q0, A + d*(Alen - 2), d, mod);
73 _n_fq_mul(q0, q0, u, ctx, t);
74
75 _nmod_vec_neg(q1, q1, d, mod);
76
77 _n_fq_mul(u, q0, B + d*0, ctx, t);
78 _n_fq_add(A + d*(-1 + Alen - Blen),
79 A + d*(-1 + Alen - Blen), u, d, mod);
80
81 for (i = 0; i < Blen - 1; i++)
82 {
83 _n_fq_mul2(t, q1, B + d*i, ctx);
84 _n_fq_madd2(t, q0, B + d*(i + 1), ctx, t + 2*d);
85 _n_fq_reduce2(u, t, ctx, t + 2*d);
86 _n_fq_add(A + d*(i + Alen - Blen),
87 A + d*(i + Alen - Blen), u, d, mod);
88 }
89
90 Alen -= 2;
91 while (Alen > 0 && _n_fq_is_zero(A + d*(Alen - 1), d))
92 Alen--;
93
94 goto again;
95 }
96 else if (Blen > Alen)
97 {
98 _n_fq_inv(u, A + d*(Alen - 1), ctx, t);
99 _n_fq_mul(q1, B + d*(Blen - 1), u, ctx, t);
100 _n_fq_mul(q0, q1, A + d*(Alen - 2), ctx, t);
101 _n_fq_sub(q0, q0, B + d*(Blen - 2), d, mod);
102 _n_fq_mul(q0, q0, u, ctx, t);
103
104 _nmod_vec_neg(q1, q1, d, mod);
105
106 i = -1;
107 _n_fq_mul(u, q0, A + d*(i + 1), ctx, t);
108 _n_fq_add(B + d*(i + Blen - Alen), B + d*(i + Blen - Alen), u, d, mod);
109
110 for (i = 0; i < Alen - 2; i++)
111 {
112 _n_fq_mul2(t, q1, A + d*i, ctx);
113 _n_fq_madd2(t, q0, A + d*(i + 1), ctx, t + 2*d);
114 _n_fq_reduce2(u, t, ctx, t + 2*d);
115 _n_fq_add(B + d*(i + Blen - Alen), B + d*(i + Blen - Alen), u, d, mod);
116 }
117
118 Blen -= 2;
119 while (Blen > 0 && _n_fq_is_zero(B + d*(Blen - 1), d))
120 Blen--;
121
122 goto again;
123 }
124 else
125 {
126 _n_fq_inv(u, B + d*(Blen - 1), ctx, t);
127 _n_fq_mul(q0, A + d*(Alen - 1), u, ctx, t);
128
129 for (i = 0; i < Blen - 1; i++)
130 {
131 _n_fq_mul(u, q0, B + d*i, ctx, t);
132 _n_fq_sub(A + d*i, A + d*i, u, d, mod);
133 }
134
135 Alen -= 1;
136 while (Alen > 0 && _n_fq_is_zero(A + d*(Alen - 1), d))
137 Alen--;
138
139 goto again;
140 }
141 }
142
143
n_fq_poly_gcd_(n_fq_poly_t G,const n_fq_poly_t A,const n_fq_poly_t B,const fq_nmod_ctx_t ctx,n_poly_stack_t St)144 void n_fq_poly_gcd_(
145 n_fq_poly_t G,
146 const n_fq_poly_t A,
147 const n_fq_poly_t B,
148 const fq_nmod_ctx_t ctx,
149 n_poly_stack_t St)
150 {
151 slong d = fq_nmod_ctx_degree(ctx);
152 slong n;
153 mp_limb_t * a, * b, * t;
154 #if FLINT_WANT_ASSERT
155 fq_nmod_poly_t GG, AA, BB;
156 fq_nmod_poly_init(GG, ctx);
157 fq_nmod_poly_init(AA, ctx);
158 fq_nmod_poly_init(BB, ctx);
159 FLINT_ASSERT(n_fq_poly_is_canonical(A, ctx));
160 FLINT_ASSERT(n_fq_poly_is_canonical(B, ctx));
161 n_fq_poly_get_fq_nmod_poly(AA, A, ctx);
162 n_fq_poly_get_fq_nmod_poly(BB, B, ctx);
163 fq_nmod_poly_gcd(GG, AA, BB, ctx);
164 #endif
165
166 n_poly_stack_fit_request(St, 3);
167
168 t = n_poly_stack_vec_init(St, 8*d);
169 a = n_poly_stack_vec_init(St, d*A->length + 1);
170 b = n_poly_stack_vec_init(St, d*B->length + 1);
171
172 _nmod_vec_set(a, A->coeffs, d*A->length);
173 _nmod_vec_set(b, B->coeffs, d*B->length);
174
175 n = _n_fq_poly_gcd_euclidean_inplace_(a, A->length, b, B->length, ctx, t);
176
177 if (n < 0)
178 {
179 n = -n - 1;
180 n_poly_fit_length(G, d*n);
181 _nmod_vec_set(G->coeffs, b, d*n);
182 G->length = n;
183 }
184 else
185 {
186 n_poly_fit_length(G, d*n);
187 _nmod_vec_set(G->coeffs, a, d*n);
188 G->length = n;
189 }
190
191 n_poly_stack_vec_clear(St);
192 n_poly_stack_vec_clear(St);
193 n_poly_stack_vec_clear(St);
194
195 #if FLINT_WANT_ASSERT
196 FLINT_ASSERT(n_fq_poly_is_canonical(G, ctx));
197 n_fq_poly_get_fq_nmod_poly(AA, G, ctx);
198 FLINT_ASSERT(fq_nmod_poly_equal(AA, GG, ctx));
199 fq_nmod_poly_clear(GG, ctx);
200 fq_nmod_poly_clear(AA, ctx);
201 fq_nmod_poly_clear(BB, ctx);
202 #endif
203 }
204
n_fq_poly_gcd(n_fq_poly_t G,const n_fq_poly_t A,const n_fq_poly_t B,const fq_nmod_ctx_t ctx)205 void n_fq_poly_gcd(
206 n_fq_poly_t G,
207 const n_fq_poly_t A,
208 const n_fq_poly_t B,
209 const fq_nmod_ctx_t ctx)
210 {
211 #if 0
212 n_poly_stack_t St;
213 n_poly_stack_init(St);
214 n_fq_poly_gcd_(G, A, B, ctx, St);
215 n_poly_stack_clear(St);
216 #else
217 fq_nmod_poly_t g, a, b;
218 fq_nmod_poly_init(g, ctx);
219 fq_nmod_poly_init(a, ctx);
220 fq_nmod_poly_init(b, ctx);
221 n_fq_poly_get_fq_nmod_poly(a, A, ctx);
222 n_fq_poly_get_fq_nmod_poly(b, B, ctx);
223 fq_nmod_poly_gcd(g, a, b, ctx);
224 n_fq_poly_set_fq_nmod_poly(G, g, ctx);
225 fq_nmod_poly_clear(g, ctx);
226 fq_nmod_poly_clear(a, ctx);
227 fq_nmod_poly_clear(b, ctx);
228 #endif
229 }
230
231