source: git/factory/canonicalform.cc @ 341696

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