source: git/factory/libfac/factor/SqrFree.cc @ be5dff

spielwiese
Last change on this file since be5dff was be5dff, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
fix extra ';' error of libfac
  • Property mode set to 100644
File size: 13.1 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2// emacs edit mode for this file is -*- C++ -*-
3/* $Id$ */
4static const char * errmsg = "\nYou found a bug!\nPlease inform singular@mathematik.uni-kl.de\n Please include above information and your input (the ideal/polynomial and characteristic) in your bug-report.\nThank you.";
5///////////////////////////////////////////////////////////////////////////////
6// FACTORY - Includes
7
8#include<factory/factory.h>
9
10#ifndef NOSTREAMIO
11#ifdef HAVE_IOSTREAM
12#include <iostream>
13#define OSTREAM std::ostream
14#define ISTREAM std::istream
15#define CERR std::cerr
16#define CIN std::cin
17#elif defined(HAVE_IOSTREAM_H)
18#include <iostream.h>
19#define OSTREAM ostream
20#define ISTREAM istream
21#define CERR cerr
22#define CIN cin
23#endif
24#endif
25// Factor - Includes
26#include "tmpl_inst.h"
27#include "helpstuff.h"
28// some CC's need this:
29#include "SqrFree.h"
30
31#ifdef SQRFREEDEBUG
32# define DEBUGOUTPUT
33#else
34# undef DEBUGOUTPUT
35#endif
36
37#include <libfac/factor/debug.h>
38#include "timing.h"
39
40TIMING_DEFINE_PRINT(squarefree_time)
41TIMING_DEFINE_PRINT(gcd_time)
42
43static inline CFFactor
44Powerup( const CFFactor & F , int exp=1)
45{
46  return CFFactor(F.factor(), exp*F.exp()) ;
47}
48
49static CFFList
50Powerup( const CFFList & Inputlist , int exp=1 )
51{
52  CFFList Outputlist;
53
54  for ( CFFListIterator i=Inputlist; i.hasItem(); i++ )
55    Outputlist.append(Powerup(i.getItem(), exp));
56  return Outputlist ;
57}
58
59///////////////////////////////////////////////////////////////
60// Compute the Pth root of a polynomial in characteristic p  //
61// f must be a polynomial which we can take the Pth root of. //
62// Domain is q=p^m , f a uni/multivariate polynomial         //
63///////////////////////////////////////////////////////////////
64static CanonicalForm
65PthRoot( const CanonicalForm & f )
66{
67  CanonicalForm RES, R = f;
68  int n= max(level(R),getNumVars(R)), p= getCharacteristic();
69
70  if (n==0)
71  { // constant
72    if (R.inExtension()) // not in prime field; f over |F(q=p^k)
73    {
74      R = power(R,Powerup(p,getGFDegree() - 1)) ;
75    }
76    // if f in prime field, do nothing
77    return R;
78  }
79  // we assume R is a Pth power here
80  RES = R.genZero();
81  Variable x(n);
82  for (int i=0; i<= (int) (degree(R,level(R))/p) ; i++)
83    RES += PthRoot( R[i*p] ) * power(x,i);
84  return RES;
85}
86
87///////////////////////////////////////////////////////////////
88// Compute the Pth root of a polynomial in characteristic p  //
89// f must be a polynomial which we can take the Pth root of. //
90// Domain is q=p^m , f a uni/multivariate polynomial         //
91///////////////////////////////////////////////////////////////
92static CanonicalForm
93PthRoot( const CanonicalForm & f ,const CanonicalForm & mipo)
94{
95  CanonicalForm RES, R = f;
96  int n= max(level(R),getNumVars(R)), p= getCharacteristic();
97  int mipodeg=-1;
98  if (f.level()==mipo.level()) mipodeg=mipo.degree();
99  else if ((f.level()==1) &&(!mipo.isZero()))
100  {
101    Variable t;
102    CanonicalForm tt=getMipo(mipo.mvar(),t);
103    mipodeg=degree(tt,t);
104  }
105
106  if ((n==0)
107  ||(mipodeg!=-1))
108  { // constant
109    if (R.inExtension()) // not in prime field; f over |F(q=p^k)
110    {
111      R = power(R,Powerup(p,getGFDegree() - 1)) ;
112    }
113    else if ((f.level()==mipo.level())
114    ||((f.level()==1) &&(!mipo.isZero())))
115    {
116      R = power(R,Powerup(p,mipodeg - 1)) ;
117      R=mod(R, mipo);
118    }
119    // if f in prime field, do nothing
120    return R;
121  }
122  // we assume R is a Pth power here
123  RES = R.genZero();
124  Variable x(n);
125  for (int i=0; i<= (int) (degree(R,level(R))/p) ; i++)
126    RES += PthRoot( R[i*p], mipo ) * power(x,i);
127  return RES;
128}
129
130///////////////////////////////////////////////////////////////
131// A uni/multivariate SqrFreeTest routine.                   //
132// Cheaper to run if all you want is a test.                 //
133// Works for charcteristic 0 and q=p^m                       //
134// Returns 1 if poly r is SqrFree, 0 if SqrFree will do some //
135// kind of factorization.                                    //
136// Would be much more effcient iff we had *good*             //
137//  uni/multivariate gcd's and/or gcdtest's                  //
138///////////////////////////////////////////////////////////////
139int
140SqrFreeTest( const CanonicalForm & r, int opt)
141{
142  CanonicalForm f=r, g;
143  int n=level(f);
144
145  if (getNumVars(f)==0) return 1 ; // a constant is SqrFree
146  if ( f.isUnivariate() )
147  {
148    g= f.deriv();
149    if ( getCharacteristic() > 0 && g.isZero() ) return 0 ;
150    // Next: it would be best to have a *univariate* gcd-test which returns
151    // 0 iff gcdtest(f,g) == 1 or a constant ( for real Polynomials )
152    g = gcd(f,g);
153    if ( g.isOne() || (-g).isOne() ) return 1;
154    else
155      if ( getNumVars(g) == 0 ) return 1;// <- totaldegree!!!
156      else return 0 ;
157  }
158  else
159  { // multivariate case
160    for ( int k=1; k<=n; k++ )
161    {
162      g = swapvar(f,k,n); g = content(g);
163      // g = 1 || -1 : sqr-free, g poly : not sqr-free, g number : opt helps
164      if ( ! (g.isOne() || (-g).isOne() || getNumVars(g)==0 ) ) {
165        if ( opt==0 ) return 0;
166        else {
167          if ( SqrFreeTest(g,1) == 0 ) return 0;
168          g = swapvar(g,k,n);
169          f /=g ;
170        }
171      }
172    }
173    // Now f is primitive
174    n = level(f); // maybe less indeterminants
175    //    if ( totaldegree(f) <= 1 ) return 1;
176
177    // Let`s look if it is a Pth root
178    if ( getCharacteristic() > 0 )
179      for (int k=1; k<=n; k++ )
180      {
181        g=swapvar(f,k,n); g=g.deriv();
182        if ( ! g.isZero() ) break ;
183        else if ( k==n) return 0 ; // really is Pth root
184      }
185    g = f.deriv() ;
186    // Next: it would be best to have a *multivariate* gcd-test which returns
187    // 0 iff gcdtest(f,g) == 1 or a constant ( for real Polynomials )
188    g= gcd(f,g);
189    if ( g.isOne() || (-g).isOne() || (g==f) || (getNumVars(g)==0) ) return 1 ;
190    else return 0 ;
191  }
192  factoryError("libfac: ERROR: SqrFreeTest: we should never fall trough here!");
193  return 0;
194}
195
196///////////////////////////////////////////////////////////////
197// A uni/multivariate SqrFree routine.Works for polynomials  //
198// which don\'t have a constant as the content.              //
199// Works for charcteristic 0 and q=p^m                       //
200// returns a list of polys each of sqrfree, but gcd(f_i,f_j) //
201// needs not to be 1 !!!!!                                   //
202///////////////////////////////////////////////////////////////
203static CFFList
204SqrFreed( const CanonicalForm & r , const CanonicalForm &mipo=0)
205{
206  CanonicalForm h, g, f = r;
207  CFFList Outputlist;
208  int n = level(f);
209
210  DEBINCLEVEL(CERR, "SqrFreed");
211  DEBOUTLN(CERR, "Called with r= ", r);
212  if (getNumVars(f)==0 )
213  { // just a constant; return it
214    Outputlist= CFFactor(f,1);
215    return Outputlist ;
216  }
217
218// We look if we do have a content; if so, SqrFreed the content
219// and continue computations with pp(f)
220  for (int k=1; k<=n; k++)
221  {
222    if ((mipo.isZero())/*||(k!=1)*/)
223    {
224      g = swapvar(f,k,n); g = content(g);
225      if ( ! (g.isOne() || (-g).isOne() || (degree(g)==0) ))
226      {
227        g = swapvar(g,k,n);
228        DEBOUTLN(CERR, "We have a content: ", g);
229        Outputlist = myUnion(SqrFreeMV(g,mipo),Outputlist); // should we add a
230                                                // SqrFreeTest(g) first ?
231        DEBOUTLN(CERR, "Outputlist is now: ", Outputlist);
232        f /=g;
233        DEBOUTLN(CERR, "f is now: ", f);
234      }
235    }
236  }
237
238// Now f is primitive; Let`s look if f is univariate
239  if ( f.isUnivariate() )
240  {
241    DEBOUTLN(CERR, "f is univariate: ", f);
242    g = content(f);
243    if ( ! (g.isOne() || (-g).isOne() ) )
244    {
245      Outputlist= myappend(Outputlist,CFFactor(g,1)) ;
246      f /= g;
247    }
248    Outputlist = Union(sqrFree(f),Outputlist) ;
249    DEBOUTLN(CERR, "Outputlist after univ. sqrFree(f) = ", Outputlist);
250    DEBDECLEVEL(CERR, "SqrFreed");
251    return Outputlist ;
252  }
253
254// Linear?
255  if ( totaldegree(f) <= 1 )
256  {
257    Outputlist= myappend(Outputlist,CFFactor(f,1)) ;
258    DEBDECLEVEL(CERR, "SqrFreed");
259    return Outputlist ;
260  }
261
262// is it Pth root?
263  n=level(f); // maybe less indeterminants
264  g= f.deriv();
265  if ( getCharacteristic() > 0 && g.isZero() )
266  {  // Pth roots only apply to char > 0
267    for (int k=1; k<=n; k++)
268    {
269      if ((mipo.isZero())/*||(k!=1)*/)
270      {
271        g=swapvar(f,k,n) ;
272        g = g.deriv();
273
274        if ( ! g.isZero() )
275        { // can`t be Pth root
276          CFFList Outputlist2= SqrFreed(swapvar(f,k,n));
277          for (CFFListIterator inter=Outputlist2; inter.hasItem(); inter++)
278          {
279            Outputlist= myappend(Outputlist, CFFactor(swapvar(inter.getItem().factor(),k,n), inter.getItem().exp()));
280          }
281          return Outputlist;
282        }
283      }
284    }
285    // really is Pth power
286    DEBOUTLN(CERR, "f is a p'th root: ", f);
287    CFMap m;
288    g = compress(f,m);
289    if (mipo.isZero())
290      f = m(PthRoot(g));
291    else
292      f = m(PthRoot(g,mipo));
293    DEBOUTLN(CERR, "  that is       : ", f);
294    // now : Outputlist union ( SqrFreed(f) )^getCharacteristic()
295    Outputlist=myUnion(Powerup(SqrFreeMV(f),getCharacteristic()),Outputlist);
296    DEBDECLEVEL(CERR, "SqrFreed");
297    return Outputlist ;
298  }
299  g = f.deriv();
300  DEBOUTLN(CERR, "calculating gcd of ", f);
301  DEBOUTLN(CERR, "               and ", g);
302  h = gcd(f,pp(g));  h /= lc(h);
303  DEBOUTLN(CERR,"gcd(f,g)= ",h);
304  if ( (h.isOne()) || ( h==f) || ((-h).isOne()) || getNumVars(h)==0 )
305  { // no common factor
306    Outputlist= myappend(Outputlist,CFFactor(f,1)) ;
307    DEBOUTLN(CERR, "Outputlist= ", Outputlist);
308    DEBDECLEVEL(CERR, "SqrFreed");
309    return Outputlist ;
310  }
311  else
312  { // we can split into two nontrivial pieces
313    f /= h; // Now we have split the poly into f and h
314    g = lc(f);
315    if ( !g.isOne() && getNumVars(g) == 0 )
316    {
317       Outputlist= myappend(Outputlist,CFFactor(g,1)) ;
318       f /= g;
319    }
320    DEBOUTLN(CERR, "Split into f= ", f);
321    DEBOUTLN(CERR, "       and h= ", h);
322    // For char > 0 the polys f and h can have Pth roots; so we need a test
323    // Test is disabled because timing is the same
324
325//    if ( SqrFreeTest(f,0) )
326//      Outputlist= myappend(Outputlist,CFFactor(f,1)) ;
327//    else
328    Outputlist=myUnion(Outputlist, SqrFreeMV(f));
329//    if ( SqrFreeTest(h,0) )
330//      Outputlist= myappend(Outputlist,CFFactor(h,1)) ;
331//    else
332    Outputlist=myUnion(Outputlist,SqrFreeMV(h));
333    DEBOUTLN(CERR, "Returning list ", Outputlist);
334    DEBDECLEVEL(CERR, "SqrFreed");
335    return Outputlist ;
336  }
337  factoryError("libfac: ERROR: SqrFreed: we should never fall trough here!");
338  DEBDECLEVEL(CERR, "SqrFreed");
339  return Outputlist; // for safety purpose
340}
341
342///////////////////////////////////////////////////////////////
343// The user front-end for the SqrFreed routine.              //
344// Input can have a constant as content                      //
345///////////////////////////////////////////////////////////////
346CFFList
347SqrFreeMV( const CanonicalForm & r , const CanonicalForm & mipo )
348{
349  CanonicalForm g=icontent(r), f = r;
350  CFFList Outputlist, Outputlist2,tmpOutputlist;
351
352  DEBINCLEVEL(CERR, "SqrFreeMV");
353  DEBOUTLN(CERR,"Called with f= ", f);
354
355  // Take care of stupid users giving us constants
356  if ( getNumVars(f) == 0 )
357  { // a constant ; Exp==1 even if f==0
358      Outputlist= myappend(Outputlist,CFFactor(f,1));
359  }
360  else
361  {
362      // Now we are sure: we have a nonconstant polynomial
363      g = lc(f);
364      while ( getNumVars(g) != 0 ) g=content(g);
365      if ( ! g.isOne() ) Outputlist= myappend(Outputlist,CFFactor(g,1)) ;
366      f /= g;
367      if ( getNumVars(f) != 0 ) // a real polynomial
368      {
369        if (!mipo.isZero())
370        {
371          #if 0
372          Variable alpha=rootOf(mipo);
373          CanonicalForm ff=replacevar(f,mipo.mvar(),alpha);
374          tmpOutputlist=SqrFreeMV(ff,0);
375          ff=replacevar(f,alpha,mipo.mvar());
376          for ( CFFListIterator i=tmpOutputlist; i.hasItem(); i++ )
377          {
378            ff=i.getItem().factor();
379            ff /= ff.Lc();
380            ff=replacevar(ff,alpha,mipo.mvar());
381            Outputlist=myappend(Outputlist,CFFactor(ff,1));
382          }
383          #else
384          Outputlist=myUnion(SqrFreed(f,mipo),Outputlist) ;
385          #endif
386        }
387        else
388          Outputlist=myUnion(SqrFreed(f),Outputlist) ;
389      }
390  }
391  DEBOUTLN(CERR,"Outputlist = ", Outputlist);
392  for ( CFFListIterator i=Outputlist; i.hasItem(); i++ )
393    if ( getNumVars(i.getItem().factor()) > 0 )
394      Outputlist2.append(i.getItem());
395
396  DEBOUTLN(CERR,"Outputlist2 = ", Outputlist2);
397  DEBDECLEVEL(CERR, "SqrFreeMV");
398  return Outputlist2 ;
399}
400
401CFFList SqrFree(const CanonicalForm & r )
402{
403  CFFList outputlist, sqrfreelist = SqrFreeMV(r);
404  CFFListIterator i;
405  CanonicalForm elem;
406  int n=totaldegree(r);
407
408  DEBINCLEVEL(CERR, "SqrFree");
409
410  if ( sqrfreelist.length() < 2 )
411  {
412    DEBDECLEVEL(CERR, "SqrFree");
413    return sqrfreelist;
414  }
415  for ( int j=1; j<=n; j++ )
416  {
417    elem =1;
418    for ( i = sqrfreelist; i.hasItem() ; i++ )
419    {
420      if ( i.getItem().exp() == j ) elem *= i.getItem().factor();
421    }
422    if ( !(elem.isOne()) ) outputlist.append(CFFactor(elem,j));
423  }
424  elem=1;
425  for ( i=outputlist; i.hasItem(); i++ )
426    if ( getNumVars(i.getItem().factor()) > 0 )
427      elem*= power(i.getItem().factor(),i.getItem().exp());
428  elem= r/elem;
429  outputlist.insert(CFFactor(elem,1));
430
431  DEBOUTLN(CERR, "SqrFree returns list:", outputlist);
432  DEBDECLEVEL(CERR, "SqrFree");
433  return outputlist;
434}
Note: See TracBrowser for help on using the repository browser.