Source code

001/*
002 * $Id$
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        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        re = ring.ring.parse(sr.trim());
124        im = ring.ring.parse(si.trim());
125    }
126
127
128    /**
129     * Get the corresponding element factory.
130     * @return factory for this Element.
131     * @see edu.jas.structure.Element#factory()
132     */
133    public ComplexRing<C> factory() {
134        return ring;
135    }
136
137
138    /**
139     * Get the real part.
140     * @return re.
141     */
142    public C getRe() {
143        return re;
144    }
145
146
147    /**
148     * Get the imaginary part.
149     * @return im.
150     */
151    public C getIm() {
152        return im;
153    }
154
155
156    /**
157     * Copy this.
158     * @see edu.jas.structure.Element#copy()
159     */
160    @Override
161    public Complex<C> copy() {
162        return new Complex<C>(ring, re, im);
163    }
164
165
166    /**
167     * Get the String representation.
168     */
169    @Override
170    public String toString() {
171        String s = re.toString();
172        if (im.isZERO()) {
173            return s;
174        }
175        s += "i" + im;
176        return s;
177    }
178
179
180    /**
181     * Get a scripting compatible string representation.
182     * @return script compatible representation for this Element.
183     * @see edu.jas.structure.Element#toScript()
184     */
185    @Override
186    public String toScript() {
187        // Python case
188        StringBuffer s = new StringBuffer();
189        if (im.isZERO()) {
190            s.append(re.toScript());
191        } else {
192            C mi = im;
193            //s.append("");
194            if (!re.isZERO()) {
195                s.append(re.toScript());
196                if (mi.signum() > 0) {
197                    s.append(" + ");
198                } else {
199                    s.append(" - ");
200                    mi = mi.negate();
201                }
202            }
203            if (mi.isONE()) {
204                s.append("I");
205            } else {
206                s.append(mi.toScript()).append(" * I");
207            }
208            s.append("");
209        }
210        return s.toString();
211    }
212
213
214    /**
215     * Get a scripting compatible string representation of the factory.
216     * @return script compatible representation for this ElemFactory.
217     * @see edu.jas.structure.Element#toScriptFactory()
218     */
219    @Override
220    public String toScriptFactory() {
221        // Python case
222        return ring.toScript();
223    }
224
225
226    /**
227     * Is Complex number zero.
228     * @return If this is 0 then true is returned, else false.
229     * @see edu.jas.structure.RingElem#isZERO()
230     */
231    public boolean isZERO() {
232        return re.isZERO() && im.isZERO();
233    }
234
235
236    /**
237     * Is Complex number one.
238     * @return If this is 1 then true is returned, else false.
239     * @see edu.jas.structure.RingElem#isONE()
240     */
241    public boolean isONE() {
242        return re.isONE() && im.isZERO();
243    }
244
245
246    /**
247     * Is Complex imaginary one.
248     * @return If this is i then true is returned, else false.
249     */
250    public boolean isIMAG() {
251        return re.isZERO() && im.isONE();
252    }
253
254
255    /**
256     * Is Complex unit element.
257     * @return If this is a unit then true is returned, else false.
258     * @see edu.jas.structure.RingElem#isUnit()
259     */
260    public boolean isUnit() {
261        if (isZERO()) {
262            return false;
263        }
264        if (ring.isField()) {
265            return true;
266        }
267        return norm().re.isUnit();
268    }
269
270
271    /**
272     * Comparison with any other object.
273     * @see java.lang.Object#equals(java.lang.Object)
274     */
275    @Override
276    @SuppressWarnings("unchecked")
277    public boolean equals(Object b) {
278        if (b == null) {
279            return false;
280        }
281        if (!(b instanceof Complex)) {
282            return false;
283        }
284        Complex<C> bc = (Complex<C>) b;
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 &gt; b.re, or re == b.re and im
306     *         &gt; b.im; -1 if re &lt; b.re, or re == b.re and im &lt; b.im
307     */
308    @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 &gt; 0, or re == 0 and im &gt;
322     *         0; -1 if re &lt; 0, or re == 0 and im &lt; 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, absolute 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.multiply(this.conjugate());
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)),
416                        re.multiply(B.im).sum(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", "cast" })
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 greatest common divisor.
529     * @param S Complex<C>.
530     * @return gcd(this,S).
531     */
532    public Complex<C> gcd(Complex<C> S) {
533        if (S == null || S.isZERO()) {
534            return this;
535        }
536        if (this.isZERO()) {
537            return S;
538        }
539        if (ring.isField()) {
540            return ring.getONE();
541        }
542        Complex<C> a = this;
543        Complex<C> b = S;
544        if (a.re.signum() < 0) {
545            a = a.negate();
546        }
547        if (b.re.signum() < 0) {
548            b = b.negate();
549        }
550        while (!b.isZERO()) {
551            if (debug) {
552                logger.info("norm(b), a, b = {}, {}, {}", b.norm(), a, b);
553            }
554            Complex<C>[] qr = a.quotientRemainder(b);
555            if (qr[0].isZERO()) {
556                System.out.println("a = " + a);
557            }
558            a = b;
559            b = qr[1];
560        }
561        if (a.re.signum() < 0) {
562            a = a.negate();
563        }
564        return a;
565    }
566
567
568    /**
569     * Complex extended greatest common divisor.
570     * @param S Complex<C>.
571     * @return [ gcd(this,S), a, b ] with a*this + b*S = gcd(this,S).
572     */
573    @SuppressWarnings("unchecked")
574    public Complex<C>[] egcd(Complex<C> S) {
575        Complex<C>[] ret = new Complex[3];
576        ret[0] = null;
577        ret[1] = null;
578        ret[2] = null;
579        if (S == null || S.isZERO()) {
580            ret[0] = this;
581            return ret;
582        }
583        if (this.isZERO()) {
584            ret[0] = S;
585            return ret;
586        }
587        if (ring.isField()) {
588            Complex<C> half = new Complex<C>(ring, ring.ring.fromInteger(1).divide(ring.ring.fromInteger(2)));
589            ret[0] = ring.getONE();
590            ret[1] = this.inverse().multiply(half);
591            ret[2] = S.inverse().multiply(half);
592            return ret;
593        }
594        Complex<C>[] qr;
595        Complex<C> q = this;
596        Complex<C> r = S;
597        Complex<C> c1 = ring.getONE();
598        Complex<C> d1 = ring.getZERO();
599        Complex<C> c2 = ring.getZERO();
600        Complex<C> d2 = ring.getONE();
601        Complex<C> x1;
602        Complex<C> x2;
603        while (!r.isZERO()) {
604            if (debug) {
605                logger.info("norm(r), q, r = {}, {}, {}", r.norm(), q, r);
606            }
607            qr = q.quotientRemainder(r);
608            q = qr[0];
609            x1 = c1.subtract(q.multiply(d1));
610            x2 = c2.subtract(q.multiply(d2));
611            c1 = d1;
612            c2 = d2;
613            d1 = x1;
614            d2 = x2;
615            q = r;
616            r = qr[1];
617        }
618        if (q.re.signum() < 0) {
619            q = q.negate();
620            c1 = c1.negate();
621            c2 = c2.negate();
622        }
623        ret[0] = q;
624        ret[1] = c1;
625        ret[2] = c2;
626        return ret;
627    }
628
629}