source: git/factory/libfac/charset/csutil.cc @ 16055bd

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