View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.numbers.complex;
19  
20  import java.io.Serializable;
21  import java.util.ArrayList;
22  import java.util.List;
23  import java.util.function.DoubleUnaryOperator;
24  
25  /**
26   * Cartesian representation of a complex number. The complex number is expressed
27   * in the form \( a + ib \) where \( a \) and \( b \) are real numbers and \( i \)
28   * is the imaginary unit which satisfies the equation \( i^2 = -1 \). For the
29   * complex number \( a + ib \), \( a \) is called the <em>real part</em> and
30   * \( b \) is called the <em>imaginary part</em>.
31   *
32   * <p>This class is immutable. All arithmetic will create a new instance for the
33   * result.</p>
34   *
35   * <p>Arithmetic in this class conforms to the C99 standard for complex numbers
36   * defined in ISO/IEC 9899, Annex G. Methods have been named using the equivalent
37   * method in ISO C99. The behavior for special cases is listed as defined in C99.</p>
38   *
39   * <p>For functions \( f \) which obey the conjugate equality \( conj(f(z)) = f(conj(z)) \),
40   * the specifications for the upper half-plane imply the specifications for the lower
41   * half-plane.</p>
42   *
43   * <p>For functions that are either odd, \( f(z) = -f(-z) \), or even, \( f(z) =  f(-z) \),
44   * the specifications for the first quadrant imply the specifications for the other three
45   * quadrants.</p>
46   *
47   * <p>Special cases of <a href="http://mathworld.wolfram.com/BranchCut.html">branch cuts</a>
48   * for multivalued functions adopt the principle value convention from C99. Specials cases
49   * from C99 that raise the "invalid" or "divide-by-zero"
50   * <a href="https://en.cppreference.com/w/c/numeric/fenv/FE_exceptions">floating-point
51   * exceptions</a> return the documented value without an explicit mechanism to notify
52   * of the exception case, that is no exceptions are thrown during computations in-line with
53   * the convention of the corresponding single-valued functions in
54   * {@link Math Math}.
55   * These cases are documented in the method special cases as "invalid" or "divide-by-zero"
56   * floating-point operation.
57   * Note: Invalid floating-point exception cases will result in a complex number where the
58   * cardinality of NaN component parts has increased as a real or imaginary part could
59   * not be computed and is set to NaN.
60   *
61   * @see <a href="http://www.open-std.org/JTC1/SC22/WG14/www/standards">
62   *    ISO/IEC 9899 - Programming languages - C</a>
63   */
64  public final class Complex implements Serializable  {
65      /**
66       * A complex number representing \( i \), the square root of \( -1 \).
67       *
68       * <p>\( (0 + i 1) \).
69       */
70      public static final Complex I = new Complex(0, 1);
71      /**
72       * A complex number representing one.
73       *
74       * <p>\( (1 + i 0) \).
75       */
76      public static final Complex ONE = new Complex(1, 0);
77      /**
78       * A complex number representing zero.
79       *
80       * <p>\( (0 + i 0) \).
81       */
82      public static final Complex ZERO = new Complex(0, 0);
83  
84      /** A complex number representing {@code NaN + i NaN}. */
85      private static final Complex NAN = new Complex(Double.NaN, Double.NaN);
86      /** &pi;/2. */
87      private static final double PI_OVER_2 = 0.5 * Math.PI;
88      /** &pi;/4. */
89      private static final double PI_OVER_4 = 0.25 * Math.PI;
90      /** Natural logarithm of 2 (ln(2)). */
91      private static final double LN_2 = Math.log(2);
92      /** Base 10 logarithm of 10 divided by 2 (log10(e)/2). */
93      private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2;
94      /** Base 10 logarithm of 2 (log10(2)). */
95      private static final double LOG10_2 = Math.log10(2);
96      /** {@code 1/2}. */
97      private static final double HALF = 0.5;
98      /** {@code sqrt(2)}. */
99      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 }