source: git/factory/canonicalform.cc @ 43568c

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