source: git/factory/cf_factor.cc @ 3fd1ff

spielwiese
Last change on this file since 3fd1ff was 3fd1ff, checked in by Hans Schönemann <hannes@…>, 16 years ago
*hannes: ezgcd_p git-svn-id: file:///usr/local/Singular/svn/trunk@10821 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 18.3 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_factor.cc,v 1.41 2008-07-01 13:45:44 Singular Exp $ */
3
4//{{{ docu
5//
6// cf_factor.cc - factorization and square free algorithms.
7//
8// Used by: fac_multivar.cc, fac_univar.cc, cf_irred.cc
9//
10// Header file: cf_algorithm.h
11//
12//}}}
13
14#include <config.h>
15
16#include "cf_gmp.h"
17
18#include "assert.h"
19
20#include "cf_defs.h"
21#include "canonicalform.h"
22#include "cf_iter.h"
23#include "fac_berlekamp.h"
24#include "fac_cantzass.h"
25#include "fac_univar.h"
26#include "fac_multivar.h"
27#include "fac_sqrfree.h"
28#include "cf_algorithm.h"
29#include "cf_map.h"
30
31#include "int_int.h"
32#ifdef HAVE_NTL
33#include "NTLconvert.h"
34#endif
35
36int getExp(); /* cf_char.cc */
37
38static bool isUnivariateBaseDomain( const CanonicalForm & f )
39{
40    CFIterator i = f;
41    bool ok = i.coeff().inBaseDomain();
42    i++;
43    while ( i.hasTerms() && ( ok = ok && i.coeff().inBaseDomain() ) ) i++;
44    return ok;
45}
46
47void find_exp(const CanonicalForm & f, int * exp_f)
48{
49  if ( ! f.inCoeffDomain() )
50  {
51    int e=f.level();
52    CFIterator i = f;
53    if (e>=0)
54    {
55      if (i.exp() > exp_f[e]) exp_f[e]=i.exp();
56    }
57    for (; i.hasTerms(); i++ )
58    {
59      find_exp(i.coeff(), exp_f);
60    }
61  }
62}
63
64int find_mvar(const CanonicalForm & f)
65{
66  int mv=f.level();
67  int *exp_f=new int[mv+1];
68  int i;
69  for(i=mv;i>0;i--) exp_f[i]=0;
70  find_exp(f,exp_f);
71  for(i=mv;i>0;i--)
72  {
73    if ((exp_f[i]>0) && (exp_f[i]<exp_f[mv]))
74    {
75      mv=i;
76    }
77  }
78  delete[] exp_f;
79  return mv;
80}
81
82#if 1
83//#ifndef NOSTREAMIO
84void out_cf(char *s1,const CanonicalForm &f,char *s2)
85{
86  printf("%s",s1);
87  if (f.isZero()) printf("+0");
88  //else if (! f.inCoeffDomain() )
89  else if (! f.inBaseDomain() )
90  {
91    int l = f.level();
92    for ( CFIterator i = f; i.hasTerms(); i++ )
93    {
94      int e=i.exp();
95      if (i.coeff().isOne())
96      {
97        printf("+");
98        if (e==0) printf("1");
99        else
100        {
101          printf("v(%d)",l);
102          if (e!=1) printf("^%d",e);
103        }
104      }
105      else
106      {
107        out_cf("+(",i.coeff(),")");
108        if (e!=0)
109        {
110          printf("*v(%d)",l);
111          if (e!=1) printf("^%d",e);
112        }
113      }
114    }
115  }
116  else
117  {
118    if ( f.isImm() )
119    {
120      printf("+%d",f.intval());
121    }
122    else
123    {
124    #ifdef NOSTREAMIO
125      #ifdef SINGULAR
126      if (f.inZ())
127      {
128        MP_INT m=gmp_numerator(f);
129        char * str = new char[mpz_sizeinbase( &m, 10 ) + 2];
130        str = mpz_get_str( str, 10, &m );
131        printf("%s",str);
132        delete[] str;
133      }
134      else if (f.inQ())
135      {
136        MP_INT m=gmp_numerator(f);
137        char * str = new char[mpz_sizeinbase( &m, 10 ) + 2];
138        str = mpz_get_str( str, 10, &m );
139        printf("%s/",str);
140        delete[] str;
141        m=gmp_denominator(f);
142        str = new char[mpz_sizeinbase( &m, 10 ) + 2];
143        str = mpz_get_str( str, 10, &m );
144        printf("%s",str);
145        delete[] str;
146      }
147      #else
148      printf("+...");
149      #endif
150    #else
151       std::cout << f;
152    #endif
153    }
154    //if (f.inZ()) printf("(Z)");
155    //else if (f.inQ()) printf("(Q)");
156    //else if (f.inFF()) printf("(FF)");
157    //else if (f.inPP()) printf("(PP)");
158    //else if (f.inGF()) printf("(PP)");
159    //else
160    if (f.inExtension()) printf("E(%d)",f.level());
161  }
162  printf("%s",s2);
163}
164void out_cff(CFFList &L)
165{
166  int n = L.length();
167  CFFListIterator J=L;
168  int j=0;
169  for ( ; J.hasItem(); J++, j++ )
170  {
171    printf("F%d",j);out_cf(":",J.getItem().factor()," ^ ");
172    printf("%d\n", J.getItem().exp());
173  }
174}
175void test_cff(CFFList &L,const CanonicalForm & f)
176{
177  int n = L.length();
178  CFFListIterator J=L;
179  CanonicalForm t=1;
180  int j=0;
181  if (!(L.getFirst().factor().inCoeffDomain()))
182    printf("first entry is not const\n");
183  for ( ; J.hasItem(); J++, j++ )
184  {
185    CanonicalForm tt=J.getItem().factor();
186    if (tt.inCoeffDomain() && (j!=0))
187      printf("other entry is const\n");
188    j=J.getItem().exp();
189    while(j>0) { t*=tt; j--; }
190  }
191  if ((f-t)!=0) { printf("problem:\n");out_cf("factor:",f," has problems\n");}
192}
193//#endif
194#endif
195
196bool isPurePoly_m(const CanonicalForm & f)
197{
198  if (f.inBaseDomain()) return true;
199  if (f.level()<0) return false;
200  for (CFIterator i=f;i.hasTerms();i++)
201  {
202    if (!isPurePoly_m(i.coeff())) return false;
203  }
204  return true;
205}
206bool isPurePoly(const CanonicalForm & f)
207{
208  if (f.level()<=0) return false;
209  for (CFIterator i=f;i.hasTerms();i++)
210  {
211    if (!(i.coeff().inBaseDomain())) return false;
212  }
213  return true;
214}
215
216
217///////////////////////////////////////////////////////////////
218// get_max_degree_Variable returns Variable with             //
219// highest degree. We assume f is *not* a constant!          //
220///////////////////////////////////////////////////////////////
221Variable
222get_max_degree_Variable(const CanonicalForm & f)
223{
224  ASSERT( ( ! f.inCoeffDomain() ), "no constants" );
225  int max=0, maxlevel=0, n=level(f);
226  for ( int i=1; i<=n; i++ )
227  {
228    if (degree(f,Variable(i)) >= max)
229    {
230      max= degree(f,Variable(i)); maxlevel= i;
231    }
232  }
233  return Variable(maxlevel);
234}
235
236///////////////////////////////////////////////////////////////
237// get_Terms: Split the polynomial in the containing terms.  //
238// getTerms: the real work is done here.                     //
239///////////////////////////////////////////////////////////////
240void
241getTerms( const CanonicalForm & f, const CanonicalForm & t, CFList & result )
242{
243  if ( getNumVars(f) == 0 ) result.append(f*t);
244  else{
245    Variable x(level(f));
246    for ( CFIterator i=f; i.hasTerms(); i++ )
247      getTerms( i.coeff(), t*power(x,i.exp()), result);
248  }
249}
250CFList
251get_Terms( const CanonicalForm & f ){
252  CFList result,dummy,dummy2;
253  CFIterator i;
254  CFListIterator j;
255
256  if ( getNumVars(f) == 0 ) result.append(f);
257  else{
258    Variable _x(level(f));
259    for ( i=f; i.hasTerms(); i++ ){
260      getTerms(i.coeff(), 1, dummy);
261      for ( j=dummy; j.hasItem(); j++ )
262        result.append(j.getItem() * power(_x, i.exp()));
263
264      dummy= dummy2; // have to initalize new
265    }
266  }
267  return result;
268}
269
270
271///////////////////////////////////////////////////////////////
272// homogenize homogenizes f with Variable x                  //
273///////////////////////////////////////////////////////////////
274
275CanonicalForm
276homogenize( const CanonicalForm & f, const Variable & x)
277{
278#if 0
279  int maxdeg=totaldegree(f), deg;
280  CFIterator i;
281  CanonicalForm elem, result(0);
282
283  for (i=f; i.hasTerms(); i++)
284  {
285    elem= i.coeff()*power(f.mvar(),i.exp());
286    deg = totaldegree(elem);
287    if ( deg < maxdeg )
288      result += elem * power(x,maxdeg-deg);
289    else
290      result+=elem;
291  }
292  return result;
293#else
294  CFList Newlist, Termlist= get_Terms(f);
295  int maxdeg=totaldegree(f), deg;
296  CFListIterator i;
297  CanonicalForm elem, result(0);
298
299  for (i=Termlist; i.hasItem(); i++)
300  {
301    elem= i.getItem();
302    deg = totaldegree(elem);
303    if ( deg < maxdeg )
304      Newlist.append(elem * power(x,maxdeg-deg));
305    else
306      Newlist.append(elem);
307  }
308  for (i=Newlist; i.hasItem(); i++) // rebuild
309    result += i.getItem();
310
311  return result;
312#endif
313}
314
315CanonicalForm
316homogenize( const CanonicalForm & f, const Variable & x, const Variable & v1, const Variable & v2)
317{
318#if 0
319  int maxdeg=totaldegree(f), deg;
320  CFIterator i;
321  CanonicalForm elem, result(0);
322
323  for (i=f; i.hasTerms(); i++)
324  {
325    elem= i.coeff()*power(f.mvar(),i.exp());
326    deg = totaldegree(elem);
327    if ( deg < maxdeg )
328      result += elem * power(x,maxdeg-deg);
329    else
330      result+=elem;
331  }
332  return result;
333#else
334  CFList Newlist, Termlist= get_Terms(f);
335  int maxdeg=totaldegree(f), deg;
336  CFListIterator i;
337  CanonicalForm elem, result(0);
338
339  for (i=Termlist; i.hasItem(); i++)
340  {
341    elem= i.getItem();
342    deg = totaldegree(elem,v1,v2);
343    if ( deg < maxdeg )
344      Newlist.append(elem * power(x,maxdeg-deg));
345    else
346      Newlist.append(elem);
347  }
348  for (i=Newlist; i.hasItem(); i++) // rebuild
349    result += i.getItem();
350
351  return result;
352#endif
353}
354
355#ifdef SINGULAR
356extern int singular_homog_flag;
357#else
358#define singular_homog_flag 1
359#endif
360int cmpCF( const CFFactor & f, const CFFactor & g )
361{
362  if (f.exp() > g.exp()) return 1;
363  if (f.exp() < g.exp()) return 0;
364  if (f.factor() > g.factor()) return 1;
365  return 0;
366}
367
368CFFList factorize ( const CanonicalForm & f, bool issqrfree )
369{
370  if ( f.inCoeffDomain() )
371        return CFFList( f );
372  int mv=f.level();
373  int org_v=mv;
374  //out_cf("factorize:",f,"==================================\n");
375  if (! f.isUnivariate() )
376  {
377    if ( singular_homog_flag && f.isHomogeneous())
378    {
379      Variable xn = get_max_degree_Variable(f);
380      int d_xn = degree(f,xn);
381      CFMap n;
382      CanonicalForm F = compress(f(1,xn),n);
383      CFFList Intermediatelist;
384      Intermediatelist = factorize(F);
385      CFFList Homoglist;
386      CFFListIterator j;
387      for ( j=Intermediatelist; j.hasItem(); j++ )
388      {
389        Homoglist.append(
390            CFFactor( n(j.getItem().factor()), j.getItem().exp()) );
391      }
392      CFFList Unhomoglist;
393      CanonicalForm unhomogelem;
394      for ( j=Homoglist; j.hasItem(); j++ )
395      {
396        unhomogelem= homogenize(j.getItem().factor(),xn);
397        Unhomoglist.append(CFFactor(unhomogelem,j.getItem().exp()));
398        d_xn -= (degree(unhomogelem,xn)*j.getItem().exp());
399      }
400      if ( d_xn != 0 ) // have to append xn^(d_xn)
401        Unhomoglist.append(CFFactor(CanonicalForm(xn),d_xn));
402      if(isOn(SW_USE_NTL_SORT)) Unhomoglist.sort(cmpCF);
403      return Unhomoglist;
404    }
405    mv=find_mvar(f);
406    if ( getCharacteristic() == 0 )
407    {
408      if (mv!=f.level())
409      {
410        swapvar(f,Variable(mv),f.mvar());
411      }
412    }
413    else
414    {
415      if (mv!=1)
416      {
417        swapvar(f,Variable(mv),Variable(1));
418        org_v=1;
419      }
420    }
421  }
422  CFFList F;
423  if ( getCharacteristic() > 0 )
424  {
425    ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
426    #ifdef HAVE_NTL
427    if (isOn(SW_USE_NTL) && (isPurePoly(f)))
428    {
429      // USE NTL
430      if (getCharacteristic()!=2)
431      {
432        if (fac_NTL_char!=getCharacteristic())
433        {
434          fac_NTL_char=getCharacteristic();
435          #ifndef NTL_ZZ
436          if (fac_NTL_char >NTL_SP_BOUND)
437          {
438            ZZ r;
439            r=getCharacteristic();
440            ZZ_pContext ccc(r);
441            ccc.restore();
442            ZZ_p::init(r);
443          }
444          else
445          #endif
446          {
447            #ifdef NTL_ZZ
448            ZZ r;
449            r=getCharacteristic();
450            ZZ_pContext ccc(r);
451            #else
452            zz_pContext ccc(getCharacteristic());
453            #endif
454            ccc.restore();
455            #ifdef NTL_ZZ
456            ZZ_p::init(r);
457            #else
458            zz_p::init(getCharacteristic());
459            #endif
460          }
461        }
462       #ifndef NTL_ZZ
463       if (fac_NTL_char >NTL_SP_BOUND)
464       {
465          // convert to NTL
466          ZZ_pX f1=convertFacCF2NTLZZpX(f);
467          ZZ_p leadcoeff = LeadCoeff(f1);
468          //make monic
469          f1=f1 / LeadCoeff(f1);
470          // factorize
471          vec_pair_ZZ_pX_long factors;
472          CanZass(factors,f1);
473          // convert back to factory
474          F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
475        }
476        else
477        #endif
478        {
479          // convert to NTL
480          #ifdef NTL_ZZ
481          ZZ_pX f1=convertFacCF2NTLZZpX(f);
482          ZZ_p leadcoeff = LeadCoeff(f1);
483          #else
484          zz_pX f1=convertFacCF2NTLzzpX(f);
485          zz_p leadcoeff = LeadCoeff(f1);
486          #endif
487          //make monic
488          f1=f1 / LeadCoeff(f1);
489          // factorize
490          #ifdef NTL_ZZ
491          vec_pair_ZZ_pX_long factors;
492          #else
493          vec_pair_zz_pX_long factors;
494          #endif
495          CanZass(factors,f1);
496          // convert back to factory
497          #ifdef NTL_ZZ
498          F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
499          #else
500          F=convertNTLvec_pair_zzpX_long2FacCFFList(factors,leadcoeff,f.mvar());
501          #endif
502        }
503        //test_cff(F,f);
504      }
505      else
506      {
507        // Specialcase characteristic==2
508        if (fac_NTL_char!=2)
509        {
510          fac_NTL_char=2;
511          zz_p::init(2);
512        }
513        // convert to NTL using the faster conversion routine for characteristic 2
514        GF2X f1=convertFacCF2NTLGF2X(f);
515        // no make monic necessary in GF2
516        //factorize
517        vec_pair_GF2X_long factors;
518        CanZass(factors,f1);
519
520        // convert back to factory again using the faster conversion routine for vectors over GF2X
521        F=convertNTLvec_pair_GF2X_long2FacCFFList(factors,LeadCoeff(f1),f.mvar());
522      }
523    }
524    else
525    #endif
526    {  // Use Factory without NTL
527      if ( isOn( SW_BERLEKAMP ) )
528         F=FpFactorizeUnivariateB( f, issqrfree );
529      else
530        F=FpFactorizeUnivariateCZ( f, issqrfree, 0, Variable(), Variable() );
531    }
532  }
533  else
534  {
535    bool on_rational = isOn(SW_RATIONAL);
536    On(SW_RATIONAL);
537    CanonicalForm cd = bCommonDen( f );
538    CanonicalForm fz = f * cd;
539    Off(SW_RATIONAL);
540    if ( f.isUnivariate() )
541    {
542      #ifdef HAVE_NTL
543      if ((isOn(SW_USE_NTL)) && (isPurePoly(f)))
544      {
545        //USE NTL
546        CanonicalForm ic=icontent(fz);
547        fz/=ic;
548        ZZ c;
549        vec_pair_ZZX_long factors;
550        //factorize the converted polynomial
551        factor(c,factors,convertFacCF2NTLZZX(fz));
552
553        //convert the result back to Factory
554        F=convertNTLvec_pair_ZZX_long2FacCFFList(factors,c,fz.mvar());
555        if ( ! ic.isOne() )
556        {
557          if ( F.getFirst().factor().inCoeffDomain() )
558          {
559            CFFactor new_first( F.getFirst().factor() * ic );
560            F.removeFirst();
561            F.insert( new_first );
562          }
563          else
564            F.insert( CFFactor( ic ) );
565        }
566        else
567        {
568          if ( !F.getFirst().factor().inCoeffDomain() )
569          {
570            CFFactor new_first( 1 );
571            F.insert( new_first );
572          }
573        }
574        //if ( F.getFirst().factor().isOne() )
575        //{
576        //  F.removeFirst();
577        //}
578        //printf("NTL:\n");out_cff(F);
579        //F=ZFactorizeUnivariate( fz, issqrfree );
580        //printf("fac.:\n");out_cff(F);
581      }
582      else
583      #endif
584      {
585        //Use Factory without NTL
586        F = ZFactorizeUnivariate( fz, issqrfree );
587      }
588    }
589    else
590    {
591      F = ZFactorizeMultivariate( fz, issqrfree );
592    }
593
594    if ( on_rational )
595      On(SW_RATIONAL);
596    if ( ! cd.isOne() )
597    {
598      if ( F.getFirst().factor().inCoeffDomain() )
599      {
600        CFFactor new_first( F.getFirst().factor() / cd );
601        F.removeFirst();
602        F.insert( new_first );
603      }
604      else
605      {
606        F.insert( CFFactor( 1/cd ) );
607      }
608    }
609  }
610
611  if ((mv!=org_v) && (! f.isUnivariate() ))
612  {
613    CFFListIterator J=F;
614    for ( ; J.hasItem(); J++)
615    {
616      swapvar(J.getItem().factor(),Variable(mv),Variable(org_v));
617    }
618    swapvar(f,Variable(mv),Variable(org_v));
619  }
620  //out_cff(F);
621  if(isOn(SW_USE_NTL_SORT)) F.sort(cmpCF);
622  return F;
623}
624
625#ifdef HAVE_NTL
626CanonicalForm fntl ( const CanonicalForm & f, int j )
627{
628  ZZX f1=convertFacCF2NTLZZX(f);
629  return convertZZ2CF(coeff(f1,j));
630}
631#endif
632
633CFFList factorize ( const CanonicalForm & f, const Variable & alpha )
634{
635  //out_cf("factorize:",f,"==================================\n");
636  //out_cf("mipo:",getMipo(alpha),"\n");
637  CFFList F;
638  ASSERT( alpha.level() < 0, "not an algebraic extension" );
639  ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
640  ASSERT( getCharacteristic() > 0, "char 0 factorization not implemented" );
641  #ifdef HAVE_NTL
642  if  (isOn(SW_USE_NTL))
643  {
644    //USE NTL
645    if (getCharacteristic()!=2)
646    {
647      // First all cases with characteristic !=2
648      // set remainder
649      if (fac_NTL_char!=getCharacteristic())
650      {
651        fac_NTL_char=getCharacteristic();
652        #ifdef NTL_ZZ
653        ZZ r;
654        r=getCharacteristic();
655        ZZ_pContext ccc(r);
656        #else
657        zz_pContext ccc(getCharacteristic());
658        #endif
659        ccc.restore();
660        #ifdef NTL_ZZ
661        ZZ_p::init(r);
662        #else
663        zz_p::init(getCharacteristic());
664        #endif
665      }
666
667      // set minimal polynomial in NTL
668      #ifdef NTL_ZZ
669      ZZ_pX minPo=convertFacCF2NTLZZpX(getMipo(alpha));
670      ZZ_pEContext c(minPo);
671      #else
672      zz_pX minPo=convertFacCF2NTLzzpX(getMipo(alpha));
673      zz_pEContext c(minPo);
674      #endif
675
676      c.restore();
677
678      // convert to NTL
679      #ifdef NTL_ZZ
680      ZZ_pEX f1=convertFacCF2NTLZZ_pEX(f,minPo);
681      ZZ_pE leadcoeff= LeadCoeff(f1);
682      #else
683      zz_pEX f1=convertFacCF2NTLzz_pEX(f,minPo);
684      zz_pE leadcoeff= LeadCoeff(f1);
685      #endif
686
687      //make monic
688      f1=f1 / leadcoeff;
689
690      // factorize using NTL
691      #ifdef NTL_ZZ
692      vec_pair_ZZ_pEX_long factors;
693      #else
694      vec_pair_zz_pEX_long factors;
695      #endif
696      CanZass(factors,f1);
697
698      // return converted result
699      F=convertNTLvec_pair_zzpEX_long2FacCFFList(factors,leadcoeff,f.mvar(),alpha);
700    }
701    else
702    {
703      // special case : GF2
704
705      // remainder is two ==> nothing to do
706      // set remainder
707      ZZ r;
708      r=getCharacteristic();
709      ZZ_pContext ccc(r);
710      ccc.restore();
711
712      // set minimal polynomial in NTL using the optimized conversion routines for characteristic 2
713      GF2X minPo=convertFacCF2NTLGF2X(getMipo(alpha,f.mvar()));
714      GF2EContext c(minPo);
715      c.restore();
716
717      // convert to NTL again using the faster conversion routines
718      GF2EX f1;
719      if (isPurePoly(f))
720      {
721        GF2X f_tmp=convertFacCF2NTLGF2X(f);
722        f1=to_GF2EX(f_tmp);
723      }
724      else
725      {
726        f1=convertFacCF2NTLGF2EX(f,minPo);
727      }
728
729      // make monic (in Z/2(a))
730      GF2E f1_coef=LeadCoeff(f1);
731      MakeMonic(f1);
732
733      // factorize using NTL
734      vec_pair_GF2EX_long factors;
735      CanZass(factors,f1);
736
737      // return converted result
738      F=convertNTLvec_pair_GF2EX_long2FacCFFList(factors,f1_coef,f.mvar(),alpha);
739    }
740
741  }
742  else
743  #endif
744  {
745    F=FpFactorizeUnivariateCZ( f, false, 1, alpha, Variable() );
746  }
747  return F;
748}
749
750CFFList sqrFree ( const CanonicalForm & f )
751{
752//    ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
753    CFFList result;
754
755    if ( getCharacteristic() == 0 )
756        result = sqrFreeZ( f );
757    else
758        result = sqrFreeFp( f );
759
760    //return ( sort ? sortCFFList( result ) : result );
761    return result;
762}
763
764bool isSqrFree ( const CanonicalForm & f )
765{
766//    ASSERT( f.isUnivariate(), "multivariate factorization not implemented" );
767    if ( getCharacteristic() == 0 )
768        return isSqrFreeZ( f );
769    else
770        return isSqrFreeFp( f );
771}
772
Note: See TracBrowser for help on using the repository browser.