source: git/factory/libfac/factor/Factor.cc @ 21263cf

jengelh-datetimespielwiese
Last change on this file since 21263cf was 21263cf, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: replaced hack by setMipo
  • Property mode set to 100644
File size: 39.5 KB
Line 
1///////////////////////////////////////////////////////////////////////////////
2static const char * errmsg = "\nYou found a bug!\nPlease inform singular@mathematik.uni-kl.de\nPlease include above information and your input (the ideal/polynomial and characteristic) in your bug-report.\nThank you.";
3///////////////////////////////////////////////////////////////////////////////
4// FACTORY - Includes
5#include <factory.h>
6#ifndef NOSTREAMIO
7#ifdef HAVE_IOSTREAM
8#include <iostream>
9#define CERR std::cerr
10#define CIN std::cin
11#elif defined(HAVE_IOSTREAM_H)
12#include <iostream.h>
13#define CERR cerr
14#define CIN cin
15#endif
16#endif
17// Factor - Includes
18#include "tmpl_inst.h"
19#include "SqrFree.h"
20#include "helpstuff.h"
21#include "MVMultiHensel.h"
22#include "Truefactor.h"
23#include "homogfactor.h"
24#include "interrupt.h"
25// some CC's need this:
26#include "Factor.h"
27
28#include "alg_factor.h"
29void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
30void out_cff(CFFList &L);
31
32
33#ifdef FACTORDEBUG
34#  define DEBUGOUTPUT
35#else
36#  undef DEBUGOUTPUT
37#endif
38
39#include <libfac/factor/debug.h>
40#include "timing.h"
41TIMING_DEFINE_PRINT(factorize_time)
42TIMING_DEFINE_PRINT(sqrfree_time)
43TIMING_DEFINE_PRINT(discr_time)
44TIMING_DEFINE_PRINT(evaluate_time)
45TIMING_DEFINE_PRINT(hensel_time)
46TIMING_DEFINE_PRINT(truefactor_time)
47
48/*
49* a wrapper for factorize over algebraic extensions:
50* does a sanity check and, if nec., a conversion
51* before calling factorize(f,alpha)
52* ( in factorize, alpha.level() must be < 0 )
53*/
54static
55CFFList factorize2 ( const CanonicalForm & f,
56                     const Variable & alpha, const CanonicalForm & mipo )
57{
58  if (alpha.level() <0)
59  {
60    return factorize(f,alpha);
61  }
62  else
63  {
64    //out_cf("f2 - factor:",f,"\n");
65    //out_cf("f2 - ext:",alpha,"\n");
66    //out_cf("f2 - mipo:",mipo,"\n");
67    Variable X=rootOf(mipo);
68    CanonicalForm F=f;
69    F=replacevar(f,alpha,X);
70    //out_cf("call - factor:",F,"\n");
71    //out_cf("call - ext:",X,"\n");
72    //out_cf("call - mipo:",getMipo(X,'A'),"\n");
73    CFFList L=factorize(F,X);
74    CFFListIterator i=L;
75    {
76      CFFList Outputlist;
77      for(;i.hasItem(); i++ )
78      {
79        Outputlist.append(CFFactor(
80        replacevar(i.getItem().factor(),X,alpha),
81        i.getItem().exp()));
82      }
83      //out_cff(Outputlist);
84      return Outputlist;
85    }
86  }
87}
88///////////////////////////////////////////////////////////////
89// Choose a main variable if the user didn`t wish a          //
90// special one. Returns level of main variable.              //
91///////////////////////////////////////////////////////////////
92static int
93choose_main_variable( const CanonicalForm & f, int Mainvar=0){
94  CanonicalForm remlc, newlc;
95  int n= level(f), mainvar= Mainvar;
96
97  if (mainvar != 0) return mainvar ; // We force use of the wished mainvar
98  remlc= LC(f,n); mainvar = n;
99  if ( totaldegree(remlc)==0 ){ remlc=f.genOne() ; }
100  DEBOUTLN(CERR, "remlc= " , remlc);
101  for ( int i=n-1; i>=1; i-- )
102  {
103    newlc= LC(f,i);
104    if ( totaldegree(newlc)==0 ){ newlc=f.genOne() ; }
105    DEBOUTLN(CERR, "newlc= " , newlc);
106    if ( (remlc.isOne()) && (newlc.isOne()) ){ // take care of the degrees
107      if ( degree(f,i) < degree(f,mainvar) ){
108        remlc= newlc;
109        mainvar= i;
110      }
111    }
112    else  if ( (! remlc.isOne() ) && ( newlc.isOne() ) ){
113      remlc= newlc;
114      mainvar= i;
115    }
116  }
117  return mainvar;
118}
119
120///////////////////////////////////////////////////////////////
121// Check if the derivative is nonzero for oldmainvar.        //
122// Returns the level of the choosen main variable.           //
123///////////////////////////////////////////////////////////////
124static int
125necessary_condition( const CanonicalForm & F, int oldmainvar){
126  CanonicalForm g;
127  int n=level(F);
128
129  g= swapvar(F,oldmainvar,n);
130  g= g.deriv();
131  if ( g.isZero() )
132  {
133    for ( int i=n; i>=1; i-- )
134    {
135      g= swapvar(F,i,n);
136      g= g.deriv();
137      if ( ! g.isZero() ) return i;
138    }
139  }
140  return oldmainvar;
141}
142
143///////////////////////////////////////////////////////////////
144// Make F monic. Return monic polynomial.                    //
145///////////////////////////////////////////////////////////////
146static CanonicalForm
147make_monic( const CanonicalForm & F, const CanonicalForm & lt)
148{
149  CanonicalForm intermediatpoly,f;
150  Variable x(level(F));
151
152  if ( degree(lt) == 0 ) f= 1/lt * F ;
153  else
154  {
155    intermediatpoly= power(lt,degree(F)-1);
156    for ( int i=0; i<=degree(F); i++ )
157      if ( ! F[i].isZero())
158        f+= (F[i] * intermediatpoly*power(x,i))/power(lt,i);
159  }
160  return f;
161}
162
163///////////////////////////////////////////////////////////////
164// Decide whether num/denum (num,denum both from the         //
165// FiniteFielddomain)  lies in the RationalDomain.           //
166// If false, return num/denum else return the zero poly from //
167// the FiniteFielddomain.                                    //
168///////////////////////////////////////////////////////////////
169static CanonicalForm
170is_rational( const CanonicalForm & num, const CanonicalForm & denum ){
171  CanonicalForm a, b;
172  int retvalue;
173
174  retvalue= mydivremt(num,denum,a,b);
175  if ( retvalue && b == num.genZero() ) // num/denum from FFdomain
176    return a;
177  else // num/denum is rational
178    return num.genZero();
179}
180
181///////////////////////////////////////////////////////////////
182// lt_is_product returns 1 iff lt is a product, 0 iff lt is  //
183// a sum.                                                    //
184///////////////////////////////////////////////////////////////
185static int
186lt_is_product( const CanonicalForm & lt ){
187  CFList result;
188
189  result= get_Terms(lt);
190  if ( result.length() > 1 ) return 0;
191  else return 1;
192}
193
194///////////////////////////////////////////////////////////////
195// Reverse the make_monic transformation.                    //
196// Return the list of factors.                               //
197///////////////////////////////////////////////////////////////
198static CFFList
199not_monic( const CFFList & TheList, const CanonicalForm & ltt, const CanonicalForm & F, int levelF)
200{
201  CFFList Returnlist,IntermediateList;
202  CFFListIterator i;
203  CanonicalForm intermediate,lt= ltt,savelc;
204  CanonicalForm numerator,denumerator,test,a,b;
205  Variable x(level(F));
206  int test1;
207
208  if ( lt.isOne() ) return TheList; // the poly was already monic
209  if ( TheList.length() <= 1 ) // only one factor to substitute back
210  {
211    if ( totaldegree(lt) == 0 ) // lt is type numeric
212      Returnlist.append( CFFactor(lt*TheList.getFirst().factor(),
213                                  TheList.getFirst().exp()) );
214    else
215    {
216      intermediate = F(x*lt, levelF)/power(lt,degree(F,levelF)-1);
217      Returnlist.append(CFFactor(intermediate,TheList.getFirst().exp()));
218    }
219  }
220  else // more then one factor
221  {
222    IntermediateList= TheList;
223    if ( totaldegree(lt) == 0 ){ // lt is type numeric;(SqrFree-use, see above)
224      Returnlist.append( CFFactor(lt*IntermediateList.getFirst().factor()
225                                  , IntermediateList.getFirst().exp()) );
226      IntermediateList.removeFirst();
227      Returnlist= Union(Returnlist,IntermediateList);
228    }
229    else // lt is a) a product or b) a sum of terms
230    {
231      if ( lt_is_product(lt) ) // case a)
232      {
233        DEBOUTLN(CERR, "lt_is_product: ", lt);
234        savelc= content(lt) ; // can we simplify to savelc= lc(lt); ?
235        while ( getNumVars(savelc) != 0 )
236          savelc= content(savelc);
237        for ( i=TheList; i.hasItem();i++ )
238        {
239          numerator= i.getItem().factor();
240          numerator= numerator(x*lt,levelF); // x <- x*lt
241          denumerator= power(lt,degree(F,levelF)-1); // == lt^(1-degree(F,x)
242          while (numerator.genZero() == is_rational(numerator, denumerator))
243            numerator*= lt;
244          intermediate= is_rational(numerator,denumerator);
245
246          Returnlist.append( CFFactor(lc(content(intermediate))*intermediate/content(intermediate), i.getItem().exp() ) );
247        }
248        // Now we add a test. If product(factors)/F is a multiple of
249        // savelc, we have to add 1/multiplicity to the factors
250        IntermediateList= Returnlist;
251        intermediate= 1;
252        for ( CFFListIterator j=IntermediateList; j.hasItem(); j++)
253          intermediate*= j.getItem().factor();
254        test1= mydivremt( intermediate,F,a,b);
255        if ( test1 && b == intermediate.genZero() ) // Yupp!
256        {
257          IntermediateList.append(CFFactor(1/a,1));
258          Returnlist= IntermediateList;
259        }
260        else { Returnlist= IntermediateList; }
261      }
262      else // case b)
263      {
264        DEBOUTLN(CERR, "lt_is_sum: ", lt);
265        CanonicalForm save_denumerator= 1;
266        for ( i=TheList; i.hasItem(); i++ )
267        {
268          numerator= i.getItem().factor();
269          numerator= numerator(x*lt,levelF); // x <- x*lt
270          denumerator= power(lt,degree(numerator,levelF)); // == lt^(-degree(numerator,x)
271          test= content(numerator,x);
272          test1= mydivremt(denumerator,test,a,b);
273          if ( test1 && b == numerator.genZero() ) // Yupp!
274          {
275            save_denumerator*= a;
276            Returnlist.append(CFFactor(numerator/test ,1));
277          }
278          else
279          {
280            factoryError("libfac: ERROR: not_monic1: case lt is a sum.");
281          }
282        }
283        // Now we add a test if we did the right thing:
284        // save_denumerator should be a multiple of the leading coeff
285        test1= mydivremt(save_denumerator,lt,a,b);
286        if ( test1 && b == save_denumerator.genZero() ) // Yupp!
287          // We have to multiply one of the factors with
288          // the multiplicity of the save_denumerator <-> lc
289          // the following will do what we want
290          Returnlist= myUnion( CFFList(CFFactor(1/a,1)),Returnlist) ;
291        else
292        {
293          factoryError("libfac: ERROR: not_monic2: case lt is a sum.");
294        }
295      }
296    }
297  }
298  DEBOUTLN(CERR,"Returnlist: ", Returnlist);
299  return Returnlist;
300}
301
302///////////////////////////////////////////////////////////////
303// Substitute the (Variable,Value)-Pair(s) from Substitution-//
304// list into the polynomial F. Returns the resulting poly.   //
305///////////////////////////////////////////////////////////////
306static CanonicalForm
307substitutePoly( const CanonicalForm & F, const SFormList & Substitutionlist){
308  CanonicalForm f= F;
309
310  for ( SFormListIterator i=Substitutionlist; i.hasItem(); i++ )
311    f= f(i.getItem().exp(),level(i.getItem().factor()));
312  return f;
313}
314
315///////////////////////////////////////////////////////////////
316// Find specialization values for the poly F. Returns 0 if   //
317// procedure failed, 1 otherwise. On success Substitutionlist//
318// holds (Variable,Value)-pairs. On failure we only have a   //
319// partitial list.                                           //
320///////////////////////////////////////////////////////////////
321//      *** This is the version with extensions ***          //
322///////////////////////////////////////////////////////////////
323
324///////////////////////////////////////////////////////////////
325// is CF g ok?                                               //
326///////////////////////////////////////////////////////////////
327static int
328various_tests( const CanonicalForm & g, int deg, int vars_left)
329{
330  CFMap m;
331
332  if ( degree(g) == deg ) // degrees match
333    if ( level(compress(g,m)) == (vars_left) ) // exactly one variable less
334      if ( SqrFreeTest(g,1) ) // poly is sqrfree
335        if ( gcd(g,g.deriv()).isOne() ) // Discriminante != 0
336           return 1;
337  return 0;
338}
339
340///////////////////////////////////////////////////////////////
341// specialize one variable over the given field.             //
342///////////////////////////////////////////////////////////////
343// substitutes in poly f of degree deg with former
344// former_nr_of_variables variables the variable nr_of_variable ;
345// this is done in the field of Char getCharacteristic() and
346// Extension given by Extgenerator.
347///////////////////////////////////////////////////////////////
348static int
349specialize_variable( CanonicalForm & f, int deg, SFormList & Substitutionlist, int nr_of_variable,
350                     int former_nr_of_variables, CFGenerator & Extgenerator ){
351  CanonicalForm g;
352  Variable x(nr_of_variable);
353
354  DEBOUTLN(CERR, "specialize_variable: called with: ", f);
355  for ( Extgenerator.reset(); Extgenerator.hasItems(); Extgenerator.next() ){
356    DEBOUTLN(CERR, "  specialize_variable: trying:  ", Extgenerator.item());
357    g= f( Extgenerator.item(), x );
358    DEBOUTLN(CERR, "  specialize_variable: resulting g= ", g);
359    if ( various_tests(g,deg,former_nr_of_variables - nr_of_variable ) ){
360      Substitutionlist.insert(SForm(x,Extgenerator.item())); // append (Var,value) pair
361      f= g;
362      return 1;
363    }
364  }
365  return 0;
366}
367static int
368specialize_agvariable( CanonicalForm & f, int deg, SFormList & Substitutionlist, int nr_of_variable,
369                     int former_nr_of_variables, AlgExtGenerator & Extgenerator ){
370  CanonicalForm g;
371  Variable x(nr_of_variable);
372
373  DEBOUTLN(CERR, "specialize_variable: called with: ", f);
374  for ( Extgenerator.reset(); Extgenerator.hasItems(); Extgenerator.next() ){
375    DEBOUTLN(CERR, "  specialize_variable: trying:  ", Extgenerator.item());
376    g= f( Extgenerator.item(), x );
377    DEBOUTLN(CERR, "  specialize_variable: resulting g= ", g);
378    if ( various_tests(g,deg,former_nr_of_variables - nr_of_variable ) ){
379      Substitutionlist.insert(SForm(x,Extgenerator.item())); // append (Var,value) pair
380      f= g;
381      return 1;
382    }
383  }
384  return 0;
385}
386
387///////////////////////////////////////////////////////////////
388// generate a minpoly of degree degree_of_Extension in the   //
389// field getCharacteristik()^Extension.                      //
390///////////////////////////////////////////////////////////////
391CanonicalForm
392generate_mipo( int degree_of_Extension , const Variable & Extension ){
393  FFRandom gen;
394  if (degree (Extension) < 0)
395    factoryError("libfac: evaluate: Extension not inFF() or inGF() !");
396  return find_irreducible( degree_of_Extension, gen, Variable(1) );
397}
398
399///////////////////////////////////////////////////////////////
400// Try to find a specialization for f over the field of char //
401// f.getCharacteristic() and (possible) extension defined by //
402// the variable Extension .                                  //
403// Returns 1 iff specialisation was found, 0 otherwise.      //
404// If 0 is returned there are variables left to substitute.  //
405// We check if Substitutionlist.length() > 0, i.e. we        //
406// previously tried to find specialization values for some   //
407// values. We take them and work with the resulting poly.    //
408///////////////////////////////////////////////////////////////
409static int
410try_specializePoly(const CanonicalForm & f, const Variable & Extension, int deg, SFormList & Substitutionlist, int ii,int j)
411{
412  int ok,i= ii;
413  CanonicalForm ff= f;
414
415  if ( Substitutionlist.length() > 0 ){ // we formerly tried to specialize
416    ff= substitutePoly(f,Substitutionlist); // substitute found values
417    i= Substitutionlist.length() + 1;
418  }
419
420  if ( degree(Extension) > 0 )
421  { // working over Extensions
422    DEBOUTLN(CERR, "try_specializePoly: working over Extensions: ", Extension);
423    if (Extension.level() > 0)
424    {
425    //  AlgExtGenerator g(Extension,minpoly );
426    //  for ( int k=i ; k<j ; k++ ) // try to find specialization for all
427    //  {                           // variables (# = k ) beginning with the
428    //                             // starting value i
429    //    ok= specialize_agvariable( ff, deg, Substitutionlist, k, j, g );
430    //    if ( ! ok ) return 0; // we failed
431    //  }
432      #ifndef NDEBUG
433        //printf("libfac: try_specializePoly: extension level >0\n");
434      #endif
435      return 0; // we failed
436    }
437    else
438    {
439      AlgExtGenerator g(Extension);
440      for ( int k=i ; k<j ; k++ ) // try to find specialization for all
441      {                           // variables (# = k ) beginning with the
442                                 // starting value i
443        ok= specialize_agvariable( ff, deg, Substitutionlist, k, j, g );
444        if ( ! ok ) return 0; // we failed
445      }
446    }
447  }
448  else{ // working over the ground-field
449    FFGenerator g;
450    DEBOUTMSG(CERR, "try_specializePoly: working over the ground-field.");
451    for ( int k=i ; k<j ; k++ ){
452      ok= specialize_variable( ff, deg, Substitutionlist, k, j, g );
453      if ( ! ok ) return 0; // we failed
454    }
455  }
456  return 1;
457}
458
459static int
460specializePoly(const CanonicalForm & f, Variable & Extension, int deg, SFormList & Substitutionlist, int i,int j)
461{
462  Variable minpoly= Extension;
463  int ok,extended= degree(Extension), working_over_extension;
464
465  // Remember if we are working over an extension-field
466  if ( extended >= 2 )    { working_over_extension = 1; }
467  else                    { working_over_extension = 0; extended = 1; }
468  // First try:
469  ok = try_specializePoly(f,minpoly,deg,Substitutionlist,i,j);
470  while ( ! ok ) // we have to extend!
471  {
472    SFormList origS=Substitutionlist;
473    extended+= 1;
474    if ( ! working_over_extension )
475    {
476      if (!hasMipo(Extension))
477        minpoly= rootOf (generate_mipo (extended, Extension));
478      else
479      {
480        setReduce (Extension, false);
481        setMipo (minpoly, generate_mipo ( extended, Extension));
482        setReduce (Extension, true);
483      }
484      Extension= minpoly;
485      ok= try_specializePoly(f,minpoly,deg,Substitutionlist,i,j);
486      if (!ok)
487        Substitutionlist=origS;
488    }
489    else
490    {
491      factoryError("libfac: spezializePoly ERROR: Working over given extension-field not yet implemented!");
492      return 0;
493    }
494  }
495  return 1;
496}
497
498
499// This is a procedure to play with: lot's of parameters!
500// returns: 0  iff no success (possibly because Extension isn't great enough
501//          >0 iff g (univariate) splits into n factors;
502// if n>0 BestEvaluationpoint contains the choice of values for the variables
503//
504// tries to find at least maxtries evaluation points
505// if g factored sametries into the same number of poly's the procedure stops
506// if we tried failtries evaluations not found valid, we stop. Perhaps
507// Extension isn't big enough!
508static int
509evaluate( int maxtries, int sametries, int failtries, const CanonicalForm &f , const Variable & Extension, const CanonicalForm &mipo, SFormList & BestEvaluationpoint, CFFList & BestFactorisation ){
510  int minfactors=degree(f),degf=degree(f),n=level(f.mvar())-1;
511  SFormList minEvaluation;
512  CFFList minFactorisation;
513  int samefactors=0, failedfactor=0, tried=0;
514  FFRandom gen;
515  CFFList unilist;
516
517  if ( degree(Extension) >0 ) GFRandom gen;
518  else
519  {
520    if ( degree(Extension) == 0 ) FFRandom gen;
521    else
522    {
523      factoryError("libfac: evaluate: Extension not inFF() or inGF() !");
524    }
525  }
526  REvaluation k(1,n,gen);
527  k.nextpoint();
528  for ( int i=1; i<=maxtries ; i++)
529  {
530    // k.nextpoint();
531    SFormList Substitutionlist;
532    for ( int j=1; j<=n; j++ )
533     Substitutionlist.insert(SForm(Variable(j),k[j]));
534    k.nextpoint();
535    CanonicalForm g=substitutePoly(f,Substitutionlist);
536    if ( various_tests(g, degf,1) )
537    { // found a valid point
538      failedfactor = 0; tried += 1;
539      if ( degree(Extension) == 0   )
540        unilist = factorize(g,1); // poly is sqr-free!
541      else
542      {
543        unilist = factorize2(g,Extension,mipo);
544      }
545      if (unilist.length() <= minfactors )
546      {
547        minfactors=unilist.length();
548        minEvaluation=Substitutionlist;
549        minFactorisation=unilist;
550      }
551      else samefactors +=1;
552
553      if (unilist.length() == 1 ) // wow! we found f is irreducible!
554      {
555        BestEvaluationpoint=minEvaluation;
556        BestFactorisation=minFactorisation;
557        return 1;
558      }
559
560      if ( samefactors >= sametries ) // now we stop ( maybe polynomial *has*
561                                      // minfactors factors? )
562      {
563        BestEvaluationpoint=minEvaluation;
564        BestFactorisation=minFactorisation;
565        return minfactors;
566      }
567
568    }
569    else
570      failedfactor += 1;
571
572    if ( failedfactor >= failtries ) // now we stop ( perhaps Extension isn't
573                                     // big enough )
574    {
575      if ( tried == 0 )
576        return 0;
577      else
578      {
579        BestEvaluationpoint=minEvaluation;
580        BestFactorisation=minFactorisation;
581        return minfactors;
582      }
583    }
584  }
585  BestEvaluationpoint=minEvaluation;
586  BestFactorisation=minFactorisation;
587  return minfactors;
588}
589
590///////////////////////////////////////////////////////////////
591// A factorization routine for a sqrfree polynomial.         //
592// Returns the list of factors.                              //
593///////////////////////////////////////////////////////////////
594CFFList
595Factorized( const CanonicalForm & F, const CanonicalForm & alpha, int Mainvar)
596{
597  CanonicalForm f,lt,ff,ffuni;
598  Variable Extension=alpha.mvar();
599  CFFList Outputlist,UnivariateFactorlist,Outputlist2;
600  SFormList Substitutionlist, Evaluationpoint;
601  CFFactor copy;
602  int mainvar=Mainvar,success,oldmainvar;
603  CFMap m;
604
605  // INTERRUPTHANDLER
606  if ( interrupt_handle() ) return CFFList() ;
607  // INTERRUPTHANDLER
608
609  if (F.inCoeffDomain()) {  return CFFList(CFFactor(F,1)); }
610  if ( F.isUnivariate() ) // could have lost one Variable elsewhere
611  {
612    if ( degree(Extension) == 0 )
613    {
614      TIMING_START(evaluate_time);
615      Outputlist = factorize(F,1); // poly is sqr-free
616      TIMING_END(evaluate_time);
617      return Outputlist;
618    }
619    else
620    {
621      if (Extension.level()<0)
622      DEBOUTLN(CERR, "Univ. Factorization over extension of degree ",
623               degree(getMipo(Extension,'x')) );
624      else
625      DEBOUTLN(CERR, "Univ. Factorization over extension of level ??",
626                Extension.level());
627      TIMING_START(evaluate_time);
628      Outputlist = factorize2(F,Extension,alpha);
629      TIMING_END(evaluate_time);
630      return Outputlist;
631    }
632  }
633
634  if ( Mainvar ) oldmainvar=Mainvar; else oldmainvar=level(F);
635  // First choose a main variable; this may be revisted in the next step
636  mainvar = choose_main_variable(F);
637  // Let`s look if @f/@mainvar is nonzero
638  mainvar = necessary_condition(F,mainvar);
639  // Now we have definetly choosen a main variable
640  // swap poly such that the mainvar has highest level
641  f=swapvar(F,mainvar,level(F));
642
643  // INTERRUPTHANDLER
644  if ( interrupt_handle() ) return CFFList() ;
645  // INTERRUPTHANDLER
646
647  if ( oldmainvar != mainvar ){
648    DEBOUTSL(CERR); DEBOUT(CERR,"Swapped poly ", F);
649    DEBOUT(CERR, " in ", f); DEBOUTNL(CERR);
650    DEBOUTSL(CERR); DEBOUT(CERR,"Swapped  ", oldmainvar );
651    DEBOUT(CERR, " <-- ", mainvar ); DEBOUT(CERR, "  Mainvar= ", f.mvar());
652    DEBOUTNL(CERR);
653    ff = f.deriv();
654    TIMING_START(discr_time);
655    ffuni = gcd(f,ff);
656    TIMING_END(discr_time);
657    if ( !(ffuni.isOne()) ){ //discriminante nonzero: split poly
658      DEBOUTLN(CERR,"Nontrivial GCD of f= ", f);
659      DEBOUTLN(CERR,"             and @f= ", ff);
660      DEBOUTLN(CERR,"          GCD(f,@f)= ", ffuni);
661      ff=f/ffuni;
662      CFFList Outputlist_a, Outputlist_b;
663      Outputlist_a = Factorized(ff,alpha);
664      DEBOUTLN(CERR, "Outputlist_a = ", Outputlist_a);
665      Outputlist_b = Factorized(ffuni,alpha);
666      DEBOUTLN(CERR, "Outputlist_b = ", Outputlist_b);
667      Outputlist = myUnion(Outputlist_a, Outputlist_b);
668      // have to back-swapvar the factors....
669      for ( CFFListIterator i=Outputlist; i.hasItem(); i++ ){
670        copy=i.getItem();
671        Outputlist2.append(CFFactor(swapvar(copy.factor(),oldmainvar,mainvar),copy.exp()));
672      }
673      DEBOUTLN(CERR, "Outputlist2 (a+b swapped) (to return) = ", Outputlist2);
674      return Outputlist2;
675    }
676  }
677
678  // Check special cases
679  for ( int i=1; i<=level(F); i++)
680  {
681    if ( degree(f,Variable(i) ) == 1 )
682    //test trivial case; only true iff F is primitiv w.r.t every variable; else check (if F=ax+b) gcd(a,b)=1 ?
683    {
684      DEBOUTLN(CERR, "Trivial case: ", F);
685      Outputlist.append(CFFactor(F,1));
686      return Outputlist;
687    }
688  }
689
690  // Look at the leading term:
691  lt = LC(f);
692  DEBOUTLN(CERR, "Leading term: ", lt);
693  //if ( lt != f.genOne() )
694  if ( !lt.isOne() )
695  {
696    // make the polynomial monic in the main variable
697    ff = make_monic(f,lt); ffuni = ff;
698    DEBOUTLN(CERR, "make_monic returned: ", ff);
699  }
700  else{ ff= f; ffuni= ff; }
701
702  TIMING_START(evaluate_time);
703  success=evaluate(min(10,max(degree(ff), 5)), min(degree(ff),3), min(degree(ff),3), ff, Extension, alpha, Substitutionlist,UnivariateFactorlist);
704  DEBOUTLN(CERR,  "Returned from evaluate: success: ", success);
705  for ( SFormListIterator ii=Substitutionlist; ii.hasItem(); ii++ )
706  {
707    DEBOUTLN(CERR, "Substituting ", ii.getItem().factor());
708    DEBOUTLN(CERR, "       with value: ", ii.getItem().exp());
709  }
710
711  if ( success==0 ) // evalute wasn't successfull
712  {
713    success= specializePoly(ffuni,Extension,degree(ff),Substitutionlist,1,getNumVars(compress(ff,m)));
714    DEBOUTLN(CERR,  "Returned from specializePoly: success: ", success);
715    if (success == 0 ) // No spezialisation could be found
716    {
717      factoryError("libfac: Factorize: ERROR: Not able to find a valid specialization!");
718      Outputlist.append(CFFactor(F,1));
719      return Outputlist;
720    }
721
722    // INTERRUPTHANDLER
723    if ( interrupt_handle() ) return CFFList() ;
724    // INTERRUPTHANDLER
725
726    ffuni = substitutePoly(ff,Substitutionlist);
727    // We now have an univariat poly; factorize that
728    if ( degree(Extension) == 0   )
729    {
730      DEBOUTMSG(CERR, "Univ. Factorization over the ground field");
731      UnivariateFactorlist = factorize(ffuni,1); // univ. poly is sqr-free!
732    }
733    else
734    {
735      DEBOUTLN(CERR, "Univ. Factorization over extension of degree ",
736               degree(getMipo(Extension,'x')) );
737      UnivariateFactorlist = factorize2(ffuni,Extension,alpha);
738    }
739  }
740  else
741  {
742    ffuni = substitutePoly(ff,Substitutionlist);
743  }
744    TIMING_END(evaluate_time);
745  if (UnivariateFactorlist.length() == 1)
746  { // poly is irreduzibel
747    DEBOUTLN(CERR, "Univ. poly is irreduzible: ", UnivariateFactorlist);
748    Outputlist.append(CFFactor(F,1));
749    return Outputlist;
750  }
751  else
752  { // we have factors
753    DEBOUTSL(CERR);
754    DEBOUT(CERR, "Univariate poly has " , UnivariateFactorlist.length());
755    DEBOUT(CERR, " factors:  ", ffuni);
756    DEBOUT(CERR, " = ", UnivariateFactorlist); DEBOUTNL(CERR);
757
758    // INTERRUPTHANDLER
759    if ( interrupt_handle() ) return CFFList() ;
760    // INTERRUPTHANDLER
761
762    TIMING_START(hensel_time);
763    Outputlist = MultiHensel(ff,UnivariateFactorlist,Substitutionlist, alpha);
764    DEBOUTLN(CERR, "Outputlist after MultiHensel: ", Outputlist);
765    TIMING_END(hensel_time);
766
767    // INTERRUPTHANDLER
768    if ( interrupt_handle() ) return CFFList() ;
769    // INTERRUPTHANDLER
770
771    TIMING_START(truefactor_time);
772    Outputlist = Truefactors(ff, level(ff), Substitutionlist, Outputlist);
773    DEBOUTLN(CERR, "Outputlist after Truefactors: ", Outputlist);
774    TIMING_END(truefactor_time);
775
776    // INTERRUPTHANDLER
777    if ( interrupt_handle() ) return CFFList() ;
778    // INTERRUPTHANDLER
779
780    //if ( lt != f.genOne() )
781    if ( !lt.isOne() )
782    {
783      Outputlist = not_monic(Outputlist,lt,ff,level(ff));
784      DEBOUTLN(CERR, "not_monic returned: ", Outputlist);
785    }
786
787    // have to back-swapvar the factors....
788    for ( CFFListIterator i=Outputlist; i.hasItem(); i++ )
789    {
790      copy=i.getItem();
791      Outputlist2.append(CFFactor(swapvar(copy.factor(),oldmainvar,mainvar),copy.exp()));
792    }
793
794    return Outputlist2;
795  }
796}
797
798int cmpCF( const CFFactor & f, const CFFactor & g );
799
800///////////////////////////////////////////////////////////////
801// The user front-end for a uni/multivariate factorization   //
802// routine. F needs not to be SqrFree.                       //
803// Option of * choosing a  main variable (n.y.i.)            //
804//           * choosing an algebraic extension (n.y.u.)      //
805//           * ensuring poly is sqrfree (n.y.i.)             //
806// use Factorize(F,alpha,is_SqrFree) if not over Zp[x]/Q[x]  //
807///////////////////////////////////////////////////////////////
808int find_mvar(const CanonicalForm &f);
809CFFList Factorize(const CanonicalForm & F, int is_SqrFree )
810{
811  //out_cf("Factorize ",F,"\n");
812  CFFList Outputlist;
813
814  // INTERRUPTHANDLER
815  if ( interrupt_handle() ) return CFFList() ;
816  // INTERRUPTHANDLER
817
818  DEBINCLEVEL(CERR, "Factorize");
819  DEBOUTLN(CERR, "Called with F= ", F);
820  if (( getCharacteristic() == 0 ) || (F.isUnivariate()))
821  { // char == 0
822    TIMING_START(factorize_time);
823    //CERR << "Factoring in char=0 of " << F << " = " << Outputlist << "\n";
824    Outputlist= factorize(F);
825    // Factorization in char=0 doesn't sometimes return at least two elements!!!
826    if ( getNumVars(Outputlist.getFirst().factor()) != 0 )
827      Outputlist.insert(CFFactor(1,1));
828    //CERR << "  Factorize in char=0: returning with: " << Outputlist << "\n";
829    TIMING_END(factorize_time);
830    DEBDECLEVEL(CERR, "Factorize");
831    TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
832    return Outputlist;
833  }
834  CFFList SqrFreeList,Intermediatelist,Outputlist2;
835  ListIterator<CFFactor> i,j;
836  CanonicalForm g=1,unit=1,r=1;
837  Variable minpoly; // dummy
838  int exp;
839  CFMap m;
840  TIMING_START(factorize_time);
841  // search an "optimal" main variavble
842  int mv=F.level();
843  if ((mv != LEVELBASE) /* && (! F.isUnivariate()) */)
844  {
845     mv=find_mvar(F);
846     if (mv!=F.level())
847     {
848       swapvar(F,Variable(mv),F.mvar());
849     }
850  }
851
852  ///////
853  // Maybe it`s better to add a sqrfree-test before?
854  // (If gcd is fast...)
855  ///////
856  //  if ( ! SqrFreeTest(F) ){
857  if ( ! is_SqrFree )
858  {
859    TIMING_START(sqrfree_time);
860    SqrFreeList = SqrFreeMV(F) ; // first sqrfree the polynomial
861    // don't use sqrFree(F), factory's internal sqrFree for multiv.
862    // Polynomials; it's wrong!! Ex.: char=p   f= x^p*(y+1);
863    // SqrFreeMV(f)= ( y+1, (x)^p ), sqrFree(f)= ( y+1 ) .
864    TIMING_END(sqrfree_time);
865
866    // INTERRUPTHANDLER
867    if ( interrupt_handle() ) return CFFList() ;
868    // INTERRUPTHANDLER
869
870  }
871  else
872    SqrFreeList.append(CFFactor(F,1));
873
874  DEBOUTLN(CERR, "SqrFreeMV= ", SqrFreeList);
875  for ( i=SqrFreeList; i.hasItem(); i++ )
876  {
877    DEBOUTLN(CERR, "Factor under consideration: ", i.getItem().factor());
878    // We need a compress on each list item ! Maybe we have less variables!
879    g =compress(i.getItem().factor(),m);
880    exp = i.getItem().exp();
881    if ( getNumVars(g) ==0 ) // a constant; Exp==1
882      Outputlist.append( CFFactor(g,1) ) ;
883    else// a real polynomial
884      if ( g.isUnivariate() )
885      {
886        //out_cf("univ. poly: ",g,"\n");
887        Intermediatelist=factorize(g,1); // poly is sqr-free!
888        for ( j=Intermediatelist; j.hasItem(); j++ )
889          //Normally j.getItem().exp() should be 1
890          Outputlist.append( CFFactor( m(j.getItem().factor()),exp*j.getItem().exp()));
891      }
892      else
893      { // multivariate polynomial
894        if ( g.isHomogeneous() )
895        {
896          DEBOUTLN(CERR, "Poly is homogeneous! : ", g);
897          // Now we can substitute one variable to 1, factorize and then
898          // look on the resulting factors and their monomials for
899          // backsubstitution of the substituted variable.
900          Intermediatelist = HomogFactor(g, minpoly, 0);
901        }
902        else // not homogeneous
903          Intermediatelist = Factorized(g, minpoly, 0);
904
905        // INTERRUPTHANDLER
906        if ( interrupt_handle() ) return CFFList() ;
907        // INTERRUPTHANDLER
908
909        for ( j=Intermediatelist; j.hasItem(); j++ )
910          //Normally j.getItem().exp() should be 1
911          Outputlist= myappend( Outputlist, CFFactor(m(j.getItem().factor()),exp*j.getItem().exp()));
912      }
913  }
914  g=1; unit=1;
915  DEBOUTLN(CERR, "Outputlist is ", Outputlist);
916  for ( i=Outputlist; i.hasItem(); i++ )
917    if ( level(i.getItem().factor()) > 0 )
918    {
919      unit = lc(i.getItem().factor());
920      if ( getNumVars(unit) == 0 )
921      { // a constant; possibly 1
922        Outputlist2.append(CFFactor(i.getItem().factor()/unit , i.getItem().exp()));
923        g *=power(i.getItem().factor()/unit,i.getItem().exp());
924      }
925      else
926      {
927        Outputlist2.append(i.getItem());
928        g *=power(i.getItem().factor(),i.getItem().exp());
929      }
930    }
931
932  r=F/g;
933  Outputlist2.insert(CFFactor(r,1));
934
935  if ((mv!=F.level()) && (! F.isUnivariate() ))
936  {
937    CFFListIterator J=Outputlist2;
938    for ( ; J.hasItem(); J++)
939    {
940      swapvar(J.getItem().factor(),Variable(mv),F.mvar());
941    }
942    swapvar(F,Variable(mv),F.mvar());
943  }
944  DEBDECLEVEL(CERR, "Factorize");
945  TIMING_END(factorize_time);
946
947  TIMING_PRINT(sqrfree_time, "\ntime used for sqrfree   : ");
948  TIMING_PRINT(discr_time, "time used for discriminante   : ");
949  TIMING_PRINT(evaluate_time, "time used for evaluation and univ. factorization  : ");
950  TIMING_PRINT(hensel_time, "time used for hensel-lift   : ");
951  TIMING_PRINT(truefactor_time, "time used for truefactors   : ");
952  TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
953
954  if(isOn(SW_USE_NTL_SORT)) Outputlist2.sort(cmpCF);
955
956  return Outputlist2;
957}
958
959///////////////////////////////////////////////////////////////
960// The user front-end for a uni/multivariate factorization   //
961// routine. F needs not to be SqrFree.                       //
962// Option of * choosing a  main variable (n.y.i.)            //
963//           * choosing an algebraic extension (n.y.u.)      //
964//           * ensuring poly is sqrfree (n.y.i.)             //
965///////////////////////////////////////////////////////////////
966static bool fdivides2(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &minpoly)
967{
968  if (!minpoly.isZero())
969  {
970  #if 0
971    Variable Alpha=minpoly.mvar();
972    Variable X=rootOf(minpoly);
973    CanonicalForm rF=replacevar(F,Alpha,X);
974    CanonicalForm rG=replacevar(G,Alpha,X);
975    return fdivides(rF,rG);;
976  #else
977    if (degree(F,F.mvar()) > degree(G,F.mvar())) return false;
978    return true;
979    //CanonicalForm a,b;
980    //mydivrem(G,F,a,b);
981    //if (b.isZero()) return true;
982    //else return false;
983  #endif
984  }
985  else
986   return fdivides(F,G);
987}
988
989CFFList
990Factorize(const CanonicalForm & F, const CanonicalForm & minpoly, int is_SqrFree )
991{
992  //out_cf("Factorize: F=",F,"\n");
993  //out_cf("           minpoly:",minpoly,"\n");
994  CFFList Outputlist,SqrFreeList,Intermediatelist,Outputlist2;
995  ListIterator<CFFactor> i,j;
996  CanonicalForm g=1,unit=1,r=1;
997  //Variable minpoly; // reserved (-> Factorisation over algebraic Extensions)
998  int exp;
999  CFMap m;
1000
1001  // INTERRUPTHANDLER
1002  if ( interrupt_handle() ) return CFFList() ;
1003  // INTERRUPTHANDLER
1004
1005  DEBINCLEVEL(CERR, "Factorize");
1006  DEBOUTLN(CERR, "Called with F= ", F);
1007  if ( getCharacteristic() == 0 )
1008  { // char == 0
1009    TIMING_START(factorize_time);
1010    //CERR << "Factoring in char=0 of " << F << " = " << Outputlist << "\n";
1011    #if 0
1012    // SHOULD: Outputlist= factorize(F,minpoly);
1013    Outputlist= factorize(F);
1014    #else
1015    if (!minpoly.isZero())
1016    {
1017      if ( F.isHomogeneous() )
1018      {
1019        DEBOUTLN(CERR, "Poly is homogeneous! : ", F);
1020        // Now we can substitute one variable to 1, factorize and then
1021        // look on the resulting factors and their monomials for
1022        // backsubstitution of the substituted variable.
1023        Outputlist=HomogFactor(F, minpoly, 0);
1024      }
1025      else
1026      {
1027        CFList as(minpoly);
1028        //CFFList sqF=sqrFree(F); // sqrFreeZ
1029        CFFList sqF=SqrFreeMV(F,minpoly);
1030        if (sqF.isEmpty()) sqF=sqrFree(F);
1031        CFFList G,H;
1032        CanonicalForm fac;
1033        int d;
1034        ListIterator<CFFactor> i,k;
1035        for ( i = sqF; i.hasItem(); ++i )
1036        {
1037          d = i.getItem().exp();
1038          fac = i.getItem().factor();
1039          int success=1;
1040          G = newfactoras( fac, as, success);
1041          for ( k = G; k.hasItem(); ++k )
1042          {
1043            fac = k.getItem().factor();
1044            int dd = k.getItem().exp();
1045            H.append( CFFactor( fac , d*dd ) );
1046          }
1047        }
1048        Outputlist = H;
1049      }
1050    }
1051    else // minpoly==0
1052      Outputlist=factorize(F);
1053    #endif
1054    // Factorization in char=0 doesn't sometimes return at least two elements!!!
1055    if ( getNumVars(Outputlist.getFirst().factor()) != 0 )
1056      Outputlist.insert(CFFactor(1,1));
1057    //CERR << "  Factorize in char=0: returning with: " << Outputlist << "\n";
1058    TIMING_END(factorize_time);
1059    DEBDECLEVEL(CERR, "Factorize");
1060    TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
1061    //out_cff(Outputlist);
1062    return Outputlist;
1063  }
1064  TIMING_START(factorize_time);
1065  // search an "optimal" main variavble
1066  int mv=F.level();
1067  if (mv != LEVELBASE && ! F.isUnivariate() )
1068  {
1069     mv=find_mvar(F);
1070     if (mv!=F.level())
1071     {
1072       swapvar(F,Variable(mv),F.mvar());
1073     }
1074  }
1075
1076  ///////
1077  // Maybe it`s better to add a sqrfree-test before?
1078  // (If gcd is fast...)
1079  ///////
1080  //  if ( ! SqrFreeTest(F) ){
1081  if ( ! is_SqrFree )
1082  {
1083    TIMING_START(sqrfree_time);
1084    SqrFreeList = SqrFreeMV(F, minpoly) ; // first sqrfree the polynomial
1085    // don't use sqrFree(F), factory's internal sqrFree for multiv.
1086    // Polynomials; it's wrong!! Ex.: char=p   f= x^p*(y+1);
1087    // SqrFreeMV(f)= ( y+1, (x)^p ), sqrFree(f)= ( y+1 ) .
1088    TIMING_END(sqrfree_time);
1089
1090    // INTERRUPTHANDLER
1091    if ( interrupt_handle() ) return CFFList() ;
1092    // INTERRUPTHANDLER
1093
1094  }
1095  else
1096    SqrFreeList.append(CFFactor(F,1));
1097  DEBOUTLN(CERR, "SqrFreeMV= ", SqrFreeList);
1098  for ( i=SqrFreeList; i.hasItem(); i++ )
1099  {
1100    DEBOUTLN(CERR, "Factor under consideration: ", i.getItem().factor());
1101    // We need a compress on each list item ! Maybe we have less variables!
1102    g =compress(i.getItem().factor(),m);
1103    exp = i.getItem().exp();
1104    if ( getNumVars(g) ==0 ) // a constant; Exp==1
1105      Outputlist.append( CFFactor(g,1) ) ;
1106    else// a real polynomial
1107    {
1108      if ( g.isUnivariate() )
1109      {
1110        Variable alpha=rootOf(minpoly);
1111        Intermediatelist=factorize2(g,alpha,minpoly); // poly is sqr-free!
1112        for ( j=Intermediatelist; j.hasItem(); j++ )
1113          //Normally j.getItem().exp() should be 1
1114          Outputlist.append(
1115           CFFactor( m(replacevar(j.getItem().factor(),alpha,minpoly.mvar())),
1116             exp*j.getItem().exp()));
1117      }
1118      else // multivariate polynomial
1119      {
1120        if ( g.isHomogeneous() )
1121        {
1122          DEBOUTLN(CERR, "Poly is homogeneous! : ", g);
1123          // Now we can substitute one variable to 1, factorize and then
1124          // look on the resulting factors and their monomials for
1125          // backsubstitution of the substituted variable.
1126          Intermediatelist = HomogFactor(g, minpoly, 0);
1127        }
1128        else // not homogeneous
1129          Intermediatelist = Factorized(g, minpoly, 0);
1130
1131        // INTERRUPTHANDLER
1132        if ( interrupt_handle() ) return CFFList() ;
1133        // INTERRUPTHANDLER
1134
1135        for ( j=Intermediatelist; j.hasItem(); j++ )
1136          //Normally j.getItem().exp() should be 1
1137          Outputlist= myappend( Outputlist, CFFactor(m(j.getItem().factor()),exp*j.getItem().exp()));
1138      }
1139    }
1140  }
1141  g=1; unit=1;
1142  DEBOUTLN(CERR, "Outputlist is ", Outputlist);
1143  for ( i=Outputlist; i.hasItem(); i++ )
1144    if ( level(i.getItem().factor()) > 0 )
1145    {
1146      unit = lc(i.getItem().factor());
1147      if ( getNumVars(unit) == 0 ){ // a constant; possibly 1
1148        Outputlist2.append(CFFactor(i.getItem().factor()/unit , i.getItem().exp()));
1149        g *=power(i.getItem().factor()/unit,i.getItem().exp());
1150      }
1151      else
1152      {
1153        Outputlist2.append(i.getItem());
1154        g *=power(i.getItem().factor(),i.getItem().exp());
1155      }
1156    }
1157
1158  r=F/g;
1159  Outputlist2.insert(CFFactor(r,1));
1160
1161  if ((mv!=F.level()) && (! F.isUnivariate() ))
1162  {
1163    CFFListIterator J=Outputlist2;
1164    for ( ; J.hasItem(); J++)
1165    {
1166      swapvar(J.getItem().factor(),Variable(mv),F.mvar());
1167    }
1168    swapvar(F,Variable(mv),F.mvar());
1169  }
1170
1171  DEBDECLEVEL(CERR, "Factorize");
1172  TIMING_END(factorize_time);
1173
1174  TIMING_PRINT(sqrfree_time, "\ntime used for sqrfree   : ");
1175  TIMING_PRINT(discr_time, "time used for discriminante   : ");
1176  TIMING_PRINT(evaluate_time, "time used for evaluation and univ. factorization  : ");
1177  TIMING_PRINT(hensel_time, "time used for hensel-lift   : ");
1178  TIMING_PRINT(truefactor_time, "time used for truefactors   : ");
1179  TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
1180
1181  if(isOn(SW_USE_NTL_SORT)) Outputlist2.sort(cmpCF);
1182
1183  //out_cff(Outputlist2);
1184  return Outputlist2;
1185}
1186
Note: See TracBrowser for help on using the repository browser.