My Project
Loading...
Searching...
No Matches
cfNTLzzpEXGCD.cc
Go to the documentation of this file.
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
20NTL -- A Library for Doing Number Theory
21Copyright (C) 1996-2009 Victor Shoup
22
23The most recent version of NTL is available at http://www.shoup.net
24
25This program is free software; you can redistribute it and/or
26modify it under the terms of the GNU General Public License
27as published by the Free Software Foundation; either version 2
28of the License, or (at your option) any later version.
29
30This program is distributed in the hope that it will be useful,
31but WITHOUT ANY WARRANTY; without even the implied warranty of
32MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33GNU General Public License for more details.
34
35You should have received a copy of the GNU General Public License
36along with this program; if not, write to the Free Software
37Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
38
39This entire copyright notice should be placed in an appropriately
40conspicuous place accompanying all distributions of software that
41make use of NTL.
42
43The above terms apply to all of the software modules distributed with NTL,
44i.e., all source files in either the ntl-xxx.tar.gz or WinNTL-xxx.zip
45distributions. In general, the individual files do not contain
46copyright notices.
47
48Note that the quad_float package is derived from the doubledouble package,
49originally developed by Keith Briggs, and also licensed unger the GNU GPL.
50The files quad_float.c and quad_float.h contain more detailed copyright
51notices.
52
53Note that the traditional long integer package used by NTL, lip.c, is derived
54from---and represents an extensive modification of---
55a package originally developed and copyrighted by Arjen Lenstra,
56who has agreed to renounce any copyright claims on the particular
57version of the long integer package appearing in NTL, so that the
58this package now is covered by the GNU GPL as well.
59
60Note that the alternative long integer package used by NTL is GMP,
61which is written by Torbjorn Granlund <tege@swox.com>.
62GMP is licensed under the terms of the GNU Lesser General Public License.
63
64Note that NTL makes use of the RSA Data Security, Inc. MD5 Message
65Digest Algorithm.
66
67Note 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
75The hope is that the GNU GPL is actually less restrictive than these
76older terms; however, in any circumstances such that GNU GPL is more
77restrictive, then the following rule is in force:
78versions prior to 4.0 may continue to be used under the old terms,
79but 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
95long 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
100static
101void 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
109void 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
167void 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
242void 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
281void 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
Conversion to and from NTL.
#define swap(_i, _j)
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
CanonicalForm b
Definition: cfModGcd.cc:4103
static void SetSize(vec_zz_pX &x, long n, long m)
void tryNTLGCD(zz_pEX &x, const zz_pEX &a, const zz_pEX &b, bool &fail)
compute the GCD x of a and b, fail is set to true if a zero divisor is encountered
long InvModStatus(zz_pE &x, const zz_pE &a)
void tryNTLXGCD(zz_pEX &d, zz_pEX &s, zz_pEX &t, const zz_pEX &a, const zz_pEX &b, bool &fail)
compute the extended GCD d=s*a+t*b, fail is set to true if a zero divisor is encountered
void tryPlainRem(zz_pEX &r, const zz_pEX &a, const zz_pEX &b, vec_zz_pX &x, bool &fail)
void tryPlainDivRem(zz_pEX &q, zz_pEX &r, const zz_pEX &a, const zz_pEX &b, bool &fail)
This file defines functions for univariate GCD and extended GCD over Z/p[t]/(f)[x] for reducible f.
const CanonicalForm int s
Definition: facAbsFact.cc:51
CFList conv(const CFFList &L)
convert a CFFList to a CFList by dropping the multiplicity
Definition: facBivar.cc:126
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
static int max(int a, int b)
Definition: fast_mult.cc:264
static BOOLEAN IsOne(number a, const coeffs)
Definition: flintcf_Q.cc:332
static BOOLEAN IsZero(number a, const coeffs)
Definition: flintcf_Q.cc:328
STATIC_VAR unsigned add[]
Definition: misc_ip.cc:107