1 // -*- c++ -*-
2 //*****************************************************************************
3 /** @file cfNTLzzpEXGCD.cc
4 *
5 * @author Martin Lee
6 *
7 * Univariate GCD and extended GCD over Z/p[t]/(f)[x] for reducible f
8 *
9 * @note the following code is slightly modified code out of
10 * lzz_pEX.c from Victor Shoup's NTL. Below is NTL's copyright notice.
11 *
12 * @par Copyright:
13 * (c) by The SINGULAR Team, see LICENSE file
14 *
15
16 COPYRIGHT NOTICE
17 for NTL 5.5
18 (modified for Singular 2-0-6 - 3-1)
19
20 NTL -- A Library for Doing Number Theory
21 Copyright (C) 1996-2009 Victor Shoup
22
23 The most recent version of NTL is available at http://www.shoup.net
24
25 This program is free software; you can redistribute it and/or
26 modify it under the terms of the GNU General Public License
27 as published by the Free Software Foundation; either version 2
28 of the License, or (at your option) any later version.
29
30 This program is distributed in the hope that it will be useful,
31 but WITHOUT ANY WARRANTY; without even the implied warranty of
32 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33 GNU General Public License for more details.
34
35 You should have received a copy of the GNU General Public License
36 along with this program; if not, write to the Free Software
37 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
38
39 This entire copyright notice should be placed in an appropriately
40 conspicuous place accompanying all distributions of software that
41 make use of NTL.
42
43 The above terms apply to all of the software modules distributed with NTL,
44 i.e., all source files in either the ntl-xxx.tar.gz or WinNTL-xxx.zip
45 distributions. In general, the individual files do not contain
46 copyright notices.
47
48 Note that the quad_float package is derived from the doubledouble package,
49 originally developed by Keith Briggs, and also licensed unger the GNU GPL.
50 The files quad_float.c and quad_float.h contain more detailed copyright
51 notices.
52
53 Note that the traditional long integer package used by NTL, lip.c, is derived
54 from---and represents an extensive modification of---
55 a package originally developed and copyrighted by Arjen Lenstra,
56 who has agreed to renounce any copyright claims on the particular
57 version of the long integer package appearing in NTL, so that the
58 this package now is covered by the GNU GPL as well.
59
60 Note that the alternative long integer package used by NTL is GMP,
61 which is written by Torbjorn Granlund <tege@swox.com>.
62 GMP is licensed under the terms of the GNU Lesser General Public License.
63
64 Note that NTL makes use of the RSA Data Security, Inc. MD5 Message
65 Digest Algorithm.
66
67 Note that prior to version 4.0, NTL was distributed under the following terms:
68 NTL is freely available for research and educational purposes.
69 I don't want to attach any legalistic licensing restrictions on
70 users of NTL.
71 However, NTL should not be linked in a commercial program
72 (although using data in a commercial
73 product produced by a program that used NTL is fine).
74
75 The hope is that the GNU GPL is actually less restrictive than these
76 older terms; however, in any circumstances such that GNU GPL is more
77 restrictive, then the following rule is in force:
78 versions prior to 4.0 may continue to be used under the old terms,
79 but users of versions 4.0 or later should adhere to the terms of the GNU GPL.
80 **/
81
82
83
84 #include "config.h"
85
86
87 #include "cfNTLzzpEXGCD.h"
88
89 #ifdef HAVE_NTL
90 #include "NTLconvert.h"
91 #endif
92
93 #ifdef HAVE_NTL
94
InvModStatus(zz_pE & x,const zz_pE & a)95 long InvModStatus (zz_pE& x, const zz_pE& a)
96 {
97 return InvModStatus(x._zz_pE__rep, a._zz_pE__rep, zz_pE::modulus());
98 }
99
100 static
SetSize(vec_zz_pX & x,long n,long m)101 void SetSize(vec_zz_pX& x, long n, long m)
102 {
103 x.SetLength(n);
104 long i;
105 for (i = 0; i < n; i++)
106 x[i].rep.SetMaxLength(m);
107 }
108
tryPlainRem(zz_pEX & r,const zz_pEX & a,const zz_pEX & b,vec_zz_pX & x,bool & fail)109 void tryPlainRem(zz_pEX& r, const zz_pEX& a, const zz_pEX& b, vec_zz_pX& x,
110 bool& fail)
111 {
112 long da, db, dq, i, j, LCIsOne;
113 const zz_pE *bp;
114 zz_pX *xp;
115
116
117 zz_pE LCInv, t;
118 zz_pX s;
119
120 da = deg(a);
121 db = deg(b);
122
123 if (db < 0) Error("zz_pEX: division by zero");
124
125 if (da < db) {
126 r = a;
127 return;
128 }
129
130 bp = b.rep.elts();
131
132 if (IsOne(bp[db]))
133 LCIsOne = 1;
134 else {
135 LCIsOne = 0;
136 fail= InvModStatus (LCInv, bp[db]);
137 //inv(LCInv, bp[db]);
138 if (fail)
139 return;
140 }
141
142 for (i = 0; i <= da; i++)
143 x[i] = rep(a.rep[i]);
144
145 xp = x.elts();
146
147 dq = da - db;
148
149 for (i = dq; i >= 0; i--) {
150 conv(t, xp[i+db]);
151 if (!LCIsOne)
152 mul(t, t, LCInv);
153 NTL::negate(t, t);
154
155 for (j = db-1; j >= 0; j--) {
156 mul(s, rep(t), rep(bp[j]));
157 add(xp[i+j], xp[i+j], s);
158 }
159 }
160
161 r.rep.SetLength(db);
162 for (i = 0; i < db; i++)
163 conv(r.rep[i], xp[i]);
164 r.normalize();
165 }
166
tryPlainDivRem(zz_pEX & q,zz_pEX & r,const zz_pEX & a,const zz_pEX & b,bool & fail)167 void tryPlainDivRem(zz_pEX& q, zz_pEX& r, const zz_pEX& a, const zz_pEX& b,
168 bool& fail)
169 {
170 long da, db, dq, i, j, LCIsOne;
171 const zz_pE *bp;
172 zz_pE *qp;
173 zz_pX *xp;
174
175
176 zz_pE LCInv, t;
177 zz_pX s;
178
179 da = deg(a);
180 db = deg(b);
181
182 if (db < 0) Error("zz_pEX: division by zero");
183
184 if (da < db) {
185 r = a;
186 clear(q);
187 return;
188 }
189
190 zz_pEX lb;
191
192 if (&q == &b) {
193 lb = b;
194 bp = lb.rep.elts();
195 }
196 else
197 bp = b.rep.elts();
198
199 if (IsOne(bp[db]))
200 LCIsOne = 1;
201 else {
202 LCIsOne = 0;
203 fail= InvModStatus (LCInv, bp[db]);
204 //inv(LCInv, bp[db]);
205 if (fail)
206 return;
207 }
208
209 vec_zz_pX x;
210
211 SetSize(x, da+1, 2*zz_pE::degree());
212
213 for (i = 0; i <= da; i++)
214 x[i] = rep(a.rep[i]);
215
216 xp = x.elts();
217
218 dq = da - db;
219 q.rep.SetLength(dq+1);
220 qp = q.rep.elts();
221
222 for (i = dq; i >= 0; i--) {
223 conv(t, xp[i+db]);
224 if (!LCIsOne)
225 mul(t, t, LCInv);
226 qp[i] = t;
227 NTL::negate(t, t);
228
229 for (j = db-1; j >= 0; j--) {
230 mul(s, rep(t), rep(bp[j]));
231 add(xp[i+j], xp[i+j], s);
232 }
233 }
234
235 r.rep.SetLength(db);
236 for (i = 0; i < db; i++)
237 conv(r.rep[i], xp[i]);
238 r.normalize();
239 }
240
241
tryNTLGCD(zz_pEX & x,const zz_pEX & a,const zz_pEX & b,bool & fail)242 void tryNTLGCD(zz_pEX& x, const zz_pEX& a, const zz_pEX& b, bool& fail)
243 {
244 zz_pE t;
245
246 if (IsZero(b))
247 x = a;
248 else if (IsZero(a))
249 x = b;
250 else {
251 long n = max(deg(a),deg(b)) + 1;
252 zz_pEX u(INIT_SIZE, n), v(INIT_SIZE, n);
253
254 vec_zz_pX tmp;
255 SetSize(tmp, n, 2*zz_pE::degree());
256
257 u = a;
258 v = b;
259 do {
260 tryPlainRem(u, u, v, tmp, fail);
261 if (fail)
262 return;
263 swap(u, v);
264 } while (!IsZero(v));
265
266 x = u;
267 }
268
269 if (IsZero(x)) return;
270 if (IsOne(LeadCoeff(x))) return;
271
272 /* make gcd monic */
273
274
275 fail= InvModStatus (t, LeadCoeff(x));
276 if (fail)
277 return;
278 mul(x, x, t);
279 }
280
tryNTLXGCD(zz_pEX & d,zz_pEX & s,zz_pEX & t,const zz_pEX & a,const zz_pEX & b,bool & fail)281 void tryNTLXGCD(zz_pEX& d, zz_pEX& s, zz_pEX& t, const zz_pEX& a,
282 const zz_pEX& b, bool& fail)
283 {
284 zz_pE z;
285
286
287 if (IsZero(b)) {
288 set(s);
289 clear(t);
290 d = a;
291 }
292 else if (IsZero(a)) {
293 clear(s);
294 set(t);
295 d = b;
296 }
297 else {
298 long e = max(deg(a), deg(b)) + 1;
299
300 zz_pEX temp(INIT_SIZE, e), u(INIT_SIZE, e), v(INIT_SIZE, e),
301 u0(INIT_SIZE, e), v0(INIT_SIZE, e),
302 u1(INIT_SIZE, e), v1(INIT_SIZE, e),
303 u2(INIT_SIZE, e), v2(INIT_SIZE, e), q(INIT_SIZE, e);
304
305
306 set(u1); clear(v1);
307 clear(u2); set(v2);
308 u = a; v = b;
309
310 do {
311 fail= InvModStatus (z, LeadCoeff(v));
312 if (fail)
313 return;
314 DivRem(q, u, u, v);
315 swap(u, v);
316 u0 = u2;
317 v0 = v2;
318 mul(temp, q, u2);
319 sub(u2, u1, temp);
320 mul(temp, q, v2);
321 sub(v2, v1, temp);
322 u1 = u0;
323 v1 = v0;
324 } while (!IsZero(v));
325
326 d = u;
327 s = u1;
328 t = v1;
329 }
330
331 if (IsZero(d)) return;
332 if (IsOne(LeadCoeff(d))) return;
333
334 /* make gcd monic */
335
336 fail= InvModStatus (z, LeadCoeff(d));
337 if (fail)
338 return;
339 mul(d, d, z);
340 mul(s, s, z);
341 mul(t, t, z);
342 }
343 #endif
344
345