source: git/factory/canonicalform.cc @ 4deddb

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