source: git/factory/canonicalform.cc @ e3198f

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