source: git/factory/canonicalform.cc @ fd80670

spielwiese
Last change on this file since fd80670 was b79ed5, checked in by Hans Schoenemann <hannes@…>, 11 years ago
removed genCoeff: unused
  • Property mode set to 100644
File size: 55.6 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3#include "config.h"
4
5#include "cf_assert.h"
6#include "cf_factory.h"
7
8#include "cf_defs.h"
9#include "cf_globals.h"
10#include "canonicalform.h"
11#include "cf_iter.h"
12#include "int_cf.h"
13#include "cf_algorithm.h"
14#include "imm.h"
15#include "gfops.h"
16#include "cf_binom.h"
17
18#include <factory/cf_gmp.h>
19
20
21#ifndef NOSTREAMIO
22CanonicalForm readCF( ISTREAM& );
23#endif /* NOSTREAMIO */
24
25//{{{ initialization
26int initializeCharacteristic();
27
28#ifdef SINGULAR
29extern int mmInit(void);
30#endif
31
32int
33initCanonicalForm( void )
34{
35    static bool initialized = false;
36    if ( ! initialized ) {
37#if defined (SINGULAR)
38        (void)mmInit();
39#endif
40
41        (void)initializeCharacteristic();
42        initPT();
43        initialized = true;
44    }
45    return 1;
46}
47//}}}
48
49//{{{ constructors, destructors, selectors
50CanonicalForm::CanonicalForm( const char * str, const int base ) : value( CFFactory::basic( str, base ) )
51{
52}
53
54InternalCF*
55CanonicalForm::getval() const
56{
57    if ( is_imm( value ) )
58        return value;
59    else
60        return value->copyObject();
61}
62
63CanonicalForm
64CanonicalForm::deepCopy() const
65{
66    if ( is_imm( value ) )
67        return *this;
68    else
69        return CanonicalForm( value->deepCopyObject() );
70}
71//}}}
72
73//{{{ predicates
74#if 0
75bool
76CanonicalForm::isImm() const
77{
78    return is_imm( value );
79}
80#endif
81
82bool
83CanonicalForm::inZ() const
84{
85    if ( is_imm( value ) == INTMARK )
86        return true;
87    else if ( is_imm( value ) )
88        return false;
89    else
90        return value->levelcoeff() == IntegerDomain;
91}
92
93bool
94CanonicalForm::inQ() const
95{
96    if ( is_imm( value ) == INTMARK )
97        return true;
98    else if ( is_imm( value ) )
99        return false;
100    else
101        return value->levelcoeff() == IntegerDomain ||
102            value->levelcoeff() == RationalDomain;
103}
104
105bool
106CanonicalForm::inFF() const
107{
108    return is_imm( value ) == FFMARK;
109}
110
111bool
112CanonicalForm::inGF() const
113{
114    return is_imm( value ) == GFMARK;
115}
116
117bool
118CanonicalForm::inPP() const
119{
120    return ! is_imm( value ) && ( value->levelcoeff() == PrimePowerDomain );
121}
122
123bool
124CanonicalForm::inBaseDomain() const
125{
126    if ( is_imm( value ) )
127        return true;
128    else
129        return value->inBaseDomain();
130}
131
132bool
133CanonicalForm::inExtension() const
134{
135    if ( is_imm( value ) )
136        return false;
137    else
138        return value->inExtension();
139}
140
141bool
142CanonicalForm::inCoeffDomain() const
143{
144    if ( is_imm( value ) )
145        return true;
146    else
147        return value->inCoeffDomain();
148}
149
150bool
151CanonicalForm::inPolyDomain() const
152{
153    if ( is_imm( value ) )
154        return false;
155    else
156        return value->inPolyDomain();
157}
158
159bool
160CanonicalForm::inQuotDomain() const
161{
162    if ( is_imm( value ) )
163        return false;
164    else
165        return value->inQuotDomain();
166}
167
168bool
169CanonicalForm::isFFinGF() const
170{
171    return is_imm( value ) == GFMARK && gf_isff( imm2int( value ) );
172}
173
174bool
175CanonicalForm::isUnivariate() const
176{
177    if ( is_imm( value ) )
178        return false;
179    else
180        return value->isUnivariate();
181}
182
183// is_homogeneous returns 1 iff f is homogeneous, 0 otherwise//
184bool
185CanonicalForm::isHomogeneous() const
186{
187  if (this->isZero()) return true;
188  else if (this->inCoeffDomain()) return true;
189  else
190  {
191#if 0
192    CFIterator i;
193    int cdeg = -2, dummy;
194    for ( i = *this; i.hasTerms(); i++ )
195    {
196      if (!(i.coeff().isHomogeneous())) return false;
197      if ( (dummy = totaldegree( i.coeff() ) + i.exp()) != cdeg )
198      {
199         if (cdeg == -2) cdeg = dummy;
200         else return false;
201      }
202    }
203    return true;
204#else
205    CFList termlist= get_Terms(*this);
206    CFListIterator i;
207    int deg= totaldegree(termlist.getFirst());
208
209    for ( i=termlist; i.hasItem(); i++ )
210      if ( totaldegree(i.getItem()) != deg ) return false;
211    return true;
212#endif
213  }
214}
215
216//}}}
217
218//{{{ conversion functions
219long
220CanonicalForm::intval() const
221{
222    if ( is_imm( value ) )
223        return imm_intval( value );
224    else
225        return value->intval();
226}
227
228CanonicalForm
229CanonicalForm::mapinto () const
230{
231    //ASSERT( is_imm( value ) ||  ! value->inExtension(), "cannot map into different Extension" );
232    if ( is_imm( value ) )
233        if ( getCharacteristic() == 0 )
234            if ( is_imm( value ) == FFMARK )
235                return CanonicalForm( int2imm( ff_symmetric( imm2int( value ) ) ) );
236            else  if ( is_imm( value ) == GFMARK )
237                return CanonicalForm( int2imm( ff_symmetric( gf_gf2ff( imm2int( value ) ) ) ) );
238            else
239                return *this;
240        else  if ( CFFactory::gettype() == PrimePowerDomain )
241            return CanonicalForm( CFFactory::basic( imm2int( value ) ) );
242        else  if ( getGFDegree() == 1 )
243            return CanonicalForm( int2imm_p( ff_norm( imm2int( value ) ) ) );
244        else
245            return CanonicalForm( int2imm_gf( gf_int2gf( imm2int( value ) ) ) );
246    else  if ( value->inBaseDomain() )
247        if ( getCharacteristic() == 0 )
248            if ( value->levelcoeff() == PrimePowerDomain )
249                return CFFactory::basic( getmpi( value, true ) );
250            else
251                return *this;
252        else  if ( CFFactory::gettype() == PrimePowerDomain )
253        {
254            ASSERT( value->levelcoeff() == PrimePowerDomain || value->levelcoeff() == IntegerDomain, "no proper map defined" );
255            if ( value->levelcoeff() == PrimePowerDomain )
256                return *this;
257            else
258                return CFFactory::basic( getmpi( value ) );
259        }
260        else
261        {
262            int val;
263            if ( value->levelcoeff() == IntegerDomain )
264                val = value->intmod( ff_prime );
265            else  if ( value->levelcoeff() == RationalDomain )
266                return num().mapinto() / den().mapinto();
267            else {
268                ASSERT( 0, "illegal domain" );
269                return 0;
270            }
271            if ( getGFDegree() > 1 )
272                return CanonicalForm( int2imm_gf( gf_int2gf( val ) ) );
273            else
274                return CanonicalForm( int2imm_p( val ) );
275        }
276    else
277    {
278        Variable x = value->variable();
279        CanonicalForm result;
280        for ( CFIterator i = *this; i.hasTerms(); i++ )
281            result += (power( x, i.exp() ) * i.coeff().mapinto());
282        return result;
283    }
284}
285//}}}
286
287//{{{ CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const
288//{{{ docu
289//
290// lc(), Lc(), LC() - leading coefficient functions.
291//
292// All methods return CO if CO is in a base domain.
293//
294// lc() returns the leading coefficient of CO with respect to
295// lexicographic ordering.  Elements in an algebraic extension
296// are considered polynomials so lc() always returns a leading
297// coefficient in a base domain.  This method is useful to get
298// the base domain over which CO is defined.
299//
300// Lc() returns the leading coefficient of CO with respect to
301// lexicographic ordering.  In contrast to lc() elements in an
302// algebraic extension are considered coefficients so Lc() always
303// returns a leading coefficient in a coefficient domain.
304//
305// LC() returns the leading coefficient of CO where CO is
306// considered a univariate polynomial in its main variable.  An
307// element of an algebraic extension is considered an univariate
308// polynomial, too.
309//
310// LC( v ) returns the leading coefficient of CO where CO is
311// considered an univariate polynomial in the polynomial variable
312// v.
313// Note: If v is less than the main variable of CO we have to
314// swap variables which may be quite expensive.
315//
316// Examples:
317// Let x < y be polynomial variables, a an algebraic variable.
318//
319// (3*a*x*y^2+y+x).lc() = 3
320// (3*a*x*y^2+y+x).Lc() = 3*a
321// (3*a*x*y^2+y+x).LC() = 3*a*x
322// (3*a*x*y^2+y+x).LC( x ) = 3*a*y^2+1
323//
324// (3*a^2+4*a).lc() = 3
325// (3*a^2+4*a).Lc() = 3*a^2+4*a
326// (3*a^2+4*a).LC() = 3
327// (3*a^2+4*a).LC( x ) = 3*a^2+4*a
328//
329// See also: InternalCF::lc(), InternalCF::Lc(), InternalCF::LC(),
330// InternalPoly::lc(), InternalPoly::Lc(), InternalPoly::LC(),
331// ::lc(), ::Lc(), ::LC(), ::LC( v )
332//
333//}}}
334CanonicalForm
335CanonicalForm::lc () const
336{
337    if ( is_imm( value ) )
338        return *this;
339    else
340        return value->lc();
341}
342
343CanonicalForm
344CanonicalForm::Lc () const
345{
346    if ( is_imm( value ) || value->inCoeffDomain() )
347        return *this;
348    else
349        return value->Lc();
350}
351
352CanonicalForm
353CanonicalForm::LC () const
354{
355    if ( is_imm( value ) )
356        return *this;
357    else
358        return value->LC();
359}
360
361CanonicalForm
362CanonicalForm::LC ( const Variable & v ) const
363{
364    if ( is_imm( value ) || value->inCoeffDomain() )
365        return *this;
366
367    Variable x = value->variable();
368    if ( v > x )
369        return *this;
370    else if ( v == x )
371        return value->LC();
372    else {
373        CanonicalForm f = swapvar( *this, v, x );
374         if ( f.mvar() == x )
375             return swapvar( f.value->LC(), v, x );
376         else
377            // v did not occur in f
378            return *this;
379    }
380}
381//}}}
382
383//{{{ int CanonicalForm::degree (), degree ( v ) const
384//{{{ docu
385//
386// degree() - degree methods.
387//
388// Both methods returns -1 for the zero polynomial and 0 if
389// CO is in a base domain.
390//
391// degree() returns the degree of CO in its main variable.
392// Elements in an algebraic extension are considered polynomials.
393//
394// degree( v ) returns the degree of CO with respect to v.
395// Elements in an algebraic extension are considered polynomials,
396// and v may be algebraic.
397//
398// See also: InternalCf::degree(), InternalPoly::degree(),
399// ::degree(), ::degree( v )
400//
401//}}}
402int
403CanonicalForm::degree() const
404{
405    int what = is_imm( value );
406    if ( what )
407        if ( what == FFMARK )
408            return imm_iszero_p( value ) ? -1 : 0;
409        else if ( what == INTMARK )
410            return imm_iszero( value ) ? -1 : 0;
411        else
412            return imm_iszero_gf( value ) ? -1 : 0;
413    else
414        return value->degree();
415}
416
417int
418CanonicalForm::degree( const Variable & v ) const
419{
420    int what = is_imm( value );
421#if 0
422    if ( what )
423        if ( what == FFMARK )
424            return imm_iszero_p( value ) ? -1 : 0;
425        else if ( what == INTMARK )
426            return imm_iszero( value ) ? -1 : 0;
427        else
428            return imm_iszero_gf( value ) ? -1 : 0;
429    else if ( value->inBaseDomain() )
430        return value->degree();
431#else
432    switch(what)
433    {
434      case FFMARK: return imm_iszero_p( value ) ? -1 : 0;
435      case INTMARK: return imm_iszero( value ) ? -1 : 0;
436      case GFMARK:  return imm_iszero_gf( value ) ? -1 : 0;
437      case 0: if ( value->inBaseDomain() )
438              return value->degree();
439              break;
440    }
441#endif
442
443    Variable x = value->variable();
444    if ( v == x )
445        return value->degree();
446    else  if ( v > x )
447        // relatively to v, f is in a coefficient ring
448        return 0;
449    else {
450        int coeffdeg, result = 0;
451        // search for maximum of coefficient degree
452        for ( CFIterator i = *this; i.hasTerms(); i++ ) {
453            coeffdeg = i.coeff().degree( v );
454            if ( coeffdeg > result )
455                result = coeffdeg;
456        }
457        return result;
458    }
459}
460//}}}
461
462//{{{ CanonicalForm CanonicalForm::tailcoeff (), int CanonicalForm::taildegree () const
463//{{{ docu
464//
465// tailcoeff(), taildegree() - return least coefficient and
466//   degree, resp.
467//
468// tailcoeff() returns the coefficient of the term with the least
469// degree in CO where CO is considered an univariate polynomial
470// in its main variable.  Elements in an algebraic extension are
471// considered coefficients.
472//
473// taildegree() returns -1 for the zero polynomial, 0 if CO is in
474// a base domain, otherwise the least degree of CO where CO is
475// considered a univariate polynomial in its main variable.  In
476// contrast to tailcoeff(), elements in an algebraic extension
477// are considered polynomials, not coefficients, and such may
478// have a taildegree larger than zero.
479//
480// tailcoeff( v ) returns the tail coefficient of CO where CO is
481// considered an univariate polynomial in the polynomial variable
482// v.
483// Note: If v is less than the main variable of CO we have to
484// swap variables which may be quite expensive.
485//
486// See also: InternalCF::tailcoeff(), InternalCF::tailcoeff(),
487// InternalPoly::tailcoeff(), InternalPoly::taildegree,
488// ::tailcoeff(), ::taildegree()
489//
490//}}}
491CanonicalForm
492CanonicalForm::tailcoeff () const
493{
494    if ( is_imm( value ) || value->inCoeffDomain() )
495        return *this;
496    else
497        return value->tailcoeff();
498}
499
500CanonicalForm
501CanonicalForm::tailcoeff (const Variable& v) const
502{
503    if ( is_imm( value ) || value->inCoeffDomain() )
504        return *this;
505
506    Variable x = value->variable();
507    if ( v > x )
508        return *this;
509    else if ( v == x )
510        return value->tailcoeff();
511    else {
512        CanonicalForm f = swapvar( *this, v, x );
513         if ( f.mvar() == x )
514             return swapvar( f.value->tailcoeff(), v, x );
515         else
516            // v did not occur in f
517            return *this;
518    }
519}
520
521int
522CanonicalForm::taildegree () const
523{
524    int what = is_imm( value );
525    if ( what )
526        if ( what == FFMARK )
527            return imm_iszero_p( value ) ? -1 : 0;
528        else if ( what == INTMARK )
529            return imm_iszero( value ) ? -1 : 0;
530        else
531            return imm_iszero_gf( value ) ? -1 : 0;
532    else
533        return value->taildegree();
534}
535//}}}
536
537//{{{ int CanonicalForm::level (), Variable CanonicalForm::mvar () const
538//{{{ docu
539//
540// level(), mvar() - return level and main variable of CO.
541//
542// level() returns the level of CO.  For a list of the levels and
543// their meanings, see cf_defs.h.
544//
545// mvar() returns the main variable of CO or Variable() if CO is
546// in a base domain.
547//
548// See also: InternalCF::level(), InternalCF::variable(),
549// InternalPoly::level(), InternalPoly::variable(), ::level(),
550// ::mvar()
551//
552//}}}
553int
554CanonicalForm::level () const
555{
556    if ( is_imm( value ) )
557        return LEVELBASE;
558    else
559        return value->level();
560}
561
562Variable
563CanonicalForm::mvar () const
564{
565    if ( is_imm( value ) )
566        return Variable();
567    else
568        return value->variable();
569}
570//}}}
571
572//{{{ CanonicalForm CanonicalForm::num (), den () const
573//{{{ docu
574//
575// num(), den() - return numinator and denominator of CO.
576//
577// num() returns the numinator of CO if CO is a rational number,
578// CO itself otherwise.
579//
580// den() returns the denominator of CO if CO is a rational
581// number, 1 (from the current domain!) otherwise.
582//
583// See also: InternalCF::num(), InternalCF::den(),
584// InternalRational::num(), InternalRational::den(), ::num(),
585// ::den()
586//
587//}}}
588CanonicalForm
589CanonicalForm::num () const
590{
591    if ( is_imm( value ) )
592        return *this;
593    else
594        return CanonicalForm( value->num() );
595}
596
597CanonicalForm
598CanonicalForm::den () const
599{
600    if ( is_imm( value ) )
601        return CanonicalForm( 1 );
602    else
603        return CanonicalForm( value->den() );
604}
605//}}}
606
607//{{{ assignment operators
608CanonicalForm &
609CanonicalForm::operator += ( const CanonicalForm & cf )
610{
611    int what = is_imm( value );
612    if ( what ) {
613        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
614        if ( (what = is_imm( cf.value )) == FFMARK )
615            value = imm_add_p( value, cf.value );
616        else  if ( what == GFMARK )
617            value = imm_add_gf( value, cf.value );
618        else  if ( what )
619            value = imm_add( value, cf.value );
620        else {
621            InternalCF * dummy = cf.value->copyObject();
622            value = dummy->addcoeff( value );
623        }
624    }
625    else  if ( is_imm( cf.value ) )
626        value = value->addcoeff( cf.value );
627    else  if ( value->level() == cf.value->level() ) {
628        if ( value->levelcoeff() == cf.value->levelcoeff() )
629            value = value->addsame( cf.value );
630        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
631            value = value->addcoeff( cf.value );
632        else {
633            InternalCF * dummy = cf.value->copyObject();
634            dummy = dummy->addcoeff( value );
635            if ( value->deleteObject() ) delete value;
636            value = dummy;
637        }
638    }
639    else  if ( level() > cf.level() )
640        value = value->addcoeff( cf.value );
641    else {
642        InternalCF * dummy = cf.value->copyObject();
643        dummy = dummy->addcoeff( value );
644        if ( value->deleteObject() ) delete value;
645        value = dummy;
646    }
647    return *this;
648}
649
650CanonicalForm &
651CanonicalForm::operator -= ( const CanonicalForm & cf )
652{
653    int what = is_imm( value );
654    if ( what ) {
655        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
656        if ( (what = is_imm( cf.value )) == FFMARK )
657            value = imm_sub_p( value, cf.value );
658        else  if ( what == GFMARK )
659            value = imm_sub_gf( value, cf.value );
660        else  if ( what )
661            value = imm_sub( value, cf.value );
662        else {
663            InternalCF * dummy = cf.value->copyObject();
664            value = dummy->subcoeff( value, true );
665        }
666    }
667    else  if ( is_imm( cf.value ) )
668        value = value->subcoeff( cf.value, false );
669    else  if ( value->level() == cf.value->level() ) {
670        if ( value->levelcoeff() == cf.value->levelcoeff() )
671            value = value->subsame( cf.value );
672        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
673            value = value->subcoeff( cf.value, false );
674        else {
675            InternalCF * dummy = cf.value->copyObject();
676            dummy = dummy->subcoeff( value, true );
677            if ( value->deleteObject() ) delete value;
678            value = dummy;
679        }
680    }
681    else  if ( level() > cf.level() )
682        value = value->subcoeff( cf.value, false );
683    else {
684        InternalCF * dummy = cf.value->copyObject();
685        dummy = dummy->subcoeff( value, true );
686        if ( value->deleteObject() ) delete value;
687        value = dummy;
688    }
689    return *this;
690}
691
692CanonicalForm &
693CanonicalForm::operator *= ( const CanonicalForm & cf )
694{
695    int what = is_imm( value );
696    if ( what ) {
697        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
698        if ( (what = is_imm( cf.value )) == FFMARK )
699            value = imm_mul_p( value, cf.value );
700        else  if ( what == GFMARK )
701            value = imm_mul_gf( value, cf.value );
702        else  if ( what )
703            value = imm_mul( value, cf.value );
704        else {
705            InternalCF * dummy = cf.value->copyObject();
706            value = dummy->mulcoeff( value );
707        }
708    }
709    else  if ( is_imm( cf.value ) )
710        value = value->mulcoeff( cf.value );
711    else  if ( value->level() == cf.value->level() ) {
712        if ( value->levelcoeff() == cf.value->levelcoeff() )
713            value = value->mulsame( cf.value );
714        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
715            value = value->mulcoeff( cf.value );
716        else {
717            InternalCF * dummy = cf.value->copyObject();
718            dummy = dummy->mulcoeff( value );
719            if ( value->deleteObject() ) delete value;
720            value = dummy;
721        }
722    }
723    else  if ( level() > cf.level() )
724        value = value->mulcoeff( cf.value );
725    else {
726        InternalCF * dummy = cf.value->copyObject();
727        dummy = dummy->mulcoeff( value );
728        if ( value->deleteObject() ) delete value;
729        value = dummy;
730    }
731    return *this;
732}
733
734CanonicalForm &
735CanonicalForm::operator /= ( const CanonicalForm & cf )
736{
737    int what = is_imm( value );
738    if ( what ) {
739        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
740        if ( (what = is_imm( cf.value )) == FFMARK )
741            value = imm_div_p( value, cf.value );
742        else  if ( what == GFMARK )
743            value = imm_div_gf( value, cf.value );
744        else  if ( what )
745            value = imm_divrat( value, cf.value );
746        else {
747            InternalCF * dummy = cf.value->copyObject();
748            value = dummy->dividecoeff( value, true );
749        }
750    }
751    else  if ( is_imm( cf.value ) )
752        value = value->dividecoeff( cf.value, false );
753    else  if ( value->level() == cf.value->level() ) {
754        if ( value->levelcoeff() == cf.value->levelcoeff() )
755            value = value->dividesame( cf.value );
756        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
757            value = value->dividecoeff( cf.value, false );
758        else {
759            InternalCF * dummy = cf.value->copyObject();
760            dummy = dummy->dividecoeff( value, true );
761            if ( value->deleteObject() ) delete value;
762            value = dummy;
763        }
764    }
765    else  if ( level() > cf.level() )
766        value = value->dividecoeff( cf.value, false );
767    else {
768        InternalCF * dummy = cf.value->copyObject();
769        dummy = dummy->dividecoeff( value, true );
770        if ( value->deleteObject() ) delete value;
771        value = dummy;
772    }
773    return *this;
774}
775
776CanonicalForm &
777CanonicalForm::div ( const CanonicalForm & cf )
778{
779    int what = is_imm( value );
780    if ( what ) {
781        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
782        if ( (what = is_imm( cf.value )) == FFMARK )
783            value = imm_div_p( value, cf.value );
784        else  if ( what == GFMARK )
785            value = imm_div_gf( value, cf.value );
786        else  if ( what )
787            value = imm_div( value, cf.value );
788        else {
789            InternalCF * dummy = cf.value->copyObject();
790            value = dummy->divcoeff( value, true );
791        }
792    }
793    else  if ( is_imm( cf.value ) )
794        value = value->divcoeff( cf.value, false );
795    else  if ( value->level() == cf.value->level() ) {
796        if ( value->levelcoeff() == cf.value->levelcoeff() )
797            value = value->divsame( cf.value );
798        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
799            value = value->divcoeff( cf.value, false );
800        else {
801            InternalCF * dummy = cf.value->copyObject();
802            dummy = dummy->divcoeff( value, true );
803            if ( value->deleteObject() ) delete value;
804            value = dummy;
805        }
806    }
807    else  if ( level() > cf.level() )
808        value = value->divcoeff( cf.value, false );
809    else {
810        InternalCF * dummy = cf.value->copyObject();
811        dummy = dummy->divcoeff( value, true );
812        if ( value->deleteObject() ) delete value;
813        value = dummy;
814    }
815    return *this;
816}
817
818//same as divremt but handles zero divisors in case we are in Z_p[x]/(f) where f is not irreducible
819CanonicalForm &
820CanonicalForm::tryDiv ( const CanonicalForm & cf, const CanonicalForm& M, bool& fail )
821{
822    ASSERT (getCharacteristic() > 0, "expected positive characteristic");
823    ASSERT (!getReduce (M.mvar()), "do not reduce modulo M");
824    fail= false;
825    int what = is_imm( value );
826    if ( what ) {
827        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
828        if ( (what = is_imm( cf.value )) == FFMARK )
829            value = imm_div_p( value, cf.value );
830        else  if ( what == GFMARK )
831            value = imm_div_gf( value, cf.value );
832        else {
833            InternalCF * dummy = cf.value->copyObject();
834            value = dummy->divcoeff( value, true );
835        }
836    }
837    else  if ( is_imm( cf.value ) )
838        value = value->tryDivcoeff (cf.value, false, M, fail);
839    else  if ( value->level() == cf.value->level() ) {
840        if ( value->levelcoeff() == cf.value->levelcoeff() )
841            value = value->tryDivsame( cf.value, M, fail );
842        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
843            value = value->tryDivcoeff( cf.value, false, M, fail );
844        else {
845            InternalCF * dummy = cf.value->copyObject();
846            dummy = dummy->tryDivcoeff( value, true, M, fail );
847            if ( value->deleteObject() ) delete value;
848            value = dummy;
849        }
850    }
851    else  if ( level() > cf.level() )
852        value = value->tryDivcoeff( cf.value, false, M, fail );
853    else {
854        InternalCF * dummy = cf.value->copyObject();
855        dummy = dummy->tryDivcoeff( value, true, M, fail );
856        if ( value->deleteObject() ) delete value;
857        value = dummy;
858    }
859    return *this;
860}
861
862CanonicalForm &
863CanonicalForm::operator %= ( const CanonicalForm & cf )
864{
865    int what = is_imm( value );
866    if ( what ) {
867        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
868        if ( (what = is_imm( cf.value )) == FFMARK )
869            value = imm_mod_p( value, cf.value );
870        else  if ( what == GFMARK )
871            value = imm_mod_gf( value, cf.value );
872        else  if ( what )
873            value = imm_mod( value, cf.value );
874        else {
875            InternalCF * dummy = cf.value->copyObject();
876            value = dummy->modulocoeff( value, true );
877        }
878    }
879    else  if ( is_imm( cf.value ) )
880        value = value->modulocoeff( cf.value, false );
881    else  if ( value->level() == cf.value->level() ) {
882        if ( value->levelcoeff() == cf.value->levelcoeff() )
883            value = value->modulosame( cf.value );
884        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
885            value = value->modulocoeff( cf.value, false );
886        else {
887            InternalCF * dummy = cf.value->copyObject();
888            dummy = dummy->modulocoeff( value, true );
889            if ( value->deleteObject() ) delete value;
890            value = dummy;
891        }
892    }
893    else  if ( level() > cf.level() )
894        value = value->modulocoeff( cf.value, false );
895    else {
896        InternalCF * dummy = cf.value->copyObject();
897        dummy = dummy->modulocoeff( value, true );
898        if ( value->deleteObject() ) delete value;
899        value = dummy;
900    }
901    return *this;
902}
903
904CanonicalForm &
905CanonicalForm::mod ( const CanonicalForm & cf )
906{
907    int what = is_imm( value );
908    if ( what ) {
909        ASSERT ( ! is_imm( cf.value ) || (what==is_imm( cf.value )), "illegal base coefficients" );
910        if ( (what = is_imm( cf.value )) == FFMARK )
911            value = imm_mod_p( value, cf.value );
912        else  if ( what == GFMARK )
913            value = imm_mod_gf( value, cf.value );
914        else  if ( what )
915            value = imm_mod( value, cf.value );
916        else {
917            InternalCF * dummy = cf.value->copyObject();
918            value = dummy->modcoeff( value, true );
919        }
920    }
921    else  if ( is_imm( cf.value ) )
922        value = value->modcoeff( cf.value, false );
923    else  if ( value->level() == cf.value->level() ) {
924        if ( value->levelcoeff() == cf.value->levelcoeff() )
925            value = value->modsame( cf.value );
926        else  if ( value->levelcoeff() > cf.value->levelcoeff() )
927            value = value->modcoeff( cf.value, false );
928        else {
929            InternalCF * dummy = cf.value->copyObject();
930            dummy = dummy->modcoeff( value, true );
931            if ( value->deleteObject() ) delete value;
932            value = dummy;
933        }
934    }
935    else  if ( level() > cf.level() )
936        value = value->modcoeff( cf.value, false );
937    else {
938        InternalCF * dummy = cf.value->copyObject();
939        dummy = dummy->modcoeff( value, true );
940        if ( value->deleteObject() ) delete value;
941        value = dummy;
942    }
943    return *this;
944}
945
946void
947divrem ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r )
948{
949    InternalCF * qq = 0, * rr = 0;
950    int what = is_imm( f.value );
951    if ( what )
952        if ( is_imm( g.value ) ) {
953            if ( what == FFMARK )
954                imm_divrem_p( f.value, g.value, qq, rr );
955            else  if ( what == GFMARK )
956                imm_divrem_gf( f.value, g.value, qq, rr );
957            else
958                imm_divrem( f.value, g.value, qq, rr );
959        }
960        else
961            g.value->divremcoeff( f.value, qq, rr, true );
962    else  if ( (what=is_imm( g.value )) )
963        f.value->divremcoeff( g.value, qq, rr, false );
964    else  if ( f.value->level() == g.value->level() )
965        if ( f.value->levelcoeff() == g.value->levelcoeff() )
966            f.value->divremsame( g.value, qq, rr );
967        else  if ( f.value->levelcoeff() > g.value->levelcoeff() )
968            f.value->divremcoeff( g.value, qq, rr, false );
969        else
970            g.value->divremcoeff( f.value, qq, rr, true );
971    else  if ( f.value->level() > g.value->level() )
972        f.value->divremcoeff( g.value, qq, rr, false );
973    else
974        g.value->divremcoeff( f.value, qq, rr, true );
975    ASSERT( qq != 0 && rr != 0, "error in divrem" );
976    q = CanonicalForm( qq );
977    r = CanonicalForm( rr );
978}
979
980bool
981divremt ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r )
982{
983    InternalCF * qq = 0, * rr = 0;
984    int what = is_imm( f.value );
985    bool result = true;
986    if ( what )
987        if ( is_imm( g.value ) ) {
988            if ( what == FFMARK )
989                imm_divrem_p( f.value, g.value, qq, rr );
990            else  if ( what == GFMARK )
991                imm_divrem_gf( f.value, g.value, qq, rr );
992            else
993                imm_divrem( f.value, g.value, qq, rr );
994        }
995        else
996            result = g.value->divremcoefft( f.value, qq, rr, true );
997    else  if ( (what=is_imm( g.value )) )
998        result = f.value->divremcoefft( g.value, qq, rr, false );
999    else  if ( f.value->level() == g.value->level() )
1000        if ( f.value->levelcoeff() == g.value->levelcoeff() )
1001            result = f.value->divremsamet( g.value, qq, rr );
1002        else  if ( f.value->levelcoeff() > g.value->levelcoeff() )
1003            result = f.value->divremcoefft( g.value, qq, rr, false );
1004        else
1005            result = g.value->divremcoefft( f.value, qq, rr, true );
1006    else  if ( f.value->level() > g.value->level() )
1007        result = f.value->divremcoefft( g.value, qq, rr, false );
1008    else
1009        result = g.value->divremcoefft( f.value, qq, rr, true );
1010    if ( result ) {
1011        ASSERT( qq != 0 && rr != 0, "error in divrem" );
1012        q = CanonicalForm( qq );
1013        r = CanonicalForm( rr );
1014    }
1015    else {
1016        q = 0; r = 0;
1017    }
1018    return result;
1019}
1020
1021//same as divremt but handles zero divisors in case we are in Z_p[x]/(f) where f is not irreducible
1022bool
1023tryDivremt ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r, const CanonicalForm& M, bool& fail )
1024{
1025    ASSERT (getCharacteristic() > 0, "expected positive characteristic");
1026    ASSERT (!getReduce (M.mvar()), "do not reduce modulo M");
1027    fail= false;
1028    InternalCF * qq = 0, * rr = 0;
1029    int what = is_imm( f.value );
1030    bool result = true;
1031    if ( what )
1032        if ( is_imm( g.value ) ) {
1033            if ( what == FFMARK )
1034                imm_divrem_p( f.value, g.value, qq, rr );
1035            else  if ( what == GFMARK )
1036                imm_divrem_gf( f.value, g.value, qq, rr );
1037        }
1038        else
1039            result = g.value->tryDivremcoefft( f.value, qq, rr, true, M, fail );
1040    else  if ( (what=is_imm( g.value )) )
1041        result = f.value->tryDivremcoefft( g.value, qq, rr, false, M, fail );
1042    else  if ( f.value->level() == g.value->level() )
1043        if ( f.value->levelcoeff() == g.value->levelcoeff() )
1044            result = f.value->tryDivremsamet( g.value, qq, rr, M, fail );
1045        else  if ( f.value->levelcoeff() > g.value->levelcoeff() )
1046            result = f.value->tryDivremcoefft( g.value, qq, rr, false, M, fail );
1047        else
1048            result = g.value->tryDivremcoefft( f.value, qq, rr, true, M, fail );
1049    else  if ( f.value->level() > g.value->level() )
1050        result = f.value->tryDivremcoefft( g.value, qq, rr, false, M, fail );
1051    else
1052        result = g.value->tryDivremcoefft( f.value, qq, rr, true, M, fail );
1053    if (fail)
1054    {
1055      q= 0;
1056      r= 0;
1057      return false;
1058    }
1059    if ( result ) {
1060        ASSERT( qq != 0 && rr != 0, "error in divrem" );
1061        q = CanonicalForm( qq );
1062        r = CanonicalForm( rr );
1063        q= reduce (q, M);
1064        r= reduce (r, M);
1065    }
1066    else {
1067        q = 0; r = 0;
1068    }
1069    return result;
1070}
1071
1072//}}}
1073
1074//{{{ CanonicalForm CanonicalForm::operator () ( f ), operator () ( f, v ) const
1075//{{{ docu
1076//
1077// operator ()() - evaluation operator.
1078//
1079// Both operators return CO if CO is in a base domain.
1080//
1081// operator () ( f ) returns CO with f inserted for the main
1082// variable.  Elements in an algebraic extension are considered
1083// polynomials.
1084//
1085// operator () ( f, v ) returns CO with f inserted for v.
1086// Elements in an algebraic extension are considered polynomials
1087// and v may be an algebraic variable.
1088//
1089//}}}
1090CanonicalForm
1091CanonicalForm::operator () ( const CanonicalForm & f ) const
1092{
1093    if ( is_imm( value ) || value->inBaseDomain() )
1094        return *this;
1095    else {
1096#if 0
1097        CFIterator i = *this;
1098        int lastExp = i.exp();
1099        CanonicalForm result = i.coeff();
1100        i++;
1101        while ( i.hasTerms() ) {
1102            if ( (lastExp - i.exp()) == 1 )
1103                result *= f;
1104            else
1105                result *= power( f, lastExp - i.exp() );
1106            result += i.coeff();
1107            lastExp = i.exp();
1108            i++;
1109        }
1110        if ( lastExp != 0 )
1111            result *= power( f, lastExp );
1112#else
1113        CFIterator i = *this;
1114        int lastExp = i.exp();
1115        CanonicalForm result = i.coeff();
1116        i++;
1117        while ( i.hasTerms() )
1118        {
1119            int i_exp=i.exp();
1120            if ( (lastExp - i_exp /* i.exp()*/) == 1 )
1121                result *= f;
1122            else
1123                result *= power( f, lastExp - i_exp /*i.exp()*/ );
1124            result += i.coeff();
1125            lastExp = i_exp /*i.exp()*/;
1126            i++;
1127        }
1128        if ( lastExp != 0 )
1129            result *= power( f, lastExp );
1130#endif
1131        return result;
1132    }
1133}
1134
1135CanonicalForm
1136CanonicalForm::operator () ( const CanonicalForm & f, const Variable & v ) const
1137{
1138    if ( is_imm( value ) || value->inBaseDomain() )
1139        return *this;
1140
1141    Variable x = value->variable();
1142    if ( v > x )
1143        return *this;
1144    else  if ( v == x )
1145        return (*this)( f );
1146    else {
1147        // v is less than main variable of f
1148        CanonicalForm result = 0;
1149        for ( CFIterator i = *this; i.hasTerms(); i++ )
1150            result += i.coeff()( f, v ) * power( x, i.exp() );
1151        return result;
1152    }
1153}
1154//}}}
1155
1156//{{{ CanonicalForm CanonicalForm::operator [] ( int i ) const
1157//{{{ docu
1158//
1159// operator []() - return i'th coefficient from CO.
1160//
1161// Returns CO if CO is in a base domain and i equals zero.
1162// Returns zero (from the current domain) if CO is in a base
1163// domain and i is larger than zero.  Otherwise, returns the
1164// coefficient to x^i in CO (if x denotes the main variable of
1165// CO) or zero if CO does not contain x^i.  Elements in an
1166// algebraic extension are considered polynomials.  i should be
1167// larger or equal zero.
1168//
1169// Note: Never use a loop like
1170//
1171// for ( int i = degree( f ); i >= 0; i-- )
1172//     foo( i, f[ i ] );
1173//
1174// which is much slower than
1175//
1176// for ( int i = degree( f ), CFIterator I = f; I.hasTerms(); I++ ) {
1177//     // fill gap with zeroes
1178//     for ( ; i > I.exp(); i-- )
1179//         foo( i, 0 );
1180//     // at this point, i == I.exp()
1181//     foo( i, i.coeff() );
1182//     i--;
1183// }
1184// // work through trailing zeroes
1185// for ( ; i >= 0; i-- )
1186//     foo( i, 0 );
1187//
1188//}}}
1189CanonicalForm
1190CanonicalForm::operator [] ( int i ) const
1191{
1192    ASSERT( i >= 0, "index to operator [] less than zero" );
1193    if ( is_imm( value ) )
1194        if ( i == 0 )
1195            return *this;
1196        else
1197            return CanonicalForm( 0 );
1198    else
1199        return value->coeff( i );
1200}
1201//}}}
1202
1203//{{{ CanonicalForm CanonicalForm::deriv (), deriv ( x )
1204//{{{ docu
1205//
1206// deriv() - return the formal derivation of CO.
1207//
1208// deriv() derives CO with respect to its main variable.  Returns
1209// zero from the current domain if f is in a coefficient domain.
1210//
1211// deriv( x ) derives CO with respect to x.  x should be a
1212// polynomial variable.  Returns zero from the current domain if
1213// f is in a coefficient domain.
1214//
1215// See also: ::deriv()
1216//
1217//}}}
1218CanonicalForm
1219CanonicalForm::deriv () const
1220{
1221    if ( is_imm( value ) || value->inCoeffDomain() )
1222        return CanonicalForm( 0 );
1223    else {
1224        CanonicalForm result = 0;
1225        Variable x = value->variable();
1226        for ( CFIterator i = *this; i.hasTerms(); i++ )
1227            if ( i.exp() > 0 )
1228                result += power( x, i.exp()-1 ) * i.coeff() * i.exp();
1229        return result;
1230    }
1231}
1232
1233CanonicalForm
1234CanonicalForm::deriv ( const Variable & x ) const
1235{
1236    ASSERT( x.level() > 0, "cannot derive with respect to algebraic variables" );
1237    if ( is_imm( value ) || value->inCoeffDomain() )
1238        return CanonicalForm( 0 );
1239
1240    Variable y = value->variable();
1241    if ( x > y )
1242        return CanonicalForm( 0 );
1243    else if ( x == y )
1244        return deriv();
1245    else {
1246        CanonicalForm result = 0;
1247        for ( CFIterator i = *this; i.hasTerms(); i++ )
1248            result += i.coeff().deriv( x ) * power( y, i.exp() );
1249        return result;
1250    }
1251}
1252//}}}
1253
1254//{{{ int CanonicalForm::sign () const
1255//{{{ docu
1256//
1257// sign() - return sign of CO.
1258//
1259// If CO is an integer or a rational number, the sign is defined
1260// as usual.  If CO is an element of a prime power domain or of
1261// FF(p) and SW_SYMMETRIC_FF is on, the sign of CO is the sign of
1262// the symmetric representation of CO.  If CO is in GF(q) or in
1263// FF(p) and SW_SYMMETRIC_FF is off, the sign of CO is zero iff
1264// CO is zero, otherwise the sign is one.
1265//
1266// If CO is a polynomial or in an extension of one of the base
1267// domains, the sign of CO is the sign of its leading
1268// coefficient.
1269//
1270// See also: InternalCF::sign(), InternalInteger::sign(),
1271// InternalPrimePower::sign(), InternalRational::sign(),
1272// InternalPoly::sign(), imm_sign(), gf_sign()
1273//
1274//}}}
1275int
1276CanonicalForm::sign () const
1277{
1278    if ( is_imm( value ) )
1279        return imm_sign( value );
1280    else
1281        return value->sign();
1282}
1283//}}}
1284
1285//{{{ CanonicalForm CanonicalForm::sqrt () const
1286//{{{ docu
1287//
1288// sqrt() - calculate integer square root.
1289//
1290// CO has to be an integer greater or equal zero.  Returns the
1291// largest integer less or equal sqrt(CO).
1292//
1293// In the immediate case, we use the newton method to find the
1294// root.  The algorithm is from H. Cohen - 'A Course in
1295// Computational Algebraic Number Theory', ch. 1.7.1.
1296//
1297// See also: InternalCF::sqrt(), InternalInteger::sqrt(), ::sqrt()
1298//
1299//}}}
1300CanonicalForm
1301CanonicalForm::sqrt () const
1302{
1303    if ( is_imm( value ) ) {
1304        ASSERT( is_imm( value ) == INTMARK, "sqrt() not implemented" );
1305        long n = imm2int( value );
1306        ASSERT( n >= 0, "arg to sqrt() less than zero" );
1307        if ( n == 0 || n == 1 )
1308            return CanonicalForm( n );
1309        else {
1310            long x, y = n;
1311            do {
1312                x = y;
1313                // the intermediate result may not fit into an
1314                // integer, but the result does
1315                y = (unsigned long)(x + n/x)/2;
1316            } while ( y < x );
1317            return CanonicalForm( x );
1318        }
1319    }
1320    else
1321        return CanonicalForm( value->sqrt() );
1322}
1323//}}}
1324
1325//{{{ int CanonicalForm::ilog2 () const
1326//{{{ docu
1327//
1328// ilog2() - integer logarithm to base 2.
1329//
1330// Returns the largest integer less or equal logarithm of CO to
1331// base 2.  CO should be a positive integer.
1332//
1333// See also: InternalCF::ilog2(), InternalInteger::ilog2(), ::ilog2()
1334//
1335//}}}
1336int
1337CanonicalForm::ilog2 () const
1338{
1339    if ( is_imm( value ) )
1340    {
1341        ASSERT( is_imm( value ) == INTMARK, "ilog2() not implemented" );
1342        long a = imm2int( value );
1343        ASSERT( a > 0, "arg to ilog2() less or equal zero" );
1344        int n = -1;
1345        while ( a > 0 )
1346        {
1347          n++;
1348          a /=2;
1349        }
1350        return n;
1351    }
1352    else
1353        return value->ilog2();
1354}
1355//}}}
1356
1357//{{{ bool operator ==, operator != ( const CanonicalForm & lhs, const CanonicalForm & rhs )
1358//{{{ docu
1359//
1360// operator ==(), operator !=() - compare canonical forms on
1361//   (in)equality.
1362//
1363// operator ==() returns true iff lhs equals rhs.
1364// operator !=() returns true iff lhs does not equal rhs.
1365//
1366// This is the point in factory where we essentially use that
1367// CanonicalForms in fact are canonical.  There must not be two
1368// different representations of the same mathematical object,
1369// otherwise, such (in)equality will not be recognized by these
1370// operators.  In other word, we rely on the fact that structural
1371// different factory objects in any case represent different
1372// mathematical objects.
1373//
1374// So we use the following procedure to test on equality (and
1375// analogously on inequality).  First, we check whether lhs.value
1376// equals rhs.value.  If so we are ready and return true.
1377// Second, if one of the operands is immediate, but the other one
1378// not, we return false.  Third, if the operand's levels differ
1379// we return false.  Fourth, if the operand's levelcoeffs differ
1380// we return false.  At last, we call the corresponding internal
1381// method to compare both operands.
1382//
1383// Both operands should have coefficients from the same base domain.
1384//
1385// Note: To compare with the zero or the unit of the current domain,
1386// you better use the methods `CanonicalForm::isZero()' or
1387// `CanonicalForm::isOne()', resp., than something like `f == 0',
1388// since the latter is quite a lot slower.
1389//
1390// See also: InternalCF::comparesame(),
1391// InternalInteger::comparesame(), InternalRational::comparesame(),
1392// InternalPrimePower::comparesame(), InternalPoly::comparesame()
1393//
1394//}}}
1395bool
1396operator == ( const CanonicalForm & lhs, const CanonicalForm & rhs )
1397{
1398    if ( lhs.value == rhs.value )
1399        return true;
1400    else if ( is_imm( rhs.value ) || is_imm( lhs.value ) ) {
1401        ASSERT( ! is_imm( rhs.value ) ||
1402                ! is_imm( lhs.value ) ||
1403                is_imm( rhs.value ) == is_imm( lhs.value ),
1404                "incompatible operands" );
1405        return false;
1406    }
1407    else  if ( lhs.value->level() != rhs.value->level() )
1408        return false;
1409    else  if ( lhs.value->levelcoeff() != rhs.value->levelcoeff() )
1410        return false;
1411    else
1412        return rhs.value->comparesame( lhs.value ) == 0;
1413}
1414
1415bool
1416operator != ( const CanonicalForm & lhs, const CanonicalForm & rhs )
1417{
1418    if ( lhs.value == rhs.value )
1419        return false;
1420    else if ( is_imm( rhs.value ) || is_imm( lhs.value ) ) {
1421        ASSERT( ! is_imm( rhs.value ) ||
1422                ! is_imm( lhs.value ) ||
1423                is_imm( rhs.value ) == is_imm( lhs.value ),
1424                "incompatible operands" );
1425        return true;
1426    }
1427    else  if ( lhs.value->level() != rhs.value->level() )
1428        return true;
1429    else  if ( lhs.value->levelcoeff() != rhs.value->levelcoeff() )
1430        return true;
1431    else        return rhs.value->comparesame( lhs.value ) != 0;
1432}
1433//}}}
1434
1435//{{{ bool operator >, operator < ( const CanonicalForm & lhs, const CanonicalForm & rhs )
1436//{{{ docu
1437//
1438// operator >(), operator <() - compare canonical forms. on size or
1439//   level.
1440//
1441// The most common and most useful application of these operators
1442// is to compare two integers or rationals, of course.  However,
1443// these operators are defined on all other base domains and on
1444// polynomials, too.  From a mathematical point of view this may
1445// seem meaningless, since there is no ordering on finite fields
1446// or on polynomials respecting the algebraic structure.
1447// Nevertheless, from a programmer's point of view it may be
1448// sensible to order these objects, e.g. to sort them.
1449//
1450// Therefore, the ordering defined by these operators in any case
1451// is a total ordering which fulfills the law of trichotomy.
1452//
1453// It is clear how this is done in the case of the integers and
1454// the rationals.  For finite fields, all you can say is that
1455// zero is the minimal element w.r.t. the ordering, the other
1456// elements are ordered in an arbitrary (but total!)  way.  For
1457// polynomials, you have an ordering derived from the
1458// lexicographical ordering of monomials.  E.g. if lm(f) < lm(g)
1459// w.r.t. lexicographic ordering, then f < g.  For more details,
1460// refer to the documentation of `InternalPoly::operator <()'.
1461//
1462// Both operands should have coefficients from the same base domain.
1463//
1464// The scheme how both operators are implemented is allmost the
1465// same as for the assignment operators (check for immediates,
1466// then check levels, then check levelcoeffs, then call the
1467// appropriate internal comparesame()/comparecoeff() method).
1468// For more information, confer to the overview for the
1469// arithmetic operators.
1470//
1471// See also: InternalCF::comparesame(),
1472// InternalInteger::comparesame(), InternalRational::comparesame(),
1473// InternalPrimePower::comparesame(), InternalPoly::comparesame(),
1474// InternalCF::comparecoeff(), InternalInteger::comparecoeff(),
1475// InternalRational::comparecoeff(),
1476// InternalPrimePower::comparecoeff(), InternalPoly::comparecoeff(),
1477// imm_cmp(), imm_cmp_p(), imm_cmp_gf()
1478//
1479//}}}
1480bool
1481operator > ( const CanonicalForm & lhs, const CanonicalForm & rhs )
1482{
1483    int what = is_imm( rhs.value );
1484    if ( is_imm( lhs.value ) ) {
1485        ASSERT( ! what || (what == is_imm( lhs.value )), "incompatible operands" );
1486        if ( what == 0 )
1487            return rhs.value->comparecoeff( lhs.value ) < 0;
1488        else if ( what == INTMARK )
1489            return imm_cmp( lhs.value, rhs.value ) > 0;
1490        else if ( what == FFMARK )
1491            return imm_cmp_p( lhs.value, rhs.value ) > 0;
1492        else
1493            return imm_cmp_gf( lhs.value, rhs.value ) > 0;
1494    }
1495    else  if ( what )
1496        return lhs.value->comparecoeff( rhs.value ) > 0;
1497    else  if ( lhs.value->level() == rhs.value->level() )
1498        if ( lhs.value->levelcoeff() == rhs.value->levelcoeff() )
1499            return lhs.value->comparesame( rhs.value ) > 0;
1500        else  if ( lhs.value->levelcoeff() > rhs.value->levelcoeff() )
1501            return lhs.value->comparecoeff( rhs.value ) > 0;
1502        else
1503            return rhs.value->comparecoeff( lhs.value ) < 0;
1504    else
1505        return lhs.value->level() > rhs.value->level();
1506}
1507
1508bool
1509operator < ( const CanonicalForm & lhs, const CanonicalForm & rhs )
1510{
1511    int what = is_imm( rhs.value );
1512    if ( is_imm( lhs.value ) ) {
1513        ASSERT( ! what || (what == is_imm( lhs.value )), "incompatible operands" );
1514        if ( what == 0 )
1515            return rhs.value->comparecoeff( lhs.value ) > 0;
1516        else if ( what == INTMARK )
1517            return imm_cmp( lhs.value, rhs.value ) < 0;
1518        else if ( what == FFMARK )
1519            return imm_cmp_p( lhs.value, rhs.value ) < 0;
1520        else
1521            return imm_cmp_gf( lhs.value, rhs.value ) < 0;
1522    }
1523    else  if ( what )
1524        return lhs.value->comparecoeff( rhs.value ) < 0;
1525    else  if ( lhs.value->level() == rhs.value->level() )
1526        if ( lhs.value->levelcoeff() == rhs.value->levelcoeff() )
1527            return lhs.value->comparesame( rhs.value ) < 0;
1528        else  if ( lhs.value->levelcoeff() > rhs.value->levelcoeff() )
1529            return lhs.value->comparecoeff( rhs.value ) < 0;
1530        else
1531            return rhs.value->comparecoeff( lhs.value ) > 0;
1532    else
1533        return lhs.value->level() < rhs.value->level();
1534}
1535//}}}
1536
1537//{{{ CanonicalForm bgcd ( const CanonicalForm & f, const CanonicalForm & g )
1538//{{{ docu
1539//
1540// bgcd() - return base coefficient gcd.
1541//
1542// If both f and g are integers and `SW_RATIONAL' is off the
1543// positive greatest common divisor of f and g is returned.
1544// Otherwise, if `SW_RATIONAL' is on or one of f and g is not an
1545// integer, the greatest common divisor is trivial: either zero
1546// if f and g equal zero or one (both from the current domain).
1547//
1548// f and g should come from one base domain which should be not
1549// the prime power domain.
1550//
1551// Implementation:
1552//
1553// CanonicalForm::bgcd() handles the immediate case with a
1554//   standard euclidean algorithm.  For the non-immediate cases
1555//   `InternalCF::bgcdsame()' or `InternalCF::bgcdcoeff()', resp. are
1556//   called following the usual level/levelcoeff approach.
1557//
1558// InternalCF::bgcdsame() and
1559// InternalCF::bgcdcoeff() throw an assertion ("not implemented")
1560//
1561// InternalInteger::bgcdsame() is a wrapper around `mpz_gcd()'
1562//   which takes some care about immediate results and the sign
1563//   of the result
1564// InternalInteger::bgcdcoeff() is a wrapper around
1565//   `mpz_gcd_ui()' which takes some care about the sign
1566//   of the result
1567//
1568// InternalRational::bgcdsame() and
1569// InternalRational::bgcdcoeff() always return one
1570//
1571//}}}
1572CanonicalForm
1573bgcd ( const CanonicalForm & f, const CanonicalForm & g )
1574{
1575    // check immediate cases
1576    int what = is_imm( g.value );
1577    if ( is_imm( f.value ) )
1578    {
1579        ASSERT( ! what || (what == is_imm( f.value )), "incompatible operands" );
1580        if ( what == 0 )
1581            return g.value->bgcdcoeff( f.value );
1582        else if ( what == INTMARK && ! cf_glob_switches.isOn( SW_RATIONAL ) )
1583        {
1584            // calculate gcd using standard integer
1585            // arithmetic
1586            long fInt = imm2int( f.value );
1587            long gInt = imm2int( g.value );
1588
1589            if ( fInt < 0 ) fInt = -fInt;
1590            if ( gInt < 0 ) gInt = -gInt;
1591            // swap fInt and gInt
1592            if ( gInt > fInt )
1593            {
1594                long swap = gInt;
1595                gInt = fInt;
1596                fInt = swap;
1597            }
1598
1599            // now, 0 <= gInt <= fInt.  Start the loop.
1600            while ( gInt )
1601            {
1602                // calculate (fInt, gInt) = (gInt, fInt%gInt)
1603                long r = fInt % gInt;
1604                fInt = gInt;
1605                gInt = r;
1606            }
1607
1608            return CanonicalForm( fInt );
1609        }
1610        else
1611            // we do not go for maximal speed for these stupid
1612            // special cases
1613            return CanonicalForm( f.isZero() && g.isZero() ? 0 : 1 );
1614    }
1615    else if ( what )
1616        return f.value->bgcdcoeff( g.value );
1617
1618    int fLevel = f.value->level();
1619    int gLevel = g.value->level();
1620
1621    // check levels
1622    if ( fLevel == gLevel )
1623    {
1624        fLevel = f.value->levelcoeff();
1625        gLevel = g.value->levelcoeff();
1626
1627        // check levelcoeffs
1628        if ( fLevel == gLevel )
1629            return f.value->bgcdsame( g.value );
1630        else if ( fLevel < gLevel )
1631            return g.value->bgcdcoeff( f.value );
1632        else
1633            return f.value->bgcdcoeff( g.value );
1634    }
1635    else if ( fLevel < gLevel )
1636        return g.value->bgcdcoeff( f.value );
1637    else
1638        return f.value->bgcdcoeff( g.value );
1639}
1640//}}}
1641
1642//{{{ CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
1643//{{{ docu
1644//
1645// bextgcd() - return base coefficient extended gcd.
1646//
1647//}}}
1648CanonicalForm
1649bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
1650{
1651    // check immediate cases
1652    int what = is_imm( g.value );
1653    if ( is_imm( f.value ) ) {
1654        ASSERT( ! what || (what == is_imm( f.value )), "incompatible operands" );
1655        if ( what == 0 )
1656            return g.value->bextgcdcoeff( f.value, b, a );
1657        else if ( what == INTMARK && ! cf_glob_switches.isOn( SW_RATIONAL ) ) {
1658            // calculate extended gcd using standard integer
1659            // arithmetic
1660            long fInt = imm2int( f.value );
1661            long gInt = imm2int( g.value );
1662
1663            // to avoid any system dpendencies with `%', we work
1664            // with positive numbers only.  To a pity, we have to
1665            // redo all the checks when assigning to a and b.
1666            if ( fInt < 0 ) fInt = -fInt;
1667            if ( gInt < 0 ) gInt = -gInt;
1668            // swap fInt and gInt
1669            if ( gInt > fInt ) {
1670                long swap = gInt;
1671                gInt = fInt;
1672                fInt = swap;
1673            }
1674
1675            long u = 1; long v = 0;
1676            long uNext = 0; long vNext = 1;
1677
1678            // at any step, we have:
1679            // fInt_0 * u + gInt_0 * v = fInt
1680            // fInt_0 * uNext + gInt_0 * vNext = gInt
1681            // where fInt_0 and gInt_0 denote the values of fint
1682            // and gInt, resp., at the beginning
1683            while ( gInt ) {
1684                long r = fInt % gInt;
1685                long q = fInt / gInt;
1686                long uSwap = u - q * uNext;
1687                long vSwap = v - q * vNext;
1688
1689                // update variables
1690                fInt = gInt;
1691                gInt = r;
1692                u = uNext; v = vNext;
1693                uNext = uSwap; vNext = vSwap;
1694            }
1695
1696            // now, assign to a and b
1697            long fTest = imm2int( f.value );
1698            long gTest = imm2int( g.value );
1699            if ( gTest > fTest ) {
1700                a = v; b = u;
1701            } else {
1702                a = u; b = v;
1703            }
1704            if ( fTest < 0 ) a = -a;
1705            if ( gTest < 0 ) b = -b;
1706            return CanonicalForm( fInt );
1707        } else
1708            // stupid special cases
1709            if ( ! f.isZero() ) {
1710                a = 1/f; b = 0; return CanonicalForm( 1 );
1711            } else if ( ! g.isZero() ) {
1712                a = 0; b = 1/g; return CanonicalForm( 1 );
1713            } else {
1714                a = 0; b = 0; return CanonicalForm( 0 );
1715            }
1716    }
1717    else if ( what )
1718        return f.value->bextgcdcoeff( g.value, a, b );
1719
1720    int fLevel = f.value->level();
1721    int gLevel = g.value->level();
1722
1723    // check levels
1724    if ( fLevel == gLevel ) {
1725        fLevel = f.value->levelcoeff();
1726        gLevel = g.value->levelcoeff();
1727
1728        // check levelcoeffs
1729        if ( fLevel == gLevel )
1730            return f.value->bextgcdsame( g.value, a, b );
1731        else if ( fLevel < gLevel )
1732            return g.value->bextgcdcoeff( f.value, b, a );
1733        else
1734            return f.value->bextgcdcoeff( g.value, a, b );
1735    }
1736    else if ( fLevel < gLevel )
1737        return g.value->bextgcdcoeff( f.value, b, a );
1738    else
1739        return f.value->bextgcdcoeff( g.value, a, b );
1740}
1741//}}}
1742
1743CanonicalForm
1744blcm ( const CanonicalForm & f, const CanonicalForm & g )
1745{
1746    if ( f.isZero() || g.isZero() )
1747        return CanonicalForm( 0 );
1748/*
1749    else if (f.isOne())
1750        return g;
1751    else if (g.isOne())
1752        return f;
1753*/
1754    else
1755        return (f / bgcd( f, g )) * g;
1756}
1757
1758//{{{ input/output
1759#ifndef NOSTREAMIO
1760void
1761CanonicalForm::print( OSTREAM & os, char * str ) const
1762{
1763    if ( is_imm( value ) )
1764        imm_print( os, value, str );
1765    else
1766        value->print( os, str );
1767}
1768
1769void
1770CanonicalForm::print( OSTREAM & os ) const
1771{
1772    if ( is_imm( value ) )
1773        imm_print( os, value, "" );
1774    else
1775        value->print( os, "" );
1776}
1777
1778OSTREAM&
1779operator << ( OSTREAM & os, const CanonicalForm & cf )
1780{
1781    cf.print( os, "" );
1782    return os;
1783}
1784
1785ISTREAM&
1786operator >> ( ISTREAM & is, CanonicalForm & cf )
1787{
1788    cf = readCF( is );
1789    return is;
1790}
1791#endif /* NOSTREAMIO */
1792//}}}
1793
1794//{{{ genOne(), genZero()
1795CanonicalForm
1796CanonicalForm::genZero() const
1797{
1798    int what = is_imm( value );
1799    if ( what == FFMARK )
1800        return CanonicalForm( CFFactory::basic( FiniteFieldDomain, 0L ) );
1801    else  if ( what == GFMARK )
1802        return CanonicalForm( CFFactory::basic( GaloisFieldDomain, 0L ) );
1803    else  if ( what )
1804        return CanonicalForm( CFFactory::basic( IntegerDomain, 0L ) );
1805    else
1806        return CanonicalForm( value->genZero() );
1807}
1808
1809CanonicalForm
1810CanonicalForm::genOne() const
1811{
1812    int what = is_imm( value );
1813    if ( what == FFMARK )
1814        return CanonicalForm( CFFactory::basic( FiniteFieldDomain, 1L ) );
1815    else  if ( what == GFMARK )
1816        return CanonicalForm( CFFactory::basic( GaloisFieldDomain, 1L ) );
1817    else  if ( what )
1818        return CanonicalForm( CFFactory::basic( IntegerDomain, 1L ) );
1819    else
1820        return CanonicalForm( value->genOne() );
1821}
1822//}}}
1823
1824//{{{ exponentiation
1825CanonicalForm
1826power ( const CanonicalForm & f, int n )
1827{
1828  ASSERT( n >= 0, "illegal exponent" );
1829  if ( f.isZero() )
1830    return 0;
1831  else  if ( f.isOne() )
1832    return f;
1833  else  if ( f == -1 )
1834  {
1835    if ( n % 2 == 0 )
1836      return 1;
1837    else
1838      return -1;
1839  }
1840  else  if ( n == 0 )
1841    return 1;
1842
1843  //else if (f.inGF())
1844  //{
1845  //}
1846  else
1847  {
1848    CanonicalForm g,h;
1849    h=f;
1850    while(n%2==0)
1851    {
1852      h*=h;
1853      n/=2;
1854    }
1855    g=h;
1856    while(1)
1857    {
1858      n/=2;
1859      if(n==0)
1860        return g;
1861      h*=h;
1862      if(n%2!=0) g*=h;
1863    }
1864  }
1865}
1866
1867CanonicalForm
1868power ( const Variable & v, int n )
1869{
1870    //ASSERT( n >= 0, "illegal exponent" );
1871    if ( n == 0 )
1872        return 1;
1873    else  if ( n == 1 )
1874        return v;
1875    else  if (( v.level() < 0 ) && (hasMipo(v)))
1876    {
1877        CanonicalForm result( v, n-1 );
1878        return result * v;
1879    }
1880    else
1881        return CanonicalForm( v, n );
1882}
1883//}}}
1884
1885//{{{ switches
1886void
1887On( int sw )
1888{
1889    cf_glob_switches.On( sw );
1890}
1891
1892void
1893Off( int sw )
1894{
1895    cf_glob_switches.Off( sw );
1896}
1897
1898bool
1899isOn( int sw )
1900{
1901    return cf_glob_switches.isOn( sw );
1902}
1903//}}}
Note: See TracBrowser for help on using the repository browser.