1 /*
2  * Licensed to the Apache Software Foundation (ASF) under one or more
3  * contributor license agreements.  See the NOTICE file distributed with
4  * this work for additional information regarding copyright ownership.
5  * The ASF licenses this file to You under the Apache License, Version 2.0
6  * (the "License"); you may not use this file except in compliance with
7  * the License.  You may obtain a copy of the License at
8  *
9  *      http://www.apache.org/licenses/LICENSE-2.0
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  */
17 
18 package org.apache.commons.math3.complex;
19 
20 import java.io.Serializable;
21 import java.util.ArrayList;
22 import java.util.List;
23 
24 import org.apache.commons.math3.FieldElement;
25 import org.apache.commons.math3.exception.NotPositiveException;
26 import org.apache.commons.math3.exception.NullArgumentException;
27 import org.apache.commons.math3.exception.util.LocalizedFormats;
28 import org.apache.commons.math3.util.FastMath;
29 import org.apache.commons.math3.util.MathUtils;
30 import org.apache.commons.math3.util.Precision;
31 
32 /**
33  * Representation of a Complex number, i.e. a number which has both a
34  * real and imaginary part.
35  * <p>
36  * Implementations of arithmetic operations handle {@code NaN} and
37  * infinite values according to the rules for {@link java.lang.Double}, i.e.
38  * {@link #equals} is an equivalence relation for all instances that have
39  * a {@code NaN} in either real or imaginary part, e.g. the following are
40  * considered equal:
41  * <ul>
42  *  <li>{@code 1 + NaNi}</li>
43  *  <li>{@code NaN + i}</li>
44  *  <li>{@code NaN + NaNi}</li>
45  * </ul><p>
46  * Note that this contradicts the IEEE-754 standard for floating
47  * point numbers (according to which the test {@code x == x} must fail if
48  * {@code x} is {@code NaN}). The method
49  * {@link org.apache.commons.math3.util.Precision#equals(double,double,int)
50  * equals for primitive double} in {@link org.apache.commons.math3.util.Precision}
51  * conforms with IEEE-754 while this class conforms with the standard behavior
52  * for Java object types.</p>
53  *
54  */
55 public class Complex implements FieldElement<Complex>, Serializable  {
56     /** The square root of -1. A number representing "0.0 + 1.0i" */
57     public static final Complex I = new Complex(0.0, 1.0);
58     // CHECKSTYLE: stop ConstantName
59     /** A complex number representing "NaN + NaNi" */
60     public static final Complex NaN = new Complex(Double.NaN, Double.NaN);
61     // CHECKSTYLE: resume ConstantName
62     /** A complex number representing "+INF + INFi" */
63     public static final Complex INF = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
64     /** A complex number representing "1.0 + 0.0i" */
65     public static final Complex ONE = new Complex(1.0, 0.0);
66     /** A complex number representing "0.0 + 0.0i" */
67     public static final Complex ZERO = new Complex(0.0, 0.0);
68 
69     /** Serializable version identifier */
70     private static final long serialVersionUID = -6195664516687396620L;
71 
72     /** The imaginary part. */
73     private final double imaginary;
74     /** The real part. */
75     private final double real;
76     /** Record whether this complex number is equal to NaN. */
77     private final transient boolean isNaN;
78     /** Record whether this complex number is infinite. */
79     private final transient boolean isInfinite;
80 
81     /**
82      * Create a complex number given only the real part.
83      *
84      * @param real Real part.
85      */
Complex(double real)86     public Complex(double real) {
87         this(real, 0.0);
88     }
89 
90     /**
91      * Create a complex number given the real and imaginary parts.
92      *
93      * @param real Real part.
94      * @param imaginary Imaginary part.
95      */
Complex(double real, double imaginary)96     public Complex(double real, double imaginary) {
97         this.real = real;
98         this.imaginary = imaginary;
99 
100         isNaN = Double.isNaN(real) || Double.isNaN(imaginary);
101         isInfinite = !isNaN &&
102             (Double.isInfinite(real) || Double.isInfinite(imaginary));
103     }
104 
105     /**
106      * Return the absolute value of this complex number.
107      * Returns {@code NaN} if either real or imaginary part is {@code NaN}
108      * and {@code Double.POSITIVE_INFINITY} if neither part is {@code NaN},
109      * but at least one part is infinite.
110      *
111      * @return the absolute value.
112      */
abs()113     public double abs() {
114         if (isNaN) {
115             return Double.NaN;
116         }
117         if (isInfinite()) {
118             return Double.POSITIVE_INFINITY;
119         }
120         if (FastMath.abs(real) < FastMath.abs(imaginary)) {
121             if (imaginary == 0.0) {
122                 return FastMath.abs(real);
123             }
124             double q = real / imaginary;
125             return FastMath.abs(imaginary) * FastMath.sqrt(1 + q * q);
126         } else {
127             if (real == 0.0) {
128                 return FastMath.abs(imaginary);
129             }
130             double q = imaginary / real;
131             return FastMath.abs(real) * FastMath.sqrt(1 + q * q);
132         }
133     }
134 
135     /**
136      * Returns a {@code Complex} whose value is
137      * {@code (this + addend)}.
138      * Uses the definitional formula
139      * <p>
140      *   {@code (a + bi) + (c + di) = (a+c) + (b+d)i}
141      * </p>
142      * If either {@code this} or {@code addend} has a {@code NaN} value in
143      * either part, {@link #NaN} is returned; otherwise {@code Infinite}
144      * and {@code NaN} values are returned in the parts of the result
145      * according to the rules for {@link java.lang.Double} arithmetic.
146      *
147      * @param  addend Value to be added to this {@code Complex}.
148      * @return {@code this + addend}.
149      * @throws NullArgumentException if {@code addend} is {@code null}.
150      */
add(Complex addend)151     public Complex add(Complex addend) throws NullArgumentException {
152         MathUtils.checkNotNull(addend);
153         if (isNaN || addend.isNaN) {
154             return NaN;
155         }
156 
157         return createComplex(real + addend.getReal(),
158                              imaginary + addend.getImaginary());
159     }
160 
161     /**
162      * Returns a {@code Complex} whose value is {@code (this + addend)},
163      * with {@code addend} interpreted as a real number.
164      *
165      * @param addend Value to be added to this {@code Complex}.
166      * @return {@code this + addend}.
167      * @see #add(Complex)
168      */
add(double addend)169     public Complex add(double addend) {
170         if (isNaN || Double.isNaN(addend)) {
171             return NaN;
172         }
173 
174         return createComplex(real + addend, imaginary);
175     }
176 
177      /**
178      * Returns the conjugate of this complex number.
179      * The conjugate of {@code a + bi} is {@code a - bi}.
180      * <p>
181      * {@link #NaN} is returned if either the real or imaginary
182      * part of this Complex number equals {@code Double.NaN}.
183      * </p><p>
184      * If the imaginary part is infinite, and the real part is not
185      * {@code NaN}, the returned value has infinite imaginary part
186      * of the opposite sign, e.g. the conjugate of
187      * {@code 1 + POSITIVE_INFINITY i} is {@code 1 - NEGATIVE_INFINITY i}.
188      * </p>
189      * @return the conjugate of this Complex object.
190      */
conjugate()191     public Complex conjugate() {
192         if (isNaN) {
193             return NaN;
194         }
195 
196         return createComplex(real, -imaginary);
197     }
198 
199     /**
200      * Returns a {@code Complex} whose value is
201      * {@code (this / divisor)}.
202      * Implements the definitional formula
203      * <pre>
204      *  <code>
205      *    a + bi          ac + bd + (bc - ad)i
206      *    ----------- = -------------------------
207      *    c + di         c<sup>2</sup> + d<sup>2</sup>
208      *  </code>
209      * </pre>
210      * but uses
211      * <a href="http://doi.acm.org/10.1145/1039813.1039814">
212      * prescaling of operands</a> to limit the effects of overflows and
213      * underflows in the computation.
214      * <p>
215      * {@code Infinite} and {@code NaN} values are handled according to the
216      * following rules, applied in the order presented:
217      * <ul>
218      *  <li>If either {@code this} or {@code divisor} has a {@code NaN} value
219      *   in either part, {@link #NaN} is returned.
220      *  </li>
221      *  <li>If {@code divisor} equals {@link #ZERO}, {@link #NaN} is returned.
222      *  </li>
223      *  <li>If {@code this} and {@code divisor} are both infinite,
224      *   {@link #NaN} is returned.
225      *  </li>
226      *  <li>If {@code this} is finite (i.e., has no {@code Infinite} or
227      *   {@code NaN} parts) and {@code divisor} is infinite (one or both parts
228      *   infinite), {@link #ZERO} is returned.
229      *  </li>
230      *  <li>If {@code this} is infinite and {@code divisor} is finite,
231      *   {@code NaN} values are returned in the parts of the result if the
232      *   {@link java.lang.Double} rules applied to the definitional formula
233      *   force {@code NaN} results.
234      *  </li>
235      * </ul>
236      *
237      * @param divisor Value by which this {@code Complex} is to be divided.
238      * @return {@code this / divisor}.
239      * @throws NullArgumentException if {@code divisor} is {@code null}.
240      */
divide(Complex divisor)241     public Complex divide(Complex divisor)
242         throws NullArgumentException {
243         MathUtils.checkNotNull(divisor);
244         if (isNaN || divisor.isNaN) {
245             return NaN;
246         }
247 
248         final double c = divisor.getReal();
249         final double d = divisor.getImaginary();
250         if (c == 0.0 && d == 0.0) {
251             return NaN;
252         }
253 
254         if (divisor.isInfinite() && !isInfinite()) {
255             return ZERO;
256         }
257 
258         if (FastMath.abs(c) < FastMath.abs(d)) {
259             double q = c / d;
260             double denominator = c * q + d;
261             return createComplex((real * q + imaginary) / denominator,
262                 (imaginary * q - real) / denominator);
263         } else {
264             double q = d / c;
265             double denominator = d * q + c;
266             return createComplex((imaginary * q + real) / denominator,
267                 (imaginary - real * q) / denominator);
268         }
269     }
270 
271     /**
272      * Returns a {@code Complex} whose value is {@code (this / divisor)},
273      * with {@code divisor} interpreted as a real number.
274      *
275      * @param  divisor Value by which this {@code Complex} is to be divided.
276      * @return {@code this / divisor}.
277      * @see #divide(Complex)
278      */
divide(double divisor)279     public Complex divide(double divisor) {
280         if (isNaN || Double.isNaN(divisor)) {
281             return NaN;
282         }
283         if (divisor == 0d) {
284             return NaN;
285         }
286         if (Double.isInfinite(divisor)) {
287             return !isInfinite() ? ZERO : NaN;
288         }
289         return createComplex(real / divisor,
290                              imaginary  / divisor);
291     }
292 
293     /** {@inheritDoc} */
reciprocal()294     public Complex reciprocal() {
295         if (isNaN) {
296             return NaN;
297         }
298 
299         if (real == 0.0 && imaginary == 0.0) {
300             return INF;
301         }
302 
303         if (isInfinite) {
304             return ZERO;
305         }
306 
307         if (FastMath.abs(real) < FastMath.abs(imaginary)) {
308             double q = real / imaginary;
309             double scale = 1. / (real * q + imaginary);
310             return createComplex(scale * q, -scale);
311         } else {
312             double q = imaginary / real;
313             double scale = 1. / (imaginary * q + real);
314             return createComplex(scale, -scale * q);
315         }
316     }
317 
318     /**
319      * Test for equality with another object.
320      * If both the real and imaginary parts of two complex numbers
321      * are exactly the same, and neither is {@code Double.NaN}, the two
322      * Complex objects are considered to be equal.
323      * The behavior is the same as for JDK's {@link Double#equals(Object)
324      * Double}:
325      * <ul>
326      *  <li>All {@code NaN} values are considered to be equal,
327      *   i.e, if either (or both) real and imaginary parts of the complex
328      *   number are equal to {@code Double.NaN}, the complex number is equal
329      *   to {@code NaN}.
330      *  </li>
331      *  <li>
332      *   Instances constructed with different representations of zero (i.e.
333      *   either "0" or "-0") are <em>not</em> considered to be equal.
334      *  </li>
335      * </ul>
336      *
337      * @param other Object to test for equality with this instance.
338      * @return {@code true} if the objects are equal, {@code false} if object
339      * is {@code null}, not an instance of {@code Complex}, or not equal to
340      * this instance.
341      */
342     @Override
equals(Object other)343     public boolean equals(Object other) {
344         if (this == other) {
345             return true;
346         }
347         if (other instanceof Complex){
348             Complex c = (Complex) other;
349             if (c.isNaN) {
350                 return isNaN;
351             } else {
352                 return MathUtils.equals(real, c.real) &&
353                     MathUtils.equals(imaginary, c.imaginary);
354             }
355         }
356         return false;
357     }
358 
359     /**
360      * Test for the floating-point equality between Complex objects.
361      * It returns {@code true} if both arguments are equal or within the
362      * range of allowed error (inclusive).
363      *
364      * @param x First value (cannot be {@code null}).
365      * @param y Second value (cannot be {@code null}).
366      * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
367      * values between the real (resp. imaginary) parts of {@code x} and
368      * {@code y}.
369      * @return {@code true} if there are fewer than {@code maxUlps} floating
370      * point values between the real (resp. imaginary) parts of {@code x}
371      * and {@code y}.
372      *
373      * @see Precision#equals(double,double,int)
374      * @since 3.3
375      */
equals(Complex x, Complex y, int maxUlps)376     public static boolean equals(Complex x, Complex y, int maxUlps) {
377         return Precision.equals(x.real, y.real, maxUlps) &&
378             Precision.equals(x.imaginary, y.imaginary, maxUlps);
379     }
380 
381     /**
382      * Returns {@code true} iff the values are equal as defined by
383      * {@link #equals(Complex,Complex,int) equals(x, y, 1)}.
384      *
385      * @param x First value (cannot be {@code null}).
386      * @param y Second value (cannot be {@code null}).
387      * @return {@code true} if the values are equal.
388      *
389      * @since 3.3
390      */
equals(Complex x, Complex y)391     public static boolean equals(Complex x, Complex y) {
392         return equals(x, y, 1);
393     }
394 
395     /**
396      * Returns {@code true} if, both for the real part and for the imaginary
397      * part, there is no double value strictly between the arguments or the
398      * difference between them is within the range of allowed error
399      * (inclusive).  Returns {@code false} if either of the arguments is NaN.
400      *
401      * @param x First value (cannot be {@code null}).
402      * @param y Second value (cannot be {@code null}).
403      * @param eps Amount of allowed absolute error.
404      * @return {@code true} if the values are two adjacent floating point
405      * numbers or they are within range of each other.
406      *
407      * @see Precision#equals(double,double,double)
408      * @since 3.3
409      */
equals(Complex x, Complex y, double eps)410     public static boolean equals(Complex x, Complex y, double eps) {
411         return Precision.equals(x.real, y.real, eps) &&
412             Precision.equals(x.imaginary, y.imaginary, eps);
413     }
414 
415     /**
416      * Returns {@code true} if, both for the real part and for the imaginary
417      * part, there is no double value strictly between the arguments or the
418      * relative difference between them is smaller or equal to the given
419      * tolerance. Returns {@code false} if either of the arguments is NaN.
420      *
421      * @param x First value (cannot be {@code null}).
422      * @param y Second value (cannot be {@code null}).
423      * @param eps Amount of allowed relative error.
424      * @return {@code true} if the values are two adjacent floating point
425      * numbers or they are within range of each other.
426      *
427      * @see Precision#equalsWithRelativeTolerance(double,double,double)
428      * @since 3.3
429      */
equalsWithRelativeTolerance(Complex x, Complex y, double eps)430     public static boolean equalsWithRelativeTolerance(Complex x, Complex y,
431                                                       double eps) {
432         return Precision.equalsWithRelativeTolerance(x.real, y.real, eps) &&
433             Precision.equalsWithRelativeTolerance(x.imaginary, y.imaginary, eps);
434     }
435 
436     /**
437      * Get a hashCode for the complex number.
438      * Any {@code Double.NaN} value in real or imaginary part produces
439      * the same hash code {@code 7}.
440      *
441      * @return a hash code value for this object.
442      */
443     @Override
hashCode()444     public int hashCode() {
445         if (isNaN) {
446             return 7;
447         }
448         return 37 * (17 * MathUtils.hash(imaginary) +
449             MathUtils.hash(real));
450     }
451 
452     /**
453      * Access the imaginary part.
454      *
455      * @return the imaginary part.
456      */
getImaginary()457     public double getImaginary() {
458         return imaginary;
459     }
460 
461     /**
462      * Access the real part.
463      *
464      * @return the real part.
465      */
getReal()466     public double getReal() {
467         return real;
468     }
469 
470     /**
471      * Checks whether either or both parts of this complex number is
472      * {@code NaN}.
473      *
474      * @return true if either or both parts of this complex number is
475      * {@code NaN}; false otherwise.
476      */
isNaN()477     public boolean isNaN() {
478         return isNaN;
479     }
480 
481     /**
482      * Checks whether either the real or imaginary part of this complex number
483      * takes an infinite value (either {@code Double.POSITIVE_INFINITY} or
484      * {@code Double.NEGATIVE_INFINITY}) and neither part
485      * is {@code NaN}.
486      *
487      * @return true if one or both parts of this complex number are infinite
488      * and neither part is {@code NaN}.
489      */
isInfinite()490     public boolean isInfinite() {
491         return isInfinite;
492     }
493 
494     /**
495      * Returns a {@code Complex} whose value is {@code this * factor}.
496      * Implements preliminary checks for {@code NaN} and infinity followed by
497      * the definitional formula:
498      * <p>
499      *   {@code (a + bi)(c + di) = (ac - bd) + (ad + bc)i}
500      * </p>
501      * Returns {@link #NaN} if either {@code this} or {@code factor} has one or
502      * more {@code NaN} parts.
503      * <p>
504      * Returns {@link #INF} if neither {@code this} nor {@code factor} has one
505      * or more {@code NaN} parts and if either {@code this} or {@code factor}
506      * has one or more infinite parts (same result is returned regardless of
507      * the sign of the components).
508      * </p><p>
509      * Returns finite values in components of the result per the definitional
510      * formula in all remaining cases.</p>
511      *
512      * @param  factor value to be multiplied by this {@code Complex}.
513      * @return {@code this * factor}.
514      * @throws NullArgumentException if {@code factor} is {@code null}.
515      */
multiply(Complex factor)516     public Complex multiply(Complex factor)
517         throws NullArgumentException {
518         MathUtils.checkNotNull(factor);
519         if (isNaN || factor.isNaN) {
520             return NaN;
521         }
522         if (Double.isInfinite(real) ||
523             Double.isInfinite(imaginary) ||
524             Double.isInfinite(factor.real) ||
525             Double.isInfinite(factor.imaginary)) {
526             // we don't use isInfinite() to avoid testing for NaN again
527             return INF;
528         }
529         return createComplex(real * factor.real - imaginary * factor.imaginary,
530                              real * factor.imaginary + imaginary * factor.real);
531     }
532 
533     /**
534      * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
535      * interpreted as a integer number.
536      *
537      * @param  factor value to be multiplied by this {@code Complex}.
538      * @return {@code this * factor}.
539      * @see #multiply(Complex)
540      */
multiply(final int factor)541     public Complex multiply(final int factor) {
542         if (isNaN) {
543             return NaN;
544         }
545         if (Double.isInfinite(real) ||
546             Double.isInfinite(imaginary)) {
547             return INF;
548         }
549         return createComplex(real * factor, imaginary * factor);
550     }
551 
552     /**
553      * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
554      * interpreted as a real number.
555      *
556      * @param  factor value to be multiplied by this {@code Complex}.
557      * @return {@code this * factor}.
558      * @see #multiply(Complex)
559      */
multiply(double factor)560     public Complex multiply(double factor) {
561         if (isNaN || Double.isNaN(factor)) {
562             return NaN;
563         }
564         if (Double.isInfinite(real) ||
565             Double.isInfinite(imaginary) ||
566             Double.isInfinite(factor)) {
567             // we don't use isInfinite() to avoid testing for NaN again
568             return INF;
569         }
570         return createComplex(real * factor, imaginary * factor);
571     }
572 
573     /**
574      * Returns a {@code Complex} whose value is {@code (-this)}.
575      * Returns {@code NaN} if either real or imaginary
576      * part of this Complex number is {@code Double.NaN}.
577      *
578      * @return {@code -this}.
579      */
negate()580     public Complex negate() {
581         if (isNaN) {
582             return NaN;
583         }
584 
585         return createComplex(-real, -imaginary);
586     }
587 
588     /**
589      * Returns a {@code Complex} whose value is
590      * {@code (this - subtrahend)}.
591      * Uses the definitional formula
592      * <p>
593      *  {@code (a + bi) - (c + di) = (a-c) + (b-d)i}
594      * </p>
595      * If either {@code this} or {@code subtrahend} has a {@code NaN]} value in either part,
596      * {@link #NaN} is returned; otherwise infinite and {@code NaN} values are
597      * returned in the parts of the result according to the rules for
598      * {@link java.lang.Double} arithmetic.
599      *
600      * @param  subtrahend value to be subtracted from this {@code Complex}.
601      * @return {@code this - subtrahend}.
602      * @throws NullArgumentException if {@code subtrahend} is {@code null}.
603      */
subtract(Complex subtrahend)604     public Complex subtract(Complex subtrahend)
605         throws NullArgumentException {
606         MathUtils.checkNotNull(subtrahend);
607         if (isNaN || subtrahend.isNaN) {
608             return NaN;
609         }
610 
611         return createComplex(real - subtrahend.getReal(),
612                              imaginary - subtrahend.getImaginary());
613     }
614 
615     /**
616      * Returns a {@code Complex} whose value is
617      * {@code (this - subtrahend)}.
618      *
619      * @param  subtrahend value to be subtracted from this {@code Complex}.
620      * @return {@code this - subtrahend}.
621      * @see #subtract(Complex)
622      */
subtract(double subtrahend)623     public Complex subtract(double subtrahend) {
624         if (isNaN || Double.isNaN(subtrahend)) {
625             return NaN;
626         }
627         return createComplex(real - subtrahend, imaginary);
628     }
629 
630     /**
631      * Compute the
632      * <a href="http://mathworld.wolfram.com/InverseCosine.html" TARGET="_top">
633      * inverse cosine</a> of this complex number.
634      * Implements the formula:
635      * <p>
636      *  {@code acos(z) = -i (log(z + i (sqrt(1 - z<sup>2</sup>))))}
637      * </p>
638      * Returns {@link Complex#NaN} if either real or imaginary part of the
639      * input argument is {@code NaN} or infinite.
640      *
641      * @return the inverse cosine of this complex number.
642      * @since 1.2
643      */
acos()644     public Complex acos() {
645         if (isNaN) {
646             return NaN;
647         }
648 
649         return this.add(this.sqrt1z().multiply(I)).log().multiply(I.negate());
650     }
651 
652     /**
653      * Compute the
654      * <a href="http://mathworld.wolfram.com/InverseSine.html" TARGET="_top">
655      * inverse sine</a> of this complex number.
656      * Implements the formula:
657      * <p>
658      *  {@code asin(z) = -i (log(sqrt(1 - z<sup>2</sup>) + iz))}
659      * </p><p>
660      * Returns {@link Complex#NaN} if either real or imaginary part of the
661      * input argument is {@code NaN} or infinite.</p>
662      *
663      * @return the inverse sine of this complex number.
664      * @since 1.2
665      */
asin()666     public Complex asin() {
667         if (isNaN) {
668             return NaN;
669         }
670 
671         return sqrt1z().add(this.multiply(I)).log().multiply(I.negate());
672     }
673 
674     /**
675      * Compute the
676      * <a href="http://mathworld.wolfram.com/InverseTangent.html" TARGET="_top">
677      * inverse tangent</a> of this complex number.
678      * Implements the formula:
679      * <p>
680      * {@code atan(z) = (i/2) log((i + z)/(i - z))}
681      * </p><p>
682      * Returns {@link Complex#NaN} if either real or imaginary part of the
683      * input argument is {@code NaN} or infinite.</p>
684      *
685      * @return the inverse tangent of this complex number
686      * @since 1.2
687      */
atan()688     public Complex atan() {
689         if (isNaN) {
690             return NaN;
691         }
692 
693         return this.add(I).divide(I.subtract(this)).log()
694                 .multiply(I.divide(createComplex(2.0, 0.0)));
695     }
696 
697     /**
698      * Compute the
699      * <a href="http://mathworld.wolfram.com/Cosine.html" TARGET="_top">
700      * cosine</a> of this complex number.
701      * Implements the formula:
702      * <p>
703      *  {@code cos(a + bi) = cos(a)cosh(b) - sin(a)sinh(b)i}
704      * </p><p>
705      * where the (real) functions on the right-hand side are
706      * {@link FastMath#sin}, {@link FastMath#cos},
707      * {@link FastMath#cosh} and {@link FastMath#sinh}.
708      * </p><p>
709      * Returns {@link Complex#NaN} if either real or imaginary part of the
710      * input argument is {@code NaN}.
711      * </p><p>
712      * Infinite values in real or imaginary parts of the input may result in
713      * infinite or NaN values returned in parts of the result.</p>
714      * <pre>
715      *  Examples:
716      *  <code>
717      *   cos(1 &plusmn; INFINITY i) = 1 \u2213 INFINITY i
718      *   cos(&plusmn;INFINITY + i) = NaN + NaN i
719      *   cos(&plusmn;INFINITY &plusmn; INFINITY i) = NaN + NaN i
720      *  </code>
721      * </pre>
722      *
723      * @return the cosine of this complex number.
724      * @since 1.2
725      */
cos()726     public Complex cos() {
727         if (isNaN) {
728             return NaN;
729         }
730 
731         return createComplex(FastMath.cos(real) * FastMath.cosh(imaginary),
732                              -FastMath.sin(real) * FastMath.sinh(imaginary));
733     }
734 
735     /**
736      * Compute the
737      * <a href="http://mathworld.wolfram.com/HyperbolicCosine.html" TARGET="_top">
738      * hyperbolic cosine</a> of this complex number.
739      * Implements the formula:
740      * <pre>
741      *  <code>
742      *   cosh(a + bi) = cosh(a)cos(b) + sinh(a)sin(b)i
743      *  </code>
744      * </pre>
745      * where the (real) functions on the right-hand side are
746      * {@link FastMath#sin}, {@link FastMath#cos},
747      * {@link FastMath#cosh} and {@link FastMath#sinh}.
748      * <p>
749      * Returns {@link Complex#NaN} if either real or imaginary part of the
750      * input argument is {@code NaN}.
751      * </p>
752      * Infinite values in real or imaginary parts of the input may result in
753      * infinite or NaN values returned in parts of the result.
754      * <pre>
755      *  Examples:
756      *  <code>
757      *   cosh(1 &plusmn; INFINITY i) = NaN + NaN i
758      *   cosh(&plusmn;INFINITY + i) = INFINITY &plusmn; INFINITY i
759      *   cosh(&plusmn;INFINITY &plusmn; INFINITY i) = NaN + NaN i
760      *  </code>
761      * </pre>
762      *
763      * @return the hyperbolic cosine of this complex number.
764      * @since 1.2
765      */
cosh()766     public Complex cosh() {
767         if (isNaN) {
768             return NaN;
769         }
770 
771         return createComplex(FastMath.cosh(real) * FastMath.cos(imaginary),
772                              FastMath.sinh(real) * FastMath.sin(imaginary));
773     }
774 
775     /**
776      * Compute the
777      * <a href="http://mathworld.wolfram.com/ExponentialFunction.html" TARGET="_top">
778      * exponential function</a> of this complex number.
779      * Implements the formula:
780      * <pre>
781      *  <code>
782      *   exp(a + bi) = exp(a)cos(b) + exp(a)sin(b)i
783      *  </code>
784      * </pre>
785      * where the (real) functions on the right-hand side are
786      * {@link FastMath#exp}, {@link FastMath#cos}, and
787      * {@link FastMath#sin}.
788      * <p>
789      * Returns {@link Complex#NaN} if either real or imaginary part of the
790      * input argument is {@code NaN}.
791      * </p>
792      * Infinite values in real or imaginary parts of the input may result in
793      * infinite or NaN values returned in parts of the result.
794      * <pre>
795      *  Examples:
796      *  <code>
797      *   exp(1 &plusmn; INFINITY i) = NaN + NaN i
798      *   exp(INFINITY + i) = INFINITY + INFINITY i
799      *   exp(-INFINITY + i) = 0 + 0i
800      *   exp(&plusmn;INFINITY &plusmn; INFINITY i) = NaN + NaN i
801      *  </code>
802      * </pre>
803      *
804      * @return <code><i>e</i><sup>this</sup></code>.
805      * @since 1.2
806      */
exp()807     public Complex exp() {
808         if (isNaN) {
809             return NaN;
810         }
811 
812         double expReal = FastMath.exp(real);
813         return createComplex(expReal *  FastMath.cos(imaginary),
814                              expReal * FastMath.sin(imaginary));
815     }
816 
817     /**
818      * Compute the
819      * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html" TARGET="_top">
820      * natural logarithm</a> of this complex number.
821      * Implements the formula:
822      * <pre>
823      *  <code>
824      *   log(a + bi) = ln(|a + bi|) + arg(a + bi)i
825      *  </code>
826      * </pre>
827      * where ln on the right hand side is {@link FastMath#log},
828      * {@code |a + bi|} is the modulus, {@link Complex#abs},  and
829      * {@code arg(a + bi) = }{@link FastMath#atan2}(b, a).
830      * <p>
831      * Returns {@link Complex#NaN} if either real or imaginary part of the
832      * input argument is {@code NaN}.
833      * </p>
834      * Infinite (or critical) values in real or imaginary parts of the input may
835      * result in infinite or NaN values returned in parts of the result.
836      * <pre>
837      *  Examples:
838      *  <code>
839      *   log(1 &plusmn; INFINITY i) = INFINITY &plusmn; (&pi;/2)i
840      *   log(INFINITY + i) = INFINITY + 0i
841      *   log(-INFINITY + i) = INFINITY + &pi;i
842      *   log(INFINITY &plusmn; INFINITY i) = INFINITY &plusmn; (&pi;/4)i
843      *   log(-INFINITY &plusmn; INFINITY i) = INFINITY &plusmn; (3&pi;/4)i
844      *   log(0 + 0i) = -INFINITY + 0i
845      *  </code>
846      * </pre>
847      *
848      * @return the value <code>ln &nbsp; this</code>, the natural logarithm
849      * of {@code this}.
850      * @since 1.2
851      */
log()852     public Complex log() {
853         if (isNaN) {
854             return NaN;
855         }
856 
857         return createComplex(FastMath.log(abs()),
858                              FastMath.atan2(imaginary, real));
859     }
860 
861     /**
862      * Returns of value of this complex number raised to the power of {@code x}.
863      * Implements the formula:
864      * <pre>
865      *  <code>
866      *   y<sup>x</sup> = exp(x&middot;log(y))
867      *  </code>
868      * </pre>
869      * where {@code exp} and {@code log} are {@link #exp} and
870      * {@link #log}, respectively.
871      * <p>
872      * Returns {@link Complex#NaN} if either real or imaginary part of the
873      * input argument is {@code NaN} or infinite, or if {@code y}
874      * equals {@link Complex#ZERO}.</p>
875      *
876      * @param  x exponent to which this {@code Complex} is to be raised.
877      * @return <code> this<sup>x</sup></code>.
878      * @throws NullArgumentException if x is {@code null}.
879      * @since 1.2
880      */
pow(Complex x)881     public Complex pow(Complex x)
882         throws NullArgumentException {
883         MathUtils.checkNotNull(x);
884         return this.log().multiply(x).exp();
885     }
886 
887     /**
888      * Returns of value of this complex number raised to the power of {@code x}.
889      *
890      * @param  x exponent to which this {@code Complex} is to be raised.
891      * @return <code>this<sup>x</sup></code>.
892      * @see #pow(Complex)
893      */
pow(double x)894      public Complex pow(double x) {
895         return this.log().multiply(x).exp();
896     }
897 
898     /**
899      * Compute the
900      * <a href="http://mathworld.wolfram.com/Sine.html" TARGET="_top">
901      * sine</a>
902      * of this complex number.
903      * Implements the formula:
904      * <pre>
905      *  <code>
906      *   sin(a + bi) = sin(a)cosh(b) - cos(a)sinh(b)i
907      *  </code>
908      * </pre>
909      * where the (real) functions on the right-hand side are
910      * {@link FastMath#sin}, {@link FastMath#cos},
911      * {@link FastMath#cosh} and {@link FastMath#sinh}.
912      * <p>
913      * Returns {@link Complex#NaN} if either real or imaginary part of the
914      * input argument is {@code NaN}.
915      * </p><p>
916      * Infinite values in real or imaginary parts of the input may result in
917      * infinite or {@code NaN} values returned in parts of the result.
918      * <pre>
919      *  Examples:
920      *  <code>
921      *   sin(1 &plusmn; INFINITY i) = 1 &plusmn; INFINITY i
922      *   sin(&plusmn;INFINITY + i) = NaN + NaN i
923      *   sin(&plusmn;INFINITY &plusmn; INFINITY i) = NaN + NaN i
924      *  </code>
925      * </pre>
926      *
927      * @return the sine of this complex number.
928      * @since 1.2
929      */
sin()930     public Complex sin() {
931         if (isNaN) {
932             return NaN;
933         }
934 
935         return createComplex(FastMath.sin(real) * FastMath.cosh(imaginary),
936                              FastMath.cos(real) * FastMath.sinh(imaginary));
937     }
938 
939     /**
940      * Compute the
941      * <a href="http://mathworld.wolfram.com/HyperbolicSine.html" TARGET="_top">
942      * hyperbolic sine</a> of this complex number.
943      * Implements the formula:
944      * <pre>
945      *  <code>
946      *   sinh(a + bi) = sinh(a)cos(b)) + cosh(a)sin(b)i
947      *  </code>
948      * </pre>
949      * where the (real) functions on the right-hand side are
950      * {@link FastMath#sin}, {@link FastMath#cos},
951      * {@link FastMath#cosh} and {@link FastMath#sinh}.
952      * <p>
953      * Returns {@link Complex#NaN} if either real or imaginary part of the
954      * input argument is {@code NaN}.
955      * </p><p>
956      * Infinite values in real or imaginary parts of the input may result in
957      * infinite or NaN values returned in parts of the result.
958      * <pre>
959      *  Examples:
960      *  <code>
961      *   sinh(1 &plusmn; INFINITY i) = NaN + NaN i
962      *   sinh(&plusmn;INFINITY + i) = &plusmn; INFINITY + INFINITY i
963      *   sinh(&plusmn;INFINITY &plusmn; INFINITY i) = NaN + NaN i
964      *  </code>
965      * </pre>
966      *
967      * @return the hyperbolic sine of {@code this}.
968      * @since 1.2
969      */
sinh()970     public Complex sinh() {
971         if (isNaN) {
972             return NaN;
973         }
974 
975         return createComplex(FastMath.sinh(real) * FastMath.cos(imaginary),
976             FastMath.cosh(real) * FastMath.sin(imaginary));
977     }
978 
979     /**
980      * Compute the
981      * <a href="http://mathworld.wolfram.com/SquareRoot.html" TARGET="_top">
982      * square root</a> of this complex number.
983      * Implements the following algorithm to compute {@code sqrt(a + bi)}:
984      * <ol><li>Let {@code t = sqrt((|a| + |a + bi|) / 2)}</li>
985      * <li><pre>if {@code  a &#8805; 0} return {@code t + (b/2t)i}
986      *  else return {@code |b|/2t + sign(b)t i }</pre></li>
987      * </ol>
988      * where <ul>
989      * <li>{@code |a| = }{@link FastMath#abs}(a)</li>
990      * <li>{@code |a + bi| = }{@link Complex#abs}(a + bi)</li>
991      * <li>{@code sign(b) =  }{@link FastMath#copySign(double,double) copySign(1d, b)}
992      * </ul>
993      * <p>
994      * Returns {@link Complex#NaN} if either real or imaginary part of the
995      * input argument is {@code NaN}.
996      * </p>
997      * Infinite values in real or imaginary parts of the input may result in
998      * infinite or NaN values returned in parts of the result.
999      * <pre>
1000      *  Examples:
1001      *  <code>
1002      *   sqrt(1 &plusmn; INFINITY i) = INFINITY + NaN i
1003      *   sqrt(INFINITY + i) = INFINITY + 0i
1004      *   sqrt(-INFINITY + i) = 0 + INFINITY i
1005      *   sqrt(INFINITY &plusmn; INFINITY i) = INFINITY + NaN i
1006      *   sqrt(-INFINITY &plusmn; INFINITY i) = NaN &plusmn; INFINITY i
1007      *  </code>
1008      * </pre>
1009      *
1010      * @return the square root of {@code this}.
1011      * @since 1.2
1012      */
sqrt()1013     public Complex sqrt() {
1014         if (isNaN) {
1015             return NaN;
1016         }
1017 
1018         if (real == 0.0 && imaginary == 0.0) {
1019             return createComplex(0.0, 0.0);
1020         }
1021 
1022         double t = FastMath.sqrt((FastMath.abs(real) + abs()) / 2.0);
1023         if (real >= 0.0) {
1024             return createComplex(t, imaginary / (2.0 * t));
1025         } else {
1026             return createComplex(FastMath.abs(imaginary) / (2.0 * t),
1027                                  FastMath.copySign(1d, imaginary) * t);
1028         }
1029     }
1030 
1031     /**
1032      * Compute the
1033      * <a href="http://mathworld.wolfram.com/SquareRoot.html" TARGET="_top">
1034      * square root</a> of <code>1 - this<sup>2</sup></code> for this complex
1035      * number.
1036      * Computes the result directly as
1037      * {@code sqrt(ONE.subtract(z.multiply(z)))}.
1038      * <p>
1039      * Returns {@link Complex#NaN} if either real or imaginary part of the
1040      * input argument is {@code NaN}.
1041      * </p>
1042      * Infinite values in real or imaginary parts of the input may result in
1043      * infinite or NaN values returned in parts of the result.
1044      *
1045      * @return the square root of <code>1 - this<sup>2</sup></code>.
1046      * @since 1.2
1047      */
sqrt1z()1048     public Complex sqrt1z() {
1049         return createComplex(1.0, 0.0).subtract(this.multiply(this)).sqrt();
1050     }
1051 
1052     /**
1053      * Compute the
1054      * <a href="http://mathworld.wolfram.com/Tangent.html" TARGET="_top">
1055      * tangent</a> of this complex number.
1056      * Implements the formula:
1057      * <pre>
1058      *  <code>
1059      *   tan(a + bi) = sin(2a)/(cos(2a)+cosh(2b)) + [sinh(2b)/(cos(2a)+cosh(2b))]i
1060      *  </code>
1061      * </pre>
1062      * where the (real) functions on the right-hand side are
1063      * {@link FastMath#sin}, {@link FastMath#cos}, {@link FastMath#cosh} and
1064      * {@link FastMath#sinh}.
1065      * <p>
1066      * Returns {@link Complex#NaN} if either real or imaginary part of the
1067      * input argument is {@code NaN}.
1068      * </p>
1069      * Infinite (or critical) values in real or imaginary parts of the input may
1070      * result in infinite or NaN values returned in parts of the result.
1071      * <pre>
1072      *  Examples:
1073      *  <code>
1074      *   tan(a &plusmn; INFINITY i) = 0 &plusmn; i
1075      *   tan(&plusmn;INFINITY + bi) = NaN + NaN i
1076      *   tan(&plusmn;INFINITY &plusmn; INFINITY i) = NaN + NaN i
1077      *   tan(&plusmn;&pi;/2 + 0 i) = &plusmn;INFINITY + NaN i
1078      *  </code>
1079      * </pre>
1080      *
1081      * @return the tangent of {@code this}.
1082      * @since 1.2
1083      */
tan()1084     public Complex tan() {
1085         if (isNaN || Double.isInfinite(real)) {
1086             return NaN;
1087         }
1088         if (imaginary > 20.0) {
1089             return createComplex(0.0, 1.0);
1090         }
1091         if (imaginary < -20.0) {
1092             return createComplex(0.0, -1.0);
1093         }
1094 
1095         double real2 = 2.0 * real;
1096         double imaginary2 = 2.0 * imaginary;
1097         double d = FastMath.cos(real2) + FastMath.cosh(imaginary2);
1098 
1099         return createComplex(FastMath.sin(real2) / d,
1100                              FastMath.sinh(imaginary2) / d);
1101     }
1102 
1103     /**
1104      * Compute the
1105      * <a href="http://mathworld.wolfram.com/HyperbolicTangent.html" TARGET="_top">
1106      * hyperbolic tangent</a> of this complex number.
1107      * Implements the formula:
1108      * <pre>
1109      *  <code>
1110      *   tan(a + bi) = sinh(2a)/(cosh(2a)+cos(2b)) + [sin(2b)/(cosh(2a)+cos(2b))]i
1111      *  </code>
1112      * </pre>
1113      * where the (real) functions on the right-hand side are
1114      * {@link FastMath#sin}, {@link FastMath#cos}, {@link FastMath#cosh} and
1115      * {@link FastMath#sinh}.
1116      * <p>
1117      * Returns {@link Complex#NaN} if either real or imaginary part of the
1118      * input argument is {@code NaN}.
1119      * </p>
1120      * Infinite values in real or imaginary parts of the input may result in
1121      * infinite or NaN values returned in parts of the result.
1122      * <pre>
1123      *  Examples:
1124      *  <code>
1125      *   tanh(a &plusmn; INFINITY i) = NaN + NaN i
1126      *   tanh(&plusmn;INFINITY + bi) = &plusmn;1 + 0 i
1127      *   tanh(&plusmn;INFINITY &plusmn; INFINITY i) = NaN + NaN i
1128      *   tanh(0 + (&pi;/2)i) = NaN + INFINITY i
1129      *  </code>
1130      * </pre>
1131      *
1132      * @return the hyperbolic tangent of {@code this}.
1133      * @since 1.2
1134      */
tanh()1135     public Complex tanh() {
1136         if (isNaN || Double.isInfinite(imaginary)) {
1137             return NaN;
1138         }
1139         if (real > 20.0) {
1140             return createComplex(1.0, 0.0);
1141         }
1142         if (real < -20.0) {
1143             return createComplex(-1.0, 0.0);
1144         }
1145         double real2 = 2.0 * real;
1146         double imaginary2 = 2.0 * imaginary;
1147         double d = FastMath.cosh(real2) + FastMath.cos(imaginary2);
1148 
1149         return createComplex(FastMath.sinh(real2) / d,
1150                              FastMath.sin(imaginary2) / d);
1151     }
1152 
1153 
1154 
1155     /**
1156      * Compute the argument of this complex number.
1157      * The argument is the angle phi between the positive real axis and
1158      * the point representing this number in the complex plane.
1159      * The value returned is between -PI (not inclusive)
1160      * and PI (inclusive), with negative values returned for numbers with
1161      * negative imaginary parts.
1162      * <p>
1163      * If either real or imaginary part (or both) is NaN, NaN is returned.
1164      * Infinite parts are handled as {@code Math.atan2} handles them,
1165      * essentially treating finite parts as zero in the presence of an
1166      * infinite coordinate and returning a multiple of pi/4 depending on
1167      * the signs of the infinite parts.
1168      * See the javadoc for {@code Math.atan2} for full details.
1169      *
1170      * @return the argument of {@code this}.
1171      */
getArgument()1172     public double getArgument() {
1173         return FastMath.atan2(getImaginary(), getReal());
1174     }
1175 
1176     /**
1177      * Computes the n-th roots of this complex number.
1178      * The nth roots are defined by the formula:
1179      * <pre>
1180      *  <code>
1181      *   z<sub>k</sub> = abs<sup>1/n</sup> (cos(phi + 2&pi;k/n) + i (sin(phi + 2&pi;k/n))
1182      *  </code>
1183      * </pre>
1184      * for <i>{@code k=0, 1, ..., n-1}</i>, where {@code abs} and {@code phi}
1185      * are respectively the {@link #abs() modulus} and
1186      * {@link #getArgument() argument} of this complex number.
1187      * <p>
1188      * If one or both parts of this complex number is NaN, a list with just
1189      * one element, {@link #NaN} is returned.
1190      * if neither part is NaN, but at least one part is infinite, the result
1191      * is a one-element list containing {@link #INF}.
1192      *
1193      * @param n Degree of root.
1194      * @return a List of all {@code n}-th roots of {@code this}.
1195      * @throws NotPositiveException if {@code n <= 0}.
1196      * @since 2.0
1197      */
nthRoot(int n)1198     public List<Complex> nthRoot(int n) throws NotPositiveException {
1199 
1200         if (n <= 0) {
1201             throw new NotPositiveException(LocalizedFormats.CANNOT_COMPUTE_NTH_ROOT_FOR_NEGATIVE_N,
1202                                            n);
1203         }
1204 
1205         final List<Complex> result = new ArrayList<Complex>();
1206 
1207         if (isNaN) {
1208             result.add(NaN);
1209             return result;
1210         }
1211         if (isInfinite()) {
1212             result.add(INF);
1213             return result;
1214         }
1215 
1216         // nth root of abs -- faster / more accurate to use a solver here?
1217         final double nthRootOfAbs = FastMath.pow(abs(), 1.0 / n);
1218 
1219         // Compute nth roots of complex number with k = 0, 1, ... n-1
1220         final double nthPhi = getArgument() / n;
1221         final double slice = 2 * FastMath.PI / n;
1222         double innerPart = nthPhi;
1223         for (int k = 0; k < n ; k++) {
1224             // inner part
1225             final double realPart = nthRootOfAbs *  FastMath.cos(innerPart);
1226             final double imaginaryPart = nthRootOfAbs *  FastMath.sin(innerPart);
1227             result.add(createComplex(realPart, imaginaryPart));
1228             innerPart += slice;
1229         }
1230 
1231         return result;
1232     }
1233 
1234     /**
1235      * Create a complex number given the real and imaginary parts.
1236      *
1237      * @param realPart Real part.
1238      * @param imaginaryPart Imaginary part.
1239      * @return a new complex number instance.
1240      * @since 1.2
1241      * @see #valueOf(double, double)
1242      */
createComplex(double realPart, double imaginaryPart)1243     protected Complex createComplex(double realPart,
1244                                     double imaginaryPart) {
1245         return new Complex(realPart, imaginaryPart);
1246     }
1247 
1248     /**
1249      * Create a complex number given the real and imaginary parts.
1250      *
1251      * @param realPart Real part.
1252      * @param imaginaryPart Imaginary part.
1253      * @return a Complex instance.
1254      */
valueOf(double realPart, double imaginaryPart)1255     public static Complex valueOf(double realPart,
1256                                   double imaginaryPart) {
1257         if (Double.isNaN(realPart) ||
1258             Double.isNaN(imaginaryPart)) {
1259             return NaN;
1260         }
1261         return new Complex(realPart, imaginaryPart);
1262     }
1263 
1264     /**
1265      * Create a complex number given only the real part.
1266      *
1267      * @param realPart Real part.
1268      * @return a Complex instance.
1269      */
valueOf(double realPart)1270     public static Complex valueOf(double realPart) {
1271         if (Double.isNaN(realPart)) {
1272             return NaN;
1273         }
1274         return new Complex(realPart);
1275     }
1276 
1277     /**
1278      * Resolve the transient fields in a deserialized Complex Object.
1279      * Subclasses will need to override {@link #createComplex} to
1280      * deserialize properly.
1281      *
1282      * @return A Complex instance with all fields resolved.
1283      * @since 2.0
1284      */
readResolve()1285     protected final Object readResolve() {
1286         return createComplex(real, imaginary);
1287     }
1288 
1289     /** {@inheritDoc} */
getField()1290     public ComplexField getField() {
1291         return ComplexField.getInstance();
1292     }
1293 
1294     /** {@inheritDoc} */
1295     @Override
toString()1296     public String toString() {
1297         return "(" + real + ", " + imaginary + ")";
1298     }
1299 
1300 }
1301