001/* 002 * $Id: GreatestCommonDivisorSimple.java 4025 2012-07-23 16:41:43Z kredel $ 003 */ 004 005package edu.jas.ufd; 006 007 008import org.apache.log4j.Logger; 009 010import edu.jas.poly.GenPolynomial; 011import edu.jas.poly.PolyUtil; 012import edu.jas.structure.GcdRingElem; 013import edu.jas.structure.RingFactory; 014import edu.jas.structure.Power; 015 016 017/** 018 * Greatest common divisor algorithms with monic polynomial remainder sequence. 019 * If C is a field, then the monic PRS (on coefficients) is computed otherwise 020 * no simplifications in the reduction are made. 021 * @author Heinz Kredel 022 */ 023 024public class GreatestCommonDivisorSimple<C extends GcdRingElem<C>> extends GreatestCommonDivisorAbstract<C> { 025 026 027 private static final Logger logger = Logger.getLogger(GreatestCommonDivisorSimple.class); 028 029 030 private final boolean debug = logger.isDebugEnabled(); 031 032 033 /** 034 * Univariate GenPolynomial greatest comon divisor. Uses pseudoRemainder for 035 * remainder. 036 * @param P univariate GenPolynomial. 037 * @param S univariate GenPolynomial. 038 * @return gcd(P,S). 039 */ 040 @Override 041 public GenPolynomial<C> baseGcd(GenPolynomial<C> P, GenPolynomial<C> S) { 042 if (S == null || S.isZERO()) { 043 return P; 044 } 045 if (P == null || P.isZERO()) { 046 return S; 047 } 048 if (P.ring.nvar > 1) { 049 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial"); 050 } 051 boolean field = P.ring.coFac.isField(); 052 long e = P.degree(0); 053 long f = S.degree(0); 054 GenPolynomial<C> q; 055 GenPolynomial<C> r; 056 if (f > e) { 057 r = P; 058 q = S; 059 long g = f; 060 f = e; 061 e = g; 062 } else { 063 q = P; 064 r = S; 065 } 066 if (debug) { 067 logger.debug("degrees: e = " + e + ", f = " + f); 068 } 069 C c; 070 if (field) { 071 r = r.monic(); 072 q = q.monic(); 073 c = P.ring.getONECoefficient(); 074 } else { 075 r = r.abs(); 076 q = q.abs(); 077 C a = baseContent(r); 078 C b = baseContent(q); 079 c = gcd(a, b); // indirection 080 r = divide(r, a); // indirection 081 q = divide(q, b); // indirection 082 } 083 if (r.isONE()) { 084 return r.multiply(c); 085 } 086 if (q.isONE()) { 087 return q.multiply(c); 088 } 089 GenPolynomial<C> x; 090 //System.out.println("q = " + q); 091 //System.out.println("r = " + r); 092 while (!r.isZERO()) { 093 x = PolyUtil.<C> baseSparsePseudoRemainder(q, r); 094 q = r; 095 if (field) { 096 r = x.monic(); 097 } else { 098 r = x; 099 } 100 //System.out.println("q = " + q); 101 //System.out.println("r = " + r); 102 } 103 q = basePrimitivePart(q); 104 return (q.multiply(c)).abs(); 105 } 106 107 108 /** 109 * Univariate GenPolynomial recursive greatest comon divisor. Uses 110 * pseudoRemainder for remainder. 111 * @param P univariate recursive GenPolynomial. 112 * @param S univariate recursive GenPolynomial. 113 * @return gcd(P,S). 114 */ 115 @Override 116 public GenPolynomial<GenPolynomial<C>> recursiveUnivariateGcd(GenPolynomial<GenPolynomial<C>> P, 117 GenPolynomial<GenPolynomial<C>> S) { 118 if (S == null || S.isZERO()) { 119 return P; 120 } 121 if (P == null || P.isZERO()) { 122 return S; 123 } 124 if (P.ring.nvar > 1) { 125 throw new IllegalArgumentException(this.getClass().getName() + " no univariate polynomial"); 126 } 127 boolean field = P.leadingBaseCoefficient().ring.coFac.isField(); 128 long e = P.degree(0); 129 long f = S.degree(0); 130 GenPolynomial<GenPolynomial<C>> q; 131 GenPolynomial<GenPolynomial<C>> r; 132 if (f > e) { 133 r = P; 134 q = S; 135 long g = f; 136 f = e; 137 e = g; 138 } else { 139 q = P; 140 r = S; 141 } 142 if (debug) { 143 logger.debug("degrees: e = " + e + ", f = " + f); 144 } 145 if (field) { 146 r = PolyUtil.<C> monic(r); 147 q = PolyUtil.<C> monic(q); 148 } else { 149 r = r.abs(); 150 q = q.abs(); 151 } 152 GenPolynomial<C> a = recursiveContent(r); 153 GenPolynomial<C> b = recursiveContent(q); 154 155 GenPolynomial<C> c = gcd(a, b); // go to recursion 156 //System.out.println("rgcd c = " + c); 157 r = PolyUtil.<C> recursiveDivide(r, a); 158 q = PolyUtil.<C> recursiveDivide(q, b); 159 if (r.isONE()) { 160 return r.multiply(c); 161 } 162 if (q.isONE()) { 163 return q.multiply(c); 164 } 165 GenPolynomial<GenPolynomial<C>> x; 166 while (!r.isZERO()) { 167 x = PolyUtil.<C> recursivePseudoRemainder(q, r); 168 q = r; 169 if (field) { 170 r = PolyUtil.<C> monic(x); 171 } else { 172 r = x; 173 } 174 } 175 q = recursivePrimitivePart(q); 176 q = q.abs().multiply(c); 177 return q; 178 } 179 180 181 /** 182 * Univariate GenPolynomial resultant. 183 * @param P univariate GenPolynomial. 184 * @param S univariate GenPolynomial. 185 * @return res(P,S). 186 */ 187 @Override 188 public GenPolynomial<C> baseResultant(GenPolynomial<C> P, GenPolynomial<C> S) { 189 if (S == null || S.isZERO()) { 190 return S; 191 } 192 if (P == null || P.isZERO()) { 193 return P; 194 } 195 if (P.ring.nvar > 1 || P.ring.nvar == 0) { 196 throw new IllegalArgumentException("no univariate polynomial"); 197 } 198 long e = P.degree(0); 199 long f = S.degree(0); 200 if (f == 0 && e == 0) { 201 return P.ring.getONE(); 202 } 203 if ( e == 0 ) { 204 return Power.<GenPolynomial<C>> power(P.ring,P,f); 205 } 206 if ( f == 0 ) { 207 return Power.<GenPolynomial<C>> power(S.ring,S,e); 208 } 209 GenPolynomial<C> q; 210 GenPolynomial<C> r; 211 int s = 0; // sign is +, 1 for sign is - 212 if (e < f) { 213 r = P; 214 q = S; 215 long t = e; 216 e = f; 217 f = t; 218 if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f) 219 s = 1; 220 } 221 } else { 222 q = P; 223 r = S; 224 } 225 RingFactory<C> cofac = P.ring.coFac; 226 boolean field = cofac.isField(); 227 C c = cofac.getONE(); 228 GenPolynomial<C> x; 229 long g; 230 do { 231 if (field) { 232 x = q.remainder(r); 233 } else { 234 x = PolyUtil.<C>baseSparsePseudoRemainder(q,r); 235 //System.out.println("x_s = " + x + ", lbcf(r) = " + r.leadingBaseCoefficient()); 236 } 237 if ( x.isZERO() ) { 238 return x; 239 } 240 //System.out.println("x = " + x); 241 e = q.degree(0); 242 f = r.degree(0); 243 if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f) 244 s = 1 - s; 245 } 246 g = x.degree(0); 247 C c2 = r.leadingBaseCoefficient(); 248 for (int i = 0; i < (e-g); i++ ) { 249 c = c.multiply(c2); 250 } 251 q = r; 252 r = x; 253 } while (g != 0); 254 C c2 = r.leadingBaseCoefficient(); 255 for (int i = 0; i < f; i++ ) { 256 c = c.multiply(c2); 257 } 258 if ( s == 1 ) { 259 c = c.negate(); 260 } 261 x = P.ring.getONE().multiply(c); 262 return x; 263 } 264 265 266 /** 267 * Univariate GenPolynomial recursive resultant. 268 * @param P univariate recursive GenPolynomial. 269 * @param S univariate recursive GenPolynomial. 270 * @return res(P,S). 271 */ 272 @Override 273 public GenPolynomial<GenPolynomial<C>> recursiveUnivariateResultant(GenPolynomial<GenPolynomial<C>> P, 274 GenPolynomial<GenPolynomial<C>> S) { 275 if (S == null || S.isZERO()) { 276 return S; 277 } 278 if (P == null || P.isZERO()) { 279 return P; 280 } 281 if (P.ring.nvar > 1 || P.ring.nvar == 0) { 282 throw new IllegalArgumentException("no recursive univariate polynomial"); 283 } 284 long e = P.degree(0); 285 long f = S.degree(0); 286 if ( f == 0 && e == 0 ) { 287 // if coeffs are multivariate (and non constant) 288 // otherwise it would be 1 289 GenPolynomial<C> t = resultant(P.leadingBaseCoefficient(), S.leadingBaseCoefficient()); 290 return P.ring.getONE().multiply(t); 291 } 292 if ( e == 0 ) { 293 return Power.<GenPolynomial<GenPolynomial<C>>> power(P.ring,P,f); 294 } 295 if ( f == 0 ) { 296 return Power.<GenPolynomial<GenPolynomial<C>>> power(S.ring,S,e); 297 } 298 GenPolynomial<GenPolynomial<C>> q; 299 GenPolynomial<GenPolynomial<C>> r; 300 int s = 0; // sign is +, 1 for sign is - 301 if (f > e) { 302 r = P; 303 q = S; 304 long g = f; 305 f = e; 306 e = g; 307 if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f) 308 s = 1; 309 } 310 } else { 311 q = P; 312 r = S; 313 } 314 GenPolynomial<GenPolynomial<C>> x; 315 RingFactory<GenPolynomial<C>> cofac = P.ring.coFac; 316 GenPolynomial<C> c = cofac.getONE(); 317 long g; 318 do { 319 x = PolyUtil.<C>recursiveSparsePseudoRemainder(q,r); 320 //x = PolyUtil.<C>recursiveDensePseudoRemainder(q,r); 321 if ( x.isZERO() ) { 322 return x; 323 } 324 //no: x = recursivePrimitivePart(x); 325 //System.out.println("x = " + x); 326 e = q.degree(0); 327 f = r.degree(0); 328 if ((e % 2 != 0) && (f % 2 != 0)) { // odd(e) && odd(f) 329 s = 1 - s; 330 } 331 g = x.degree(0); 332 GenPolynomial<C> c2 = r.leadingBaseCoefficient(); 333 for (int i = 0; i < (e-g); i++ ) { 334 c = c.multiply(c2); 335 } 336 q = r; 337 r = x; 338 } while (g != 0); 339 GenPolynomial<C> c2 = r.leadingBaseCoefficient(); 340 for (int i = 0; i < f; i++ ) { 341 c = c.multiply(c2); 342 } 343 if ( s == 1 ) { 344 c = c.negate(); 345 } 346 x = P.ring.getONE().multiply(c); 347 return x; 348 } 349 350}