1 /*
2 Copyright (C) 2011 William Hart
3 Copyright (C) 2011, 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 "nmod_poly.h"
15 #include "mpn_extras.h"
16
17 #define __set(B, lenB, A, lenA) \
18 do { \
19 _nmod_vec_set((B), (A), (lenA)); \
20 (lenB) = (lenA); \
21 } while (0)
22
23 #define __rem(R, lenR, A, lenA, B, lenB) \
24 do { \
25 if ((lenA) >= (lenB)) \
26 { \
27 _nmod_poly_rem((R), (A), (lenA), (B), (lenB), mod); \
28 (lenR) = (lenB) - 1; \
29 MPN_NORM((R), (lenR)); \
30 } \
31 else \
32 { \
33 _nmod_vec_set((R), (A), (lenA)); \
34 (lenR) = (lenA); \
35 } \
36 } while (0)
37
38 /*
39 XXX: Incidentally, this implementation currently supports aliasing.
40 But since this may change in the future, no function other than
41 nmod_poly_gcd_hgcd() should rely on this.
42 */
43
_nmod_poly_gcd_hgcd(mp_ptr G,mp_srcptr A,slong lenA,mp_srcptr B,slong lenB,nmod_t mod)44 slong _nmod_poly_gcd_hgcd(mp_ptr G, mp_srcptr A, slong lenA,
45 mp_srcptr B, slong lenB, nmod_t mod)
46 {
47 const slong cutoff = FLINT_BIT_COUNT(mod.n) <= 8 ?
48 NMOD_POLY_SMALL_GCD_CUTOFF : NMOD_POLY_GCD_CUTOFF;
49
50 mp_ptr J = _nmod_vec_init(2 * lenB);
51 mp_ptr R = J + lenB;
52
53 slong lenG, lenJ, lenR;
54
55 __rem(R, lenR, A, lenA, B, lenB);
56
57 if (lenR == 0)
58 {
59 __set(G, lenG, B, lenB);
60 }
61 else
62 {
63 _nmod_poly_hgcd(NULL, NULL, G, &(lenG), J, &(lenJ), B, lenB, R, lenR, mod);
64
65 while (lenJ != 0)
66 {
67 __rem(R, lenR, G, lenG, J, lenJ);
68
69 if (lenR == 0)
70 {
71 __set(G, lenG, J, lenJ);
72 break;
73 }
74 if (lenJ < cutoff)
75 {
76 lenG = _nmod_poly_gcd_euclidean(G, J, lenJ, R, lenR, mod);
77 break;
78 }
79
80 _nmod_poly_hgcd(NULL, NULL, G, &(lenG), J, &(lenJ), J, lenJ, R, lenR, mod);
81 }
82 }
83 _nmod_vec_clear(J);
84
85 return lenG;
86 }
87
nmod_poly_gcd_hgcd(nmod_poly_t G,const nmod_poly_t A,const nmod_poly_t B)88 void nmod_poly_gcd_hgcd(nmod_poly_t G,
89 const nmod_poly_t A, const nmod_poly_t B)
90 {
91 if (A->length < B->length)
92 {
93 nmod_poly_gcd_hgcd(G, B, A);
94 }
95 else /* lenA >= lenB >= 0 */
96 {
97 slong lenA = A->length, lenB = B->length, lenG;
98 nmod_poly_t tG;
99 mp_ptr g;
100
101 if (lenA == 0) /* lenA = lenB = 0 */
102 {
103 nmod_poly_zero(G);
104 }
105 else if (lenB == 0) /* lenA > lenB = 0 */
106 {
107 nmod_poly_make_monic(G, A);
108 }
109 else /* lenA >= lenB >= 1 */
110 {
111 if (G == A || G == B)
112 {
113 nmod_poly_init2(tG, A->mod.n, FLINT_MIN(lenA, lenB));
114 g = tG->coeffs;
115 }
116 else
117 {
118 nmod_poly_fit_length(G, FLINT_MIN(lenA, lenB));
119 g = G->coeffs;
120 }
121
122 lenG = _nmod_poly_gcd_hgcd(g, A->coeffs, lenA,
123 B->coeffs, lenB, A->mod);
124
125 if (G == A || G == B)
126 {
127 nmod_poly_swap(tG, G);
128 nmod_poly_clear(tG);
129 }
130 G->length = lenG;
131
132 if (G->length == 1)
133 G->coeffs[0] = 1;
134 else
135 nmod_poly_make_monic(G, G);
136 }
137 }
138 }
139
140 #undef __set
141 #undef __rem
142
143