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 | **/ |
---|
41 | CFArray |
---|
42 | subResChain ( 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 | **/ |
---|
145 | static CanonicalForm |
---|
146 | trivialResultant ( 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 | **/ |
---|
172 | CanonicalForm |
---|
173 | resultant ( 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 | } |
---|