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 
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 
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
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 
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 
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 
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 
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 
const CanonicalForm int s
Definition: facAbsFact.cc:55
int j
Definition: facHensel.cc:105
Conversion to and from NTL.
CFList conv(const CFFList &L)
convert a CFFList to a CFList by dropping the multiplicity
Definition: facBivar.cc:126
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 ...
CanonicalForm b
Definition: cfModGcd.cc:4044
static int max(int a, int b)
Definition: fast_mult.cc:264
int m
Definition: cfEzgcd.cc:121
int i
Definition: cfEzgcd.cc:125
This file defines functions for univariate GCD and extended GCD over Z/p[t]/(f)[x] for reducible f...
void tryPlainDivRem(zz_pEX &q, zz_pEX &r, const zz_pEX &a, const zz_pEX &b, bool &fail)
static BOOLEAN IsOne(number a, const coeffs r)
Definition: flintcf_Q.cc:330
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
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 ...
#define swap(_i, _j)
void tryPlainRem(zz_pEX &r, const zz_pEX &a, const zz_pEX &b, vec_zz_pX &x, bool &fail)
Variable x
Definition: cfModGcd.cc:4023
STATIC_VAR unsigned add[]
Definition: misc_ip.cc:117
static BOOLEAN IsZero(number a, const coeffs r)
Definition: flintcf_Q.cc:326
long InvModStatus(zz_pE &x, const zz_pE &a)
int degree(const CanonicalForm &f)
static void SetSize(vec_zz_pX &x, long n, long m)