source: git/factory/cf_resultant.cc @ 362fc67

jengelh-datetimespielwiese
Last change on this file since 362fc67 was 362fc67, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: remove $Id$
  • Property mode set to 100644
File size: 6.6 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3//{{{ docu
4//
5// cf_resultant.cc - algorithms for calculating resultants.
6//
7// Header file: cf_algorithm.h
8//
9//}}}
10
11#include "config.h"
12
13#include "cf_assert.h"
14
15#include "canonicalform.h"
16#include "variable.h"
17#include "cf_algorithm.h"
18
19//{{{ CFArray subResChain ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
20//{{{ docu
21//
22// subResChain() - caculate extended subresultant chain.
23//
24// The chain is calculated from f and g with respect to variable
25// x which should not be an algebraic variable.  If f or q equals
26// zero, an array consisting of one zero entry is returned.
27//
28// Note: this is not the standard subresultant chain but the
29// *extended* chain!
30//
31// This algorithm is from the article of R. Loos - 'Generalized
32// Polynomial Remainder Sequences' in B. Buchberger - 'Computer
33// Algebra - Symbolic and Algebraic Computation' with some
34// necessary extensions concerning the calculation of the first
35// step.
36//
37//}}}
38CFArray
39subResChain ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
40{
41    ASSERT( x.level() > 0, "cannot calculate subresultant sequence with respect to algebraic variables" );
42
43    CFArray trivialResult( 0, 0 );
44    CanonicalForm F, G;
45    Variable X;
46
47    // some checks on triviality
48    if ( f.isZero() || g.isZero() ) {
49        trivialResult[0] = 0;
50        return trivialResult;
51    }
52
53    // make x main variable
54    if ( f.mvar() > x || g.mvar() > x ) {
55        if ( f.mvar() > g.mvar() )
56            X = f.mvar();
57        else
58            X = g.mvar();
59        F = swapvar( f, X, x );
60        G = swapvar( g, X, x );
61    }
62    else {
63        X = x;
64        F = f;
65        G = g;
66    }
67    // at this point, we have to calculate the sequence of F and
68    // G in respect to X where X is equal to or greater than the
69    // main variables of F and G
70
71    // initialization of chain
72    int m = degree( F, X );
73    int n = degree( G, X );
74
75    int j = (m <= n) ? n : m-1;
76    int r;
77
78    CFArray S( 0, j+1 );
79    CanonicalForm R;
80    S[j+1] = F; S[j] = G;
81
82    // make sure that S[j+1] is regular and j < n
83    if ( m == n && j > 0 ) {
84        S[j-1] = LC( S[j], X ) * psr( S[j+1], S[j], X );
85        j--;
86    } else if ( m < n ) {
87        S[j-1] = LC( S[j], X ) * LC( S[j], X ) * S[j+1];
88        j--;
89    } else if ( m > n && j > 0 ) {
90        // calculate first step
91        r = degree( S[j], X );
92        R = LC( S[j+1], X );
93
94        // if there was a gap calculate similar polynomial
95        if ( j > r && r >= 0 )
96            S[r] = power( LC( S[j], X ), j - r ) * S[j] * power( R, j - r );
97
98        if ( r > 0 ) {
99            // calculate remainder
100            S[r-1] = psr( S[j+1], S[j], X ) * power( -R, j - r );
101            j = r-1;
102        }
103    }
104
105    while ( j > 0 ) {
106        // at this point, 0 < j < n and S[j+1] is regular
107        r = degree( S[j], X );
108        R = LC( S[j+1], X );
109
110        // if there was a gap calculate similar polynomial
111        if ( j > r && r >= 0 )
112            S[r] = (power( LC( S[j], X ), j - r ) * S[j]) / power( R, j - r );
113
114        if ( r <= 0 ) break;
115        // calculate remainder
116        S[r-1] = psr( S[j+1], S[j], X ) / power( -R, j - r + 2 );
117
118        j = r-1;
119        // again 0 <= j < r <= jOld and S[j+1] is regular
120    }
121
122    for ( j = 0; j <= S.max(); j++ ) {
123        // reswap variables if necessary
124        if ( X != x ) {
125            S[j] = swapvar( S[j], X, x );
126        }
127    }
128
129    return S;
130}
131//}}}
132
133//{{{ static CanonicalForm trivialResultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
134//{{{ docu
135//
136// trivialResultant - calculate trivial resultants.
137//
138// x's level should be larger than f's and g's levels.  Either f
139// or g should be constant or both linear.
140//
141// Used by resultant().
142//
143//}}}
144static CanonicalForm
145trivialResultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
146{
147    // f or g in R
148    if ( degree( f, x ) == 0 )
149        return power( f, degree( g, x ) );
150    if ( degree( g, x ) == 0 )
151        return power( g, degree( f, x ) );
152
153    // f and g are linear polynomials
154    return LC( f, x ) * g - LC( g, x ) * f;
155}
156//}}}
157
158//{{{ CanonicalForm resultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
159//{{{ docu
160//
161// resultant() - return resultant of f and g with respect to x.
162//
163// The chain is calculated from f and g with respect to variable
164// x which should not be an algebraic variable.  If f or q equals
165// zero, zero is returned.  If f is a coefficient with respect to
166// x, f^degree(g, x) is returned, analogously for g.
167//
168// This algorithm serves as a wrapper around other resultant
169// algorithms which do the real work.  Here we use standard
170// properties of resultants only.
171//
172//}}}
173CanonicalForm
174resultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
175{
176    ASSERT( x.level() > 0, "cannot calculate resultant with respect to algebraic variables" );
177
178    // some checks on triviality.  We will not use degree( v )
179    // here because this may involve variable swapping.
180    if ( f.isZero() || g.isZero() )
181        return 0;
182    if ( f.mvar() < x )
183        return power( f, g.degree( x ) );
184    if ( g.mvar() < x )
185        return power( g, f.degree( x ) );
186
187    // make x main variale
188    CanonicalForm F, G;
189    Variable X;
190    if ( f.mvar() > x || g.mvar() > x ) {
191        if ( f.mvar() > g.mvar() )
192            X = f.mvar();
193        else
194            X = g.mvar();
195        F = swapvar( f, X, x );
196        G = swapvar( g, X, x );
197    }
198    else {
199        X = x;
200        F = f;
201        G = g;
202    }
203    // at this point, we have to calculate resultant( F, G, X )
204    // where X is equal to or greater than the main variables
205    // of F and G
206
207    int m = degree( F, X );
208    int n = degree( G, X );
209    // catch trivial cases
210    if ( m+n <= 2 || m == 0 || n == 0 )
211        return swapvar( trivialResultant( F, G, X ), X, x );
212
213    // exchange F and G if necessary
214    int flipFactor;
215    if ( m < n ) {
216        CanonicalForm swap = F;
217        F = G; G = swap;
218        int degswap = m;
219        m = n; n = degswap;
220        if ( m & 1 && n & 1 )
221            flipFactor = -1;
222        else
223            flipFactor = 1;
224    } else
225        flipFactor = 1;
226
227    // this is not an effective way to calculate the resultant!
228    CanonicalForm extFactor;
229    if ( m == n ) {
230        if ( n & 1 )
231            extFactor = -LC( G, X );
232        else
233            extFactor = LC( G, X );
234    } else
235        extFactor = power( LC( F, X ), m-n-1 );
236
237    CanonicalForm result;
238    result = subResChain( F, G, X )[0] / extFactor;
239
240    return swapvar( result, X, x ) * flipFactor;
241}
242//}}}
Note: See TracBrowser for help on using the repository browser.