source: git/factory/facAlgFunc.cc @ 5355082

fieker-DuValspielwiese
Last change on this file since 5355082 was 5355082, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: moved helper functions for facAlgFunc from libfac to facAlgFunc
  • Property mode set to 100644
File size: 42.1 KB
Line 
1////////////////////////////////////////////////////////////
2// emacs edit mode for this file is -*- C++ -*-
3////////////////////////////////////////////////////////////
4
5#include "config.h"
6
7#include "cf_assert.h"
8#include "debug.h"
9
10#include "algext.h"
11#include "cf_random.h"
12#include "cf_generator.h"
13#include "cf_irred.h"
14#include "cf_iter.h"
15#include "cf_util.h"
16#include "cf_algorithm.h"
17#include "cf_map.h"
18#include "cfModResultant.h"
19#include "cfCharSets.h"
20#include "facAlgFunc.h"
21
22void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
23
24static
25CanonicalForm
26Prem (const CanonicalForm& F, const CanonicalForm& G)
27{
28  CanonicalForm f, g, l, test, lu, lv, t, retvalue;
29  int degF, degG, levelF, levelG;
30  bool reord;
31  Variable v, vg= G.mvar();
32
33  if ( (levelF= F.level()) < (levelG= G.level()))
34    return F;
35  else
36  {
37    if ( levelF == levelG )
38    {
39      f= F;
40      g= G;
41      reord= false;
42      v= F.mvar();
43    }
44    else
45    {
46      v= Variable (levelF + 1);
47      f= swapvar (F, vg, v);
48      g= swapvar (G, vg, v);
49      reord= true;
50    }
51    degG= degree (g, v );
52    degF= degree (f, v );
53    if (degG <= degF)
54    {
55      l= LC (g);
56      g= g - l*power (v, degG);
57    }
58    else
59      l= 1;
60    while ((degG <= degF) && (!f.isZero()))
61    {
62      test= gcd (l, LC(f));
63      lu= l / test;
64      lv= LC(f) / test;
65
66      t= power (v, degF - degG)*g*lv;
67
68      if (degF == 0)
69        f= 0;
70      else
71        f= f - LC (f)*power (v, degF);
72
73      f= lu*f - t;
74      degF= degree (f, v);
75    }
76
77    if (reord)
78      retvalue= swapvar (f, vg, v);
79    else
80      retvalue= f;
81
82    return retvalue;
83  }
84}
85
86static
87CanonicalForm
88Prem( const CanonicalForm &f, const CFList &L)
89{
90  CanonicalForm rem= f;
91  CFListIterator i= L;
92
93  for (i.lastItem(); i.hasItem(); i--)
94    rem= normalize (Prem (rem, i.getItem()));
95
96  return rem;
97}
98
99static
100CanonicalForm
101Prem (const CanonicalForm &f, const CFList &L, const CFList &as)
102{
103  CanonicalForm rem = f;
104  CFListIterator i = L;
105  for ( i.lastItem(); i.hasItem(); i-- )
106    rem = Prem( rem, i.getItem() );
107  return normalize (rem); //TODO better normalize in case as.length() == 1 && as.getFirst().level() == 1 ???
108}
109
110CFFList
111myappend( const CFFList & Inputlist, const CFFactor & TheFactor)
112{
113  CFFList Outputlist ;
114  CFFactor copy;
115  CFFListIterator i;
116  int exp=0;
117
118  for ( i=Inputlist ; i.hasItem() ; i++ )
119  {
120    copy = i.getItem();
121    if ( copy.factor() == TheFactor.factor() )
122      exp += copy.exp();
123    else
124      Outputlist.append(copy);
125  }
126  Outputlist.append( CFFactor(TheFactor.factor(), exp + TheFactor.exp()));
127  return Outputlist;
128}
129
130CFFList
131myUnion(const CFFList & Inputlist1,const CFFList & Inputlist2)
132{
133  CFFList Outputlist;
134  CFFListIterator i;
135
136  for ( i=Inputlist1 ; i.hasItem() ; i++ )
137    Outputlist = myappend(Outputlist, i.getItem() );
138  for ( i=Inputlist2 ; i.hasItem() ; i++ )
139    Outputlist = myappend(Outputlist, i.getItem() );
140
141  return Outputlist;
142}
143
144///////////////////////////////////////////////////////////////
145// generate a minpoly of degree degree_of_Extension in the   //
146// field getCharacteristik()^Extension.                      //
147///////////////////////////////////////////////////////////////
148CanonicalForm
149generate_mipo( int degree_of_Extension , const Variable & Extension )
150{
151  FFRandom gen;
152  /*if (degree (Extension) < 0)
153    factoryError("libfac: evaluate: Extension not inFF() or inGF() !");*/
154  return find_irreducible( degree_of_Extension, gen, Variable(1) );
155}
156
157static Varlist
158Var_is_in_AS(const Varlist & uord, const CFList & Astar);
159
160////////////////////////////////////////////////////////////////////////
161// This implements the algorithm of Trager for factorization of
162// (multivariate) polynomials over algebraic extensions and so called
163// function field extensions.
164////////////////////////////////////////////////////////////////////////
165
166// // missing class: IntGenerator:
167bool IntGenerator::hasItems() const
168{
169    return 1;
170}
171
172CanonicalForm IntGenerator::item() const
173//int IntGenerator::item() const
174{
175  //return current; //CanonicalForm( current );
176  return mapinto(CanonicalForm( current ));
177}
178
179void IntGenerator::next()
180{
181    current++;
182}
183
184int hasAlgVar(const CanonicalForm &f, const Variable &v)
185{
186  if (f.inBaseDomain()) return 0;
187  if (f.inCoeffDomain())
188  {
189    if (f.mvar()==v) return 1;
190    return hasAlgVar(f.LC(),v);
191  }
192  if (f.inPolyDomain())
193  {
194    if (hasAlgVar(f.LC(),v)) return 1;
195    for( CFIterator i=f; i.hasTerms(); i++)
196    {
197      if (hasAlgVar(i.coeff(),v)) return 1;
198    }
199  }
200  return 0;
201}
202
203int hasVar(const CanonicalForm &f, const Variable &v)
204{
205  if (f.inBaseDomain()) return 0;
206  if (f.inCoeffDomain())
207  {
208    if (f.mvar()==v) return 1;
209    return hasAlgVar(f.LC(),v);
210  }
211  if (f.inPolyDomain())
212  {
213    if (f.mvar()==v) return 1;
214    if (hasVar(f.LC(),v)) return 1;
215    for( CFIterator i=f; i.hasTerms(); i++)
216    {
217      if (hasVar(i.coeff(),v)) return 1;
218    }
219  }
220  return 0;
221}
222
223int hasAlgVar(const CanonicalForm &f)
224{
225  if (f.inBaseDomain()) return 0;
226  if (f.inCoeffDomain())
227  {
228    if (f.level()!=0)
229      return 1;
230    return hasAlgVar(f.LC());
231  }
232  if (f.inPolyDomain())
233  {
234    if (hasAlgVar(f.LC())) return 1;
235    for( CFIterator i=f; i.hasTerms(); i++)
236    {
237      if (hasAlgVar(i.coeff())) return 1;
238    }
239  }
240  return 0;
241}
242
243/// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r
244void
245psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, 
246      CanonicalForm & r, CanonicalForm& multiplier, const Variable& x)
247{
248    ASSERT( x.level() > 0, "type error: polynomial variable expected" );
249    ASSERT( ! g.isZero(), "math error: division by zero" );
250
251    // swap variables such that x's level is larger or equal
252    // than both f's and g's levels.
253    Variable X;
254    if (f.level() > g.level())
255      X= f.mvar();
256    else
257      X= g.mvar();
258    if (X.level() < x.level())
259      X= x;
260    CanonicalForm F= swapvar (f, x, X);
261    CanonicalForm G= swapvar (g, x, X);
262
263    // now, we have to calculate the pseudo remainder of F and G
264    // w.r.t. X
265    int fDegree= degree (F, X);
266    int gDegree= degree (G, X);
267    if (fDegree < 0 || fDegree < gDegree)
268    {
269      q= 0;
270      r= f;
271    }
272    else
273    {
274      CanonicalForm LCG= LC (G, X);
275      multiplier= power (LCG, fDegree - gDegree + 1);
276      divrem (multiplier*F, G, q, r);
277      q= swapvar (q, x, X);
278      r= swapvar (r, x, X);
279    }
280}
281
282static CanonicalForm
283Sprem ( const CanonicalForm &f, const CanonicalForm &g, CanonicalForm & m, CanonicalForm & q )
284{
285  CanonicalForm ff, gg, l, test, retvalue;
286  int df, dg,n;
287  bool reord;
288  Variable vf, vg, v;
289
290  if ( (vf = f.mvar()) < (vg = g.mvar()) )
291  {
292    m=CanonicalForm(0); q=CanonicalForm(0);
293    return f;
294  }
295  else
296  {
297    if ( vf == vg )
298    {
299      ff = f; gg = g;
300      reord = false;
301      v = vg; // == x
302    }
303    else
304    {
305      v = Variable(level(f.mvar()) + 1);
306      ff = swapvar(f,vg,v); // == r
307      gg = swapvar(g,vg,v); // == v
308      reord=true;
309    }
310    dg = degree( gg, v ); // == dv
311    df = degree( ff, v ); // == dr
312    if (dg <= df) {l=LC(gg); gg = gg -LC(gg)*power(v,dg);}
313    else { l = 1; }
314    n= 0;
315    while ( ( dg <= df  ) && ( !ff.isZero()) )
316    {
317      test= power(v,df-dg) * gg * LC(ff);
318      if ( df == 0 ){ff= ff.genZero();}
319      else {ff= ff - LC(ff)*power(v,df);}
320      ff = l*ff-test;
321      df= degree(ff,v);
322      n++;
323    }
324    if ( reord )
325    {
326      retvalue= swapvar( ff, vg, v );
327    }
328    else
329    {
330      retvalue= ff;
331    }
332    m= power(l,n);
333    if ( fdivides(g,m*f-retvalue) )
334      q= (m*f-retvalue)/g;
335    else
336      q= CanonicalForm(0);
337    return retvalue;
338  }
339}
340
341CanonicalForm
342divide( const CanonicalForm & ff, const CanonicalForm & f, const CFList & as)
343{
344  CanonicalForm r,m,q;
345
346  if (f.inCoeffDomain())
347  {
348    bool isRat=isOn(SW_RATIONAL);
349    if (getCharacteristic() == 0)
350      On(SW_RATIONAL);
351    q=ff/f;
352    if (!isRat && getCharacteristic() == 0)
353      Off(SW_RATIONAL);
354  }
355  else
356    r= Sprem(ff,f,m,q);
357
358  r= Prem(q,as);
359  return r;
360}
361
362CanonicalForm
363alg_content (const CanonicalForm& f, const CFList& as)
364{
365  if (!f.inCoeffDomain())
366  {
367    CFIterator i= f;
368    CanonicalForm result= abs (i.coeff());
369    i++;
370    while (i.hasTerms() && !result.isOne())
371    {
372      result= alg_gcd (i.coeff(), result, as);
373      i++;
374    }
375    return result;
376  }
377
378  return abs (f);
379}
380
381CanonicalForm alg_gcd(const CanonicalForm & fff, const CanonicalForm &ggg,
382                      const CFList &as)
383{
384  if (fff.inCoeffDomain() || ggg.inCoeffDomain())
385    return 1;
386  CanonicalForm f=fff;
387  CanonicalForm g=ggg;
388  f=Prem(f,as);
389  g=Prem(g,as);
390  if ( f.isZero() )
391  {
392    if ( g.lc().sign() < 0 ) return -g;
393    else                     return g;
394  }
395  else  if ( g.isZero() )
396  {
397    if ( f.lc().sign() < 0 ) return -f;
398    else                     return f;
399  }
400
401  int v= as.getLast().level();
402  if (f.level() <= v || g.level() <= v)
403    return 1;
404
405  CanonicalForm res;
406
407  // does as appear in f and g ?
408  bool has_alg_var=false;
409  for ( CFListIterator j=as;j.hasItem(); j++ )
410  {
411    Variable v=j.getItem().mvar();
412    if (hasVar (f, v))
413      has_alg_var=true;
414    if (hasVar (g, v))
415      has_alg_var=true;
416  }
417  if (!has_alg_var)
418  {
419    /*if ((hasAlgVar(f))
420    || (hasAlgVar(g)))
421    {
422      Varlist ord;
423      for ( CFListIterator j=as;j.hasItem(); j++ )
424        ord.append(j.getItem().mvar());
425      res=algcd(f,g,as,ord);
426    }
427    else*/
428    if (!hasAlgVar (f) && !hasAlgVar (g))
429      res=gcd(f,g);
430
431    return res;
432  }
433
434  int mvf=f.level();
435  int mvg=g.level();
436  if (mvg > mvf)
437  {
438    CanonicalForm tmp=f; f=g; g=tmp;
439    int tmp2=mvf; mvf=mvg; mvg=tmp2;
440  }
441  if (g.inBaseDomain() || f.inBaseDomain())
442    return CanonicalForm(1);
443
444  CanonicalForm c_f= alg_content (f, as);
445
446  if (mvf != mvg)
447  {
448    res= alg_gcd (g, c_f, as);
449    return res;
450  }
451  Variable x= f.mvar();
452
453  // now: mvf==mvg, f.level()==g.level()
454  CanonicalForm c_g= alg_content (g, as);
455
456  int delta= degree (f) - degree (g);
457
458  f= divide (f, c_f, as);
459  g= divide (g, c_g, as);
460
461  // gcd of contents
462  CanonicalForm c_gcd= alg_gcd (c_f, c_g, as);
463  CanonicalForm tmp;
464
465  if (delta < 0)
466  {
467    tmp= f;
468    f= g;
469    g= tmp;
470    delta= -delta;
471  }
472
473  CanonicalForm r= 1;
474
475  while (degree (g, x) > 0)
476  {
477    r= Prem (f, g, as);
478    r= Prem (r, as);
479    if (!r.isZero())
480    {
481      r= divide (r, alg_content (r,as), as);
482      r /= vcontent (r,Variable (v+1));
483    }
484    f= g;
485    g= r;
486  }
487
488  if (degree (g, x) == 0)
489    return c_gcd;
490
491  c_f= alg_content (f, as);
492
493  f= divide (f, c_f, as);
494
495  f *= c_gcd;
496  f /= vcontent (f, Variable (v+1));
497
498  return f;
499}
500
501static CanonicalForm
502resultante( const CanonicalForm & f, const CanonicalForm& g, const Variable & v )
503{
504  bool on_rational = isOn(SW_RATIONAL);
505  if (!on_rational && getCharacteristic() == 0)
506    On(SW_RATIONAL);
507  CanonicalForm cd = bCommonDen( f );
508  CanonicalForm fz = f * cd;
509  cd = bCommonDen( g );
510  CanonicalForm gz = g * cd;
511  if (!on_rational && getCharacteristic() == 0)
512    Off(SW_RATIONAL);
513  CanonicalForm result;
514  if (getCharacteristic() == 0)
515    result= resultantZ (fz, gz,v);
516  else
517    result= resultant (fz,gz,v);
518
519  return result;
520}
521
522// sqr-free routine for algebraic extensions
523// we need it! Ex.: f=c^2+2*a*c-1; as=[a^2+1]; f=(c+a)^2
524//static CFFList alg_sqrfree( const CanonicalForm & f )
525//{
526//  CFFList L;
527//
528//  L.append(CFFactor(f,1));
529//  return L;
530//}
531
532// Calculates a square free norm
533// Input: f(x, alpha) a square free polynomial over K(alpha),
534// alpha is defined by the minimal polynomial Palpha
535// K has more than S elements (S is defined in thesis; look getextension)
536static void
537sqrf_norm_sub( const CanonicalForm & f, const CanonicalForm & PPalpha,
538           CFGenerator & myrandom, CanonicalForm & s,  CanonicalForm & g,
539           CanonicalForm & R)
540{
541  Variable y=PPalpha.mvar(),vf=f.mvar();
542  CanonicalForm temp, Palpha=PPalpha, t;
543  int sqfreetest=0;
544  CFFList testlist;
545  CFFListIterator i;
546
547  myrandom.reset();   s=myrandom.item();   g=f;
548  R= CanonicalForm(0);
549
550  // Norm, resultante taken with respect to y
551  while ( !sqfreetest )
552  {
553    R = resultante(Palpha, g, y); R= R* bCommonDen(R);
554    R /= content (R);
555    // sqfree check ; R is a polynomial in K[x]
556    if ( getCharacteristic() == 0 )
557    {
558      temp= gcd(R, R.deriv(vf));
559      if (degree(temp,vf) != 0 || temp == temp.genZero() )
560        sqfreetest= 0;
561      else
562        sqfreetest= 1;
563    }
564    else
565    {
566      Variable X;
567      if (hasFirstAlgVar(R,X))
568        testlist=factorize( R, X );
569      else
570        testlist= factorize(R);
571
572      if (testlist.getFirst().factor().inCoeffDomain())
573        testlist.removeFirst();
574      sqfreetest=1;
575      for ( i=testlist; i.hasItem(); i++)
576      {
577        if ( i.getItem().exp() > 1 && degree(i.getItem().factor(), R.mvar()) > 0)
578        {
579          sqfreetest=0;
580          break;
581        }
582      }
583    }
584    if ( ! sqfreetest )
585    {
586      myrandom.next();
587      if ( getCharacteristic() == 0 )
588        t= CanonicalForm(mapinto(myrandom.item()));
589      else
590        t= CanonicalForm(myrandom.item());
591      s= t;
592      g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
593    }
594  }
595}
596static void
597sqrf_agnorm_sub( const CanonicalForm & f, const CanonicalForm & PPalpha,
598           AlgExtGenerator & myrandom, CanonicalForm & s,  CanonicalForm & g,
599           CanonicalForm & R)
600{
601  Variable y=PPalpha.mvar(),vf=f.mvar();
602  CanonicalForm temp, Palpha=PPalpha, t;
603  int sqfreetest=0;
604  CFFList testlist;
605  CFFListIterator i;
606
607  myrandom.reset();   s=myrandom.item();   g=f;
608  R= CanonicalForm(0);
609
610  // Norm, resultante taken with respect to y
611  while ( !sqfreetest )
612  {
613    R = resultante(Palpha, g, y); R= R* bCommonDen(R);
614    R /= content (R);
615    // sqfree check ; R is a polynomial in K[x]
616    if ( getCharacteristic() == 0 )
617    {
618      temp= gcd(R, R.deriv(vf));
619      if (degree(temp,vf) != 0 || temp == temp.genZero() )
620        sqfreetest= 0;
621      else
622        sqfreetest= 1;
623    }
624    else
625    {
626      Variable X;
627      if (hasFirstAlgVar(R,X))
628        testlist= factorize( R, X );
629      else
630        testlist= factorize(R);
631      if (testlist.getFirst().factor().inCoeffDomain())
632        testlist.removeFirst();
633      sqfreetest=1;
634      for ( i=testlist; i.hasItem(); i++)
635      {
636        if ( i.getItem().exp() > 1 && degree(i.getItem().factor(), R.mvar()) > 0)
637        {
638          sqfreetest=0;
639          break;
640        }
641      }
642    }
643    if ( ! sqfreetest )
644    {
645      myrandom.next();
646      if ( getCharacteristic() == 0 )
647        t= CanonicalForm(mapinto(myrandom.item()));
648      else
649        t= CanonicalForm(myrandom.item());
650      s= t;
651      g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
652    }
653  }
654}
655
656static void
657sqrf_norm( const CanonicalForm & f, const CanonicalForm & PPalpha,
658           const Variable & Extension, CanonicalForm & s,  CanonicalForm & g,
659           CanonicalForm & R)
660{
661  if ( getCharacteristic() == 0 )
662  {
663    IntGenerator myrandom;
664    sqrf_norm_sub(f,PPalpha, myrandom, s,g,R);
665  }
666  else if ( degree(Extension) > 0 ) // working over Extensions
667  {
668    AlgExtGenerator myrandom(Extension);
669    sqrf_agnorm_sub(f,PPalpha, myrandom, s,g,R);
670  }
671  else
672  {
673    FFGenerator myrandom;
674    sqrf_norm_sub(f,PPalpha, myrandom, s,g,R);
675  }
676}
677
678static Varlist
679Var_is_in_AS(const Varlist & uord, const CFList & Astar){
680  Varlist output;
681  CanonicalForm elem;
682  Variable x;
683
684  for ( VarlistIterator i=uord; i.hasItem(); i++)
685  {
686    x=i.getItem();
687    for ( CFListIterator j=Astar; j.hasItem(); j++ )
688    {
689      elem= j.getItem();
690      if ( degree(elem,x) > 0 ) // x actually occures in Astar
691      {
692        output.append(x);
693        break;
694      }
695    }
696  }
697  return output;
698}
699
700// Look if Minimalpolynomials in Astar define seperable Extensions
701// Must be a power of p: i.e. y^{p^e}-x
702static int
703inseperable(const CFList & Astar)
704{
705  CanonicalForm elem;
706  int Counter= 1;
707
708  if ( Astar.length() == 0 ) return 0;
709  for ( CFListIterator i=Astar; i.hasItem(); i++)
710  {
711    elem= i.getItem();
712    if ( elem.deriv() == elem.genZero() ) return Counter;
713    else Counter += 1;
714  }
715  return 0;
716}
717
718// calculate big enough extension for finite fields
719// Idea: first calculate k, such that q^k > S (->thesis, -> getextension)
720// Second, search k with gcd(k,m_i)=1, where m_i is the degree of the i'th
721// minimal polynomial. Then the minpoly f_i remains irrd. over q^k and we
722// have enough elements to plug in.
723static int
724getextension( IntList & degreelist, int n)
725{
726  int charac= getCharacteristic();
727  setCharacteristic(0); // need it for k !
728  int k=1, m=1, length=degreelist.length();
729  IntListIterator i;
730
731  for (i=degreelist; i.hasItem(); i++) m= m*i.getItem();
732  int q=charac;
733  while (q <= ((n*m)*(n*m)/2)) { k= k+1; q= q*charac;}
734  int l=0;
735  do {
736    for (i=degreelist; i.hasItem(); i++){
737      l= l+1;
738      if ( igcd(k,i.getItem()) == 1){
739        if ( l==length ){ setCharacteristic(charac);  return k; }
740      }
741      else { break; }
742    }
743    k= k+1; l=0;
744  }
745  while ( 1 );
746}
747
748CanonicalForm
749QuasiInverse (const CanonicalForm& f, const CanonicalForm& g,
750              const Variable& x)
751{
752  CanonicalForm pi, pi1, q, t0, t1, Hi, bi, pi2;
753  bool isRat= isOn (SW_RATIONAL);
754  CanonicalForm m,tmp;
755  if (isRat && getCharacteristic() == 0)
756    Off (SW_RATIONAL);
757  pi= f/content (f,x);
758  pi1= g/content (g,x);
759
760  t0= 0;
761  t1= 1;
762  bi= 1;
763
764  int delta= degree (f, x) - degree (g, x);
765  Hi= power (LC (pi1, x), delta);
766  if ( (delta+1) % 2 )
767      bi = 1;
768  else
769      bi = -1;
770
771  while (degree (pi1,x) > 0)
772  {
773    psqr( pi, pi1, q, pi2, m, x);
774    pi2 /= bi;
775
776    tmp= t1;
777    t1= t0*m - t1*q;
778    t0= tmp;
779    t1 /= bi;
780    pi = pi1; pi1 = pi2;
781    if ( degree ( pi1, x ) > 0 )
782    {
783      delta = degree( pi, x ) - degree( pi1, x );
784      if ( (delta+1) % 2 )
785        bi = LC( pi, x ) * power( Hi, delta );
786      else
787        bi = -LC( pi, x ) * power( Hi, delta );
788      Hi = power( LC( pi1, x ), delta ) / power( Hi, delta-1 );
789    }
790  }
791  t1 /= gcd (pi1, t1);
792  if (!isRat && getCharacteristic() == 0)
793    Off (SW_RATIONAL);
794  return t1;
795}
796
797CanonicalForm
798evaluate (const CanonicalForm& f, const CanonicalForm& g, const CanonicalForm& h, const CanonicalForm& powH)
799{
800  if (f.inCoeffDomain())
801    return f;
802  CFIterator i= f;
803  int lastExp = i.exp();
804  CanonicalForm result = i.coeff()*powH;
805  i++;
806  while (i.hasTerms())
807  {
808    int i_exp= i.exp();
809    if ((lastExp - i_exp) == 1)
810    {
811      result *= g;
812      result /= h;
813    }
814    else
815    {
816      result *= power (g, lastExp - i_exp);
817      result /= power (h, lastExp - i_exp);
818    }
819    result += i.coeff()*powH;
820    lastExp = i_exp;
821    i++;
822  }
823  if (lastExp != 0)
824  {
825    result *= power (g, lastExp);
826    result /= power (h, lastExp);
827  }
828  return result;
829}
830
831
832/// evaluate f at g/h at v such that powH*f is integral i.e. powH is assumed to be h^degree(f,v)
833CanonicalForm
834evaluate (const CanonicalForm& f, const CanonicalForm& g,
835          const CanonicalForm& h, const CanonicalForm& powH,
836          const Variable& v)
837{
838  if (f.inCoeffDomain())
839  {
840    return f*powH;
841  }
842
843  Variable x = f.mvar();
844  if ( v > x )
845    return f*powH;
846  else  if ( v == x )
847    return evaluate (f, g, h, powH);
848
849  // v is less than main variable of f
850  CanonicalForm result= 0;
851  for (CFIterator i= f; i.hasTerms(); i++)
852    result += evaluate (i.coeff(), g, h, powH, v)*power (x, i.exp());
853  return result;
854}
855
856// calculate a "primitive element"
857// K must have more than S elements (->thesis, -> getextension)
858static CFList
859simpleextension(CFList& backSubst, const CFList & Astar,
860                const Variable & Extension, bool& isFunctionField,
861                CanonicalForm & R)
862{
863  CFList Returnlist, Bstar=Astar;
864  CanonicalForm s, g, ra, rb, oldR, h, denra, denrb=1;
865  Variable alpha;
866  CFList tmp;
867
868  bool isRat= isOn (SW_RATIONAL);
869
870  CFListIterator j;
871  if (Astar.length() == 1)
872  {
873    R= Astar.getFirst();
874    rb= R.mvar();
875  }
876  else
877  {
878    R=Bstar.getFirst();
879    Bstar.removeFirst();
880    for (CFListIterator i=Bstar; i.hasItem(); i++)
881    {
882      j= i;
883      j++;
884      if (getCharacteristic() == 0)
885        Off (SW_RATIONAL);
886      R /= icontent (R);
887      if (getCharacteristic() == 0)
888        On (SW_RATIONAL);
889      oldR= R;
890      sqrf_norm (i.getItem(), R, Extension, s, g, R);
891
892      backSubst.insert (s);
893
894      if (getCharacteristic() == 0)
895        Off (SW_RATIONAL);
896      R /= icontent (R);
897
898      if (getCharacteristic() == 0)
899        On (SW_RATIONAL);
900
901      if (!isFunctionField)
902      {
903        alpha= rootOf (R);
904        h= replacevar (g, g.mvar(), alpha);
905        if (getCharacteristic() == 0)
906          On (SW_RATIONAL); //needed for GCD
907        h= gcd (h, oldR);
908        h /= Lc (h);
909        ra= -h[0];
910        ra= replacevar(ra, alpha, g.mvar());
911        rb= R.mvar()-s*ra;
912        for (; j.hasItem(); j++)
913        {
914          j.getItem()= j.getItem() (ra, oldR.mvar());
915          j.getItem()= j.getItem() (rb, i.getItem().mvar());
916        }
917      }
918      else
919      {
920        if (getCharacteristic() == 0)
921          On (SW_RATIONAL);
922        h= swapvar (g, g.mvar(), oldR.mvar());
923        tmp= CFList (swapvar (R, g.mvar(), oldR.mvar()));
924        h= alg_gcd (h, swapvar (oldR, g.mvar(), oldR.mvar()), tmp);
925        CanonicalForm hh= replacevar (h, oldR.mvar(), alpha);
926
927        CanonicalForm numinv, deninv;
928        numinv= QuasiInverse (tmp.getFirst(), LC (h), tmp.getFirst().mvar());
929
930        if (getCharacteristic() == 0)
931          Off (SW_RATIONAL);
932        h *= numinv;
933        h= reduce (h, tmp.getFirst());
934        deninv= LC(h);
935
936        ra= -h[0];
937        denra= gcd (ra, deninv);
938        ra /= denra;
939        denra= deninv/denra;
940        denra= replacevar (denra, ra.mvar(), g.mvar());
941        ra= replacevar(ra, ra.mvar(), g.mvar());
942        rb= R.mvar()*denra-s*ra;
943        denrb= denra;
944        for (; j.hasItem(); j++)
945        {
946          CanonicalForm powdenra= power (denra, degree (j.getItem(), oldR.mvar()));
947          j.getItem()= evaluate (j.getItem(),ra, denra, powdenra, oldR.mvar());
948          powdenra= power (denra, degree (j.getItem(), i.getItem().mvar()));
949          j.getItem()= evaluate (j.getItem(), rb, denrb, powdenra, i.getItem().mvar());
950        }
951      }
952
953      Returnlist.append (ra);
954      if (isFunctionField)
955        Returnlist.append (denra);
956    }
957  }
958  Returnlist.append (rb);
959  if (isFunctionField)
960    Returnlist.append (denrb);
961
962  if (isRat && getCharacteristic() == 0)
963    On (SW_RATIONAL);
964  else if (!isRat && getCharacteristic() == 0)
965    Off (SW_RATIONAL);
966
967  return Returnlist;
968}
969
970CanonicalForm alg_lc(const CanonicalForm &f)
971{
972  if (f.level()>0)
973  {
974    return alg_lc(f.LC());
975  }
976
977  return f;
978}
979
980CanonicalForm alg_LC (const CanonicalForm& f, int lev)
981{
982  CanonicalForm result= f;
983  while (result.level() > lev)
984    result= LC (result);
985  return result;
986}
987
988CanonicalForm
989subst (const CanonicalForm& f, const CFList& a, const CFList& b,
990       const CanonicalForm& Rstar, bool isFunctionField)
991{
992  if (isFunctionField)
993    ASSERT (2*a.length() == b.length(), "wrong length of lists");
994  else
995    ASSERT (a.length() == b.length(), "lists of equal length expected");
996  CFListIterator j= b;
997  CanonicalForm result= f, tmp, powj;
998  for (CFListIterator i= a; i.hasItem() && j.hasItem(); i++, j++)
999  {
1000    if (!isFunctionField)
1001      result= result (j.getItem(), i.getItem().mvar());
1002    else
1003    {
1004      tmp= j.getItem();
1005      j++;
1006      powj= power (j.getItem(), degree (result, i.getItem().mvar()));
1007      result= evaluate (result, tmp, j.getItem(), powj, i.getItem().mvar());
1008
1009      if (fdivides (powj, result, tmp))
1010        result= tmp;
1011
1012      result /= vcontent (result, Variable (i.getItem().level() + 1));
1013    }
1014  }
1015  result= reduce (result, Rstar);
1016  result /= vcontent (result, Variable (Rstar.level() + 1));
1017  return result;
1018}
1019
1020CanonicalForm
1021backSubst (const CanonicalForm& F, const CFList& a, const CFList& b)
1022{
1023  ASSERT (a.length() == b.length() - 1, "wrong length of lists in backSubst");
1024  CanonicalForm result= F;
1025  Variable tmp;
1026  CFList tmp2= b;
1027  tmp= tmp2.getLast().mvar();
1028  tmp2.removeLast();
1029  for (CFListIterator iter= a; iter.hasItem(); iter++)
1030  {
1031    result= result (tmp+iter.getItem()*tmp2.getLast().mvar(), tmp);
1032    tmp= tmp2.getLast().mvar();
1033    tmp2.removeLast();
1034  }
1035  return result;
1036}
1037
1038void deflateDegree (const CanonicalForm & F, int & pExp, int n)
1039{
1040  if (n == 0 || n > F.level())
1041  {
1042    pExp= -1;
1043    return;
1044  }
1045  if (F.level() == n)
1046  {
1047    ASSERT (F.deriv().isZero(), "derivative of F is not zero");
1048    int termCount=0;
1049    CFIterator i= F;
1050    for (; i.hasTerms(); i++)
1051    {
1052      if (i.exp() != 0)
1053        termCount++;
1054    }
1055
1056    int j= 1;
1057    i= F;
1058    for (;j < termCount; j++, i++)
1059      ;
1060
1061    int exp= i.exp();
1062    int count= 0;
1063    int p= getCharacteristic();
1064    while ((exp >= p) && (exp != 0) && (exp % p == 0))
1065    {
1066      exp /= p;
1067      count++;
1068    }
1069    pExp= count;
1070  }
1071  else
1072  {
1073    CFIterator i= F;
1074    deflateDegree (i.coeff(), pExp, n);
1075    i++;
1076    int tmp= pExp;
1077    for (; i.hasTerms(); i++)
1078    {
1079      deflateDegree (i.coeff(), pExp, n);
1080      if (tmp == -1)
1081        tmp= pExp;
1082      else if (tmp != -1 && pExp != -1)
1083        pExp= (pExp < tmp) ? pExp : tmp;
1084      else
1085        pExp= tmp;
1086    }
1087  }
1088}
1089
1090CanonicalForm deflatePoly (const CanonicalForm & F, int exp)
1091{
1092  if (exp == 0)
1093    return F;
1094  int p= getCharacteristic();
1095  int pToExp= ipower (p, exp);
1096  Variable x=F.mvar();
1097  CanonicalForm result= 0;
1098  for (CFIterator i= F; i.hasTerms(); i++)
1099    result += i.coeff()*power (x, i.exp()/pToExp);
1100  return result;
1101}
1102
1103CanonicalForm deflatePoly (const CanonicalForm & F, int exps, int n)
1104{
1105  if (n == 0 || exps <= 0 || F.level() < n)
1106    return F;
1107  if (F.level() == n)
1108    return deflatePoly (F, exps);
1109  else
1110  {
1111    CanonicalForm result= 0;
1112    for (CFIterator i= F; i.hasTerms(); i++)
1113      result += deflatePoly (i.coeff(), exps, n)*power(F.mvar(), i.exp());
1114    return result;
1115  }
1116}
1117
1118CanonicalForm inflatePoly (const CanonicalForm & F, int exp)
1119{
1120  if (exp == 0)
1121    return F;
1122  int p= getCharacteristic();
1123  int pToExp= ipower (p, exp);
1124  Variable x=F.mvar();
1125  CanonicalForm result= 0;
1126  for (CFIterator i= F; i.hasTerms(); i++)
1127    result += i.coeff()*power (x, i.exp()*pToExp);
1128  return result;
1129}
1130
1131CanonicalForm inflatePoly (const CanonicalForm & F, int exps, int n)
1132{
1133  if (n == 0 || exps <= 0 || F.level() < n)
1134    return F;
1135  if (F.level() == n)
1136    return inflatePoly (F, exps);
1137  else
1138  {
1139    CanonicalForm result= 0;
1140    for (CFIterator i= F; i.hasTerms(); i++)
1141      result += inflatePoly (i.coeff(), exps, n)*power(F.mvar(), i.exp());
1142    return result;
1143  }
1144}
1145
1146// the heart of the algorithm: the one from Trager
1147static CFFList
1148alg_factor( const CanonicalForm & F, const CFList & Astar, const Variable & vminpoly, const Varlist & oldord, const CFList & as, bool isFunctionField)
1149{
1150  bool isRat= isOn (SW_RATIONAL);
1151  CFFList L, Factorlist;
1152  CanonicalForm R, Rstar, s, g, h, f= F;
1153  CFList substlist, backSubsts;
1154
1155  substlist= simpleextension (backSubsts, Astar, vminpoly, isFunctionField, Rstar);
1156
1157  f= subst (f, Astar, substlist, Rstar, isFunctionField);
1158
1159  Variable alpha;
1160  if (!isFunctionField)
1161  {
1162    alpha= rootOf (Rstar);
1163    g= replacevar (f, Rstar.mvar(), alpha);
1164
1165    if (!isRat && getCharacteristic() == 0)
1166      On (SW_RATIONAL);
1167    Factorlist= factorize (g, alpha);
1168
1169    for (CFFListIterator i= Factorlist; i.hasItem(); i++)
1170    {
1171      h= i.getItem().factor();
1172      if (!h.inCoeffDomain())
1173      {
1174        h= replacevar (h, alpha, Rstar.mvar());
1175        h *= bCommonDen(h);
1176        h= backSubst (h, backSubsts, Astar);
1177        h= Prem (h, as);
1178        L.append (CFFactor (h, i.getItem().exp()));
1179      }
1180    }
1181    if (!isRat && getCharacteristic() == 0)
1182      Off (SW_RATIONAL);
1183    return L;
1184  }
1185  // after here we are over an extension of a function field
1186
1187  // make quasi monic
1188  CFList Rstarlist= CFList (Rstar);
1189  int algExtLevel= Astar.getLast().level(); //highest level of algebraic variables
1190  CanonicalForm numinv;
1191  if (!isRat && getCharacteristic() == 0)
1192    On (SW_RATIONAL);
1193  numinv= QuasiInverse (Rstar, alg_LC(f, algExtLevel), Rstar.mvar());
1194
1195  f *= numinv;
1196  f= Prem (f, Rstarlist);
1197  f /= vcontent (f, Rstar.mvar());
1198  // end quasi monic
1199
1200  if (degree (f) == 1)
1201  {
1202    f= backSubst (f, backSubsts, Astar);
1203    f *= bCommonDen (f);
1204    f= Prem (f, as);
1205    f /= vcontent (f, as.getFirst().mvar());
1206
1207    return CFFList (CFFactor (f, 1));
1208  }
1209
1210  sqrf_norm(f, Rstar, vminpoly, s, g, R );
1211
1212  Variable X;
1213  if (hasFirstAlgVar(R,X))
1214  {
1215    // factorize R over alg.extension with X
1216    Factorlist =  factorize( R, X );
1217  }
1218  else
1219  {
1220    // factor R over k
1221    Factorlist = factorize(R);
1222  }
1223
1224  if ( !Factorlist.getFirst().factor().inCoeffDomain() )
1225    Factorlist.insert(CFFactor(1,1));
1226  if ( Factorlist.length() == 2 && Factorlist.getLast().exp()== 1)
1227  { // irreduzibel (first entry is a constant)
1228    f= backSubst (f, backSubsts, Astar);
1229    f *= bCommonDen (f);
1230    f= Prem (f, as);
1231    f /= vcontent (f, as.getFirst().mvar());
1232
1233    L.append(CFFactor(f,1));
1234  }
1235  else
1236  {
1237    g= f;
1238    for ( CFFListIterator i=Factorlist; i.hasItem(); i++)
1239    {
1240      CanonicalForm fnew=i.getItem().factor();
1241      if (fnew.level() < Rstar.level()) //factor is a constant from the function field
1242        continue;
1243      else
1244      {
1245        fnew= fnew (g.mvar()+s*Rstar.mvar(), g.mvar());
1246        fnew= reduce (fnew, Rstar);
1247      }
1248
1249      h= alg_gcd (g, fnew, Rstarlist);
1250      numinv= QuasiInverse(Rstar, alg_LC(h, algExtLevel), Rstar.mvar());
1251      h *= numinv;
1252      h= Prem (h, Rstarlist);
1253      h /= vcontent (h, Rstar.mvar());
1254
1255      if (h.level() >= Rstar.level())
1256      {
1257        g= divide (g, h, Rstarlist);
1258        h= backSubst (h, backSubsts, Astar);
1259        h= Prem (h, as);
1260        h *= bCommonDen (h);
1261        h /= vcontent (h, as.getFirst().mvar());
1262        L.append (CFFactor (h, 1));
1263      }
1264    }
1265    // we are not interested in a
1266    // constant (over K_r, which can be a polynomial!)
1267    if (degree(g, f.mvar())>0){ L.append(CFFactor(g,1)); }
1268  }
1269  CFFList LL;
1270  if (getCharacteristic()>0)
1271  {
1272    CFFListIterator i=L;
1273    CanonicalForm c_fac=1;
1274    CanonicalForm c;
1275    for(;i.hasItem(); i++ )
1276    {
1277      CanonicalForm ff=i.getItem().factor();
1278      c=alg_lc(ff);
1279      int e=i.getItem().exp();
1280      ff/=c;
1281      if (!ff.isOne()) LL.append(CFFactor(ff,e));
1282      while (e>0) { c_fac*=c;e--; }
1283    }
1284    if (!c_fac.isOne()) LL.insert(CFFactor(c_fac,1));
1285  }
1286  else
1287  {
1288    LL=L;
1289  }
1290
1291  if (!isRat && getCharacteristic() == 0)
1292    Off (SW_RATIONAL);
1293
1294  return LL;
1295}
1296
1297/// algorithm of A. Steel described in "Conquering Inseparability: Primary
1298/// decomposition and multivariate factorization over algebraic function fields
1299/// of positive characteristic" with the following modifications: we use
1300/// characteristic sets in IdealOverSubfield and Trager's primitive element
1301/// algorithm instead of FGLM
1302CFFList
1303SteelTrager (const CanonicalForm & f, const CFList & AS, const Varlist & uord)
1304{
1305  CanonicalForm F= f, lcmVars= 1;
1306  CFList asnew, as= AS;
1307  CFListIterator i, ii;
1308
1309  bool derivZeroF= false, derivZero= false;
1310  int j, exp=0, expF= 0, tmpExp;
1311  CFFList varsMapLevel;
1312  CFFListIterator iter;
1313
1314  // check if F is separable
1315  if (F.deriv().isZero())
1316  {
1317    derivZeroF= true;
1318    deflateDegree (F, expF, F.level());
1319  }
1320
1321  CanonicalForm varsF= getVars (F);
1322  varsF /= F.mvar();
1323
1324  lcmVars= lcm (varsF, lcmVars);
1325
1326  if (derivZeroF)
1327    as.append (F);
1328
1329  CFFList varsGMapLevel;
1330  CFFList tmp;
1331  CFFList * varsGMap= new CFFList [as.length()];
1332  for (j= 0; j < as.length(); j++)
1333    varsGMap[j]= CFFList();
1334
1335  CanonicalForm varsG;
1336  j= 0;
1337  bool recurse= false;
1338  i= as;
1339  // make minimal polys and F separable
1340  while (i.hasItem())
1341  {
1342    derivZero= false;
1343
1344    if (i.getItem().deriv() == 0)
1345    {
1346      deflateDegree (i.getItem(), exp, i.getItem().level());
1347      i.getItem()= deflatePoly (i.getItem(), exp, i.getItem().level());
1348     
1349      varsG= getVars (i.getItem());
1350      varsG /= i.getItem().mvar();
1351
1352      lcmVars= lcm (varsG, lcmVars);
1353
1354      while (!varsG.isOne())
1355      {
1356        if (i.getItem().deriv (varsG.level()).isZero())
1357        {
1358          deflateDegree (i.getItem(), tmpExp, varsG.level());
1359          if (exp >= tmpExp)
1360          {
1361            if (exp == tmpExp)
1362              i.getItem()= deflatePoly (i.getItem(), exp, varsG.level());
1363            else
1364            {
1365              if (j != 0)
1366                recurse= true;
1367              i.getItem()= deflatePoly (i.getItem(), tmpExp, varsG.level());
1368            }
1369            varsGMapLevel.insert (CFFactor (varsG.mvar(), exp - tmpExp));
1370          }
1371          else
1372          {
1373            i.getItem()= deflatePoly (i.getItem(), exp, varsG.level());
1374            varsGMapLevel.insert (CFFactor (varsG.mvar(), 0));
1375          }
1376        }
1377        else
1378        {
1379          if (j != 0)
1380            recurse= true;
1381          varsGMapLevel.insert (CFFactor (varsG.mvar(),exp));
1382        }
1383        varsG /= varsG.mvar();
1384      }
1385
1386      if (recurse)
1387      {
1388        ii=as;
1389        for (; ii.hasItem(); ii++)
1390        {
1391          if (ii.getItem() == i.getItem())
1392            continue;
1393          for (iter= varsGMapLevel; iter.hasItem(); iter++)
1394            ii.getItem()= inflatePoly (ii.getItem(), iter.getItem().exp(),
1395                                       iter.getItem().factor().level());
1396        }
1397      }
1398      else
1399      {
1400        ii= i;
1401        ii++;
1402        for (; ii.hasItem(); ii++)
1403        {
1404          for (iter= varsGMapLevel; iter.hasItem(); iter++)
1405          {
1406            ii.getItem()= inflatePoly (ii.getItem(), iter.getItem().exp(),
1407                                       iter.getItem().factor().level());
1408          }
1409        }
1410      }
1411      if (varsGMap[j].isEmpty())
1412        varsGMap[j]= varsGMapLevel;
1413      else
1414      {
1415        if (!varsGMapLevel.isEmpty())
1416        {
1417          tmp= varsGMap[j];
1418          CFFListIterator iter2= varsGMapLevel;
1419          ASSERT (tmp.length() == varsGMapLevel.length(), "wrong length of lists");
1420          for (iter=tmp; iter.hasItem(); iter++, iter2++)
1421            iter.getItem()= CFFactor (iter.getItem().factor(),
1422                                  iter.getItem().exp() + iter2.getItem().exp());
1423          varsGMap[j]= tmp;
1424        }
1425      }
1426      varsGMapLevel= CFFList();
1427    }
1428    asnew.append (i.getItem());
1429    if (!recurse)
1430    {
1431      i++;
1432      j++;
1433    }
1434    else
1435    {
1436      i= as;
1437      j= 0;
1438      recurse= false;
1439      asnew= CFList();
1440    }
1441  }
1442
1443  while (!lcmVars.isOne())
1444  {
1445    varsMapLevel.insert (CFFactor (lcmVars.mvar(), 0));
1446    lcmVars /= lcmVars.mvar();
1447  }
1448
1449  // compute how to map variables in F
1450  for (j= 0; j < as.length(); j++)
1451  {
1452    if (varsGMap[j].isEmpty())
1453      continue;
1454
1455    for (CFFListIterator iter2= varsGMap[j]; iter2.hasItem(); iter2++)
1456    {
1457      for (iter= varsMapLevel; iter.hasItem(); iter++)
1458      {
1459        if (iter.getItem().factor() == iter2.getItem().factor())
1460        {
1461          iter.getItem()= CFFactor (iter.getItem().factor(),
1462                                  iter.getItem().exp() + iter2.getItem().exp());
1463        }
1464      }
1465    }
1466  }
1467
1468  delete [] varsGMap;
1469
1470  if (derivZeroF)
1471  {
1472    as.removeLast();
1473    asnew.removeLast();
1474    F= deflatePoly (F, expF, F.level());
1475  }
1476
1477  // map variables in F
1478  for (iter= varsMapLevel; iter.hasItem(); iter++)
1479  {
1480    if (expF > 0)
1481      tmpExp= iter.getItem().exp() - expF;
1482    else
1483      tmpExp= iter.getItem().exp();
1484
1485    if (tmpExp > 0)
1486      F= inflatePoly (F, tmpExp, iter.getItem().factor().level());
1487    else if (tmpExp < 0)
1488      F= deflatePoly (F, -tmpExp, iter.getItem().factor().level());
1489  }
1490
1491  // factor F with minimal polys given in asnew
1492  asnew.append (F);
1493  asnew= charSetViaModCharSet (asnew); // TODO use modCharSet
1494
1495  F= asnew.getLast();
1496  F /= content (F); // should be obsolete if we use modCharSet
1497
1498  asnew.removeLast();
1499  for (CFListIterator i= asnew; i.hasItem(); i++) // should be obsolete if we use modCharSet
1500    i.getItem() /= content (i.getItem());
1501
1502  j= 0;
1503  tmp= newcfactor (F, asnew, j);
1504
1505  // transform back
1506  j= 0;
1507  int p= getCharacteristic();
1508  CFList transBack;
1509  CFMap M;
1510  CanonicalForm g;
1511
1512  for (iter= varsMapLevel; iter.hasItem(); iter++)
1513  {
1514    if (iter.getItem().exp() > 0)
1515    {
1516      j++;
1517      g= power (Variable (f.level() + j), ipower (p, iter.getItem().exp())) -
1518         iter.getItem().factor().mvar();
1519      transBack.append (g);
1520      M.newpair (iter.getItem().factor().mvar(), Variable (f.level() + j));
1521    }
1522  }
1523
1524  for (i= asnew; i.hasItem(); i++)
1525    transBack.insert (M (i.getItem()));
1526
1527  if (expF > 0)
1528    tmpExp= ipower (p, expF);
1529
1530  CFFList result;
1531  CFList transform;
1532
1533  for (CFFListIterator k= tmp; k.hasItem(); k++)
1534  {
1535    transform= transBack;
1536    CanonicalForm factor= k.getItem().factor();
1537    factor= M (factor);
1538    transform.append (factor);
1539    transform= charSetViaModCharSet (transform); //TODO use modCharSet
1540    for (i= transform; i.hasItem(); i++)
1541    {
1542      if (degree (i.getItem(), f.mvar()) > 0)
1543      {
1544        factor= i.getItem();
1545        break;
1546      }
1547    }
1548
1549    factor /= content (factor); // should be superflous if we use modCharSet instead of CharSet
1550
1551    if (expF > 0)
1552    {
1553      int mult= tmpExp/(degree (factor)/degree (k.getItem().factor()));
1554      result.append (CFFactor (factor, k.getItem().exp()*mult));
1555    }
1556    else
1557      result.append (CFFactor (factor, k.getItem().exp()));
1558  }
1559
1560  return result;
1561}
1562
1563void
1564multiplicity (CFFList& factors, const CanonicalForm& F, const CFList& as)
1565{
1566  CanonicalForm G= F;
1567  Variable x= F.mvar();
1568  CanonicalForm q, r;
1569  int count= -1;
1570  for (CFFListIterator iter=factors; iter.hasItem(); iter++)
1571  {
1572    count= -1;
1573    if (iter.getItem().factor().inCoeffDomain())
1574      continue;
1575    while (1)
1576    {
1577      psqr (G, iter.getItem().factor(), q, r, x);
1578
1579      q= Prem (q, as);
1580      r= Prem (r, as);
1581      if (!r.isZero())
1582        break;
1583      count++;
1584      G= q;
1585    }
1586    iter.getItem()= CFFactor (iter.getItem().factor(), iter.getItem().exp()+count);
1587  }
1588}
1589
1590// 1) prepares data
1591// 2) for char=p we distinguish 3 cases:
1592//           no transcendentals, seperable and inseperable extensions
1593CFFList
1594newfactoras( const CanonicalForm & f, const CFList & as, int &success)
1595{
1596  bool isRat= isOn (SW_RATIONAL);
1597  if (!isRat && getCharacteristic() == 0)
1598    On (SW_RATIONAL);
1599  Variable vf=f.mvar();
1600  CFListIterator i;
1601  CFFListIterator jj;
1602  CFList reduceresult;
1603  CFFList result;
1604
1605  success=1;
1606
1607// F1: [Test trivial cases]
1608// 1) first trivial cases:
1609  if (vf.level() <= as.getLast().level())
1610  {
1611// ||( (as.length()==1) && (degree(f,vf)==3) && (degree(as.getFirst()==2)) )
1612    if (!isRat && getCharacteristic() == 0)
1613      Off (SW_RATIONAL);
1614    return CFFList(CFFactor(f,1));
1615  }
1616
1617// 2) List of variables:
1618// 2a) Setup list of those polys in AS having degree(AS[i], AS[i].mvar()) > 1
1619// 2b) Setup variableordering
1620  CFList Astar;
1621  Variable x;
1622  CanonicalForm elem;
1623  Varlist ord, uord,oldord;
1624  for ( int ii=1; ii< level(vf) ; ii++ ) { uord.append(Variable(ii));  }
1625  oldord= uord; oldord.append(vf);
1626
1627  for ( i=as; i.hasItem(); i++ ){
1628    elem= i.getItem();
1629    x= elem.mvar();
1630    if ( degree(elem,x) > 1){ // otherwise it's not an extension
1631      Astar.append(elem);
1632      ord.append(x);
1633    }
1634  }
1635  uord= Difference(uord,ord);
1636
1637// 3) second trivial cases: we already proved irr. of f over no extensions
1638  if ( Astar.length() == 0 ){
1639    if (!isRat && getCharacteristic() == 0)
1640      Off (SW_RATIONAL);
1641    return CFFList(CFFactor(f,1));
1642  }
1643
1644// 4) Try to obtain a partial factorization using prop2 and prop3
1645//    Use with caution! We have to proof these propositions first!
1646  // Not yet implemented
1647
1648// 5) Look if elements in uord actually occure in any of the minimal
1649//    polynomials. If no element of uord occures in any of the minimal
1650//   polynomials, we don't have transzendentals.
1651  Varlist newuord=Var_is_in_AS(uord,Astar);
1652
1653  CFFList Factorlist;
1654  Varlist gcdord= Union(ord,newuord); gcdord.append(f.mvar());
1655  bool isFunctionField= (newuord.length() > 0);
1656
1657  // This is for now. we need alg_sqrfree implemented!
1658  CanonicalForm Fgcd= 0;
1659  if (isFunctionField)
1660    Fgcd= alg_gcd(f,f.deriv(),Astar);
1661
1662  bool derivZero= f.deriv().isZero();
1663  if (isFunctionField && (degree (Fgcd, f.mvar()) > 0) && !derivZero)
1664  {
1665    CanonicalForm Ggcd= divide(f, Fgcd,Astar);
1666    if (getCharacteristic() == 0)
1667    {
1668      CFFList result= newfactoras (Ggcd,as,success); //Ggcd is the squarefree part of f
1669      multiplicity (result, f, Astar);
1670      if (!isRat && getCharacteristic() == 0)
1671        Off (SW_RATIONAL);
1672      return result;
1673    }
1674
1675    Fgcd= pp(Fgcd); Ggcd= pp(Ggcd);
1676    if (!isRat && getCharacteristic() == 0)
1677      Off (SW_RATIONAL);
1678    return myUnion(newfactoras(Fgcd,as,success) , newfactoras(Ggcd,as,success));
1679  }
1680
1681  if ( getCharacteristic() > 0 )
1682  {
1683    // First look for extension!
1684    IntList degreelist;
1685    Variable vminpoly;
1686    for (i=Astar; i.hasItem(); i++){degreelist.append(degree(i.getItem()));}
1687    int extdeg= getextension(degreelist, degree(f));
1688
1689    // Now the real stuff!
1690    if ( newuord.length() == 0 ) // no transzendentals
1691    {
1692      if ( extdeg > 1 ){
1693        CanonicalForm MIPO= generate_mipo( extdeg, vminpoly);
1694        vminpoly= rootOf(MIPO);
1695      }
1696      Factorlist= alg_factor(f, Astar, vminpoly, oldord, as, isFunctionField);
1697      return Factorlist;
1698    }
1699    else if ( inseperable(Astar) > 0 || derivZero) // Look if extensions are seperable
1700    {
1701      Factorlist= SteelTrager(f, Astar, newuord);
1702      return Factorlist;
1703    }
1704    else{ // we are on the save side: Use trager
1705      if (extdeg > 1 ){
1706        CanonicalForm MIPO=generate_mipo(extdeg, vminpoly );
1707        vminpoly= rootOf(MIPO);
1708      }
1709      Factorlist= alg_factor(f, Astar, vminpoly, oldord, as, isFunctionField);
1710      return Factorlist;
1711    }
1712  }
1713  else{ // char=0 apply trager directly
1714    Variable vminpoly;
1715    Factorlist= alg_factor(f, Astar, vminpoly, oldord, as, isFunctionField);
1716    if (!isRat && getCharacteristic() == 0)
1717      Off (SW_RATIONAL);
1718    return Factorlist;
1719  }
1720
1721  return CFFList(CFFactor(f,1));
1722}
1723
1724CFFList
1725newcfactor(const CanonicalForm & f, const CFList & as, int & success )
1726{
1727  bool isRat= isOn (SW_RATIONAL);
1728  if (!isRat && getCharacteristic() == 0)
1729    On (SW_RATIONAL);
1730  CFFList Output, output, Factors= factorize(f);
1731  if (Factors.getFirst().factor().inCoeffDomain())
1732    Factors.removeFirst();
1733
1734  if ( as.length() == 0 )
1735  {
1736    if (!isRat && getCharacteristic() == 0)
1737      Off (SW_RATIONAL);
1738    success=1;
1739    return Factors;
1740  }
1741  if (f.level() <= as.getLast().level())
1742  {
1743    if (!isRat && getCharacteristic() == 0)
1744      Off (SW_RATIONAL);
1745    success=1;
1746    return Factors;
1747  }
1748
1749  success=1;
1750  for ( CFFListIterator i=Factors; i.hasItem(); i++ )
1751  {
1752    if (i.getItem().factor().level() > as.getLast().level())
1753    {
1754      output=newfactoras(i.getItem().factor(),as, success);
1755      for ( CFFListIterator j=output; j.hasItem(); j++)
1756        Output = myappend(Output,CFFactor(j.getItem().factor(),j.getItem().exp()*i.getItem().exp()));
1757    }
1758  }
1759
1760  if (!isRat && getCharacteristic() == 0)
1761    Off (SW_RATIONAL);
1762  return Output;
1763}
Note: See TracBrowser for help on using the repository browser.