source: git/libfac/factor/SqrFree.cc @ 052e4f

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