source: git/libfac/factor/Factor.cc @ b87513c

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