source: git/factory/cf_ops.cc @ 24c6177

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