001 /*
002 * $Id: Complex.java 3364 2010-10-24 12:56:06Z kredel $
003 */
004
005 package edu.jas.poly;
006
007
008 import org.apache.log4j.Logger;
009
010 import edu.jas.arith.BigComplex;
011 import edu.jas.arith.BigDecimal;
012 import edu.jas.arith.BigInteger;
013 import edu.jas.arith.BigRational;
014 import edu.jas.structure.GcdRingElem;
015 import edu.jas.structure.RingElem;
016 import 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.
023 * @author Heinz Kredel
024 */
025 public 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> clone() {
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 s.append("");
195 if (!re.isZERO()) {
196 s.append(re.toScript());
197 s.append(" + ");
198 }
199 if (im.isONE()) {
200 s.append("I");
201 } else {
202 s.append(im.toScript()).append(" * I");
203 }
204 s.append("");
205 }
206 return s.toString();
207 }
208
209
210 /**
211 * Get a scripting compatible string representation of the factory.
212 * @return script compatible representation for this ElemFactory.
213 * @see edu.jas.structure.Element#toScriptFactory()
214 */
215 //JAVA6only: @Override
216 public String toScriptFactory() {
217 // Python case
218 return ring.toScript();
219 }
220
221
222 /**
223 * Is Complex number zero.
224 * @return If this is 0 then true is returned, else false.
225 * @see edu.jas.structure.RingElem#isZERO()
226 */
227 public boolean isZERO() {
228 return re.isZERO() && im.isZERO();
229 }
230
231
232 /**
233 * Is Complex number one.
234 * @return If this is 1 then true is returned, else false.
235 * @see edu.jas.structure.RingElem#isONE()
236 */
237 public boolean isONE() {
238 return re.isONE() && im.isZERO();
239 }
240
241
242 /**
243 * Is Complex imaginary one.
244 * @return If this is i then true is returned, else false.
245 */
246 public boolean isIMAG() {
247 return re.isZERO() && im.isONE();
248 }
249
250
251 /**
252 * Is Complex unit element.
253 * @return If this is a unit then true is returned, else false.
254 * @see edu.jas.structure.RingElem#isUnit()
255 */
256 public boolean isUnit() {
257 if (isZERO()) {
258 return false;
259 }
260 if (ring.isField()) {
261 return true;
262 }
263 return norm().re.isUnit();
264 }
265
266
267 /**
268 * Comparison with any other object.
269 * @see java.lang.Object#equals(java.lang.Object)
270 */
271 @Override
272 @SuppressWarnings("unchecked")
273 public boolean equals(Object b) {
274 if (!(b instanceof Complex)) {
275 return false;
276 }
277 Complex<C> bc = null;
278 try {
279 bc = (Complex<C>) b;
280 } catch (ClassCastException e) {
281 }
282 if (bc == null) {
283 return false;
284 }
285 if (!ring.equals(bc.ring)) {
286 return false;
287 }
288 return re.equals(bc.re) && im.equals(bc.im);
289 }
290
291
292 /**
293 * Hash code for this Complex.
294 * @see java.lang.Object#hashCode()
295 */
296 @Override
297 public int hashCode() {
298 return 37 * re.hashCode() + im.hashCode();
299 }
300
301
302 /**
303 * Since complex numbers are unordered, we use lexicographical order of re
304 * and im.
305 * @return 0 if this is equal to b; 1 if re > b.re, or re == b.re and im >
306 * b.im; -1 if re < b.re, or re == b.re and im < b.im
307 */
308 //JAVA6only: @Override
309 public int compareTo(Complex<C> b) {
310 int s = re.compareTo(b.re);
311 if (s != 0) {
312 return s;
313 }
314 return im.compareTo(b.im);
315 }
316
317
318 /**
319 * Since complex numbers are unordered, we use lexicographical order of re
320 * and im.
321 * @return 0 if this is equal to 0; 1 if re > 0, or re == 0 and im > 0; -1
322 * if re < 0, or re == 0 and im < 0
323 * @see edu.jas.structure.RingElem#signum()
324 */
325 public int signum() {
326 int s = re.signum();
327 if (s != 0) {
328 return s;
329 }
330 return im.signum();
331 }
332
333
334 /* arithmetic operations: +, -, -
335 */
336
337 /**
338 * Complex number summation.
339 * @param B a Complex<C> number.
340 * @return this+B.
341 */
342 public Complex<C> sum(Complex<C> B) {
343 return new Complex<C>(ring, re.sum(B.re), im.sum(B.im));
344 }
345
346
347 /**
348 * Complex number subtract.
349 * @param B a Complex<C> number.
350 * @return this-B.
351 */
352 public Complex<C> subtract(Complex<C> B) {
353 return new Complex<C>(ring, re.subtract(B.re), im.subtract(B.im));
354 }
355
356
357 /**
358 * Complex number negative.
359 * @return -this.
360 * @see edu.jas.structure.RingElem#negate()
361 */
362 public Complex<C> negate() {
363 return new Complex<C>(ring, re.negate(), im.negate());
364 }
365
366
367 /* arithmetic operations: conjugate, absolut value
368 */
369
370 /**
371 * Complex number conjugate.
372 * @return the complex conjugate of this.
373 */
374 public Complex<C> conjugate() {
375 return new Complex<C>(ring, re, im.negate());
376 }
377
378
379 /**
380 * Complex number norm.
381 * @see edu.jas.structure.StarRingElem#norm()
382 * @return ||this||.
383 */
384 public Complex<C> norm() {
385 // this.conjugate().multiply(this);
386 C v = re.multiply(re);
387 v = v.sum(im.multiply(im));
388 return new Complex<C>(ring, v);
389 }
390
391
392 /**
393 * Complex number absolute value.
394 * @see edu.jas.structure.RingElem#abs()
395 * @return |this|^2. <b>Note:</b> The square root is not jet implemented.
396 */
397 public Complex<C> abs() {
398 Complex<C> n = norm();
399 logger.error("abs() square root missing");
400 // n = n.sqrt();
401 return n;
402 }
403
404
405 /* arithmetic operations: *, inverse, /
406 */
407
408
409 /**
410 * Complex number product.
411 * @param B is a complex number.
412 * @return this*B.
413 */
414 public Complex<C> multiply(Complex<C> B) {
415 return new Complex<C>(ring, re.multiply(B.re).subtract(im.multiply(B.im)), re.multiply(B.im).sum(
416 im.multiply(B.re)));
417 }
418
419
420 /**
421 * Complex number inverse.
422 * @return S with S*this = 1, if it is defined.
423 * @see edu.jas.structure.RingElem#inverse()
424 */
425 public Complex<C> inverse() {
426 C a = norm().re.inverse();
427 return new Complex<C>(ring, re.multiply(a), im.multiply(a.negate()));
428 }
429
430
431 /**
432 * Complex number remainder.
433 * @param S is a complex number.
434 * @return 0.
435 */
436 public Complex<C> remainder(Complex<C> S) {
437 if (ring.isField()) {
438 return ring.getZERO();
439 }
440 return quotientRemainder(S)[1];
441 }
442
443
444 /**
445 * Complex number divide.
446 * @param B is a complex number, non-zero.
447 * @return this/B.
448 */
449 public Complex<C> divide(Complex<C> B) {
450 if (ring.isField()) {
451 return this.multiply(B.inverse());
452 }
453 return quotientRemainder(B)[0];
454 }
455
456
457 /**
458 * Complex number quotient and remainder.
459 * @param S Complex.
460 * @return Complex[] { q, r } with q = this/S and r = rem(this,S).
461 */
462 @SuppressWarnings("unchecked")
463 public Complex<C>[] quotientRemainder(Complex<C> S) {
464 Complex<C>[] ret = new Complex[2];
465 C n = S.norm().re;
466 Complex<C> Sp = this.multiply(S.conjugate()); // == this*inv(S)*n
467 C qr = Sp.re.divide(n);
468 C rr = Sp.re.remainder(n);
469 C qi = Sp.im.divide(n);
470 C ri = Sp.im.remainder(n);
471 C rr1 = rr;
472 C ri1 = ri;
473 if (rr.signum() < 0) {
474 rr = rr.negate();
475 }
476 if (ri.signum() < 0) {
477 ri = ri.negate();
478 }
479 C one = n.factory().fromInteger(1);
480 if (rr.sum(rr).compareTo(n) > 0) { // rr > n/2
481 if (rr1.signum() < 0) {
482 qr = qr.subtract(one);
483 } else {
484 qr = qr.sum(one);
485 }
486 }
487 if (ri.sum(ri).compareTo(n) > 0) { // ri > n/2
488 if (ri1.signum() < 0) {
489 qi = qi.subtract(one);
490 } else {
491 qi = qi.sum(one);
492 }
493 }
494 Sp = new Complex<C>(ring, qr, qi);
495 Complex<C> Rp = this.subtract(Sp.multiply(S));
496 if (debug && n.compareTo(Rp.norm().re) < 0) {
497 System.out.println("n = " + n);
498 System.out.println("qr = " + qr);
499 System.out.println("qi = " + qi);
500 System.out.println("rr = " + rr);
501 System.out.println("ri = " + ri);
502 System.out.println("rr1 = " + rr1);
503 System.out.println("ri1 = " + ri1);
504 System.out.println("this = " + this);
505 System.out.println("S = " + S);
506 System.out.println("Sp = " + Sp);
507 BigInteger tr = (BigInteger) (Object) this.re;
508 BigInteger ti = (BigInteger) (Object) this.im;
509 BigInteger sr = (BigInteger) (Object) S.re;
510 BigInteger si = (BigInteger) (Object) S.im;
511 BigComplex tc = new BigComplex(new BigRational(tr), new BigRational(ti));
512 BigComplex sc = new BigComplex(new BigRational(sr), new BigRational(si));
513 BigComplex qc = tc.divide(sc);
514 System.out.println("qc = " + qc);
515 BigDecimal qrd = new BigDecimal(qc.getRe());
516 BigDecimal qid = new BigDecimal(qc.getIm());
517 System.out.println("qrd = " + qrd);
518 System.out.println("qid = " + qid);
519 throw new ArithmeticException("QR norm not decreasing " + Rp + ", " + Rp.norm());
520 }
521 ret[0] = Sp;
522 ret[1] = Rp;
523 return ret;
524 }
525
526
527 /**
528 * Complex number quotient and remainder.
529 * @param S Complex.
530 * @return Complex[] { q, r } with q = this/S and r = rem(this,S).
531 * @deprecated use quotientRemainder()
532 */
533 @Deprecated
534 public Complex<C>[] divideAndRemainder(Complex<C> S) {
535 return quotientRemainder(S);
536 }
537
538
539 /**
540 * Complex number greatest common divisor.
541 * @param S Complex<C>.
542 * @return gcd(this,S).
543 */
544 public Complex<C> gcd(Complex<C> S) {
545 if (S == null || S.isZERO()) {
546 return this;
547 }
548 if (this.isZERO()) {
549 return S;
550 }
551 if (ring.isField()) {
552 return ring.getONE();
553 }
554 Complex<C> a = this;
555 Complex<C> b = S;
556 if (a.re.signum() < 0) {
557 a = a.negate();
558 }
559 if (b.re.signum() < 0) {
560 b = b.negate();
561 }
562 while (!b.isZERO()) {
563 if (debug) {
564 logger.info("norm(b), a, b = " + b.norm() + ", " + a + ", " + b);
565 }
566 Complex<C>[] qr = a.quotientRemainder(b);
567 if (qr[0].isZERO()) {
568 System.out.println("a = " + a);
569 }
570 a = b;
571 b = qr[1];
572 }
573 if (a.re.signum() < 0) {
574 a = a.negate();
575 }
576 return a;
577 }
578
579
580 /**
581 * Complex extended greatest common divisor.
582 * @param S Complex<C>.
583 * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S).
584 */
585 @SuppressWarnings("unchecked")
586 public Complex<C>[] egcd(Complex<C> S) {
587 Complex<C>[] ret = new Complex[3];
588 ret[0] = null;
589 ret[1] = null;
590 ret[2] = null;
591 if (S == null || S.isZERO()) {
592 ret[0] = this;
593 return ret;
594 }
595 if (this.isZERO()) {
596 ret[0] = S;
597 return ret;
598 }
599 if (ring.isField()) {
600 Complex<C> half = new Complex<C>(ring, ring.ring.fromInteger(1).divide(ring.ring.fromInteger(2)));
601 ret[0] = ring.getONE();
602 ret[1] = this.inverse().multiply(half);
603 ret[2] = S.inverse().multiply(half);
604 return ret;
605 }
606 Complex<C>[] qr;
607 Complex<C> q = this;
608 Complex<C> r = S;
609 Complex<C> c1 = ring.getONE();
610 Complex<C> d1 = ring.getZERO();
611 Complex<C> c2 = ring.getZERO();
612 Complex<C> d2 = ring.getONE();
613 Complex<C> x1;
614 Complex<C> x2;
615 while (!r.isZERO()) {
616 if (debug) {
617 logger.info("norm(r), q, r = " + r.norm() + ", " + q + ", " + r);
618 }
619 qr = q.quotientRemainder(r);
620 q = qr[0];
621 x1 = c1.subtract(q.multiply(d1));
622 x2 = c2.subtract(q.multiply(d2));
623 c1 = d1;
624 c2 = d2;
625 d1 = x1;
626 d2 = x2;
627 q = r;
628 r = qr[1];
629 }
630 if (q.re.signum() < 0) {
631 q = q.negate();
632 c1 = c1.negate();
633 c2 = c2.negate();
634 }
635 ret[0] = q;
636 ret[1] = c1;
637 ret[2] = c2;
638 return ret;
639 }
640
641 }