source: git/factory/canonicalform.cc @ 9f7665

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