source: git/libfac/charset/csutil.cc @ 924d8f

spielwiese
Last change on this file since 924d8f was 924d8f, checked in by Hans Schoenemann <hannes@…>, 14 years ago
svn does not support $ git-svn-id: file:///usr/local/Singular/svn/trunk@12915 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 22.6 KB
RevLine 
[1a80b4]1////////////////////////////////////////////////////////////
2// emacs edit mode for this file is -*- C++ -*-
[341696]3/* $Id$ */
[1a80b4]4/////////////////////////////////////////////////////////////
5// FACTORY - Includes
6#include <factory.h>
7// Factor - Includes
8#include <tmpl_inst.h>
9#include <Factor.h>
10#include <SqrFree.h>
[4a81ec]11#include <helpstuff.h>
12#include <homogfactor.h>
[1a80b4]13// Charset - Includes
14#include "csutil.h"
[9b25e7e]15extern void out_cf(char *s1,const CanonicalForm &f,char *s2);
[639047e]16extern CanonicalForm alg_lc(const CanonicalForm &f);
17extern int hasAlgVar(const CanonicalForm &f);
[1a80b4]18
19static bool
[4a81ec]20lowerRank ( const CanonicalForm & f, const CanonicalForm & g, int & ind )
[1a80b4]21{
22  int df, dg;
23  Variable vf = f.mvar(), vg = g.mvar();
[0be2bc]24
[56c5ce7]25  if ( f.inCoeffDomain() )
26  {
[0be2bc]27    if ( g.inCoeffDomain() ) ind= 1;
[4a81ec]28    return true;//( vg > vf );
29  }
30  else if ( g.inCoeffDomain() ) return false;
31  else if ( vf < vg ) return true;
[56c5ce7]32  else if ( vf == vg )
33  {
[4a81ec]34    df = degree( f ); dg = degree( g );
35    if ( df < dg ) return true;
36    else if ( df == dg ) return lowerRank( LC( f ), LC( g ) , ind);
37    else return false;
38  }
[1a80b4]39  return false;
40}
41
42CanonicalForm
43lowestRank( const CFList & F )
44{
45  CFListIterator i = F;
46  CanonicalForm f;
[4a81ec]47  int ind=0;
[0be2bc]48  if ( ! i.hasItem() )        return f;
[1a80b4]49  f = i.getItem(); ++i;
[56c5ce7]50  while ( i.hasItem() )
51  {
[e2ca88]52    //CERR << "comparing " << f << "  and " << i.getItem()
53    // << " == " << lowerRank( i.getItem(), f, ind ) << "\n";
[56c5ce7]54    if ( lowerRank( i.getItem(), f, ind ) )
55    {
56      if ( ind )
57      {
[0be2bc]58        CFList Itemlist= get_Terms(i.getItem());
59        CFList Flist= get_Terms(f);
60
61        // Have to further compare number of terms!
[e2ca88]62        //CERR << "compare terms! f= " << Flist.length() << "  item= "
63        //     << Itemlist.length() <<"\n";
[0be2bc]64        if ( Itemlist.length() < Flist.length()) f = i.getItem();
65        ind=0;
[4a81ec]66      }
[56c5ce7]67      else
68      {
[0be2bc]69        f = i.getItem();
[4a81ec]70      }
71    }
[1a80b4]72    ++i;
73  }
74  return f;
75}
76
77// old version
78// CanonicalForm
79// prem ( const CanonicalForm &f, const CanonicalForm &g )
80// {
81//   CanonicalForm ff, gg, cg;
82//   int df, dg;
83//   bool reord;
84//   Variable vf, vg, v;
85//
86//   if ( (vf = f.mvar()) < (vg = g.mvar()) ) return f;
87//   else {
88//     if ( vf == vg ) {
89//       ff = f; gg = g;
90//       reord = false;
91//       v = vg;
92//     }
[0be2bc]93//     else {
[1a80b4]94//       v = Variable(level(f.mvar()) + 1);
95//       ff = swapvar(f,vg,v);
96//       gg = swapvar(g,vg,v);
97//       reord=true;
98//     }
99//     cg = ini( gg, v );
100//     dg = degree( gg, v );
101//     while ( ( df = degree( ff, v ) ) >= dg )
102//       ff = cg * ff - power( v, df-dg ) * gg * LC( ff, v );
103//     if ( reord ) {
104//       return swapvar( ff, vg, v );
105//     }
106//     else
107//       return ff;
108//   }
109// }
110
111CanonicalForm
[56c5ce7]112Prem ( const CanonicalForm &f, const CanonicalForm &g )
113{
[4a81ec]114  CanonicalForm ff, gg, l, test, lu, lv, t, retvalue;
[1a80b4]115  int df, dg;
116  bool reord;
117  Variable vf, vg, v;
[4a81ec]118
[1a80b4]119  if ( (vf = f.mvar()) < (vg = g.mvar()) ) return f;
[56c5ce7]120  else
121  {
122    if ( vf == vg )
123    {
[1a80b4]124      ff = f; gg = g;
125      reord = false;
126      v = vg;
127    }
[56c5ce7]128    else
129    {
[1a80b4]130      v = Variable(level(f.mvar()) + 1);
131      ff = swapvar(f,vg,v);
132      gg = swapvar(g,vg,v);
133      reord=true;
134    }
135    dg = degree( gg, v );
136    df = degree( ff, v );
137    if (dg <= df) {l=LC(gg); gg = gg -LC(gg)*power(v,dg);}
138    else { l = 1; }
[6f801b]139    while ( ( dg <= df  ) && ( !ff.isZero()) )
[56c5ce7]140    {
[e2ca88]141      // CERR << "Start gcd..." << "\n";
[1a80b4]142      test = gcd(l,LC(ff));
[e2ca88]143      //CERR << "gcd(" << l << "," << LC(ff) << ")= " << test << "\n";
[1a80b4]144      lu = l/test; lv = LC(ff)/test;
145      t = power(v,df-dg) * gg * lv;
146      if ( df == 0 ){ ff = ff.genZero(); }
147      else { ff = ff - LC(ff)*power(v,df); }
148      ff = lu*ff - t;
149      df = degree( ff, v );
150    }
[56c5ce7]151    if ( reord )
152    {
153      retvalue= swapvar( ff, vg, v );
[4a81ec]154    }
[56c5ce7]155    else
156    {
157      retvalue= ff;
[4a81ec]158    }
159    return retvalue;
160  }
161}
162
163static CanonicalForm
[56c5ce7]164Sprem ( const CanonicalForm &f, const CanonicalForm &g, CanonicalForm & m, CanonicalForm & q )
165{
[4a81ec]166  CanonicalForm ff, gg, l, test, retvalue;
167  int df, dg,n;
168  bool reord;
169  Variable vf, vg, v;
170
[56c5ce7]171  if ( (vf = f.mvar()) < (vg = g.mvar()) )
172  {
[0be2bc]173    m=CanonicalForm(0); q=CanonicalForm(0);
[4a81ec]174    return f;
175  }
[56c5ce7]176  else
177  {
178    if ( vf == vg )
179    {
[4a81ec]180      ff = f; gg = g;
181      reord = false;
182      v = vg; // == x
183    }
[56c5ce7]184    else
185    {
[4a81ec]186      v = Variable(level(f.mvar()) + 1);
187      ff = swapvar(f,vg,v); // == r
188      gg = swapvar(g,vg,v); // == v
189      reord=true;
190    }
[0be2bc]191    dg = degree( gg, v ); // == dv
[4a81ec]192    df = degree( ff, v ); // == dr
193    if (dg <= df) {l=LC(gg); gg = gg -LC(gg)*power(v,dg);}
194    else { l = 1; }
195    n= 0;
[6f801b]196    while ( ( dg <= df  ) && ( !ff.isZero()) )
[56c5ce7]197    {
[4a81ec]198      test= power(v,df-dg) * gg * LC(ff);
199      if ( df == 0 ){ff= ff.genZero();}
200      else {ff= ff - LC(ff)*power(v,df);}
201      ff = l*ff-test;
202      df= degree(ff,v);
203      n++;
204    }
[56c5ce7]205    if ( reord )
206    {
207      retvalue= swapvar( ff, vg, v );
[4a81ec]208    }
[56c5ce7]209    else
210    {
211      retvalue= ff;
[4a81ec]212    }
213    m= power(l,n);
[dbcf42a]214    if ( fdivides(g,m*f-retvalue) )
[4a81ec]215      q= (m*f-retvalue)/g;
[56c5ce7]216    else
[4a81ec]217      q= CanonicalForm(0);
218    return retvalue;
[1a80b4]219  }
220}
221
[4a81ec]222CanonicalForm
[38e7b3]223divide( const CanonicalForm & ff, const CanonicalForm & f, const CFList & as)
224{
[4a81ec]225  CanonicalForm r,m,q;
226
[fb4f62e]227  //out_cf("divide f=",ff,"\n");
228  //out_cf("divide g=",f,"\n");
229  if (f.inCoeffDomain())
230  {
[56c5ce7]231    bool b=!isOn(SW_RATIONAL);
232    On(SW_RATIONAL);
[fb4f62e]233    q=ff/f;
234    if (b) Off(SW_RATIONAL);
[56c5ce7]235  }
[fb4f62e]236  else
237    r= Sprem(ff,f,m,q); //result in q, ignore r,m
[e2ca88]238  //CERR << "r= " << r << "  , m= " << m << "  , q= " << q << "\n";
[4a81ec]239  r= Prem(q,as);
[e2ca88]240  //CERR << "r= " << r << "\n";
[fb4f62e]241  //out_cf(" ->",r,"\n");
[4a81ec]242  return r;
243}
244
[56c5ce7]245// This function allows as to be empty; in that case, it is equivalent
246// to the previous version (treating no variables as algebraic).
[4a81ec]247static CanonicalForm
[56c5ce7]248myfitting( const CanonicalForm &f, const CFList &as )
[38e7b3]249{
[4a81ec]250 CanonicalForm rem=f;
251
[38e7b3]252 if ( !(rem.isZero()) )
253 {
[4a81ec]254   if ( getCharacteristic() > 0 )
255     return num((rem/lc(rem)));
[38e7b3]256   else
257   {
[4a81ec]258     On(SW_RATIONAL);
259     CanonicalForm temp= mapinto(rem);
[e2ca88]260//      CERR << "temp= " << temp << "\n";
261//      CERR << "lc(temp)= " << lc(temp) << "\n";
262//      CERR << "temp/lc(temp)= " << temp/lc(temp) << "\n";
263//      CERR << "num(rem/lc(rem))= " << num(rem/lc(rem)) << "\n";
[56c5ce7]264
265     // If as is of length 1, and its polynomial is level 1, then
266     // we treat the first variable as algebraic and invert the leading
267     // coefficient where this variable is part of the coefficient domain.
268
269     if (as.length() == 1 && level(as.getFirst()) == 1)
270     {
271       CanonicalForm lcoeff = temp;
272       while (level(lcoeff) > 1)
273       {
274         lcoeff = LC(lcoeff);
275       }
276       // out_cf("myfitting: lcoeff = ", lcoeff, "\n");
277
278       CanonicalForm p = as.getFirst();
279       // out_cf("myfitting: p = ", p, "\n");
280
281       CanonicalForm unused, inverse;
282
283       extgcd(lcoeff, p, inverse, unused);
284       // out_cf("myfitting: inverse = ", inverse, "\n");
285
286       // This may leave temp with non-integral coefficients,
287       // which will be cleaned up below.
288       temp = temp * inverse;
289       // out_cf("myfitting: temp = ", temp, "\n");
290     }
291
[17ce9d]292     temp= bCommonDen(temp/lc(temp))*(temp/lc(temp));
[4a81ec]293     Off(SW_RATIONAL);
294     rem= mapinto(temp);
295     return rem;
296   }
297 }
298 else
299   return rem;
300}
301
[1a80b4]302CanonicalForm
[56c5ce7]303Prem( const CanonicalForm &f, const CFList &L )
304{
[1a80b4]305  CanonicalForm rem = f;
306  CFListIterator i = L;
[56c5ce7]307  for ( i.lastItem(); i.hasItem(); i-- )
308  {
309    //CERR << "   PREM: Prem(" << rem << "," ;
[1a80b4]310    rem = Prem( rem, i.getItem() );
[56c5ce7]311    //CERR << "   PREM: Prem(" << rem << "," << i.getItem() << ")  = " << rem << "\n";
[1a80b4]312  }
[56c5ce7]313  return myfitting(rem, CFList());
314}
315
316CanonicalForm
317Prem( const CanonicalForm &f, const CFList &L, const CFList &as )
318{
319  CanonicalForm rem = f;
320  CFListIterator i = L;
321  for ( i.lastItem(); i.hasItem(); i-- )
322  {
323    //CERR << "   PREM: Prem(" << rem << "," ;
324    rem = Prem( rem, i.getItem() );
325    //CERR << "   PREM: Prem(" << rem << "," << i.getItem() << ")  = " << rem << "\n";
326  }
327  return myfitting(rem, as);
[1a80b4]328}
329
330CFList
[56c5ce7]331Prem( const CFList &AS, const CFList &L )
332{
[1a80b4]333  CFList Output;
[4a81ec]334
[1a80b4]335  for ( CFListIterator i=AS; i.hasItem(); i++ )
336    Output = Union(CFList(Prem(i.getItem(),L)), Output);
[4a81ec]337
[1a80b4]338  return Output;
339}
340
[4a81ec]341static CanonicalForm
[38e7b3]342premasb( const CanonicalForm & f, const CFList & as)
343{
[4a81ec]344  CanonicalForm remd=f;
345  CFList AS=as;
346
[38e7b3]347  if ( as.length() > 1 )
348  {
[4a81ec]349    AS.removeFirst(); // get rid of first elem
350    CanonicalForm elem;
[38e7b3]351    while ( ! AS.isEmpty() )
352    { // thats true for at least the first iteration
[4a81ec]353      elem= AS.getLast();
354      remd= Prem(remd,elem);
355      AS.removeLast();
356    }
357  }
358  CanonicalForm a,b;
359  if ( mydivremt(remd, as.getFirst(), a,b )){ remd= remd.genZero();}
360  else { remd= Prem(remd, as.getFirst()); }
361
362  return remd;
363}
364
365CFList
[38e7b3]366remsetb( const CFList & ps, const CFList & as)
367{
[4a81ec]368  CFList output;
369  CanonicalForm elem;
[38e7b3]370  for (CFListIterator i=ps; i.hasItem(); i++)
371  {
[4a81ec]372    elem= premasb(i.getItem(),as);
[6f801b]373    if ( !elem.isZero() ) output.append(elem);
[4a81ec]374  }
375  return output;
376}
377
[1a80b4]378// for characteristic sets
379//////////////////////////////////
380// replace the power of factors of polys in as by 1 if any
381static CFList
[38e7b3]382nopower( const CanonicalForm & init )
383{
[1a80b4]384  CFFList sqrfreelist;// = Factorize(init);//SqrFree(init);
385  CFList output;
386  CanonicalForm elem;
387  int count=0;
[0be2bc]388
[1a80b4]389  for ( CFIterator j=init; j.hasTerms(); j++ )
[56c5ce7]390    if (!(j.coeff().isOne()) ) count++;
[1a80b4]391  //  if ( init != 1 ){
[e2ca88]392  //  CERR << "nopower: f is " << init << "\n";
393  //  CERR << "nopower: count is " << count << "\n";}
[1a80b4]394  if ( count > 1 ) sqrfreelist = CFFList( CFFactor(init,1));
[38e7b3]395  else
396  {
[1a80b4]397    sqrfreelist = Factorize(init);
398    //sqrfreelist.removeFirst();
399  }
[38e7b3]400  for ( CFFListIterator i=sqrfreelist; i.hasItem(); i++ )
401  {
[1a80b4]402    elem=i.getItem().factor();
403    if ( cls(elem) > 0 ) output.append(elem);
404  }
405  return output;
406}
407
408// remove the content of polys in PS; add the removed content to
409// Remembern.FS2 ( the set of removed factors )
[0be2bc]410CFList
[d704fa]411removecontent ( const CFList & PS, PremForm & Remembern )
412{
413  CFListIterator i=PS;
414  if ((!i.hasItem()) || ( cls(PS.getFirst()) == 0 )) return PS;
415
[1a80b4]416  CFList output;
417  CanonicalForm cc,elem;
418
[d704fa]419  for (; i.hasItem(); i++)
420  {
[1a80b4]421    elem = i.getItem();
422    cc = content(elem, elem.mvar());
[d704fa]423    if ( cls(cc) > 0 )
424    {
[0be2bc]425      output.append(elem/cc);
426      Remembern.FS2 = Union(CFList(cc), Remembern.FS2);
[1a80b4]427    }
428    else{ output.append(elem); }
429  }
430  return output;
431}
432
433// remove possible factors in Remember.FS1 from poly r
434// Remember.FS2 contains all factors removed before
[0be2bc]435void
[56c5ce7]436removefactor( CanonicalForm & r , PremForm & Remembern)
437{
[1a80b4]438  int test;
439  CanonicalForm a,b,testelem;
440  CFList testlist;
441  int n=level(r);
442  CFListIterator j ;
443
[38e7b3]444  for ( int J=1; J<= n ; J++ )
445  {
[1a80b4]446    testlist.append(CanonicalForm(Variable(J)));
447  }
[0be2bc]448
[4a81ec]449  //  testlist = Union(Remembern.FS1, testlist); // add candidates
[1a80b4]450
451  // remove already removed factors
[38e7b3]452  for ( j = Remembern.FS2 ; j.hasItem(); j++ )
453  {
[1a80b4]454    testelem = j.getItem();
[38e7b3]455    while ( 1 )
456    {
[4a81ec]457      test = mydivremt(r,testelem,a,b);
[6f801b]458      if ( test && b.isZero() ) r = a;
[1a80b4]459      else break;
460    }
461  }
[4a81ec]462
[1a80b4]463  // Let's look if we have other canditates to remove
[38e7b3]464  for ( j = testlist ; j.hasItem(); j++ )
465  {
[1a80b4]466    testelem = j.getItem();
[4a81ec]467//    if ( testelem != r && testelem != r.mvar() ){
[38e7b3]468    if ( testelem != r )
469    {
470      while ( 1 )
471      {
[0be2bc]472        test = divremt(r,testelem,a,b);
[6f801b]473        if ( test && b.isZero() )
[38e7b3]474        {
[0be2bc]475          Remembern.FS2= Union(Remembern.FS2, CFList(testelem));
476          r = a;
477          if ( r == 1 ) break;
478        }
479        else break;
[1a80b4]480      }
481    }
482  }
[e2ca88]483  //  CERR << "Remembern.FS1 = " << Remembern.FS1 << "\n";
484  //  CERR << "Remembern.FS2 = " << Remembern.FS2 << "\n";
[4a81ec]485  //  Remembern.FS1 = Difference(Remembern.FS1, Remembern.FS2);
[e2ca88]486  //  CERR << "  New Remembern.FS1 = " << Remembern.FS1 << "\n";
[1a80b4]487}
488
489
490// all irreducible nonconstant factors of a set of polynomials
491CFList
[38e7b3]492factorps( const CFList &ps )
493{
[1a80b4]494  CFList qs;
495  CFFList q;
496  CanonicalForm elem;
[0be2bc]497
[38e7b3]498  for ( CFListIterator i=ps; i. hasItem(); i++ )
499  {
[1a80b4]500    q=Factorize(i.getItem());
501    q.removeFirst();
[0be2bc]502    // Next can be simplified ( first (already removed) elem in q is the only constant
[38e7b3]503    for ( CFFListIterator j=q; j.hasItem(); j++ )
504    {
[1a80b4]505      elem = j.getItem().factor();
506      if ( getNumVars(elem) > 0 )
[56c5ce7]507        qs= Union(qs, CFList(myfitting(elem, CFList())));
[1a80b4]508    }
509  }
510  return qs;
511}
512
513// the initial of poly f wrt to the order of the variables
514static CanonicalForm
[38e7b3]515inital( const CanonicalForm &f )
516{
[1a80b4]517  CanonicalForm leadcoeff;
518
519  if ( cls(f) == 0 ) {return f.genOne(); }
[38e7b3]520  else
521  {
[1a80b4]522    leadcoeff = LC(f,lvar(f));
[0be2bc]523    //    if ( leadcoeff != 0 )
[56c5ce7]524    return myfitting(leadcoeff, CFList()); //num(leadcoeff/lc(leadcoeff));
[4a81ec]525    //    else return leadcoeff;
[1a80b4]526  }
527}
528
529// the set of all nonconstant factors of initals of polys in as
530// CFList
531// initalset(const CFList &as){
532//   CanonicalForm elem;
533//   CFList is, iss,iniset;
534
535//   for ( CFListIterator i=as ; i.hasItem(); i++ ){
536//     elem = inital(i.getItem());
537//     if ( cls(elem) > 0 ) is.append(elem);
538//   }
539//   iss = factorps(is);
540//   for ( CFListIterator j=iss; j.hasItem();j++ ){
541//     elem = j.getItem();
542//     if ( cls(elem) > 0 ) iniset.append(num(elem/lc(elem)));
543//   }
544//   return iniset;
545// }
546
[fdf28d]547// the set of nonconstant initials of Cset
[1a80b4]548// with certain repeated factors cancelled
549CFList
[56c5ce7]550initalset1(const CFList & Cset)
551{
[1a80b4]552  CFList temp;
553  CFList initals;
554  CanonicalForm init;
555
[56c5ce7]556  for ( CFListIterator i = Cset ; i.hasItem(); i++ )
557  {
[1a80b4]558    initals= nopower( inital(i.getItem()) );
559    //    init= inital(i.getItem());
[56c5ce7]560    for ( CFListIterator j = initals; j.hasItem(); j++)
561    {
[1a80b4]562      init = j.getItem();
563      if ( cls(init) > 0 )
[0be2bc]564        temp= Union(temp, CFList(init));
[1a80b4]565    }
566  }
567  return temp;
568}
569
[fdf28d]570// the set of nonconstant initials of Cset of those polys
[1a80b4]571// not having their cls higher than reducible
572// with certain repeated factors cancelled
573CFList
[56c5ce7]574initalset2(const CFList & Cset, const CanonicalForm & reducible)
575{
[1a80b4]576  CFList temp;
577  CFList initals;
578  CanonicalForm init;
579  int clsred = cls(reducible);
580
[56c5ce7]581  for ( CFListIterator i = Cset ; i.hasItem(); i++ )
582  {
[1a80b4]583    init = i.getItem();
[56c5ce7]584    if ( cls(init) < clsred )
585    {
[1a80b4]586      initals= nopower( inital(init) );
587      //    init= inital(i.getItem());
[56c5ce7]588      for ( CFListIterator j = initals; j.hasItem(); j++)
589      {
[0be2bc]590        init = j.getItem();
591        if ( cls(init) > 0 )
592          temp= Union(temp, CFList(init));
[1a80b4]593      }
594    }
595  }
596  return temp;
597}
598
599//replace the power of factors of poly in CF init by 1 if any
600// and return the sqrfree poly
601// static CanonicalForm
602// nopower1( const CanonicalForm & init ){
603//   CFFList returnlist=Factorize(init);
604//   CanonicalForm elem, test=init.genOne();
605//   for ( CFFListIterator i= returnlist; i.hasItem(); i++){
606//     elem = i.getItem().factor();
607//     if ( cls(elem)>0 ) test *= elem;
608//   }
609//   return test;
610// }
611
612// the sequence of distinct factors of f
[0be2bc]613//CF pfactor( ..... )
[1a80b4]614
[4a81ec]615// //////////////////////////////////////////
616// // for IrrCharSeries
[1a80b4]617
618#ifdef IRRCHARSERIESDEBUG
619#  define DEBUGOUTPUT
620#else
621#  undef DEBUGOUTPUT
622#endif
623#include "debug.h"
624// examine the irreducibility of as for IrrCharSeries
625int
[56c5ce7]626irreducible( const CFList & AS)
627{
[1a80b4]628// AS is given by AS = { A1, A2, .. Ar }, d_i = degree(Ai)
629
630// 1) we test: if d_i > 1, d_j =1 for all j<>i, then AS is irreducible.
[56c5ce7]631  bool deg1=true;
632  for ( CFListIterator i = AS ; i.hasItem(); i++ )
633  {
634    if ( degree(i.getItem()) > 1 )
635    {
[1a80b4]636      if ( deg1 ) deg1=0;
637      else return 0; // found 2nd poly with deg > 1
638    }
639  }
640  return 1;
641}
642
643
644// select an item from PS for irras
645CFList
[56c5ce7]646select( const ListCFList & PS)
647{
[1a80b4]648  return PS.getFirst();
649}
650
651// divide list ppi in elems having length <= and > length
652void
[56c5ce7]653select( const ListCFList & ppi, int length, ListCFList & ppi1, ListCFList & ppi2)
654{
[1a80b4]655  CFList elem;
[56c5ce7]656  for ( ListCFListIterator i=ppi ; i.hasItem(); i++ )
657  {
[1a80b4]658    elem = i.getItem();
659    if ( ! elem.isEmpty() )
660      if ( length <= elem.length() ){ ppi2.append(elem); }
661      else { ppi1.append(elem); }
662  }
663}
664
665
666//////////////////////////////////////////////////////////////
667// help-functions for sets
668
669// is f in F ?
670static bool
[56c5ce7]671member( const CanonicalForm &f, const CFList &F)
672{
[1a80b4]673  for ( CFListIterator i=F; i.hasItem(); i++ )
674    if ( i.getItem() == f ) return 1;
675  return 0;
676}
677
678// are list A and B the same?
679bool
[56c5ce7]680same( const CFList &A, const CFList &B )
681{
[1a80b4]682  CFListIterator i;
683
684  for (i = A; i.hasItem(); i++)
685    if (! member(i.getItem(), B) )  return 0;
686  for (i = B; i.hasItem(); i++)
687    if (! member(i.getItem(), A) )  return 0;
688  return 1;
689}
690
691
692// is List cs contained in List of lists pi?
693bool
[56c5ce7]694member( const CFList & cs, const ListCFList & pi )
695{
[1a80b4]696  ListCFListIterator i;
697  CFList elem;
698
[56c5ce7]699  for ( i=pi; i.hasItem(); i++)
700  {
[1a80b4]701    elem = i.getItem();
702    if ( same(cs,elem) ) return 1;
703  }
704  return 0;
705}
706
[fdf28d]707// is PS a subset of Cset ?
[1a80b4]708bool
[56c5ce7]709subset( const CFList &PS, const CFList &Cset )
710{
[fdf28d]711  //  CERR << "subset: called with: " << PS << "   " << Cset << "\n";
[1a80b4]712  for ( CFListIterator i=PS; i.hasItem(); i++ )
[56c5ce7]713    if ( ! member(i.getItem(), Cset) )
714    {
[fdf28d]715      //      CERR << "subset: " << i.getItem() << "  is not a member of " << Cset << "\n";
[1a80b4]716      return 0;
717    }
718  return 1;
719}
720
721// Union of two List of Lists
722ListCFList
[56c5ce7]723MyUnion( const ListCFList & a, const ListCFList &b )
724{
[1a80b4]725  ListCFList output;
726  ListCFListIterator i;
727  CFList elem;
[0be2bc]728
[56c5ce7]729  for ( i = a ; i.hasItem(); i++ )
730  {
[1a80b4]731    elem=i.getItem();
[56c5ce7]732    if ( (! elem.isEmpty()) && ( ! member(elem,output)) )
733    {
[1a80b4]734      output.append(elem);
735    }
736  }
[0be2bc]737
[56c5ce7]738  for ( i = b ; i.hasItem(); i++ )
739  {
[1a80b4]740    elem=i.getItem();
[56c5ce7]741    if ( (! elem.isEmpty()) && ( ! member(elem,output)) )
742    {
[1a80b4]743      output.append(elem);
744    }
745  }
746  return output;
747}
748
749//if list b is member of the list of lists remove b and return the rest
750ListCFList
[56c5ce7]751MyDifference( const ListCFList & a, const CFList &b)
752{
[1a80b4]753  ListCFList output;
754  ListCFListIterator i;
755  CFList elem;
756
[56c5ce7]757  for ( i = a ; i.hasItem(); i++ )
758  {
[1a80b4]759    elem=i.getItem();
[56c5ce7]760    if ( (! elem.isEmpty()) && ( ! same(elem,b)) )
761    {
[1a80b4]762      output.append(elem);
763    }
764  }
765return output;
766}
767
768// remove all elements of b from list of lists a and return the rest
769ListCFList
[56c5ce7]770Minus( const ListCFList & a, const ListCFList & b)
771{
[1a80b4]772  ListCFList output=a;
773
774  for ( ListCFListIterator i=b; i.hasItem(); i++ )
775    output = MyDifference(output, i.getItem() );
776
777  return output;
778}
779
[639047e]780#if 0
781static CanonicalForm alg_lc(const CanonicalForm &f, const CFList as)
782{
783  for(CFListIterator i=as; i.hasItem(); i++)
784  {
785    if (f.mvar()==i.getItem().mvar()) return f;
[9b25e7e]786  }
[639047e]787  if (f.level()>0)
788  {
789    return alg_lc(f.LC(),as);
790  }
791  return CanonicalForm(1);
792}
793#endif
794
795CanonicalForm alg_gcd(const CanonicalForm & fff, const CanonicalForm &ggg,
796                      const CFList &as)
797{
798  CanonicalForm f=fff;
799  CanonicalForm g=ggg;
[9b25e7e]800  f=Prem(f,as);
801  g=Prem(g,as);
802  if ( f.isZero() )
803  {
804    if ( g.lc().sign() < 0 ) return -g;
805    else                     return g;
806  }
807  else  if ( g.isZero() )
808  {
809    if ( f.lc().sign() < 0 ) return -f;
810    else                     return f;
811  }
812  //out_cf("alg_gcd(",fff," , ");
813  //out_cf("",ggg,")\n");
[639047e]814  CanonicalForm res;
815  // does as appear in f and g ?
816  bool has_alg_var=false;
817  for ( CFListIterator j=as;j.hasItem(); j++ )
818  {
819    Variable v=j.getItem().mvar();
[9b25e7e]820    if (hasVar(f,v)) {has_alg_var=true; /*break;*/}
821    if (hasVar(g,v)) {has_alg_var=true; /*break;*/}
822    //out_cf("as:",j.getItem(),"\n");
823  }
[639047e]824  if (!has_alg_var)
825  {
826    if ((hasAlgVar(f))
827    || (hasAlgVar(g)))
828    {
829      Varlist ord;
830      for ( CFListIterator j=as;j.hasItem(); j++ )
831        ord.append(j.getItem().mvar());
832      res=algcd(f,g,as,ord);
833    }
834    else
835      res=gcd(f,g);
836    //out_cf("gcd=",res,"\n");
[9b25e7e]837    //out_cf("of f=",fff," , ");
838    //out_cf("and g=",ggg,"\n");
839
[639047e]840    return res;
[9b25e7e]841  }
[639047e]842
843  int mvf=f.level();
844  int mvg=g.level();
845  if (mvg > mvf)
[9b25e7e]846  {
[639047e]847    CanonicalForm tmp=f; f=g; g=tmp;
848    int tmp2=mvf; mvf=mvg; mvg=tmp2;
[9b25e7e]849  }
[639047e]850  if (g.inBaseDomain() || f.inBaseDomain())
851  {
852    //printf("const\n");
[9b25e7e]853    //out_cf("of f=",fff," , ");
854    //out_cf("and g=",ggg,"\n");
[639047e]855    return CanonicalForm(1);
[9b25e7e]856  }
857
[639047e]858  // gcd of all coefficients:
859  CFIterator i=f;
860  CanonicalForm c_gcd=i.coeff(); i++;
861  while( i.hasTerms())
862  {
863    c_gcd=alg_gcd(i.coeff(),c_gcd,as);
864    if (c_gcd.inBaseDomain()) break;
865    i++;
866  }
867  //printf("f.mvar=%d (%d), g.mvar=%d (%d)\n",f.level(),mvf,g.level(),mvg);
868  if (mvf!=mvg) // => mvf > mvg
869  {
870    res=alg_gcd(g,c_gcd,as);
871    //out_cf("alg_gcd1=",res,"\n");
[9b25e7e]872    //out_cf("of f=",fff," , ");
873    //out_cf("and g=",ggg,"\n");
[639047e]874    return res;
875  }
[9b25e7e]876  // now: mvf==mvg, f.level()==g.level()
[639047e]877  if (!c_gcd.inBaseDomain())
878  {
879    i=g;
880    while( i.hasTerms())
881    {
882      c_gcd=alg_gcd(i.coeff(),c_gcd,as);
883      if (c_gcd.inBaseDomain()) break;
884      i++;
885    }
[9b25e7e]886  }
[639047e]887
[fb4f62e]888  //f/=c_gcd;
889  //g/=c_gcd;
890  if (!c_gcd.isOne())
[56c5ce7]891  {
[fb4f62e]892    f=divide(f,c_gcd,as);
893    g=divide(g,c_gcd,as);
894  }
[639047e]895
896  CFList gg;
[9b25e7e]897  CanonicalForm r=1;
898  while (1)
[639047e]899  {
900    //printf("f.mvar=%d, g.mvar=%d\n",f.level(),g.level());
901    gg=as;
[9b25e7e]902    if (!g.inCoeffDomain()) gg.append(g);
[639047e]903    //out_cf("Prem(",f," , ");
904    //out_cf("",g,")\n");
[9b25e7e]905    if (g.inCoeffDomain()|| g.isZero())
906    {
907      //printf("in coeff domain:");
908      if (g.isZero()) { //printf("0\n");
909        i=f;
910        CanonicalForm f_gcd=i.coeff(); i++;
911        while( i.hasTerms())
912        {
913          f_gcd=alg_gcd(i.coeff(),f_gcd,as);
914          if (f_gcd.inBaseDomain()) break;
915          i++;
[56c5ce7]916        }
917        //out_cf("g=0 -> f:",f,"\n");
918        //out_cf("f_gcd:",f_gcd,"\n");
919        //out_cf("c_gcd:",c_gcd,"\n");
920        //f/=f_gcd;
921        f=divide(f,f_gcd,as);
922        //out_cf("f/f_gcd:",f,"\n");
923        f*=c_gcd;
924        //out_cf("f*c_gcd:",f,"\n");
[fb4f62e]925        CanonicalForm r_lc=alg_lc(f);
[56c5ce7]926        //out_cf("r_lc:",r_lc,"\n");
927        //f/=r_lc;
928        f=divide(f,r_lc,as);
929        //out_cf(" -> gcd:",f,"\n");
[fb4f62e]930        return f;
[9b25e7e]931      }
[56c5ce7]932      else { //printf("c\n");
933        return c_gcd;}
934    }
935    else if (g.level()==f.level()) r=Prem(f,gg,as);
[639047e]936    else
937    {
[9b25e7e]938      //printf("diff. level: %d, %d\n", f.level(), g.level());
[639047e]939      // g and f have different levels
[9b25e7e]940      //out_cf("c_gcd=",c_gcd,"\n");
941    //out_cf("of f=",fff," , ");
942    //out_cf("and g=",ggg,"\n");
943      return c_gcd;
[639047e]944    }
945    //out_cf("->",r,"\n");
946    f=g;
947    g=r;
948  }
[9b25e7e]949  if (degree(f,Variable(mvg))==0)
[639047e]950  {
951  // remark: mvf == mvg == f.level() ==g.level()
952    //out_cf("c_gcd=",c_gcd,"\n");
[9b25e7e]953    //out_cf("of f=",fff," , ");
954    //out_cf("and g=",ggg,"\n");
[639047e]955    return c_gcd;
[9b25e7e]956  }
[639047e]957  //out_cf("f=",f,"\n");
958  i=f;
959  CanonicalForm k=i.coeff(); i++;
960  //out_cf("k=",k,"\n");
961  while( i.hasTerms())
962  {
963    if (k.inCoeffDomain()) break;
964    k=alg_gcd(i.coeff(),k,as);
965    //out_cf("k=",k,"\n");
966    i++;
967  }
[9b25e7e]968  //out_cf("end-k=",k,"\n");
[639047e]969  f/=k;
[9b25e7e]970  //out_cf("f/k=",f,"\n");
[639047e]971  res=c_gcd*f;
972  CanonicalForm r_lc=alg_lc(res);
973  res/=r_lc;
974  //CanonicalForm r_lc=alg_lc(res,as);
975  //res/=r_lc;
976  //out_cf("alg_gcd2=",res,"\n");
[9b25e7e]977  //  out_cf("of f=",fff," , ");
978  //  out_cf("and g=",ggg,"\n");
979  //return res;
980  //if (res.level()<fff.level() && res.level() < ggg.level())
981  //  return alg_gcd(alg_gcd(fff,res,as),ggg,as);
982  //else
983    return res;
[639047e]984}
Note: See TracBrowser for help on using the repository browser.