source: git/libfac/factor/Factor.cc @ 3e55bc

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