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