source: git/factory/facAlgFunc.cc @ ae8a3a

spielwiese
Last change on this file since ae8a3a was ae8a3a, checked in by Martin Lee <martinlee84@…>, 9 years ago
fix: bug in substitution
  • Property mode set to 100644
File size: 27.0 KB
RevLine 
[611df0]1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facAlgFunc.cc
5 *
6 * This file provides functions to factorize polynomials over alg. function
7 * fields
[bbf054]8 *
9 * @note some of the code is code from libfac or derived from code from libfac.
10 * Libfac is written by M. Messollen. See also COPYING for license information
11 * and README for general information on characteristic sets.
[fea494]12 *
[611df0]13 * ABSTRACT: Descriptions can be found in B. Trager "Algebraic Factoring and
14 * Rational Function Integration" and A. Steel "Conquering Inseparability:
15 * Primary decomposition and multivariate factorization over algebraic function
16 * fields of positive characteristic"
17 *
[bbf054]18 * @author Martin Lee
[611df0]19 *
20 **/
21/*****************************************************************************/
[9d6bf4]22
[5355082]23#include "config.h"
24
25#include "cf_assert.h"
26#include "debug.h"
27
28#include "cf_generator.h"
29#include "cf_iter.h"
30#include "cf_util.h"
31#include "cf_algorithm.h"
[36dfca]32#include "templates/ftmpl_functions.h"
[5355082]33#include "cf_map.h"
34#include "cfModResultant.h"
35#include "cfCharSets.h"
36#include "facAlgFunc.h"
[c514f7]37#include "facAlgFuncUtil.h"
[0479e09]38
[070c1b4]39void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
[639047e]40
[5355082]41CanonicalForm
42alg_content (const CanonicalForm& f, const CFList& as)
43{
44  if (!f.inCoeffDomain())
45  {
46    CFIterator i= f;
47    CanonicalForm result= abs (i.coeff());
48    i++;
49    while (i.hasTerms() && !result.isOne())
50    {
51      result= alg_gcd (i.coeff(), result, as);
52      i++;
53    }
54    return result;
55  }
56
57  return abs (f);
58}
59
[b17fa5]60CanonicalForm
61alg_gcd (const CanonicalForm & fff, const CanonicalForm &ggg, const CFList &as)
[5355082]62{
63  if (fff.inCoeffDomain() || ggg.inCoeffDomain())
64    return 1;
65  CanonicalForm f=fff;
66  CanonicalForm g=ggg;
67  f=Prem(f,as);
68  g=Prem(g,as);
69  if ( f.isZero() )
70  {
71    if ( g.lc().sign() < 0 ) return -g;
72    else                     return g;
73  }
74  else  if ( g.isZero() )
75  {
76    if ( f.lc().sign() < 0 ) return -f;
77    else                     return f;
78  }
79
80  int v= as.getLast().level();
81  if (f.level() <= v || g.level() <= v)
82    return 1;
83
84  CanonicalForm res;
85
86  // does as appear in f and g ?
87  bool has_alg_var=false;
88  for ( CFListIterator j=as;j.hasItem(); j++ )
89  {
90    Variable v=j.getItem().mvar();
91    if (hasVar (f, v))
92      has_alg_var=true;
93    if (hasVar (g, v))
94      has_alg_var=true;
95  }
96  if (!has_alg_var)
97  {
98    /*if ((hasAlgVar(f))
99    || (hasAlgVar(g)))
100    {
101      Varlist ord;
102      for ( CFListIterator j=as;j.hasItem(); j++ )
103        ord.append(j.getItem().mvar());
104      res=algcd(f,g,as,ord);
105    }
106    else*/
107    if (!hasAlgVar (f) && !hasAlgVar (g))
[4758c9]108      return res=gcd(f,g);
[5355082]109  }
110
111  int mvf=f.level();
112  int mvg=g.level();
113  if (mvg > mvf)
114  {
115    CanonicalForm tmp=f; f=g; g=tmp;
116    int tmp2=mvf; mvf=mvg; mvg=tmp2;
117  }
118  if (g.inBaseDomain() || f.inBaseDomain())
119    return CanonicalForm(1);
120
121  CanonicalForm c_f= alg_content (f, as);
122
123  if (mvf != mvg)
124  {
125    res= alg_gcd (g, c_f, as);
126    return res;
127  }
128  Variable x= f.mvar();
129
130  // now: mvf==mvg, f.level()==g.level()
131  CanonicalForm c_g= alg_content (g, as);
132
133  int delta= degree (f) - degree (g);
134
135  f= divide (f, c_f, as);
136  g= divide (g, c_g, as);
137
138  // gcd of contents
139  CanonicalForm c_gcd= alg_gcd (c_f, c_g, as);
140  CanonicalForm tmp;
141
142  if (delta < 0)
143  {
144    tmp= f;
145    f= g;
146    g= tmp;
147    delta= -delta;
148  }
149
150  CanonicalForm r= 1;
151
152  while (degree (g, x) > 0)
153  {
[c514f7]154    r= Prem (f, g);
[5355082]155    r= Prem (r, as);
156    if (!r.isZero())
157    {
158      r= divide (r, alg_content (r,as), as);
159      r /= vcontent (r,Variable (v+1));
160    }
161    f= g;
162    g= r;
163  }
164
165  if (degree (g, x) == 0)
166    return c_gcd;
167
168  c_f= alg_content (f, as);
169
170  f= divide (f, c_f, as);
171
172  f *= c_gcd;
173  f /= vcontent (f, Variable (v+1));
174
175  return f;
176}
177
[0479e09]178static CanonicalForm
[ea83e2]179resultante (const CanonicalForm & f, const CanonicalForm& g, const Variable & v)
[9d6bf4]180{
[625105]181  bool on_rational = isOn(SW_RATIONAL);
[be5bf9]182  if (!on_rational && getCharacteristic() == 0)
183    On(SW_RATIONAL);
[625105]184  CanonicalForm cd = bCommonDen( f );
185  CanonicalForm fz = f * cd;
186  cd = bCommonDen( g );
187  CanonicalForm gz = g * cd;
[be5bf9]188  if (!on_rational && getCharacteristic() == 0)
189    Off(SW_RATIONAL);
[a199a3]190  CanonicalForm result;
[1936fb]191#ifdef HAVE_NTL
[a199a3]192  if (getCharacteristic() == 0)
193    result= resultantZ (fz, gz,v);
194  else
[1936fb]195#endif
[a199a3]196    result= resultant (fz,gz,v);
[625105]197
[a199a3]198  return result;
[0479e09]199}
200
[0ac77c4]201/// compute the norm R of f over PPalpha, g= f (x-s*alpha)
202/// if proof==true, R is squarefree and if in addition getCharacteristic() > 0
203/// the squarefree factors of R are returned.
204/// Based on Trager's sqrf_norm algorithm.
[868240]205static CFFList
[0ac77c4]206norm (const CanonicalForm & f, const CanonicalForm & PPalpha,
207      CFGenerator & myrandom, CanonicalForm & s, CanonicalForm & g,
208      CanonicalForm & R, bool proof)
[9d6bf4]209{
[c9f357]210  Variable y= PPalpha.mvar(), vf= f.mvar();
[868240]211  CanonicalForm temp, Palpha= PPalpha, t;
212  int sqfreetest= 0;
[0479e09]213  CFFList testlist;
214  CFFListIterator i;
[0be2bc]215
[c9f357]216  if (proof)
217  {
218    myrandom.reset();
219    s= myrandom.item();
220    g= f;
221    R= CanonicalForm(0);
222  }
223  else
224  {
225    if (getCharacteristic() == 0)
226      t= CanonicalForm (mapinto (myrandom.item()));
227    else
228      t= CanonicalForm (myrandom.item());
229    s= t;
230    g= f (vf - t*Palpha.mvar(), vf);
231  }
[0479e09]232
233  // Norm, resultante taken with respect to y
[868240]234  while (!sqfreetest)
[9d6bf4]235  {
[c9f357]236    R= resultante (Palpha, g, y);
[868240]237    R= R* bCommonDen(R);
[88e95f]238    R /= content (R);
[c9f357]239    if (proof)
[405ebc]240    {
[c9f357]241      // sqfree check ; R is a polynomial in K[x]
242      if (getCharacteristic() == 0)
243      {
244        temp= gcd (R, R.deriv (vf));
245        if (degree(temp,vf) != 0 || temp == temp.genZero())
246          sqfreetest= 0;
247        else
248          sqfreetest= 1;
249      }
[ba6139]250      else
251      {
[c9f357]252        Variable X;
253        testlist= sqrFree (R);
254
255        if (testlist.getFirst().factor().inCoeffDomain())
256          testlist.removeFirst();
257        sqfreetest= 1;
258        for (i= testlist; i.hasItem(); i++)
[ba6139]259        {
[c9f357]260          if (i.getItem().exp() > 1 && degree (i.getItem().factor(),R.mvar()) > 0)
261          {
262            sqfreetest= 0;
263            break;
264          }
[ba6139]265        }
266      }
[c9f357]267      if (!sqfreetest)
268      {
269        myrandom.next();
270        if (getCharacteristic() == 0)
271          t= CanonicalForm (mapinto (myrandom.item()));
272        else
273          t= CanonicalForm (myrandom.item());
274        s= t;
275        g= f (vf - t*Palpha.mvar(), vf);
276      }
[b96e07]277    }
[c9f357]278    else
279      break;
[b96e07]280  }
[868240]281  return testlist;
[b96e07]282}
[0479e09]283
[0ac77c4]284/// see @a norm, R is guaranteed to be squarefree
285/// Based on Trager's sqrf_norm algorithm.
[868240]286static CFFList
[0ac77c4]287sqrfNorm (const CanonicalForm & f, const CanonicalForm & PPalpha,
288          const Variable & Extension, CanonicalForm & s,  CanonicalForm & g,
289          CanonicalForm & R)
[9d6bf4]290{
[868240]291  CFFList result;
[3ff2cf]292  if (getCharacteristic() == 0)
[9d6bf4]293  {
[0be2bc]294    IntGenerator myrandom;
[0ac77c4]295    result= norm (f, PPalpha, myrandom, s, g, R, true);
[0be2bc]296  }
[3ff2cf]297  else if (degree (Extension) > 0)
[9d6bf4]298  {
[868240]299    AlgExtGenerator myrandom (Extension);
[0ac77c4]300    result= norm (f, PPalpha, myrandom, s, g, R, true);
[0479e09]301  }
[9d6bf4]302  else
303  {
[90e403f]304    FFGenerator myrandom;
[0ac77c4]305    result= norm (f, PPalpha, myrandom, s, g, R, true);
[0479e09]306  }
[868240]307  return result;
[0479e09]308}
309
310// calculate a "primitive element"
[3ff2cf]311// K must have more than S elements (-> getDegOfExt)
[0be2bc]312static CFList
[ea83e2]313simpleExtension (CFList& backSubst, const CFList & Astar,
314                 const Variable & Extension, bool& isFunctionField,
315                 CanonicalForm & R)
[dbcaa4]316{
[abde36]317  CFList Returnlist, Bstar= Astar;
318  CanonicalForm s, g, ra, rb, oldR, h, denra, denrb= 1;
[dbcaa4]319  Variable alpha;
320  CFList tmp;
321
322  bool isRat= isOn (SW_RATIONAL);
[0479e09]323
[dbcaa4]324  CFListIterator j;
325  if (Astar.length() == 1)
326  {
327    R= Astar.getFirst();
328    rb= R.mvar();
[ae8a3a]329    Returnlist.append (rb);
330    if (isFunctionField)
331      Returnlist.append (denrb);
[dbcaa4]332  }
333  else
334  {
335    R=Bstar.getFirst();
336    Bstar.removeFirst();
[3ff2cf]337    for (CFListIterator i= Bstar; i.hasItem(); i++)
[dbcaa4]338    {
339      j= i;
340      j++;
[be5bf9]341      if (getCharacteristic() == 0)
342        Off (SW_RATIONAL);
[dbcaa4]343      R /= icontent (R);
[be5bf9]344      if (getCharacteristic() == 0)
345        On (SW_RATIONAL);
[dbcaa4]346      oldR= R;
[868240]347      //TODO normalize i.getItem over K(R)?
[0ac77c4]348      (void) sqrfNorm (i.getItem(), R, Extension, s, g, R);
[dbcaa4]349
[ea0a9d]350      backSubst.insert (s);
351
[be5bf9]352      if (getCharacteristic() == 0)
353        Off (SW_RATIONAL);
[dbcaa4]354      R /= icontent (R);
355
[be5bf9]356      if (getCharacteristic() == 0)
357        On (SW_RATIONAL);
[dbcaa4]358
359      if (!isFunctionField)
360      {
361        alpha= rootOf (R);
362        h= replacevar (g, g.mvar(), alpha);
[be5bf9]363        if (getCharacteristic() == 0)
364          On (SW_RATIONAL); //needed for GCD
[dbcaa4]365        h= gcd (h, oldR);
366        h /= Lc (h);
367        ra= -h[0];
368        ra= replacevar(ra, alpha, g.mvar());
369        rb= R.mvar()-s*ra;
370        for (; j.hasItem(); j++)
371        {
372          j.getItem()= j.getItem() (rb, i.getItem().mvar());
[4ab95b6]373          j.getItem()= j.getItem() (ra, oldR.mvar());
[dbcaa4]374        }
[03f640]375        prune (alpha);
[dbcaa4]376      }
377      else
378      {
[be5bf9]379        if (getCharacteristic() == 0)
380          On (SW_RATIONAL);
[36dfca]381        Variable v= Variable (tmax (g.level(), oldR.level()) + 1);
382        h= swapvar (g, oldR.mvar(), v);
383        tmp= CFList (R);
384        h= alg_gcd (h, swapvar (oldR, oldR.mvar(), v), tmp);
[dbcaa4]385
[ac2c30]386        CanonicalForm numinv, deninv;
387        numinv= QuasiInverse (tmp.getFirst(), LC (h), tmp.getFirst().mvar());
388        h *= numinv;
[6ae00c]389        h= Prem (h, tmp);
[ac2c30]390        deninv= LC(h);
[dbcaa4]391
392        ra= -h[0];
[ac2c30]393        denra= gcd (ra, deninv);
[dbcaa4]394        ra /= denra;
[ac2c30]395        denra= deninv/denra;
[dbcaa4]396        rb= R.mvar()*denra-s*ra;
397        denrb= denra;
398        for (; j.hasItem(); j++)
399        {
[ea83e2]400          CanonicalForm powdenra= power (denra, degree (j.getItem(),
[4ab95b6]401                                         i.getItem().mvar()));
[ea83e2]402          j.getItem()= evaluate (j.getItem(), rb, denrb, powdenra,
403                                 i.getItem().mvar());
[4ab95b6]404          powdenra= power (denra, degree (j.getItem(), oldR.mvar()));
405          j.getItem()= evaluate (j.getItem(),ra, denra, powdenra, oldR.mvar());
406
[dbcaa4]407        }
408      }
409
410      Returnlist.append (ra);
411      if (isFunctionField)
412        Returnlist.append (denra);
[ae8a3a]413      Returnlist.append (rb);
414      if (isFunctionField)
415        Returnlist.append (denrb);
[0479e09]416    }
417  }
[dbcaa4]418
[be5bf9]419  if (isRat && getCharacteristic() == 0)
[dbcaa4]420    On (SW_RATIONAL);
[be5bf9]421  else if (!isRat && getCharacteristic() == 0)
[dbcaa4]422    Off (SW_RATIONAL);
[0479e09]423
424  return Returnlist;
425}
426
[c514f7]427/// Trager's algorithm, i.e. convert to one field extension and
428/// factorize over this field extension
[0479e09]429static CFFList
[c514f7]430Trager (const CanonicalForm & F, const CFList & Astar,
[2ed604a]431        const Variable & vminpoly, const CFList & as, bool isFunctionField)
[38e7b3]432{
[be5bf9]433  bool isRat= isOn (SW_RATIONAL);
[9e74d36]434  CFFList L, normFactors, tmp;
[868240]435  CFFListIterator iter, iter2;
[0aa5d5]436  CanonicalForm R, Rstar, s, g, h, f= F;
[ea0a9d]437  CFList substlist, backSubsts;
[0479e09]438
[ea83e2]439  substlist= simpleExtension (backSubsts, Astar, vminpoly, isFunctionField,
440                              Rstar);
[0be2bc]441
[0aa5d5]442  f= subst (f, Astar, substlist, Rstar, isFunctionField);
443
444  Variable alpha;
445  if (!isFunctionField)
446  {
447    alpha= rootOf (Rstar);
448    g= replacevar (f, Rstar.mvar(), alpha);
449
[be5bf9]450    if (!isRat && getCharacteristic() == 0)
451      On (SW_RATIONAL);
[9e74d36]452    tmp= factorize (g, alpha); // factorization over number field
[0aa5d5]453
[9e74d36]454    for (iter= tmp; iter.hasItem(); iter++)
[0aa5d5]455    {
[868240]456      h= iter.getItem().factor();
[0aa5d5]457      if (!h.inCoeffDomain())
458      {
459        h= replacevar (h, alpha, Rstar.mvar());
460        h *= bCommonDen(h);
[ea0a9d]461        h= backSubst (h, backSubsts, Astar);
[0aa5d5]462        h= Prem (h, as);
[868240]463        L.append (CFFactor (h, iter.getItem().exp()));
[0aa5d5]464      }
465    }
[be5bf9]466    if (!isRat && getCharacteristic() == 0)
467      Off (SW_RATIONAL);
[03f640]468    prune (alpha);
[0aa5d5]469    return L;
470  }
[956535]471  // after here we are over an extension of a function field
472
473  // make quasi monic
474  CFList Rstarlist= CFList (Rstar);
[36aaf8]475  int algExtLevel= Astar.getLast().level(); //highest level of algebraic variables
[ac2c30]476  CanonicalForm numinv;
[be5bf9]477  if (!isRat && getCharacteristic() == 0)
478    On (SW_RATIONAL);
[abde36]479  numinv= QuasiInverse (Rstar, alg_LC (f, algExtLevel), Rstar.mvar());
[956535]480
[ac2c30]481  f *= numinv;
[956535]482  f= Prem (f, Rstarlist);
483  f /= vcontent (f, Rstar.mvar());
484  // end quasi monic
[0aa5d5]485
[e031fa1]486  if (degree (f) == 1)
[204e3c]487  {
488    f= backSubst (f, backSubsts, Astar);
489    f *= bCommonDen (f);
490    f= Prem (f, as);
491    f /= vcontent (f, as.getFirst().mvar());
492
[e031fa1]493    return CFFList (CFFactor (f, 1));
[204e3c]494  }
[e031fa1]495
[86acc34]496  CFGenerator * Gen;
497  if (getCharacteristic() == 0)
[67294f1]498    Gen= CFGenFactory::generate();
[86acc34]499  else if (degree (vminpoly) > 0)
500    Gen= AlgExtGenerator (vminpoly).clone();
[0be2bc]501  else
[67294f1]502    Gen= CFGenFactory::generate();
[115a9a]503
[86acc34]504  CFFList LL= CFFList (CFFactor (f, 1));
[204e3c]505
[86acc34]506  Variable X;
507  do
[38e7b3]508  {
[86acc34]509    tmp= CFFList();
510    for (iter= LL; iter.hasItem(); iter++)
[38e7b3]511    {
[86acc34]512      f= iter.getItem().factor();
[9e74d36]513      (void) norm (f, Rstar, *Gen, s, g, R, false);
[86acc34]514
515      if (hasFirstAlgVar (R, X))
[9e74d36]516        normFactors= factorize (R, X);
[574650]517      else
[9e74d36]518        normFactors= factorize (R);
[86acc34]519
[9e74d36]520      if (normFactors.getFirst().factor().inCoeffDomain())
521        normFactors.removeFirst();
[6ae00c]522      if (normFactors.length() < 1 || (normFactors.length() == 1 && normFactors.getLast().exp() == 1))
[574650]523      {
[86acc34]524        f= backSubst (f, backSubsts, Astar);
525        f *= bCommonDen (f);
526        f= Prem (f, as);
527        f /= vcontent (f, as.getFirst().mvar());
528
529        L.append (CFFactor (f, 1));
[0479e09]530      }
[86acc34]531      else
532      {
533        g= f;
[9e74d36]534        int count= 0;
535        for (iter2= normFactors; iter2.hasItem(); iter2++)
[86acc34]536        {
537          CanonicalForm fnew= iter2.getItem().factor();
[9e74d36]538          if (fnew.level() <= Rstar.level()) //factor is a constant from the function field
[86acc34]539            continue;
540          else
541          {
542            fnew= fnew (g.mvar() + s*Rstar.mvar(), g.mvar());
[6ae00c]543            fnew= Prem (fnew, CFList (Rstar));
[86acc34]544          }
[574650]545
[86acc34]546          h= alg_gcd (g, fnew, Rstarlist);
547          numinv= QuasiInverse (Rstar, alg_LC (h, algExtLevel), Rstar.mvar());
548          h *= numinv;
549          h= Prem (h, Rstarlist);
550          h /= vcontent (h, Rstar.mvar());
[574650]551
[9e74d36]552          if (h.level() > Rstar.level())
[86acc34]553          {
554            g= divide (g, h, Rstarlist);
555            if (degree (h) == 1 || iter2.getItem().exp() == 1)
556            {
557              h= backSubst (h, backSubsts, Astar);
558              h= Prem (h, as);
559              h *= bCommonDen (h);
560              h /= vcontent (h, as.getFirst().mvar());
561              L.append (CFFactor (h, 1));
562            }
563            else
564              tmp.append (CFFactor (h, iter2.getItem().exp()));
565          }
[9e74d36]566          count++;
567          if (g.level() <= Rstar.level())
568            break;
569          if (count == normFactors.length() - 1)
570          {
571            if (normFactors.getLast().exp() == 1 && g.level() > Rstar.level())
572            {
573              g= backSubst (g, backSubsts, Astar);
574              g= Prem (g, as);
575              g *= bCommonDen (g);
576              g /= vcontent (g, as.getFirst().mvar());
577              L.append (CFFactor (g, 1));
578            }
579            else if (normFactors.getLast().exp() > 1 &&
580                     g.level() > Rstar.level())
581            {
582              g /= vcontent (g, Rstar.mvar());
583              tmp.append (CFFactor (g, normFactors.getLast().exp()));
584            }
585            break;
586          }
[86acc34]587        }
[0479e09]588      }
589    }
[86acc34]590    LL= tmp;
591    (*Gen).next();
[639047e]592  }
[86acc34]593  while (!LL.isEmpty());
[5355082]594
[be5bf9]595  if (!isRat && getCharacteristic() == 0)
596    Off (SW_RATIONAL);
597
[86acc34]598  delete Gen;
599
600  return L;
[0479e09]601}
602
[23f3c87]603
[73e699]604/// map elements in @a AS into a PIE \f$ L \f$ and record where the variables
605/// are mapped to in @a varsMapLevel, i.e @a varsMapLevel contains a list of
606/// pairs of variables \f$ v_i \f$ and integers \f$ e_i \f$ such that
607/// \f$ L=K(\sqrt[p^{e_1}]{v_1}, \ldots, \sqrt[p^{e_n}]{v_n}) \f$
[23f3c87]608CFList
609mapIntoPIE (CFFList& varsMapLevel, CanonicalForm& lcmVars, const CFList & AS)
[9d6bf4]610{
[23f3c87]611  CanonicalForm varsG;
612  int j, exp= 0, tmpExp;
613  bool recurse= false;
[ebf7d8]614  CFList asnew, as= AS;
[23f3c87]615  CFListIterator i= as, ii;
616  CFFList varsGMapLevel, tmp;
[ebf7d8]617  CFFListIterator iter;
618  CFFList * varsGMap= new CFFList [as.length()];
619  for (j= 0; j < as.length(); j++)
620    varsGMap[j]= CFFList();
621  j= 0;
622  while (i.hasItem())
[9d6bf4]623  {
[ebf7d8]624    if (i.getItem().deriv() == 0)
[9d6bf4]625    {
[ebf7d8]626      deflateDegree (i.getItem(), exp, i.getItem().level());
627      i.getItem()= deflatePoly (i.getItem(), exp, i.getItem().level());
[23f3c87]628
[ebf7d8]629      varsG= getVars (i.getItem());
630      varsG /= i.getItem().mvar();
631
632      lcmVars= lcm (varsG, lcmVars);
633
634      while (!varsG.isOne())
635      {
636        if (i.getItem().deriv (varsG.level()).isZero())
637        {
638          deflateDegree (i.getItem(), tmpExp, varsG.level());
639          if (exp >= tmpExp)
640          {
641            if (exp == tmpExp)
642              i.getItem()= deflatePoly (i.getItem(), exp, varsG.level());
643            else
644            {
645              if (j != 0)
646                recurse= true;
647              i.getItem()= deflatePoly (i.getItem(), tmpExp, varsG.level());
648            }
649            varsGMapLevel.insert (CFFactor (varsG.mvar(), exp - tmpExp));
650          }
651          else
652          {
653            i.getItem()= deflatePoly (i.getItem(), exp, varsG.level());
654            varsGMapLevel.insert (CFFactor (varsG.mvar(), 0));
655          }
656        }
657        else
658        {
659          if (j != 0)
660            recurse= true;
661          varsGMapLevel.insert (CFFactor (varsG.mvar(),exp));
662        }
663        varsG /= varsG.mvar();
664      }
665
666      if (recurse)
[9d6bf4]667      {
[868240]668        ii= as;
[ebf7d8]669        for (; ii.hasItem(); ii++)
670        {
671          if (ii.getItem() == i.getItem())
672            continue;
673          for (iter= varsGMapLevel; iter.hasItem(); iter++)
674            ii.getItem()= inflatePoly (ii.getItem(), iter.getItem().exp(),
675                                       iter.getItem().factor().level());
676        }
[0479e09]677      }
[ebf7d8]678      else
679      {
680        ii= i;
681        ii++;
682        for (; ii.hasItem(); ii++)
683        {
684          for (iter= varsGMapLevel; iter.hasItem(); iter++)
685          {
686            ii.getItem()= inflatePoly (ii.getItem(), iter.getItem().exp(),
687                                       iter.getItem().factor().level());
688          }
689        }
690      }
691      if (varsGMap[j].isEmpty())
692        varsGMap[j]= varsGMapLevel;
693      else
694      {
695        if (!varsGMapLevel.isEmpty())
696        {
697          tmp= varsGMap[j];
698          CFFListIterator iter2= varsGMapLevel;
[ea83e2]699          ASSERT (tmp.length() == varsGMapLevel.length(),
700                  "wrong length of lists");
[ebf7d8]701          for (iter=tmp; iter.hasItem(); iter++, iter2++)
702            iter.getItem()= CFFactor (iter.getItem().factor(),
703                                  iter.getItem().exp() + iter2.getItem().exp());
704          varsGMap[j]= tmp;
705        }
706      }
707      varsGMapLevel= CFFList();
708    }
709    asnew.append (i.getItem());
710    if (!recurse)
711    {
712      i++;
713      j++;
714    }
715    else
716    {
717      i= as;
718      j= 0;
719      recurse= false;
720      asnew= CFList();
721    }
722  }
723
724  while (!lcmVars.isOne())
725  {
726    varsMapLevel.insert (CFFactor (lcmVars.mvar(), 0));
727    lcmVars /= lcmVars.mvar();
728  }
729
730  for (j= 0; j < as.length(); j++)
731  {
732    if (varsGMap[j].isEmpty())
733      continue;
734
735    for (CFFListIterator iter2= varsGMap[j]; iter2.hasItem(); iter2++)
736    {
737      for (iter= varsMapLevel; iter.hasItem(); iter++)
[9d6bf4]738      {
[ebf7d8]739        if (iter.getItem().factor() == iter2.getItem().factor())
740        {
741          iter.getItem()= CFFactor (iter.getItem().factor(),
742                                  iter.getItem().exp() + iter2.getItem().exp());
[0be2bc]743        }
[0479e09]744      }
745    }
746  }
[ebf7d8]747
748  delete [] varsGMap;
749
[23f3c87]750  return asnew;
751}
752
753/// algorithm of A. Steel described in "Conquering Inseparability: Primary
754/// decomposition and multivariate factorization over algebraic function fields
755/// of positive characteristic" with the following modifications: we use
756/// characteristic sets in IdealOverSubfield and Trager's primitive element
757/// algorithm instead of FGLM
758CFFList
759SteelTrager (const CanonicalForm & f, const CFList & AS)
760{
761  CanonicalForm F= f, lcmVars= 1;
762  CFList asnew, as= AS;
763  CFListIterator i, ii;
764
765  bool derivZeroF= false;
766  int j, expF= 0, tmpExp;
767  CFFList varsMapLevel, tmp;
768  CFFListIterator iter;
769
770  // check if F is separable
771  if (F.deriv().isZero())
772  {
773    derivZeroF= true;
774    deflateDegree (F, expF, F.level());
775  }
776
777  CanonicalForm varsF= getVars (F);
778  varsF /= F.mvar();
779
780  lcmVars= lcm (varsF, lcmVars);
781
782  if (derivZeroF)
783    as.append (F);
784
785  asnew= mapIntoPIE (varsMapLevel, lcmVars, as);
786
[ebf7d8]787  if (derivZeroF)
788  {
789    asnew.removeLast();
790    F= deflatePoly (F, expF, F.level());
791  }
792
793  // map variables in F
794  for (iter= varsMapLevel; iter.hasItem(); iter++)
795  {
796    if (expF > 0)
797      tmpExp= iter.getItem().exp() - expF;
798    else
799      tmpExp= iter.getItem().exp();
800
801    if (tmpExp > 0)
802      F= inflatePoly (F, tmpExp, iter.getItem().factor().level());
803    else if (tmpExp < 0)
804      F= deflatePoly (F, -tmpExp, iter.getItem().factor().level());
805  }
806
807  // factor F with minimal polys given in asnew
808  asnew.append (F);
[e3cdb4]809  asnew= charSetViaModCharSet (asnew, false);
[ebf7d8]810
811  F= asnew.getLast();
[bbf054]812  F /= content (F);
[ebf7d8]813
814  asnew.removeLast();
[868240]815  for (i= asnew; i.hasItem(); i++)
[ebf7d8]816    i.getItem() /= content (i.getItem());
817
[c514f7]818  tmp= facAlgFunc (F, asnew);
[ebf7d8]819
820  // transform back
821  j= 0;
822  int p= getCharacteristic();
823  CFList transBack;
824  CFMap M;
825  CanonicalForm g;
826
827  for (iter= varsMapLevel; iter.hasItem(); iter++)
[9d6bf4]828  {
[ebf7d8]829    if (iter.getItem().exp() > 0)
830    {
831      j++;
832      g= power (Variable (f.level() + j), ipower (p, iter.getItem().exp())) -
833         iter.getItem().factor().mvar();
834      transBack.append (g);
835      M.newpair (iter.getItem().factor().mvar(), Variable (f.level() + j));
836    }
837  }
838
839  for (i= asnew; i.hasItem(); i++)
840    transBack.insert (M (i.getItem()));
841
842  if (expF > 0)
843    tmpExp= ipower (p, expF);
844
845  CFFList result;
846  CFList transform;
847
[da9cd4]848  bool found;
[868240]849  for (iter= tmp; iter.hasItem(); iter++)
[ebf7d8]850  {
[da9cd4]851    found= false;
[ebf7d8]852    transform= transBack;
[868240]853    CanonicalForm factor= iter.getItem().factor();
[ebf7d8]854    factor= M (factor);
855    transform.append (factor);
[cf6604]856    transform= modCharSet (transform, false);
[da9cd4]857
858retry:
859    if (transform.isEmpty())
860    {
861      transform= transBack;
862      transform.append (factor);
863      transform= charSetViaCharSetN (transform);
864    }
[ebf7d8]865    for (i= transform; i.hasItem(); i++)
866    {
867      if (degree (i.getItem(), f.mvar()) > 0)
868      {
[da9cd4]869        if (i.getItem().level() > f.level())
870          break;
871        found= true;
[ebf7d8]872        factor= i.getItem();
873        break;
874      }
875    }
876
[da9cd4]877    if (!found)
878    {
879      found= false;
880      transform= CFList();
881      goto retry;
882    }
883
[bbf054]884    factor /= content (factor);
[ebf7d8]885
886    if (expF > 0)
[9d6bf4]887    {
[868240]888      int mult= tmpExp/(degree (factor)/degree (iter.getItem().factor()));
889      result.append (CFFactor (factor, iter.getItem().exp()*mult));
[0479e09]890    }
[ebf7d8]891    else
[868240]892      result.append (CFFactor (factor, iter.getItem().exp()));
[0479e09]893  }
894
[ebf7d8]895  return result;
[0479e09]896}
897
[c514f7]898/// factorize a polynomial that is irreducible over the ground field modulo an
899/// extension given by an irreducible characteristic set
[e0fbbeb]900
[0479e09]901// 1) prepares data
[0be2bc]902// 2) for char=p we distinguish 3 cases:
[c514f7]903//           no transcendentals, separable and inseparable extensions
[0479e09]904CFFList
[c514f7]905facAlgFunc2 (const CanonicalForm & f, const CFList & as)
[9d6bf4]906{
[863b03]907  bool isRat= isOn (SW_RATIONAL);
[be5bf9]908  if (!isRat && getCharacteristic() == 0)
[863b03]909    On (SW_RATIONAL);
[0479e09]910  Variable vf=f.mvar();
911  CFListIterator i;
912  CFFListIterator jj;
913  CFList reduceresult;
914  CFFList result;
[0be2bc]915
[0479e09]916// F1: [Test trivial cases]
917// 1) first trivial cases:
[5355082]918  if (vf.level() <= as.getLast().level())
[e031fa1]919  {
[be5bf9]920    if (!isRat && getCharacteristic() == 0)
[863b03]921      Off (SW_RATIONAL);
[0479e09]922    return CFFList(CFFactor(f,1));
923  }
924
[abde36]925// 2) Setup list of those polys in AS having degree > 1
[0479e09]926  CFList Astar;
927  Variable x;
928  CanonicalForm elem;
[2ed604a]929  Varlist ord, uord;
[abde36]930  for (int ii= 1; ii < level (vf); ii++)
931    uord.append (Variable (ii));
[0479e09]932
[abde36]933  for (i= as; i.hasItem(); i++)
934  {
[0479e09]935    elem= i.getItem();
936    x= elem.mvar();
[abde36]937    if (degree (elem, x) > 1)  // otherwise it's not an extension
938    {
939      Astar.append (elem);
940      ord.append (x);
[0479e09]941    }
942  }
[abde36]943  uord= Difference (uord, ord);
[0479e09]944
[ebf7d8]945// 3) second trivial cases: we already proved irr. of f over no extensions
[abde36]946  if (Astar.length() == 0)
947  {
[be5bf9]948    if (!isRat && getCharacteristic() == 0)
[863b03]949      Off (SW_RATIONAL);
[abde36]950    return CFFList (CFFactor (f, 1));
[0479e09]951  }
952
[abde36]953// 4) Look if elements in uord actually occur in any of the minimal
[0be2bc]954//    polynomials. If no element of uord occures in any of the minimal
[abde36]955//    polynomials the field is an alg. number field not an alg. function field
[c514f7]956  Varlist newuord= varsInAs (uord, Astar);
[0479e09]957
958  CFFList Factorlist;
[abde36]959  Varlist gcdord= Union (ord, newuord);
960  gcdord.append (f.mvar());
[0aa5d5]961  bool isFunctionField= (newuord.length() > 0);
962
[abde36]963  // TODO alg_sqrfree?
[0aa5d5]964  CanonicalForm Fgcd= 0;
965  if (isFunctionField)
[c514f7]966    Fgcd= alg_gcd (f, f.deriv(), Astar);
[0aa5d5]967
[ebf7d8]968  bool derivZero= f.deriv().isZero();
969  if (isFunctionField && (degree (Fgcd, f.mvar()) > 0) && !derivZero)
[0aa5d5]970  {
[0479e09]971    CanonicalForm Ggcd= divide(f, Fgcd,Astar);
[0aa5d5]972    if (getCharacteristic() == 0)
973    {
[c514f7]974      CFFList result= facAlgFunc2 (Ggcd, as); //Ggcd is the squarefree part of f
[0aa5d5]975      multiplicity (result, f, Astar);
[be5bf9]976      if (!isRat && getCharacteristic() == 0)
[863b03]977        Off (SW_RATIONAL);
[0aa5d5]978      return result;
979    }
[5355082]980
[abde36]981    Fgcd= pp (Fgcd);
982    Ggcd= pp (Ggcd);
[be5bf9]983    if (!isRat && getCharacteristic() == 0)
[863b03]984      Off (SW_RATIONAL);
[c514f7]985    return merge (facAlgFunc2 (Fgcd, as), facAlgFunc2 (Ggcd, as));
[0479e09]986  }
[ebf7d8]987
[abde36]988  if (getCharacteristic() > 0)
[0aa5d5]989  {
[0479e09]990    IntList degreelist;
991    Variable vminpoly;
[abde36]992    for (i= Astar; i.hasItem(); i++)
993      degreelist.append (degree (i.getItem()));
[0479e09]994
[abde36]995    int extdeg= getDegOfExt (degreelist, degree (f));
996
997    if (newuord.length() == 0) // no parameters
[ebf7d8]998    {
[abde36]999      if (extdeg > 1)
1000      {
[c514f7]1001        CanonicalForm MIPO= generateMipo (extdeg);
[0be2bc]1002        vminpoly= rootOf(MIPO);
[0479e09]1003      }
[2ed604a]1004      Factorlist= Trager(f, Astar, vminpoly, as, isFunctionField);
[03f640]1005      if (extdeg > 1)
1006        prune (vminpoly);
[0479e09]1007      return Factorlist;
1008    }
[abde36]1009    else if (isInseparable(Astar) || derivZero) // inseparable case
[ebf7d8]1010    {
[abde36]1011      Factorlist= SteelTrager (f, Astar);
[ebf7d8]1012      return Factorlist;
[0479e09]1013    }
[abde36]1014    else // separable case
1015    {
1016      if (extdeg > 1)
1017      {
[c514f7]1018        CanonicalForm MIPO=generateMipo (extdeg);
[abde36]1019        vminpoly= rootOf (MIPO);
[0479e09]1020      }
[abde36]1021      Factorlist= Trager (f, Astar, vminpoly, as, isFunctionField);
[03f640]1022      if (extdeg > 1)
1023        prune (vminpoly);
[0479e09]1024      return Factorlist;
1025    }
1026  }
[abde36]1027  else // char 0
1028  {
[0479e09]1029    Variable vminpoly;
[abde36]1030    Factorlist= Trager (f, Astar, vminpoly, as, isFunctionField);
[be5bf9]1031    if (!isRat && getCharacteristic() == 0)
[863b03]1032      Off (SW_RATIONAL);
1033    return Factorlist;
[0479e09]1034  }
1035
[c514f7]1036  return CFFList (CFFactor(f,1));
[0479e09]1037}
1038
[c514f7]1039
1040/// factorize a polynomial modulo an extension given by an irreducible
1041/// characteristic set
[0479e09]1042CFFList
[c514f7]1043facAlgFunc (const CanonicalForm & f, const CFList & as)
[9d6bf4]1044{
[863b03]1045  bool isRat= isOn (SW_RATIONAL);
[be5bf9]1046  if (!isRat && getCharacteristic() == 0)
[863b03]1047    On (SW_RATIONAL);
[0aa5d5]1048  CFFList Output, output, Factors= factorize(f);
1049  if (Factors.getFirst().factor().inCoeffDomain())
1050    Factors.removeFirst();
[0479e09]1051
[abde36]1052  if (as.length() == 0)
[0aa5d5]1053  {
[be5bf9]1054    if (!isRat && getCharacteristic() == 0)
[863b03]1055      Off (SW_RATIONAL);
[0aa5d5]1056    return Factors;
1057  }
[5355082]1058  if (f.level() <= as.getLast().level())
[0aa5d5]1059  {
[be5bf9]1060    if (!isRat && getCharacteristic() == 0)
[863b03]1061      Off (SW_RATIONAL);
[0aa5d5]1062    return Factors;
1063  }
[0479e09]1064
[abde36]1065  for (CFFListIterator i=Factors; i.hasItem(); i++)
[9d6bf4]1066  {
[ebf7d8]1067    if (i.getItem().factor().level() > as.getLast().level())
1068    {
[abde36]1069      output= facAlgFunc2 (i.getItem().factor(), as);
1070      for (CFFListIterator j= output; j.hasItem(); j++)
1071        Output= append (Output, CFFactor (j.getItem().factor(),
1072                                          j.getItem().exp()*i.getItem().exp()));
[ebf7d8]1073    }
[0479e09]1074  }
[863b03]1075
[be5bf9]1076  if (!isRat && getCharacteristic() == 0)
[863b03]1077    Off (SW_RATIONAL);
[0479e09]1078  return Output;
1079}
Note: See TracBrowser for help on using the repository browser.