001/* 002 * $Id: Complex.java 4125 2012-08-19 19:05:22Z kredel $ 003 */ 004 005package edu.jas.poly; 006 007 008import org.apache.log4j.Logger; 009 010import edu.jas.arith.BigComplex; 011import edu.jas.arith.BigDecimal; 012import edu.jas.arith.BigInteger; 013import edu.jas.arith.BigRational; 014import edu.jas.structure.RingElem; 015import edu.jas.structure.GcdRingElem; 016import edu.jas.structure.StarRingElem; 017 018 019/** 020 * Generic Complex class implementing the RingElem interface. Objects of this 021 * class are immutable. 022 * @param <C> base type of RingElem (for complex polynomials). 023 * @author Heinz Kredel 024 */ 025public class Complex<C extends RingElem<C>> implements StarRingElem<Complex<C>>, GcdRingElem<Complex<C>> { 026 027 028 private static final Logger logger = Logger.getLogger(Complex.class); 029 030 031 private static final boolean debug = logger.isDebugEnabled(); 032 033 034 /** 035 * Complex class factory data structure. 036 */ 037 public final ComplexRing<C> ring; 038 039 040 /** 041 * Real part of the data structure. 042 */ 043 protected final C re; 044 045 046 /** 047 * Imaginary part of the data structure. 048 */ 049 protected final C im; 050 051 052 /** 053 * The constructor creates a Complex object from two C objects as real and 054 * imaginary part. 055 * @param ring factory for Complex objects. 056 * @param r real part. 057 * @param i imaginary part. 058 */ 059 public Complex(ComplexRing<C> ring, C r, C i) { 060 this.ring = ring; 061 re = r; 062 im = i; 063 } 064 065 066 /** 067 * The constructor creates a Complex object from a C object as real part, 068 * the imaginary part is set to 0. 069 * @param r real part. 070 */ 071 public Complex(ComplexRing<C> ring, C r) { 072 this(ring, r, ring.ring.getZERO()); 073 } 074 075 076 /** 077 * The constructor creates a Complex object from a long element as real 078 * part, the imaginary part is set to 0. 079 * @param r real part. 080 */ 081 public Complex(ComplexRing<C> ring, long r) { 082 this(ring, ring.ring.fromInteger(r)); 083 } 084 085 086 /** 087 * The constructor creates a Complex object with real part 0 and imaginary 088 * part 0. 089 */ 090 public Complex(ComplexRing<C> ring) { 091 this(ring, ring.ring.getZERO()); 092 } 093 094 095 /** 096 * The constructor creates a Complex object from a String representation. 097 * @param s string of a Complex. 098 * @throws NumberFormatException 099 */ 100 public Complex(ComplexRing<C> ring, String s) throws NumberFormatException { 101 this.ring = ring; 102 if (s == null || s.length() == 0) { 103 re = ring.ring.getZERO(); 104 im = ring.ring.getZERO(); 105 return; 106 } 107 s = s.trim(); 108 int i = s.indexOf("i"); 109 if (i < 0) { 110 re = ring.ring.parse(s); 111 im = ring.ring.getZERO(); 112 return; 113 } 114 //logger.warn("String constructor not done"); 115 String sr = ""; 116 if (i > 0) { 117 sr = s.substring(0, i); 118 } 119 String si = ""; 120 if (i < s.length()) { 121 si = s.substring(i + 1, s.length()); 122 } 123 //int j = sr.indexOf("+"); 124 re = ring.ring.parse(sr.trim()); 125 im = ring.ring.parse(si.trim()); 126 } 127 128 129 /** 130 * Get the corresponding element factory. 131 * @return factory for this Element. 132 * @see edu.jas.structure.Element#factory() 133 */ 134 public ComplexRing<C> factory() { 135 return ring; 136 } 137 138 139 /** 140 * Get the real part. 141 * @return re. 142 */ 143 public C getRe() { 144 return re; 145 } 146 147 148 /** 149 * Get the imaginary part. 150 * @return im. 151 */ 152 public C getIm() { 153 return im; 154 } 155 156 157 /** 158 * Clone this. 159 * @see java.lang.Object#clone() 160 */ 161 @Override 162 public Complex<C> copy() { 163 return new Complex<C>(ring, re, im); 164 } 165 166 167 /** 168 * Get the String representation. 169 */ 170 @Override 171 public String toString() { 172 String s = re.toString(); 173 //logger.info("compareTo "+im+" ? 0 = "+i); 174 if (im.isZERO()) { 175 return s; 176 } 177 s += "i" + im; 178 return s; 179 } 180 181 182 /** 183 * Get a scripting compatible string representation. 184 * @return script compatible representation for this Element. 185 * @see edu.jas.structure.Element#toScript() 186 */ 187 //JAVA6only: @Override 188 public String toScript() { 189 // Python case 190 StringBuffer s = new StringBuffer(); 191 if (im.isZERO()) { 192 s.append(re.toScript()); 193 } else { 194 C mi = im; 195 //s.append(""); 196 if (!re.isZERO()) { 197 s.append(re.toScript()); 198 if (mi.signum() > 0) { 199 s.append(" + "); 200 } else { 201 s.append(" - "); 202 mi = mi.negate(); 203 } 204 } 205 if (mi.isONE()) { 206 s.append("I"); 207 } else { 208 s.append(mi.toScript()).append(" * I"); 209 } 210 s.append(""); 211 } 212 return s.toString(); 213 } 214 215 216 /** 217 * Get a scripting compatible string representation of the factory. 218 * @return script compatible representation for this ElemFactory. 219 * @see edu.jas.structure.Element#toScriptFactory() 220 */ 221 //JAVA6only: @Override 222 public String toScriptFactory() { 223 // Python case 224 return ring.toScript(); 225 } 226 227 228 /** 229 * Is Complex number zero. 230 * @return If this is 0 then true is returned, else false. 231 * @see edu.jas.structure.RingElem#isZERO() 232 */ 233 public boolean isZERO() { 234 return re.isZERO() && im.isZERO(); 235 } 236 237 238 /** 239 * Is Complex number one. 240 * @return If this is 1 then true is returned, else false. 241 * @see edu.jas.structure.RingElem#isONE() 242 */ 243 public boolean isONE() { 244 return re.isONE() && im.isZERO(); 245 } 246 247 248 /** 249 * Is Complex imaginary one. 250 * @return If this is i then true is returned, else false. 251 */ 252 public boolean isIMAG() { 253 return re.isZERO() && im.isONE(); 254 } 255 256 257 /** 258 * Is Complex unit element. 259 * @return If this is a unit then true is returned, else false. 260 * @see edu.jas.structure.RingElem#isUnit() 261 */ 262 public boolean isUnit() { 263 if (isZERO()) { 264 return false; 265 } 266 if (ring.isField()) { 267 return true; 268 } 269 return norm().re.isUnit(); 270 } 271 272 273 /** 274 * Comparison with any other object. 275 * @see java.lang.Object#equals(java.lang.Object) 276 */ 277 @Override 278 @SuppressWarnings("unchecked") 279 public boolean equals(Object b) { 280 if (!(b instanceof Complex)) { 281 return false; 282 } 283 Complex<C> bc = null; 284 try { 285 bc = (Complex<C>) b; 286 } catch (ClassCastException e) { 287 } 288 if (bc == null) { 289 return false; 290 } 291 if (!ring.equals(bc.ring)) { 292 return false; 293 } 294 return re.equals(bc.re) && im.equals(bc.im); 295 } 296 297 298 /** 299 * Hash code for this Complex. 300 * @see java.lang.Object#hashCode() 301 */ 302 @Override 303 public int hashCode() { 304 return 37 * re.hashCode() + im.hashCode(); 305 } 306 307 308 /** 309 * Since complex numbers are unordered, we use lexicographical order of re 310 * and im. 311 * @return 0 if this is equal to b; 1 if re > b.re, or re == b.re and im > 312 * b.im; -1 if re < b.re, or re == b.re and im < b.im 313 */ 314 //JAVA6only: @Override 315 public int compareTo(Complex<C> b) { 316 int s = re.compareTo(b.re); 317 if (s != 0) { 318 return s; 319 } 320 return im.compareTo(b.im); 321 } 322 323 324 /** 325 * Since complex numbers are unordered, we use lexicographical order of re 326 * and im. 327 * @return 0 if this is equal to 0; 1 if re > 0, or re == 0 and im > 0; -1 328 * if re < 0, or re == 0 and im < 0 329 * @see edu.jas.structure.RingElem#signum() 330 */ 331 public int signum() { 332 int s = re.signum(); 333 if (s != 0) { 334 return s; 335 } 336 return im.signum(); 337 } 338 339 340 /* arithmetic operations: +, -, - 341 */ 342 343 /** 344 * Complex number summation. 345 * @param B a Complex<C> number. 346 * @return this+B. 347 */ 348 public Complex<C> sum(Complex<C> B) { 349 return new Complex<C>(ring, re.sum(B.re), im.sum(B.im)); 350 } 351 352 353 /** 354 * Complex number subtract. 355 * @param B a Complex<C> number. 356 * @return this-B. 357 */ 358 public Complex<C> subtract(Complex<C> B) { 359 return new Complex<C>(ring, re.subtract(B.re), im.subtract(B.im)); 360 } 361 362 363 /** 364 * Complex number negative. 365 * @return -this. 366 * @see edu.jas.structure.RingElem#negate() 367 */ 368 public Complex<C> negate() { 369 return new Complex<C>(ring, re.negate(), im.negate()); 370 } 371 372 373 /* arithmetic operations: conjugate, absolut value 374 */ 375 376 /** 377 * Complex number conjugate. 378 * @return the complex conjugate of this. 379 */ 380 public Complex<C> conjugate() { 381 return new Complex<C>(ring, re, im.negate()); 382 } 383 384 385 /** 386 * Complex number norm. 387 * @see edu.jas.structure.StarRingElem#norm() 388 * @return ||this||. 389 */ 390 public Complex<C> norm() { 391 // this.conjugate().multiply(this); 392 C v = re.multiply(re); 393 v = v.sum(im.multiply(im)); 394 return new Complex<C>(ring, v); 395 } 396 397 398 /** 399 * Complex number absolute value. 400 * @see edu.jas.structure.RingElem#abs() 401 * @return |this|^2. <b>Note:</b> The square root is not jet implemented. 402 */ 403 public Complex<C> abs() { 404 Complex<C> n = norm(); 405 logger.error("abs() square root missing"); 406 // n = n.sqrt(); 407 return n; 408 } 409 410 411 /* arithmetic operations: *, inverse, / 412 */ 413 414 415 /** 416 * Complex number product. 417 * @param B is a complex number. 418 * @return this*B. 419 */ 420 public Complex<C> multiply(Complex<C> B) { 421 return new Complex<C>(ring, re.multiply(B.re).subtract(im.multiply(B.im)), re.multiply(B.im).sum( 422 im.multiply(B.re))); 423 } 424 425 426 /** 427 * Complex number inverse. 428 * @return S with S*this = 1, if it is defined. 429 * @see edu.jas.structure.RingElem#inverse() 430 */ 431 public Complex<C> inverse() { 432 C a = norm().re.inverse(); 433 return new Complex<C>(ring, re.multiply(a), im.multiply(a.negate())); 434 } 435 436 437 /** 438 * Complex number remainder. 439 * @param S is a complex number. 440 * @return 0. 441 */ 442 public Complex<C> remainder(Complex<C> S) { 443 if (ring.isField()) { 444 return ring.getZERO(); 445 } 446 return quotientRemainder(S)[1]; 447 } 448 449 450 /** 451 * Complex number divide. 452 * @param B is a complex number, non-zero. 453 * @return this/B. 454 */ 455 public Complex<C> divide(Complex<C> B) { 456 if (ring.isField()) { 457 return this.multiply(B.inverse()); 458 } 459 return quotientRemainder(B)[0]; 460 } 461 462 463 /** 464 * Complex number quotient and remainder. 465 * @param S Complex. 466 * @return Complex[] { q, r } with q = this/S and r = rem(this,S). 467 */ 468 @SuppressWarnings("unchecked") 469 public Complex<C>[] quotientRemainder(Complex<C> S) { 470 Complex<C>[] ret = new Complex[2]; 471 C n = S.norm().re; 472 Complex<C> Sp = this.multiply(S.conjugate()); // == this*inv(S)*n 473 C qr = Sp.re.divide(n); 474 C rr = Sp.re.remainder(n); 475 C qi = Sp.im.divide(n); 476 C ri = Sp.im.remainder(n); 477 C rr1 = rr; 478 C ri1 = ri; 479 if (rr.signum() < 0) { 480 rr = rr.negate(); 481 } 482 if (ri.signum() < 0) { 483 ri = ri.negate(); 484 } 485 C one = n.factory().fromInteger(1); 486 if (rr.sum(rr).compareTo(n) > 0) { // rr > n/2 487 if (rr1.signum() < 0) { 488 qr = qr.subtract(one); 489 } else { 490 qr = qr.sum(one); 491 } 492 } 493 if (ri.sum(ri).compareTo(n) > 0) { // ri > n/2 494 if (ri1.signum() < 0) { 495 qi = qi.subtract(one); 496 } else { 497 qi = qi.sum(one); 498 } 499 } 500 Sp = new Complex<C>(ring, qr, qi); 501 Complex<C> Rp = this.subtract(Sp.multiply(S)); 502 if (debug && n.compareTo(Rp.norm().re) < 0) { 503 System.out.println("n = " + n); 504 System.out.println("qr = " + qr); 505 System.out.println("qi = " + qi); 506 System.out.println("rr = " + rr); 507 System.out.println("ri = " + ri); 508 System.out.println("rr1 = " + rr1); 509 System.out.println("ri1 = " + ri1); 510 System.out.println("this = " + this); 511 System.out.println("S = " + S); 512 System.out.println("Sp = " + Sp); 513 BigInteger tr = (BigInteger) (Object) this.re; 514 BigInteger ti = (BigInteger) (Object) this.im; 515 BigInteger sr = (BigInteger) (Object) S.re; 516 BigInteger si = (BigInteger) (Object) S.im; 517 BigComplex tc = new BigComplex(new BigRational(tr), new BigRational(ti)); 518 BigComplex sc = new BigComplex(new BigRational(sr), new BigRational(si)); 519 BigComplex qc = tc.divide(sc); 520 System.out.println("qc = " + qc); 521 BigDecimal qrd = new BigDecimal(qc.getRe()); 522 BigDecimal qid = new BigDecimal(qc.getIm()); 523 System.out.println("qrd = " + qrd); 524 System.out.println("qid = " + qid); 525 throw new ArithmeticException("QR norm not decreasing " + Rp + ", " + Rp.norm()); 526 } 527 ret[0] = Sp; 528 ret[1] = Rp; 529 return ret; 530 } 531 532 533 /** 534 * Complex number quotient and remainder. 535 * @param S Complex. 536 * @return Complex[] { q, r } with q = this/S and r = rem(this,S). 537 * @deprecated use quotientRemainder() 538 */ 539 @Deprecated 540 public Complex<C>[] divideAndRemainder(Complex<C> S) { 541 return quotientRemainder(S); 542 } 543 544 545 /** 546 * Complex number greatest common divisor. 547 * @param S Complex<C>. 548 * @return gcd(this,S). 549 */ 550 public Complex<C> gcd(Complex<C> S) { 551 if (S == null || S.isZERO()) { 552 return this; 553 } 554 if (this.isZERO()) { 555 return S; 556 } 557 if (ring.isField()) { 558 return ring.getONE(); 559 } 560 Complex<C> a = this; 561 Complex<C> b = S; 562 if (a.re.signum() < 0) { 563 a = a.negate(); 564 } 565 if (b.re.signum() < 0) { 566 b = b.negate(); 567 } 568 while (!b.isZERO()) { 569 if (debug) { 570 logger.info("norm(b), a, b = " + b.norm() + ", " + a + ", " + b); 571 } 572 Complex<C>[] qr = a.quotientRemainder(b); 573 if (qr[0].isZERO()) { 574 System.out.println("a = " + a); 575 } 576 a = b; 577 b = qr[1]; 578 } 579 if (a.re.signum() < 0) { 580 a = a.negate(); 581 } 582 return a; 583 } 584 585 586 /** 587 * Complex extended greatest common divisor. 588 * @param S Complex<C>. 589 * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S). 590 */ 591 @SuppressWarnings("unchecked") 592 public Complex<C>[] egcd(Complex<C> S) { 593 Complex<C>[] ret = new Complex[3]; 594 ret[0] = null; 595 ret[1] = null; 596 ret[2] = null; 597 if (S == null || S.isZERO()) { 598 ret[0] = this; 599 return ret; 600 } 601 if (this.isZERO()) { 602 ret[0] = S; 603 return ret; 604 } 605 if (ring.isField()) { 606 Complex<C> half = new Complex<C>(ring, ring.ring.fromInteger(1).divide(ring.ring.fromInteger(2))); 607 ret[0] = ring.getONE(); 608 ret[1] = this.inverse().multiply(half); 609 ret[2] = S.inverse().multiply(half); 610 return ret; 611 } 612 Complex<C>[] qr; 613 Complex<C> q = this; 614 Complex<C> r = S; 615 Complex<C> c1 = ring.getONE(); 616 Complex<C> d1 = ring.getZERO(); 617 Complex<C> c2 = ring.getZERO(); 618 Complex<C> d2 = ring.getONE(); 619 Complex<C> x1; 620 Complex<C> x2; 621 while (!r.isZERO()) { 622 if (debug) { 623 logger.info("norm(r), q, r = " + r.norm() + ", " + q + ", " + r); 624 } 625 qr = q.quotientRemainder(r); 626 q = qr[0]; 627 x1 = c1.subtract(q.multiply(d1)); 628 x2 = c2.subtract(q.multiply(d2)); 629 c1 = d1; 630 c2 = d2; 631 d1 = x1; 632 d2 = x2; 633 q = r; 634 r = qr[1]; 635 } 636 if (q.re.signum() < 0) { 637 q = q.negate(); 638 c1 = c1.negate(); 639 c2 = c2.negate(); 640 } 641 ret[0] = q; 642 ret[1] = c1; 643 ret[2] = c2; 644 return ret; 645 } 646 647}