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