source: git/libfac/factor/Factor.cc @ 5e1226

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