source: git/factory/facNTLzzpEXGCD.cc @ 70c40f

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