source: git/factory/cf_map.cc @ 9752db

spielwiese
Last change on this file since 9752db was e88604, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: deleted sparse modular gcd from sm_sparsemod.*
  • Property mode set to 100644
File size: 10.1 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4//{{{ docu
5//
6// cf_map.cc - definition of class CFMap.
7//
8// Used by: cf_gcd.cc, fac_multivar.cc
9//
10//}}}
11
12#include "config.h"
13
14#include "canonicalform.h"
15#include "cf_map.h"
16#include "cf_iter.h"
17#include "templates/ftmpl_functions.h"
18
19//{{{ MapPair & MapPair::operator = ( const MapPair & p )
20//{{{ docu
21//
22// MapPair::operator = - assignment operator.
23//
24//}}}
25MapPair &
26MapPair::operator = ( const MapPair & p )
27{
28    if ( this != &p ) {
29        V = p.V;
30        S = p.S;
31    }
32    return *this;
33}
34//}}}
35
36#ifndef NOSTREAMIO
37//{{{ OSTREAM & operator << ( OSTREAM & s, const MapPair & p )
38//{{{ docu
39//
40// operator << - print a map pair ("V -> S").
41//
42//}}}
43OSTREAM &
44operator << ( OSTREAM & s, const MapPair & p )
45{
46    s << p.var() << " -> " << p.subst();
47    return s;
48}
49//}}}
50
51void MapPair::print( OSTREAM&) const
52{
53}
54#endif /* NOSTREAMIO */
55
56//{{{ CFMap::CFMap ( const CFList & L )
57//{{{ docu
58//
59// CFMap::CFMap() - construct a CFMap from a CFList.
60//
61// Variable[i] will be mapped to CFList[i] under the resulting
62// map.
63//
64//}}}
65CFMap::CFMap ( const CFList & L )
66{
67    CFListIterator i;
68    int j;
69    for ( i = L, j = 1; i.hasItem(); i++, j++ )
70        P.insert( MapPair( Variable(j), i.getItem() ) );
71}
72//}}}
73
74//{{{ CFMap & CFMap::operator = ( const CFMap & m )
75//{{{ docu
76//
77// CFMap::operator = - assignment operator.
78//
79//}}}
80CFMap &
81CFMap::operator = ( const CFMap & m )
82{
83    if ( this != &m )
84        P = m.P;
85    return *this;
86}
87//}}}
88
89//{{{ static int cmpfunc ( const MapPair & p1, const MapPair & p2 )
90//{{{ docu
91//
92// cmpfunc() - compare two map pairs.
93//
94// Return -1 if p2's variable is less than p1's, 0 if they are
95// equal, 1 if p2's level is greater than p1's.
96//
97//}}}
98static int
99cmpfunc ( const MapPair & p1, const MapPair & p2 )
100{
101    if ( p1.var() > p2.var() ) return -1;
102    else if ( p1.var() == p2.var() ) return 0;
103    else return 1;
104}
105//}}}
106
107//{{{ static void insfunc ( MapPair & orgp, const MapPair & newp )
108//{{{ docu
109//
110// insfunc() - assign newp to orgp.
111//
112// cmpfunc() and insfunc() are used as functions for inserting a
113// map pair into a map by CFMap::newpair().
114//
115//}}}
116static void
117insfunc ( MapPair & orgp, const MapPair & newp )
118{
119    orgp = newp;
120}
121//}}}
122
123//{{{ void CFMap::newpair ( const Variable & v, const CanonicalForm & s )
124//{{{ docu
125//
126// CFMap::newpair() - insert a MapPair into a CFMap.
127//
128//}}}
129void
130CFMap::newpair ( const Variable & v, const CanonicalForm & s )
131{
132    P.insert( MapPair( v, s ), cmpfunc, insfunc );
133}
134//}}}
135
136//{{{ static CanonicalForm subsrec ( const CanonicalForm & f, const MPListIterator & i )
137//{{{ docu
138//
139// subsrec() - recursively apply the substitutions in i to f.
140//
141// Substitutes algebraic variables, too.  The substituted
142// expression are not subject to further substitutions.
143//
144// Used by: CFMap::operator ()().
145//
146//}}}
147static CanonicalForm
148subsrec ( const CanonicalForm & f, const MPListIterator & i )
149{
150    if ( f.inBaseDomain() ) return f;
151    MPListIterator j = i;
152
153    // skip MapPairs larger than the main variable of f
154    while ( j.hasItem() && j.getItem().var() > f.mvar() ) j++;
155
156    if ( j.hasItem() )
157        if ( j.getItem().var() != f.mvar() ) {
158            // simply descend if the current MapPair variable is
159            // not the main variable of f
160            CanonicalForm result = 0;
161            CFIterator I;
162            for ( I = f; I.hasTerms(); I++ )
163                result += power( f.mvar(), I.exp() ) * subsrec( I.coeff(), j );
164            return result;
165        }
166        else {
167            // replace the main variable of f with the image of
168            // the current variable under MapPair
169            CanonicalForm result = 0;
170            CanonicalForm s = j.getItem().subst();
171            CFIterator I;
172            // move on to the next MapPair
173            j++;
174            for ( I = f; I.hasTerms(); I++ )
175                result += subsrec( I.coeff(), j ) * power( s, I.exp() );
176            return result;
177        }
178    else
179        return f;
180}
181//}}}
182
183//{{{ CanonicalForm CFMap::operator () ( const CanonicalForm & f ) const
184//{{{ docu
185//
186// CFMap::operator () - apply CO to f.
187//
188// See subsrec() for more detailed information.
189//
190//}}}
191CanonicalForm
192CFMap::operator () ( const CanonicalForm & f ) const
193{
194    MPListIterator i = P;
195    return subsrec( f, i );
196}
197//}}}
198
199#ifndef NOSTREAMIO
200//{{{ OSTREAM & operator << ( OSTREAM & s, const CFMap & m )
201//{{{ docu
202//
203// operator << - print a CFMap ("( V[1] -> S[1], ..., V[n] -> // S[n] )".
204//
205//}}}
206OSTREAM &
207operator << ( OSTREAM & s, const CFMap & m )
208{
209    m.P.print(s);
210    return s;
211}
212//}}}
213#endif /* NOSTREAMIO */
214
215//{{{ CanonicalForm compress ( const CanonicalForm & f, CFMap & m )
216//{{{ docu
217//
218// compress() - compress the canonical form f.
219//
220// Compress the polynomial f such that the levels of its
221// polynomial variables are ordered without any gaps starting
222// from level 1.  Return the compressed polynomial and a map m to
223// undo the compression.  That is, if f' = compress(f, m), than f
224// = m(f').
225//
226//}}}
227CanonicalForm
228compress ( const CanonicalForm & f, CFMap & m )
229{
230    CanonicalForm result = f;
231    int i, n;
232    int * degs = degrees( f );
233
234    m = CFMap();
235    n = i = 1;
236    while ( i <= level( f ) ) {
237        while( degs[i] == 0 ) i++;
238        if ( i != n ) {
239            // swap variables and remember the swap in the map
240            m.newpair( Variable( n ), Variable( i ) );
241            result = swapvar( result, Variable( i ), Variable( n ) );
242        }
243        n++; i++;
244    }
245    delete [] degs;
246    return result;
247}
248//}}}
249
250//{{{ void compress ( const CFArray & a, CFMap & M, CFMap & N )
251//{{{ docu
252//
253// compress() - compress the variables occuring in an a.
254//
255// Compress the polynomial variables occuring in a so that their
256// levels are ordered without any gaps starting from level 1.
257// Return the CFMap M to realize the compression and its inverse,
258// the CFMap N.  Note that if you compress a member of a using M
259// the result of the compression is not necessarily compressed,
260// since the map is constructed using all variables occuring in
261// a.
262//
263//}}}
264void
265compress ( const CFArray & a, CFMap & M, CFMap & N )
266{
267    M = N = CFMap();
268    if ( a.size() == 0 )
269        return;
270    int maxlevel = level( a[a.min()] );
271    int i, j;
272
273    // get the maximum of levels in a
274    for ( i = a.min() + 1; i <= a.max(); i++ )
275        if ( level( a[i] ) > maxlevel )
276            maxlevel = level( a[i] );
277    if ( maxlevel <= 0 )
278        return;
279
280    int * degs = new int[maxlevel+1];
281    int * tmp = new int[maxlevel+1];
282    for ( i = 1; i <= maxlevel; i++ )
283        degs[i] = 0;
284
285    // calculate the union of all levels occuring in a
286    for ( i = a.min(); i <= a.max(); i++ ) {
287        tmp = degrees( a[i], tmp );
288        for ( j = 1; j <= level( a[i] ); j++ )
289            if ( tmp[j] != 0 )
290                degs[j] = 1;
291    }
292
293    // create the maps
294    i = 1; j = 1;
295    while ( i <= maxlevel ) {
296        if ( degs[i] != 0 ) {
297            M.newpair( Variable(i), Variable(j) );
298            N.newpair( Variable(j), Variable(i) );
299            j++;
300        }
301        i++;
302    }
303    delete [] tmp;
304    delete [] degs;
305}
306//}}}
307
308/*
309*  compute positions p1 and pe of optimal variables:
310*    pe is used in "ezgcd" and
311*    p1 in "gcd_poly1"
312*/
313static
314void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
315{
316    int i, o1, oe;
317    i = p1 = pe = 0;
318    do
319    {
320        i++;
321        if ( i > n ) return;
322    } while ( ( df[i] == 0 ) || ( dg[i] == 0 ) );
323    p1 = pe = i;
324    if ( df[i] > dg[i] )
325    {
326        o1 = df[i]; oe = dg[i];
327    }
328    else
329    {
330        o1 = dg[i]; oe = df[i];
331    }
332    while ( i < n )
333    {
334        i++;
335        if ( ( df[i] != 0 ) && ( dg[i] != 0 ) )
336        {
337            if ( df[i] > dg[i] )
338            {
339                if ( o1 >= df[i]) { o1 = df[i]; p1 = i; }
340                if ( oe < dg[i]) { oe = dg[i]; pe = i; }
341            }
342            else
343            {
344                if ( o1 >= dg[i]) { o1 = dg[i]; p1 = i; }
345                if ( oe < df[i]) { oe = df[i]; pe = i; }
346            }
347        }
348    }
349}
350
351
352//{{{ void compress ( const CanonicalForm & f, const CanonicalForm & g, CFMap & M, CFMap & N )
353//{{{ docu
354//
355// compress() - compress the variables occurring in f and g with respect
356// to optimal variables
357//
358// Compress the polynomial variables occurring in f and g so that
359// the levels of variables common to f and g are ordered without
360// any gaps starting from level 1, whereas the variables occuring
361// in only one of f or g are moved to levels higher than the
362// levels of the common variables.  Return the CFMap M to realize
363// the compression and its inverse, the CFMap N.
364// N needs only variables common to f and g.
365//
366//}}}
367void
368compress ( const CanonicalForm & f, const CanonicalForm & g, CFMap & M, CFMap & N )
369{
370    int n = tmax( f.level(), g.level() );
371    int i, k, p1, pe;
372    int * degsf = new int[n+1];
373    int * degsg = new int[n+1];
374
375    for ( i = 0; i <= n; i++ )
376    {
377        degsf[i] = degsg[i] = 0;
378    }
379
380    degsf = degrees( f, degsf );
381    degsg = degrees( g, degsg );
382    optvalues( degsf, degsg, n, p1, pe );
383
384    i = 1; k = 1;
385    if ( pe > 1 )
386    {
387        M.newpair( Variable(pe), Variable(k) );
388        N.newpair( Variable(k), Variable(pe) );
389        k++;
390    }
391    while ( i <= n )
392    {
393        if ( degsf[i] > 0 && degsg[i] > 0 )
394        {
395            if ( ( i != k ) && ( i != pe ) && ( i != p1 ) )
396            {
397                M.newpair( Variable(i), Variable(k) );
398                N.newpair( Variable(k), Variable(i) );
399            }
400            k++;
401        }
402        i++;
403    }
404    if ( p1 != pe )
405    {
406        M.newpair( Variable(p1), Variable(k) );
407        N.newpair( Variable(k), Variable(p1) );
408        k++;
409    }
410    i = 1;
411    while ( i <= n )
412    {
413        if ( degsf[i] > 0 && degsg[i] == 0 ) {
414            if ( i != k )
415            {
416                M.newpair( Variable(i), Variable(k) );
417                k++;
418            }
419        }
420        else if ( degsf[i] == 0 && degsg[i] > 0 )
421        {
422            if ( i != k )
423            {
424                M.newpair( Variable(i), Variable(k) );
425                k++;
426            }
427        }
428        i++;
429    }
430
431    delete [] degsf;
432    delete [] degsg;
433}
434//}}}
Note: See TracBrowser for help on using the repository browser.