source: git/factory/cf_factor.cc @ 341696

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