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

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