1 /* emacs edit mode for this file is -*- C++ -*- */
2
3 /**
4 *
5 * @file cf_ops.cc
6 *
7 * simple structural algorithms.
8 *
9 * A 'structural' algorithm is an algorithm which gives
10 * structural information on polynomials in contrast to a
11 * 'mathematical' algorithm which calculates some mathematical
12 * function.
13 *
14 * Compare these functions with the functions in cf_algorithm.cc,
15 * which are mathematical algorithms.
16 *
17 *
18 * Header file: canonicalform.h
19 *
20 **/
21
22
23 #include "config.h"
24
25
26 #include "cf_assert.h"
27
28 #include "canonicalform.h"
29 #include "variable.h"
30 #include "cf_iter.h"
31
32 /** static Variable sv_x1, sv_x2;
33 *
34 * sv_x1, sv_x2 - variables to swap by swapvar() and replacevar.
35 *
36 * These variables are initialized by swapvar() such that sv_x1 <
37 * sv_x2. They are used by swapvar_between() and swapvar_rec()
38 * to swap variables efficiently.
39 * Furthermore, sv_x1 and sv_x2 are used by replacevar() and
40 * replacevar_between().
41 *
42 **/
43 STATIC_INST_VAR Variable sv_x1, sv_x2;
44
45 /** static void swapvar_between ( const CanonicalForm & f, CanonicalForm & result, const CanonicalForm & term, int expx2 )
46 *
47 * swapvar_between() - replace occurences of sv_x1 in f with sv_x2.
48 *
49 * If Psi denotes the map which maps sv_x1 to sv_x2, this
50 * function returns
51 *
52 * result + Psi(f) * term * sv_x1^expx2
53 *
54 * Used by: swapvar()
55 *
56 **/
57 static void
swapvar_between(const CanonicalForm & f,CanonicalForm & result,const CanonicalForm & term,int expx2)58 swapvar_between ( const CanonicalForm & f, CanonicalForm & result, const CanonicalForm & term, int expx2 )
59 {
60 if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
61 // in this case, we do not have to replace anything
62 result += term * power( sv_x1, expx2 ) * f;
63 else if ( f.mvar() == sv_x1 )
64 // this is where the real work is done: this iterator
65 // replaces sv_x1 with sv_x2
66 for ( CFIterator i = f; i.hasTerms(); i++ )
67 result += power( sv_x2, i.exp() ) * term * power( sv_x1, expx2 ) * i.coeff();
68 else
69 // f's level is larger than sv_x1: descend down
70 for ( CFIterator i = f; i.hasTerms(); i++ )
71 swapvar_between( i.coeff(), result, term * power( f.mvar(), i.exp() ), expx2 );
72 }
73 #if 0
74 static CanonicalForm
75 swapvar_between1 ( const CanonicalForm & f )
76 {
77 if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
78 // in this case, we do not have to replace anything
79 return f;
80 else if ( f.mvar() == sv_x1 )
81 {
82 // this is where the real work is done: this iterator
83 // replaces sv_x1 with sv_x2
84 CanonicalForm result;
85 for ( CFIterator i = f; i.hasTerms(); i++ )
86 result += power( sv_x2, i.exp() ) * i.coeff();
87 return result;
88 }
89 else
90 {
91 // f's level is larger than sv_x1: descend down
92 CanonicalForm result;
93 for ( CFIterator i = f; i.hasTerms(); i++ )
94 result += swapvar_between1( i.coeff() ) * power( f.mvar(), i.exp() );
95 return result;
96 }
97 }
98 #endif
99
100 /**
101 *
102 * swapvar_between() - swap occurences of sv_x1 and sv_x2 in f.
103 *
104 * If Psi denotes the map which swaps sv_x1 and sv_x2, this
105 * function returns
106 *
107 * result + Psi(f) * term
108 *
109 * Used by: swapvar()
110 *
111 **/
112 static void
swapvar_rec(const CanonicalForm & f,CanonicalForm & result,const CanonicalForm & term)113 swapvar_rec ( const CanonicalForm & f, CanonicalForm & result, const CanonicalForm & term )
114 {
115 if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
116 // in this case, we do not have to swap anything
117 result += term * f;
118 else if ( f.mvar() == sv_x2 )
119 // this is where the real work is done: this iterator
120 // replaces sv_x1 with sv_x2 in the coefficients of f and
121 // remembers the exponents of sv_x2 in the last argument
122 // of the call to swapvar_between()
123 for ( CFIterator i = f; i.hasTerms(); i++ )
124 swapvar_between( i.coeff(), result, term, i.exp() );
125 else if ( f.mvar() < sv_x2 )
126 // sv_x2 does not occur in f, but sv_x1 does. Replace it.
127 swapvar_between( f, result, term, 0 );
128 else
129 // f's level is larger than sv_x2: descend down
130 for ( CFIterator i = f; i.hasTerms(); i++ )
131 swapvar_rec( i.coeff(), result, term * power( f.mvar(), i.exp() ) );
132 }
133 #if 0
134 static CanonicalForm
135 swapvar_rec1 ( const CanonicalForm & f )
136 {
137 if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
138 return f;
139 else if ( f.mvar() == sv_x2 )
140 {
141 CanonicalForm result;
142 for ( CFIterator i = f; i.hasTerms(); i++ )
143 result += swapvar_between1( i.coeff() ) * power( sv_x1, i.exp() );
144 return result;
145 }
146 else if ( f.mvar() < sv_x2 )
147 return swapvar_between1( f );
148 else
149 {
150 CanonicalForm result;
151 for ( CFIterator i = f; i.hasTerms(); i++ )
152 result += swapvar_rec1( i.coeff() ) * power( f.mvar(), i.exp() );
153 return result;
154 }
155 }
156 #endif
157
158 /**
159 *
160 * swapvar() - swap variables x1 and x2 in f.
161 *
162 * Returns the image of f under the map which maps x1 to x2 and
163 * x2 to x1. This is done quite efficiently because it is used
164 * really often. x1 and x2 should be polynomial variables.
165 *
166 **/
167 CanonicalForm
swapvar(const CanonicalForm & f,const Variable & x1,const Variable & x2)168 swapvar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
169 {
170 ASSERT( x1.level() > 0 && x2.level() > 0, "cannot swap algebraic Variables" );
171 if ( f.inCoeffDomain() || x1 == x2 || ( x1 > f.mvar() && x2 > f.mvar() ) )
172 return f;
173 else
174 {
175 CanonicalForm result = 0;
176 if ( x1 > x2 )
177 {
178 sv_x1 = x2; sv_x2 = x1;
179 }
180 else
181 {
182 sv_x1 = x1; sv_x2 = x2;
183 }
184 if ( f.mvar() < sv_x2 )
185 // we only have to replace sv_x1 by sv_x2
186 swapvar_between( f, result, 1, 0 );
187 else
188 // we really have to swap variables
189 swapvar_rec( f, result, 1 );
190 return result;
191 }
192 }
193 #if 0
194 CanonicalForm
195 swapvar1 ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
196 {
197 ASSERT( x1.level() > 0 && x2.level() > 0, "cannot swap algebraic variables" );
198 if ( f.inCoeffDomain() || x1 == x2 || ( x1 > f.mvar() && x2 > f.mvar() ) )
199 return f;
200 else
201 {
202 CanonicalForm result = 0;
203 if ( x1 > x2 ) {
204 sv_x1 = x2; sv_x2 = x1;
205 }
206 else
207 {
208 sv_x1 = x1; sv_x2 = x2;
209 }
210 if ( f.mvar() < sv_x2 )
211 // we only have to replace sv_x1 by sv_x2
212 return swapvar_between1( f );
213 else
214 // we really have to swap variables
215 return swapvar_rec1( f );
216 }
217 }
218 #endif
219
220 /**
221 *
222 * replacevar_between() - replace occurences of sv_x1 in f with sv_x2.
223 *
224 * This is allmost the same as swapvar_between() except that
225 * sv_x1 may be an algebraic variable, so we have to test on
226 * 'f.inBaseDomain()' instead of 'f.inCoeffDomain()' in the
227 * beginning.
228 *
229 * Used by: replacevar()
230 *
231 **/
232 static CanonicalForm
replacevar_between(const CanonicalForm & f)233 replacevar_between ( const CanonicalForm & f )
234 {
235 if ( f.inBaseDomain() )
236 return f;
237
238 Variable x = f.mvar();
239
240 if ( x < sv_x1 )
241 // in this case, we do not have to replace anything
242 return f;
243 else if ( x == sv_x1 )
244 {
245 // this is where the real work is done: this iterator
246 // replaces sv_x1 with sv_x2
247 CanonicalForm result;
248 for ( CFIterator i = f; i.hasTerms(); i++ )
249 result += power( sv_x2, i.exp() ) * i.coeff();
250 return result;
251 }
252 else
253 {
254 // f's level is larger than sv_x1: descend down
255 CanonicalForm result;
256 for ( CFIterator i = f; i.hasTerms(); i++ )
257 result += replacevar_between( i.coeff() ) * power( x, i.exp() );
258 return result;
259 }
260 }
261
262 /** CanonicalForm replacevar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
263 *
264 * replacevar() - replace all occurences of x1 in f by x2.
265 *
266 * In contrast to swapvar(), x1 may be an algebraic variable, but
267 * x2 must be a polynomial variable.
268 *
269 **/
270 CanonicalForm
replacevar(const CanonicalForm & f,const Variable & x1,const Variable & x2)271 replacevar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
272 {
273 //ASSERT( x2.level() > 0, "cannot replace with algebraic variable" );
274 if ( f.inBaseDomain() || x1 == x2 || ( x1 > f.mvar() ) )
275 return f;
276 else
277 {
278 sv_x1 = x1;
279 sv_x2 = x2;
280 return replacevar_between( f );
281 }
282 }
283
284 /** static void fillVarsRec ( const CanonicalForm & f, int * vars )
285 *
286 * fillVarsRec - fill array describing occurences of variables in f.
287 *
288 * Only polynomial variables are looked up. The information is
289 * stored in the arrary vars. vars should be large enough to
290 * hold all information, i.e. larger than the level of f.
291 *
292 * Used by getVars() and getNumVars().
293 *
294 **/
295 static void
fillVarsRec(const CanonicalForm & f,int * vars)296 fillVarsRec ( const CanonicalForm & f, int * vars )
297 {
298 int n;
299 if ( (n = f.level()) > 0 )
300 {
301 vars[n] = 1;
302 CFIterator i;
303 for ( i = f; i.hasTerms(); ++i )
304 fillVarsRec( i.coeff(), vars );
305 }
306 }
307
308 /** int getNumVars ( const CanonicalForm & f )
309 *
310 * getNumVars() - get number of polynomial variables in f.
311 *
312 **/
313 int
getNumVars(const CanonicalForm & f)314 getNumVars ( const CanonicalForm & f )
315 {
316 int n;
317 if ( f.inCoeffDomain() )
318 return 0;
319 else if ( (n = f.level()) == 1 )
320 return 1;
321 else
322 {
323 int * vars = NEW_ARRAY(int, n+1);
324 int i;
325 for ( i = n-1; i >=0; i-- ) vars[i] = 0;
326
327 // look for variables
328 for ( CFIterator I = f; I.hasTerms(); ++I )
329 fillVarsRec( I.coeff(), vars );
330
331 // count them
332 int m = 0;
333 for ( i = 1; i < n; i++ )
334 if ( vars[i] != 0 ) m++;
335
336 DELETE_ARRAY(vars);
337 // do not forget to count our own variable
338 return m+1;
339 }
340 }
341
342 /** CanonicalForm getVars ( const CanonicalForm & f )
343 *
344 * getVars() - get polynomial variables of f.
345 *
346 * Return the product of all of them, 1 if there are not any.
347 *
348 **/
349 CanonicalForm
getVars(const CanonicalForm & f)350 getVars ( const CanonicalForm & f )
351 {
352 int n;
353 if ( f.inCoeffDomain() )
354 return 1;
355 else if ( (n = f.level()) == 1 )
356 return Variable( 1 );
357 else
358 {
359 int * vars = NEW_ARRAY(int, n+1);
360 int i;
361 for ( i = n; i >= 0; i-- ) vars[i] = 0;
362
363 // look for variables
364 for ( CFIterator I = f; I.hasTerms(); ++I )
365 fillVarsRec( I.coeff(), vars );
366
367 // multiply them all
368 CanonicalForm result = 1;
369 for ( i = n; i > 0; i-- )
370 if ( vars[i] != 0 ) result *= Variable( i );
371
372 DELETE_ARRAY(vars);
373 // do not forget our own variable
374 return f.mvar() * result;
375 }
376 }
377
378 /** CanonicalForm apply ( const CanonicalForm & f, void (*mf)( CanonicalForm &, int & ) )
379 *
380 * apply() - apply mf to terms of f.
381 *
382 * Calls mf( f[i], i ) for each term f[i]*x^i of f and builds a
383 * new term from the result. If f is in a coefficient domain,
384 * mf( f, i ) should result in an i == 0, since otherwise it is
385 * not clear which variable to use for the resulting term.
386 *
387 * An example:
388 *
389 ~~~~~~~~~~~~~~~~~~~~~{.c}
390 void
391 diff( CanonicalForm & f, int & i )
392 {
393 f = f * i;
394 if ( i > 0 ) i--;
395 }
396 ~~~~~~~~~~~~~~~~~~~~~
397 * Then apply( f, diff ) is differentation of f with respect to the
398 * main variable of f.
399 *
400 **/
401 CanonicalForm
apply(const CanonicalForm & f,void (* mf)(CanonicalForm &,int &))402 apply ( const CanonicalForm & f, void (*mf)( CanonicalForm &, int & ) )
403 {
404 if ( f.inCoeffDomain() )
405 {
406 int exp = 0;
407 CanonicalForm result = f;
408 mf( result, exp );
409 ASSERT( exp == 0, "illegal result, do not know what variable to use" );
410 return result;
411 }
412 else
413 {
414 CanonicalForm result, coeff;
415 CFIterator i;
416 int exp;
417 Variable x = f.mvar();
418 for ( i = f; i.hasTerms(); i++ )
419 {
420 coeff = i.coeff();
421 exp = i.exp();
422 mf( coeff, exp );
423 if ( ! coeff.isZero() )
424 result += power( x, exp ) * coeff;
425 }
426 return result;
427 }
428 }
429
430 /** CanonicalForm mapdomain ( const CanonicalForm & f, CanonicalForm (*mf)( const CanonicalForm & ) )
431 *
432 * mapdomain() - map all coefficients of f through mf.
433 *
434 * Recursively descends down through f to the coefficients which
435 * are in a coefficient domain mapping each such coefficient
436 * through mf and returns the result.
437 *
438 **/
439 CanonicalForm
mapdomain(const CanonicalForm & f,CanonicalForm (* mf)(const CanonicalForm &))440 mapdomain ( const CanonicalForm & f, CanonicalForm (*mf)( const CanonicalForm & ) )
441 {
442 if ( f.inBaseDomain() )
443 return mf( f );
444 else
445 {
446 CanonicalForm result = 0;
447 CFIterator i;
448 Variable x = f.mvar();
449 for ( i = f; i.hasTerms(); i++ )
450 result += power( x, i.exp() ) * mapdomain( i.coeff(), mf );
451 return result;
452 }
453 }
454
455 /** static void degreesRec ( const CanonicalForm & f, int * degs )
456 *
457 * degreesRec() - recursively get degrees of f.
458 *
459 * Used by degrees().
460 *
461 **/
462 static void
degreesRec(const CanonicalForm & f,int * degs)463 degreesRec ( const CanonicalForm & f, int * degs )
464 {
465 if ( ! f.inCoeffDomain() )
466 {
467 int level = f.level();
468 int deg = f.degree();
469 // calculate the maximum degree of all coefficients which
470 // are in the same level
471 if ( degs[level] < deg )
472 degs[level] = f.degree();
473 for ( CFIterator i = f; i.hasTerms(); i++ )
474 degreesRec( i.coeff(), degs );
475 }
476 }
477
478 /** int * degrees ( const CanonicalForm & f, int * degs )
479 *
480 * degress() - return the degrees of all polynomial variables in f.
481 *
482 * Returns 0 if f is in a coefficient domain, the degrees of f in
483 * all its polynomial variables in an array of int otherwise:
484 *
485 * degrees( f, 0 )[i] = degree( f, Variable(i) )
486 *
487 * If degs is not the zero pointer the degrees are stored in this
488 * array. In this case degs should be larger than the level of
489 * f. If degs is the zero pointer, an array of sufficient size
490 * is allocated automatically.
491 *
492 **/
degrees(const CanonicalForm & f,int * degs)493 int * degrees ( const CanonicalForm & f, int * degs )
494 {
495 if ( f.inCoeffDomain() )
496 {
497 if (degs != 0)
498 return degs;
499 else
500 return 0;
501 }
502 else
503 {
504 int level = f.level();
505 if ( degs == NULL )
506 degs = NEW_ARRAY(int,level+1);
507 for ( int i = level; i >= 0; i-- )
508 degs[i] = 0;
509 degreesRec( f, degs );
510 return degs;
511 }
512 }
513
514 /** int totaldegree ( const CanonicalForm & f )
515 *
516 * totaldegree() - return the total degree of f.
517 *
518 * If f is zero, return -1. If f is in a coefficient domain,
519 * return 0. Otherwise return the total degree of f in all
520 * polynomial variables.
521 *
522 **/
totaldegree(const CanonicalForm & f)523 int totaldegree ( const CanonicalForm & f )
524 {
525 if ( f.isZero() )
526 return -1;
527 else if ( f.inCoeffDomain() )
528 return 0;
529 else
530 {
531 CFIterator i;
532 int cdeg = 0, dummy;
533 // calculate maximum over all coefficients of f, taking
534 // in account our own exponent
535 for ( i = f; i.hasTerms(); i++ )
536 if ( (dummy = totaldegree( i.coeff() ) + i.exp()) > cdeg )
537 cdeg = dummy;
538 return cdeg;
539 }
540 }
541
542 /** int totaldegree ( const CanonicalForm & f, const Variable & v1, const Variable & v2 )
543 *
544 * totaldegree() - return the total degree of f as a polynomial
545 * in the polynomial variables between v1 and v2 (inclusively).
546 *
547 * If f is zero, return -1. If f is in a coefficient domain,
548 * return 0. Also, return 0 if v1 > v2. Otherwise, take f to be
549 * a polynomial in the polynomial variables between v1 and v2 and
550 * return its total degree.
551 *
552 **/
553 int
totaldegree(const CanonicalForm & f,const Variable & v1,const Variable & v2)554 totaldegree ( const CanonicalForm & f, const Variable & v1, const Variable & v2 )
555 {
556 if ( f.isZero() )
557 return -1;
558 else if ( v1 > v2 )
559 return 0;
560 else if ( f.inCoeffDomain() )
561 return 0;
562 else if ( f.mvar() < v1 )
563 return 0;
564 else if ( f.mvar() == v1 )
565 return f.degree();
566 else if ( f.mvar() > v2 )
567 {
568 // v2's level is larger than f's level, descend down
569 CFIterator i;
570 int cdeg = 0, dummy;
571 // calculate maximum over all coefficients of f
572 for ( i = f; i.hasTerms(); i++ )
573 if ( (dummy = totaldegree( i.coeff(), v1, v2 )) > cdeg )
574 cdeg = dummy;
575 return cdeg;
576 }
577 else
578 {
579 // v1 < f.mvar() <= v2
580 CFIterator i;
581 int cdeg = 0, dummy;
582 // calculate maximum over all coefficients of f, taking
583 // in account our own exponent
584 for ( i = f; i.hasTerms(); i++ )
585 if ( (dummy = totaldegree( i.coeff(), v1, v2 ) + i.exp()) > cdeg )
586 cdeg = dummy;
587 return cdeg;
588 }
589 }
590
591 /** int size ( const CanonicalForm & f, const Variable & v )
592 *
593 * size() - count number of monomials of f with level higher
594 * or equal than level of v.
595 *
596 * Returns one if f is in an base domain.
597 *
598 **/
599 int
size(const CanonicalForm & f,const Variable & v)600 size ( const CanonicalForm & f, const Variable & v )
601 {
602 if ( f.inBaseDomain() )
603 return 1;
604
605 if ( f.mvar() < v )
606 // polynomials with level < v1 are counted as coefficients
607 return 1;
608 else
609 {
610 CFIterator i;
611 int result = 0;
612 // polynomials with level > v2 are not counted al all
613 for ( i = f; i.hasTerms(); i++ )
614 result += size( i.coeff(), v );
615 return result;
616 }
617 }
618
619 /** int size ( const CanonicalForm & f )
620 *
621 * size() - return number of monomials in f which are in an
622 * coefficient domain.
623 *
624 * Returns one if f is in an coefficient domain.
625 *
626 **/
size(const CanonicalForm & f)627 int size ( const CanonicalForm & f )
628 {
629 if ( f.inCoeffDomain() )
630 return 1;
631 else
632 {
633 int result = 0;
634 CFIterator i;
635 for ( i = f; i.hasTerms(); i++ )
636 result += size( i.coeff() );
637 return result;
638 }
639 }
640
size_maxexp(const CanonicalForm & f,int & maxexp)641 int size_maxexp ( const CanonicalForm & f, int& maxexp )
642 {
643 if ( f.inCoeffDomain() )
644 return 1;
645 else
646 {
647 if (f.degree()>maxexp) maxexp=f.degree();
648 int result = 0;
649 CFIterator i;
650 for ( i = f; i.hasTerms(); i++ )
651 result += size_maxexp( i.coeff(), maxexp );
652 return result;
653 }
654 }
655
656 /** polynomials in M.mvar() are considered coefficients
657 * M univariate monic polynomial
658 * the coefficients of f are reduced modulo M
659 **/
reduce(const CanonicalForm & f,const CanonicalForm & M)660 CanonicalForm reduce(const CanonicalForm & f, const CanonicalForm & M)
661 {
662 if(f.inBaseDomain() || f.level() < M.level())
663 return f;
664 if(f.level() == M.level())
665 {
666 if(f.degree() < M.degree())
667 return f;
668 CanonicalForm tmp = mod (f, M);
669 return tmp;
670 }
671 // here: f.level() > M.level()
672 CanonicalForm result = 0;
673 for(CFIterator i=f; i.hasTerms(); i++)
674 result += reduce(i.coeff(),M) * power(f.mvar(),i.exp());
675 return result;
676 }
677
678 /** check if poly f contains an algebraic variable a **/
hasFirstAlgVar(const CanonicalForm & f,Variable & a)679 bool hasFirstAlgVar( const CanonicalForm & f, Variable & a )
680 {
681 if( f.inBaseDomain() ) // f has NO alg. variable
682 return false;
683 if( f.level()<0 ) // f has only alg. vars, so take the first one
684 {
685 a = f.mvar();
686 return true;
687 }
688 for(CFIterator i=f; i.hasTerms(); i++)
689 if( hasFirstAlgVar( i.coeff(), a ))
690 return true; // 'a' is already set
691 return false;
692 }
693
694 /** left shift the main variable of F by n
695 * @return if x is the main variable of F the result is F(x^n)
696 **/
leftShift(const CanonicalForm & F,int n)697 CanonicalForm leftShift (const CanonicalForm& F, int n)
698 {
699 ASSERT (n >= 0, "cannot left shift by negative number");
700 if (F.inBaseDomain())
701 return F;
702 if (n == 0)
703 return F;
704 Variable x=F.mvar();
705 CanonicalForm result= 0;
706 for (CFIterator i= F; i.hasTerms(); i++)
707 result += i.coeff()*power (x, i.exp()*n);
708 return result;
709 }
710