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}