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