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}