source: git/factory/facAlgFuncUtil.cc @ aebdf8

fieker-DuValspielwiese
Last change on this file since aebdf8 was aebdf8, checked in by Martin Lee <martinlee84@…>, 10 years ago
fix: wrong handling of SW_RATIONAL
  • Property mode set to 100644
File size: 15.4 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facAlgFuncUtil.cc
5 *
6 * This file provides utility functions to factorize polynomials over alg.
7 * function fields
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.
12 *
13 * @author Martin Lee
14 *
15 **/
16/*****************************************************************************/
17
18#include "config.h"
19
20#include "cf_assert.h"
21
22#include "canonicalform.h"
23#include "facAlgFuncUtil.h"
24#include "cfCharSetsUtil.h"
25#include "cf_random.h"
26#include "cf_irred.h"
27#include "cf_algorithm.h"
28#include "cf_util.h"
29#include "cf_iter.h"
30
31CFFList
32append (const CFFList & Inputlist, const CFFactor & TheFactor)
33{
34  CFFList Outputlist;
35  CFFactor copy;
36  CFFListIterator i;
37  int exp=0;
38
39  for (i= Inputlist; i.hasItem() ; i++)
40  {
41    copy= i.getItem();
42    if (copy.factor() == TheFactor.factor())
43      exp += copy.exp();
44    else
45      Outputlist.append(copy);
46  }
47  Outputlist.append (CFFactor (TheFactor.factor(), exp + TheFactor.exp()));
48  return Outputlist;
49}
50
51CFFList
52merge (const CFFList & Inputlist1, const CFFList & Inputlist2)
53{
54  CFFList Outputlist;
55  CFFListIterator i;
56
57  for (i= Inputlist1; i.hasItem(); i++)
58    Outputlist= append (Outputlist, i.getItem());
59  for (i= Inputlist2; i.hasItem() ; i++)
60    Outputlist= append (Outputlist, i.getItem());
61
62  return Outputlist;
63}
64
65Varlist
66varsInAs (const Varlist & uord, const CFList & Astar)
67{
68  Varlist output;
69  CanonicalForm elem;
70  Variable x;
71
72  for (VarlistIterator i= uord; i.hasItem(); i++)
73  {
74    x= i.getItem();
75    for (CFListIterator j= Astar; j.hasItem(); j++ )
76    {
77      elem= j.getItem();
78      if (degree (elem, x) > 0) // x actually occures in Astar
79      {
80        output.append (x);
81        break;
82      }
83    }
84  }
85  return output;
86}
87
88// generate an irreducible poly of degree degOfExt over F_p
89CanonicalForm
90generateMipo (int degOfExt)
91{
92#ifndef HAVE_NTL
93  FFRandom gen;
94  return find_irreducible (degOfExt, gen, Variable (1));
95#else
96  return randomIrredpoly (degOfExt, Variable (1));
97#endif
98}
99
100CanonicalForm alg_lc (const CanonicalForm & f)
101{
102  if (f.level()>0)
103  {
104    return alg_lc(f.LC());
105  }
106
107  return f;
108}
109
110CanonicalForm alg_LC (const CanonicalForm& f, int lev)
111{
112  CanonicalForm result= f;
113  while (result.level() > lev)
114    result= LC (result);
115  return result;
116}
117
118
119CanonicalForm
120subst (const CanonicalForm& f, const CFList& a, const CFList& b,
121       const CanonicalForm& Rstar, bool isFunctionField)
122{
123  if (isFunctionField)
124    ASSERT (2*a.length() == b.length(), "wrong length of lists");
125  else
126    ASSERT (a.length() == b.length(), "lists of equal length expected");
127  CFListIterator j= b;
128  CanonicalForm result= f, tmp, powj;
129  CFListIterator i= a;
130  int length= a.length();
131  int count= 0;
132  for (; i.hasItem() && j.hasItem(); i++, j++, count++)
133  {
134    if (length - count == 2)
135    {
136      if (!isFunctionField)
137      {
138        result= result (b.getLast(), a.getLast().mvar());
139        result= result (j.getItem(), i.getItem().mvar());
140        break;
141      }
142      else
143      {
144        tmp= b.getLast();
145        j++;
146        j++;
147        powj= power (tmp, degree (result, a.getLast().mvar()));
148        result= evaluate (result, j.getItem(), tmp, powj, a.getLast().mvar());
149
150        if (fdivides (powj, result, tmp))
151          result= tmp;
152
153        result /= vcontent (result, Variable (a.getLast().level() + 1));
154        j--;
155        j--;
156        tmp= j.getItem();
157        j++;
158        powj= power (j.getItem(), degree (result, i.getItem().mvar()));
159        result= evaluate (result, tmp, j.getItem(), powj, i.getItem().mvar());
160
161        if (fdivides (powj, result, tmp))
162          result= tmp;
163
164        result /= vcontent (result, Variable (i.getItem().level() + 1));
165        break;
166      }
167    }
168    else
169    {
170      if (!isFunctionField)
171        result= result (j.getItem(), i.getItem().mvar());
172      else
173      {
174        tmp= j.getItem();
175        j++;
176        powj= power (j.getItem(), degree (result, i.getItem().mvar()));
177        result= evaluate (result, tmp, j.getItem(), powj, i.getItem().mvar());
178
179        if (fdivides (powj, result, tmp))
180          result= tmp;
181
182        result /= vcontent (result, Variable (i.getItem().level() + 1));
183      }
184    }
185  }
186  result= reduce (result, Rstar);
187  result /= vcontent (result, Variable (Rstar.level() + 1));
188  return result;
189}
190
191CanonicalForm
192backSubst (const CanonicalForm& F, const CFList& a, const CFList& b)
193{
194  ASSERT (a.length() == b.length() - 1, "wrong length of lists in backSubst");
195  CanonicalForm result= F;
196  Variable tmp;
197  CFList tmp2= b;
198  tmp= tmp2.getLast().mvar();
199  tmp2.removeLast();
200  for (CFListIterator iter= a; iter.hasItem(); iter++)
201  {
202    result= result (tmp+iter.getItem()*tmp2.getLast().mvar(), tmp);
203    tmp= tmp2.getLast().mvar();
204    tmp2.removeLast();
205  }
206  return result;
207}
208
209void deflateDegree (const CanonicalForm & F, int & pExp, int n)
210{
211  if (n == 0 || n > F.level())
212  {
213    pExp= -1;
214    return;
215  }
216  if (F.level() == n)
217  {
218    ASSERT (F.deriv().isZero(), "derivative of F is not zero");
219    int termCount=0;
220    CFIterator i= F;
221    for (; i.hasTerms(); i++)
222    {
223      if (i.exp() != 0)
224        termCount++;
225    }
226
227    int j= 1;
228    i= F;
229    for (;j < termCount; j++, i++)
230      ;
231
232    int exp= i.exp();
233    int count= 0;
234    int p= getCharacteristic();
235    while ((exp >= p) && (exp != 0) && (exp % p == 0))
236    {
237      exp /= p;
238      count++;
239    }
240    pExp= count;
241  }
242  else
243  {
244    CFIterator i= F;
245    deflateDegree (i.coeff(), pExp, n);
246    i++;
247    int tmp= pExp;
248    for (; i.hasTerms(); i++)
249    {
250      deflateDegree (i.coeff(), pExp, n);
251      if (tmp == -1)
252        tmp= pExp;
253      else if (tmp != -1 && pExp != -1)
254        pExp= (pExp < tmp) ? pExp : tmp;
255      else
256        pExp= tmp;
257    }
258  }
259}
260
261CanonicalForm deflatePoly (const CanonicalForm & F, int exp)
262{
263  if (exp == 0)
264    return F;
265  int p= getCharacteristic();
266  int pToExp= ipower (p, exp);
267  Variable x=F.mvar();
268  CanonicalForm result= 0;
269  for (CFIterator i= F; i.hasTerms(); i++)
270    result += i.coeff()*power (x, i.exp()/pToExp);
271  return result;
272}
273
274CanonicalForm deflatePoly (const CanonicalForm & F, int exps, int n)
275{
276  if (n == 0 || exps <= 0 || F.level() < n)
277    return F;
278  if (F.level() == n)
279    return deflatePoly (F, exps);
280  else
281  {
282    CanonicalForm result= 0;
283    for (CFIterator i= F; i.hasTerms(); i++)
284      result += deflatePoly (i.coeff(), exps, n)*power(F.mvar(), i.exp());
285    return result;
286  }
287}
288
289CanonicalForm inflatePoly (const CanonicalForm & F, int exp)
290{
291  if (exp == 0)
292    return F;
293  int p= getCharacteristic();
294  int pToExp= ipower (p, exp);
295  Variable x=F.mvar();
296  CanonicalForm result= 0;
297  for (CFIterator i= F; i.hasTerms(); i++)
298    result += i.coeff()*power (x, i.exp()*pToExp);
299  return result;
300}
301
302CanonicalForm inflatePoly (const CanonicalForm & F, int exps, int n)
303{
304  if (n == 0 || exps <= 0 || F.level() < n)
305    return F;
306  if (F.level() == n)
307    return inflatePoly (F, exps);
308  else
309  {
310    CanonicalForm result= 0;
311    for (CFIterator i= F; i.hasTerms(); i++)
312      result += inflatePoly (i.coeff(), exps, n)*power(F.mvar(), i.exp());
313    return result;
314  }
315}
316
317void
318multiplicity (CFFList& factors, const CanonicalForm& F, const CFList& as)
319{
320  CanonicalForm G= F;
321  Variable x= F.mvar();
322  CanonicalForm q, r;
323  int count= -1;
324  for (CFFListIterator iter=factors; iter.hasItem(); iter++)
325  {
326    count= -1;
327    if (iter.getItem().factor().inCoeffDomain())
328      continue;
329    while (1)
330    {
331      psqr (G, iter.getItem().factor(), q, r, x);
332
333      q= Prem (q, as);
334      r= Prem (r, as);
335      if (!r.isZero())
336        break;
337      count++;
338      G= q;
339    }
340    iter.getItem()= CFFactor (iter.getItem().factor(),
341                              iter.getItem().exp() + count);
342  }
343}
344
345int hasAlgVar (const CanonicalForm &f, const Variable &v)
346{
347  if (f.inBaseDomain())
348    return 0;
349  if (f.inCoeffDomain())
350  {
351    if (f.mvar() == v)
352      return 1;
353    return hasAlgVar (f.LC(), v);
354  }
355  if (f.inPolyDomain())
356  {
357    if (hasAlgVar(f.LC(), v))
358      return 1;
359    for (CFIterator i= f; i.hasTerms(); i++)
360    {
361      if (hasAlgVar (i.coeff(), v))
362        return 1;
363    }
364  }
365  return 0;
366}
367
368int hasVar (const CanonicalForm &f, const Variable &v)
369{
370  if (f.inBaseDomain())
371    return 0;
372  if (f.inCoeffDomain())
373  {
374    if (f.mvar() == v)
375      return 1;
376    return hasAlgVar (f.LC(), v);
377  }
378  if (f.inPolyDomain())
379  {
380    if (f.mvar() == v)
381      return 1;
382    if (hasVar (f.LC(), v))
383      return 1;
384    for (CFIterator i= f; i.hasTerms(); i++)
385    {
386      if (hasVar (i.coeff(), v))
387        return 1;
388    }
389  }
390  return 0;
391}
392
393int hasAlgVar (const CanonicalForm &f)
394{
395  if (f.inBaseDomain())
396    return 0;
397  if (f.inCoeffDomain())
398  {
399    if (f.level() != 0)
400      return 1;
401    return hasAlgVar(f.LC());
402  }
403  if (f.inPolyDomain())
404  {
405    if (hasAlgVar(f.LC()))
406      return 1;
407    for (CFIterator i= f; i.hasTerms(); i++)
408    {
409      if (hasAlgVar (i.coeff()))
410        return 1;
411    }
412  }
413  return 0;
414}
415
416/// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r
417void
418psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, 
419      CanonicalForm & r, CanonicalForm& multiplier, const Variable& x)
420{
421    ASSERT( x.level() > 0, "type error: polynomial variable expected" );
422    ASSERT( ! g.isZero(), "math error: division by zero" );
423
424    // swap variables such that x's level is larger or equal
425    // than both f's and g's levels.
426    Variable X;
427    if (f.level() > g.level())
428      X= f.mvar();
429    else
430      X= g.mvar();
431    if (X.level() < x.level())
432      X= x;
433    CanonicalForm F= swapvar (f, x, X);
434    CanonicalForm G= swapvar (g, x, X);
435
436    // now, we have to calculate the pseudo remainder of F and G
437    // w.r.t. X
438    int fDegree= degree (F, X);
439    int gDegree= degree (G, X);
440    if (fDegree < 0 || fDegree < gDegree)
441    {
442      q= 0;
443      r= f;
444    }
445    else
446    {
447      CanonicalForm LCG= LC (G, X);
448      multiplier= power (LCG, fDegree - gDegree + 1);
449      divrem (multiplier*F, G, q, r);
450      q= swapvar (q, x, X);
451      r= swapvar (r, x, X);
452    }
453}
454
455CanonicalForm
456Sprem (const CanonicalForm &f, const CanonicalForm &g, CanonicalForm & m,
457       CanonicalForm & q )
458{
459  CanonicalForm ff, gg, l, test, retvalue;
460  int df, dg,n;
461  bool reord;
462  Variable vf, vg, v;
463
464  if ((vf = f.mvar()) < (vg = g.mvar()))
465  {
466    m= 0;
467    q= 0;
468    return f;
469  }
470  else
471  {
472    if ( vf == vg )
473    {
474      ff= f;
475      gg= g;
476      reord= false;
477      v= vg; // == x
478    }
479    else
480    {
481      v= Variable (f.level() + 1);
482      ff= swapvar (f, vg, v); // == r
483      gg= swapvar (g, vg, v); // == v
484      reord=true;
485    }
486    dg= degree (gg, v); // == dv
487    df= degree (ff, v); // == dr
488    if (dg <= df)
489    {
490      l= LC (gg);
491      gg= gg - LC(gg)*power(v,dg);
492    }
493    else
494      l = 1;
495    n= 0;
496    while ((dg <= df) && (!ff.isZero()))
497    {
498      test= gg*LC (ff)*power (v, df - dg);
499      if (df == 0)
500        ff= 0;
501      else
502        ff= ff - LC(ff)*power(v,df);
503      ff= l*ff - test;
504      df= degree (ff, v);
505      n++;
506    }
507
508    if (reord)
509      retvalue= swapvar (ff, vg, v);
510    else
511      retvalue= ff;
512
513    m= power (l, n);
514    if (fdivides (g, m*f - retvalue))
515      q= (m*f - retvalue)/g;
516    else
517      q= 0;
518    return retvalue;
519  }
520}
521
522CanonicalForm
523divide (const CanonicalForm & ff, const CanonicalForm & f, const CFList & as)
524{
525  CanonicalForm r, m, q;
526
527  if (f.inCoeffDomain())
528  {
529    bool isRat= isOn(SW_RATIONAL);
530    if (getCharacteristic() == 0)
531      On(SW_RATIONAL);
532    q= ff/f;
533    if (!isRat && getCharacteristic() == 0)
534      Off(SW_RATIONAL);
535  }
536  else
537    r= Sprem (ff, f, m, q);
538
539  r= Prem (q, as);
540  return r;
541}
542
543// check if polynomials in Astar define a separable extension
544bool
545isInseparable (const CFList & Astar)
546{
547  CanonicalForm elem;
548
549  if (Astar.length() == 0)
550    return false;
551  for (CFListIterator i= Astar; i.hasItem(); i++)
552  {
553    elem= i.getItem();
554    if (elem.deriv().isZero())
555      return true;
556  }
557  return false;
558}
559
560// calculate big enough extension for finite fields
561// Idea: first calculate k, such that q^k > S (->thesis)
562// Second, search k with gcd(k,m_i)=1, where m_i is the degree of the i'th
563// minimal polynomial. Then the minpoly f_i remains irrd. over q^k and we
564// have enough elements to plug in.
565int
566getDegOfExt (IntList & degreelist, int n)
567{
568  int charac= getCharacteristic();
569  setCharacteristic(0); // need it for k !
570  int k= 1, m= 1, length= degreelist.length();
571  IntListIterator i;
572
573  for (i= degreelist; i.hasItem(); i++)
574     m= m*i.getItem();
575  int q= charac;
576  while (q <= ((n*m)*(n*m)/2))
577  {
578    k= k+1;
579    q= q*charac;
580  }
581  int l= 0;
582  do
583  {
584    for (i= degreelist; i.hasItem(); i++)
585    {
586      l= l + 1;
587      if (igcd (k, i.getItem()) == 1)
588      {
589        if (l == length)
590        {
591          setCharacteristic (charac);
592          return k;
593        }
594      }
595      else
596        break;
597    }
598    k= k + 1;
599    l= 0;
600  }
601  while (1);
602}
603
604CanonicalForm
605QuasiInverse (const CanonicalForm& f, const CanonicalForm& g,
606              const Variable& x)
607{
608  CanonicalForm pi, pi1, q, t0, t1, Hi, bi, pi2;
609  bool isRat= isOn (SW_RATIONAL);
610  pi= f;
611  pi1= g;
612  if (isRat)
613  {
614    pi *= bCommonDen (pi);
615    pi1 *= bCommonDen (pi1);
616  }
617  CanonicalForm m,tmp;
618  if (isRat && getCharacteristic() == 0)
619    Off (SW_RATIONAL);
620
621  pi= pi/content (pi,x);
622  pi1= pi1/content (pi1,x);
623
624  t0= 0;
625  t1= 1;
626  bi= 1;
627
628  int delta= degree (f, x) - degree (g, x);
629  Hi= power (LC (pi1, x), delta);
630  if ( (delta+1) % 2 )
631      bi = 1;
632  else
633      bi = -1;
634
635  while (degree (pi1,x) > 0)
636  {
637    psqr (pi, pi1, q, pi2, m, x);
638    pi2 /= bi;
639
640    tmp= t1;
641    t1= t0*m - t1*q;
642    t0= tmp;
643    t1 /= bi;
644    pi= pi1;
645    pi1= pi2;
646    if (degree (pi1, x) > 0)
647    {
648      delta= degree (pi, x) - degree (pi1, x);
649      if ((delta + 1) % 2)
650        bi= LC (pi, x)*power (Hi, delta);
651      else
652        bi= -LC (pi, x)*power (Hi, delta);
653      Hi= power (LC (pi1, x), delta)/power (Hi, delta - 1);
654    }
655  }
656  t1 /= gcd (pi1, t1);
657  if (isRat && getCharacteristic() == 0)
658    On (SW_RATIONAL);
659  return t1;
660}
661
662CanonicalForm
663evaluate (const CanonicalForm& f, const CanonicalForm& g,
664          const CanonicalForm& h, const CanonicalForm& powH)
665{
666  if (f.inCoeffDomain())
667    return f;
668  CFIterator i= f;
669  int lastExp = i.exp();
670  CanonicalForm result = i.coeff()*powH;
671  i++;
672  while (i.hasTerms())
673  {
674    int i_exp= i.exp();
675    if ((lastExp - i_exp) == 1)
676    {
677      result *= g;
678      result /= h;
679    }
680    else
681    {
682      result *= power (g, lastExp - i_exp);
683      result /= power (h, lastExp - i_exp);
684    }
685    result += i.coeff()*powH;
686    lastExp = i_exp;
687    i++;
688  }
689  if (lastExp != 0)
690  {
691    result *= power (g, lastExp);
692    result /= power (h, lastExp);
693  }
694  return result;
695}
696
697
698/// evaluate f at g/h at v such that powH*f is integral i.e. powH is assumed to be h^degree(f,v)
699CanonicalForm
700evaluate (const CanonicalForm& f, const CanonicalForm& g,
701          const CanonicalForm& h, const CanonicalForm& powH,
702          const Variable& v)
703{
704  if (f.inCoeffDomain())
705  {
706    return f*powH;
707  }
708
709  Variable x = f.mvar();
710  if ( v > x )
711    return f*powH;
712  else  if ( v == x )
713    return evaluate (f, g, h, powH);
714
715  // v is less than main variable of f
716  CanonicalForm result= 0;
717  for (CFIterator i= f; i.hasTerms(); i++)
718    result += evaluate (i.coeff(), g, h, powH, v)*power (x, i.exp());
719  return result;
720}
721
722
Note: See TracBrowser for help on using the repository browser.