001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.numbers.complex;
019
020import java.io.Serializable;
021import java.util.ArrayList;
022import java.util.List;
023import java.util.function.DoubleUnaryOperator;
024
025/**
026 * Cartesian representation of a complex number. The complex number is expressed
027 * in the form \( a + ib \) where \( a \) and \( b \) are real numbers and \( i \)
028 * is the imaginary unit which satisfies the equation \( i^2 = -1 \). For the
029 * complex number \( a + ib \), \( a \) is called the <em>real part</em> and
030 * \( b \) is called the <em>imaginary part</em>.
031 *
032 * <p>This class is immutable. All arithmetic will create a new instance for the
033 * result.</p>
034 *
035 * <p>Arithmetic in this class conforms to the C99 standard for complex numbers
036 * defined in ISO/IEC 9899, Annex G. Methods have been named using the equivalent
037 * method in ISO C99. The behavior for special cases is listed as defined in C99.</p>
038 *
039 * <p>For functions \( f \) which obey the conjugate equality \( conj(f(z)) = f(conj(z)) \),
040 * the specifications for the upper half-plane imply the specifications for the lower
041 * half-plane.</p>
042 *
043 * <p>For functions that are either odd, \( f(z) = -f(-z) \), or even, \( f(z) =  f(-z) \),
044 * the specifications for the first quadrant imply the specifications for the other three
045 * quadrants.</p>
046 *
047 * <p>Special cases of <a href="http://mathworld.wolfram.com/BranchCut.html">branch cuts</a>
048 * for multivalued functions adopt the principle value convention from C99. Specials cases
049 * from C99 that raise the "invalid" or "divide-by-zero"
050 * <a href="https://en.cppreference.com/w/c/numeric/fenv/FE_exceptions">floating-point
051 * exceptions</a> return the documented value without an explicit mechanism to notify
052 * of the exception case, that is no exceptions are thrown during computations in-line with
053 * the convention of the corresponding single-valued functions in
054 * {@link Math Math}.
055 * These cases are documented in the method special cases as "invalid" or "divide-by-zero"
056 * floating-point operation.
057 * Note: Invalid floating-point exception cases will result in a complex number where the
058 * cardinality of NaN component parts has increased as a real or imaginary part could
059 * not be computed and is set to NaN.
060 *
061 * @see <a href="http://www.open-std.org/JTC1/SC22/WG14/www/standards">
062 *    ISO/IEC 9899 - Programming languages - C</a>
063 */
064public final class Complex implements Serializable  {
065    /**
066     * A complex number representing \( i \), the square root of \( -1 \).
067     *
068     * <p>\( (0 + i 1) \).
069     */
070    public static final Complex I = new Complex(0, 1);
071    /**
072     * A complex number representing one.
073     *
074     * <p>\( (1 + i 0) \).
075     */
076    public static final Complex ONE = new Complex(1, 0);
077    /**
078     * A complex number representing zero.
079     *
080     * <p>\( (0 + i 0) \).
081     */
082    public static final Complex ZERO = new Complex(0, 0);
083
084    /** A complex number representing {@code NaN + i NaN}. */
085    private static final Complex NAN = new Complex(Double.NaN, Double.NaN);
086    /** &pi;/2. */
087    private static final double PI_OVER_2 = 0.5 * Math.PI;
088    /** &pi;/4. */
089    private static final double PI_OVER_4 = 0.25 * Math.PI;
090    /** Natural logarithm of 2 (ln(2)). */
091    private static final double LN_2 = Math.log(2);
092    /** Base 10 logarithm of 10 divided by 2 (log10(e)/2). */
093    private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2;
094    /** Base 10 logarithm of 2 (log10(2)). */
095    private static final double LOG10_2 = Math.log10(2);
096    /** {@code 1/2}. */
097    private static final double HALF = 0.5;
098    /** {@code sqrt(2)}. */
099    private static final double ROOT2 = 1.4142135623730951;
100    /** {@code 1.0 / sqrt(2)}.
101     * This is pre-computed to the closest double from the exact result.
102     * It is 1 ULP different from 1.0 / Math.sqrt(2) but equal to Math.sqrt(2) / 2.
103     */
104    private static final double ONE_OVER_ROOT2 = 0.7071067811865476;
105    /** The bit representation of {@code -0.0}. */
106    private static final long NEGATIVE_ZERO_LONG_BITS = Double.doubleToLongBits(-0.0);
107    /** Exponent offset in IEEE754 representation. */
108    private static final int EXPONENT_OFFSET = 1023;
109    /**
110     * Largest double-precision floating-point number such that
111     * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
112     * bound on the relative error due to rounding real numbers to double
113     * precision floating-point numbers.
114     *
115     * <p>In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
116     * Copied from o.a.c.numbers.Precision.
117     *
118     * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
119     */
120    private static final double EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
121    /** Mask to remove the sign bit from a long. */
122    private static final long UNSIGN_MASK = 0x7fff_ffff_ffff_ffffL;
123    /** Mask to extract the 52-bit mantissa from a long representation of a double. */
124    private static final long MANTISSA_MASK = 0x000f_ffff_ffff_ffffL;
125    /** The multiplier used to split the double value into hi and low parts. This must be odd
126     * and a value of 2^s + 1 in the range {@code p/2 <= s <= p-1} where p is the number of
127     * bits of precision of the floating point number. Here {@code s = 27}.*/
128    private static final double MULTIPLIER = 1.34217729E8;
129
130    /**
131     * Crossover point to switch computation for asin/acos factor A.
132     * This has been updated from the 1.5 value used by Hull et al to 10
133     * as used in boost::math::complex.
134     * @see <a href="https://svn.boost.org/trac/boost/ticket/7290">Boost ticket 7290</a>
135     */
136    private static final double A_CROSSOVER = 10.0;
137    /** Crossover point to switch computation for asin/acos factor B. */
138    private static final double B_CROSSOVER = 0.6471;
139    /**
140     * The safe maximum double value {@code x} to avoid loss of precision in asin/acos.
141     * Equal to sqrt(M) / 8 in Hull, et al (1997) with M the largest normalised floating-point value.
142     */
143    private static final double SAFE_MAX = Math.sqrt(Double.MAX_VALUE) / 8;
144    /**
145     * The safe minimum double value {@code x} to avoid loss of precision/underflow in asin/acos.
146     * Equal to sqrt(u) * 4 in Hull, et al (1997) with u the smallest normalised floating-point value.
147     */
148    private static final double SAFE_MIN = Math.sqrt(Double.MIN_NORMAL) * 4;
149    /**
150     * The safe maximum double value {@code x} to avoid loss of precision in atanh.
151     * Equal to sqrt(M) / 2 with M the largest normalised floating-point value.
152     */
153    private static final double SAFE_UPPER = Math.sqrt(Double.MAX_VALUE) / 2;
154    /**
155     * The safe minimum double value {@code x} to avoid loss of precision/underflow in atanh.
156     * Equal to sqrt(u) * 2 with u the smallest normalised floating-point value.
157     */
158    private static final double SAFE_LOWER = Math.sqrt(Double.MIN_NORMAL) * 2;
159    /** The safe maximum double value {@code x} to avoid overflow in sqrt. */
160    private static final double SQRT_SAFE_UPPER = Double.MAX_VALUE / 8;
161    /**
162     * A safe maximum double value {@code m} where {@code e^m} is not infinite.
163     * This can be used when functions require approximations of sinh(x) or cosh(x)
164     * when x is large using exp(x):
165     * <pre>
166     * sinh(x) = (e^x - e^-x) / 2 = sign(x) * e^|x| / 2
167     * cosh(x) = (e^x + e^-x) / 2 = e^|x| / 2 </pre>
168     *
169     * <p>This value can be used to approximate e^x using a product:
170     *
171     * <pre>
172     * e^x = product_n (e^m) * e^(x-nm)
173     * n = (int) x/m
174     * e.g. e^2000 = e^m * e^m * e^(2000 - 2m) </pre>
175     *
176     * <p>The value should be below ln(max_value) ~ 709.783.
177     * The value m is set to an integer for less error when subtracting m and chosen as
178     * even (m=708) as it is used as a threshold in tanh with m/2.
179     *
180     * <p>The value is used to compute e^x multiplied by a small number avoiding
181     * overflow (sinh/cosh) or a small number divided by e^x without underflow due to
182     * infinite e^x (tanh). The following conditions are used:
183     * <pre>
184     * 0.5 * e^m * Double.MIN_VALUE * e^m * e^m = Infinity
185     * 2.0 / e^m / e^m = 0.0 </pre>
186     */
187    private static final double SAFE_EXP = 708;
188    /**
189     * The value of Math.exp(SAFE_EXP): e^708.
190     * To be used in overflow/underflow safe products of e^m to approximate e^x where {@code x > m}.
191     */
192    private static final double EXP_M = Math.exp(SAFE_EXP);
193
194    /** 54 shifted 20-bits to align with the exponent of the upper 32-bits of a double. */
195    private static final int EXP_54 = 0x36_00000;
196    /** Represents an exponent of 500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */
197    private static final int EXP_500 = 0x5f3_00000;
198    /** Represents an exponent of 1024 in unbiased form (infinite or nan)
199     * shifted 20-bits to align with the upper 32-bits of a double. */
200    private static final int EXP_1024 = 0x7ff_00000;
201    /** Represents an exponent of -500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */
202    private static final int EXP_NEG_500 = 0x20b_00000;
203    /** 2^600. */
204    private static final double TWO_POW_600 = 0x1.0p+600;
205    /** 2^-600. */
206    private static final double TWO_POW_NEG_600 = 0x1.0p-600;
207
208    /** Serializable version identifier. */
209    private static final long serialVersionUID = 20180201L;
210
211    /**
212     * The size of the buffer for {@link #toString()}.
213     *
214     * <p>The longest double will require a sign, a maximum of 17 digits, the decimal place
215     * and the exponent, e.g. for max value this is 24 chars: -1.7976931348623157e+308.
216     * Set the buffer size to twice this and round up to a power of 2 thus
217     * allowing for formatting characters. The size is 64.
218     */
219    private static final int TO_STRING_SIZE = 64;
220    /** The minimum number of characters in the format. This is 5, e.g. {@code "(0,0)"}. */
221    private static final int FORMAT_MIN_LEN = 5;
222    /** {@link #toString() String representation}. */
223    private static final char FORMAT_START = '(';
224    /** {@link #toString() String representation}. */
225    private static final char FORMAT_END = ')';
226    /** {@link #toString() String representation}. */
227    private static final char FORMAT_SEP = ',';
228    /** The minimum number of characters before the separator. This is 2, e.g. {@code "(0"}. */
229    private static final int BEFORE_SEP = 2;
230
231    /** The imaginary part. */
232    private final double imaginary;
233    /** The real part. */
234    private final double real;
235
236    /**
237     * Define a constructor for a Complex.
238     * This is used in functions that implement trigonomic identities.
239     */
240    @FunctionalInterface
241    private interface ComplexConstructor {
242        /**
243         * Create a complex number given the real and imaginary parts.
244         *
245         * @param real Real part.
246         * @param imaginary Imaginary part.
247         * @return {@code Complex} object.
248         */
249        Complex create(double real, double imaginary);
250    }
251
252    /**
253     * Private default constructor.
254     *
255     * @param real Real part.
256     * @param imaginary Imaginary part.
257     */
258    private Complex(double real, double imaginary) {
259        this.real = real;
260        this.imaginary = imaginary;
261    }
262
263    /**
264     * Create a complex number given the real and imaginary parts.
265     *
266     * @param real Real part.
267     * @param imaginary Imaginary part.
268     * @return {@code Complex} number.
269     */
270    public static Complex ofCartesian(double real, double imaginary) {
271        return new Complex(real, imaginary);
272    }
273
274    /**
275     * Creates a complex number from its polar representation using modulus {@code rho} (\( \rho \))
276     * and phase angle {@code theta} (\( \theta \)).
277     *
278     * \[ \begin{aligned}
279     *    x &amp;= \rho \cos(\theta) \\
280     *    y &amp;= \rho \sin(\theta) \end{aligned} \]
281     *
282     * <p>Requires that {@code rho} is non-negative and non-NaN and {@code theta} is finite;
283     * otherwise returns a complex with NaN real and imaginary parts. A {@code rho} value of
284     * {@code -0.0} is considered negative and an invalid modulus.
285     *
286     * <p>A non-NaN complex number constructed using this method will satisfy the following
287     * to within floating-point error when {@code theta} is in the range
288     * \( -\pi\ \lt \theta \leq \pi \):
289     *
290     * <pre>
291     *  Complex.ofPolar(rho, theta).abs() == rho
292     *  Complex.ofPolar(rho, theta).arg() == theta</pre>
293     *
294     * <p>If {@code rho} is infinite then the resulting parts may be infinite or NaN
295     * following the rules for double arithmetic, for example:</p>
296     *
297     * <ul>
298     * <li>{@code ofPolar(}\( -0.0 \){@code , }\( 0 \){@code ) = }\( \text{NaN} + i \text{NaN} \)
299     * <li>{@code ofPolar(}\( 0.0 \){@code , }\( 0 \){@code ) = }\( 0 + i 0 \)
300     * <li>{@code ofPolar(}\( 1 \){@code , }\( 0 \){@code ) = }\( 1 + i 0 \)
301     * <li>{@code ofPolar(}\( 1 \){@code , }\( \pi \){@code ) = }\( -1 + i \sin(\pi) \)
302     * <li>{@code ofPolar(}\( \infty \){@code , }\( \pi \){@code ) = }\( -\infty + i \infty \)
303     * <li>{@code ofPolar(}\( \infty \){@code , }\( 0 \){@code ) = }\( -\infty + i \text{NaN} \)
304     * <li>{@code ofPolar(}\( \infty \){@code , }\( -\frac{\pi}{4} \){@code ) = }\( \infty - i \infty \)
305     * <li>{@code ofPolar(}\( \infty \){@code , }\( 5\frac{\pi}{4} \){@code ) = }\( -\infty - i \infty \)
306     * </ul>
307     *
308     * <p>This method is the functional equivalent of the C++ method {@code std::polar}.
309     *
310     * @param rho The modulus of the complex number.
311     * @param theta The argument of the complex number.
312     * @return {@code Complex} number.
313     * @see <a href="http://mathworld.wolfram.com/PolarCoordinates.html">Polar Coordinates</a>
314     */
315    public static Complex ofPolar(double rho, double theta) {
316        // Require finite theta and non-negative, non-nan rho
317        if (!Double.isFinite(theta) || negative(rho) || Double.isNaN(rho)) {
318            return NAN;
319        }
320        final double x = rho * Math.cos(theta);
321        final double y = rho * Math.sin(theta);
322        return new Complex(x, y);
323    }
324
325    /**
326     * Create a complex cis number. This is also known as the complex exponential:
327     *
328     * \[ \text{cis}(x) = e^{ix} = \cos(x) + i \sin(x) \]
329     *
330     * @param x {@code double} to build the cis number.
331     * @return {@code Complex} cis number.
332     * @see <a href="http://mathworld.wolfram.com/Cis.html">Cis</a>
333     */
334    public static Complex ofCis(double x) {
335        return new Complex(Math.cos(x), Math.sin(x));
336    }
337
338    /**
339     * Returns a {@code Complex} instance representing the specified string {@code s}.
340     *
341     * <p>If {@code s} is {@code null}, then a {@code NullPointerException} is thrown.
342     *
343     * <p>The string must be in a format compatible with that produced by
344     * {@link #toString() Complex.toString()}.
345     * The format expects a start and end parentheses surrounding two numeric parts split
346     * by a separator. Leading and trailing spaces are allowed around each numeric part.
347     * Each numeric part is parsed using {@link Double#parseDouble(String)}. The parts
348     * are interpreted as the real and imaginary parts of the complex number.
349     *
350     * <p>Examples of valid strings and the equivalent {@code Complex} are shown below:
351     *
352     * <pre>
353     * "(0,0)"             = Complex.ofCartesian(0, 0)
354     * "(0.0,0.0)"         = Complex.ofCartesian(0, 0)
355     * "(-0.0, 0.0)"       = Complex.ofCartesian(-0.0, 0)
356     * "(-1.23, 4.56)"     = Complex.ofCartesian(-1.23, 4.56)
357     * "(1e300,-1.1e-2)"   = Complex.ofCartesian(1e300, -1.1e-2)</pre>
358     *
359     * @param s String representation.
360     * @return {@code Complex} number.
361     * @throws NullPointerException if the string is null.
362     * @throws NumberFormatException if the string does not contain a parsable complex number.
363     * @see Double#parseDouble(String)
364     * @see #toString()
365     */
366    public static Complex parse(String s) {
367        final int len = s.length();
368        if (len < FORMAT_MIN_LEN) {
369            throw new NumberFormatException(
370                parsingExceptionMsg("Input too short, expected format",
371                                    FORMAT_START + "x" + FORMAT_SEP + "y" + FORMAT_END, s));
372        }
373
374        // Confirm start: '('
375        if (s.charAt(0) != FORMAT_START) {
376            throw new NumberFormatException(
377                parsingExceptionMsg("Expected start delimiter", FORMAT_START, s));
378        }
379
380        // Confirm end: ')'
381        if (s.charAt(len - 1) != FORMAT_END) {
382            throw new NumberFormatException(
383                parsingExceptionMsg("Expected end delimiter", FORMAT_END, s));
384        }
385
386        // Confirm separator ',' is between at least 2 characters from
387        // either end: "(x,x)"
388        // Count back from the end ignoring the last 2 characters.
389        final int sep = s.lastIndexOf(FORMAT_SEP, len - 3);
390        if (sep < BEFORE_SEP) {
391            throw new NumberFormatException(
392                parsingExceptionMsg("Expected separator between two numbers", FORMAT_SEP, s));
393        }
394
395        // Should be no more separators
396        if (s.indexOf(FORMAT_SEP, sep + 1) != -1) {
397            throw new NumberFormatException(
398                parsingExceptionMsg("Incorrect number of parts, expected only 2 using separator",
399                                    FORMAT_SEP, s));
400        }
401
402        // Try to parse the parts
403
404        final String rePart = s.substring(1, sep);
405        final double re;
406        try {
407            re = Double.parseDouble(rePart);
408        } catch (final NumberFormatException ex) {
409            throw new NumberFormatException(
410                parsingExceptionMsg("Could not parse real part", rePart, s));
411        }
412
413        final String imPart = s.substring(sep + 1, len - 1);
414        final double im;
415        try {
416            im = Double.parseDouble(imPart);
417        } catch (final NumberFormatException ex) {
418            throw new NumberFormatException(
419                parsingExceptionMsg("Could not parse imaginary part", imPart, s));
420        }
421
422        return ofCartesian(re, im);
423    }
424
425    /**
426     * Creates an exception message.
427     *
428     * @param message Message prefix.
429     * @param error Input that caused the error.
430     * @param s String representation.
431     * @return A message.
432     */
433    private static String parsingExceptionMsg(String message,
434                                              Object error,
435                                              String s) {
436        final StringBuilder sb = new StringBuilder(100)
437            .append(message)
438            .append(" '").append(error)
439            .append("' for input \"").append(s).append('"');
440        return sb.toString();
441    }
442
443    /**
444     * Gets the real part \( a \) of this complex number \( (a + i b) \).
445     *
446     * @return The real part.
447     */
448    public double getReal() {
449        return real;
450    }
451
452    /**
453     * Gets the real part \( a \) of this complex number \( (a + i b) \).
454     *
455     * <p>This method is the equivalent of the C++ method {@code std::complex::real}.
456     *
457     * @return The real part.
458     * @see #getReal()
459     */
460    public double real() {
461        return getReal();
462    }
463
464    /**
465     * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
466     *
467     * @return The imaginary part.
468     */
469    public double getImaginary() {
470        return imaginary;
471    }
472
473    /**
474     * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
475     *
476     * <p>This method is the equivalent of the C++ method {@code std::complex::imag}.
477     *
478     * @return The imaginary part.
479     * @see #getImaginary()
480     */
481    public double imag() {
482        return getImaginary();
483    }
484
485    /**
486     * Returns the absolute value of this complex number. This is also called complex norm, modulus,
487     * or magnitude.
488     *
489     * <p>\[ \text{abs}(x + i y) = \sqrt{(x^2 + y^2)} \]
490     *
491     * <p>Special cases:
492     *
493     * <ul>
494     * <li>{@code abs(x + iy) == abs(y + ix) == abs(x - iy)}.
495     * <li>If {@code z} is ±∞ + iy for any y, returns +∞.
496     * <li>If {@code z} is x + iNaN for non-infinite x, returns NaN.
497     * <li>If {@code z} is x + i0, returns |x|.
498     * </ul>
499     *
500     * <p>The cases ensure that if either component is infinite then the result is positive
501     * infinity. If either component is NaN and this is not {@link #isInfinite() infinite} then
502     * the result is NaN.
503     *
504     * <p>This method follows the
505     * <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
506     * in calculating the returned value without intermediate overflow or underflow.
507     *
508     * <p>The computed result will be within 1 ulp of the exact result.
509     *
510     * @return The absolute value.
511     * @see #isInfinite()
512     * @see #isNaN()
513     * @see <a href="http://mathworld.wolfram.com/ComplexModulus.html">Complex modulus</a>
514     */
515    public double abs() {
516        return abs(real, imaginary);
517    }
518
519    /**
520     * Returns the absolute value of the complex number.
521     * <pre>abs(x + i y) = sqrt(x^2 + y^2)</pre>
522     *
523     * <p>This should satisfy the special cases of the hypot function in ISO C99 F.9.4.3:
524     * "The hypot functions compute the square root of the sum of the squares of x and y,
525     * without undue overflow or underflow."
526     *
527     * <ul>
528     * <li>hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent.
529     * <li>hypot(x, ±0) is equivalent to |x|.
530     * <li>hypot(±∞, y) returns +∞, even if y is a NaN.
531     * </ul>
532     *
533     * <p>This method is called by all methods that require the absolute value of the complex
534     * number, e.g. abs(), sqrt() and log().
535     *
536     * @param real Real part.
537     * @param imaginary Imaginary part.
538     * @return The absolute value.
539     */
540    private static double abs(double real, double imaginary) {
541        // Specialised implementation of hypot.
542        // See NUMBERS-143
543        return hypot(real, imaginary);
544    }
545
546    /**
547     * Returns the argument of this complex number.
548     *
549     * <p>The argument is the angle phi between the positive real axis and
550     * the point representing this number in the complex plane.
551     * The value returned is between \( -\pi \) (not inclusive)
552     * and \( \pi \) (inclusive), with negative values returned for numbers with
553     * negative imaginary parts.
554     *
555     * <p>If either real or imaginary part (or both) is NaN, then the result is NaN.
556     * Infinite parts are handled as {@linkplain Math#atan2} handles them,
557     * essentially treating finite parts as zero in the presence of an
558     * infinite coordinate and returning a multiple of \( \frac{\pi}{4} \) depending on
559     * the signs of the infinite parts.
560     *
561     * <p>This code follows the
562     * <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
563     * in calculating the returned value using the {@code atan2(y, x)} method for complex
564     * \( x + iy \).
565     *
566     * @return The argument of this complex number.
567     * @see Math#atan2(double, double)
568     */
569    public double arg() {
570        // Delegate
571        return Math.atan2(imaginary, real);
572    }
573
574    /**
575     * Returns the squared norm value of this complex number. This is also called the absolute
576     * square.
577     *
578     * <p>\[ \text{norm}(x + i y) = x^2 + y^2 \]
579     *
580     * <p>If either component is infinite then the result is positive infinity. If either
581     * component is NaN and this is not {@link #isInfinite() infinite} then the result is NaN.
582     *
583     * <p>Note: This method may not return the same value as the square of {@link #abs()} as
584     * that method uses an extended precision computation.
585     *
586     * <p>{@code norm()} can be used as a faster alternative than {@code abs()} for ranking by
587     * magnitude. If used for ranking any overflow to infinity will create an equal ranking for
588     * values that may be still distinguished by {@code abs()}.
589     *
590     * @return The square norm value.
591     * @see #isInfinite()
592     * @see #isNaN()
593     * @see #abs()
594     * @see <a href="http://mathworld.wolfram.com/AbsoluteSquare.html">Absolute square</a>
595     */
596    public double norm() {
597        if (isInfinite()) {
598            return Double.POSITIVE_INFINITY;
599        }
600        return real * real + imaginary * imaginary;
601    }
602
603    /**
604     * Returns {@code true} if either the real <em>or</em> imaginary component of the complex number is NaN
605     * <em>and</em> the complex number is not infinite.
606     *
607     * <p>Note that:
608     * <ul>
609     *   <li>There is more than one complex number that can return {@code true}.
610     *   <li>Different representations of NaN can be distinguished by the
611     *       {@link #equals(Object) Complex.equals(Object)} method.
612     * </ul>
613     *
614     * @return {@code true} if this instance contains NaN and no infinite parts.
615     * @see Double#isNaN(double)
616     * @see #isInfinite()
617     * @see #equals(Object) Complex.equals(Object)
618     */
619    public boolean isNaN() {
620        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
621            return !isInfinite();
622        }
623        return false;
624    }
625
626    /**
627     * Returns {@code true} if either real or imaginary component of the complex number is infinite.
628     *
629     * <p>Note: A complex number with at least one infinite part is regarded
630     * as an infinity (even if its other part is a NaN).
631     *
632     * @return {@code true} if this instance contains an infinite value.
633     * @see Double#isInfinite(double)
634     */
635    public boolean isInfinite() {
636        return Double.isInfinite(real) || Double.isInfinite(imaginary);
637    }
638
639    /**
640     * Returns {@code true} if both real and imaginary component of the complex number are finite.
641     *
642     * @return {@code true} if this instance contains finite values.
643     * @see Double#isFinite(double)
644     */
645    public boolean isFinite() {
646        return Double.isFinite(real) && Double.isFinite(imaginary);
647    }
648
649    /**
650     * Returns the
651     * <a href="http://mathworld.wolfram.com/ComplexConjugate.html">conjugate</a>
652     * \( \overline{z} \) of this complex number \( z \).
653     *
654     * <p>\[ \begin{aligned}
655     *                z  &amp;= a + i b \\
656     *      \overline{z} &amp;= a - i b \end{aligned}\]
657     *
658     * @return The conjugate (\( \overline{z} \)) of this complex number.
659     */
660    public Complex conj() {
661        return new Complex(real, -imaginary);
662    }
663
664    /**
665     * Returns a {@code Complex} whose value is the negation of both the real and imaginary parts
666     * of complex number \( z \).
667     *
668     * <p>\[ \begin{aligned}
669     *       z  &amp;=  a + i b \\
670     *      -z  &amp;= -a - i b \end{aligned} \]
671     *
672     * @return \( -z \).
673     */
674    public Complex negate() {
675        return new Complex(-real, -imaginary);
676    }
677
678    /**
679     * Returns the projection of this complex number onto the Riemann sphere.
680     *
681     * <p>\( z \) projects to \( z \), except that all complex infinities (even those
682     * with one infinite part and one NaN part) project to positive infinity on the real axis.
683     *
684     * If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to:
685     *
686     * <pre>return Complex.ofCartesian(Double.POSITIVE_INFINITY, Math.copySign(0.0, z.imag());</pre>
687     *
688     * @return \( z \) projected onto the Riemann sphere.
689     * @see #isInfinite()
690     * @see <a href="http://pubs.opengroup.org/onlinepubs/9699919799/functions/cproj.html">
691     * IEEE and ISO C standards: cproj</a>
692     */
693    public Complex proj() {
694        if (isInfinite()) {
695            return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0.0, imaginary));
696        }
697        return this;
698    }
699
700    /**
701     * Returns a {@code Complex} whose value is {@code (this + addend)}.
702     * Implements the formula:
703     *
704     * <p>\[ (a + i b) + (c + i d) = (a + c) + i (b + d) \]
705     *
706     * @param  addend Value to be added to this complex number.
707     * @return {@code this + addend}.
708     * @see <a href="http://mathworld.wolfram.com/ComplexAddition.html">Complex Addition</a>
709     */
710    public Complex add(Complex addend) {
711        return new Complex(real + addend.real,
712                           imaginary + addend.imaginary);
713    }
714
715    /**
716     * Returns a {@code Complex} whose value is {@code (this + addend)},
717     * with {@code addend} interpreted as a real number.
718     * Implements the formula:
719     *
720     * <p>\[ (a + i b) + c = (a + c) + i b \]
721     *
722     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
723     * real-only and complex numbers.</p>
724     *
725     * <p>Note: This method preserves the sign of the imaginary component \( b \) if it is {@code -0.0}.
726     * The sign would be lost if adding \( (c + i 0) \) using
727     * {@link #add(Complex) add(Complex.ofCartesian(addend, 0))} since
728     * {@code -0.0 + 0.0 = 0.0}.
729     *
730     * @param addend Value to be added to this complex number.
731     * @return {@code this + addend}.
732     * @see #add(Complex)
733     * @see #ofCartesian(double, double)
734     */
735    public Complex add(double addend) {
736        return new Complex(real + addend, imaginary);
737    }
738
739    /**
740     * Returns a {@code Complex} whose value is {@code (this + addend)},
741     * with {@code addend} interpreted as an imaginary number.
742     * Implements the formula:
743     *
744     * <p>\[ (a + i b) + i d = a + i (b + d) \]
745     *
746     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
747     * imaginary-only and complex numbers.</p>
748     *
749     * <p>Note: This method preserves the sign of the real component \( a \) if it is {@code -0.0}.
750     * The sign would be lost if adding \( (0 + i d) \) using
751     * {@link #add(Complex) add(Complex.ofCartesian(0, addend))} since
752     * {@code -0.0 + 0.0 = 0.0}.
753     *
754     * @param addend Value to be added to this complex number.
755     * @return {@code this + addend}.
756     * @see #add(Complex)
757     * @see #ofCartesian(double, double)
758     */
759    public Complex addImaginary(double addend) {
760        return new Complex(real, imaginary + addend);
761    }
762
763    /**
764     * Returns a {@code Complex} whose value is {@code (this - subtrahend)}.
765     * Implements the formula:
766     *
767     * <p>\[ (a + i b) - (c + i d) = (a - c) + i (b - d) \]
768     *
769     * @param  subtrahend Value to be subtracted from this complex number.
770     * @return {@code this - subtrahend}.
771     * @see <a href="http://mathworld.wolfram.com/ComplexSubtraction.html">Complex Subtraction</a>
772     */
773    public Complex subtract(Complex subtrahend) {
774        return new Complex(real - subtrahend.real,
775                           imaginary - subtrahend.imaginary);
776    }
777
778    /**
779     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
780     * with {@code subtrahend} interpreted as a real number.
781     * Implements the formula:
782     *
783     * <p>\[ (a + i b) - c = (a - c) + i b \]
784     *
785     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
786     * real-only and complex numbers.</p>
787     *
788     * @param  subtrahend Value to be subtracted from this complex number.
789     * @return {@code this - subtrahend}.
790     * @see #subtract(Complex)
791     */
792    public Complex subtract(double subtrahend) {
793        return new Complex(real - subtrahend, imaginary);
794    }
795
796    /**
797     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
798     * with {@code subtrahend} interpreted as an imaginary number.
799     * Implements the formula:
800     *
801     * <p>\[ (a + i b) - i d = a + i (b - d) \]
802     *
803     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
804     * imaginary-only and complex numbers.</p>
805     *
806     * @param  subtrahend Value to be subtracted from this complex number.
807     * @return {@code this - subtrahend}.
808     * @see #subtract(Complex)
809     */
810    public Complex subtractImaginary(double subtrahend) {
811        return new Complex(real, imaginary - subtrahend);
812    }
813
814    /**
815     * Returns a {@code Complex} whose value is {@code (minuend - this)},
816     * with {@code minuend} interpreted as a real number.
817     * Implements the formula:
818     * \[ c - (a + i b) = (c - a) - i b \]
819     *
820     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
821     * real-only and complex numbers.</p>
822     *
823     * <p>Note: This method inverts the sign of the imaginary component \( b \) if it is {@code 0.0}.
824     * The sign would not be inverted if subtracting from \( c + i 0 \) using
825     * {@link #subtract(Complex) Complex.ofCartesian(minuend, 0).subtract(this)} since
826     * {@code 0.0 - 0.0 = 0.0}.
827     *
828     * @param  minuend Value this complex number is to be subtracted from.
829     * @return {@code minuend - this}.
830     * @see #subtract(Complex)
831     * @see #ofCartesian(double, double)
832     */
833    public Complex subtractFrom(double minuend) {
834        return new Complex(minuend - real, -imaginary);
835    }
836
837    /**
838     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
839     * with {@code minuend} interpreted as an imaginary number.
840     * Implements the formula:
841     * \[ i d - (a + i b) = -a + i (d - b) \]
842     *
843     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
844     * imaginary-only and complex numbers.</p>
845     *
846     * <p>Note: This method inverts the sign of the real component \( a \) if it is {@code 0.0}.
847     * The sign would not be inverted if subtracting from \( 0 + i d \) using
848     * {@link #subtract(Complex) Complex.ofCartesian(0, minuend).subtract(this)} since
849     * {@code 0.0 - 0.0 = 0.0}.
850     *
851     * @param  minuend Value this complex number is to be subtracted from.
852     * @return {@code this - subtrahend}.
853     * @see #subtract(Complex)
854     * @see #ofCartesian(double, double)
855     */
856    public Complex subtractFromImaginary(double minuend) {
857        return new Complex(-real, minuend - imaginary);
858    }
859
860    /**
861     * Returns a {@code Complex} whose value is {@code this * factor}.
862     * Implements the formula:
863     *
864     * <p>\[ (a + i b)(c + i d) = (ac - bd) + i (ad + bc) \]
865     *
866     * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
867     *
868     * @param  factor Value to be multiplied by this complex number.
869     * @return {@code this * factor}.
870     * @see <a href="http://mathworld.wolfram.com/ComplexMultiplication.html">Complex Muliplication</a>
871     */
872    public Complex multiply(Complex factor) {
873        return multiply(real, imaginary, factor.real, factor.imaginary);
874    }
875
876    /**
877     * Returns a {@code Complex} whose value is:
878     * <pre>
879     *  (a + i b)(c + i d) = (ac - bd) + i (ad + bc)</pre>
880     *
881     * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
882     *
883     * @param re1 Real component of first number.
884     * @param im1 Imaginary component of first number.
885     * @param re2 Real component of second number.
886     * @param im2 Imaginary component of second number.
887     * @return (a + b i)(c + d i).
888     */
889    private static Complex multiply(double re1, double im1, double re2, double im2) {
890        double a = re1;
891        double b = im1;
892        double c = re2;
893        double d = im2;
894        final double ac = a * c;
895        final double bd = b * d;
896        final double ad = a * d;
897        final double bc = b * c;
898        double x = ac - bd;
899        double y = ad + bc;
900
901        // --------------
902        // NaN can occur if:
903        // - any of (a,b,c,d) are NaN (for NaN or Infinite complex numbers)
904        // - a multiplication of infinity by zero (ac,bd,ad,bc).
905        // - a subtraction of infinity from infinity (e.g. ac - bd)
906        //   Note that (ac,bd,ad,bc) can be infinite due to overflow.
907        //
908        // Detect a NaN result and perform correction.
909        //
910        // Modification from the listing in ISO C99 G.5.1 (6)
911        // Do not correct infinity multiplied by zero. This is left as NaN.
912        // --------------
913
914        if (Double.isNaN(x) && Double.isNaN(y)) {
915            // Recover infinities that computed as NaN+iNaN ...
916            boolean recalc = false;
917            if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
918                isNotZero(c, d)) {
919                // This complex is infinite.
920                // "Box" the infinity and change NaNs in the other factor to 0.
921                a = boxInfinity(a);
922                b = boxInfinity(b);
923                c = changeNaNtoZero(c);
924                d = changeNaNtoZero(d);
925                recalc = true;
926            }
927            if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
928                isNotZero(a, b)) {
929                // The other complex is infinite.
930                // "Box" the infinity and change NaNs in the other factor to 0.
931                c = boxInfinity(c);
932                d = boxInfinity(d);
933                a = changeNaNtoZero(a);
934                b = changeNaNtoZero(b);
935                recalc = true;
936            }
937            if (!recalc && (Double.isInfinite(ac) || Double.isInfinite(bd) ||
938                            Double.isInfinite(ad) || Double.isInfinite(bc))) {
939                // The result overflowed to infinity.
940                // Recover infinities from overflow by changing NaNs to 0 ...
941                a = changeNaNtoZero(a);
942                b = changeNaNtoZero(b);
943                c = changeNaNtoZero(c);
944                d = changeNaNtoZero(d);
945                recalc = true;
946            }
947            if (recalc) {
948                x = Double.POSITIVE_INFINITY * (a * c - b * d);
949                y = Double.POSITIVE_INFINITY * (a * d + b * c);
950            }
951        }
952        return new Complex(x, y);
953    }
954
955    /**
956     * Box values for the real or imaginary component of an infinite complex number.
957     * Any infinite value will be returned as one. Non-infinite values will be returned as zero.
958     * The sign is maintained.
959     *
960     * <pre>
961     *  inf  =  1
962     * -inf  = -1
963     *  x    =  0
964     * -x    = -0
965     * </pre>
966     *
967     * @param component the component
968     * @return The boxed value
969     */
970    private static double boxInfinity(double component) {
971        return Math.copySign(Double.isInfinite(component) ? 1.0 : 0.0, component);
972    }
973
974    /**
975     * Checks if the complex number is not zero.
976     *
977     * @param real the real component
978     * @param imaginary the imaginary component
979     * @return true if the complex is not zero
980     */
981    private static boolean isNotZero(double real, double imaginary) {
982        // The use of equals is deliberate.
983        // This method must distinguish NaN from zero thus ruling out:
984        // (real != 0.0 || imaginary != 0.0)
985        return !(real == 0.0 && imaginary == 0.0);
986    }
987
988    /**
989     * Change NaN to zero preserving the sign; otherwise return the value.
990     *
991     * @param value the value
992     * @return The new value
993     */
994    private static double changeNaNtoZero(double value) {
995        return Double.isNaN(value) ? Math.copySign(0.0, value) : value;
996    }
997
998    /**
999     * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
1000     * interpreted as a real number.
1001     * Implements the formula:
1002     *
1003     * <p>\[ (a + i b) c =  (ac) + i (bc) \]
1004     *
1005     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1006     * real-only and complex numbers.</p>
1007     *
1008     * <p>Note: This method should be preferred over using
1009     * {@link #multiply(Complex) multiply(Complex.ofCartesian(factor, 0))}. Multiplication
1010     * can generate signed zeros if either {@code this} complex has zeros for the real
1011     * and/or imaginary component, or if the factor is zero. The summation of signed zeros
1012     * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
1013     * from the equivalent call to multiply by a real-only number.
1014     *
1015     * @param  factor Value to be multiplied by this complex number.
1016     * @return {@code this * factor}.
1017     * @see #multiply(Complex)
1018     */
1019    public Complex multiply(double factor) {
1020        return new Complex(real * factor, imaginary * factor);
1021    }
1022
1023    /**
1024     * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
1025     * interpreted as an imaginary number.
1026     * Implements the formula:
1027     *
1028     * <p>\[ (a + i b) id = (-bd) + i (ad) \]
1029     *
1030     * <p>This method can be used to compute the multiplication of this complex number \( z \)
1031     * by \( i \) using a factor with magnitude 1.0. This should be used in preference to
1032     * {@link #multiply(Complex) multiply(Complex.I)} with or without {@link #negate() negation}:</p>
1033     *
1034     * \[ \begin{aligned}
1035     *    iz &amp;= (-b + i a) \\
1036     *   -iz &amp;= (b - i a) \end{aligned} \]
1037     *
1038     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1039     * imaginary-only and complex numbers.</p>
1040     *
1041     * <p>Note: This method should be preferred over using
1042     * {@link #multiply(Complex) multiply(Complex.ofCartesian(0, factor))}. Multiplication
1043     * can generate signed zeros if either {@code this} complex has zeros for the real
1044     * and/or imaginary component, or if the factor is zero. The summation of signed zeros
1045     * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
1046     * from the equivalent call to multiply by an imaginary-only number.
1047     *
1048     * @param  factor Value to be multiplied by this complex number.
1049     * @return {@code this * factor}.
1050     * @see #multiply(Complex)
1051     */
1052    public Complex multiplyImaginary(double factor) {
1053        return new Complex(-imaginary * factor, real * factor);
1054    }
1055
1056    /**
1057     * Returns a {@code Complex} whose value is {@code (this / divisor)}.
1058     * Implements the formula:
1059     *
1060     * <p>\[ \frac{a + i b}{c + i d} = \frac{(ac + bd) + i (bc - ad)}{c^2+d^2} \]
1061     *
1062     * <p>Re-calculates NaN result values to recover infinities as specified in C99 standard G.5.1.
1063     *
1064     * @param divisor Value by which this complex number is to be divided.
1065     * @return {@code this / divisor}.
1066     * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
1067     */
1068    public Complex divide(Complex divisor) {
1069        return divide(real, imaginary, divisor.real, divisor.imaginary);
1070    }
1071
1072    /**
1073     * Returns a {@code Complex} whose value is:
1074     * <pre>
1075     * <code>
1076     *   a + i b     (ac + bd) + i (bc - ad)
1077     *   -------  =  -----------------------
1078     *   c + i d            c<sup>2</sup> + d<sup>2</sup>
1079     * </code>
1080     * </pre>
1081     *
1082     * <p>Recalculates to recover infinities as specified in C99
1083     * standard G.5.1. Method is fully in accordance with
1084     * C++11 standards for complex numbers.</p>
1085     *
1086     * <p>Note: In the event of divide by zero this method produces the same result
1087     * as dividing by a real-only zero using {@link #divide(double)}.
1088     *
1089     * @param re1 Real component of first number.
1090     * @param im1 Imaginary component of first number.
1091     * @param re2 Real component of second number.
1092     * @param im2 Imaginary component of second number.
1093     * @return (a + i b) / (c + i d).
1094     * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
1095     * @see #divide(double)
1096     */
1097    private static Complex divide(double re1, double im1, double re2, double im2) {
1098        double a = re1;
1099        double b = im1;
1100        double c = re2;
1101        double d = im2;
1102        int ilogbw = 0;
1103        // Get the exponent to scale the divisor parts to the range [1, 2).
1104        final int exponent = getScale(c, d);
1105        if (exponent <= Double.MAX_EXPONENT) {
1106            ilogbw = exponent;
1107            c = Math.scalb(c, -ilogbw);
1108            d = Math.scalb(d, -ilogbw);
1109        }
1110        final double denom = c * c + d * d;
1111
1112        // Note: Modification from the listing in ISO C99 G.5.1 (8):
1113        // Avoid overflow if a or b are very big.
1114        // Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
1115        // when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
1116        // (bc - ad) with large negative values.
1117        // Use the maximum exponent as an approximation to the magnitude.
1118        if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
1119            ilogbw -= 2;
1120            a /= 4;
1121            b /= 4;
1122        }
1123
1124        double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
1125        double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
1126        // Recover infinities and zeros that computed as NaN+iNaN
1127        // the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
1128        if (Double.isNaN(x) && Double.isNaN(y)) {
1129            if (denom == 0.0 &&
1130                    (!Double.isNaN(a) || !Double.isNaN(b))) {
1131                // nonzero/zero
1132                // This case produces the same result as divide by a real-only zero
1133                // using Complex.divide(+/-0.0)
1134                x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
1135                y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
1136            } else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
1137                    Double.isFinite(c) && Double.isFinite(d)) {
1138                // infinite/finite
1139                a = boxInfinity(a);
1140                b = boxInfinity(b);
1141                x = Double.POSITIVE_INFINITY * (a * c + b * d);
1142                y = Double.POSITIVE_INFINITY * (b * c - a * d);
1143            } else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
1144                    Double.isFinite(a) && Double.isFinite(b)) {
1145                // finite/infinite
1146                c = boxInfinity(c);
1147                d = boxInfinity(d);
1148                x = 0.0 * (a * c + b * d);
1149                y = 0.0 * (b * c - a * d);
1150            }
1151        }
1152        return new Complex(x, y);
1153    }
1154
1155    /**
1156     * Returns a {@code Complex} whose value is {@code (this / divisor)},
1157     * with {@code divisor} interpreted as a real number.
1158     * Implements the formula:
1159     *
1160     * <p>\[ \frac{a + i b}{c} = \frac{a}{c} + i \frac{b}{c} \]
1161     *
1162     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1163     * real-only and complex numbers.</p>
1164     *
1165     * <p>Note: This method should be preferred over using
1166     * {@link #divide(Complex) divide(Complex.ofCartesian(divisor, 0))}. Division
1167     * can generate signed zeros if {@code this} complex has zeros for the real
1168     * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
1169     * in {@link #divide(Complex)} may create zeros in the result that differ in sign
1170     * from the equivalent call to divide by a real-only number.
1171     *
1172     * @param  divisor Value by which this complex number is to be divided.
1173     * @return {@code this / divisor}.
1174     * @see #divide(Complex)
1175     */
1176    public Complex divide(double divisor) {
1177        return new Complex(real / divisor, imaginary / divisor);
1178    }
1179
1180    /**
1181     * Returns a {@code Complex} whose value is {@code (this / divisor)},
1182     * with {@code divisor} interpreted as an imaginary number.
1183     * Implements the formula:
1184     *
1185     * <p>\[ \frac{a + i b}{id} = \frac{b}{d} - i \frac{a}{d} \]
1186     *
1187     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
1188     * imaginary-only and complex numbers.</p>
1189     *
1190     * <p>Note: This method should be preferred over using
1191     * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))}. Division
1192     * can generate signed zeros if {@code this} complex has zeros for the real
1193     * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
1194     * in {@link #divide(Complex)} may create zeros in the result that differ in sign
1195     * from the equivalent call to divide by an imaginary-only number.
1196     *
1197     * <p>Warning: This method will generate a different result from
1198     * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))} if the divisor is zero.
1199     * In this case the divide method using a zero-valued Complex will produce the same result
1200     * as dividing by a real-only zero. The output from dividing by imaginary zero will create
1201     * infinite and NaN values in the same component parts as the output from
1202     * {@code this.divide(Complex.ZERO).multiplyImaginary(1)}, however the sign
1203     * of some infinite values may be negated.
1204     *
1205     * @param  divisor Value by which this complex number is to be divided.
1206     * @return {@code this / divisor}.
1207     * @see #divide(Complex)
1208     * @see #divide(double)
1209     */
1210    public Complex divideImaginary(double divisor) {
1211        return new Complex(imaginary / divisor, -real / divisor);
1212    }
1213
1214    /**
1215     * Returns the
1216     * <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
1217     * exponential function</a> of this complex number.
1218     *
1219     * <p>\[ \exp(z) = e^z \]
1220     *
1221     * <p>The exponential function of \( z \) is an entire function in the complex plane.
1222     * Special cases:
1223     *
1224     * <ul>
1225     * <li>{@code z.conj().exp() == z.exp().conj()}.
1226     * <li>If {@code z} is ±0 + i0, returns 1 + i0.
1227     * <li>If {@code z} is x + i∞ for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1228     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1229     * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
1230     * <li>If {@code z} is −∞ + iy for finite y, returns +0 cis(y) (see {@link #ofCis(double)}).
1231     * <li>If {@code z} is +∞ + iy for finite nonzero y, returns +∞ cis(y).
1232     * <li>If {@code z} is −∞ + i∞, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
1233     * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
1234     * <li>If {@code z} is −∞ + iNaN, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
1235     * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
1236     * <li>If {@code z} is NaN + i0, returns NaN + i0.
1237     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
1238     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1239     * </ul>
1240     *
1241     * <p>Implements the formula:
1242     *
1243     * <p>\[ \exp(x + iy) = e^x (\cos(y) + i \sin(y)) \]
1244     *
1245     * @return The exponential of this complex number.
1246     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</a>
1247     */
1248    public Complex exp() {
1249        if (Double.isInfinite(real)) {
1250            // Set the scale factor applied to cis(y)
1251            final double zeroOrInf;
1252            if (real < 0) {
1253                if (!Double.isFinite(imaginary)) {
1254                    // (−∞ + i∞) or (−∞ + iNaN) returns (±0 ± i0) (where the signs of the
1255                    // real and imaginary parts of the result are unspecified).
1256                    // Here we preserve the conjugate equality.
1257                    return new Complex(0, Math.copySign(0, imaginary));
1258                }
1259                // (−∞ + iy) returns +0 cis(y), for finite y
1260                zeroOrInf = 0;
1261            } else {
1262                // (+∞ + i0) returns +∞ + i0.
1263                if (imaginary == 0) {
1264                    return this;
1265                }
1266                // (+∞ + i∞) or (+∞ + iNaN) returns (±∞ + iNaN) and raises the invalid
1267                // floating-point exception (where the sign of the real part of the
1268                // result is unspecified).
1269                if (!Double.isFinite(imaginary)) {
1270                    return new Complex(real, Double.NaN);
1271                }
1272                // (+∞ + iy) returns (+∞ cis(y)), for finite nonzero y.
1273                zeroOrInf = real;
1274            }
1275            return new Complex(zeroOrInf * Math.cos(imaginary),
1276                               zeroOrInf * Math.sin(imaginary));
1277        } else if (Double.isNaN(real)) {
1278            // (NaN + i0) returns (NaN + i0)
1279            // (NaN + iy) returns (NaN + iNaN) and optionally raises the invalid floating-point exception
1280            // (NaN + iNaN) returns (NaN + iNaN)
1281            return imaginary == 0 ? this : NAN;
1282        } else if (!Double.isFinite(imaginary)) {
1283            // (x + i∞) or (x + iNaN) returns (NaN + iNaN) and raises the invalid
1284            // floating-point exception, for finite x.
1285            return NAN;
1286        }
1287        // real and imaginary are finite.
1288        // Compute e^a * (cos(b) + i sin(b)).
1289
1290        // Special case:
1291        // (±0 + i0) returns (1 + i0)
1292        final double exp = Math.exp(real);
1293        if (imaginary == 0) {
1294            return new Complex(exp, imaginary);
1295        }
1296        return new Complex(exp * Math.cos(imaginary),
1297                           exp * Math.sin(imaginary));
1298    }
1299
1300    /**
1301     * Returns the
1302     * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html">
1303     * natural logarithm</a> of this complex number.
1304     *
1305     * <p>The natural logarithm of \( z \) is unbounded along the real axis and
1306     * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
1307     * natural logarithm has a branch cut along the negative real axis \( (-infty,0] \).
1308     * Special cases:
1309     *
1310     * <ul>
1311     * <li>{@code z.conj().log() == z.log().conj()}.
1312     * <li>If {@code z} is −0 + i0, returns −∞ + iπ ("divide-by-zero" floating-point operation).
1313     * <li>If {@code z} is +0 + i0, returns −∞ + i0 ("divide-by-zero" floating-point operation).
1314     * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
1315     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1316     * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +∞ + iπ.
1317     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
1318     * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
1319     * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
1320     * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
1321     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
1322     * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
1323     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1324     * </ul>
1325     *
1326     * <p>Implements the formula:
1327     *
1328     * <p>\[ \ln(z) = \ln |z| + i \arg(z) \]
1329     *
1330     * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
1331     *
1332     * <p>The implementation is based on the method described in:</p>
1333     * <blockquote>
1334     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
1335     * Implementing complex elementary functions using exception handling.
1336     * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
1337     * </blockquote>
1338     *
1339     * @return The natural logarithm of this complex number.
1340     * @see Math#log(double)
1341     * @see #abs()
1342     * @see #arg()
1343     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
1344     */
1345    public Complex log() {
1346        return log(Math::log, HALF, LN_2, Complex::ofCartesian);
1347    }
1348
1349    /**
1350     * Returns the base 10
1351     * <a href="http://mathworld.wolfram.com/CommonLogarithm.html">
1352     * common logarithm</a> of this complex number.
1353     *
1354     * <p>The common logarithm of \( z \) is unbounded along the real axis and
1355     * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
1356     * common logarithm has a branch cut along the negative real axis \( (-infty,0] \).
1357     * Special cases are as defined in the {@link #log() natural logarithm}:
1358     *
1359     * <p>Implements the formula:
1360     *
1361     * <p>\[ \log_{10}(z) = \log_{10} |z| + i \arg(z) \]
1362     *
1363     * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
1364     *
1365     * @return The base 10 logarithm of this complex number.
1366     * @see Math#log10(double)
1367     * @see #abs()
1368     * @see #arg()
1369     */
1370    public Complex log10() {
1371        return log(Math::log10, LOG_10E_O_2, LOG10_2, Complex::ofCartesian);
1372    }
1373
1374    /**
1375     * Returns the logarithm of this complex number using the provided function.
1376     * Implements the formula:
1377     *
1378     * <pre>
1379     *   log(x + i y) = log(|x + i y|) + i arg(x + i y)</pre>
1380     *
1381     * <p>Warning: The argument {@code logOf2} must be equal to {@code log(2)} using the
1382     * provided log function otherwise scaling using powers of 2 in the case of overflow
1383     * will be incorrect. This is provided as an internal optimisation.
1384     *
1385     * @param log Log function.
1386     * @param logOfeOver2 The log function applied to e, then divided by 2.
1387     * @param logOf2 The log function applied to 2.
1388     * @param constructor Constructor for the returned complex.
1389     * @return The logarithm of this complex number.
1390     * @see #abs()
1391     * @see #arg()
1392     */
1393    private Complex log(DoubleUnaryOperator log, double logOfeOver2, double logOf2, ComplexConstructor constructor) {
1394        // Handle NaN
1395        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
1396            // Return NaN unless infinite
1397            if (isInfinite()) {
1398                return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
1399            }
1400            // No-use of the input constructor
1401            return NAN;
1402        }
1403
1404        // Returns the real part:
1405        // log(sqrt(x^2 + y^2))
1406        // log(x^2 + y^2) / 2
1407
1408        // Compute with positive values
1409        double x = Math.abs(real);
1410        double y = Math.abs(imaginary);
1411
1412        // Find the larger magnitude.
1413        if (x < y) {
1414            final double tmp = x;
1415            x = y;
1416            y = tmp;
1417        }
1418
1419        if (x == 0) {
1420            // Handle zero: raises the ‘‘divide-by-zero’’ floating-point exception.
1421            return constructor.create(Double.NEGATIVE_INFINITY,
1422                                      negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary);
1423        }
1424
1425        double re;
1426
1427        // This alters the implementation of Hull et al (1994) which used a standard
1428        // precision representation of |z|: sqrt(x*x + y*y).
1429        // This formula should use the same definition of the magnitude returned
1430        // by Complex.abs() which is a high precision computation with scaling.
1431        // The checks for overflow thus only require ensuring the output of |z|
1432        // will not overflow or underflow.
1433
1434        if (x > HALF && x < ROOT2) {
1435            // x^2+y^2 close to 1. Use log1p(x^2+y^2 - 1) / 2.
1436            re = Math.log1p(x2y2m1(x, y)) * logOfeOver2;
1437        } else {
1438            // Check for over/underflow in |z|
1439            // When scaling:
1440            // log(a / b) = log(a) - log(b)
1441            // So initialize the result with the log of the scale factor.
1442            re = 0;
1443            if (x > Double.MAX_VALUE / 2) {
1444                // Potential overflow.
1445                if (isPosInfinite(x)) {
1446                    // Handle infinity
1447                    return constructor.create(x, arg());
1448                }
1449                // Scale down.
1450                x /= 2;
1451                y /= 2;
1452                // log(2)
1453                re = logOf2;
1454            } else if (y < Double.MIN_NORMAL) {
1455                // Potential underflow.
1456                if (y == 0) {
1457                    // Handle real only number
1458                    return constructor.create(log.applyAsDouble(x), arg());
1459                }
1460                // Scale up sub-normal numbers to make them normal by scaling by 2^54,
1461                // i.e. more than the mantissa digits.
1462                x *= 0x1.0p54;
1463                y *= 0x1.0p54;
1464                // log(2^-54) = -54 * log(2)
1465                re = -54 * logOf2;
1466            }
1467            re += log.applyAsDouble(abs(x, y));
1468        }
1469
1470        // All ISO C99 edge cases for the imaginary are satisfied by the Math library.
1471        return constructor.create(re, arg());
1472    }
1473
1474    /**
1475     * Returns the complex power of this complex number raised to the power of {@code x}.
1476     * Implements the formula:
1477     *
1478     * <p>\[ z^x = e^{x \ln(z)} \]
1479     *
1480     * <p>If this complex number is zero then this method returns zero if {@code x} is positive
1481     * in the real component and zero in the imaginary component;
1482     * otherwise it returns NaN + iNaN.
1483     *
1484     * @param  x The exponent to which this complex number is to be raised.
1485     * @return This complex number raised to the power of {@code x}.
1486     * @see #log()
1487     * @see #multiply(Complex)
1488     * @see #exp()
1489     * @see <a href="http://mathworld.wolfram.com/ComplexExponentiation.html">Complex exponentiation</a>
1490     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
1491     */
1492    public Complex pow(Complex x) {
1493        if (real == 0 &&
1494            imaginary == 0) {
1495            // This value is zero. Test the other.
1496            if (x.real > 0 &&
1497                x.imaginary == 0) {
1498                // 0 raised to positive number is 0
1499                return ZERO;
1500            }
1501            // 0 raised to anything else is NaN
1502            return NAN;
1503        }
1504        return log().multiply(x).exp();
1505    }
1506
1507    /**
1508     * Returns the complex power of this complex number raised to the power of {@code x},
1509     * with {@code x} interpreted as a real number.
1510     * Implements the formula:
1511     *
1512     * <p>\[ z^x = e^{x \ln(z)} \]
1513     *
1514     * <p>If this complex number is zero then this method returns zero if {@code x} is positive;
1515     * otherwise it returns NaN + iNaN.
1516     *
1517     * @param  x The exponent to which this complex number is to be raised.
1518     * @return This complex number raised to the power of {@code x}.
1519     * @see #log()
1520     * @see #multiply(double)
1521     * @see #exp()
1522     * @see #pow(Complex)
1523     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
1524     */
1525    public Complex pow(double x) {
1526        if (real == 0 &&
1527            imaginary == 0) {
1528            // This value is zero. Test the other.
1529            if (x > 0) {
1530                // 0 raised to positive number is 0
1531                return ZERO;
1532            }
1533            // 0 raised to anything else is NaN
1534            return NAN;
1535        }
1536        return log().multiply(x).exp();
1537    }
1538
1539    /**
1540     * Returns the
1541     * <a href="http://mathworld.wolfram.com/SquareRoot.html">
1542     * square root</a> of this complex number.
1543     *
1544     * <p>\[ \sqrt{x + iy} = \frac{1}{2} \sqrt{2} \left( \sqrt{ \sqrt{x^2 + y^2} + x } + i\ \text{sgn}(y) \sqrt{ \sqrt{x^2 + y^2} - x } \right) \]
1545     *
1546     * <p>The square root of \( z \) is in the range \( [0, +\infty) \) along the real axis and
1547     * is unbounded along the imaginary axis. The imaginary part of the square root has a
1548     * branch cut along the negative real axis \( (-infty,0) \). Special cases:
1549     *
1550     * <ul>
1551     * <li>{@code z.conj().sqrt() == z.sqrt().conj()}.
1552     * <li>If {@code z} is ±0 + i0, returns +0 + i0.
1553     * <li>If {@code z} is x + i∞ for all x (including NaN), returns +∞ + i∞.
1554     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
1555     * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +0 + i∞.
1556     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
1557     * <li>If {@code z} is −∞ + iNaN, returns NaN ± i∞ (where the sign of the imaginary part of the result is unspecified).
1558     * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
1559     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
1560     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1561     * </ul>
1562     *
1563     * <p>Implements the following algorithm to compute \( \sqrt{x + iy} \):
1564     * <ol>
1565     * <li>Let \( t = \sqrt{2 (|x| + |x + iy|)} \)
1566     * <li>if \( x \geq 0 \) return \( \frac{t}{2} + i \frac{y}{t} \)
1567     * <li>else return \( \frac{|y|}{t} + i\ \text{sgn}(y) \frac{t}{2} \)
1568     * </ol>
1569     * where:
1570     * <ul>
1571     * <li>\( |x| =\ \){@link Math#abs(double) abs}(x)
1572     * <li>\( |x + y i| =\ \){@link Complex#abs}
1573     * <li>\( \text{sgn}(y) =\ \){@link Math#copySign(double,double) copySign}(1.0, y)
1574     * </ul>
1575     *
1576     * <p>The implementation is overflow and underflow safe based on the method described in:</p>
1577     * <blockquote>
1578     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
1579     * Implementing complex elementary functions using exception handling.
1580     * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
1581     * </blockquote>
1582     *
1583     * @return The square root of this complex number.
1584     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sqrt/">Sqrt</a>
1585     */
1586    public Complex sqrt() {
1587        return sqrt(real, imaginary);
1588    }
1589
1590    /**
1591     * Returns the square root of the complex number {@code sqrt(x + i y)}.
1592     *
1593     * @param real Real component.
1594     * @param imaginary Imaginary component.
1595     * @return The square root of the complex number.
1596     */
1597    private static Complex sqrt(double real, double imaginary) {
1598        // Handle NaN
1599        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
1600            // Check for infinite
1601            if (Double.isInfinite(imaginary)) {
1602                return new Complex(Double.POSITIVE_INFINITY, imaginary);
1603            }
1604            if (Double.isInfinite(real)) {
1605                if (real == Double.NEGATIVE_INFINITY) {
1606                    return new Complex(Double.NaN, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
1607                }
1608                return new Complex(Double.POSITIVE_INFINITY, Double.NaN);
1609            }
1610            return NAN;
1611        }
1612
1613        // Compute with positive values and determine sign at the end
1614        final double x = Math.abs(real);
1615        final double y = Math.abs(imaginary);
1616
1617        // Compute
1618        final double t;
1619
1620        // This alters the implementation of Hull et al (1994) which used a standard
1621        // precision representation of |z|: sqrt(x*x + y*y).
1622        // This formula should use the same definition of the magnitude returned
1623        // by Complex.abs() which is a high precision computation with scaling.
1624        // Worry about overflow if 2 * (|z| + |x|) will overflow.
1625        // Worry about underflow if |z| or |x| are sub-normal components.
1626
1627        if (inRegion(x, y, Double.MIN_NORMAL, SQRT_SAFE_UPPER)) {
1628            // No over/underflow
1629            t = Math.sqrt(2 * (abs(x, y) + x));
1630        } else {
1631            // Potential over/underflow. First check infinites and real/imaginary only.
1632
1633            // Check for infinite
1634            if (isPosInfinite(y)) {
1635                return new Complex(Double.POSITIVE_INFINITY, imaginary);
1636            } else if (isPosInfinite(x)) {
1637                if (real == Double.NEGATIVE_INFINITY) {
1638                    return new Complex(0, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
1639                }
1640                return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0, imaginary));
1641            } else if (y == 0) {
1642                // Real only
1643                final double sqrtAbs = Math.sqrt(x);
1644                if (real < 0) {
1645                    return new Complex(0, Math.copySign(sqrtAbs, imaginary));
1646                }
1647                return new Complex(sqrtAbs, imaginary);
1648            } else if (x == 0) {
1649                // Imaginary only. This sets the two components to the same magnitude.
1650                // Note: In polar coordinates this does not happen:
1651                // real = sqrt(abs()) * Math.cos(arg() / 2)
1652                // imag = sqrt(abs()) * Math.sin(arg() / 2)
1653                // arg() / 2 = pi/4 and cos and sin should both return sqrt(2)/2 but
1654                // are different by 1 ULP.
1655                final double sqrtAbs = Math.sqrt(y) * ONE_OVER_ROOT2;
1656                return new Complex(sqrtAbs, Math.copySign(sqrtAbs, imaginary));
1657            } else {
1658                // Over/underflow.
1659                // Full scaling is not required as this is done in the hypotenuse function.
1660                // Keep the number as big as possible for maximum precision in the second sqrt.
1661                // Note if we scale by an even power of 2, we can re-scale by sqrt of the number.
1662                // a = sqrt(b)
1663                // a = sqrt(b/4) * sqrt(4)
1664
1665                final double rescale;
1666                final double sx;
1667                final double sy;
1668                if (Math.max(x, y) > SQRT_SAFE_UPPER) {
1669                    // Overflow. Scale down by 16 and rescale by sqrt(16).
1670                    sx = x / 16;
1671                    sy = y / 16;
1672                    rescale = 4;
1673                } else {
1674                    // Sub-normal numbers. Make them normal by scaling by 2^54,
1675                    // i.e. more than the mantissa digits, and rescale by sqrt(2^54) = 2^27.
1676                    sx = x * 0x1.0p54;
1677                    sy = y * 0x1.0p54;
1678                    rescale = 0x1.0p-27;
1679                }
1680                t = rescale * Math.sqrt(2 * (abs(sx, sy) + sx));
1681            }
1682        }
1683
1684        if (real >= 0) {
1685            return new Complex(t / 2, imaginary / t);
1686        }
1687        return new Complex(y / t, Math.copySign(t / 2, imaginary));
1688    }
1689
1690    /**
1691     * Returns the
1692     * <a href="http://mathworld.wolfram.com/Sine.html">
1693     * sine</a> of this complex number.
1694     *
1695     * <p>\[ \sin(z) = \frac{1}{2} i \left( e^{-iz} - e^{iz} \right) \]
1696     *
1697     * <p>This is an odd function: \( \sin(z) = -\sin(-z) \).
1698     * The sine is an entire function and requires no branch cuts.
1699     *
1700     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
1701     *
1702     * <p>\[ \sin(x + iy) = \sin(x)\cosh(y) + i \cos(x)\sinh(y) \]
1703     *
1704     * <p>As per the C99 standard this function is computed using the trigonomic identity:
1705     *
1706     * <p>\[ \sin(z) = -i \sinh(iz) \]
1707     *
1708     * @return The sine of this complex number.
1709     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sin/">Sin</a>
1710     */
1711    public Complex sin() {
1712        // Define in terms of sinh
1713        // sin(z) = -i sinh(iz)
1714        // Multiply this number by I, compute sinh, then multiply by back
1715        return sinh(-imaginary, real, Complex::multiplyNegativeI);
1716    }
1717
1718    /**
1719     * Returns the
1720     * <a href="http://mathworld.wolfram.com/Cosine.html">
1721     * cosine</a> of this complex number.
1722     *
1723     * <p>\[ \cos(z) = \frac{1}{2} \left( e^{iz} + e^{-iz} \right) \]
1724     *
1725     * <p>This is an even function: \( \cos(z) = \cos(-z) \).
1726     * The cosine is an entire function and requires no branch cuts.
1727     *
1728     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
1729     *
1730     * <p>\[ \cos(x + iy) = \cos(x)\cosh(y) - i \sin(x)\sinh(y) \]
1731     *
1732     * <p>As per the C99 standard this function is computed using the trigonomic identity:
1733     *
1734     * <p>\[ cos(z) = cosh(iz) \]
1735     *
1736     * @return The cosine of this complex number.
1737     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cos/">Cos</a>
1738     */
1739    public Complex cos() {
1740        // Define in terms of cosh
1741        // cos(z) = cosh(iz)
1742        // Multiply this number by I and compute cosh.
1743        return cosh(-imaginary, real, Complex::ofCartesian);
1744    }
1745
1746    /**
1747     * Returns the
1748     * <a href="http://mathworld.wolfram.com/Tangent.html">
1749     * tangent</a> of this complex number.
1750     *
1751     * <p>\[ \tan(z) = \frac{i(e^{-iz} - e^{iz})}{e^{-iz} + e^{iz}} \]
1752     *
1753     * <p>This is an odd function: \( \tan(z) = -\tan(-z) \).
1754     * The tangent is an entire function and requires no branch cuts.
1755     *
1756     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:</p>
1757     * \[ \tan(x + iy) = \frac{\sin(2x)}{\cos(2x)+\cosh(2y)} + i \frac{\sinh(2y)}{\cos(2x)+\cosh(2y)} \]
1758     *
1759     * <p>As per the C99 standard this function is computed using the trigonomic identity:</p>
1760     * \[ \tan(z) = -i \tanh(iz) \]
1761     *
1762     * @return The tangent of this complex number.
1763     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tan/">Tangent</a>
1764     */
1765    public Complex tan() {
1766        // Define in terms of tanh
1767        // tan(z) = -i tanh(iz)
1768        // Multiply this number by I, compute tanh, then multiply by back
1769        return tanh(-imaginary, real, Complex::multiplyNegativeI);
1770    }
1771
1772    /**
1773     * Returns the
1774     * <a href="http://mathworld.wolfram.com/InverseSine.html">
1775     * inverse sine</a> of this complex number.
1776     *
1777     * <p>\[ \sin^{-1}(z) = - i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
1778     *
1779     * <p>The inverse sine of \( z \) is unbounded along the imaginary axis and
1780     * in the range \( [-\pi, \pi] \) along the real axis. Special cases are handled
1781     * as if the operation is implemented using \( \sin^{-1}(z) = -i \sinh^{-1}(iz) \).
1782     *
1783     * <p>The inverse sine is a multivalued function and requires a branch cut in
1784     * the complex plane; the cut is conventionally placed at the line segments
1785     * \( (\infty,-1) \) and \( (1,\infty) \) of the real axis.
1786     *
1787     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
1788     *
1789     * <p>\[ \begin{aligned}
1790     *   \sin^{-1}(z) &amp;= \sin^{-1}(B) + i\ \text{sgn}(y)\ln \left(A + \sqrt{A^2-1} \right) \\
1791     *   A &amp;= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
1792     *   B &amp;= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \end{aligned} \]
1793     *
1794     * <p>where \( \text{sgn}(y) \) is the sign function implemented using
1795     * {@link Math#copySign(double,double) copySign(1.0, y)}.
1796     *
1797     * <p>The implementation is based on the method described in:</p>
1798     * <blockquote>
1799     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
1800     * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
1801     * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
1802     * </blockquote>
1803     *
1804     * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
1805     * {@code c++} implementation {@code <boost/math/complex/asin.hpp>}.
1806     *
1807     * @return The inverse sine of this complex number.
1808     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSin/">ArcSin</a>
1809     */
1810    public Complex asin() {
1811        return asin(real, imaginary, Complex::ofCartesian);
1812    }
1813
1814    /**
1815     * Returns the inverse sine of the complex number.
1816     *
1817     * <p>This function exists to allow implementation of the identity
1818     * {@code asinh(z) = -i asin(iz)}.
1819     *
1820     * <p>Adapted from {@code <boost/math/complex/asin.hpp>}. This method only (and not
1821     * invoked methods within) is distributed under the Boost Software License V1.0.
1822     * The original notice is shown below and the licence is shown in full in LICENSE:
1823     * <pre>
1824     * (C) Copyright John Maddock 2005.
1825     * Distributed under the Boost Software License, Version 1.0. (See accompanying
1826     * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
1827     * </pre>
1828     *
1829     * @param real Real part.
1830     * @param imaginary Imaginary part.
1831     * @param constructor Constructor.
1832     * @return The inverse sine of this complex number.
1833     */
1834    private static Complex asin(final double real, final double imaginary,
1835                                final ComplexConstructor constructor) {
1836        // Compute with positive values and determine sign at the end
1837        final double x = Math.abs(real);
1838        final double y = Math.abs(imaginary);
1839        // The result (without sign correction)
1840        final double re;
1841        final double im;
1842
1843        // Handle C99 special cases
1844        if (Double.isNaN(x)) {
1845            if (isPosInfinite(y)) {
1846                re = x;
1847                im = y;
1848            } else {
1849                // No-use of the input constructor
1850                return NAN;
1851            }
1852        } else if (Double.isNaN(y)) {
1853            if (x == 0) {
1854                re = 0;
1855                im = y;
1856            } else if (isPosInfinite(x)) {
1857                re = y;
1858                im = x;
1859            } else {
1860                // No-use of the input constructor
1861                return NAN;
1862            }
1863        } else if (isPosInfinite(x)) {
1864            re = isPosInfinite(y) ? PI_OVER_4 : PI_OVER_2;
1865            im = x;
1866        } else if (isPosInfinite(y)) {
1867            re = 0;
1868            im = y;
1869        } else {
1870            // Special case for real numbers:
1871            if (y == 0 && x <= 1) {
1872                return constructor.create(Math.asin(real), imaginary);
1873            }
1874
1875            final double xp1 = x + 1;
1876            final double xm1 = x - 1;
1877
1878            if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
1879                final double yy = y * y;
1880                final double r = Math.sqrt(xp1 * xp1 + yy);
1881                final double s = Math.sqrt(xm1 * xm1 + yy);
1882                final double a = 0.5 * (r + s);
1883                final double b = x / a;
1884
1885                if (b <= B_CROSSOVER) {
1886                    re = Math.asin(b);
1887                } else {
1888                    final double apx = a + x;
1889                    if (x <= 1) {
1890                        re = Math.atan(x / Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))));
1891                    } else {
1892                        re = Math.atan(x / (y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))));
1893                    }
1894                }
1895
1896                if (a <= A_CROSSOVER) {
1897                    final double am1;
1898                    if (x < 1) {
1899                        am1 = 0.5 * (yy / (r + xp1) + yy / (s - xm1));
1900                    } else {
1901                        am1 = 0.5 * (yy / (r + xp1) + (s + xm1));
1902                    }
1903                    im = Math.log1p(am1 + Math.sqrt(am1 * (a + 1)));
1904                } else {
1905                    im = Math.log(a + Math.sqrt(a * a - 1));
1906                }
1907            } else {
1908                // Hull et al: Exception handling code from figure 4
1909                if (y <= (EPSILON * Math.abs(xm1))) {
1910                    if (x < 1) {
1911                        re = Math.asin(x);
1912                        im = y / Math.sqrt(xp1 * (1 - x));
1913                    } else {
1914                        re = PI_OVER_2;
1915                        if ((Double.MAX_VALUE / xp1) > xm1) {
1916                            // xp1 * xm1 won't overflow:
1917                            im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
1918                        } else {
1919                            im = LN_2 + Math.log(x);
1920                        }
1921                    }
1922                } else if (y <= SAFE_MIN) {
1923                    // Hull et al: Assume x == 1.
1924                    // True if:
1925                    // E^2 > 8*sqrt(u)
1926                    //
1927                    // E = Machine epsilon: (1 + epsilon) = 1
1928                    // u = Double.MIN_NORMAL
1929                    re = PI_OVER_2 - Math.sqrt(y);
1930                    im = Math.sqrt(y);
1931                } else if (EPSILON * y - 1 >= x) {
1932                    // Possible underflow:
1933                    re = x / y;
1934                    im = LN_2 + Math.log(y);
1935                } else if (x > 1) {
1936                    re = Math.atan(x / y);
1937                    final double xoy = x / y;
1938                    im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
1939                } else {
1940                    final double a = Math.sqrt(1 + y * y);
1941                    // Possible underflow:
1942                    re = x / a;
1943                    im = 0.5 * Math.log1p(2 * y * (y + a));
1944                }
1945            }
1946        }
1947
1948        return constructor.create(changeSign(re, real),
1949                                  changeSign(im, imaginary));
1950    }
1951
1952    /**
1953     * Returns the
1954     * <a href="http://mathworld.wolfram.com/InverseCosine.html">
1955     * inverse cosine</a> of this complex number.
1956     *
1957     * <p>\[ \cos^{-1}(z) = \frac{\pi}{2} + i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
1958     *
1959     * <p>The inverse cosine of \( z \) is in the range \( [0, \pi) \) along the real axis and
1960     * unbounded along the imaginary axis. Special cases:
1961     *
1962     * <ul>
1963     * <li>{@code z.conj().acos() == z.acos().conj()}.
1964     * <li>If {@code z} is ±0 + i0, returns π/2 − i0.
1965     * <li>If {@code z} is ±0 + iNaN, returns π/2 + iNaN.
1966     * <li>If {@code z} is x + i∞ for finite x, returns π/2 − i∞.
1967     * <li>If {@code z} is x + iNaN, returns NaN + iNaN ("invalid" floating-point operation).
1968     * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns π − i∞.
1969     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +0 − i∞.
1970     * <li>If {@code z} is −∞ + i∞, returns 3π/4 − i∞.
1971     * <li>If {@code z} is +∞ + i∞, returns π/4 − i∞.
1972     * <li>If {@code z} is ±∞ + iNaN, returns NaN ± i∞ where the sign of the imaginary part of the result is unspecified.
1973     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
1974     * <li>If {@code z} is NaN + i∞, returns NaN − i∞.
1975     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
1976     * </ul>
1977     *
1978     * <p>The inverse cosine is a multivalued function and requires a branch cut in
1979     * the complex plane; the cut is conventionally placed at the line segments
1980     * \( (-\infty,-1) \) and \( (1,\infty) \) of the real axis.
1981     *
1982     * <p>This function is implemented using real \( x \) and imaginary \( y \) parts:
1983     *
1984     * <p>\[ \begin{aligned}
1985     *   \cos^{-1}(z) &amp;= \cos^{-1}(B) - i\ \text{sgn}(y) \ln\left(A + \sqrt{A^2-1}\right) \\
1986     *   A &amp;= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
1987     *   B &amp;= \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \end{aligned} \]
1988     *
1989     * <p>where \( \text{sgn}(y) \) is the sign function implemented using
1990     * {@link Math#copySign(double,double) copySign(1.0, y)}.
1991     *
1992     * <p>The implementation is based on the method described in:</p>
1993     * <blockquote>
1994     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
1995     * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
1996     * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
1997     * </blockquote>
1998     *
1999     * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
2000     * {@code c++} implementation {@code <boost/math/complex/acos.hpp>}.
2001     *
2002     * @return The inverse cosine of this complex number.
2003     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCos/">ArcCos</a>
2004     */
2005    public Complex acos() {
2006        return acos(real, imaginary, Complex::ofCartesian);
2007    }
2008
2009    /**
2010     * Returns the inverse cosine of the complex number.
2011     *
2012     * <p>This function exists to allow implementation of the identity
2013     * {@code acosh(z) = +-i acos(z)}.
2014     *
2015     * <p>Adapted from {@code <boost/math/complex/acos.hpp>}. This method only (and not
2016     * invoked methods within) is distributed under the Boost Software License V1.0.
2017     * The original notice is shown below and the licence is shown in full in LICENSE:
2018     * <pre>
2019     * (C) Copyright John Maddock 2005.
2020     * Distributed under the Boost Software License, Version 1.0. (See accompanying
2021     * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
2022     * </pre>
2023     *
2024     * @param real Real part.
2025     * @param imaginary Imaginary part.
2026     * @param constructor Constructor.
2027     * @return The inverse cosine of the complex number.
2028     */
2029    private static Complex acos(final double real, final double imaginary,
2030                                final ComplexConstructor constructor) {
2031        // Compute with positive values and determine sign at the end
2032        final double x = Math.abs(real);
2033        final double y = Math.abs(imaginary);
2034        // The result (without sign correction)
2035        final double re;
2036        final double im;
2037
2038        // Handle C99 special cases
2039        if (isPosInfinite(x)) {
2040            if (isPosInfinite(y)) {
2041                re = PI_OVER_4;
2042                im = y;
2043            } else if (Double.isNaN(y)) {
2044                // sign of the imaginary part of the result is unspecified
2045                return constructor.create(imaginary, real);
2046            } else {
2047                re = 0;
2048                im = Double.POSITIVE_INFINITY;
2049            }
2050        } else if (Double.isNaN(x)) {
2051            if (isPosInfinite(y)) {
2052                return constructor.create(x, -imaginary);
2053            }
2054            // No-use of the input constructor
2055            return NAN;
2056        } else if (isPosInfinite(y)) {
2057            re = PI_OVER_2;
2058            im = y;
2059        } else if (Double.isNaN(y)) {
2060            return constructor.create(x == 0 ? PI_OVER_2 : y, y);
2061        } else {
2062            // Special case for real numbers:
2063            if (y == 0 && x <= 1) {
2064                return constructor.create(x == 0 ? PI_OVER_2 : Math.acos(real), -imaginary);
2065            }
2066
2067            final double xp1 = x + 1;
2068            final double xm1 = x - 1;
2069
2070            if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
2071                final double yy = y * y;
2072                final double r = Math.sqrt(xp1 * xp1 + yy);
2073                final double s = Math.sqrt(xm1 * xm1 + yy);
2074                final double a = 0.5 * (r + s);
2075                final double b = x / a;
2076
2077                if (b <= B_CROSSOVER) {
2078                    re = Math.acos(b);
2079                } else {
2080                    final double apx = a + x;
2081                    if (x <= 1) {
2082                        re = Math.atan(Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))) / x);
2083                    } else {
2084                        re = Math.atan((y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))) / x);
2085                    }
2086                }
2087
2088                if (a <= A_CROSSOVER) {
2089                    final double am1;
2090                    if (x < 1) {
2091                        am1 = 0.5 * (yy / (r + xp1) + yy / (s - xm1));
2092                    } else {
2093                        am1 = 0.5 * (yy / (r + xp1) + (s + xm1));
2094                    }
2095                    im = Math.log1p(am1 + Math.sqrt(am1 * (a + 1)));
2096                } else {
2097                    im = Math.log(a + Math.sqrt(a * a - 1));
2098                }
2099            } else {
2100                // Hull et al: Exception handling code from figure 6
2101                if (y <= (EPSILON * Math.abs(xm1))) {
2102                    if (x < 1) {
2103                        re = Math.acos(x);
2104                        im = y / Math.sqrt(xp1 * (1 - x));
2105                    } else {
2106                        // This deviates from Hull et al's paper as per
2107                        // https://svn.boost.org/trac/boost/ticket/7290
2108                        if ((Double.MAX_VALUE / xp1) > xm1) {
2109                            // xp1 * xm1 won't overflow:
2110                            re = y / Math.sqrt(xm1 * xp1);
2111                            im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
2112                        } else {
2113                            re = y / x;
2114                            im = LN_2 + Math.log(x);
2115                        }
2116                    }
2117                } else if (y <= SAFE_MIN) {
2118                    // Hull et al: Assume x == 1.
2119                    // True if:
2120                    // E^2 > 8*sqrt(u)
2121                    //
2122                    // E = Machine epsilon: (1 + epsilon) = 1
2123                    // u = Double.MIN_NORMAL
2124                    re = Math.sqrt(y);
2125                    im = Math.sqrt(y);
2126                } else if (EPSILON * y - 1 >= x) {
2127                    re = PI_OVER_2;
2128                    im = LN_2 + Math.log(y);
2129                } else if (x > 1) {
2130                    re = Math.atan(y / x);
2131                    final double xoy = x / y;
2132                    im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
2133                } else {
2134                    re = PI_OVER_2;
2135                    final double a = Math.sqrt(1 + y * y);
2136                    im = 0.5 * Math.log1p(2 * y * (y + a));
2137                }
2138            }
2139        }
2140
2141        return constructor.create(negative(real) ? Math.PI - re : re,
2142                                  negative(imaginary) ? im : -im);
2143    }
2144
2145    /**
2146     * Returns the
2147     * <a href="http://mathworld.wolfram.com/InverseTangent.html">
2148     * inverse tangent</a> of this complex number.
2149     *
2150     * <p>\[ \tan^{-1}(z) = \frac{i}{2} \ln \left( \frac{i + z}{i - z} \right) \]
2151     *
2152     * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the imaginary axis and
2153     * in the range \( [-\pi/2, \pi/2] \) along the real axis.
2154     *
2155     * <p>The inverse tangent is a multivalued function and requires a branch cut in
2156     * the complex plane; the cut is conventionally placed at the line segments
2157     * \( (i \infty,-i] \) and \( [i,i \infty) \) of the imaginary axis.
2158     *
2159     * <p>As per the C99 standard this function is computed using the trigonomic identity:
2160     * \[ \tan^{-1}(z) = -i \tanh^{-1}(iz) \]
2161     *
2162     * @return The inverse tangent of this complex number.
2163     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTan/">ArcTan</a>
2164     */
2165    public Complex atan() {
2166        // Define in terms of atanh
2167        // atan(z) = -i atanh(iz)
2168        // Multiply this number by I, compute atanh, then multiply by back
2169        return atanh(-imaginary, real, Complex::multiplyNegativeI);
2170    }
2171
2172    /**
2173     * Returns the
2174     * <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
2175     * hyperbolic sine</a> of this complex number.
2176     *
2177     * <p>\[ \sinh(z) = \frac{1}{2} \left( e^{z} - e^{-z} \right) \]
2178     *
2179     * <p>The hyperbolic sine of \( z \) is an entire function in the complex plane
2180     * and is periodic with respect to the imaginary component with period \( 2\pi i \).
2181     * Special cases:
2182     *
2183     * <ul>
2184     * <li>{@code z.conj().sinh() == z.sinh().conj()}.
2185     * <li>This is an odd function: \( \sinh(z) = -\sinh(-z) \).
2186     * <li>If {@code z} is +0 + i0, returns +0 + i0.
2187     * <li>If {@code z} is +0 + i∞, returns ±0 + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
2188     * <li>If {@code z} is +0 + iNaN, returns ±0 + iNaN (where the sign of the real part of the result is unspecified).
2189     * <li>If {@code z} is x + i∞ for positive finite x, returns NaN + iNaN ("invalid" floating-point operation).
2190     * <li>If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
2191     * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
2192     * <li>If {@code z} is +∞ + iy for positive finite y, returns +∞ cis(y) (see {@link #ofCis(double)}.
2193     * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
2194     * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
2195     * <li>If {@code z} is NaN + i0, returns NaN + i0.
2196     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
2197     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2198     * </ul>
2199     *
2200     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
2201     *
2202     * <p>\[ \sinh(x + iy) = \sinh(x)\cos(y) + i \cosh(x)\sin(y) \]
2203     *
2204     * @return The hyperbolic sine of this complex number.
2205     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sinh/">Sinh</a>
2206     */
2207    public Complex sinh() {
2208        return sinh(real, imaginary, Complex::ofCartesian);
2209    }
2210
2211    /**
2212     * Returns the hyperbolic sine of the complex number.
2213     *
2214     * <p>This function exists to allow implementation of the identity
2215     * {@code sin(z) = -i sinh(iz)}.<p>
2216     *
2217     * @param real Real part.
2218     * @param imaginary Imaginary part.
2219     * @param constructor Constructor.
2220     * @return The hyperbolic sine of the complex number.
2221     */
2222    private static Complex sinh(double real, double imaginary, ComplexConstructor constructor) {
2223        if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
2224            return constructor.create(real, Double.NaN);
2225        }
2226        if (real == 0) {
2227            // Imaginary-only sinh(iy) = i sin(y).
2228            if (Double.isFinite(imaginary)) {
2229                // Maintain periodic property with respect to the imaginary component.
2230                // sinh(+/-0.0) * cos(+/-x) = +/-0 * cos(x)
2231                return constructor.create(changeSign(real, Math.cos(imaginary)),
2232                                          Math.sin(imaginary));
2233            }
2234            // If imaginary is inf/NaN the sign of the real part is unspecified.
2235            // Returning the same real value maintains the conjugate equality.
2236            // It is not possible to also maintain the odd function (hence the unspecified sign).
2237            return constructor.create(real, Double.NaN);
2238        }
2239        if (imaginary == 0) {
2240            // Real-only sinh(x).
2241            return constructor.create(Math.sinh(real), imaginary);
2242        }
2243        final double x = Math.abs(real);
2244        if (x > SAFE_EXP) {
2245            // Approximate sinh/cosh(x) using exp^|x| / 2
2246            return coshsinh(x, real, imaginary, true, constructor);
2247        }
2248        // No overflow of sinh/cosh
2249        return constructor.create(Math.sinh(real) * Math.cos(imaginary),
2250                                  Math.cosh(real) * Math.sin(imaginary));
2251    }
2252
2253    /**
2254     * Returns the
2255     * <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
2256     * hyperbolic cosine</a> of this complex number.
2257     *
2258     * <p>\[ \cosh(z) = \frac{1}{2} \left( e^{z} + e^{-z} \right) \]
2259     *
2260     * <p>The hyperbolic cosine of \( z \) is an entire function in the complex plane
2261     * and is periodic with respect to the imaginary component with period \( 2\pi i \).
2262     * Special cases:
2263     *
2264     * <ul>
2265     * <li>{@code z.conj().cosh() == z.cosh().conj()}.
2266     * <li>This is an even function: \( \cosh(z) = \cosh(-z) \).
2267     * <li>If {@code z} is +0 + i0, returns 1 + i0.
2268     * <li>If {@code z} is +0 + i∞, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified; "invalid" floating-point operation).
2269     * <li>If {@code z} is +0 + iNaN, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified).
2270     * <li>If {@code z} is x + i∞ for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
2271     * <li>If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
2272     * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
2273     * <li>If {@code z} is +∞ + iy for finite nonzero y, returns +∞ cis(y) (see {@link #ofCis(double)}).
2274     * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
2275     * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
2276     * <li>If {@code z} is NaN + i0, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified).
2277     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
2278     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2279     * </ul>
2280     *
2281     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
2282     *
2283     * <p>\[ \cosh(x + iy) = \cosh(x)\cos(y) + i \sinh(x)\sin(y) \]
2284     *
2285     * @return The hyperbolic cosine of this complex number.
2286     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cosh/">Cosh</a>
2287     */
2288    public Complex cosh() {
2289        return cosh(real, imaginary, Complex::ofCartesian);
2290    }
2291
2292    /**
2293     * Returns the hyperbolic cosine of the complex number.
2294     *
2295     * <p>This function exists to allow implementation of the identity
2296     * {@code cos(z) = cosh(iz)}.<p>
2297     *
2298     * @param real Real part.
2299     * @param imaginary Imaginary part.
2300     * @param constructor Constructor.
2301     * @return The hyperbolic cosine of the complex number.
2302     */
2303    private static Complex cosh(double real, double imaginary, ComplexConstructor constructor) {
2304        // ISO C99: Preserve the even function by mapping to positive
2305        // f(z) = f(-z)
2306        if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
2307            return constructor.create(Math.abs(real), Double.NaN);
2308        }
2309        if (real == 0) {
2310            // Imaginary-only cosh(iy) = cos(y).
2311            if (Double.isFinite(imaginary)) {
2312                // Maintain periodic property with respect to the imaginary component.
2313                // sinh(+/-0.0) * sin(+/-x) = +/-0 * sin(x)
2314                return constructor.create(Math.cos(imaginary),
2315                                          changeSign(real, Math.sin(imaginary)));
2316            }
2317            // If imaginary is inf/NaN the sign of the imaginary part is unspecified.
2318            // Although not required by C99 changing the sign maintains the conjugate equality.
2319            // It is not possible to also maintain the even function (hence the unspecified sign).
2320            return constructor.create(Double.NaN, changeSign(real, imaginary));
2321        }
2322        if (imaginary == 0) {
2323            // Real-only cosh(x).
2324            // Change sign to preserve conjugate equality and even function.
2325            // sin(+/-0) * sinh(+/-x) = +/-0 * +/-a (sinh is monotonic and same sign)
2326            // => change the sign of imaginary using real. Handles special case of infinite real.
2327            // If real is NaN the sign of the imaginary part is unspecified.
2328            return constructor.create(Math.cosh(real), changeSign(imaginary, real));
2329        }
2330        final double x = Math.abs(real);
2331        if (x > SAFE_EXP) {
2332            // Approximate sinh/cosh(x) using exp^|x| / 2
2333            return coshsinh(x, real, imaginary, false, constructor);
2334        }
2335        // No overflow of sinh/cosh
2336        return constructor.create(Math.cosh(real) * Math.cos(imaginary),
2337                                  Math.sinh(real) * Math.sin(imaginary));
2338    }
2339
2340    /**
2341     * Compute cosh or sinh when the absolute real component |x| is large. In this case
2342     * cosh(x) and sinh(x) can be approximated by exp(|x|) / 2:
2343     *
2344     * <pre>
2345     * cosh(x+iy) real = (e^|x| / 2) * cos(y)
2346     * cosh(x+iy) imag = (e^|x| / 2) * sin(y) * sign(x)
2347     * sinh(x+iy) real = (e^|x| / 2) * cos(y) * sign(x)
2348     * sinh(x+iy) imag = (e^|x| / 2) * sin(y)
2349     * </pre>
2350     *
2351     * @param x Absolute real component |x|.
2352     * @param real Real part (x).
2353     * @param imaginary Imaginary part (y).
2354     * @param sinh Set to true to compute sinh, otherwise cosh.
2355     * @param constructor Constructor.
2356     * @return The hyperbolic sine/cosine of the complex number.
2357     */
2358    private static Complex coshsinh(double x, double real, double imaginary, boolean sinh,
2359                                    ComplexConstructor constructor) {
2360        // Always require the cos and sin.
2361        double re = Math.cos(imaginary);
2362        double im = Math.sin(imaginary);
2363        // Compute the correct function
2364        if (sinh) {
2365            re = changeSign(re, real);
2366        } else {
2367            im = changeSign(im, real);
2368        }
2369        // Multiply by (e^|x| / 2).
2370        // Overflow safe computation since sin/cos can be very small allowing a result
2371        // when e^x overflows: e^x / 2 = (e^m / 2) * e^m * e^(x-2m)
2372        if (x > SAFE_EXP * 3) {
2373            // e^x > e^m * e^m * e^m
2374            // y * (e^m / 2) * e^m * e^m will overflow when starting with Double.MIN_VALUE.
2375            // Note: Do not multiply by +inf to safeguard against sin(y)=0.0 which
2376            // will create 0 * inf = nan.
2377            re *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
2378            im *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
2379        } else {
2380            // Initial part of (e^x / 2) using (e^m / 2)
2381            re *= EXP_M / 2;
2382            im *= EXP_M / 2;
2383            final double xm;
2384            if (x > SAFE_EXP * 2) {
2385                // e^x = e^m * e^m * e^(x-2m)
2386                re *= EXP_M;
2387                im *= EXP_M;
2388                xm = x - SAFE_EXP * 2;
2389            } else {
2390                // e^x = e^m * e^(x-m)
2391                xm = x - SAFE_EXP;
2392            }
2393            final double exp = Math.exp(xm);
2394            re *= exp;
2395            im *= exp;
2396        }
2397        return constructor.create(re, im);
2398    }
2399
2400    /**
2401     * Returns the
2402     * <a href="http://mathworld.wolfram.com/HyperbolicTangent.html">
2403     * hyperbolic tangent</a> of this complex number.
2404     *
2405     * <p>\[ \tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}} \]
2406     *
2407     * <p>The hyperbolic tangent of \( z \) is an entire function in the complex plane
2408     * and is periodic with respect to the imaginary component with period \( \pi i \)
2409     * and has poles of the first order along the imaginary line, at coordinates
2410     * \( (0, \pi(\frac{1}{2} + n)) \).
2411     * Note that the {@code double} floating-point representation is unable to exactly represent
2412     * \( \pi/2 \) and there is no value for which a pole error occurs. Special cases:
2413     *
2414     * <ul>
2415     * <li>{@code z.conj().tanh() == z.tanh().conj()}.
2416     * <li>This is an odd function: \( \tanh(z) = -\tanh(-z) \).
2417     * <li>If {@code z} is +0 + i0, returns +0 + i0.
2418     * <li>If {@code z} is 0 + i∞, returns 0 + iNaN.
2419     * <li>If {@code z} is x + i∞ for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
2420     * <li>If {@code z} is 0 + iNaN, returns 0 + iNAN.
2421     * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
2422     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns 1 + i0 sin(2y).
2423     * <li>If {@code z} is +∞ + i∞, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
2424     * <li>If {@code z} is +∞ + iNaN, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
2425     * <li>If {@code z} is NaN + i0, returns NaN + i0.
2426     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
2427     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2428     * </ul>
2429     *
2430     * <p>Special cases include the technical corrigendum
2431     * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
2432     * DR 471: Complex math functions cacosh and ctanh</a>.
2433     *
2434     * <p>This is defined using real \( x \) and imaginary \( y \) parts:
2435     *
2436     * <p>\[ \tan(x + iy) = \frac{\sinh(2x)}{\cosh(2x)+\cos(2y)} + i \frac{\sin(2y)}{\cosh(2x)+\cos(2y)} \]
2437     *
2438     * <p>The implementation uses double-angle identities to avoid overflow of {@code 2x}
2439     * and {@code 2y}.
2440     *
2441     * @return The hyperbolic tangent of this complex number.
2442     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
2443     */
2444    public Complex tanh() {
2445        return tanh(real, imaginary, Complex::ofCartesian);
2446    }
2447
2448    /**
2449     * Returns the hyperbolic tangent of this complex number.
2450     *
2451     * <p>This function exists to allow implementation of the identity
2452     * {@code tan(z) = -i tanh(iz)}.<p>
2453     *
2454     * @param real Real part.
2455     * @param imaginary Imaginary part.
2456     * @param constructor Constructor.
2457     * @return The hyperbolic tangent of the complex number.
2458     */
2459    private static Complex tanh(double real, double imaginary, ComplexConstructor constructor) {
2460        // Cache the absolute real value
2461        final double x = Math.abs(real);
2462
2463        // Handle inf or nan.
2464        if (!isPosFinite(x) || !Double.isFinite(imaginary)) {
2465            if (isPosInfinite(x)) {
2466                if (Double.isFinite(imaginary)) {
2467                    // The sign is copied from sin(2y)
2468                    // The identity sin(2a) = 2 sin(a) cos(a) is used for consistency
2469                    // with the computation below. Only the magnitude is important
2470                    // so drop the 2. When |y| is small sign(sin(2y)) = sign(y).
2471                    final double sign = Math.abs(imaginary) < PI_OVER_2 ?
2472                                        imaginary :
2473                                        Math.sin(imaginary) * Math.cos(imaginary);
2474                    return constructor.create(Math.copySign(1, real),
2475                                              Math.copySign(0, sign));
2476                }
2477                // imaginary is infinite or NaN
2478                return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
2479            }
2480            // Remaining cases:
2481            // (0 + i inf), returns (0 + i NaN)
2482            // (0 + i NaN), returns (0 + i NaN)
2483            // (x + i inf), returns (NaN + i NaN) for non-zero x (including infinite)
2484            // (x + i NaN), returns (NaN + i NaN) for non-zero x (including infinite)
2485            // (NaN + i 0), returns (NaN + i 0)
2486            // (NaN + i y), returns (NaN + i NaN) for non-zero y (including infinite)
2487            // (NaN + i NaN), returns (NaN + i NaN)
2488            return constructor.create(real == 0 ? real : Double.NaN,
2489                                      imaginary == 0 ? imaginary : Double.NaN);
2490        }
2491
2492        // Finite components
2493        // tanh(x+iy) = (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y))
2494
2495        if (real == 0) {
2496            // Imaginary-only tanh(iy) = i tan(y)
2497            // Identity: sin 2y / (1 + cos 2y) = tan(y)
2498            return constructor.create(real, Math.tan(imaginary));
2499        }
2500        if (imaginary == 0) {
2501            // Identity: sinh 2x / (1 + cosh 2x) = tanh(x)
2502            return constructor.create(Math.tanh(real), imaginary);
2503        }
2504
2505        // The double angles can be avoided using the identities:
2506        // sinh(2x) = 2 sinh(x) cosh(x)
2507        // sin(2y) = 2 sin(y) cos(y)
2508        // cosh(2x) = 2 sinh^2(x) + 1
2509        // cos(2y) = 2 cos^2(y) - 1
2510        // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
2511        // To avoid a junction when swapping between the double angles and the identities
2512        // the identities are used in all cases.
2513
2514        if (x > SAFE_EXP / 2) {
2515            // Potential overflow in sinh/cosh(2x).
2516            // Approximate sinh/cosh using exp^x.
2517            // Ignore cos^2(y) in the divisor as it is insignificant.
2518            // real = sinh(x)cosh(x) / sinh^2(x) = +/-1
2519            final double re = Math.copySign(1, real);
2520            // imag = sin(2y) / 2 sinh^2(x)
2521            // sinh(x) -> sign(x) * e^|x| / 2 when x is large.
2522            // sinh^2(x) -> e^2|x| / 4 when x is large.
2523            // imag = sin(2y) / 2 (e^2|x| / 4) = 2 sin(2y) / e^2|x|
2524            //      = 4 * sin(y) cos(y) / e^2|x|
2525            // Underflow safe divide as e^2|x| may overflow:
2526            // imag = 4 * sin(y) cos(y) / e^m / e^(2|x| - m)
2527            // (|im| is a maximum of 2)
2528            double im = Math.sin(imaginary) * Math.cos(imaginary);
2529            if (x > SAFE_EXP) {
2530                // e^2|x| > e^m * e^m
2531                // This will underflow 2.0 / e^m / e^m
2532                im = Math.copySign(0.0, im);
2533            } else {
2534                // e^2|x| = e^m * e^(2|x| - m)
2535                im = 4 * im / EXP_M / Math.exp(2 * x - SAFE_EXP);
2536            }
2537            return constructor.create(re, im);
2538        }
2539
2540        // No overflow of sinh(2x) and cosh(2x)
2541
2542        // Note: This does not use the definitional formula but uses the identity:
2543        // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
2544
2545        final double sinhx = Math.sinh(real);
2546        final double coshx = Math.cosh(real);
2547        final double siny = Math.sin(imaginary);
2548        final double cosy = Math.cos(imaginary);
2549        final double divisor = sinhx * sinhx + cosy * cosy;
2550        return constructor.create(sinhx * coshx / divisor,
2551                                  siny * cosy / divisor);
2552    }
2553
2554    /**
2555     * Returns the
2556     * <a href="http://mathworld.wolfram.com/InverseHyperbolicSine.html">
2557     * inverse hyperbolic sine</a> of this complex number.
2558     *
2559     * <p>\[ \sinh^{-1}(z) = \ln \left(z + \sqrt{1 + z^2} \right) \]
2560     *
2561     * <p>The inverse hyperbolic sine of \( z \) is unbounded along the real axis and
2562     * in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
2563     *
2564     * <ul>
2565     * <li>{@code z.conj().asinh() == z.asinh().conj()}.
2566     * <li>This is an odd function: \( \sinh^{-1}(z) = -\sinh^{-1}(-z) \).
2567     * <li>If {@code z} is +0 + i0, returns 0 + i0.
2568     * <li>If {@code z} is x + i∞ for positive-signed finite x, returns +∞ + iπ/2.
2569     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
2570     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
2571     * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
2572     * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
2573     * <li>If {@code z} is NaN + i0, returns NaN + i0.
2574     * <li>If {@code z} is NaN + iy for finite nonzero y, returns NaN + iNaN ("invalid" floating-point operation).
2575     * <li>If {@code z} is NaN + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
2576     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2577     * </ul>
2578     *
2579     * <p>The inverse hyperbolic sine is a multivalued function and requires a branch cut in
2580     * the complex plane; the cut is conventionally placed at the line segments
2581     * \( (-i \infty,-i) \) and \( (i,i \infty) \) of the imaginary axis.
2582     *
2583     * <p>This function is computed using the trigonomic identity:
2584     *
2585     * <p>\[ \sinh^{-1}(z) = -i \sin^{-1}(iz) \]
2586     *
2587     * @return The inverse hyperbolic sine of this complex number.
2588     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSinh/">ArcSinh</a>
2589     */
2590    public Complex asinh() {
2591        // Define in terms of asin
2592        // asinh(z) = -i asin(iz)
2593        // Note: This is the opposite to the identity defined in the C99 standard:
2594        // asin(z) = -i asinh(iz)
2595        // Multiply this number by I, compute asin, then multiply by back
2596        return asin(-imaginary, real, Complex::multiplyNegativeI);
2597    }
2598
2599    /**
2600     * Returns the
2601     * <a href="http://mathworld.wolfram.com/InverseHyperbolicCosine.html">
2602     * inverse hyperbolic cosine</a> of this complex number.
2603     *
2604     * <p>\[ \cosh^{-1}(z) = \ln \left(z + \sqrt{z + 1} \sqrt{z - 1} \right) \]
2605     *
2606     * <p>The inverse hyperbolic cosine of \( z \) is in the range \( [0, \infty) \) along the
2607     * real axis and in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
2608     *
2609     * <ul>
2610     * <li>{@code z.conj().acosh() == z.acosh().conj()}.
2611     * <li>If {@code z} is ±0 + i0, returns +0 + iπ/2.
2612     * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
2613     * <li>If {@code z} is 0 + iNaN, returns NaN + iπ/2 <sup>[1]</sup>.
2614     * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
2615     * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns +∞ + iπ.
2616     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
2617     * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
2618     * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
2619     * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
2620     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
2621     * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
2622     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2623     * </ul>
2624     *
2625     * <p>Special cases include the technical corrigendum
2626     * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
2627     * DR 471: Complex math functions cacosh and ctanh</a>.
2628     *
2629     * <p>The inverse hyperbolic cosine is a multivalued function and requires a branch cut in
2630     * the complex plane; the cut is conventionally placed at the line segment
2631     * \( (-\infty,-1) \) of the real axis.
2632     *
2633     * <p>This function is computed using the trigonomic identity:
2634     *
2635     * <p>\[ \cosh^{-1}(z) = \pm i \cos^{-1}(z) \]
2636     *
2637     * <p>The sign of the multiplier is chosen to give {@code z.acosh().real() >= 0}
2638     * and compatibility with the C99 standard.
2639     *
2640     * @return The inverse hyperbolic cosine of this complex number.
2641     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCosh/">ArcCosh</a>
2642     */
2643    public Complex acosh() {
2644        // Define in terms of acos
2645        // acosh(z) = +-i acos(z)
2646        // Note the special case:
2647        // acos(+-0 + iNaN) = π/2 + iNaN
2648        // acosh(0 + iNaN) = NaN + iπ/2
2649        // will not appropriately multiply by I to maintain positive imaginary if
2650        // acos() imaginary computes as NaN. So do this explicitly.
2651        if (Double.isNaN(imaginary) && real == 0) {
2652            return new Complex(Double.NaN, PI_OVER_2);
2653        }
2654        return acos(real, imaginary, (re, im) ->
2655            // Set the sign appropriately for real >= 0
2656            negative(im) ?
2657                // Multiply by I
2658                new Complex(-im, re) :
2659                // Multiply by -I
2660                new Complex(im, -re)
2661        );
2662    }
2663
2664    /**
2665     * Returns the
2666     * <a href="http://mathworld.wolfram.com/InverseHyperbolicTangent.html">
2667     * inverse hyperbolic tangent</a> of this complex number.
2668     *
2669     * <p>\[ \tanh^{-1}(z) = \frac{1}{2} \ln \left( \frac{1 + z}{1 - z} \right) \]
2670     *
2671     * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the real axis and
2672     * in the range \( [-\pi/2, \pi/2] \) along the imaginary axis. Special cases:
2673     *
2674     * <ul>
2675     * <li>{@code z.conj().atanh() == z.atanh().conj()}.
2676     * <li>This is an odd function: \( \tanh^{-1}(z) = -\tanh^{-1}(-z) \).
2677     * <li>If {@code z} is +0 + i0, returns +0 + i0.
2678     * <li>If {@code z} is +0 + iNaN, returns +0 + iNaN.
2679     * <li>If {@code z} is +1 + i0, returns +∞ + i0 ("divide-by-zero" floating-point operation).
2680     * <li>If {@code z} is x + i∞ for finite positive-signed x, returns +0 + iπ/2.
2681     * <li>If {@code z} is x+iNaN for nonzero finite x, returns NaN+iNaN ("invalid" floating-point operation).
2682     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +0 + iπ/2.
2683     * <li>If {@code z} is +∞ + i∞, returns +0 + iπ/2.
2684     * <li>If {@code z} is +∞ + iNaN, returns +0 + iNaN.
2685     * <li>If {@code z} is NaN+iy for finite y, returns NaN+iNaN ("invalid" floating-point operation).
2686     * <li>If {@code z} is NaN + i∞, returns ±0 + iπ/2 (where the sign of the real part of the result is unspecified).
2687     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
2688     * </ul>
2689     *
2690     * <p>The inverse hyperbolic tangent is a multivalued function and requires a branch cut in
2691     * the complex plane; the cut is conventionally placed at the line segments
2692     * \( (\infty,-1] \) and \( [1,\infty) \) of the real axis.
2693     *
2694     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
2695     *
2696     * <p>\[ \tanh^{-1}(z) = \frac{1}{4} \ln \left(1 + \frac{4x}{(1-x)^2+y^2} \right) + \\
2697     *                     i \frac{1}{2} \left( \tan^{-1} \left(\frac{2y}{1-x^2-y^2} \right) + \frac{\pi}{2} \left(\text{sgn}(x^2+y^2-1)+1 \right) \text{sgn}(y) \right) \]
2698     *
2699     * <p>The imaginary part is computed using {@link Math#atan2(double, double)} to ensure the
2700     * correct quadrant is returned from \( \tan^{-1} \left(\frac{2y}{1-x^2-y^2} \right) \).
2701     *
2702     * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
2703     * {@code c++} implementation {@code <boost/math/complex/atanh.hpp>}.
2704     *
2705     * @return The inverse hyperbolic tangent of this complex number.
2706     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTanh/">ArcTanh</a>
2707     */
2708    public Complex atanh() {
2709        return atanh(real, imaginary, Complex::ofCartesian);
2710    }
2711
2712    /**
2713     * Returns the inverse hyperbolic tangent of this complex number.
2714     *
2715     * <p>This function exists to allow implementation of the identity
2716     * {@code atan(z) = -i atanh(iz)}.
2717     *
2718     * <p>Adapted from {@code <boost/math/complex/atanh.hpp>}. This method only (and not
2719     * invoked methods within) is distributed under the Boost Software License V1.0.
2720     * The original notice is shown below and the licence is shown in full in LICENSE:
2721     * <pre>
2722     * (C) Copyright John Maddock 2005.
2723     * Distributed under the Boost Software License, Version 1.0. (See accompanying
2724     * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
2725     * </pre>
2726     *
2727     * @param real Real part.
2728     * @param imaginary Imaginary part.
2729     * @param constructor Constructor.
2730     * @return The inverse hyperbolic tangent of the complex number.
2731     */
2732    private static Complex atanh(final double real, final double imaginary,
2733                                 final ComplexConstructor constructor) {
2734        // Compute with positive values and determine sign at the end
2735        double x = Math.abs(real);
2736        double y = Math.abs(imaginary);
2737        // The result (without sign correction)
2738        double re;
2739        double im;
2740
2741        // Handle C99 special cases
2742        if (Double.isNaN(x)) {
2743            if (isPosInfinite(y)) {
2744                // The sign of the real part of the result is unspecified
2745                return constructor.create(0, Math.copySign(PI_OVER_2, imaginary));
2746            }
2747            // No-use of the input constructor.
2748            // Optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
2749            return NAN;
2750        } else if (Double.isNaN(y)) {
2751            if (isPosInfinite(x)) {
2752                return constructor.create(Math.copySign(0, real), Double.NaN);
2753            }
2754            if (x == 0) {
2755                return constructor.create(real, Double.NaN);
2756            }
2757            // No-use of the input constructor
2758            return NAN;
2759        } else {
2760            // x && y are finite or infinite.
2761
2762            // Check the safe region.
2763            // The lower and upper bounds have been copied from boost::math::atanh.
2764            // They are different from the safe region for asin and acos.
2765            // x >= SAFE_UPPER: (1-x) == -x
2766            // x <= SAFE_LOWER: 1 - x^2 = 1
2767
2768            if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
2769                // Normal computation within a safe region.
2770
2771                // minus x plus 1: (-x+1)
2772                final double mxp1 = 1 - x;
2773                final double yy = y * y;
2774                // The definition of real component is:
2775                // real = log( ((x+1)^2+y^2) / ((1-x)^2+y^2) ) / 4
2776                // This simplifies by adding 1 and subtracting 1 as a fraction:
2777                //      = log(1 + ((x+1)^2+y^2) / ((1-x)^2+y^2) - ((1-x)^2+y^2)/((1-x)^2+y^2) ) / 4
2778                //
2779                // real(atanh(z)) == log(1 + 4*x / ((1-x)^2+y^2)) / 4
2780                // imag(atanh(z)) == tan^-1 (2y, (1-x)(1+x) - y^2) / 2
2781                // imag(atanh(z)) == tan^-1 (2y, (1 - x^2 - y^2) / 2
2782                // The division is done at the end of the function.
2783                re = Math.log1p(4 * x / (mxp1 * mxp1 + yy));
2784                // Modified from boost which does not switch the magnitude of x and y.
2785                // The denominator for atan2 is 1 - x^2 - y^2.
2786                // This can be made more precise if |x| > |y|.
2787                final double numerator = 2 * y;
2788                final double denominator;
2789                if (x < y) {
2790                    final double tmp = x;
2791                    x = y;
2792                    y = tmp;
2793                }
2794                // 1 - x is precise if |x| >= 1
2795                if (x >= 1) {
2796                    denominator = (1 - x) * (1 + x) - y * y;
2797                } else {
2798                    // |x| < 1: Use high precision if possible:
2799                    // 1 - x^2 - y^2 = -(x^2 + y^2 - 1)
2800                    // Modified from boost to use the custom high precision method.
2801                    denominator = -x2y2m1(x, y);
2802                }
2803                im = Math.atan2(numerator, denominator);
2804            } else {
2805                // This section handles exception cases that would normally cause
2806                // underflow or overflow in the main formulas.
2807
2808                // C99. G.7: Special case for imaginary only numbers
2809                if (x == 0) {
2810                    if (imaginary == 0) {
2811                        return constructor.create(real, imaginary);
2812                    }
2813                    // atanh(iy) = i atan(y)
2814                    return constructor.create(real, Math.atan(imaginary));
2815                }
2816
2817                // Real part:
2818                // real = Math.log1p(4x / ((1-x)^2 + y^2))
2819                // real = Math.log1p(4x / (1 - 2x + x^2 + y^2))
2820                // real = Math.log1p(4x / (1 + x(x-2) + y^2))
2821                // without either overflow or underflow in the squared terms.
2822                if (x >= SAFE_UPPER) {
2823                    // (1-x) = -x to machine precision:
2824                    // log1p(4x / (x^2 + y^2))
2825                    if (isPosInfinite(x) || isPosInfinite(y)) {
2826                        re = 0;
2827                    } else if (y >= SAFE_UPPER) {
2828                        // Big x and y: divide by x*y
2829                        re = Math.log1p((4 / y) / (x / y + y / x));
2830                    } else if (y > 1) {
2831                        // Big x: divide through by x:
2832                        re = Math.log1p(4 / (x + y * y / x));
2833                    } else {
2834                        // Big x small y, as above but neglect y^2/x:
2835                        re = Math.log1p(4 / x);
2836                    }
2837                } else if (y >= SAFE_UPPER) {
2838                    if (x > 1) {
2839                        // Big y, medium x, divide through by y:
2840                        final double mxp1 = 1 - x;
2841                        re = Math.log1p((4 * x / y) / (mxp1 * mxp1 / y + y));
2842                    } else {
2843                        // Big y, small x, as above but neglect (1-x)^2/y:
2844                        // Note: log1p(v) == v - v^2/2 + v^3/3 ... Taylor series when v is small.
2845                        // Here v is so small only the first term matters.
2846                        re = 4 * x / y / y;
2847                    }
2848                } else if (x == 1) {
2849                    // x = 1, small y:
2850                    // Special case when x == 1 as (1-x) is invalid.
2851                    // Simplify the following formula:
2852                    // real = log( sqrt((x+1)^2+y^2) ) / 2 - log( sqrt((1-x)^2+y^2) ) / 2
2853                    //      = log( sqrt(4+y^2) ) / 2 - log(y) / 2
2854                    // if: 4+y^2 -> 4
2855                    //      = log( 2 ) / 2 - log(y) / 2
2856                    //      = (log(2) - log(y)) / 2
2857                    // Multiply by 2 as it will be divided by 4 at the end.
2858                    // C99: if y=0 raises the ‘‘divide-by-zero’’ floating-point exception.
2859                    re = 2 * (LN_2 - Math.log(y));
2860                } else {
2861                    // Modified from boost which checks y > SAFE_LOWER.
2862                    // if y*y -> 0 it will be ignored so always include it.
2863                    final double mxp1 = 1 - x;
2864                    re = Math.log1p((4 * x) / (mxp1 * mxp1 + y * y));
2865                }
2866
2867                // Imaginary part:
2868                // imag = atan2(2y, (1-x)(1+x) - y^2)
2869                // if x or y are large, then the formula:
2870                //   atan2(2y, (1-x)(1+x) - y^2)
2871                // evaluates to +(PI - theta) where theta is negligible compared to PI.
2872                if (x >= SAFE_UPPER || y >= SAFE_UPPER) {
2873                    im = Math.PI;
2874                } else if (x <= SAFE_LOWER) {
2875                    // (1-x)^2 -> 1
2876                    if (y <= SAFE_LOWER) {
2877                        // 1 - y^2 -> 1
2878                        im = Math.atan2(2 * y, 1);
2879                    } else {
2880                        im = Math.atan2(2 * y, 1 - y * y);
2881                    }
2882                } else {
2883                    // Medium x, small y.
2884                    // Modified from boost which checks (y == 0) && (x == 1) and sets re = 0.
2885                    // This is same as the result from calling atan2(0, 0) so exclude this case.
2886                    // 1 - y^2 = 1 so ignore subtracting y^2
2887                    im = Math.atan2(2 * y, (1 - x) * (1 + x));
2888                }
2889            }
2890        }
2891
2892        re /= 4;
2893        im /= 2;
2894        return constructor.create(changeSign(re, real),
2895                                  changeSign(im, imaginary));
2896    }
2897
2898    /**
2899     * Compute {@code x^2 + y^2 - 1} in high precision.
2900     * Assumes that the values x and y can be multiplied without overflow; that
2901     * {@code x >= y}; and both values are positive.
2902     *
2903     * @param x the x value
2904     * @param y the y value
2905     * @return {@code x^2 + y^2 - 1}.
2906     */
2907    private static double x2y2m1(double x, double y) {
2908        // Hull et al used (x-1)*(x+1)+y*y.
2909        // From the paper on page 236:
2910
2911        // If x == 1 there is no cancellation.
2912
2913        // If x > 1, there is also no cancellation, but the argument is now accurate
2914        // only to within a factor of 1 + 3 EPSILSON (note that x – 1 is exact),
2915        // so that error = 3 EPSILON.
2916
2917        // If x < 1, there can be serious cancellation:
2918
2919        // If 4 y^2 < |x^2 – 1| the cancellation is not serious ... the argument is accurate
2920        // only to within a factor of 1 + 4 EPSILSON so that error = 4 EPSILON.
2921
2922        // Otherwise there can be serious cancellation and the relative error in the real part
2923        // could be enormous.
2924
2925        final double xx = x * x;
2926        final double yy = y * y;
2927        // Modify to use high precision before the threshold set by Hull et al.
2928        // This is to preserve the monotonic output of the computation at the switch.
2929        // Set the threshold when x^2 + y^2 is above 0.5 thus subtracting 1 results in a number
2930        // that can be expressed with a higher precision than any number in the range 0.5-1.0
2931        // due to the variable exponent used below 0.5.
2932        if (x < 1 && xx + yy > 0.5) {
2933            // Large relative error.
2934            // This does not use o.a.c.numbers.LinearCombination.value(x, x, y, y, 1, -1).
2935            // It is optimised knowing that:
2936            // - the products are squares
2937            // - the final term is -1 (which does not require split multiplication and addition)
2938            // - The answer will not be NaN as the terms are not NaN components
2939            // - The order is known to be 1 > |x| >= |y|
2940            // The squares are computed using a split multiply algorithm and
2941            // the summation using an extended precision summation algorithm.
2942
2943            // Split x and y as one 26 bits number and one 27 bits number
2944            final double xHigh = splitHigh(x);
2945            final double xLow  = x - xHigh;
2946            final double yHigh = splitHigh(y);
2947            final double yLow  = y - yHigh;
2948
2949            // Accurate split multiplication x * x and y * y
2950            final double x2Low = squareLow(xLow, xHigh, xx);
2951            final double y2Low = squareLow(yLow, yHigh, yy);
2952
2953            return sumx2y2m1(xx, x2Low, yy, y2Low);
2954        }
2955        return (x - 1) * (x + 1) + yy;
2956    }
2957
2958    /**
2959     * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) create
2960     * a big value from which to derive the two split parts.
2961     * <pre>
2962     * c = (2^s + 1) * a
2963     * a_big = c - a
2964     * a_hi = c - a_big
2965     * a_lo = a - a_hi
2966     * a = a_hi + a_lo
2967     * </pre>
2968     *
2969     * <p>The multiplicand must be odd allowing a p-bit value to be split into
2970     * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.
2971     * Combined they have (p􏰔-1) bits of significand but the sign bit of {@code a_lo}
2972     * contains a bit of information.
2973     *
2974     * @param a Value.
2975     * @return the high part of the value.
2976     * @see <a href="https://doi.org/10.1007/BF01397083">
2977     * Dekker (1971) A floating-point technique for extending the available precision</a>
2978     */
2979    private static double splitHigh(double a) {
2980        final double c = MULTIPLIER * a;
2981        return c - (c - a);
2982    }
2983
2984    /**
2985     * Compute the round-off from the square of a split number with {@code low} and {@code high}
2986     * components. Uses Dekker's algorithm for split multiplication modified for a square product.
2987     *
2988     * <p>Note: This is candidate to be replaced with {@code Math.fma(x, x, -x * x)} to compute
2989     * the round-off from the square product {@code x * x}. This would remove the requirement
2990     * to compute the split number and make this method redundant. {@code Math.fma} requires
2991     * JDK 9 and FMA hardware support.
2992     *
2993     * @param low Low part of number.
2994     * @param high High part of number.
2995     * @param square Square of the number.
2996     * @return <code>low * low - (((product - high * high) - low * high) - high * low)</code>
2997     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
2998     * Shewchuk (1997) Theorum 18</a>
2999     */
3000    private static double squareLow(double low, double high, double square) {
3001        final double lh = low * high;
3002        return low * low - (((square - high * high) - lh) - lh);
3003    }
3004
3005    /**
3006     * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
3007     * Dekker's two-sum algorithm. The values are required to be ordered by magnitude:
3008     * {@code |a| >= |b|}.
3009     *
3010     * @param a First part of sum.
3011     * @param b Second part of sum.
3012     * @param x Sum.
3013     * @return <code>b - (x - a)</code>
3014     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
3015     * Shewchuk (1997) Theorum 6</a>
3016     */
3017    private static double fastSumLow(double a, double b, double x) {
3018        // x = a + b
3019        // bVirtual = x - a
3020        // y = b - bVirtual
3021        return b - (x - a);
3022    }
3023
3024    /**
3025     * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
3026     * Knuth's two-sum algorithm. The values are not required to be ordered by magnitude.
3027     *
3028     * @param a First part of sum.
3029     * @param b Second part of sum.
3030     * @param x Sum.
3031     * @return <code>(a - (x - (x - a))) + (b - (x - a))</code>
3032     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
3033     * Shewchuk (1997) Theorum 7</a>
3034     */
3035    private static double sumLow(double a, double b, double x) {
3036        // x = a + b
3037        // bVirtual = x - a
3038        // aVirtual = x - bVirtual
3039        // bRoundoff = b - bVirtual
3040        // aRoundoff = a - aVirtual
3041        // y = aRoundoff + bRoundoff
3042        final double bVirtual = x - a;
3043        return (a - (x - bVirtual)) + (b - bVirtual);
3044    }
3045
3046    /**
3047     * Sum x^2 + y^2 - 1. It is assumed that {@code y <= x < 1}.
3048     *
3049     * <p>Implement Shewchuk's expansion-sum algorithm: [x2Low, x2High] + [-1] + [y2Low, y2High].
3050     *
3051     * @param x2High High part of x^2.
3052     * @param x2Low Low part of x^2.
3053     * @param y2High High part of y^2.
3054     * @param y2Low Low part of y^2.
3055     * @return x^2 + y^2 - 1
3056     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
3057     * Shewchuk (1997) Theorum 12</a>
3058     */
3059    private static double sumx2y2m1(double x2High, double x2Low, double y2High, double y2Low) {
3060        // Let e and f be non-overlapping expansions of components of length m and n.
3061        // The following algorithm will produce a non-overlapping expansion h where the
3062        // sum h_i = e + f and components of h are in increasing order of magnitude.
3063
3064        // Expansion-sum proceeds by a grow-expansion of the first part from one expansion
3065        // into the other, extending its length by 1. The process repeats for the next part
3066        // but the grow-expansion starts at the previous merge position + 1.
3067        // Thus expansion-sum requires mn two-sum operations to merge length m into length n
3068        // resulting in length m+n-1.
3069
3070        // Variables numbered from 1 as per Figure 7 (p.12). The output expansion h is placed
3071        // into e increasing its length for each grow expansion.
3072
3073        // We have two expansions for x^2 and y^2 and the whole number -1.
3074        // Expecting (x^2 + y^2) close to 1 we generate first the intermediate expansion
3075        // (x^2 - 1) moving the result away from 1 where there are sparse floating point
3076        // representations. This is then added to a similar magnitude y^2. Leaving the -1
3077        // until last suffers from 1 ulp rounding errors more often and the requirement
3078        // for a distillation sum to reduce rounding error frequency.
3079
3080        // Note: Do not use the alternative fast-expansion-sum of the parts sorted by magnitude.
3081        // The parts can be ordered with a single comparison into:
3082        // [y2Low, (y2High|x2Low), x2High, -1]
3083        // The fast-two-sum saves 1 fast-two-sum and 3 two-sum operations (21 additions) and
3084        // adds a penalty of a single branch condition.
3085        // However the order in not "strongly non-overlapping" and the fast-expansion-sum
3086        // output will not be strongly non-overlapping. The sum of the output has 1 ulp error
3087        // on random cis numbers approximately 1 in 160 events. This can be removed by a
3088        // distillation two-sum pass over the final expansion as a cost of 1 fast-two-sum and
3089        // 3 two-sum operations! So we use the expansion sum with the same operations and
3090        // no branches.
3091
3092        // q=running sum
3093        double q = x2Low - 1;
3094        double e1 = fastSumLow(-1, x2Low, q);
3095        double e3 = q + x2High;
3096        double e2 = sumLow(q, x2High, e3);
3097
3098        final double f1 = y2Low;
3099        final double f2 = y2High;
3100
3101        // Grow expansion of f1 into e
3102        q = f1 + e1;
3103        e1 = sumLow(f1, e1, q);
3104        double p = q + e2;
3105        e2 = sumLow(q, e2, p);
3106        double e4 = p + e3;
3107        e3 = sumLow(p, e3, e4);
3108
3109        // Grow expansion of f2 into e (only required to start at e2)
3110        q = f2 + e2;
3111        e2 = sumLow(f2, e2, q);
3112        p = q + e3;
3113        e3 = sumLow(q, e3, p);
3114        final double e5 = p + e4;
3115        e4 = sumLow(p, e4, e5);
3116
3117        // Final summation:
3118        // The sum of the parts is within 1 ulp of the true expansion value e:
3119        // |e - sum| < ulp(sum).
3120        // To achieve the exact result requires iteration of a distillation two-sum through
3121        // the expansion until convergence, i.e. no smaller term changes higher terms.
3122        // This requires (n-1) iterations for length n. Here we neglect this as
3123        // although the method is not ensured to be exact is it robust on random
3124        // cis numbers.
3125        return e1 + e2 + e3 + e4 + e5;
3126    }
3127
3128    /**
3129     * Returns the n-th roots of this complex number.
3130     * The nth roots are defined by the formula:
3131     *
3132     * <p>\[ z_k = |z|^{\frac{1}{n}} \left( \cos \left(\phi + \frac{2\pi k}{n} \right) + i \sin \left(\phi + \frac{2\pi k}{n} \right) \right) \]
3133     *
3134     * <p>for \( k=0, 1, \ldots, n-1 \), where \( |z| \) and \( \phi \)
3135     * are respectively the {@link #abs() modulus} and
3136     * {@link #arg() argument} of this complex number.
3137     *
3138     * <p>If one or both parts of this complex number is NaN, a list with all
3139     * all elements set to {@code NaN + i NaN} is returned.</p>
3140     *
3141     * @param n Degree of root.
3142     * @return A list of all {@code n}-th roots of this complex number.
3143     * @throws IllegalArgumentException if {@code n} is zero.
3144     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Root/">Root</a>
3145     */
3146    public List<Complex> nthRoot(int n) {
3147        if (n == 0) {
3148            throw new IllegalArgumentException("cannot compute zeroth root");
3149        }
3150
3151        final List<Complex> result = new ArrayList<>();
3152
3153        // nth root of abs -- faster / more accurate to use a solver here?
3154        final double nthRootOfAbs = Math.pow(abs(), 1.0 / n);
3155
3156        // Compute nth roots of complex number with k = 0, 1, ... n-1
3157        final double nthPhi = arg() / n;
3158        final double slice = 2 * Math.PI / n;
3159        double innerPart = nthPhi;
3160        for (int k = 0; k < Math.abs(n); k++) {
3161            // inner part
3162            final double realPart = nthRootOfAbs *  Math.cos(innerPart);
3163            final double imaginaryPart = nthRootOfAbs *  Math.sin(innerPart);
3164            result.add(ofCartesian(realPart, imaginaryPart));
3165            innerPart += slice;
3166        }
3167
3168        return result;
3169    }
3170
3171    /**
3172     * Test for equality with another object. If the other object is a {@code Complex} then a
3173     * comparison is made of the real and imaginary parts; otherwise {@code false} is returned.
3174     *
3175     * <p>If both the real and imaginary parts of two complex numbers
3176     * are exactly the same the two {@code Complex} objects are considered to be equal.
3177     * For this purpose, two {@code double} values are considered to be
3178     * the same if and only if the method {@link Double #doubleToLongBits(double)}
3179     * returns the identical {@code long} value when applied to each.
3180     *
3181     * <p>Note that in most cases, for two instances of class
3182     * {@code Complex}, {@code c1} and {@code c2}, the
3183     * value of {@code c1.equals(c2)} is {@code true} if and only if
3184     *
3185     * <pre>
3186     *  {@code c1.getReal() == c2.getReal() && c1.getImaginary() == c2.getImaginary()}</pre>
3187     *
3188     * <p>also has the value {@code true}. However, there are exceptions:
3189     *
3190     * <ul>
3191     *  <li>
3192     *   Instances that contain {@code NaN} values in the same part
3193     *   are considered to be equal for that part, even though {@code Double.NaN == Double.NaN}
3194     *   has the value {@code false}.
3195     *  </li>
3196     *  <li>
3197     *   Instances that share a {@code NaN} value in one part
3198     *   but have different values in the other part are <em>not</em> considered equal.
3199     *  </li>
3200     *  <li>
3201     *   Instances that contain different representations of zero in the same part
3202     *   are <em>not</em> considered to be equal for that part, even though {@code -0.0 == 0.0}
3203     *   has the value {@code true}.
3204     *  </li>
3205     * </ul>
3206     *
3207     * <p>The behavior is the same as if the components of the two complex numbers were passed
3208     * to {@link java.util.Arrays#equals(double[], double[]) Arrays.equals(double[], double[])}:
3209     *
3210     * <pre>
3211     *  Arrays.equals(new double[]{c1.getReal(), c1.getImaginary()},
3212     *                new double[]{c2.getReal(), c2.getImaginary()}); </pre>
3213     *
3214     * @param other Object to test for equality with this instance.
3215     * @return {@code true} if the objects are equal, {@code false} if object
3216     * is {@code null}, not an instance of {@code Complex}, or not equal to
3217     * this instance.
3218     * @see Double#doubleToLongBits(double)
3219     * @see java.util.Arrays#equals(double[], double[])
3220     */
3221    @Override
3222    public boolean equals(Object other) {
3223        if (this == other) {
3224            return true;
3225        }
3226        if (other instanceof Complex) {
3227            final Complex c = (Complex) other;
3228            return equals(real, c.real) &&
3229                equals(imaginary, c.imaginary);
3230        }
3231        return false;
3232    }
3233
3234    /**
3235     * Gets a hash code for the complex number.
3236     *
3237     * <p>The behavior is the same as if the components of the complex number were passed
3238     * to {@link java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])}:
3239     *
3240     * <pre>
3241     *  {@code Arrays.hashCode(new double[] {getReal(), getImaginary()})}</pre>
3242     *
3243     * @return A hash code value for this object.
3244     * @see java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])
3245     */
3246    @Override
3247    public int hashCode() {
3248        return 31 * (31 + Double.hashCode(real)) + Double.hashCode(imaginary);
3249    }
3250
3251    /**
3252     * Returns a string representation of the complex number.
3253     *
3254     * <p>The string will represent the numeric values of the real and imaginary parts.
3255     * The values are split by a separator and surrounded by parentheses.
3256     * The string can be {@link #parse(String) parsed} to obtain an instance with the same value.
3257     *
3258     * <p>The format for complex number \( x + i y \) is {@code "(x,y)"}, with \( x \) and
3259     * \( y \) converted as if using {@link Double#toString(double)}.
3260     *
3261     * @return A string representation of the complex number.
3262     * @see #parse(String)
3263     * @see Double#toString(double)
3264     */
3265    @Override
3266    public String toString() {
3267        return new StringBuilder(TO_STRING_SIZE)
3268            .append(FORMAT_START)
3269            .append(real).append(FORMAT_SEP)
3270            .append(imaginary)
3271            .append(FORMAT_END)
3272            .toString();
3273    }
3274
3275    /**
3276     * Returns {@code true} if the values are equal according to semantics of
3277     * {@link Double#equals(Object)}.
3278     *
3279     * @param x Value
3280     * @param y Value
3281     * @return {@code Double.valueof(x).equals(Double.valueOf(y))}.
3282     */
3283    private static boolean equals(double x, double y) {
3284        return Double.doubleToLongBits(x) == Double.doubleToLongBits(y);
3285    }
3286
3287    /**
3288     * Check that a value is negative. It must meet all the following conditions:
3289     * <ul>
3290     *  <li>it is not {@code NaN},</li>
3291     *  <li>it is negative signed,</li>
3292     * </ul>
3293     *
3294     * <p>Note: This is true for negative zero.</p>
3295     *
3296     * @param d Value.
3297     * @return {@code true} if {@code d} is negative.
3298     */
3299    private static boolean negative(double d) {
3300        return d < 0 || Double.doubleToLongBits(d) == NEGATIVE_ZERO_LONG_BITS;
3301    }
3302
3303    /**
3304     * Check that a value is positive infinity. Used to replace {@link Double#isInfinite()}
3305     * when the input value is known to be positive (i.e. in the case where it has been
3306     * set using {@link Math#abs(double)}).
3307     *
3308     * @param d Value.
3309     * @return {@code true} if {@code d} is +inf.
3310     */
3311    private static boolean isPosInfinite(double d) {
3312        return d == Double.POSITIVE_INFINITY;
3313    }
3314
3315    /**
3316     * Check that an absolute value is finite. Used to replace {@link Double#isFinite(double)}
3317     * when the input value is known to be positive (i.e. in the case where it has been
3318     * set using {@link Math#abs(double)}).
3319     *
3320     * @param d Value.
3321     * @return {@code true} if {@code d} is +finite.
3322     */
3323    private static boolean isPosFinite(double d) {
3324        return d <= Double.MAX_VALUE;
3325    }
3326
3327    /**
3328     * Create a complex number given the real and imaginary parts, then multiply by {@code -i}.
3329     * This is used in functions that implement trigonomic identities. It is the functional
3330     * equivalent of:
3331     *
3332     * <pre>
3333     *  z = new Complex(real, imaginary).multiplyImaginary(-1);</pre>
3334     *
3335     * @param real Real part.
3336     * @param imaginary Imaginary part.
3337     * @return {@code Complex} object.
3338     */
3339    private static Complex multiplyNegativeI(double real, double imaginary) {
3340        return new Complex(imaginary, -real);
3341    }
3342
3343    /**
3344     * Change the sign of the magnitude based on the signed value.
3345     *
3346     * <p>If the signed value is negative then the result is {@code -magnitude}; otherwise
3347     * return {@code magnitude}.
3348     *
3349     * <p>A signed value of {@code -0.0} is treated as negative. A signed value of {@code NaN}
3350     * is treated as positive.
3351     *
3352     * <p>This is not the same as {@link Math#copySign(double, double)} as this method
3353     * will change the sign based on the signed value rather than copy the sign.
3354     *
3355     * @param magnitude the magnitude
3356     * @param signedValue the signed value
3357     * @return magnitude or -magnitude.
3358     * @see #negative(double)
3359     */
3360    private static double changeSign(double magnitude, double signedValue) {
3361        return negative(signedValue) ? -magnitude : magnitude;
3362    }
3363
3364    /**
3365     * Returns a scale suitable for use with {@link Math#scalb(double, int)} to normalise
3366     * the number to the interval {@code [1, 2)}.
3367     *
3368     * <p>The scale is typically the largest unbiased exponent used in the representation of the
3369     * two numbers. In contrast to {@link Math#getExponent(double)} this handles
3370     * sub-normal numbers by computing the number of leading zeros in the mantissa
3371     * and shifting the unbiased exponent. The result is that for all finite, non-zero,
3372     * numbers {@code a, b}, the magnitude of {@code scalb(x, -getScale(a, b))} is
3373     * always in the range {@code [1, 2)}, where {@code x = max(|a|, |b|)}.
3374     *
3375     * <p>This method is a functional equivalent of the c function ilogb(double) adapted for
3376     * two input arguments.
3377     *
3378     * <p>The result is to be used to scale a complex number using {@link Math#scalb(double, int)}.
3379     * Hence the special case of both zero arguments is handled using the return value for NaN
3380     * as zero cannot be scaled. This is different from {@link Math#getExponent(double)}
3381     * or {@link #getMaxExponent(double, double)}.
3382     *
3383     * <p>Special cases:
3384     *
3385     * <ul>
3386     * <li>If either argument is NaN or infinite, then the result is
3387     * {@link Double#MAX_EXPONENT} + 1.
3388     * <li>If both arguments are zero, then the result is
3389     * {@link Double#MAX_EXPONENT} + 1.
3390     * </ul>
3391     *
3392     * @param a the first value
3393     * @param b the second value
3394     * @return The maximum unbiased exponent of the values to be used for scaling
3395     * @see Math#getExponent(double)
3396     * @see Math#scalb(double, int)
3397     * @see <a href="http://www.cplusplus.com/reference/cmath/ilogb/">ilogb</a>
3398     */
3399    private static int getScale(double a, double b) {
3400        // Only interested in the exponent and mantissa so remove the sign bit
3401        final long x = Double.doubleToRawLongBits(a) & UNSIGN_MASK;
3402        final long y = Double.doubleToRawLongBits(b) & UNSIGN_MASK;
3403        // Only interested in the maximum
3404        final long bits = Math.max(x, y);
3405        // Get the unbiased exponent
3406        int exp = ((int) (bits >>> 52)) - EXPONENT_OFFSET;
3407
3408        // No case to distinguish nan/inf
3409        // Handle sub-normal numbers
3410        if (exp == Double.MIN_EXPONENT - 1) {
3411            // Special case for zero, return as nan/inf to indicate scaling is not possible
3412            if (bits == 0) {
3413                return Double.MAX_EXPONENT + 1;
3414            }
3415            // A sub-normal number has an exponent below -1022. The amount below
3416            // is defined by the number of shifts of the most significant bit in
3417            // the mantissa that is required to get a 1 at position 53 (i.e. as
3418            // if it were a normal number with assumed leading bit)
3419            final long mantissa = bits & MANTISSA_MASK;
3420            exp -= Long.numberOfLeadingZeros(mantissa << 12);
3421        }
3422        return exp;
3423    }
3424
3425    /**
3426     * Returns the largest unbiased exponent used in the representation of the
3427     * two numbers. Special cases:
3428     *
3429     * <ul>
3430     * <li>If either argument is NaN or infinite, then the result is
3431     * {@link Double#MAX_EXPONENT} + 1.
3432     * <li>If both arguments are zero or subnormal, then the result is
3433     * {@link Double#MIN_EXPONENT} -1.
3434     * </ul>
3435     *
3436     * <p>This is used by {@link #divide(double, double, double, double)} as
3437     * a simple detection that a number may overflow if multiplied
3438     * by a value in the interval [1, 2).
3439     *
3440     * @param a the first value
3441     * @param b the second value
3442     * @return The maximum unbiased exponent of the values.
3443     * @see Math#getExponent(double)
3444     * @see #divide(double, double, double, double)
3445     */
3446    private static int getMaxExponent(double a, double b) {
3447        // This could return:
3448        // Math.getExponent(Math.max(Math.abs(a), Math.abs(b)))
3449        // A speed test is required to determine performance.
3450        return Math.max(Math.getExponent(a), Math.getExponent(b));
3451    }
3452
3453    /**
3454     * Checks if both x and y are in the region defined by the minimum and maximum.
3455     *
3456     * @param x x value.
3457     * @param y y value.
3458     * @param min the minimum (exclusive).
3459     * @param max the maximum (exclusive).
3460     * @return true if inside the region.
3461     */
3462    private static boolean inRegion(double x, double y, double min, double max) {
3463        return x < max && x > min && y < max && y > min;
3464    }
3465
3466    /**
3467     * Returns {@code sqrt(x^2 + y^2)} without intermediate overflow or underflow.
3468     *
3469     * <p>Special cases:
3470     * <ul>
3471     * <li>If either argument is infinite, then the result is positive infinity.
3472     * <li>If either argument is NaN and neither argument is infinite, then the result is NaN.
3473     * </ul>
3474     *
3475     * <p>The computed result is expected to be within 1 ulp of the exact result.
3476     *
3477     * <p>This method is a replacement for {@link Math#hypot(double, double)}. There
3478     * will be differences between this method and {@code Math.hypot(double, double)} due
3479     * to the use of a different algorithm to compute the high precision sum of
3480     * {@code x^2 + y^2}. This method has been tested to have a lower maximum error from
3481     * the exact result; any differences are expected to be 1 ULP indicating a rounding
3482     * change in the sum.
3483     *
3484     * <p>JDK9 ported the hypot function to Java for bug JDK-7130085 due to the slow performance
3485     * of the method as a native function. Benchmarks of the Complex class for functions that
3486     * use hypot confirm this is slow pre-Java 9. This implementation outperforms the new faster
3487     * {@code Math.hypot(double, double)} on JDK 11 (LTS). See the Commons numbers examples JMH
3488     * module for benchmarks. Comparisons with alternative implementations indicate
3489     * performance gains are related to edge case handling and elimination of an unpredictable
3490     * branch in the computation of {@code x^2 + y^2}.
3491     *
3492     * <p>This port was adapted from the "Freely Distributable Math Library" hypot function.
3493     * This method only (and not invoked methods within) is distributed under the terms of the
3494     * original notice as shown below:
3495     * <pre>
3496     * ====================================================
3497     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
3498     *
3499     * Developed at SunSoft, a Sun Microsystems, Inc. business.
3500     * Permission to use, copy, modify, and distribute this
3501     * software is freely granted, provided that this notice
3502     * is preserved.
3503     * ====================================================
3504     * </pre>
3505     *
3506     * <p>Note: The fdlibm c code makes use of the language ability to read and write directly
3507     * to the upper and lower 32-bits of the 64-double. The function performs
3508     * checking on the upper 32-bits for the magnitude of the two numbers by accessing
3509     * the exponent and 20 most significant bits of the mantissa. These upper bits
3510     * are manipulated during scaling and then used to perform extended precision
3511     * computation of the sum {@code x^2 + y^2} where the high part of the number has 20-bit
3512     * precision. Manipulation of direct bits has no equivalent in Java
3513     * other than use of {@link Double#doubleToLongBits(double)} and
3514     * {@link Double#longBitsToDouble(long)}. To avoid conversion to and from long and double
3515     * representations this implementation only scales the double representation. The high
3516     * and low parts of a double for the extended precision computation are extracted
3517     * using the method of Dekker (1971) to create two 26-bit numbers. This works for sub-normal
3518     * numbers and reduces the maximum error in comparison to fdlibm hypot which does not
3519     * use a split number algorithm for sub-normal numbers.
3520     *
3521     * @param x Value x
3522     * @param y Value y
3523     * @return sqrt(x^2 + y^2)
3524     * @see Math#hypot(double, double)
3525     * @see <a href="https://www.netlib.org/fdlibm/e_hypot.c">fdlibm e_hypot.c</a>
3526     * @see <a href="https://bugs.java.com/bugdatabase/view_bug.do?bug_id=7130085">JDK-7130085 : Port fdlibm hypot to Java</a>
3527     */
3528    private static double hypot(double x, double y) {
3529        // Differences to the fdlibm reference:
3530        //
3531        // 1. fdlibm orders the two parts using the magnitude of the upper 32-bits.
3532        // This incorrectly orders numbers which differ only in the lower 32-bits.
3533        // This invalidates hypot(x, y) = hypot(y, x) for small sub-normal numbers and a minority
3534        // of cases of normal numbers. This implementation forces the |x| >= |y| order
3535        // using the entire 63-bits of the unsigned doubles to ensure the function
3536        // is commutative.
3537        //
3538        // 2. fdlibm computed scaling by directly writing changes to the exponent bits
3539        // and maintained the high part (ha) during scaling for use in the high
3540        // precision sum x^2 + y^2. Since exponent scaling cannot be applied to sub-normals
3541        // the original version dropped the split number representation for sub-normals
3542        // and can produce maximum errors above 1 ULP for sub-normal numbers.
3543        // This version uses Dekker's method to split the number. This can be applied to
3544        // sub-normals and allows dropping the condition to check for sub-normal numbers
3545        // since all small numbers are handled with a single scaling factor.
3546        // The effect is increased precision for the majority of sub-normal cases where
3547        // the implementations compute a different result.
3548        //
3549        // 3. An alteration is done here to add an 'else if' instead of a second
3550        // 'if' statement. Thus you cannot scale down and up at the same time.
3551        //
3552        // 4. There is no use of the absolute double value. The magnitude comparison is
3553        // performed using the long bit representation. The computation x^2+y^2 is
3554        // insensitive to the sign bit. Thus use of Math.abs(double) is only in edge-case
3555        // branches.
3556        //
3557        // 5. The exponent different to ignore the smaller component has changed from 60 to 54.
3558        //
3559        // Original comments from fdlibm are in c style: /* */
3560        // Extra comments added for reference.
3561        //
3562        // Note that the high 32-bits are compared to constants.
3563        // The lowest 20-bits are the upper bits of the 52-bit mantissa.
3564        // The next 11-bits are the biased exponent. The sign bit has been cleared.
3565        // Scaling factors are powers of two for exact scaling.
3566        // For clarity the values have been refactored to named constants.
3567
3568        // The mask is used to remove the sign bit.
3569        final long xbits = Double.doubleToRawLongBits(x) & UNSIGN_MASK;
3570        final long ybits = Double.doubleToRawLongBits(y) & UNSIGN_MASK;
3571
3572        // Order by magnitude: |a| >= |b|
3573        double a;
3574        double b;
3575        /* High word of x & y */
3576        final int ha;
3577        final int hb;
3578        if (ybits > xbits) {
3579            a = y;
3580            b = x;
3581            ha = (int) (ybits >>> 32);
3582            hb = (int) (xbits >>> 32);
3583        } else {
3584            a = x;
3585            b = y;
3586            ha = (int) (xbits >>> 32);
3587            hb = (int) (ybits >>> 32);
3588        }
3589
3590        // Check if the smaller part is significant.
3591        // a^2 is computed in extended precision for an effective mantissa of 106-bits.
3592        // An exponent difference of 54 is where b^2 will not overlap a^2.
3593        if ((ha - hb) > EXP_54) {
3594            /* a/b > 2**54 */
3595            // or a is Inf or NaN.
3596            // No addition of a + b for sNaN.
3597            return Math.abs(a);
3598        }
3599
3600        double rescale = 1.0;
3601        if (ha > EXP_500) {
3602            /* a > 2^500 */
3603            if (ha >= EXP_1024) {
3604                /* Inf or NaN */
3605                // Check b is infinite for the IEEE754 result.
3606                // No addition of a + b for sNaN.
3607                return Math.abs(b) == Double.POSITIVE_INFINITY ?
3608                    Double.POSITIVE_INFINITY :
3609                    Math.abs(a);
3610            }
3611            /* scale a and b by 2^-600 */
3612            // Before scaling: a in [2^500, 2^1023].
3613            // After scaling: a in [2^-100, 2^423].
3614            // After scaling: b in [2^-154, 2^423].
3615            a *= TWO_POW_NEG_600;
3616            b *= TWO_POW_NEG_600;
3617            rescale = TWO_POW_600;
3618        } else if (hb < EXP_NEG_500) {
3619            // No special handling of sub-normals.
3620            // These do not matter when we do not manipulate the exponent bits
3621            // for scaling the split representation.
3622
3623            // Intentional comparison with zero.
3624            if (b == 0) {
3625                return Math.abs(a);
3626            }
3627
3628            /* scale a and b by 2^600 */
3629            // Effective min exponent of a sub-normal = -1022 - 52 = -1074.
3630            // Before scaling: b in [2^-1074, 2^-501].
3631            // After scaling: b in [2^-474, 2^99].
3632            // After scaling: a in [2^-474, 2^153].
3633            a *= TWO_POW_600;
3634            b *= TWO_POW_600;
3635            rescale = TWO_POW_NEG_600;
3636        }
3637
3638        // High precision x^2 + y^2
3639        return Math.sqrt(x2y2(a, b)) * rescale;
3640    }
3641
3642    /**
3643     * Return {@code x^2 + y^2} with high accuracy.
3644     *
3645     * <p>It is assumed that {@code 2^500 > |x| >= |y| > 2^-500}. Thus there will be no
3646     * overflow or underflow of the result. The inputs are not assumed to be unsigned.
3647     *
3648     * <p>The computation is performed using Dekker's method for extended precision
3649     * multiplication of x and y and then summation of the extended precision squares.
3650     *
3651     * @param x Value x.
3652     * @param y Value y
3653     * @return x^2 + y^2
3654     * @see <a href="https://doi.org/10.1007/BF01397083">
3655     * Dekker (1971) A floating-point technique for extending the available precision</a>
3656     */
3657    private static double x2y2(double x, double y) {
3658        // Note:
3659        // This method is different from the high-accuracy summation used in fdlibm for hypot.
3660        // The summation could be any valid computation of x^2+y^2. However since this follows
3661        // the re-scaling logic in hypot(x, y) the use of high precision has relatively
3662        // less performance overhead than if used without scaling.
3663        // The Dekker algorithm is branchless for better performance
3664        // than the fdlibm method with a maximum ULP error of approximately 0.86.
3665        //
3666        // See NUMBERS-143 for analysis.
3667
3668        // Do a Dekker summation of double length products x*x and y*y
3669        // (10 multiply and 20 additions).
3670        final double xx = x * x;
3671        final double yy = y * y;
3672        // Compute the round-off from the products.
3673        // With FMA hardware support in JDK 9+ this can be replaced with the much faster:
3674        // xxLow = Math.fma(x, x, -xx)
3675        // yyLow = Math.fma(y, y, -yy)
3676        // Dekker mul12
3677        final double xHigh = splitHigh(x);
3678        final double xLow = x - xHigh;
3679        final double xxLow = squareLow(xLow, xHigh, xx);
3680        // Dekker mul12
3681        final double yHigh = splitHigh(y);
3682        final double yLow = y - yHigh;
3683        final double yyLow = squareLow(yLow, yHigh, yy);
3684        // Dekker add2
3685        final double r = xx + yy;
3686        // Note: The order is important. Assume xx > yy and drop Dekker's conditional
3687        // check for which is the greater magnitude.
3688        // s = xx - r + yy + yyLow + xxLow
3689        // z = r + s
3690        // zz = r - z + s
3691        // Here we compute z inline and ignore computing the round-off zz.
3692        // Note: The round-off could be used with Dekker's sqrt2 method.
3693        // That adds 7 multiply, 1 division and 19 additions doubling the cost
3694        // and reducing error to < 0.5 ulp for the final sqrt.
3695        return xx - r + yy + yyLow + xxLow + r;
3696    }
3697}