[75744d] | 1 | /* emacs edit mode for this file is -*- C++ -*- */ |
---|
[341696] | 2 | /* $Id$ */ |
---|
[75744d] | 3 | |
---|
| 4 | //{{{ docu |
---|
| 5 | // |
---|
| 6 | // cf_resultant.cc - algorithms for calculating resultants. |
---|
| 7 | // |
---|
| 8 | // Header file: cf_algorithm.h |
---|
| 9 | // |
---|
| 10 | //}}} |
---|
| 11 | |
---|
[e4fe2b] | 12 | #include "config.h" |
---|
[75744d] | 13 | |
---|
[650f2d8] | 14 | #include "cf_assert.h" |
---|
[75744d] | 15 | |
---|
| 16 | #include "canonicalform.h" |
---|
| 17 | #include "variable.h" |
---|
[72dd6e] | 18 | #include "cf_algorithm.h" |
---|
[75744d] | 19 | |
---|
[54c17f] | 20 | //{{{ CFArray subResChain ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) |
---|
[75744d] | 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 | //}}} |
---|
| 39 | CFArray |
---|
[54c17f] | 40 | subResChain ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) |
---|
[75744d] | 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() ) { |
---|
[806c18] | 50 | trivialResult[0] = 0; |
---|
| 51 | return trivialResult; |
---|
[75744d] | 52 | } |
---|
| 53 | |
---|
| 54 | // make x main variable |
---|
| 55 | if ( f.mvar() > x || g.mvar() > x ) { |
---|
[806c18] | 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 ); |
---|
[75744d] | 62 | } |
---|
| 63 | else { |
---|
[806c18] | 64 | X = x; |
---|
| 65 | F = f; |
---|
| 66 | G = g; |
---|
[75744d] | 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 ) { |
---|
[806c18] | 85 | S[j-1] = LC( S[j], X ) * psr( S[j+1], S[j], X ); |
---|
| 86 | j--; |
---|
[75744d] | 87 | } else if ( m < n ) { |
---|
[806c18] | 88 | S[j-1] = LC( S[j], X ) * LC( S[j], X ) * S[j+1]; |
---|
| 89 | j--; |
---|
[75744d] | 90 | } else if ( m > n && j > 0 ) { |
---|
[806c18] | 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 | } |
---|
[75744d] | 104 | } |
---|
| 105 | |
---|
| 106 | while ( j > 0 ) { |
---|
[806c18] | 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 |
---|
[75744d] | 121 | } |
---|
| 122 | |
---|
| 123 | for ( j = 0; j <= S.max(); j++ ) { |
---|
[806c18] | 124 | // reswap variables if necessary |
---|
| 125 | if ( X != x ) { |
---|
| 126 | S[j] = swapvar( S[j], X, x ); |
---|
| 127 | } |
---|
[75744d] | 128 | } |
---|
| 129 | |
---|
| 130 | return S; |
---|
| 131 | } |
---|
| 132 | //}}} |
---|
[54c17f] | 133 | |
---|
[7e4442] | 134 | //{{{ static CanonicalForm trivialResultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) |
---|
[54c17f] | 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 | // |
---|
[f0596e] | 142 | // Used by resultant(). |
---|
| 143 | // |
---|
[54c17f] | 144 | //}}} |
---|
| 145 | static CanonicalForm |
---|
[7e4442] | 146 | trivialResultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) |
---|
[54c17f] | 147 | { |
---|
| 148 | // f or g in R |
---|
| 149 | if ( degree( f, x ) == 0 ) |
---|
[806c18] | 150 | return power( f, degree( g, x ) ); |
---|
[54c17f] | 151 | if ( degree( g, x ) == 0 ) |
---|
[806c18] | 152 | return power( g, degree( f, x ) ); |
---|
[54c17f] | 153 | |
---|
| 154 | // f and g are linear polynomials |
---|
| 155 | return LC( f, x ) * g - LC( g, x ) * f; |
---|
| 156 | } |
---|
| 157 | //}}} |
---|
| 158 | |
---|
[7e4442] | 159 | //{{{ CanonicalForm resultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) |
---|
[54c17f] | 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 | //}}} |
---|
| 174 | CanonicalForm |
---|
[7e4442] | 175 | resultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x ) |
---|
[54c17f] | 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() ) |
---|
[806c18] | 182 | return 0; |
---|
[54c17f] | 183 | if ( f.mvar() < x ) |
---|
[806c18] | 184 | return power( f, g.degree( x ) ); |
---|
[54c17f] | 185 | if ( g.mvar() < x ) |
---|
[806c18] | 186 | return power( g, f.degree( x ) ); |
---|
[54c17f] | 187 | |
---|
| 188 | // make x main variale |
---|
| 189 | CanonicalForm F, G; |
---|
| 190 | Variable X; |
---|
| 191 | if ( f.mvar() > x || g.mvar() > x ) { |
---|
[806c18] | 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 ); |
---|
[54c17f] | 198 | } |
---|
| 199 | else { |
---|
[806c18] | 200 | X = x; |
---|
| 201 | F = f; |
---|
| 202 | G = g; |
---|
[54c17f] | 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 ) |
---|
[806c18] | 212 | return swapvar( trivialResultant( F, G, X ), X, x ); |
---|
[54c17f] | 213 | |
---|
| 214 | // exchange F and G if necessary |
---|
| 215 | int flipFactor; |
---|
| 216 | if ( m < n ) { |
---|
[806c18] | 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; |
---|
[54c17f] | 225 | } else |
---|
[806c18] | 226 | flipFactor = 1; |
---|
[54c17f] | 227 | |
---|
| 228 | // this is not an effective way to calculate the resultant! |
---|
| 229 | CanonicalForm extFactor; |
---|
| 230 | if ( m == n ) { |
---|
[806c18] | 231 | if ( n & 1 ) |
---|
| 232 | extFactor = -LC( G, X ); |
---|
| 233 | else |
---|
| 234 | extFactor = LC( G, X ); |
---|
[54c17f] | 235 | } else |
---|
[806c18] | 236 | extFactor = power( LC( F, X ), m-n-1 ); |
---|
[54c17f] | 237 | |
---|
| 238 | CanonicalForm result; |
---|
| 239 | result = subResChain( F, G, X )[0] / extFactor; |
---|
| 240 | |
---|
| 241 | return swapvar( result, X, x ) * flipFactor; |
---|
| 242 | } |
---|
| 243 | //}}} |
---|