source: git/libfac/factor/SqrFree.cc @ aa7480c

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