Complex.java

  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.  */

  18. package org.apache.commons.numbers.complex;

  20. import java.io.Serializable;
  21. import java.util.ArrayList;
  22. import java.util.List;
  23. import java.util.function.DoubleUnaryOperator;

  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);

  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;

  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);

  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;

  208.     /** Serializable version identifier. */
  209.     private static final long serialVersionUID = 20180201L;

  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;

  231.     /** The imaginary part. */
  232.     private final double imaginary;
  233.     /** The real part. */
  234.     private final double real;

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.         }

  374.         // Confirm start: '('
  375.         if (s.charAt(0) != FORMAT_START) {
  376.             throw new NumberFormatException(
  377.                 parsingExceptionMsg("Expected start delimiter", FORMAT_START, s));
  378.         }

  380.         // Confirm end: ')'
  381.         if (s.charAt(len - 1) != FORMAT_END) {
  382.             throw new NumberFormatException(
  383.                 parsingExceptionMsg("Expected end delimiter", FORMAT_END, s));
  384.         }

  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.         }

  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.         }

  402.         // Try to parse the parts

  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.         }

  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.         }

  422.         return ofCartesian(re, im);
  423.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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;

  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.         // --------------

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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;

  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.         }

  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.     }

  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.     }

  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.     }

  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)).

  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.     }

  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.     }

  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.     }

  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.         }

  1404.         // Returns the real part:
  1405.         // log(sqrt(x^2 + y^2))
  1406.         // log(x^2 + y^2) / 2

  1408.         // Compute with positive values
  1409.         double x = Math.abs(real);
  1410.         double y = Math.abs(imaginary);

  1412.         // Find the larger magnitude.
  1413.         if (x < y) {
  1414.             final double tmp = x;
  1415.             x = y;
  1416.             y = tmp;
  1417.         }

  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.         }

  1425.         double re;

  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.

  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.         }

  1470.         // All ISO C99 edge cases for the imaginary are satisfied by the Math library.
  1471.         return constructor.create(re, arg());
  1472.     }

  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.     }

  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.     }

  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.     }

  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.         }

  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);

  1617.         // Compute
  1618.         final double t;

  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.

  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.

  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)

  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.         }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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;

  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.             }

  1875.             final double xp1 = x + 1;
  1876.             final double xm1 = x - 1;

  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;

  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.                 }

  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.         }

  1948.         return constructor.create(changeSign(re, real),
  1949.                                   changeSign(im, imaginary));
  1950.     }

  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.     }

  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;

  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.             }

  2067.             final double xp1 = x + 1;
  2068.             final double xm1 = x - 1;

  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;

  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.                 }

  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.         }

  2141.         return constructor.create(negative(real) ? Math.PI - re : re,
  2142.                                   negative(imaginary) ? im : -im);
  2143.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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);

  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.         }

  2492.         // Finite components
  2493.         // tanh(x+iy) = (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y))

  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.         }

  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.

  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.         }

  2540.         // No overflow of sinh(2x) and cosh(2x)

  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))

  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.     }

  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.     }

  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.     }

  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.     }

  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;

  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.

  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

  2768.             if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
  2769.                 // Normal computation within a safe region.

  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.

  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.                 }

  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.                 }

  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.         }

  2892.         re /= 4;
  2893.         im /= 2;
  2894.         return constructor.create(changeSign(re, real),
  2895.                                   changeSign(im, imaginary));
  2896.     }

  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:

  2911.         // If x == 1 there is no cancellation.

  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.

  2917.         // If x < 1, there can be serious cancellation:

  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.

  2922.         // Otherwise there can be serious cancellation and the relative error in the real part
  2923.         // could be enormous.

  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.

  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;

  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);

  2953.             return sumx2y2m1(xx, x2Low, yy, y2Low);
  2954.         }
  2955.         return (x - 1) * (x + 1) + yy;
  2956.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.

  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.

  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.

  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.

  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.

  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);

  3098.         final double f1 = y2Low;
  3099.         final double f2 = y2High;

  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);

  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);

  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.     }

  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.         }

  3151.         final List<Complex> result = new ArrayList<>();

  3153.         // nth root of abs -- faster / more accurate to use a solver here?
  3154.         final double nthRootOfAbs = Math.pow(abs(), 1.0 / n);

  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.         }

  3168.         return result;
  3169.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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;

  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.     }

  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.     }

  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.     }

  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.

  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;

  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.         }

  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.         }

  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.

  3623.             // Intentional comparison with zero.
  3624.             if (b == 0) {
  3625.                 return Math.abs(a);
  3626.             }

  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.         }

  3638.         // High precision x^2 + y^2
  3639.         return Math.sqrt(x2y2(a, b)) * rescale;
  3640.     }

  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.

  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. }
Created with JaCoCo 0.8.12.202403310830