source: git/factory/cf_resultant.cc

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