source: git/factory/libfac/factor/Factor.cc @ 6ce030f

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