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}