source: git/factory/cf_factor.cc @ a81073

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