source: git/libfac/factor/Factor.cc @ 070c1b4

jengelh-datetimespielwiese
Last change on this file since 070c1b4 was 070c1b4, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: error about choosing field extensions (Tst/Short/bug_charseries2) update version in libfac and factory
  • Property mode set to 100644
File size: 39.6 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 ) GFRandom gen;
395  else {
396    if ( degree(Extension) == 0 ) FFRandom gen;
397    else
398    {
399      factoryError("libfac: evaluate: Extension not inFF() or inGF() !");
400    }
401  }
402  return find_irreducible( degree_of_Extension, gen, Variable(1) );
403}
404
405///////////////////////////////////////////////////////////////
406// Try to find a specialization for f over the field of char //
407// f.getCharacteristic() and (possible) extension defined by //
408// the variable Extension .                                  //
409// Returns 1 iff specialisation was found, 0 otherwise.      //
410// If 0 is returned there are variables left to substitute.  //
411// We check if Substitutionlist.length() > 0, i.e. we        //
412// previously tried to find specialization values for some   //
413// values. We take them and work with the resulting poly.    //
414///////////////////////////////////////////////////////////////
415static int
416try_specializePoly(const CanonicalForm & f, const Variable & Extension, int deg, SFormList & Substitutionlist, int ii,int j)
417{
418  int ok,i= ii;
419  CanonicalForm ff= f;
420
421  if ( Substitutionlist.length() > 0 ){ // we formerly tried to specialize
422    ff= substitutePoly(f,Substitutionlist); // substitute found values
423    i= Substitutionlist.length() + 1;
424  }
425
426  if ( degree(Extension) > 0 )
427  { // working over Extensions
428    DEBOUTLN(CERR, "try_specializePoly: working over Extensions: ", Extension);
429    if (Extension.level() > 0)
430    {
431    //  AlgExtGenerator g(Extension,minpoly );
432    //  for ( int k=i ; k<j ; k++ ) // try to find specialization for all
433    //  {                           // variables (# = k ) beginning with the
434    //                             // starting value i
435    //    ok= specialize_agvariable( ff, deg, Substitutionlist, k, j, g );
436    //    if ( ! ok ) return 0; // we failed
437    //  }
438      #ifndef NDEBUG
439        //printf("libfac: try_specializePoly: extension level >0\n");
440      #endif
441      return 0; // we failed
442    }
443    else
444    {
445      AlgExtGenerator g(Extension);
446      for ( int k=i ; k<j ; k++ ) // try to find specialization for all
447      {                           // variables (# = k ) beginning with the
448                                 // starting value i
449        ok= specialize_agvariable( ff, deg, Substitutionlist, k, j, g );
450        if ( ! ok ) return 0; // we failed
451      }
452    }
453  }
454  else{ // working over the ground-field
455    FFGenerator g;
456    DEBOUTMSG(CERR, "try_specializePoly: working over the ground-field.");
457    for ( int k=i ; k<j ; k++ ){
458      ok= specialize_variable( ff, deg, Substitutionlist, k, j, g );
459      if ( ! ok ) return 0; // we failed
460    }
461  }
462  return 1;
463}
464
465static int
466specializePoly(const CanonicalForm & f, Variable & Extension, int deg, SFormList & Substitutionlist, int i,int j)
467{
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  {
478    SFormList origS=Substitutionlist;
479    extended+= 1;
480    if ( ! working_over_extension )
481    {
482      minpoly= rootOf(generate_mipo( extended,Extension ));
483      Extension= minpoly;
484      ok= try_specializePoly(f,minpoly,deg,Substitutionlist,i,j);
485      if (!ok)
486      {
487        Substitutionlist=origS;
488        // very bad hack: TODO
489        // we want to remove the newly created variable minpoly:
490        extern char *var_names_ext;
491        var_names_ext[strlen(var_names_ext)]='\0';
492      }
493    }
494    else
495    {
496      factoryError("libfac: spezializePoly ERROR: Working over given extension-field not yet implemented!");
497      return 0;
498    }
499  }
500  return 1;
501}
502
503
504// This is a procedure to play with: lot's of parameters!
505// returns: 0  iff no success (possibly because Extension isn't great enough
506//          >0 iff g (univariate) splits into n factors;
507// if n>0 BestEvaluationpoint contains the choice of values for the variables
508//
509// tries to find at least maxtries evaluation points
510// if g factored sametries into the same number of poly's the procedure stops
511// if we tried failtries evaluations not found valid, we stop. Perhaps
512// Extension isn't big enough!
513static int
514evaluate( int maxtries, int sametries, int failtries, const CanonicalForm &f , const Variable & Extension, const CanonicalForm &mipo, SFormList & BestEvaluationpoint, CFFList & BestFactorisation ){
515  int minfactors=degree(f),degf=degree(f),n=level(f.mvar())-1;
516  SFormList minEvaluation;
517  CFFList minFactorisation;
518  int samefactors=0, failedfactor=0, tried=0;
519  FFRandom gen;
520  CFFList unilist;
521
522  if ( degree(Extension) >0 ) GFRandom gen;
523  else
524  {
525    if ( degree(Extension) == 0 ) FFRandom gen;
526    else
527    {
528      factoryError("libfac: evaluate: Extension not inFF() or inGF() !");
529    }
530  }
531  REvaluation k(1,n,gen);
532  k.nextpoint();
533  for ( int i=1; i<=maxtries ; i++)
534  {
535    // k.nextpoint();
536    SFormList Substitutionlist;
537    for ( int j=1; j<=n; j++ )
538     Substitutionlist.insert(SForm(Variable(j),k[j]));
539    k.nextpoint();
540    CanonicalForm g=substitutePoly(f,Substitutionlist);
541    if ( various_tests(g, degf,1) )
542    { // found a valid point
543      failedfactor = 0; tried += 1;
544      if ( degree(Extension) == 0   )
545        unilist = factorize(g,1); // poly is sqr-free!
546      else
547      {
548        unilist = factorize2(g,Extension,mipo);
549      }
550      if (unilist.length() <= minfactors )
551      {
552        minfactors=unilist.length();
553        minEvaluation=Substitutionlist;
554        minFactorisation=unilist;
555      }
556      else samefactors +=1;
557
558      if (unilist.length() == 1 ) // wow! we found f is irreducible!
559      {
560        BestEvaluationpoint=minEvaluation;
561        BestFactorisation=minFactorisation;
562        return 1;
563      }
564
565      if ( samefactors >= sametries ) // now we stop ( maybe polynomial *has*
566                                      // minfactors factors? )
567      {
568        BestEvaluationpoint=minEvaluation;
569        BestFactorisation=minFactorisation;
570        return minfactors;
571      }
572
573    }
574    else
575      failedfactor += 1;
576
577    if ( failedfactor >= failtries ) // now we stop ( perhaps Extension isn't
578                                     // big enough )
579    {
580      if ( tried == 0 )
581        return 0;
582      else
583      {
584        BestEvaluationpoint=minEvaluation;
585        BestFactorisation=minFactorisation;
586        return minfactors;
587      }
588    }
589  }
590  BestEvaluationpoint=minEvaluation;
591  BestFactorisation=minFactorisation;
592  return minfactors;
593}
594
595///////////////////////////////////////////////////////////////
596// A factorization routine for a sqrfree polynomial.         //
597// Returns the list of factors.                              //
598///////////////////////////////////////////////////////////////
599CFFList
600Factorized( const CanonicalForm & F, const CanonicalForm & alpha, int Mainvar)
601{
602  CanonicalForm f,lt,ff,ffuni;
603  Variable Extension=alpha.mvar();
604  CFFList Outputlist,UnivariateFactorlist,Outputlist2;
605  SFormList Substitutionlist, Evaluationpoint;
606  CFFactor copy;
607  int mainvar=Mainvar,success,oldmainvar;
608  CFMap m;
609
610  // INTERRUPTHANDLER
611  if ( interrupt_handle() ) return CFFList() ;
612  // INTERRUPTHANDLER
613
614  if (F.inCoeffDomain()) {  return CFFList(CFFactor(F,1)); }
615  if ( F.isUnivariate() ) // could have lost one Variable elsewhere
616  {
617    if ( degree(Extension) == 0 )
618    {
619      TIMING_START(evaluate_time);
620      Outputlist = factorize(F,1); // poly is sqr-free
621      TIMING_END(evaluate_time);
622      return Outputlist;
623    }
624    else
625    {
626      if (Extension.level()<0)
627      DEBOUTLN(CERR, "Univ. Factorization over extension of degree ",
628               degree(getMipo(Extension,'x')) );
629      else
630      DEBOUTLN(CERR, "Univ. Factorization over extension of level ??",
631                Extension.level());
632      TIMING_START(evaluate_time);
633      Outputlist = factorize2(F,Extension,alpha);
634      TIMING_END(evaluate_time);
635      return Outputlist;
636    }
637  }
638
639  if ( Mainvar ) oldmainvar=Mainvar; else oldmainvar=level(F);
640  // First choose a main variable; this may be revisted in the next step
641  mainvar = choose_main_variable(F);
642  // Let`s look if @f/@mainvar is nonzero
643  mainvar = necessary_condition(F,mainvar);
644  // Now we have definetly choosen a main variable
645  // swap poly such that the mainvar has highest level
646  f=swapvar(F,mainvar,level(F));
647
648  // INTERRUPTHANDLER
649  if ( interrupt_handle() ) return CFFList() ;
650  // INTERRUPTHANDLER
651
652  if ( oldmainvar != mainvar ){
653    DEBOUTSL(CERR); DEBOUT(CERR,"Swapped poly ", F);
654    DEBOUT(CERR, " in ", f); DEBOUTNL(CERR);
655    DEBOUTSL(CERR); DEBOUT(CERR,"Swapped  ", oldmainvar );
656    DEBOUT(CERR, " <-- ", mainvar ); DEBOUT(CERR, "  Mainvar= ", f.mvar());
657    DEBOUTNL(CERR);
658    ff = f.deriv();
659    TIMING_START(discr_time);
660    ffuni = gcd(f,ff);
661    TIMING_END(discr_time);
662    if ( !(ffuni.isOne()) ){ //discriminante nonzero: split poly
663      DEBOUTLN(CERR,"Nontrivial GCD of f= ", f);
664      DEBOUTLN(CERR,"             and @f= ", ff);
665      DEBOUTLN(CERR,"          GCD(f,@f)= ", ffuni);
666      ff=f/ffuni;
667      CFFList Outputlist_a, Outputlist_b;
668      Outputlist_a = Factorized(ff,alpha);
669      DEBOUTLN(CERR, "Outputlist_a = ", Outputlist_a);
670      Outputlist_b = Factorized(ffuni,alpha);
671      DEBOUTLN(CERR, "Outputlist_b = ", Outputlist_b);
672      Outputlist = myUnion(Outputlist_a, Outputlist_b);
673      // have to back-swapvar the factors....
674      for ( CFFListIterator i=Outputlist; i.hasItem(); i++ ){
675        copy=i.getItem();
676        Outputlist2.append(CFFactor(swapvar(copy.factor(),oldmainvar,mainvar),copy.exp()));
677      }
678      DEBOUTLN(CERR, "Outputlist2 (a+b swapped) (to return) = ", Outputlist2);
679      return Outputlist2;
680    }
681  }
682
683  // Check special cases
684  for ( int i=1; i<=level(F); i++)
685  {
686    if ( degree(f,Variable(i) ) == 1 )
687    //test trivial case; only true iff F is primitiv w.r.t every variable; else check (if F=ax+b) gcd(a,b)=1 ?
688    {
689      DEBOUTLN(CERR, "Trivial case: ", F);
690      Outputlist.append(CFFactor(F,1));
691      return Outputlist;
692    }
693  }
694
695  // Look at the leading term:
696  lt = LC(f);
697  DEBOUTLN(CERR, "Leading term: ", lt);
698  //if ( lt != f.genOne() )
699  if ( !lt.isOne() )
700  {
701    // make the polynomial monic in the main variable
702    ff = make_monic(f,lt); ffuni = ff;
703    DEBOUTLN(CERR, "make_monic returned: ", ff);
704  }
705  else{ ff= f; ffuni= ff; }
706
707  TIMING_START(evaluate_time);
708  success=evaluate(min(10,max(degree(ff), 5)), min(degree(ff),3), min(degree(ff),3), ff, Extension, alpha, Substitutionlist,UnivariateFactorlist);
709  DEBOUTLN(CERR,  "Returned from evaluate: success: ", success);
710  for ( SFormListIterator ii=Substitutionlist; ii.hasItem(); ii++ )
711  {
712    DEBOUTLN(CERR, "Substituting ", ii.getItem().factor());
713    DEBOUTLN(CERR, "       with value: ", ii.getItem().exp());
714  }
715
716  if ( success==0 ) // evalute wasn't successfull
717  {
718    success= specializePoly(ffuni,Extension,degree(ff),Substitutionlist,1,getNumVars(compress(ff,m)));
719    DEBOUTLN(CERR,  "Returned from specializePoly: success: ", success);
720    if (success == 0 ) // No spezialisation could be found
721    {
722      factoryError("libfac: Factorize: ERROR: Not able to find a valid specialization!");
723      Outputlist.append(CFFactor(F,1));
724      return Outputlist;
725    }
726
727    // INTERRUPTHANDLER
728    if ( interrupt_handle() ) return CFFList() ;
729    // INTERRUPTHANDLER
730
731    ffuni = substitutePoly(ff,Substitutionlist);
732    // We now have an univariat poly; factorize that
733    if ( degree(Extension) == 0   )
734    {
735      DEBOUTMSG(CERR, "Univ. Factorization over the ground field");
736      UnivariateFactorlist = factorize(ffuni,1); // univ. poly is sqr-free!
737    }
738    else
739    {
740      DEBOUTLN(CERR, "Univ. Factorization over extension of degree ",
741               degree(getMipo(Extension,'x')) );
742      UnivariateFactorlist = factorize2(ffuni,Extension,alpha);
743    }
744  }
745  else
746  {
747    ffuni = substitutePoly(ff,Substitutionlist);
748  }
749    TIMING_END(evaluate_time);
750  if (UnivariateFactorlist.length() == 1)
751  { // poly is irreduzibel
752    DEBOUTLN(CERR, "Univ. poly is irreduzible: ", UnivariateFactorlist);
753    Outputlist.append(CFFactor(F,1));
754    return Outputlist;
755  }
756  else
757  { // we have factors
758    DEBOUTSL(CERR);
759    DEBOUT(CERR, "Univariate poly has " , UnivariateFactorlist.length());
760    DEBOUT(CERR, " factors:  ", ffuni);
761    DEBOUT(CERR, " = ", UnivariateFactorlist); DEBOUTNL(CERR);
762
763    // INTERRUPTHANDLER
764    if ( interrupt_handle() ) return CFFList() ;
765    // INTERRUPTHANDLER
766
767    TIMING_START(hensel_time);
768    Outputlist = MultiHensel(ff,UnivariateFactorlist,Substitutionlist, alpha);
769    DEBOUTLN(CERR, "Outputlist after MultiHensel: ", Outputlist);
770    TIMING_END(hensel_time);
771
772    // INTERRUPTHANDLER
773    if ( interrupt_handle() ) return CFFList() ;
774    // INTERRUPTHANDLER
775
776    TIMING_START(truefactor_time);
777    Outputlist = Truefactors(ff, level(ff), Substitutionlist, Outputlist);
778    DEBOUTLN(CERR, "Outputlist after Truefactors: ", Outputlist);
779    TIMING_END(truefactor_time);
780
781    // INTERRUPTHANDLER
782    if ( interrupt_handle() ) return CFFList() ;
783    // INTERRUPTHANDLER
784
785    //if ( lt != f.genOne() )
786    if ( !lt.isOne() )
787    {
788      Outputlist = not_monic(Outputlist,lt,ff,level(ff));
789      DEBOUTLN(CERR, "not_monic returned: ", Outputlist);
790    }
791
792    // have to back-swapvar the factors....
793    for ( CFFListIterator i=Outputlist; i.hasItem(); i++ )
794    {
795      copy=i.getItem();
796      Outputlist2.append(CFFactor(swapvar(copy.factor(),oldmainvar,mainvar),copy.exp()));
797    }
798
799    return Outputlist2;
800  }
801}
802
803int cmpCF( const CFFactor & f, const CFFactor & g );
804
805///////////////////////////////////////////////////////////////
806// The user front-end for a uni/multivariate factorization   //
807// routine. F needs not to be SqrFree.                       //
808// Option of * choosing a  main variable (n.y.i.)            //
809//           * choosing an algebraic extension (n.y.u.)      //
810//           * ensuring poly is sqrfree (n.y.i.)             //
811// use Factorize(F,alpha,is_SqrFree) if not over Zp[x]/Q[x]  //
812///////////////////////////////////////////////////////////////
813int find_mvar(const CanonicalForm &f);
814CFFList Factorize(const CanonicalForm & F, int is_SqrFree )
815{
816  //out_cf("Factorize ",F,"\n");
817  CFFList Outputlist;
818
819  // INTERRUPTHANDLER
820  if ( interrupt_handle() ) return CFFList() ;
821  // INTERRUPTHANDLER
822
823  DEBINCLEVEL(CERR, "Factorize");
824  DEBOUTLN(CERR, "Called with F= ", F);
825  if (( getCharacteristic() == 0 ) || (F.isUnivariate()))
826  { // char == 0
827    TIMING_START(factorize_time);
828    //CERR << "Factoring in char=0 of " << F << " = " << Outputlist << "\n";
829    Outputlist= factorize(F);
830    // Factorization in char=0 doesn't sometimes return at least two elements!!!
831    if ( getNumVars(Outputlist.getFirst().factor()) != 0 )
832      Outputlist.insert(CFFactor(1,1));
833    //CERR << "  Factorize in char=0: returning with: " << Outputlist << "\n";
834    TIMING_END(factorize_time);
835    DEBDECLEVEL(CERR, "Factorize");
836    TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
837    return Outputlist;
838  }
839  CFFList SqrFreeList,Intermediatelist,Outputlist2;
840  ListIterator<CFFactor> i,j;
841  CanonicalForm g=1,unit=1,r=1;
842  Variable minpoly; // dummy
843  int exp;
844  CFMap m;
845  TIMING_START(factorize_time);
846  // search an "optimal" main variavble
847  int mv=F.level();
848  if ((mv != LEVELBASE) /* && (! F.isUnivariate()) */)
849  {
850     mv=find_mvar(F);
851     if (mv!=F.level())
852     {
853       swapvar(F,Variable(mv),F.mvar());
854     }
855  }
856
857  ///////
858  // Maybe it`s better to add a sqrfree-test before?
859  // (If gcd is fast...)
860  ///////
861  //  if ( ! SqrFreeTest(F) ){
862  if ( ! is_SqrFree )
863  {
864    TIMING_START(sqrfree_time);
865    SqrFreeList = SqrFreeMV(F) ; // first sqrfree the polynomial
866    // don't use sqrFree(F), factory's internal sqrFree for multiv.
867    // Polynomials; it's wrong!! Ex.: char=p   f= x^p*(y+1);
868    // SqrFreeMV(f)= ( y+1, (x)^p ), sqrFree(f)= ( y+1 ) .
869    TIMING_END(sqrfree_time);
870
871    // INTERRUPTHANDLER
872    if ( interrupt_handle() ) return CFFList() ;
873    // INTERRUPTHANDLER
874
875  }
876  else
877    SqrFreeList.append(CFFactor(F,1));
878
879  DEBOUTLN(CERR, "SqrFreeMV= ", SqrFreeList);
880  for ( i=SqrFreeList; i.hasItem(); i++ )
881  {
882    DEBOUTLN(CERR, "Factor under consideration: ", i.getItem().factor());
883    // We need a compress on each list item ! Maybe we have less variables!
884    g =compress(i.getItem().factor(),m);
885    exp = i.getItem().exp();
886    if ( getNumVars(g) ==0 ) // a constant; Exp==1
887      Outputlist.append( CFFactor(g,1) ) ;
888    else// a real polynomial
889      if ( g.isUnivariate() )
890      {
891        //out_cf("univ. poly: ",g,"\n");
892        Intermediatelist=factorize(g,1); // poly is sqr-free!
893        for ( j=Intermediatelist; j.hasItem(); j++ )
894          //Normally j.getItem().exp() should be 1
895          Outputlist.append( CFFactor( m(j.getItem().factor()),exp*j.getItem().exp()));
896      }
897      else
898      { // multivariate polynomial
899        if ( g.isHomogeneous() )
900        {
901          DEBOUTLN(CERR, "Poly is homogeneous! : ", g);
902          // Now we can substitute one variable to 1, factorize and then
903          // look on the resulting factors and their monomials for
904          // backsubstitution of the substituted variable.
905          Intermediatelist = HomogFactor(g, minpoly, 0);
906        }
907        else // not homogeneous
908          Intermediatelist = Factorized(g, minpoly, 0);
909
910        // INTERRUPTHANDLER
911        if ( interrupt_handle() ) return CFFList() ;
912        // INTERRUPTHANDLER
913
914        for ( j=Intermediatelist; j.hasItem(); j++ )
915          //Normally j.getItem().exp() should be 1
916          Outputlist= myappend( Outputlist, CFFactor(m(j.getItem().factor()),exp*j.getItem().exp()));
917      }
918  }
919  g=1; unit=1;
920  DEBOUTLN(CERR, "Outputlist is ", Outputlist);
921  for ( i=Outputlist; i.hasItem(); i++ )
922    if ( level(i.getItem().factor()) > 0 )
923    {
924      unit = lc(i.getItem().factor());
925      if ( getNumVars(unit) == 0 )
926      { // a constant; possibly 1
927        Outputlist2.append(CFFactor(i.getItem().factor()/unit , i.getItem().exp()));
928        g *=power(i.getItem().factor()/unit,i.getItem().exp());
929      }
930      else
931      {
932        Outputlist2.append(i.getItem());
933        g *=power(i.getItem().factor(),i.getItem().exp());
934      }
935    }
936
937  r=F/g;
938  Outputlist2.insert(CFFactor(r,1));
939
940  if ((mv!=F.level()) && (! F.isUnivariate() ))
941  {
942    CFFListIterator J=Outputlist2;
943    for ( ; J.hasItem(); J++)
944    {
945      swapvar(J.getItem().factor(),Variable(mv),F.mvar());
946    }
947    swapvar(F,Variable(mv),F.mvar());
948  }
949  DEBDECLEVEL(CERR, "Factorize");
950  TIMING_END(factorize_time);
951
952  TIMING_PRINT(sqrfree_time, "\ntime used for sqrfree   : ");
953  TIMING_PRINT(discr_time, "time used for discriminante   : ");
954  TIMING_PRINT(evaluate_time, "time used for evaluation and univ. factorization  : ");
955  TIMING_PRINT(hensel_time, "time used for hensel-lift   : ");
956  TIMING_PRINT(truefactor_time, "time used for truefactors   : ");
957  TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
958
959  if(isOn(SW_USE_NTL_SORT)) Outputlist2.sort(cmpCF);
960
961  return Outputlist2;
962}
963
964///////////////////////////////////////////////////////////////
965// The user front-end for a uni/multivariate factorization   //
966// routine. F needs not to be SqrFree.                       //
967// Option of * choosing a  main variable (n.y.i.)            //
968//           * choosing an algebraic extension (n.y.u.)      //
969//           * ensuring poly is sqrfree (n.y.i.)             //
970///////////////////////////////////////////////////////////////
971static bool fdivides2(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &minpoly)
972{
973  if (!minpoly.isZero())
974  {
975  #if 0
976    Variable Alpha=minpoly.mvar();
977    Variable X=rootOf(minpoly);
978    CanonicalForm rF=replacevar(F,Alpha,X);
979    CanonicalForm rG=replacevar(G,Alpha,X);
980    return fdivides(rF,rG);;
981  #else
982    if (degree(F,F.mvar()) > degree(G,F.mvar())) return false;
983    return true;
984    //CanonicalForm a,b;
985    //mydivrem(G,F,a,b);
986    //if (b.isZero()) return true;
987    //else return false;
988  #endif
989  }
990  else
991   return fdivides(F,G);
992}
993
994CFFList
995Factorize(const CanonicalForm & F, const CanonicalForm & minpoly, int is_SqrFree )
996{
997  //out_cf("Factorize: F=",F,"\n");
998  //out_cf("           minpoly:",minpoly,"\n");
999  CFFList Outputlist,SqrFreeList,Intermediatelist,Outputlist2;
1000  ListIterator<CFFactor> i,j;
1001  CanonicalForm g=1,unit=1,r=1;
1002  //Variable minpoly; // reserved (-> Factorisation over algebraic Extensions)
1003  int exp;
1004  CFMap m;
1005
1006  // INTERRUPTHANDLER
1007  if ( interrupt_handle() ) return CFFList() ;
1008  // INTERRUPTHANDLER
1009
1010  DEBINCLEVEL(CERR, "Factorize");
1011  DEBOUTLN(CERR, "Called with F= ", F);
1012  if ( getCharacteristic() == 0 )
1013  { // char == 0
1014    TIMING_START(factorize_time);
1015    //CERR << "Factoring in char=0 of " << F << " = " << Outputlist << "\n";
1016    #if 0
1017    // SHOULD: Outputlist= factorize(F,minpoly);
1018    Outputlist= factorize(F);
1019    #else
1020    if (!minpoly.isZero())
1021    {
1022      if ( F.isHomogeneous() )
1023      {
1024        DEBOUTLN(CERR, "Poly is homogeneous! : ", F);
1025        // Now we can substitute one variable to 1, factorize and then
1026        // look on the resulting factors and their monomials for
1027        // backsubstitution of the substituted variable.
1028        Outputlist=HomogFactor(F, minpoly, 0);
1029      }
1030      else
1031      {
1032        CFList as(minpoly);
1033        //CFFList sqF=sqrFree(F); // sqrFreeZ
1034        CFFList sqF=SqrFreeMV(F,minpoly);
1035        if (sqF.isEmpty()) sqF=sqrFree(F);
1036        CFFList G,H;
1037        CanonicalForm fac;
1038        int d;
1039        ListIterator<CFFactor> i,k;
1040        for ( i = sqF; i.hasItem(); ++i )
1041        {
1042          d = i.getItem().exp();
1043          fac = i.getItem().factor();
1044          int success=1;
1045          G = newfactoras( fac, as, success);
1046          for ( k = G; k.hasItem(); ++k )
1047          {
1048            fac = k.getItem().factor();
1049            int dd = k.getItem().exp();
1050            H.append( CFFactor( fac , d*dd ) );
1051          }
1052        }
1053        Outputlist = H;
1054      }
1055    }
1056    else // minpoly==0
1057      Outputlist=factorize(F);
1058    #endif
1059    // Factorization in char=0 doesn't sometimes return at least two elements!!!
1060    if ( getNumVars(Outputlist.getFirst().factor()) != 0 )
1061      Outputlist.insert(CFFactor(1,1));
1062    //CERR << "  Factorize in char=0: returning with: " << Outputlist << "\n";
1063    TIMING_END(factorize_time);
1064    DEBDECLEVEL(CERR, "Factorize");
1065    TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
1066    //out_cff(Outputlist);
1067    return Outputlist;
1068  }
1069  TIMING_START(factorize_time);
1070  // search an "optimal" main variavble
1071  int mv=F.level();
1072  if (mv != LEVELBASE && ! F.isUnivariate() )
1073  {
1074     mv=find_mvar(F);
1075     if (mv!=F.level())
1076     {
1077       swapvar(F,Variable(mv),F.mvar());
1078     }
1079  }
1080
1081  ///////
1082  // Maybe it`s better to add a sqrfree-test before?
1083  // (If gcd is fast...)
1084  ///////
1085  //  if ( ! SqrFreeTest(F) ){
1086  if ( ! is_SqrFree )
1087  {
1088    TIMING_START(sqrfree_time);
1089    SqrFreeList = SqrFreeMV(F, minpoly) ; // first sqrfree the polynomial
1090    // don't use sqrFree(F), factory's internal sqrFree for multiv.
1091    // Polynomials; it's wrong!! Ex.: char=p   f= x^p*(y+1);
1092    // SqrFreeMV(f)= ( y+1, (x)^p ), sqrFree(f)= ( y+1 ) .
1093    TIMING_END(sqrfree_time);
1094
1095    // INTERRUPTHANDLER
1096    if ( interrupt_handle() ) return CFFList() ;
1097    // INTERRUPTHANDLER
1098
1099  }
1100  else
1101    SqrFreeList.append(CFFactor(F,1));
1102  DEBOUTLN(CERR, "SqrFreeMV= ", SqrFreeList);
1103  for ( i=SqrFreeList; i.hasItem(); i++ )
1104  {
1105    DEBOUTLN(CERR, "Factor under consideration: ", i.getItem().factor());
1106    // We need a compress on each list item ! Maybe we have less variables!
1107    g =compress(i.getItem().factor(),m);
1108    exp = i.getItem().exp();
1109    if ( getNumVars(g) ==0 ) // a constant; Exp==1
1110      Outputlist.append( CFFactor(g,1) ) ;
1111    else// a real polynomial
1112    {
1113      if ( g.isUnivariate() )
1114      {
1115        Variable alpha=rootOf(minpoly);
1116        Intermediatelist=factorize2(g,alpha,minpoly); // poly is sqr-free!
1117        for ( j=Intermediatelist; j.hasItem(); j++ )
1118          //Normally j.getItem().exp() should be 1
1119          Outputlist.append(
1120           CFFactor( m(replacevar(j.getItem().factor(),alpha,minpoly.mvar())),
1121             exp*j.getItem().exp()));
1122      }
1123      else // multivariate polynomial
1124      {
1125        if ( g.isHomogeneous() )
1126        {
1127          DEBOUTLN(CERR, "Poly is homogeneous! : ", g);
1128          // Now we can substitute one variable to 1, factorize and then
1129          // look on the resulting factors and their monomials for
1130          // backsubstitution of the substituted variable.
1131          Intermediatelist = HomogFactor(g, minpoly, 0);
1132        }
1133        else // not homogeneous
1134          Intermediatelist = Factorized(g, minpoly, 0);
1135
1136        // INTERRUPTHANDLER
1137        if ( interrupt_handle() ) return CFFList() ;
1138        // INTERRUPTHANDLER
1139
1140        for ( j=Intermediatelist; j.hasItem(); j++ )
1141          //Normally j.getItem().exp() should be 1
1142          Outputlist= myappend( Outputlist, CFFactor(m(j.getItem().factor()),exp*j.getItem().exp()));
1143      }
1144    }
1145  }
1146  g=1; unit=1;
1147  DEBOUTLN(CERR, "Outputlist is ", Outputlist);
1148  for ( i=Outputlist; i.hasItem(); i++ )
1149    if ( level(i.getItem().factor()) > 0 )
1150    {
1151      unit = lc(i.getItem().factor());
1152      if ( getNumVars(unit) == 0 ){ // a constant; possibly 1
1153        Outputlist2.append(CFFactor(i.getItem().factor()/unit , i.getItem().exp()));
1154        g *=power(i.getItem().factor()/unit,i.getItem().exp());
1155      }
1156      else
1157      {
1158        Outputlist2.append(i.getItem());
1159        g *=power(i.getItem().factor(),i.getItem().exp());
1160      }
1161    }
1162
1163  r=F/g;
1164  Outputlist2.insert(CFFactor(r,1));
1165
1166  if ((mv!=F.level()) && (! F.isUnivariate() ))
1167  {
1168    CFFListIterator J=Outputlist2;
1169    for ( ; J.hasItem(); J++)
1170    {
1171      swapvar(J.getItem().factor(),Variable(mv),F.mvar());
1172    }
1173    swapvar(F,Variable(mv),F.mvar());
1174  }
1175
1176  DEBDECLEVEL(CERR, "Factorize");
1177  TIMING_END(factorize_time);
1178
1179  TIMING_PRINT(sqrfree_time, "\ntime used for sqrfree   : ");
1180  TIMING_PRINT(discr_time, "time used for discriminante   : ");
1181  TIMING_PRINT(evaluate_time, "time used for evaluation and univ. factorization  : ");
1182  TIMING_PRINT(hensel_time, "time used for hensel-lift   : ");
1183  TIMING_PRINT(truefactor_time, "time used for truefactors   : ");
1184  TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
1185
1186  if(isOn(SW_USE_NTL_SORT)) Outputlist2.sort(cmpCF);
1187
1188  //out_cff(Outputlist2);
1189  return Outputlist2;
1190}
1191
Note: See TracBrowser for help on using the repository browser.