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