source: git/libfac/factor/Factor.cc @ 14b1e65

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