source: git/libfac/factor/SqrFree.cc @ 639047e

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