source: git/libfac/factor/SqrFree.cc @ 8444b0

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