source: git/libfac/factor/SqrFree.cc @ 18500b

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