source: git/factory/canonicalform.cc @ 362fc67

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