001/* 002 * $Id: Complex.java 4956 2014-10-16 22:45:10Z 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.GcdRingElem; 015import edu.jas.structure.RingElem; 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 * Copy this. 159 * @see edu.jas.structure.Element#copy() 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 @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 @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 == null) { 281 return false; 282 } 283 if (!(b instanceof Complex)) { 284 return false; 285 } 286 Complex<C> bc = (Complex<C>) b; 287 if (!ring.equals(bc.ring)) { 288 return false; 289 } 290 return re.equals(bc.re) && im.equals(bc.im); 291 } 292 293 294 /** 295 * Hash code for this Complex. 296 * @see java.lang.Object#hashCode() 297 */ 298 @Override 299 public int hashCode() { 300 return 37 * re.hashCode() + im.hashCode(); 301 } 302 303 304 /** 305 * Since complex numbers are unordered, we use lexicographical order of re 306 * and im. 307 * @return 0 if this is equal to b; 1 if re > b.re, or re == b.re and im > 308 * b.im; -1 if re < b.re, or re == b.re and im < b.im 309 */ 310 @Override 311 public int compareTo(Complex<C> b) { 312 int s = re.compareTo(b.re); 313 if (s != 0) { 314 return s; 315 } 316 return im.compareTo(b.im); 317 } 318 319 320 /** 321 * Since complex numbers are unordered, we use lexicographical order of re 322 * and im. 323 * @return 0 if this is equal to 0; 1 if re > 0, or re == 0 and im > 0; -1 324 * if re < 0, or re == 0 and im < 0 325 * @see edu.jas.structure.RingElem#signum() 326 */ 327 public int signum() { 328 int s = re.signum(); 329 if (s != 0) { 330 return s; 331 } 332 return im.signum(); 333 } 334 335 336 /* arithmetic operations: +, -, - 337 */ 338 339 /** 340 * Complex number summation. 341 * @param B a Complex<C> number. 342 * @return this+B. 343 */ 344 public Complex<C> sum(Complex<C> B) { 345 return new Complex<C>(ring, re.sum(B.re), im.sum(B.im)); 346 } 347 348 349 /** 350 * Complex number subtract. 351 * @param B a Complex<C> number. 352 * @return this-B. 353 */ 354 public Complex<C> subtract(Complex<C> B) { 355 return new Complex<C>(ring, re.subtract(B.re), im.subtract(B.im)); 356 } 357 358 359 /** 360 * Complex number negative. 361 * @return -this. 362 * @see edu.jas.structure.RingElem#negate() 363 */ 364 public Complex<C> negate() { 365 return new Complex<C>(ring, re.negate(), im.negate()); 366 } 367 368 369 /* arithmetic operations: conjugate, absolut value 370 */ 371 372 /** 373 * Complex number conjugate. 374 * @return the complex conjugate of this. 375 */ 376 public Complex<C> conjugate() { 377 return new Complex<C>(ring, re, im.negate()); 378 } 379 380 381 /** 382 * Complex number norm. 383 * @see edu.jas.structure.StarRingElem#norm() 384 * @return ||this||. 385 */ 386 public Complex<C> norm() { 387 // this.conjugate().multiply(this); 388 C v = re.multiply(re); 389 v = v.sum(im.multiply(im)); 390 return new Complex<C>(ring, v); 391 } 392 393 394 /** 395 * Complex number absolute value. 396 * @see edu.jas.structure.RingElem#abs() 397 * @return |this|^2. <b>Note:</b> The square root is not jet implemented. 398 */ 399 public Complex<C> abs() { 400 Complex<C> n = norm(); 401 logger.error("abs() square root missing"); 402 // n = n.sqrt(); 403 return n; 404 } 405 406 407 /* arithmetic operations: *, inverse, / 408 */ 409 410 411 /** 412 * Complex number product. 413 * @param B is a complex number. 414 * @return this*B. 415 */ 416 public Complex<C> multiply(Complex<C> B) { 417 return new Complex<C>(ring, re.multiply(B.re).subtract(im.multiply(B.im)), re.multiply(B.im).sum( 418 im.multiply(B.re))); 419 } 420 421 422 /** 423 * Complex number inverse. 424 * @return S with S*this = 1, if it is defined. 425 * @see edu.jas.structure.RingElem#inverse() 426 */ 427 public Complex<C> inverse() { 428 C a = norm().re.inverse(); 429 return new Complex<C>(ring, re.multiply(a), im.multiply(a.negate())); 430 } 431 432 433 /** 434 * Complex number remainder. 435 * @param S is a complex number. 436 * @return 0. 437 */ 438 public Complex<C> remainder(Complex<C> S) { 439 if (ring.isField()) { 440 return ring.getZERO(); 441 } 442 return quotientRemainder(S)[1]; 443 } 444 445 446 /** 447 * Complex number divide. 448 * @param B is a complex number, non-zero. 449 * @return this/B. 450 */ 451 public Complex<C> divide(Complex<C> B) { 452 if (ring.isField()) { 453 return this.multiply(B.inverse()); 454 } 455 return quotientRemainder(B)[0]; 456 } 457 458 459 /** 460 * Complex number quotient and remainder. 461 * @param S Complex. 462 * @return Complex[] { q, r } with q = this/S and r = rem(this,S). 463 */ 464 @SuppressWarnings({ "unchecked", "cast" }) 465 public Complex<C>[] quotientRemainder(Complex<C> S) { 466 Complex<C>[] ret = new Complex[2]; 467 C n = S.norm().re; 468 Complex<C> Sp = this.multiply(S.conjugate()); // == this*inv(S)*n 469 C qr = Sp.re.divide(n); 470 C rr = Sp.re.remainder(n); 471 C qi = Sp.im.divide(n); 472 C ri = Sp.im.remainder(n); 473 C rr1 = rr; 474 C ri1 = ri; 475 if (rr.signum() < 0) { 476 rr = rr.negate(); 477 } 478 if (ri.signum() < 0) { 479 ri = ri.negate(); 480 } 481 C one = n.factory().fromInteger(1); 482 if (rr.sum(rr).compareTo(n) > 0) { // rr > n/2 483 if (rr1.signum() < 0) { 484 qr = qr.subtract(one); 485 } else { 486 qr = qr.sum(one); 487 } 488 } 489 if (ri.sum(ri).compareTo(n) > 0) { // ri > n/2 490 if (ri1.signum() < 0) { 491 qi = qi.subtract(one); 492 } else { 493 qi = qi.sum(one); 494 } 495 } 496 Sp = new Complex<C>(ring, qr, qi); 497 Complex<C> Rp = this.subtract(Sp.multiply(S)); 498 if (debug && n.compareTo(Rp.norm().re) < 0) { 499 System.out.println("n = " + n); 500 System.out.println("qr = " + qr); 501 System.out.println("qi = " + qi); 502 System.out.println("rr = " + rr); 503 System.out.println("ri = " + ri); 504 System.out.println("rr1 = " + rr1); 505 System.out.println("ri1 = " + ri1); 506 System.out.println("this = " + this); 507 System.out.println("S = " + S); 508 System.out.println("Sp = " + Sp); 509 BigInteger tr = (BigInteger) (Object) this.re; 510 BigInteger ti = (BigInteger) (Object) this.im; 511 BigInteger sr = (BigInteger) (Object) S.re; 512 BigInteger si = (BigInteger) (Object) S.im; 513 BigComplex tc = new BigComplex(new BigRational(tr), new BigRational(ti)); 514 BigComplex sc = new BigComplex(new BigRational(sr), new BigRational(si)); 515 BigComplex qc = tc.divide(sc); 516 System.out.println("qc = " + qc); 517 BigDecimal qrd = new BigDecimal(qc.getRe()); 518 BigDecimal qid = new BigDecimal(qc.getIm()); 519 System.out.println("qrd = " + qrd); 520 System.out.println("qid = " + qid); 521 throw new ArithmeticException("QR norm not decreasing " + Rp + ", " + Rp.norm()); 522 } 523 ret[0] = Sp; 524 ret[1] = Rp; 525 return ret; 526 } 527 528 529 /** 530 * Complex number quotient and remainder. 531 * @param S Complex. 532 * @return Complex[] { q, r } with q = this/S and r = rem(this,S). 533 * @deprecated use quotientRemainder() 534 */ 535 @Deprecated 536 public Complex<C>[] divideAndRemainder(Complex<C> S) { 537 return quotientRemainder(S); 538 } 539 540 541 /** 542 * Complex number greatest common divisor. 543 * @param S Complex<C>. 544 * @return gcd(this,S). 545 */ 546 public Complex<C> gcd(Complex<C> S) { 547 if (S == null || S.isZERO()) { 548 return this; 549 } 550 if (this.isZERO()) { 551 return S; 552 } 553 if (ring.isField()) { 554 return ring.getONE(); 555 } 556 Complex<C> a = this; 557 Complex<C> b = S; 558 if (a.re.signum() < 0) { 559 a = a.negate(); 560 } 561 if (b.re.signum() < 0) { 562 b = b.negate(); 563 } 564 while (!b.isZERO()) { 565 if (debug) { 566 logger.info("norm(b), a, b = " + b.norm() + ", " + a + ", " + b); 567 } 568 Complex<C>[] qr = a.quotientRemainder(b); 569 if (qr[0].isZERO()) { 570 System.out.println("a = " + a); 571 } 572 a = b; 573 b = qr[1]; 574 } 575 if (a.re.signum() < 0) { 576 a = a.negate(); 577 } 578 return a; 579 } 580 581 582 /** 583 * Complex extended greatest common divisor. 584 * @param S Complex<C>. 585 * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S). 586 */ 587 @SuppressWarnings("unchecked") 588 public Complex<C>[] egcd(Complex<C> S) { 589 Complex<C>[] ret = new Complex[3]; 590 ret[0] = null; 591 ret[1] = null; 592 ret[2] = null; 593 if (S == null || S.isZERO()) { 594 ret[0] = this; 595 return ret; 596 } 597 if (this.isZERO()) { 598 ret[0] = S; 599 return ret; 600 } 601 if (ring.isField()) { 602 Complex<C> half = new Complex<C>(ring, ring.ring.fromInteger(1).divide(ring.ring.fromInteger(2))); 603 ret[0] = ring.getONE(); 604 ret[1] = this.inverse().multiply(half); 605 ret[2] = S.inverse().multiply(half); 606 return ret; 607 } 608 Complex<C>[] qr; 609 Complex<C> q = this; 610 Complex<C> r = S; 611 Complex<C> c1 = ring.getONE(); 612 Complex<C> d1 = ring.getZERO(); 613 Complex<C> c2 = ring.getZERO(); 614 Complex<C> d2 = ring.getONE(); 615 Complex<C> x1; 616 Complex<C> x2; 617 while (!r.isZERO()) { 618 if (debug) { 619 logger.info("norm(r), q, r = " + r.norm() + ", " + q + ", " + r); 620 } 621 qr = q.quotientRemainder(r); 622 q = qr[0]; 623 x1 = c1.subtract(q.multiply(d1)); 624 x2 = c2.subtract(q.multiply(d2)); 625 c1 = d1; 626 c2 = d2; 627 d1 = x1; 628 d2 = x2; 629 q = r; 630 r = qr[1]; 631 } 632 if (q.re.signum() < 0) { 633 q = q.negate(); 634 c1 = c1.negate(); 635 c2 = c2.negate(); 636 } 637 ret[0] = q; 638 ret[1] = c1; 639 ret[2] = c2; 640 return ret; 641 } 642 643}