source: git/factory/cf_resultant.cc @ 650f2d8

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