source: git/libfac/factor/Factor.cc @ 7d889c

spielwiese
Last change on this file since 7d889c was 7d889c, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: moved WerrorS from C++ to C (Factor.cc MVMultiHensel.cc SqrFree.cc Truefactor.cc) git-svn-id: file:///usr/local/Singular/svn/trunk@910 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 29.1 KB
Line 
1/* Copyright 1996 Michael Messollen. All rights reserved. */
2///////////////////////////////////////////////////////////////////////////////
3static char * rcsid = "$Id: Factor.cc,v 1.4 1997-11-18 16:39:04 Singular Exp $ ";
4static char * errmsg = "\nYou found a bug!\nPlease inform (Michael Messollen) michael@math.uni-sb.de \nPlease include above information and your input (the ideal/polynomial and characteristic) in your bug-report.\nThank you.";
5///////////////////////////////////////////////////////////////////////////////
6// FACTORY - Includes
7#include <factory.h>
8// Factor - Includes
9#include "tmpl_inst.h"
10#include "SqrFree.h"
11#include "helpstuff.h"
12#include "MVMultiHensel.h"
13#include "Truefactor.h"
14#include "homogfactor.h"
15#include "interrupt.h"
16// some CC's need this:
17#include "Factor.h"
18
19#ifdef SINGULAR
20#  define HAVE_SINGULAR
21   extern "C" { void WerrorS(char *); }
22#endif
23
24#ifdef FACTORDEBUG
25#  define DEBUGOUTPUT
26#else
27#  undef DEBUGOUTPUT
28#endif
29
30#include "debug.h"
31#include "timing.h"
32TIMING_DEFINE_PRINT(factorize_time);
33TIMING_DEFINE_PRINT(sqrfree_time);
34TIMING_DEFINE_PRINT(discr_time);
35TIMING_DEFINE_PRINT(evaluate_time);
36TIMING_DEFINE_PRINT(hensel_time);
37TIMING_DEFINE_PRINT(truefactor_time);
38
39
40///////////////////////////////////////////////////////////////
41// Choose a main variable if the user didn`t wish a          //
42// special one. Returns level of main variable.              //
43///////////////////////////////////////////////////////////////
44static int 
45choose_main_variable( const CanonicalForm & f, int Mainvar=0){
46  CanonicalForm remlc, newlc;
47  int n= level(f), mainvar= Mainvar;
48
49  if (mainvar != 0) return mainvar ; // We force use of the wished mainvar
50  remlc= LC(f,n); mainvar = n;
51  if ( totaldegree(remlc)==0 ){ remlc=f.genOne() ; }
52  DEBOUTLN(cout, "remlc= " , remlc);
53  for ( int i=n-1; i>=1; i-- ){
54    newlc= LC(f,i);
55    if ( totaldegree(newlc)==0 ){ newlc=f.genOne() ; }
56    DEBOUTLN(cout, "newlc= " , newlc);
57    if ( (remlc.isOne()) && (newlc.isOne()) ){ // take care of the degrees
58      if ( degree(f,i) < degree(f,mainvar) ){
59        remlc= newlc; 
60        mainvar= i;
61      }
62    }
63    else  if ( (! remlc.isOne() ) && ( newlc.isOne() ) ){
64      remlc= newlc; 
65      mainvar= i;
66    }
67  }
68  return mainvar;
69}
70
71///////////////////////////////////////////////////////////////
72// Check if the derivative is nonzero for oldmainvar.        //
73// Returns the level of the choosen main variable.           //
74///////////////////////////////////////////////////////////////
75static int 
76necessary_condition( const CanonicalForm & F, int oldmainvar){
77  CanonicalForm g;
78  int n=level(F);
79
80  g= swapvar(F,oldmainvar,n); 
81  g= g.deriv();
82  if ( g.isZero() ) 
83    for ( int i=n; i>=1; i-- ){
84      g= swapvar(F,i,n); 
85      g= g.deriv();
86      if ( ! g.isZero() ) return i;
87    }
88  return oldmainvar;
89}
90
91///////////////////////////////////////////////////////////////
92// Make F monic. Return monic polynomial.                    //
93///////////////////////////////////////////////////////////////
94static CanonicalForm
95make_monic( const CanonicalForm & F, const CanonicalForm & lt){
96  CanonicalForm intermediatpoly,f;
97  Variable x(level(F));
98
99  if ( degree(lt) == 0 ) f= 1/lt * F ;
100  else {
101    intermediatpoly= power(lt,degree(F)-1);
102    for ( int i=0; i<=degree(F); i++ )
103      if ( ! F[i].isZero())
104        f+= (F[i] * intermediatpoly*power(x,i))/power(lt,i);
105  }
106  return f;
107}
108
109///////////////////////////////////////////////////////////////
110// Decide whether num/denum (num,denum both from the         //
111// FiniteFielddomain)  lies in the RationalDomain.           //
112// If false, return num/denum else return the zero poly from //
113// the FiniteFielddomain.                                    //
114///////////////////////////////////////////////////////////////
115static CanonicalForm
116is_rational( const CanonicalForm & num, const CanonicalForm & denum ){
117  CanonicalForm a, b;
118  int retvalue;
119
120  retvalue= mydivremt(num,denum,a,b);
121  if ( retvalue && b == num.genZero() ) // num/denum from FFdomain
122    return a;
123  else // num/denum is rational
124    return num.genZero();
125}
126
127///////////////////////////////////////////////////////////////
128// lt_is_product returns 1 iff lt is a product, 0 iff lt is  //
129// a sum.                                                    //
130///////////////////////////////////////////////////////////////
131static int 
132lt_is_product( const CanonicalForm & lt ){
133  CFList result;
134
135  result= get_Terms(lt);
136  if ( result.length() > 1 ) return 0;
137  else return 1;
138}
139
140///////////////////////////////////////////////////////////////
141// Reverse the make_monic transformation.                    //
142// Return the list of factors.                               //
143///////////////////////////////////////////////////////////////
144static CFFList
145not_monic( const CFFList & TheList, const CanonicalForm & ltt, const CanonicalForm & F, int levelF){
146  CFFList Returnlist,IntermediateList;
147  CFFListIterator i;
148  CanonicalForm intermediate,lt= ltt,savelc;
149  CanonicalForm numerator,denumerator,test,a,b;
150  Variable x(level(F));
151  int test1;
152
153  if ( lt == lt.genOne() ) return TheList; // the poly was already monic
154  if ( TheList.length() <= 1 ){ // only one factor to substitute back
155    if ( totaldegree(lt) == 0 ) // lt is type numeric
156      Returnlist.append( CFFactor(lt*TheList.getFirst().factor(),
157                                  TheList.getFirst().exp()) );
158    else { 
159      intermediate = F(x*lt, levelF)/power(lt,degree(F,levelF)-1);
160      Returnlist.append(CFFactor(intermediate,TheList.getFirst().exp()));
161    }
162  }
163  else { // more then one factor
164    IntermediateList= TheList;
165    if ( totaldegree(lt) == 0 ){ // lt is type numeric;(SqrFree-use, see above)
166      Returnlist.append( CFFactor(lt*IntermediateList.getFirst().factor()
167                                  , IntermediateList.getFirst().exp()) );
168      IntermediateList.removeFirst();
169      Returnlist= Union(Returnlist,IntermediateList);
170    }
171    else{ // lt is a) a product or b) a sum of terms
172      if ( lt_is_product(lt) ){ // case a)
173        DEBOUTLN(cout, "lt_is_product: ", lt);
174        savelc= content(lt) ; // can we simplify to savelc= lc(lt); ?
175        while ( getNumVars(savelc) != 0 )
176          savelc= content(savelc);
177        for ( i=TheList; i.hasItem();i++ ){
178          numerator= i.getItem().factor();
179          numerator= numerator(x*lt,levelF); // x <- x*lt
180          denumerator= power(lt,degree(F,levelF)-1); // == lt^(1-degree(F,x)
181          while (numerator.genZero() == is_rational(numerator, denumerator))
182            numerator*= lt;
183          intermediate= is_rational(numerator,denumerator);
184
185          Returnlist.append( CFFactor(lc(content(intermediate))*intermediate/content(intermediate), i.getItem().exp() ) );
186        }
187        // Now we add a test. If product(factors)/F is a multiple of
188        // savelc, we have to add 1/multiplicity to the factors
189        IntermediateList= Returnlist;
190        intermediate= 1;
191        for ( CFFListIterator j=IntermediateList; j.hasItem(); j++)
192          intermediate*= j.getItem().factor();
193        test1= mydivremt( intermediate,F,a,b);
194        if ( test1 && b == intermediate.genZero() ) { // Yupp!
195          IntermediateList.append(CFFactor(1/a,1));
196          Returnlist= IntermediateList;
197        }
198        else { Returnlist= IntermediateList; }
199      }
200      else{ // case b)
201        DEBOUTLN(cout, "lt_is_sum: ", lt);
202        CanonicalForm save_denumerator= 1;
203        for ( i=TheList; i.hasItem(); i++ ){
204          numerator= i.getItem().factor();
205          numerator= numerator(x*lt,levelF); // x <- x*lt
206          denumerator= power(lt,degree(numerator,levelF)); // == lt^(-degree(numerator,x)
207          test= content(numerator,x); 
208          test1= mydivremt(denumerator,test,a,b);
209          if ( test1 && b == numerator.genZero() ){ // Yupp!
210            save_denumerator*= a;
211            Returnlist.append(CFFactor(numerator/test ,1));
212          }
213          else { 
214#ifdef HAVE_SINGULAR
215            WerrorS("libfac: ERROR: not_monic1: case lt is a sum.");
216#else
217            cerr << "libfac: ERROR: not_monic1: case lt is a sum.\n" 
218                 << rcsid << errmsg << endl;
219#endif
220          } 
221        }
222        // Now we add a test if we did the right thing:
223        // save_denumerator should be a multiple of the leading coeff
224        test1= mydivremt(save_denumerator,lt,a,b);
225        if ( test1 && b == save_denumerator.genZero() ) // Yupp!
226          // We have to multiply one of the factors with
227          // the multiplicity of the save_denumerator <-> lc
228          // the following will do what we want
229          Returnlist= myUnion( CFFList(CFFactor(1/a,1)),Returnlist) ;
230        else {
231#ifdef HAVE_SINGULAR
232          WerrorS("libfac: ERROR: not_monic2: case lt is a sum.");
233#else
234          cerr << "libfac: ERROR: not_monic2: case lt is a sum.\n" 
235               << rcsid << errmsg << endl;
236#endif
237        } 
238      }
239    }
240  }
241  DEBOUTLN(cout,"Returnlist: ", Returnlist);
242  return Returnlist;
243}
244
245///////////////////////////////////////////////////////////////
246// Substitute the (Variable,Value)-Pair(s) from Substitution-//
247// list into the polynomial F. Returns the resulting poly.   //
248///////////////////////////////////////////////////////////////
249static CanonicalForm
250substitutePoly( const CanonicalForm & F, const SFormList & Substitutionlist){
251  CanonicalForm f= F;
252
253  for ( SFormListIterator i=Substitutionlist; i.hasItem(); i++ )
254    f= f(i.getItem().exp(),level(i.getItem().factor()));
255  return f;
256}
257
258///////////////////////////////////////////////////////////////
259// Find specialization values for the poly F. Returns 0 if   //
260// procedure failed, 1 otherwise. On success Substitutionlist//
261// holds (Variable,Value)-pairs. On failure we only have a   //
262// partitial list.                                           //
263///////////////////////////////////////////////////////////////
264//      *** This is the version with extensions ***          //
265///////////////////////////////////////////////////////////////
266
267///////////////////////////////////////////////////////////////
268// is CF g ok?                                               //
269///////////////////////////////////////////////////////////////
270static int
271various_tests( const CanonicalForm & g, int deg, int vars_left){
272  CFMap m;
273
274  if ( degree(g) == deg ) // degrees match
275    if ( level(compress(g,m)) == (vars_left) ) // exactly one variable less
276      if ( SqrFreeTest(g,1) ) // poly is sqrfree
277        if ( mygcd(g,g.deriv()) == 1 ) // Discriminante != 0
278           return 1;
279  return 0;
280}
281
282///////////////////////////////////////////////////////////////
283// specialize one variable over the given field.             //
284///////////////////////////////////////////////////////////////
285// substitutes in poly f of degree deg with former
286// former_nr_of_variables variables the variable nr_of_variable ;
287// this is done in the field of Char getCharacteristic() and
288// Extension given by Extgenerator.
289///////////////////////////////////////////////////////////////
290static int
291specialize_variable( CanonicalForm & f, int deg, SFormList & Substitutionlist, int nr_of_variable, int former_nr_of_variables, CFGenerator & Extgenerator ){
292  CanonicalForm g;
293  Variable x(nr_of_variable);
294
295  DEBOUTLN(cout, "specialize_variable: called with: ", f);
296  for ( Extgenerator.reset(); Extgenerator.hasItems(); Extgenerator.next() ){
297    DEBOUTLN(cout, "  specialize_variable: trying:  ", Extgenerator.item());
298    g= f( Extgenerator.item(), x );
299    DEBOUTLN(cout, "  specialize_variable: resulting g= ", g);
300    if ( various_tests(g,deg,former_nr_of_variables - nr_of_variable ) ){ 
301      Substitutionlist.insert(SForm(x,Extgenerator.item())); // append (Var,value) pair
302      f= g;
303      return 1;
304    }
305  }
306  return 0;
307}
308
309///////////////////////////////////////////////////////////////
310// generate a minpoly of degree degree_of_Extension in the   //
311// field getCharacteristik()^Extension.                      //
312///////////////////////////////////////////////////////////////
313CanonicalForm
314generate_mipo( int degree_of_Extension , const Variable & Extension ){
315  FFRandom gen; 
316  if ( degree(Extension) > 0 ) GFRandom gen; 
317  else {
318    if ( degree(Extension) == 0 ) FFRandom gen;
319    else { 
320#ifdef HAVE_SINGULAR
321    WerrorS("libfac: evaluate: Extension not inFF() or inGF() !");
322#else
323    cerr << "libfac: evaluate: " << Extension << " not inFF() or inGF() !" 
324         << endl;
325#endif
326    FFRandom gen;
327    }
328  }
329  return find_irreducible( degree_of_Extension, gen, Variable(1) );
330}
331
332///////////////////////////////////////////////////////////////
333// Try to find a specialization for f over the field of char //
334// f.getCharacteristic() and (possible) extension defined by //
335// the variable Extension .                                  //
336// Returns 1 iff specialisation was found, 0 otherwise.      //
337// If 0 is returned there are variables left to substitute.  //
338// We check if Substitutionlist.length() > 0, i.e. we        //
339// previously tried to find specialization values for some   //
340// values. We take them and work with the resulting poly.    //
341///////////////////////////////////////////////////////////////
342static int
343try_specializePoly(const CanonicalForm & f, const Variable & Extension, int deg, SFormList & Substitutionlist, int ii,int j){
344  int ok,i= ii;
345  CanonicalForm ff= f;
346
347  if ( Substitutionlist.length() > 0 ){ // we formerly tried to specialize
348    ff= substitutePoly(f,Substitutionlist); // substitute found values
349    i= Substitutionlist.length() + 1;
350  }
351
352  if ( degree(Extension) > 0 ){ // working over Extensions
353    DEBOUTLN(cout, "try_specializePoly: working over Extensions: ", Extension);
354    AlgExtGenerator g(Extension);
355    for ( int k=i ; k<j ; k++ ){ // try to find specialization for all
356                                 // variables (# = k ) beginning with the
357                                 // starting value i
358      ok= specialize_variable( ff, deg, Substitutionlist, k, j, g );
359      if ( ! ok ) return 0; // we failed
360    }
361  }
362  else{ // working over the ground-field
363    FFGenerator g;
364    DEBOUTMSG(cout, "try_specializePoly: working over the ground-field.");
365    for ( int k=i ; k<j ; k++ ){
366      ok= specialize_variable( ff, deg, Substitutionlist, k, j, g );
367      if ( ! ok ) return 0; // we failed
368    }
369  }
370  return 1;
371}
372
373static int
374specializePoly(const CanonicalForm & f, Variable & Extension, int deg, SFormList & Substitutionlist, int i,int j){
375  Variable minpoly= Extension;
376  int ok,extended= degree(Extension), working_over_extension;
377
378  // Remember if we are working over an extension-field
379  if ( extended >= 2 )    { working_over_extension = 1; }
380  else                    { working_over_extension = 0; extended = 1; }
381  // First try:
382  ok = try_specializePoly(f,minpoly,deg,Substitutionlist,i,j);
383  while ( ! ok ){ // we have to extend!
384    extended+= 1;
385    if ( ! working_over_extension ){
386      minpoly= rootOf(generate_mipo( extended,Extension ));
387      Extension= minpoly;
388      ok= try_specializePoly(f,minpoly,deg,Substitutionlist,i,j);
389    }
390    else {
391#ifdef HAVE_SINGULAR
392      WerrorS("libfac: spezializePoly ERROR: Working over given extension-field not yet implemented!");
393#else
394      cerr << "libfac: spezializePoly ERROR: Working over given extension-field not yet implemented!\n" 
395           << rcsid << errmsg << endl;
396#endif
397      return 0;
398    }
399  }
400  return 1;
401}
402
403
404// This is a procedure to play with: lot's of parameters!
405// returns: 0  iff no success (possibly because Extension isn't great enough
406//          >0 iff g (univariate) splits into n factors;
407// if n>0 BestEvaluationpoint contains the choice of values for the variables
408//
409// tries to find at least maxtries evaluation points
410// if g factored sametries into the same number of poly's the procedure stops
411// if we tried failtries evaluations not found valid, we stop. Perhaps
412// Extension isn't big enough!
413static int
414evaluate( int maxtries, int sametries, int failtries, const CanonicalForm &f , const Variable & Extension, SFormList & BestEvaluationpoint, CFFList & BestFactorisation ){
415  int minfactors=degree(f),degf=degree(f),n=level(f.mvar())-1;
416  SFormList minEvaluation;
417  CFFList minFactorisation;
418  int samefactors=0, failedfactor=0, tried=0;
419  FFRandom gen;
420  CFFList unilist;
421
422  if ( degree(Extension) >0 ) GFRandom gen; 
423  else { if ( degree(Extension) == 0 ) FFRandom gen;
424  else { 
425#ifdef HAVE_SINGULAR
426    WerrorS("libfac: evaluate: Extension not inFF() or inGF() !");
427#else
428    cerr << "libfac: evaluate: " << Extension << " not inFF() or inGF() !" 
429         << endl;
430#endif
431    FFRandom gen; }}
432  REvaluation k(1,n,gen);
433  k.nextpoint();
434  for ( int i=1; i<=maxtries ; i++){
435    // k.nextpoint();
436    SFormList Substitutionlist;
437    for ( int j=1; j<=n; j++ )
438     Substitutionlist.insert(SForm(Variable(j),k[j])); 
439    k.nextpoint();
440    CanonicalForm g=substitutePoly(f,Substitutionlist);
441    if ( various_tests(g, degf,1) ){ // found a valid point
442      failedfactor = 0; tried += 1;
443      if ( degree(Extension) == 0   )
444        unilist = factorize(g,1); // poly is sqr-free!
445      else
446        unilist = factorize(g,Extension);
447      if (unilist.length() <= minfactors ) {
448        minfactors=unilist.length();
449        minEvaluation=Substitutionlist;
450        minFactorisation=unilist;
451      }
452      else samefactors +=1;
453
454      if (unilist.length() == 1 ){ // wow! we found f is irreducible!
455        BestEvaluationpoint=minEvaluation;
456        BestFactorisation=minFactorisation;
457        return 1;
458      }
459
460      if ( samefactors >= sametries ){ // now we stop ( maybe polynomial *has*
461                                       // minfactors factors? )
462        BestEvaluationpoint=minEvaluation;
463        BestFactorisation=minFactorisation;
464        return minfactors;
465      }
466
467    }
468    else failedfactor += 1;
469
470    if ( failedfactor >= failtries ){ // now we stop ( perhaps Extension isn't
471                                      // big enough )
472      if ( tried == 0 )
473        return 0;
474      else{
475        BestEvaluationpoint=minEvaluation;
476        BestFactorisation=minFactorisation;
477        return minfactors;
478      }
479    }
480
481  }
482  BestEvaluationpoint=minEvaluation;
483  BestFactorisation=minFactorisation;
484  return minfactors;
485}
486
487#ifdef EXPERIMENTAL
488static int
489find_evaluation(int maxtries, int sametries, int failtries, const CanonicalForm &f , const Variable & Extension, SFormList & BestEvaluationpoint, CFFList & BestFactorisation ){
490  int success;
491
492  success=evaluate( maxtries, sametries, failtries, f , Extension, BestEvaluationpoint, BestFactorisation );
493  return success;
494}
495#endif
496
497///////////////////////////////////////////////////////////////
498// A factorization routine for a sqrfree polynomial.         //
499// Returns the list of factors.                              //
500///////////////////////////////////////////////////////////////
501CFFList
502Factorized( const CanonicalForm & F, const Variable & alpha, int Mainvar){
503  CanonicalForm f,lt,ff,ffuni;
504  Variable Extension=alpha;
505  CFFList Outputlist,UnivariateFactorlist,Outputlist2;
506  SFormList Substitutionlist, Evaluationpoint;
507  CFFactor copy;
508  int mainvar=Mainvar,success,oldmainvar;
509  CFMap m;
510
511  // INTERRUPTHANDLER
512  if ( interrupt_handle() ) return CFFList() ;
513  // INTERRUPTHANDLER
514
515  if ( F.isUnivariate() ){ // could have lost one Variable elsewhere
516    if ( degree(Extension) == 0 ){
517      TIMING_START(evaluate_time);
518      Outputlist = factorize(F,1); // poly is sqr-free
519      TIMING_END(evaluate_time);
520      return Outputlist;
521    }
522    else{
523      TIMING_START(evaluate_time);
524      Outputlist = factorize(F,Extension);
525      TIMING_END(evaluate_time);
526      return Outputlist;
527    }
528  }
529
530  if ( Mainvar ) oldmainvar=Mainvar; else oldmainvar=level(F);
531  // First choose a main variable; this may be revisted in the next step
532  mainvar = choose_main_variable(F);
533  // Let`s look if @f/@mainvar is nonzero
534  mainvar = necessary_condition(F,mainvar);
535  // Now we have definetly choosen a main variable
536  // swap poly such that the mainvar has highest level
537  f=swapvar(F,mainvar,level(F));
538 
539  // INTERRUPTHANDLER
540  if ( interrupt_handle() ) return CFFList() ;
541  // INTERRUPTHANDLER
542
543  if ( oldmainvar != mainvar ){
544    DEBOUTSL(cout); DEBOUT(cout,"Swapped poly ", F);
545    DEBOUT(cout, " in ", f); DEBOUTNL(cout);
546    DEBOUTSL(cout); DEBOUT(cout,"Swapped  ", oldmainvar );
547    DEBOUT(cout, " <-- ", mainvar ); DEBOUT(cout, "  Mainvar= ", f.mvar());
548    DEBOUTNL(cout);
549    ff = f.deriv();
550    TIMING_START(discr_time);
551    ffuni = mygcd(f,ff);
552    TIMING_END(discr_time);
553    if ( ffuni != 1 ){ //discriminante nonzero: split poly
554      DEBOUTLN(cout,"Nontrivial GCD of f= ", f); 
555      DEBOUTLN(cout,"             and @f= ", ff);
556      DEBOUTLN(cout,"          GCD(f,@f)= ", ffuni);
557      ff=f/ffuni;
558      CFFList Outputlist_a, Outputlist_b;
559      Outputlist_a = Factorized(ff,alpha);
560      DEBOUTLN(cout, "Outputlist_a = ", Outputlist_a);
561      Outputlist_b = Factorized(ffuni,alpha);
562      DEBOUTLN(cout, "Outputlist_b = ", Outputlist_b);
563      Outputlist = myUnion(Outputlist_a, Outputlist_b);
564      // have to back-swapvar the factors....
565      for ( CFFListIterator i=Outputlist; i.hasItem(); i++ ){
566        copy=i.getItem();
567        Outputlist2.append(CFFactor(swapvar(copy.factor(),oldmainvar,mainvar),copy.exp()));
568      }
569      DEBOUTLN(cout, "Outputlist2 (a+b swapped) (to return) = ", Outputlist2);
570      return Outputlist2;
571    }
572  }
573
574  // Check special cases
575  for ( int i=1; i<=level(F); i++)
576    if ( degree(f,Variable(i) ) == 1 ) { //test trivial case; only true iff F is primitiv w.r.t every variable; else check (if F=ax+b) gcd(a,b)=1 ?
577      DEBOUTLN(cout, "Trivial case: ", F);
578      Outputlist.append(CFFactor(F,1));
579      return Outputlist;
580    }
581
582  // Look at the leading term:
583  lt = LC(f);
584  DEBOUTLN(cout, "Leading term: ", lt);
585  if ( lt != f.genOne() ){
586    // make the polynomial monic in the main variable
587    ff = make_monic(f,lt); ffuni = ff;
588    DEBOUTLN(cout, "make_monic returned: ", ff);
589  }
590  else{ ff= f; ffuni= ff; }
591
592  TIMING_START(evaluate_time);
593  success=evaluate(min(10,max(degree(ff), 5)), min(degree(ff),3), min(degree(ff),3), ff, Extension, Substitutionlist,UnivariateFactorlist);
594  DEBOUTLN(cout,  "Returned from evaluate: success: ", success);
595  for ( SFormListIterator ii=Substitutionlist; ii.hasItem(); ii++ ){
596    DEBOUTLN(cout, "Substituting ", ii.getItem().factor());
597    DEBOUTLN(cout, "       with value: ", ii.getItem().exp());
598  }
599
600  if ( success==0 ){ // evalute wasn't successfull
601    success= specializePoly(ffuni,Extension,degree(ff),Substitutionlist,1,getNumVars(compress(ff,m)));
602    DEBOUTLN(cout,  "Returned from specializePoly: success: ", success);
603    if (success == 0 ){ // No spezialisation could be found
604#ifdef HAVE_SINGULAR
605      WerrorS("libfac: Factorize: ERROR: Not able to find a valid specialization!");   
606#else
607      cerr << "libfac: Factorize: ERROR: Not able to find a valid specialization!\n" 
608           << rcsid << errmsg << endl; 
609#endif
610      Outputlist.append(CFFactor(F,1));
611      return Outputlist;
612    }
613
614    // INTERRUPTHANDLER
615    if ( interrupt_handle() ) return CFFList() ;
616    // INTERRUPTHANDLER
617
618    ffuni = substitutePoly(ff,Substitutionlist);
619    // We now have an univariat poly; factorize that
620    if ( degree(Extension) == 0   ){
621      DEBOUTMSG(cout, "Univ. Factorization over the ground field");
622      UnivariateFactorlist = factorize(ffuni,1); // univ. poly is sqr-free!
623    }
624    else{
625      DEBOUTLN(cout, "Univ. Factorization over extension of degree ",
626               degree(getMipo(Extension,'x')) );
627      UnivariateFactorlist = factorize(ffuni,Extension);
628    }
629  }
630  else{
631    ffuni = substitutePoly(ff,Substitutionlist); 
632  }
633    TIMING_END(evaluate_time);
634  if (UnivariateFactorlist.length() == 1){ // poly is irreduzibel
635    DEBOUTLN(cout, "Univ. poly is irreduzible: ", UnivariateFactorlist);
636    Outputlist.append(CFFactor(F,1));
637    return Outputlist;
638  }
639  else{ // we have factors
640    DEBOUTSL(cout); 
641    DEBOUT(cout, "Univariate poly has " , UnivariateFactorlist.length());
642    DEBOUT(cout, " factors:  ", ffuni); 
643    DEBOUT(cout, " = ", UnivariateFactorlist); DEBOUTNL(cout);
644
645    // INTERRUPTHANDLER
646    if ( interrupt_handle() ) return CFFList() ;
647    // INTERRUPTHANDLER
648
649    TIMING_START(hensel_time);
650    Outputlist = MultiHensel(ff,UnivariateFactorlist,Substitutionlist);
651    DEBOUTLN(cout, "Outputlist after MultiHensel: ", Outputlist);
652    TIMING_END(hensel_time);
653
654    // INTERRUPTHANDLER
655    if ( interrupt_handle() ) return CFFList() ;
656    // INTERRUPTHANDLER
657
658    TIMING_START(truefactor_time);
659    Outputlist = Truefactors(ff, level(ff), Substitutionlist, Outputlist);
660    DEBOUTLN(cout, "Outputlist after Truefactors: ", Outputlist);
661    TIMING_END(truefactor_time);
662
663    // INTERRUPTHANDLER
664    if ( interrupt_handle() ) return CFFList() ;
665    // INTERRUPTHANDLER
666
667    if ( lt != f.genOne() ){
668      Outputlist = not_monic(Outputlist,lt,ff,level(ff));
669      DEBOUTLN(cout, "not_monic returned: ", Outputlist);
670    }
671
672    // have to back-swapvar the factors....
673    for ( CFFListIterator i=Outputlist; i.hasItem(); i++ ){
674        copy=i.getItem();
675        Outputlist2.append(CFFactor(swapvar(copy.factor(),oldmainvar,mainvar),copy.exp()));
676    }
677
678    return Outputlist2;
679  }
680}
681
682///////////////////////////////////////////////////////////////
683// The user front-end for a uni/multivariate factorization   //
684// routine. F needs not to be SqrFree.                       //
685// Option of * choosing a  main variable (n.y.i.)            //
686//           * choosing an algebraic extension (n.y.u.)      //
687//           * ensuring poly is sqrfree (n.y.i.)             //
688///////////////////////////////////////////////////////////////
689CFFList
690Factorize( const CanonicalForm & F, int is_SqrFree ){
691  CFFList Outputlist,SqrFreeList,Intermediatelist,Outputlist2;
692  ListIterator<CFFactor> i,j;
693  CanonicalForm g=1,unit=1,r=1; 
694  Variable minpoly; // reserved (-> Factorisation over algebraic Extensions)
695  int exp;
696  CFMap m;
697
698  // INTERRUPTHANDLER
699  if ( interrupt_handle() ) return CFFList() ;
700  // INTERRUPTHANDLER
701
702  DEBINCLEVEL(cout, "Factorize");
703  DEBOUTMSG(cout, rcsid);
704  if ( getCharacteristic() == 0 ) { // char == 0
705    TIMING_START(factorize_time);
706    //cout << "Factoring in char=0 of " << F << " = " << Outputlist << endl;
707    Outputlist= factorize(F);
708    // Factorization in char=0 doesn't sometimes return at least two elements!!!
709    if ( getNumVars(Outputlist.getFirst().factor()) != 0 ) 
710      Outputlist.insert(CFFactor(1,1));
711    //cout << "  Factorize in char=0: returning with: " << Outputlist << endl;
712    TIMING_END(factorize_time);
713    DEBDECLEVEL(cout, "Factorize");
714    TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
715    return Outputlist;
716  }
717  TIMING_START(factorize_time);
718  ///////
719  // Maybe it`s better to add a sqrfree-test before?
720  // (If gcd is fast...)
721  ///////
722  //  if ( ! SqrFreeTest(F) ){
723  if ( ! is_SqrFree ){
724    TIMING_START(sqrfree_time);
725    SqrFreeList = InternalSqrFree(F) ; // first sqrfree the polynomial
726    // don't use sqrFree(F), factory's internal sqrFree for multiv.
727    // Polynomials; it's wrong!! Ex.: char=p   f= x^p*(y+1);
728    // InternalSqrFree(f)= ( y+1, (x)^p ), sqrFree(f)= ( y+1 ) .
729    TIMING_END(sqrfree_time);
730
731    // INTERRUPTHANDLER
732    if ( interrupt_handle() ) return CFFList() ;
733    // INTERRUPTHANDLER
734
735  }
736  else 
737    SqrFreeList.append(CFFactor(F,1));
738  DEBOUTLN(cout, "InternalSqrFreeList= ", SqrFreeList);
739  for ( i=SqrFreeList; i.hasItem(); i++ ){
740    DEBOUTLN(cout, "Factor under consideration: ", i.getItem().factor());
741    // We need a compress on each list item ! Maybe we have less variables!
742    g =compress(i.getItem().factor(),m); 
743    exp = i.getItem().exp();
744    if ( getNumVars(g) ==0 ) // a constant; Exp==1
745      Outputlist.append( CFFactor(g,1) ) ;
746    else// a real polynomial
747      if ( g.isUnivariate() ){
748        Intermediatelist=factorize(g,1); // poly is sqr-free!
749        for ( j=Intermediatelist; j.hasItem(); j++ )
750          //Normally j.getItem().exp() should be 1
751          Outputlist.append( CFFactor( m(j.getItem().factor()),exp*j.getItem().exp()));
752      }
753      else{ // multivariate polynomial
754        if ( is_homogeneous(g) ){
755          DEBOUTLN(cout, "Poly is homogeneous! : ", g);
756          // Now we can substitute one variable to 1, factorize and then
757          // look on the resulting factors and their monomials for
758          // backsubstitution of the substituted variable.
759          Intermediatelist = HomogFactor(g, minpoly, 0);
760        }
761        else // not homogeneous
762          Intermediatelist = Factorized(g, minpoly, 0);
763
764        // INTERRUPTHANDLER
765        if ( interrupt_handle() ) return CFFList() ;
766        // INTERRUPTHANDLER
767
768        for ( j=Intermediatelist; j.hasItem(); j++ )
769          //Normally j.getItem().exp() should be 1
770          Outputlist= myappend( Outputlist, CFFactor(m(j.getItem().factor()),exp*j.getItem().exp()));
771      }
772  }
773  g=1; unit=1;
774  DEBOUTLN(cout, "Outputlist is ", Outputlist);
775  for ( i=Outputlist; i.hasItem(); i++ )
776    if ( level(i.getItem().factor()) > 0 ){
777      unit = lc(i.getItem().factor());
778      if ( getNumVars(unit) == 0 ){ // a constant; possibly 1
779        Outputlist2.append(CFFactor(i.getItem().factor()/unit , i.getItem().exp()));
780        g *=power(i.getItem().factor()/unit,i.getItem().exp());
781      }
782      else{
783        Outputlist2.append(i.getItem());
784        g *=power(i.getItem().factor(),i.getItem().exp());
785      }
786    }
787 
788  r=F/g; 
789  Outputlist2.insert(CFFactor(r,1));
790 
791  DEBDECLEVEL(cout, "Factorize");
792  TIMING_END(factorize_time);
793
794  TIMING_PRINT(sqrfree_time, "\ntime used for sqrfree   : ");
795  TIMING_PRINT(discr_time, "time used for discriminante   : ");
796  TIMING_PRINT(evaluate_time, "time used for evaluation and univ. factorization  : ");
797  TIMING_PRINT(hensel_time, "time used for hensel-lift   : ");
798  TIMING_PRINT(truefactor_time, "time used for truefactors   : ");
799  TIMING_PRINT(factorize_time, "\ntime used for factorization   : ");
800  return Outputlist2;
801}
802
803/*
804$Log: not supported by cvs2svn $
805Revision 1.3  1997/09/12 07:19:46  Singular
806* hannes/michael: libfac-0.3.0
807
808Revision 1.6  1997/04/25 22:18:40  michael
809changed lc == 1 to lc == unit in choose_mainvar
810changed cerr and cout messages for use with Singular
811Version for libfac-0.2.1
812
813Revision 1.5  1997/01/17 05:04:03  michael
814* added support for homogenous polynomials
815* exported functions to homogfactor.cc
816
817*/
Note: See TracBrowser for help on using the repository browser.