source: git/factory/cf_ops.cc @ 85bcd6

spielwiese
Last change on this file since 85bcd6 was 362fc67, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: remove $Id$
  • Property mode set to 100644
File size: 18.1 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3//{{{ docu
4//
5// cf_ops.cc - simple structural algorithms.
6//
7// A 'structural' algorithm is an algorithm which gives
8// structural information on polynomials in contrast to a
9// 'mathematical' algorithm which calculates some mathematical
10// function.
11//
12// Compare these functions with the functions in cf_algorithm.cc,
13// which are mathematical algorithms.
14//
15// Used by: allmost everywhere
16//
17// Header file: canonicalform.h
18//
19//}}}
20
21#include "config.h"
22
23#include "cf_assert.h"
24
25#include "canonicalform.h"
26#include "variable.h"
27#include "cf_iter.h"
28
29//{{{ static Variable sv_x1, sv_x2;
30//{{{ docu
31//
32// sv_x1, sv_x2 - variables to swap by swapvar() and replacevar.
33//
34// These variables are initialized by swapvar() such that sv_x1 <
35// sv_x2.  They are used by swapvar_between() and swapvar_rec()
36// to swap variables efficiently.
37// Furthermore, sv_x1 and sv_x2 are used by replacevar() and
38// replacevar_between().
39//
40//}}}
41static Variable sv_x1, sv_x2;
42//}}}
43
44//{{{ static void swapvar_between ( const CanonicalForm & f, CanonicalForm & result, const CanonicalForm & term, int expx2 )
45//{{{ docu
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//}}}
57static void
58swapvar_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}
73static CanonicalForm
74swapvar_between1 ( const CanonicalForm & f )
75{
76    if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
77        // in this case, we do not have to replace anything
78        return f;
79    else  if ( f.mvar() == sv_x1 )
80    {
81        // this is where the real work is done: this iterator
82        // replaces sv_x1 with sv_x2
83        CanonicalForm result;
84        for ( CFIterator i = f; i.hasTerms(); i++ )
85            result += power( sv_x2, i.exp() ) * i.coeff();
86        return result;
87    }
88    else
89    {
90        // f's level is larger than sv_x1: descend down
91        CanonicalForm result;
92        for ( CFIterator i = f; i.hasTerms(); i++ )
93            result += swapvar_between1( i.coeff() ) * power( f.mvar(), i.exp() );
94        return result;
95    }
96}
97//}}}
98
99//{{{ static void swapvar_rec ( const CanonicalForm & f, CanonicalForm & result, const CanonicalForm & term )
100//{{{ docu
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//}}}
112static void
113swapvar_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}
133static CanonicalForm
134swapvar_rec1 ( const CanonicalForm & f )
135{
136    if ( f.inCoeffDomain() || f.mvar() < sv_x1 )
137        return f;
138    else  if ( f.mvar() == sv_x2 )
139    {
140        CanonicalForm result;
141        for ( CFIterator i = f; i.hasTerms(); i++ )
142            result += swapvar_between1( i.coeff() ) * power( sv_x1, i.exp() );
143        return result;
144    }
145    else  if ( f.mvar() < sv_x2 )
146        return swapvar_between1( f );
147    else
148    {
149        CanonicalForm result;
150        for ( CFIterator i = f; i.hasTerms(); i++ )
151            result += swapvar_rec1( i.coeff() ) * power( f.mvar(), i.exp() );
152        return result;
153    }
154}
155//}}}
156
157//{{{ CanonicalForm swapvar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
158//{{{ docu
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//}}}
167CanonicalForm
168swapvar ( 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}
193CanonicalForm
194swapvar1 ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
195{
196    ASSERT( x1.level() > 0 && x2.level() > 0, "cannot swap algebraic variables" );
197    if ( f.inCoeffDomain() || x1 == x2 || ( x1 > f.mvar() && x2 > f.mvar() ) )
198        return f;
199    else
200    {
201        CanonicalForm result = 0;
202        if ( x1 > x2 ) {
203            sv_x1 = x2; sv_x2 = x1;
204        }
205        else
206        {
207            sv_x1 = x1; sv_x2 = x2;
208        }
209        if ( f.mvar() < sv_x2 )
210            // we only have to replace sv_x1 by sv_x2
211            return swapvar_between1( f );
212        else
213            // we really have to swap variables
214            return swapvar_rec1( f );
215    }
216}
217//}}}
218
219//{{{ static CanonicalForm replacevar_between ( const CanonicalForm & f )
220//{{{ docu
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//}}}
232static CanonicalForm
233replacevar_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
263//{{{ CanonicalForm replacevar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
264//{{{ docu
265//
266// replacevar() - replace all occurences of x1 in f by x2.
267//
268// In contrast to swapvar(), x1 may be an algebraic variable, but
269// x2 must be a polynomial variable.
270//
271//}}}
272CanonicalForm
273replacevar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
274{
275    //ASSERT( x2.level() > 0, "cannot replace with algebraic variable" );
276    if ( f.inBaseDomain() || x1 == x2 || ( x1 > f.mvar() ) )
277        return f;
278    else
279    {
280        sv_x1 = x1;
281        sv_x2 = x2;
282        return replacevar_between( f );
283    }
284}
285//}}}
286
287//{{{ static void fillVarsRec ( const CanonicalForm & f, int * vars )
288//{{{ docu
289//
290// fillVarsRec - fill array describing occurences of variables in f.
291//
292// Only polynomial variables are looked up.  The information is
293// stored in the arrary vars.  vars should be large enough to
294// hold all information, i.e. larger than the level of f.
295//
296// Used by getVars() and getNumVars().
297//
298//}}}
299static void
300fillVarsRec ( const CanonicalForm & f, int * vars )
301{
302    int n;
303    if ( (n = f.level()) > 0 )
304    {
305        vars[n] = 1;
306        CFIterator i;
307        for ( i = f; i.hasTerms(); ++i )
308            fillVarsRec( i.coeff(), vars );
309    }
310}
311//}}}
312
313//{{{ int getNumVars ( const CanonicalForm & f )
314//{{{ docu
315//
316// getNumVars() - get number of polynomial variables in f.
317//
318//}}}
319int
320getNumVars ( const CanonicalForm & f )
321{
322    int n;
323    if ( f.inCoeffDomain() )
324        return 0;
325    else  if ( (n = f.level()) == 1 )
326        return 1;
327    else
328    {
329        int * vars = new int[ n+1 ];
330        int i;
331        for ( i = 0; i < n; i++ ) vars[i] = 0;
332
333        // look for variables
334        for ( CFIterator I = f; I.hasTerms(); ++I )
335            fillVarsRec( I.coeff(), vars );
336
337        // count them
338        int m = 0;
339        for ( i = 1; i < n; i++ )
340            if ( vars[i] != 0 ) m++;
341
342        delete [] vars;
343        // do not forget to count our own variable
344        return m+1;
345    }
346}
347//}}}
348
349//{{{ CanonicalForm getVars ( const CanonicalForm & f )
350//{{{ docu
351//
352// getVars() - get polynomial variables of f.
353//
354// Return the product of all of them, 1 if there are not any.
355//
356//}}}
357CanonicalForm
358getVars ( const CanonicalForm & f )
359{
360    int n;
361    if ( f.inCoeffDomain() )
362        return 1;
363    else  if ( (n = f.level()) == 1 )
364        return Variable( 1 );
365    else
366    {
367        int * vars = new int[ n+1 ];
368        int i;
369        for ( i = 0; i <= n; i++ ) vars[i] = 0;
370
371        // look for variables
372        for ( CFIterator I = f; I.hasTerms(); ++I )
373            fillVarsRec( I.coeff(), vars );
374
375        // multiply them all
376        CanonicalForm result = 1;
377        for ( i = n; i > 0; i-- )
378            if ( vars[i] != 0 ) result *= Variable( i );
379
380        delete [] vars;
381        // do not forget our own variable
382        return f.mvar() * result;
383    }
384}
385//}}}
386
387//{{{ CanonicalForm apply ( const CanonicalForm & f, void (*mf)( CanonicalForm &, int & ) )
388//{{{ docu
389//
390// apply() - apply mf to terms of f.
391//
392// Calls mf( f[i], i ) for each term f[i]*x^i of f and builds a
393// new term from the result.  If f is in a coefficient domain,
394// mf( f, i ) should result in an i == 0, since otherwise it is
395// not clear which variable to use for the resulting term.
396//
397// An example:
398//
399// void
400// diff( CanonicalForm & f, int & i )
401// {
402//     f = f * i;
403//     if ( i > 0 ) i--;
404// }
405//
406// Then apply( f, diff ) is differentation of f with respect to the
407// main variable of f.
408//
409//}}}
410CanonicalForm
411apply ( const CanonicalForm & f, void (*mf)( CanonicalForm &, int & ) )
412{
413    if ( f.inCoeffDomain() )
414    {
415        int exp = 0;
416        CanonicalForm result = f;
417        mf( result, exp );
418        ASSERT( exp == 0, "illegal result, do not know what variable to use" );
419        return result;
420    }
421    else
422    {
423        CanonicalForm result, coeff;
424        CFIterator i;
425        int exp;
426        Variable x = f.mvar();
427        for ( i = f; i.hasTerms(); i++ )
428        {
429            coeff = i.coeff();
430            exp = i.exp();
431            mf( coeff, exp );
432            if ( ! coeff.isZero() )
433                result += power( x, exp ) * coeff;
434        }
435        return result;
436    }
437}
438//}}}
439
440//{{{ CanonicalForm mapdomain ( const CanonicalForm & f, CanonicalForm (*mf)( const CanonicalForm & ) )
441//{{{ docu
442//
443// mapdomain() - map all coefficients of f through mf.
444//
445// Recursively descends down through f to the coefficients which
446// are in a coefficient domain mapping each such coefficient
447// through mf and returns the result.
448//
449//}}}
450CanonicalForm
451mapdomain ( const CanonicalForm & f, CanonicalForm (*mf)( const CanonicalForm & ) )
452{
453    if ( f.inBaseDomain() )
454        return mf( f );
455    else
456    {
457        CanonicalForm result = 0;
458        CFIterator i;
459        Variable x = f.mvar();
460        for ( i = f; i.hasTerms(); i++ )
461            result += power( x, i.exp() ) * mapdomain( i.coeff(), mf );
462        return result;
463    }
464}
465//}}}
466
467//{{{ static void degreesRec ( const CanonicalForm & f, int * degs )
468//{{{ docu
469//
470// degreesRec() - recursively get degrees of f.
471//
472// Used by degrees().
473//
474//}}}
475static void
476degreesRec ( const CanonicalForm & f, int * degs )
477{
478    if ( ! f.inCoeffDomain() )
479    {
480        int level = f.level();
481        int deg = f.degree();
482        // calculate the maximum degree of all coefficients which
483        // are in the same level
484        if ( degs[level] < deg )
485            degs[level] = f.degree();
486        for ( CFIterator i = f; i.hasTerms(); i++ )
487            degreesRec( i.coeff(), degs );
488    }
489}
490//}}}
491
492//{{{ int * degrees ( const CanonicalForm & f, int * degs )
493//{{{ docu
494//
495// degress() - return the degrees of all polynomial variables in f.
496//
497// Returns 0 if f is in a coefficient domain, the degrees of f in
498// all its polynomial variables in an array of int otherwise:
499//
500//   degrees( f, 0 )[i] = degree( f, Variable(i) )
501//
502// If degs is not the zero pointer the degrees are stored in this
503// array.  In this case degs should be larger than the level of
504// f.  If degs is the zero pointer, an array of sufficient size
505// is allocated automatically.
506//
507//}}}
508int * degrees ( const CanonicalForm & f, int * degs )
509{
510    if ( f.inCoeffDomain() )
511        return 0;
512    else
513    {
514        int level = f.level();
515        if ( degs == 0 )
516            degs = new int[level+1];
517        for ( int i = 0; i <= level; i++ )
518            degs[i] = 0;
519        degreesRec( f, degs );
520        return degs;
521    }
522}
523//}}}
524
525//{{{ int totaldegree ( const CanonicalForm & f )
526//{{{ docu
527//
528// totaldegree() - return the total degree of f.
529//
530// If f is zero, return -1.  If f is in a coefficient domain,
531// return 0.  Otherwise return the total degree of f in all
532// polynomial variables.
533//
534//}}}
535int totaldegree ( const CanonicalForm & f )
536{
537    if ( f.isZero() )
538        return -1;
539    else if ( f.inCoeffDomain() )
540        return 0;
541    else
542    {
543        CFIterator i;
544        int cdeg = 0, dummy;
545        // calculate maximum over all coefficients of f, taking
546        // in account our own exponent
547        for ( i = f; i.hasTerms(); i++ )
548            if ( (dummy = totaldegree( i.coeff() ) + i.exp()) > cdeg )
549                cdeg = dummy;
550        return cdeg;
551    }
552}
553//}}}
554
555//{{{ int totaldegree ( const CanonicalForm & f, const Variable & v1, const Variable & v2 )
556//{{{ docu
557//
558// totaldegree() - return the total degree of f as a polynomial
559//   in the polynomial variables between v1 and v2 (inclusively).
560//
561// If f is zero, return -1.  If f is in a coefficient domain,
562// return 0.  Also, return 0 if v1 > v2.  Otherwise, take f to be
563// a polynomial in the polynomial variables between v1 and v2 and
564// return its total degree.
565//
566//}}}
567int
568totaldegree ( const CanonicalForm & f, const Variable & v1, const Variable & v2 )
569{
570    if ( f.isZero() )
571        return -1;
572    else if ( v1 > v2 )
573        return 0;
574    else if ( f.inCoeffDomain() )
575        return 0;
576    else if ( f.mvar() < v1 )
577        return 0;
578    else if ( f.mvar() == v1 )
579        return f.degree();
580    else if ( f.mvar() > v2 )
581    {
582        // v2's level is larger than f's level, descend down
583        CFIterator i;
584        int cdeg = 0, dummy;
585        // calculate maximum over all coefficients of f
586        for ( i = f; i.hasTerms(); i++ )
587            if ( (dummy = totaldegree( i.coeff(), v1, v2 )) > cdeg )
588                cdeg = dummy;
589        return cdeg;
590    }
591    else
592    {
593        // v1 < f.mvar() <= v2
594        CFIterator i;
595        int cdeg = 0, dummy;
596        // calculate maximum over all coefficients of f, taking
597        // in account our own exponent
598        for ( i = f; i.hasTerms(); i++ )
599            if ( (dummy = totaldegree( i.coeff(), v1, v2 ) + i.exp()) > cdeg )
600                cdeg = dummy;
601        return cdeg;
602    }
603}
604//}}}
605
606//{{{ int size ( const CanonicalForm & f, const Variable & v )
607//{{{ docu
608//
609// size() - count number of monomials of f with level higher
610//   or equal than level of v.
611//
612// Returns one if f is in an base domain.
613//
614//}}}
615int
616size ( const CanonicalForm & f, const Variable & v )
617{
618    if ( f.inBaseDomain() )
619        return 1;
620
621    if ( f.mvar() < v )
622        // polynomials with level < v1 are counted as coefficients
623        return 1;
624    else
625    {
626        CFIterator i;
627        int result = 0;
628        // polynomials with level > v2 are not counted al all
629        for ( i = f; i.hasTerms(); i++ )
630            result += size( i.coeff(), v );
631        return result;
632    }
633}
634//}}}
635
636//{{{ int size ( const CanonicalForm & f )
637//{{{ docu
638//
639// size() - return number of monomials in f which are in an
640//   coefficient domain.
641//
642// Returns one if f is in an coefficient domain.
643//
644//}}}
645int
646size ( const CanonicalForm & f )
647{
648    if ( f.inCoeffDomain() )
649        return 1;
650    else
651    {
652        int result = 0;
653        CFIterator i;
654        for ( i = f; i.hasTerms(); i++ )
655            result += size( i.coeff() );
656        return result;
657    }
658}
659//}}}
660
661CanonicalForm reduce(const CanonicalForm & f, const CanonicalForm & M)
662{ // polynomials in M.mvar() are considered coefficients
663  // M univariate monic polynomial
664  // the coefficients of f are reduced modulo M
665  if(f.inBaseDomain() || f.level() < M.level())
666    return f;
667  if(f.level() == M.level())
668  {
669    if(f.degree() < M.degree())
670      return f;
671    CanonicalForm tmp = mod (f, M);
672    return tmp;
673  }
674  // here: f.level() > M.level()
675  CanonicalForm result = 0;
676  for(CFIterator i=f; i.hasTerms(); i++)
677    result += reduce(i.coeff(),M) * power(f.mvar(),i.exp());
678  return result;
679}
Note: See TracBrowser for help on using the repository browser.