source: git/libfac/factor/SqrFree.cc @ 14b1e65

spielwiese
Last change on this file since 14b1e65 was 14b1e65, checked in by Hans Schönemann <hannes@…>, 22 years ago
*hannes: Dan's NOSTREAMIO changes git-svn-id: file:///usr/local/Singular/svn/trunk@5572 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 11.6 KB
Line 
1/* Copyright 1996 Michael Messollen. All rights reserved. */
2///////////////////////////////////////////////////////////////////////////////
3// emacs edit mode for this file is -*- C++ -*-
4static char * rcsid = "$Id: SqrFree.cc,v 1.5 2001-08-08 11:59:13 Singular Exp $";
5static char * errmsg = "\nYou found a bug!\nPlease inform (Michael Messollen) michael@math.uni-sb.de .\n Please include above information and your input (the ideal/polynomial and characteristic) in your bug-report.\nThank you.";
6///////////////////////////////////////////////////////////////////////////////
7// FACTORY - Includes
8#include<factory.h>
9#ifndef NOSTREAMIO
10#include <iostream.h>
11#endif
12// Factor - Includes
13#include "tmpl_inst.h"
14#include "helpstuff.h"
15// some CC's need this:
16#include "SqrFree.h"
17
18#ifdef SINGULAR
19#  define HAVE_SINGULAR
20   extern "C" { void WerrorS(char *); }
21#endif
22
23#ifdef SQRFREEDEBUG
24# define DEBUGOUTPUT
25#else
26# undef DEBUGOUTPUT
27#endif
28
29#include "debug.h"
30#include "timing.h"
31TIMING_DEFINE_PRINT(squarefree_time);
32TIMING_DEFINE_PRINT(gcd_time);
33
34
35static inline CFFactor
36Powerup( const CFFactor & F , int exp=1){
37  return CFFactor(F.factor(), exp*F.exp()) ;
38}
39
40static CFFList
41Powerup( const CFFList & Inputlist , int exp=1 ){
42  CFFList Outputlist;
43
44  for ( CFFListIterator i=Inputlist; i.hasItem(); i++ )
45    Outputlist.append(Powerup(i.getItem(), exp));
46  return Outputlist ;
47}
48
49///////////////////////////////////////////////////////////////
50// Compute the Pth root of a polynomial in characteristic p  //
51// f must be a polynomial which we can take the Pth root of. //
52// Domain is q=p^m , f a uni/multivariate polynomial         //
53///////////////////////////////////////////////////////////////
54static CanonicalForm
55PthRoot( const CanonicalForm & f ){
56  CanonicalForm RES, R = f;
57  int n= max(level(R),getNumVars(R)), p= getCharacteristic();
58
59  if (n==0){ // constant
60    if (R.inExtension()) // not in prime field; f over |F(q=p^k)
61      R = power(R,Powerup(p,getGFDegree() - 1)) ;
62    // if f in prime field, do nothing
63    return R;
64  }
65  // we assume R is a Pth power here
66  RES = R.genZero();
67  Variable x(n);
68  for (int i=0; i<= (int) (degree(R,level(R))/p) ; i++)
69    RES += PthRoot( R[i*p] ) * power(x,i);
70  return RES;
71}
72
73///////////////////////////////////////////////////////////////
74// A uni/multivariate SqrFreeTest routine.                   //
75// Cheaper to run if all you want is a test.                 //
76// Works for charcteristic 0 and q=p^m                       //
77// Returns 1 if poly r is SqrFree, 0 if SqrFree will do some //
78// kind of factorization.                                    //
79// Would be much more effcient iff we had *good*             //
80//  uni/multivariate gcd's and/or gcdtest's                  //
81///////////////////////////////////////////////////////////////
82int
83SqrFreeTest( const CanonicalForm & r, int opt){
84  CanonicalForm f=r, g;
85  int n=level(f);
86
87  if (getNumVars(f)==0) return 1 ; // a constant is SqrFree
88  if ( f.isUnivariate() ) {
89    g= f.deriv();
90    if ( getCharacteristic() > 0 && g.isZero() ) return 0 ;
91    // Next: it would be best to have a *univariate* gcd-test which returns
92    // 0 iff gcdtest(f,g) == 1 or a constant ( for real Polynomials )
93    g = mygcd(f,g);
94    if ( g.isOne() || (-g).isOne() ) return 1;
95    else
96      if ( getNumVars(g) == 0 ) return 1;// <- totaldegree!!!
97      else return 0 ;
98  }
99  else { // multivariate case
100    for ( int k=1; k<=n; k++ ) {
101      g = swapvar(f,k,n); g = content(g);
102      // g = 1 || -1 : sqr-free, g poly : not sqr-free, g number : opt helps
103      if ( ! (g.isOne() || (-g).isOne() || getNumVars(g)==0 ) ) {
104        if ( opt==0 ) return 0;
105        else {
106          if ( SqrFreeTest(g,1) == 0 ) return 0;
107          g = swapvar(g,k,n);
108          f /=g ;
109        }
110      }
111    }
112    // Now f is primitive
113    n = level(f); // maybe less indeterminants
114    //    if ( totaldegree(f) <= 1 ) return 1;
115
116    // Let`s look if it is a Pth root
117    if ( getCharacteristic() > 0 )
118      for (int k=1; k<=n; k++ ) {
119        g=swapvar(f,k,n); g=g.deriv();
120        if ( ! g.isZero() ) break ;
121        else if ( k==n) return 0 ; // really is Pth root
122      }
123    g = f.deriv() ;
124    // Next: it would be best to have a *multivariate* gcd-test which returns
125    // 0 iff gcdtest(f,g) == 1 or a constant ( for real Polynomials )
126    g= mygcd(f,g);
127    if ( g.isOne() || (-g).isOne() || (g==f) || (getNumVars(g)==0) ) return 1 ;
128    else return 0 ;
129  }
130#ifdef HAVE_SINGULAR
131  WerrorS("libfac: ERROR: SqrFreeTest: we should never fall trough here!");
132#else
133  cerr << "\nlibfac: ERROR: SqrFreeTest: we should never fall trough here!\n"
134       << rcsid << errmsg << endl;
135#endif
136  return 0;
137}
138
139///////////////////////////////////////////////////////////////
140// A uni/multivariate SqrFree routine.Works for polynomials  //
141// which don\'t have a constant as the content.              //
142// Works for charcteristic 0 and q=p^m                       //
143// returns a list of polys each of sqrfree, but gcd(f_i,f_j) //
144// needs not to be 1 !!!!!                                   //
145///////////////////////////////////////////////////////////////
146static CFFList
147SqrFreed( const CanonicalForm & r ){
148  CanonicalForm h, g, f = r;
149  CFFList Outputlist;
150  int n = level(f);
151
152  DEBINCLEVEL(cout, "SqrFreed");
153  DEBOUTLN(cout, "Called with r= ", r);
154  if (getNumVars(f)==0 ) { // just a constant; return it
155    Outputlist= myappend(Outputlist,CFFactor(f,1));
156    return Outputlist ;
157  }
158
159// We look if we do have a content; if so, SqrFreed the content
160// and continue computations with pp(f)
161  for (int k=1; k<=n; k++) {
162    g = swapvar(f,k,n); g = content(g);
163    if ( ! (g.isOne() || (-g).isOne() || degree(g)==0 )) {
164      g = swapvar(g,k,n);
165      DEBOUTLN(cout, "We have a content: ", g);
166      Outputlist = myUnion(InternalSqrFree(g),Outputlist); // should we add a
167                                                // SqrFreeTest(g) first ?
168      DEBOUTLN(cout, "Outputlist is now: ", Outputlist);
169      f /=g;
170      DEBOUTLN(cout, "f is now: ", f);
171    }
172  }
173
174// Now f is primitive; Let`s look if f is univariate
175  if ( f.isUnivariate() ) {
176    DEBOUTLN(cout, "f is univariate: ", f);
177    g = content(g);
178    if ( ! (g.isOne() || (-g).isOne() ) ){
179      Outputlist= myappend(Outputlist,CFFactor(g,1)) ;
180      f /= g;
181    }
182    Outputlist = Union(sqrFree(f),Outputlist) ;
183    DEBOUTLN(cout, "Outputlist after univ. sqrFree(f) = ", Outputlist);
184    DEBDECLEVEL(cout, "SqrFreed");
185    return Outputlist ;
186  }
187
188// Linear?
189  if ( totaldegree(f) <= 1 ) {
190    Outputlist= myappend(Outputlist,CFFactor(f,1)) ;
191    DEBDECLEVEL(cout, "SqrFreed");
192    return Outputlist ;
193  }
194
195// is it Pth root?
196  n=level(f); // maybe less indeterminants
197  g= f.deriv();
198  if ( getCharacteristic() > 0 && g.isZero() ){  // Pth roots only apply to char > 0
199    for (int k=1; k<=n; k++) {
200      g=swapvar(f,k,n) ; g = g.deriv();
201      if ( ! g.isZero() ){ // can`t be Pth root
202        CFFList Outputlist2= SqrFreed(swapvar(f,k,n));
203        for (CFFListIterator inter=Outputlist2; inter.hasItem(); inter++){
204          Outputlist= myappend(Outputlist, CFFactor(swapvar(inter.getItem().factor(),k,n), inter.getItem().exp()));
205        }
206        return Outputlist;
207      }
208      else
209        if ( k==n ) { // really is Pth power
210          DEBOUTLN(cout, "f is a p'th root: ", f);
211          CFMap m;
212          g = compress(f,m);
213          f = m(PthRoot(g));
214          DEBOUTLN(cout, "  that is       : ", f);
215          // now : Outputlist union ( SqrFreed(f) )^getCharacteristic()
216          Outputlist=myUnion(Powerup(InternalSqrFree(f),getCharacteristic()),Outputlist);
217          DEBDECLEVEL(cout, "SqrFreed");
218          return Outputlist ;
219        }
220    }
221  }
222  g = f.deriv();
223  DEBOUTLN(cout, "calculating mygcd of ", f);
224  DEBOUTLN(cout, "               and ", g);
225  h = mygcd(f,pp(g));  h /= lc(h);
226  DEBOUTLN(cout,"mygcd(f,g)= ",h);
227  if ( (h.isOne()) || ( h==f) || ((-h).isOne()) || getNumVars(h)==0 ) { // no common factor
228    Outputlist= myappend(Outputlist,CFFactor(f,1)) ;
229    DEBOUTLN(cout, "Outputlist= ", Outputlist);
230    DEBDECLEVEL(cout, "SqrFreed");
231    return Outputlist ;
232  }
233  else { // we can split into two nontrivial pieces
234    f /= h; // Now we have split the poly into f and h
235    g = lc(f);
236    if ( g != f.genOne() && getNumVars(g) == 0 ){
237       Outputlist= myappend(Outputlist,CFFactor(g,1)) ;
238       f /= g;
239    }
240    DEBOUTLN(cout, "Split into f= ", f);
241    DEBOUTLN(cout, "       and h= ", h);
242    // For char > 0 the polys f and h can have Pth roots; so we need a test
243    // Test is disabled because timing is the same
244//    if ( SqrFreeTest(f,0) )
245//      Outputlist= myappend(Outputlist,CFFactor(f,1)) ;
246//    else
247    Outputlist=myUnion(Outputlist, InternalSqrFree(f));
248//    if ( SqrFreeTest(h,0) )
249//      Outputlist= myappend(Outputlist,CFFactor(h,1)) ;
250//    else
251    Outputlist=myUnion(Outputlist,InternalSqrFree(h));
252    DEBOUTLN(cout, "Returning list ", Outputlist);
253    DEBDECLEVEL(cout, "SqrFreed");
254    return Outputlist ;
255  }
256#ifdef HAVE_SINGULAR
257  WerrorS("libfac: ERROR: SqrFreed: we should never fall trough here!");
258#else
259  cerr << "\nlibfac: ERROR: SqrFreed: we should never fall trough here!\n"
260       << rcsid << errmsg << endl;
261#endif
262  DEBDECLEVEL(cout, "SqrFreed");
263  return Outputlist; // for safety purpose
264}
265
266///////////////////////////////////////////////////////////////
267// The user front-end for the SqrFreed routine.              //
268// Input can have a constant as content                      //
269///////////////////////////////////////////////////////////////
270CFFList
271InternalSqrFree( const CanonicalForm & r ){
272  CanonicalForm g=icontent(r), f = r;
273  CFFList Outputlist, Outputlist2;
274
275  DEBINCLEVEL(cout, "InternalSqrFree");
276  DEBOUTMSG(cout, rcsid);
277  DEBOUTLN(cout,"Called with f= ", f);
278
279  // Take care of stupid users giving us constants
280  if ( getNumVars(f) == 0 ) { // a constant ; Exp==1 even if f==0
281      Outputlist= myappend(Outputlist,CFFactor(f,1));
282  }
283  else{
284      // Now we are sure: we have a nonconstant polynomial
285      g = lc(f);
286      while ( getNumVars(g) != 0 ) g=content(g);
287      if ( ! g.isOne() ) Outputlist= myappend(Outputlist,CFFactor(g,1)) ;
288      f /= g;
289      if ( getNumVars(f) != 0 ) // a real polynomial
290        Outputlist=myUnion(SqrFreed(f),Outputlist) ;
291  }
292  DEBOUTLN(cout,"Outputlist = ", Outputlist);
293  for ( CFFListIterator i=Outputlist; i.hasItem(); i++ )
294    if ( getNumVars(i.getItem().factor()) > 0 )
295      Outputlist2.append(i.getItem());
296
297  DEBOUTLN(cout,"Outputlist2 = ", Outputlist2);
298  DEBDECLEVEL(cout, "InternalSqrFree");
299  return Outputlist2 ;
300}
301
302CFFList
303SqrFree(const CanonicalForm & r ){
304  CFFList outputlist, sqrfreelist = InternalSqrFree(r);
305  CFFListIterator i;
306  CanonicalForm elem;
307  int n=totaldegree(r);
308
309  DEBINCLEVEL(cout, "SqrFree");
310
311  if ( sqrfreelist.length() < 2 ){
312    DEBDECLEVEL(cout, "SqrFree");
313    return sqrfreelist;
314  }
315  for ( int j=1; j<=n; j++ ){
316    elem =1;
317    for ( i = sqrfreelist; i.hasItem() ; i++ ){
318      if ( i.getItem().exp() == j ) elem *= i.getItem().factor();
319    }
320    if ( elem != 1 ) outputlist.append(CFFactor(elem,j));
321  }
322  elem=1;
323  for ( i=outputlist; i.hasItem(); i++ )
324    if ( getNumVars(i.getItem().factor()) > 0 )
325      elem*= power(i.getItem().factor(),i.getItem().exp());
326  elem= r/elem;
327  outputlist.insert(CFFactor(elem,1));
328
329  DEBOUTLN(cout, "SqrFree returns list:", outputlist);
330  DEBDECLEVEL(cout, "SqrFree");
331  return outputlist;
332}
333
334/*
335$Log: not supported by cvs2svn $
336Revision 1.4  1997/11/18 16:39:06  Singular
337* hannes: moved WerrorS from C++ to C
338     (Factor.cc MVMultiHensel.cc SqrFree.cc Truefactor.cc)
339
340Revision 1.3  1997/09/12 07:19:50  Singular
341* hannes/michael: libfac-0.3.0
342
343Revision 1.4  1997/04/25 22:19:46  michael
344changed cerr and cout messages for use with Singular
345Version for libfac-0.2.1
346
347*/
Note: See TracBrowser for help on using the repository browser.