source: git/libfac/charset/alg_factor.cc @ 6f801b

spielwiese
Last change on this file since 6f801b was 6f801b, checked in by Hans Schönemann <hannes@…>, 16 years ago
hannes: genZero ->isZero git-svn-id: file:///usr/local/Singular/svn/trunk@10594 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 27.9 KB
Line 
1/* Copyright 1997 Michael Messollen. All rights reserved. */
2////////////////////////////////////////////////////////////
3// emacs edit mode for this file is -*- C++ -*-
4////////////////////////////////////////////////////////////
5static char * rcsid = "$Id: alg_factor.cc,v 1.20 2008-02-22 12:16:02 Singular Exp $";
6////////////////////////////////////////////////////////////
7// FACTORY - Includes
8#include <factory.h>
9// Factor - Includes
10#include <tmpl_inst.h>
11#include <Factor.h>
12#include <SqrFree.h>
13#include <helpstuff.h>
14// Charset - Includes
15#include "csutil.h"
16#include "charset.h"
17#include "reorder.h"
18#include "algfactor.h"
19// some CC's need this:
20#include "alg_factor.h"
21
22void out_cf(char *s1,const CanonicalForm &f,char *s2);
23
24#ifdef ALGFACTORDEBUG
25#  define DEBUGOUTPUT
26#else
27#  undef DEBUGOUTPUT
28#endif
29
30#include "debug.h"
31#include "timing.h"
32TIMING_DEFINE_PRINT(newfactoras_time);
33
34static Varlist
35Var_is_in_AS(const Varlist & uord, const CFList & Astar);
36
37int getAlgVar(const CanonicalForm &f, Variable &X)
38{
39  if (f.inBaseDomain()) return 0;
40  if (f.inCoeffDomain())
41  {
42    if (f.level()!=0)
43    {
44      X= f.mvar();
45      return 1;
46    }
47    return getAlgVar(f.LC(),X);
48  }
49  if (f.inPolyDomain())
50  {
51    if (getAlgVar(f.LC(),X)) return 1;
52    for( CFIterator i=f; i.hasTerms(); i++)
53    {
54      if (getAlgVar(i.coeff(),X)) return 1;
55    }
56  }
57  return 0;
58}
59
60////////////////////////////////////////////////////////////////////////
61// This implements the algorithm of Trager for factorization of
62// (multivariate) polynomials over algebraic extensions and so called
63// function field extensions.
64////////////////////////////////////////////////////////////////////////
65
66// // missing class: IntGenerator:
67bool IntGenerator::hasItems() const
68{
69    return 1;
70}
71
72CanonicalForm IntGenerator::item() const
73//int IntGenerator::item() const
74{
75  //return current; //CanonicalForm( current );
76  return mapinto(CanonicalForm( current ));
77}
78
79void IntGenerator::next()
80{
81    current++;
82}
83
84// replacement for factory's broken psr
85static CanonicalForm
86mypsr ( const CanonicalForm &rr, const CanonicalForm &vv, const Variable & x ){
87  CanonicalForm r=rr, v=vv, l, test, lu, lv, t, retvalue;
88  int dr, dv, d,n=0;
89
90
91  dr = degree( r, x );
92  dv = degree( v, x );
93  if (dv <= dr) {l=LC(v,x); v = v -l*power(x,dv);}
94  else { l = 1; }
95  d= dr-dv+1;
96  while ( ( dv <= dr  ) && ( !r.isZero()) )
97  {
98    test = power(x,dr-dv)*v*LC(r,x);
99    if ( dr == 0 ) { r= CanonicalForm(0); }
100    else { r= r - LC(r,x)*power(x,dr); }
101    r= l*r -test;
102    dr= degree(r,x);
103    n+=1;
104  }
105  r= power(l, d-n)*r;
106  return r;
107}
108
109// replacement for factory's broken resultant
110static CanonicalForm
111resultante( const CanonicalForm & f, const CanonicalForm& g, const Variable & v ){
112  bool on_rational = isOn(SW_RATIONAL);
113  On(SW_RATIONAL);
114  CanonicalForm cd = bCommonDen( f );
115  CanonicalForm fz = f * cd;
116  cd = bCommonDen( g );
117  CanonicalForm gz = g * cd;
118  if (!on_rational)  Off(SW_RATIONAL);
119
120return resultant(fz,gz,v);
121
122  CanonicalForm h, beta, help, F, G;
123  int delta;
124
125  DEBOUTLN( CERR, "resultante: called  f= ", f);
126  DEBOUTLN( CERR, "resultante: called  g= ", g);
127  DEBOUTLN( CERR, "resultante: called  v= ", v);
128  if ( f.mvar() < v || g.mvar() < v ){
129    DEBOUTMSG(CERR, "resultante: f.mvar() < v || g.mvar() < v");
130    return 1;
131  }
132
133  if ( f.degree( v ) < 1 || g.degree( v ) < 1 ){
134    DEBOUTMSG(CERR, "resultante: f.degree( v ) < 1 || g.degree( v ) < 1");
135    // If deg(F,v) == 0 , then resultante(F,G,v) = F^n, where n=deg(G,v)
136    if ( f.degree( v ) < 1 ) return power(f,degree(g,v));
137    else return power(g,degree(f,v));
138  }
139
140   if ( f.degree( v ) >= g.degree( v ) ) { F = f; G = g; }
141   else { G = f; F = g; }
142
143  h = CanonicalForm(1);
144  while ( !G.isZero() )
145  {
146     delta= degree(F,v) -degree(G,v);
147     beta = power(CanonicalForm(-1), delta+1) * LC(F,v)* power(h, delta);
148     h= (h * power(LC(G,v), delta)) / power(h, delta);
149     help= G;
150     G= mypsr(F,G,v);
151     G= G/beta;
152     F=help;
153   }
154   if ( degree(F,v) != 0 ) F= CanonicalForm(0);
155   return F;
156}
157
158// sqr-free routine for algebraic extensions
159// we need it! Ex.: f=c^2+2*a*c-1; as=[a^2+1]; f=(c+a)^2
160static CFFList
161alg_sqrfree( const CanonicalForm & f ){
162  CFFList L;
163
164  L.append(CFFactor(f,1));
165  return L;
166}
167
168// Calculates a square free norm
169// Input: f(x, alpha) a square free polynomial over K(alpha),
170// alpha is defined by the minimal polynomial Palpha
171// K has more than S elements (S is defined in thesis; look getextension)
172static void
173sqrf_norm_sub( const CanonicalForm & f, const CanonicalForm & PPalpha,
174           CFGenerator & myrandom, CanonicalForm & s,  CanonicalForm & g,
175           CanonicalForm & R){
176  Variable y=PPalpha.mvar(),vf=f.mvar();
177  CanonicalForm temp, Palpha=PPalpha, t;
178  int sqfreetest=0;
179  CFFList testlist;
180  CFFListIterator i;
181
182  DEBOUTLN(CERR, "sqrf_norm_sub:      f= ", f);
183  DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
184  myrandom.reset();   s=f.mvar()-myrandom.item()*Palpha.mvar();   g=f;
185  R= CanonicalForm(0);
186  DEBOUTLN(CERR, "sqrf_norm_sub: myrandom s= ", s);
187
188  // Norm, resultante taken with respect to y
189  while ( !sqfreetest ){
190    DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
191    R = resultante(Palpha, g, y); R= R* bCommonDen(R);
192    DEBOUTLN(CERR, "sqrf_norm_sub: R= ", R);
193    // sqfree check ; R is a polynomial in K[x]
194    if ( getCharacteristic() == 0 )
195    {
196      temp= gcd(R, R.deriv(vf));
197      DEBOUTLN(CERR, "sqrf_norm_sub: temp= ", temp);
198      if (degree(temp,vf) != 0 || temp.isZero() ){ sqfreetest= 0; }
199      else { sqfreetest= 1; }
200      DEBOUTLN(CERR, "sqrf_norm_sub: sqfreetest= ", sqfreetest);
201    }
202    else{
203      DEBOUTMSG(CERR, "Starting isSqrFree(R)!");
204      // Look at isSqrFree!
205      // (z+a^5+w)^4 with z<w<a should not give sqfreetest=1 !
206      // for now we use this workaround with Factorize...
207      // ...but it should go away soon!!!!
208      Variable X;
209      if (getAlgVar(R,X))
210      {
211        if (R.isUnivariate())
212          testlist=factorize( R, X );
213        else
214          testlist= Factorize(R, X, 0);
215      }
216      else
217        testlist= Factorize(R);
218      DEBOUTLN(CERR, "testlist= ", testlist);
219      testlist.removeFirst();
220      sqfreetest=1;
221      for ( i=testlist; i.hasItem(); i++)
222        if ( i.getItem().exp() > 1 && degree(i.getItem().factor(), R.mvar()) > 0) { sqfreetest=0; break; }
223      DEBOUTLN(CERR, "isSqrFree(R)= ", sqfreetest);
224    }
225    if ( ! sqfreetest ){
226      myrandom.next();
227      DEBOUTLN(CERR, "sqrf_norm_sub generated new myrandom item: ", myrandom.item());
228      if ( getCharacteristic() == 0 ) t= CanonicalForm(mapinto(myrandom.item()));
229      else t= CanonicalForm(myrandom.item());
230      s= f.mvar()+t*Palpha.mvar(); // s defines backsubstitution
231      DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);
232      g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
233      DEBOUTLN(CERR, "             gives g= ", g);
234    }
235  }
236}
237static void
238sqrf_agnorm_sub( const CanonicalForm & f, const CanonicalForm & PPalpha,
239           AlgExtGenerator & myrandom, CanonicalForm & s,  CanonicalForm & g,
240           CanonicalForm & R){
241  Variable y=PPalpha.mvar(),vf=f.mvar();
242  CanonicalForm temp, Palpha=PPalpha, t;
243  int sqfreetest=0;
244  CFFList testlist;
245  CFFListIterator i;
246
247  DEBOUTLN(CERR, "sqrf_norm_sub:      f= ", f);
248  DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
249  myrandom.reset();   s=f.mvar()-myrandom.item()*Palpha.mvar();   g=f;
250  R= CanonicalForm(0);
251  DEBOUTLN(CERR, "sqrf_norm_sub: myrandom s= ", s);
252
253  // Norm, resultante taken with respect to y
254  while ( !sqfreetest ){
255    DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
256    R = resultante(Palpha, g, y); R= R* bCommonDen(R);
257    DEBOUTLN(CERR, "sqrf_norm_sub: R= ", R);
258    // sqfree check ; R is a polynomial in K[x]
259    if ( getCharacteristic() == 0 )
260    {
261      temp= gcd(R, R.deriv(vf));
262      DEBOUTLN(CERR, "sqrf_norm_sub: temp= ", temp);
263      if (degree(temp,vf) != 0 || temp.isZero() ){ sqfreetest= 0; }
264      else { sqfreetest= 1; }
265      DEBOUTLN(CERR, "sqrf_norm_sub: sqfreetest= ", sqfreetest);
266    }
267    else{
268      DEBOUTMSG(CERR, "Starting isSqrFree(R)!");
269      // Look at isSqrFree!
270      // (z+a^5+w)^4 with z<w<a should not give sqfreetest=1 !
271      // for now we use this workaround with Factorize...
272      // ...but it should go away soon!!!!
273      Variable X;
274      if (getAlgVar(R,X))
275      {
276        if (R.isUnivariate())
277          testlist=factorize( R, X );
278        else
279          testlist= Factorize(R, X, 0);
280      }
281      else
282        testlist= Factorize(R);
283      DEBOUTLN(CERR, "testlist= ", testlist);
284      testlist.removeFirst();
285      sqfreetest=1;
286      for ( i=testlist; i.hasItem(); i++)
287        if ( i.getItem().exp() > 1 && degree(i.getItem().factor(), R.mvar()) > 0) { sqfreetest=0; break; }
288      DEBOUTLN(CERR, "isSqrFree(R)= ", sqfreetest);
289    }
290    if ( ! sqfreetest ){
291      myrandom.next();
292      DEBOUTLN(CERR, "sqrf_norm_sub generated new myrandom item: ", myrandom.item());
293      if ( getCharacteristic() == 0 ) t= CanonicalForm(mapinto(myrandom.item()));
294      else t= CanonicalForm(myrandom.item());
295      s= f.mvar()+t*Palpha.mvar(); // s defines backsubstitution
296      DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);
297      g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
298      DEBOUTLN(CERR, "             gives g= ", g);
299    }
300  }
301}
302
303static void
304sqrf_norm( const CanonicalForm & f, const CanonicalForm & PPalpha,
305           const Variable & Extension, CanonicalForm & s,  CanonicalForm & g,
306           CanonicalForm & R){
307
308  DEBOUTLN(CERR, "sqrf_norm:      f= ", f);
309  DEBOUTLN(CERR, "sqrf_norm: Palpha= ", PPalpha);
310  if ( getCharacteristic() == 0 ) {
311    IntGenerator myrandom;
312    DEBOUTMSG(CERR, "sqrf_norm: no extension, char=0");
313    sqrf_norm_sub(f,PPalpha, myrandom, s,g,R);
314    DEBOUTLN(CERR, "sqrf_norm:      f= ", f);
315    DEBOUTLN(CERR, "sqrf_norm: Palpha= ", PPalpha);
316    DEBOUTLN(CERR, "sqrf_norm:      s= ", s);
317    DEBOUTLN(CERR, "sqrf_norm:      g= ", g);
318    DEBOUTLN(CERR, "sqrf_norm:      R= ", R);
319  }
320  else if ( degree(Extension) > 0 ){ // working over Extensions
321    DEBOUTLN(CERR, "sqrf_norm: degree of extension is ", degree(Extension));
322    AlgExtGenerator myrandom(Extension);
323    sqrf_agnorm_sub(f,PPalpha, myrandom, s,g,R);
324  }
325  else{
326    FFGenerator myrandom;
327    DEBOUTMSG(CERR, "sqrf_norm: degree of extension is 0");
328    sqrf_norm_sub(f,PPalpha, myrandom, s,g,R);
329  }
330}
331
332static Varlist
333Var_is_in_AS(const Varlist & uord, const CFList & Astar){
334  Varlist output;
335  CanonicalForm elem;
336  Variable x;
337
338  for ( VarlistIterator i=uord; i.hasItem(); i++){
339    x=i.getItem();
340    for ( CFListIterator j=Astar; j.hasItem(); j++ ){
341      elem= j.getItem();
342      if ( degree(elem,x) > 0 ){ // x actually occures in Astar
343        output.append(x);
344        break;
345      }
346    }
347  }
348  return output;
349}
350
351// Look if Minimalpolynomials in Astar define seperable Extensions
352// Must be a power of p: i.e. y^{p^e}-x
353static int
354inseperable(const CFList & Astar){
355  CanonicalForm elem;
356  int Counter= 1;
357
358  if ( Astar.length() == 0 ) return 0;
359  for ( CFListIterator i=Astar; i.hasItem(); i++)
360  {
361    elem= i.getItem();
362    if ( elem.deriv().isZero() ) return Counter;
363    else Counter += 1;
364  }
365  return 0;
366}
367
368// calculate gcd of f and g in char=0
369static CanonicalForm
370gcd0( CanonicalForm f, CanonicalForm g ){
371  int charac= getCharacteristic();
372  int save=isOn(SW_RATIONAL);
373  setCharacteristic(0); Off(SW_RATIONAL);
374  CanonicalForm ff= mapinto(f), gg= mapinto(g);
375  CanonicalForm result= gcd(ff,gg);
376  setCharacteristic(charac); 
377  if (save) On(SW_RATIONAL);
378  return mapinto(result);
379}
380
381// calculate big enough extension for finite fields
382// Idea: first calculate k, such that q^k > S (->thesis, -> getextension)
383// Second, search k with gcd(k,m_i)=1, where m_i is the degree of the i'th
384// minimal polynomial. Then the minpoly f_i remains irrd. over q^k and we
385// have enough elements to plug in.
386static int
387getextension( IntList & degreelist, int n){
388  int charac= getCharacteristic();
389  setCharacteristic(0); // need it for k !
390  int k=1, m=1, length=degreelist.length();
391  IntListIterator i;
392
393  for (i=degreelist; i.hasItem(); i++) m= m*i.getItem();
394  int q=charac;
395  while (q <= ((n*m)*(n*m)/2)) { k= k+1; q= q*charac;}
396  int l=0;
397  do {
398    for (i=degreelist; i.hasItem(); i++){
399      l= l+1;
400      if ( gcd0(k,i.getItem()) == 1){
401        DEBOUTLN(CERR, "getextension: gcd == 1, l=",l);
402        if ( l==length ){ setCharacteristic(charac);  return k; }
403      }
404      else { DEBOUTMSG(CERR, "getextension: Next iteration"); break; }
405    }
406    k= k+1; l=0;
407  }
408  while ( 1 );
409}
410
411// calculate a "primitive element"
412// K must have more than S elements (->thesis, -> getextension)
413static CFList
414simpleextension(const CFList & Astar, const Variable & Extension,
415                CanonicalForm & R){
416  CFList Returnlist, Bstar=Astar;
417  CanonicalForm s, g;
418
419  DEBOUTLN(CERR, "simpleextension: Astar= ", Astar);
420  DEBOUTLN(CERR, "simpleextension:     R= ", R);
421  DEBOUTLN(CERR, "simpleextension: Extension= ", Extension);
422  if ( Astar.length() == 1 ){ R= Astar.getFirst();}
423  else{
424    R=Bstar.getFirst(); Bstar.removeFirst();
425    for ( CFListIterator i=Bstar; i.hasItem(); i++){
426      DEBOUTLN(CERR, "simpleextension: f(x)= ", i.getItem());
427      DEBOUTLN(CERR, "simpleextension: P(x)= ", R);
428      sqrf_norm(i.getItem(), R, Extension, s, g, R);
429      // spielt die Repraesentation eine Rolle?
430      // muessen wir die Nachfolger aendern, wenn s != 0 ?
431      DEBOUTLN(CERR, "simpleextension: g= ", g);
432      if ( s != 0 ) DEBOUTLN(CERR, "simpleextension: s= ", s);
433      else DEBOUTLN(CERR, "simpleextension: s= ", s);
434      DEBOUTLN(CERR, "simpleextension: R= ", R);
435      Returnlist.insert(s);
436    }
437  }
438
439  return Returnlist;
440}
441
442CanonicalForm alg_lc(const CanonicalForm &f)
443{
444  if (f.inCoeffDomain()) return f;
445  if (f.level()>0)
446  {
447    return alg_lc(f.LC());
448  }
449}
450
451// the heart of the algorithm: the one from Trager
452static CFFList
453alg_factor( const CanonicalForm & f, const CFList & Astar, const Variable & vminpoly, const Varlist & oldord, const CFList & as)
454{
455  CFFList L, Factorlist;
456  CanonicalForm R, Rstar, s, g, h;
457  CFList substlist;
458
459  DEBINCLEVEL(CERR,"alg_factor");
460  DEBOUTLN(CERR, "alg_factor: f= ", f);
461
462  //out_cf("start alg_factor:",f,"\n");
463  substlist= simpleextension(Astar, vminpoly, Rstar);
464  DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);
465  DEBOUTLN(CERR, "alg_factor: minpoly Rstar= ", Rstar);
466  DEBOUTLN(CERR, "alg_factor: vminpoly= ", vminpoly);
467
468  sqrf_norm(f, Rstar, vminpoly, s, g, R );
469  //out_cf("sqrf_norm R:",R,"\n");
470  //out_cf("sqrf_norm s:",s,"\n");
471  //out_cf("sqrf_norm g:",g,"\n");
472  DEBOUTLN(CERR, "alg_factor: g= ", g);
473  DEBOUTLN(CERR, "alg_factor: s= ", s);
474  DEBOUTLN(CERR, "alg_factor: R= ", R);
475  Off(SW_RATIONAL);
476  Variable X;
477  if (getAlgVar(R,X))
478  {
479    // factorize R over alg.extension with X
480//CERR << "alg: "<< X << " mipo=" << getMipo(X,Variable('X')) <<"\n";
481    if (R.isUnivariate())
482    {
483      DEBOUTLN(CERR, "alg_factor: factorize( ", R);
484      Factorlist =  factorize( R, X );
485    }
486    else
487    {
488      #if 1
489      Variable XX;
490      CanonicalForm mipo=getMipo(X,XX);
491      CFList as(mipo);
492      DEBOUTLN(CERR, "alg_factor: newfactoras( ", R);
493      Factorlist = newfactoras(R, as , 1);
494      #else
495      // factor R over k
496      DEBOUTLN(CERR, "alg_factor: Factorize( ", R);
497      Factorlist = Factorize(R);
498      #endif
499    }
500  }
501  else
502  {
503    // factor R over k
504    DEBOUTLN(CERR, "alg_factor: Factorize( ", R);
505    Factorlist = Factorize(R);
506  }
507  On(SW_RATIONAL);
508  DEBOUTLN(CERR, "alg_factor: Factorize(R)= ", Factorlist);
509  if ( !Factorlist.getFirst().factor().inCoeffDomain() )
510    Factorlist.insert(CFFactor(1,1));
511  if ( Factorlist.length() == 2 && Factorlist.getLast().exp()== 1)
512  { // irreduzibel (first entry is a constant)
513    L.append(CFFactor(f,1));
514  }
515  else
516  {
517    DEBOUTLN(CERR, "alg_factor: g= ", g);
518    CanonicalForm gnew= g(s,s.mvar());
519    DEBOUTLN(CERR, "alg_factor: gnew= ", gnew);
520    g=gnew;
521    for ( CFFListIterator i=Factorlist; i.hasItem(); i++)
522    {
523      CanonicalForm fnew=i.getItem().factor();
524      fnew= fnew(s,s.mvar());
525      DEBOUTLN(CERR, "alg_factor: fnew= ", fnew);
526      DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);
527      for ( CFListIterator ii=substlist; ii.hasItem(); ii++){
528        DEBOUTLN(CERR, "alg_factor: item= ", ii.getItem());
529        fnew= fnew(ii.getItem(), ii.getItem().mvar());
530        DEBOUTLN(CERR, "alg_factor: fnew= ", fnew);
531      }
532      if (degree(i.getItem().factor()) > 0 )
533      {
534        // undo linear transformation!!!! and then gcd!
535        //CERR << "algcd(" << g << "," << fnew << ",as" << as << ")" << "\n";
536        //out_cf("g=",g,"\n");
537        //out_cf("fnew=",fnew,"\n");
538        //h= algcd(g,fnew, as, oldord);
539        //if (as.length() >1)
540        //  h= algcd(g,fnew, as, oldord);
541        //else
542          h=alg_gcd(g,fnew,as);
543        //out_cf(" -> algcd=",algcd(g,fnew, as, oldord),"\n");
544        //out_cf(" -> alg_gcd=",alg_gcd(g,fnew,as),"\n");
545        //CERR << "algcd result:" << h << "\n";
546        DEBOUTLN(CERR, "  alg_factor: h= ", h);
547        DEBOUTLN(CERR, "  alg_factor: oldord= ", oldord);
548        if ( degree(h) > 0 )
549        { //otherwise it's a constant
550          //CanonicalForm c=LC(h,g.mvar());
551          //out_cf("new lc h:",c,"\n");
552          //h= divide(h,c,as);
553          //out_cf("new factor h/c:",h,"\n");
554          g= divide(g, h,as);
555          DEBOUTLN(CERR, "alg_factor: g/h= ", g);
556          DEBOUTLN(CERR, "alg_factor: s= ", s);
557          DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);
558          //out_cf("new g:",g,"\n");
559          L.append(CFFactor(h,1));
560        }
561        //else printf("degree <=1\n");
562      }
563    }
564    // we are not interested in a
565    // constant (over K_r, which can be a polynomial!)
566    if (degree(g, f.mvar())>0){ L.append(CFFactor(g,1)); }
567  }
568  CFFList LL;
569  if (getCharacteristic()>0)
570  {
571    CFFListIterator i=L;
572    CanonicalForm c_fac=1;
573    CanonicalForm c;
574    for(;i.hasItem(); i++ )
575    {
576      CanonicalForm ff=i.getItem().factor();
577      c=alg_lc(ff);
578      int e=i.getItem().exp();
579      ff/=c;
580      if (!ff.isOne()) LL.append(CFFactor(ff,e));
581      while (e>0) { c_fac*=c;e--; }
582    }
583    if (!c_fac.isOne()) LL.insert(CFFactor(c_fac,1));
584  }
585  else
586  {
587    LL=L;
588  }
589  //CFFListIterator i=LL;
590  //for(;i.hasItem(); i++ )
591  //  out_cf("end alg_f:",i.getItem().factor(),"\n");
592  //printf("end alg_factor\n");
593  DEBOUTLN(CERR, "alg_factor: L= ", LL);
594  DEBDECLEVEL(CERR,"alg_factor");
595  return LL;
596}
597
598static CFFList
599endler( const CanonicalForm & f, const CFList & AS, const Varlist & uord ){
600  CanonicalForm F=f, g, q,r;
601  CFFList Output;
602  CFList One, Two, asnew, as=AS;
603  CFListIterator i,ii;
604  VarlistIterator j;
605  Variable vg;
606
607  for (i=as; i.hasItem(); i++){
608    g= i.getItem();
609    if (g.deriv() == 0 ){
610      DEBOUTLN(CERR, "Inseperable extension detected: ", g);
611      for (j=uord; j.hasItem(); j++){
612        if ( degree(g,j.getItem()) > 0 ) vg= j.getItem();
613      }
614      // Now we have the highest transzendental in vg;
615      DEBOUTLN(CERR, "Transzendental is ", vg);
616      CanonicalForm gg=-1*g[0];
617      divrem(gg,vg,q,r); r= gg-q*vg;   gg= gg-r;
618      //DEBOUTLN(CERR, "q= ", q); DEBOUTLN(CERR, "r= ", r);
619      DEBOUTLN(CERR, "  that is ", gg);
620      DEBOUTLN(CERR, "  maps to ", g+gg);
621      One.insert(gg); Two.insert(g+gg);
622      // Now transform all remaining polys in as:
623      int x=0;
624      for (ii=i; ii.hasItem(); ii++){
625        if ( x != 0 ){
626          divrem(ii.getItem(), gg, q,r);
627//          CERR << ii.getItem() << " divided by " << gg << "\n";
628          DEBOUTLN(CERR, "q= ", q); DEBOUTLN(CERR, "r= ", r);
629          ii.append(ii.getItem()+q*g); ii.remove(1);
630          DEBOUTLN(CERR, "as= ", as);
631        }
632        x+= 1;
633      }
634      // Now transform F:
635      divrem(F, gg, q,r);
636      F= F+q*g;
637      DEBOUTLN(CERR, "new F= ", F);
638    }
639    else{ asnew.append(i.getItem());  }// just the identity
640  }
641  // factor F with minimal polys given in asnew:
642  DEBOUTLN(CERR, "Factor F=  ", F);
643  DEBOUTLN(CERR, "  with as= ", asnew);
644  int success=0;
645  CFFList factorlist= newcfactor(F,asnew, success);
646  DEBOUTLN(CERR, "  gives =  ", factorlist);
647  DEBOUTLN(CERR, "One= ", One);
648  DEBOUTLN(CERR, "Two= ", Two);
649
650  // Transform back:
651  for ( CFFListIterator k=factorlist; k.hasItem(); k++){
652    CanonicalForm factor= k.getItem().factor();
653    ii=One;
654    for (i=Two; i.hasItem(); i++){
655      DEBOUTLN(CERR, "Mapping ", i.getItem());
656      DEBOUTLN(CERR, "     to ", ii.getItem());
657      DEBOUTLN(CERR, "     in ", factor);
658      divrem(factor,i.getItem(),q,r); r=factor -q*i.getItem();
659      DEBOUTLN(CERR, "q= ", q); DEBOUTLN(CERR, "r= ", r);
660      factor= ii.getItem()*q +r; //
661      ii++;
662    }
663    Output.append(CFFactor(factor,k.getItem().exp()));
664  }
665
666  return Output;
667}
668
669
670// 1) prepares data
671// 2) for char=p we distinguish 3 cases:
672//           no transcendentals, seperable and inseperable extensions
673CFFList
674newfactoras( const CanonicalForm & f, const CFList & as, int success){
675  Variable vf=f.mvar();
676  CFListIterator i;
677  CFFListIterator jj;
678  CFList reduceresult;
679  CFFList result;
680
681  success=1;
682  DEBINCLEVEL(CERR, "newfactoras");
683  DEBOUTMSG(CERR, rcsid);
684  DEBOUTLN(CERR, "newfactoras called with f= ", f);
685  DEBOUTLN(CERR, "               content(f)= ", content(f));
686  DEBOUTLN(CERR, "                       as= ", as);
687  DEBOUTLN(CERR, "newfactoras: cls(vf)= ", cls(vf));
688  DEBOUTLN(CERR, "newfactoras: cls(as.getLast())= ", cls(as.getLast()));
689  DEBOUTLN(CERR, "newfactoras: degree(f,vf)= ", degree(f,vf));
690
691// F1: [Test trivial cases]
692// 1) first trivial cases:
693  if ( (cls(vf) <= cls(as.getLast())) ||  degree(f,vf)<=1 ){
694// ||( (as.length()==1) && (degree(f,vf)==3) && (degree(as.getFirst()==2)) )
695    DEBDECLEVEL(CERR,"newfactoras");
696    return CFFList(CFFactor(f,1));
697  }
698
699// 2) List of variables:
700// 2a) Setup list of those polys in AS having degree(AS[i], AS[i].mvar()) > 1
701// 2b) Setup variableordering
702  CFList Astar;
703  Variable x;
704  CanonicalForm elem;
705  Varlist ord, uord,oldord;
706  for ( int ii=1; ii< level(vf) ; ii++ ) { uord.append(Variable(ii));  }
707  oldord= uord; oldord.append(vf);
708
709  for ( i=as; i.hasItem(); i++ ){
710    elem= i.getItem();
711    x= elem.mvar();
712    if ( degree(elem,x) > 1){ // otherwise it's not an extension
713      Astar.append(elem);
714      ord.append(x);
715    }
716  }
717  uord= Difference(uord,ord);
718  DEBOUTLN(CERR, "Astar is: ", Astar);
719  DEBOUTLN(CERR, "ord is: ", ord);
720  DEBOUTLN(CERR, "uord is: ", uord);
721
722// 3) second trivial cases: we already prooved irr. of f over no extensions
723  if ( Astar.length() == 0 ){
724    DEBDECLEVEL(CERR,"newfactoras");
725    return CFFList(CFFactor(f,1));
726  }
727
728// 4) Try to obtain a partial factorization using prop2 and prop3
729//    Use with caution! We have to proof these propositions first!
730  // Not yet implemented
731
732// 5) Look if elements in uord actually occure in any of the minimal
733//    polynomials. If no element of uord occures in any of the minimal
734//   polynomials, we don't have transzendentals.
735  Varlist newuord=Var_is_in_AS(uord,Astar);
736  DEBOUTLN(CERR, "newuord is: ", newuord);
737
738  CFFList Factorlist;
739  Varlist gcdord= Union(ord,newuord); gcdord.append(f.mvar());
740  // This is for now. we need alg_sqrfree implemented!
741  //CERR << "algcd(" << f << "," << f.deriv() << " as:" << Astar <<"\n";
742  //CanonicalForm Fgcd= algcd(f,f.deriv(),Astar,gcdord);
743  CanonicalForm Fgcd;
744        //if (Astar.length() >1)
745        //  Fgcd= algcd(f,f.deriv(),Astar,gcdord);
746        //else
747          Fgcd= alg_gcd(f,f.deriv(),Astar);
748        //out_cf("algcd:",algcd(f,f.deriv(),Astar,gcdord),"\n");
749        //out_cf("alg_gcd:",alg_gcd(f,f.deriv(),Astar),"\n");
750 // CERR << "algcd result:"  << Fgcd << "\n";
751  if ( Fgcd == 0 ) DEBOUTMSG(CERR, "WARNING: p'th root ?");
752  if (( degree(Fgcd, f.mvar()) > 0) && (!(f.deriv().isZero())) ){
753    DEBOUTLN(CERR, "Nontrivial GCD found of ", f);
754    CanonicalForm Ggcd= divide(f, Fgcd,Astar);
755    DEBOUTLN(CERR, "  split into ", Fgcd);
756    DEBOUTLN(CERR, "         and ", Ggcd);
757    Fgcd= pp(Fgcd); Ggcd= pp(Ggcd);
758    DEBDECLEVEL(CERR,"newfactoras");
759    return UnionCFFL(newfactoras(Fgcd,as,success) , newfactoras(Ggcd,as,success));
760  }
761  if ( getCharacteristic() > 0 ){
762
763    // First look for extension!
764    IntList degreelist;
765    Variable vminpoly;
766    for (i=Astar; i.hasItem(); i++){degreelist.append(degree(i.getItem()));}
767    int extdeg= getextension(degreelist, degree(f));
768    DEBOUTLN(CERR, "Extension needed of degree ", extdeg);
769
770    // Now the real stuff!
771    if ( newuord.length() == 0 ){ // no transzendentals
772      DEBOUTMSG(CERR, "No transzendentals!");
773      if ( extdeg > 1 ){
774        CanonicalForm MIPO= generate_mipo( extdeg, vminpoly);
775        DEBOUTLN(CERR, "Minpoly produced ", MIPO);
776        vminpoly= rootOf(MIPO);
777      }
778      Factorlist= alg_factor(f, Astar, vminpoly, oldord, as);
779      DEBDECLEVEL(CERR,"newfactoras");
780      return Factorlist;
781    }
782    else if ( inseperable(Astar) > 0 ){ // Look if extensions are seperable
783      // a) Use Endler
784      DEBOUTMSG(CERR, "Inseperable extensions! Using Endler!");
785      CFFList templist= endler(f,Astar, newuord);
786      DEBOUTLN(CERR, "Endler gives: ", templist);
787      return templist;
788    }
789    else{ // we are on the save side: Use trager
790      DEBOUTMSG(CERR, "Only seperable extensions!");
791      if (extdeg > 1 ){
792        CanonicalForm MIPO=generate_mipo(extdeg, vminpoly );
793        vminpoly= rootOf(MIPO);
794        DEBOUTLN(CERR, "Minpoly generated: ", MIPO);
795        DEBOUTLN(CERR, "vminpoly= ", vminpoly);
796        DEBOUTLN(CERR, "degree(vminpoly)= ", degree(vminpoly));
797      }
798      Factorlist= alg_factor(f, Astar, vminpoly, oldord, as);
799      DEBDECLEVEL(CERR,"newfactoras");
800      return Factorlist;
801    }
802  }
803  else{ // char=0 apply trager directly
804    DEBOUTMSG(CERR, "Char=0! Apply Trager!");
805    Variable vminpoly;
806    Factorlist= alg_factor(f, Astar, vminpoly, oldord, as);
807      DEBDECLEVEL(CERR,"newfactoras");
808      return Factorlist;
809  }
810
811  DEBDECLEVEL(CERR,"newfactoras");
812  return CFFList(CFFactor(f,1));
813}
814
815CFFList
816newcfactor(const CanonicalForm & f, const CFList & as, int success ){
817  Off(SW_RATIONAL);
818  CFFList Output, output, Factors=Factorize(f); On(SW_RATIONAL);
819  Factors.removeFirst();
820
821  if ( as.length() == 0 ){ success=1; return Factors;}
822  if ( cls(f) <= cls(as.getLast()) ) { success=1; return Factors;}
823
824  success=1;
825  for ( CFFListIterator i=Factors; i.hasItem(); i++ ){
826    output=newfactoras(i.getItem().factor(),as, success);
827    for ( CFFListIterator j=output; j.hasItem(); j++)
828      Output = appendCFFL(Output,CFFactor(j.getItem().factor(),j.getItem().exp()*i.getItem().exp()));
829  }
830  return Output;
831}
832
833/*
834$Log: not supported by cvs2svn $
835Revision 1.19  2008/01/25 14:19:39  Singular
836*hannes: SqrFreeTest -> isSqrFree
837
838Revision 1.18  2008/01/22 09:51:36  Singular
839*hannes: sqrFree/InternalSqrFree -> factory
840
841Revision 1.17  2007/05/15 14:46:48  Singular
842*hannes: factorize in Zp(a)[x...]
843
844Revision 1.16  2006/05/16 14:46:48  Singular
845*hannes: gcc 4.1 fixes
846
847Revision 1.15  2005/10/17 13:16:18  Singular
848*hannes: code cleanup
849
850Revision 1.14  2005/07/08 09:18:15  Singular
851*hannes: fixed call of resultant
852
853Revision 1.13  2004/12/10 10:15:05  Singular
854*pohl: AlgExtGenerator etc.
855
856Revision 1.12  2003/02/18 11:09:25  Singular
857* hannes: alg_gcd(f,f'=0) get a special handling
858
859Revision 1.11  2002/10/24 17:22:21  Singular
860* hannes: factoring in alg.ext., alg_gcd, NTL stuff
861
862Revision 1.10  2002/08/19 11:11:29  Singular
863* hannes/pfister: alg_gcd etc.
864
865Revision 1.9  2002/07/30 15:16:19  Singular
866*hannes: fix for alg. extension
867
868Revision 1.8  2001/08/06 08:32:53  Singular
869* hannes: code cleanup
870
871Revision 1.7  2001/06/27 13:58:05  Singular
872*hannes/GP: debug newfactoras, char_series, ...
873
874Revision 1.6  2001/06/21 14:57:04  Singular
875*hannes/GP: Factorize, newfactoras, ...
876
877Revision 1.5  2001/06/18 08:44:39  pfister
878* hannes/GP/michael: factory debug, Factorize
879
880Revision 1.4  1998/05/25 18:32:38  obachman
8811998-05-25  Olaf Bachmann  <obachman@mathematik.uni-kl.de>
882
883        * charset/alg_factor.cc: replaced identifiers 'random' by
884        'myrandom' to avoid name conflicts with the built-in (stdlib)
885        library function 'random' which might be a macro -- and, actually
886        is  under gcc v 2.6.3
887
888Revision 1.3  1998/03/12 12:34:24  schmidt
889        * charset/csutil.cc, charset/alg_factor.cc: all references to
890          `common_den()' replaced by `bCommonDen()'
891
892Revision 1.2  1997/09/12 07:19:37  Singular
893* hannes/michael: libfac-0.3.0
894
895*/
Note: See TracBrowser for help on using the repository browser.