commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From aherb...@apache.org
Subject [commons-numbers] 05/08: Rearrange Complex source code.
Date Wed, 01 Apr 2020 10:58:59 GMT
This is an automated email from the ASF dual-hosted git repository.

aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git

commit 50da4295f5fd4f357797d28b855366c504aaedd5
Author: aherbert <aherbert@apache.org>
AuthorDate: Wed Apr 1 10:31:52 2020 +0100

    Rearrange Complex source code.
    
    Constructors
    
    Factory Constructor methods (including parse())
    
    Properties
    real(), imag()
    
    Properties which are methods:
    abs(), arg()
    
    Other member functions:
    double-valued: norm()
    boolean-valued: isNaN, isInfinite, isFinite
    complex-valued: conj(), proj()
    
    Math:
    
    Arithmetic
    add/subtract,multiply,divide
    
    Exponential and Log functions
    
    Power functions (pow, sqrt)
    
    Trigonomic functions
    
    Hyperbolic functions
    
    Misc functions:
    nthRoot
    
    Standard object stuff:
    toString(), equals(), hashCode()
    
    Common private helper functions. Specialised helper functions may be
    located with the only method that uses them if they are not commonly
    used.
---
 .../apache/commons/numbers/complex/Complex.java    | 3912 ++++++++++----------
 1 file changed, 1956 insertions(+), 1956 deletions(-)

diff --git a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java
index 5713d65..786e904 100644
--- a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java
+++ b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java
@@ -412,71 +412,67 @@ public final class Complex implements Serializable  {
     }
 
     /**
-     * Returns {@code true} if either the real <em>or</em> imaginary component of the complex number is NaN
-     * <em>and</em> the complex number is not infinite.
-     *
-     * <p>Note that:
-     * <ul>
-     *   <li>There is more than one complex number that can return {@code true}.
-     *   <li>Different representations of NaN can be distinguished by the
-     *       {@link #equals(Object) Complex.equals(Object)} method.
-     * </ul>
+     * Creates an exception.
      *
-     * @return {@code true} if this instance contains NaN and no infinite parts.
-     * @see Double#isNaN(double)
-     * @see #isInfinite()
-     * @see #equals(Object) Complex.equals(Object)
+     * @param message Message prefix.
+     * @param error Input that caused the error.
+     * @param cause Underlying exception (if any).
+     * @return A new instance.
      */
-    public boolean isNaN() {
-        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
-            return !isInfinite();
+    private static NumberFormatException parsingException(String message,
+                                                          Object error,
+                                                          Throwable cause) {
+        // Not called with a null message or error
+        final StringBuilder sb = new StringBuilder(100)
+            .append(message)
+            .append(" '").append(error).append('\'');
+        if (cause != null) {
+            sb.append(": ").append(cause.getMessage());
         }
-        return false;
+
+        return new NumberFormatException(sb.toString());
     }
 
     /**
-     * Returns {@code true} if either real or imaginary component of the complex number is infinite.
-     *
-     * <p>Note: A complex number with at least one infinite part is regarded
-     * as an infinity (even if its other part is a NaN).
+     * Gets the real part \( a \) of this complex number \( (a + i b) \).
      *
-     * @return {@code true} if this instance contains an infinite value.
-     * @see Double#isInfinite(double)
+     * @return The real part.
      */
-    public boolean isInfinite() {
-        return Double.isInfinite(real) || Double.isInfinite(imaginary);
+    public double getReal() {
+        return real;
     }
 
     /**
-     * Returns {@code true} if both real and imaginary component of the complex number are finite.
+     * Gets the real part \( a \) of this complex number \( (a + i b) \).
      *
-     * @return {@code true} if this instance contains finite values.
-     * @see Double#isFinite(double)
+     * <p>This method is the equivalent of the C++ method {@code std::complex::real}.
+     *
+     * @return The real part.
+     * @see #getReal()
      */
-    public boolean isFinite() {
-        return Double.isFinite(real) && Double.isFinite(imaginary);
+    public double real() {
+        return getReal();
     }
 
     /**
-     * Returns the projection of this complex number onto the Riemann sphere.
-     *
-     * <p>\( z \) projects to \( z \), except that all complex infinities (even those
-     * with one infinite part and one NaN part) project to positive infinity on the real axis.
+     * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
      *
-     * If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to:
+     * @return The imaginary part.
+     */
+    public double getImaginary() {
+        return imaginary;
+    }
+
+    /**
+     * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
      *
-     * <pre>return Complex.ofCartesian(Double.POSITIVE_INFINITY, Math.copySign(0.0, z.imag());</pre>
+     * <p>This method is the equivalent of the C++ method {@code std::complex::imag}.
      *
-     * @return \( z \) projected onto the Riemann sphere.
-     * @see #isInfinite()
-     * @see <a href="http://pubs.opengroup.org/onlinepubs/9699919799/functions/cproj.html">
-     * IEEE and ISO C standards: cproj</a>
+     * @return The imaginary part.
+     * @see #getImaginary()
      */
-    public Complex proj() {
-        if (isInfinite()) {
-            return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0.0, imaginary));
-        }
-        return this;
+    public double imag() {
+        return getImaginary();
     }
 
     /**
@@ -541,6 +537,34 @@ public final class Complex implements Serializable  {
     }
 
     /**
+     * Returns the argument of this complex number.
+     *
+     * <p>The argument is the angle phi between the positive real axis and
+     * the point representing this number in the complex plane.
+     * The value returned is between \( -\pi \) (not inclusive)
+     * and \( \pi \) (inclusive), with negative values returned for numbers with
+     * negative imaginary parts.
+     *
+     * <p>If either real or imaginary part (or both) is NaN, then the result is NaN.
+     * Infinite parts are handled as {@linkplain Math#atan2} handles them,
+     * essentially treating finite parts as zero in the presence of an
+     * infinite coordinate and returning a multiple of \( \frac{\pi}{4} \) depending on
+     * the signs of the infinite parts.
+     *
+     * <p>This code follows the
+     * <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
+     * in calculating the returned value using the {@code atan2(y, x)} method for complex
+     * \( x + iy \).
+     *
+     * @return The argument of this complex number.
+     * @see Math#atan2(double, double)
+     */
+    public double arg() {
+        // Delegate
+        return Math.atan2(imaginary, real);
+    }
+
+    /**
      * Returns the squared norm value of this complex number. This is also called the absolute
      * square.
      *
@@ -570,6 +594,101 @@ public final class Complex implements Serializable  {
     }
 
     /**
+     * Returns {@code true} if either the real <em>or</em> imaginary component of the complex number is NaN
+     * <em>and</em> the complex number is not infinite.
+     *
+     * <p>Note that:
+     * <ul>
+     *   <li>There is more than one complex number that can return {@code true}.
+     *   <li>Different representations of NaN can be distinguished by the
+     *       {@link #equals(Object) Complex.equals(Object)} method.
+     * </ul>
+     *
+     * @return {@code true} if this instance contains NaN and no infinite parts.
+     * @see Double#isNaN(double)
+     * @see #isInfinite()
+     * @see #equals(Object) Complex.equals(Object)
+     */
+    public boolean isNaN() {
+        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+            return !isInfinite();
+        }
+        return false;
+    }
+
+    /**
+     * Returns {@code true} if either real or imaginary component of the complex number is infinite.
+     *
+     * <p>Note: A complex number with at least one infinite part is regarded
+     * as an infinity (even if its other part is a NaN).
+     *
+     * @return {@code true} if this instance contains an infinite value.
+     * @see Double#isInfinite(double)
+     */
+    public boolean isInfinite() {
+        return Double.isInfinite(real) || Double.isInfinite(imaginary);
+    }
+
+    /**
+     * Returns {@code true} if both real and imaginary component of the complex number are finite.
+     *
+     * @return {@code true} if this instance contains finite values.
+     * @see Double#isFinite(double)
+     */
+    public boolean isFinite() {
+        return Double.isFinite(real) && Double.isFinite(imaginary);
+    }
+
+    /**
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/ComplexConjugate.html">conjugate</a>
+     * \( \overline{z} \) of this complex number \( z \).
+     *
+     * <p>\[ z           = a + i b \\
+     *      \overline{z} = a - i b \]
+     *
+     * @return The conjugate (\( \overline{z} \)) of this complex number.
+     */
+    public Complex conj() {
+        return new Complex(real, -imaginary);
+    }
+
+    /**
+     * Returns a {@code Complex} whose value is the negation of both the real and imaginary parts
+     * of complex number \( z \).
+     *
+     * <p>\[ z  =  a + i b \\
+     *      -z  = -a - i b \]
+     *
+     * @return \( -z \).
+     */
+    public Complex negate() {
+        return new Complex(-real, -imaginary);
+    }
+
+    /**
+     * Returns the projection of this complex number onto the Riemann sphere.
+     *
+     * <p>\( z \) projects to \( z \), except that all complex infinities (even those
+     * with one infinite part and one NaN part) project to positive infinity on the real axis.
+     *
+     * If \( z \) has an infinite part, then {@code z.proj()} shall be equivalent to:
+     *
+     * <pre>return Complex.ofCartesian(Double.POSITIVE_INFINITY, Math.copySign(0.0, z.imag());</pre>
+     *
+     * @return \( z \) projected onto the Riemann sphere.
+     * @see #isInfinite()
+     * @see <a href="http://pubs.opengroup.org/onlinepubs/9699919799/functions/cproj.html">
+     * IEEE and ISO C standards: cproj</a>
+     */
+    public Complex proj() {
+        if (isInfinite()) {
+            return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0.0, imaginary));
+        }
+        return this;
+    }
+
+    /**
      * Returns a {@code Complex} whose value is {@code (this + addend)}.
      * Implements the formula:
      *
@@ -633,869 +752,1087 @@ public final class Complex implements Serializable  {
     }
 
     /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/ComplexConjugate.html">conjugate</a>
-     * \( \overline{z} \) of this complex number \( z \).
+     * Returns a {@code Complex} whose value is {@code (this - subtrahend)}.
+     * Implements the formula:
      *
-     * <p>\[ z           = a + i b \\
-     *      \overline{z} = a - i b \]
+     * <p>\[ (a + i b) - (c + i d) = (a - c) + i (b - d) \]
      *
-     * @return The conjugate (\( \overline{z} \)) of this complex number.
+     * @param  subtrahend Value to be subtracted from this complex number.
+     * @return {@code this - subtrahend}.
+     * @see <a href="http://mathworld.wolfram.com/ComplexSubtraction.html">Complex Subtraction</a>
      */
-    public Complex conj() {
-        return new Complex(real, -imaginary);
+    public Complex subtract(Complex subtrahend) {
+        return new Complex(real - subtrahend.real,
+                           imaginary - subtrahend.imaginary);
     }
 
     /**
-     * Returns a {@code Complex} whose value is {@code (this / divisor)}.
+     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
+     * with {@code subtrahend} interpreted as a real number.
      * Implements the formula:
      *
-     * <p>\[ \frac{a + i b}{c + i d} = \frac{(ac + bd) + i (bc - ad)}{c^2+d^2} \]
+     * <p>\[ (a + i b) - c = (a - c) + i b \]
      *
-     * <p>Re-calculates NaN result values to recover infinities as specified in C99 standard G.5.1.
+     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
+     * real-only and complex numbers.</p>
      *
-     * @param divisor Value by which this complex number is to be divided.
-     * @return {@code this / divisor}.
-     * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
+     * @param  subtrahend Value to be subtracted from this complex number.
+     * @return {@code this - subtrahend}.
+     * @see #subtract(Complex)
      */
-    public Complex divide(Complex divisor) {
-        return divide(real, imaginary, divisor.real, divisor.imaginary);
+    public Complex subtract(double subtrahend) {
+        return new Complex(real - subtrahend, imaginary);
     }
 
     /**
-     * Returns a {@code Complex} whose value is:
-     * <pre>
-     * <code>
-     *   a + i b     (ac + bd) + i (bc - ad)
-     *   -------  =  -----------------------
-     *   c + i d            c<sup>2</sup> + d<sup>2</sup>
-     * </code>
-     * </pre>
+     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
+     * with {@code subtrahend} interpreted as an imaginary number.
+     * Implements the formula:
      *
-     * <p>Recalculates to recover infinities as specified in C99
-     * standard G.5.1. Method is fully in accordance with
-     * C++11 standards for complex numbers.</p>
+     * <p>\[ (a + i b) - i d = a + i (b - d) \]
      *
-     * <p>Note: In the event of divide by zero this method produces the same result
-     * as dividing by a real-only zero using {@link #divide(double)}.
+     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
+     * imaginary-only and complex numbers.</p>
+     *
+     * @param  subtrahend Value to be subtracted from this complex number.
+     * @return {@code this - subtrahend}.
+     * @see #subtract(Complex)
+     */
+    public Complex subtractImaginary(double subtrahend) {
+        return new Complex(real, imaginary - subtrahend);
+    }
+
+    /**
+     * Returns a {@code Complex} whose value is {@code (minuend - this)},
+     * with {@code minuend} interpreted as a real number.
+     * Implements the formula:
+     * \[ c - (a + i b) = (c - a) - i b \]
+     *
+     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
+     * real-only and complex numbers.</p>
+     *
+     * <p>Note: This method inverts the sign of the imaginary component \( b \) if it is {@code 0.0}.
+     * The sign would not be inverted if subtracting from \( c + i 0 \) using
+     * {@link #subtract(Complex) Complex.ofCartesian(minuend, 0).subtract(this)} since
+     * {@code 0.0 - 0.0 = 0.0}.
+     *
+     * @param  minuend Value this complex number is to be subtracted from.
+     * @return {@code minuend - this}.
+     * @see #subtract(Complex)
+     * @see #ofCartesian(double, double)
+     */
+    public Complex subtractFrom(double minuend) {
+        return new Complex(minuend - real, -imaginary);
+    }
+
+    /**
+     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
+     * with {@code minuend} interpreted as an imaginary number.
+     * Implements the formula:
+     * \[ i d - (a + i b) = -a + i (d - b) \]
+     *
+     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
+     * imaginary-only and complex numbers.</p>
+     *
+     * <p>Note: This method inverts the sign of the real component \( a \) if it is {@code 0.0}.
+     * The sign would not be inverted if subtracting from \( 0 + i d \) using
+     * {@link #subtract(Complex) Complex.ofCartesian(0, minuend).subtract(this)} since
+     * {@code 0.0 - 0.0 = 0.0}.
+     *
+     * @param  minuend Value this complex number is to be subtracted from.
+     * @return {@code this - subtrahend}.
+     * @see #subtract(Complex)
+     * @see #ofCartesian(double, double)
+     */
+    public Complex subtractFromImaginary(double minuend) {
+        return new Complex(-real, minuend - imaginary);
+    }
+
+    /**
+     * Returns a {@code Complex} whose value is {@code this * factor}.
+     * Implements the formula:
+     *
+     * <p>\[ (a + i b)(c + i d) = (ac - bd) + i (ad + bc) \]
+     *
+     * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
+     *
+     * @param  factor Value to be multiplied by this complex number.
+     * @return {@code this * factor}.
+     * @see <a href="http://mathworld.wolfram.com/ComplexMultiplication.html">Complex Muliplication</a>
+     */
+    public Complex multiply(Complex factor) {
+        return multiply(real, imaginary, factor.real, factor.imaginary);
+    }
+
+    /**
+     * Returns a {@code Complex} whose value is:
+     * <pre>
+     *  (a + i b)(c + i d) = (ac - bd) + i (ad + bc)</pre>
+     *
+     * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
      *
      * @param re1 Real component of first number.
      * @param im1 Imaginary component of first number.
      * @param re2 Real component of second number.
      * @param im2 Imaginary component of second number.
-     * @return (a + i b) / (c + i d).
-     * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
-     * @see #divide(double)
+     * @return (a + b i)(c + d i).
      */
-    private static Complex divide(double re1, double im1, double re2, double im2) {
+    private static Complex multiply(double re1, double im1, double re2, double im2) {
         double a = re1;
         double b = im1;
         double c = re2;
         double d = im2;
-        int ilogbw = 0;
-        // Get the exponent to scale the divisor parts to the range [1, 2).
-        final int exponent = getScale(c, d);
-        if (exponent <= Double.MAX_EXPONENT) {
-            ilogbw = exponent;
-            c = Math.scalb(c, -ilogbw);
-            d = Math.scalb(d, -ilogbw);
-        }
-        final double denom = c * c + d * d;
+        final double ac = a * c;
+        final double bd = b * d;
+        final double ad = a * d;
+        final double bc = b * c;
+        double x = ac - bd;
+        double y = ad + bc;
 
-        // Note: Modification from the listing in ISO C99 G.5.1 (8):
-        // Avoid overflow if a or b are very big.
-        // Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
-        // when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
-        // (bc - ad) with large negative values.
-        // Use the maximum exponent as an approximation to the magnitude.
-        if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
-            ilogbw -= 2;
-            a /= 4;
-            b /= 4;
-        }
+        // --------------
+        // NaN can occur if:
+        // - any of (a,b,c,d) are NaN (for NaN or Infinite complex numbers)
+        // - a multiplication of infinity by zero (ac,bd,ad,bc).
+        // - a subtraction of infinity from infinity (e.g. ac - bd)
+        //   Note that (ac,bd,ad,bc) can be infinite due to overflow.
+        //
+        // Detect a NaN result and perform correction.
+        //
+        // Modification from the listing in ISO C99 G.5.1 (6)
+        // Do not correct infinity multiplied by zero. This is left as NaN.
+        // --------------
 
-        double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
-        double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
-        // Recover infinities and zeros that computed as NaN+iNaN
-        // the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
         if (Double.isNaN(x) && Double.isNaN(y)) {
-            if ((denom == 0.0) &&
-                    (!Double.isNaN(a) || !Double.isNaN(b))) {
-                // nonzero/zero
-                // This case produces the same result as divide by a real-only zero
-                // using Complex.divide(+/-0.0)
-                x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
-                y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
-            } else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
-                    Double.isFinite(c) && Double.isFinite(d)) {
-                // infinite/finite
+            // Recover infinities that computed as NaN+iNaN ...
+            boolean recalc = false;
+            if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
+                isNotZero(c, d)) {
+                // This complex is infinite.
+                // "Box" the infinity and change NaNs in the other factor to 0.
                 a = boxInfinity(a);
                 b = boxInfinity(b);
-                x = Double.POSITIVE_INFINITY * (a * c + b * d);
-                y = Double.POSITIVE_INFINITY * (b * c - a * d);
-            } else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
-                    Double.isFinite(a) && Double.isFinite(b)) {
-                // finite/infinite
+                c = changeNaNtoZero(c);
+                d = changeNaNtoZero(d);
+                recalc = true;
+            }
+            if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
+                isNotZero(a, b)) {
+                // The other complex is infinite.
+                // "Box" the infinity and change NaNs in the other factor to 0.
                 c = boxInfinity(c);
                 d = boxInfinity(d);
-                x = 0.0 * (a * c + b * d);
-                y = 0.0 * (b * c - a * d);
+                a = changeNaNtoZero(a);
+                b = changeNaNtoZero(b);
+                recalc = true;
+            }
+            if (!recalc && (Double.isInfinite(ac) || Double.isInfinite(bd) ||
+                            Double.isInfinite(ad) || Double.isInfinite(bc))) {
+                // The result overflowed to infinity.
+                // Recover infinities from overflow by changing NaNs to 0 ...
+                a = changeNaNtoZero(a);
+                b = changeNaNtoZero(b);
+                c = changeNaNtoZero(c);
+                d = changeNaNtoZero(d);
+                recalc = true;
+            }
+            if (recalc) {
+                x = Double.POSITIVE_INFINITY * (a * c - b * d);
+                y = Double.POSITIVE_INFINITY * (a * d + b * c);
             }
         }
         return new Complex(x, y);
     }
 
     /**
-     * Returns a {@code Complex} whose value is {@code (this / divisor)},
-     * with {@code divisor} interpreted as a real number.
-     * Implements the formula:
+     * Box values for the real or imaginary component of an infinite complex number.
+     * Any infinite value will be returned as one. Non-infinite values will be returned as zero.
+     * The sign is maintained.
      *
-     * <p>\[ \frac{a + i b}{c} = \frac{a}{c} + i \frac{b}{c} \]
+     * <pre>
+     *  inf  =  1
+     * -inf  = -1
+     *  x    =  0
+     * -x    = -0
+     * </pre>
      *
-     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
-     * real-only and complex numbers.</p>
+     * @param component the component
+     * @return The boxed value
+     */
+    private static double boxInfinity(double component) {
+        return Math.copySign(Double.isInfinite(component) ? 1.0 : 0.0, component);
+    }
+
+    /**
+     * Checks if the complex number is not zero.
      *
-     * <p>Note: This method should be preferred over using
-     * {@link #divide(Complex) divide(Complex.ofCartesian(divisor, 0))}. Division
-     * can generate signed zeros if {@code this} complex has zeros for the real
-     * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
-     * in {@link #divide(Complex)} may create zeros in the result that differ in sign
-     * from the equivalent call to divide by a real-only number.
+     * @param real the real component
+     * @param imaginary the imaginary component
+     * @return true if the complex is not zero
+     */
+    private static boolean isNotZero(double real, double imaginary) {
+        // The use of equals is deliberate.
+        // This method must distinguish NaN from zero thus ruling out:
+        // (real != 0.0 || imaginary != 0.0)
+        return !(real == 0.0 && imaginary == 0.0);
+    }
+
+    /**
+     * Change NaN to zero preserving the sign; otherwise return the value.
      *
-     * @param  divisor Value by which this complex number is to be divided.
-     * @return {@code this / divisor}.
-     * @see #divide(Complex)
+     * @param value the value
+     * @return The new value
      */
-    public Complex divide(double divisor) {
-        return new Complex(real / divisor, imaginary / divisor);
+    private static double changeNaNtoZero(double value) {
+        return Double.isNaN(value) ? Math.copySign(0.0, value) : value;
     }
 
     /**
-     * Returns a {@code Complex} whose value is {@code (this / divisor)},
-     * with {@code divisor} interpreted as an imaginary number.
+     * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
+     * interpreted as a real number.
      * Implements the formula:
      *
-     * <p>\[ \frac{a + i b}{id} = \frac{b}{d} - i \frac{a}{d} \]
+     * <p>\[ (a + i b) c =  (ac) + i (bc) \]
      *
      * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
-     * imaginary-only and complex numbers.</p>
+     * real-only and complex numbers.</p>
      *
      * <p>Note: This method should be preferred over using
-     * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))}. Division
-     * can generate signed zeros if {@code this} complex has zeros for the real
-     * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
-     * in {@link #divide(Complex)} may create zeros in the result that differ in sign
-     * from the equivalent call to divide by an imaginary-only number.
+     * {@link #multiply(Complex) multiply(Complex.ofCartesian(factor, 0))}. Multiplication
+     * can generate signed zeros if either {@code this} complex has zeros for the real
+     * and/or imaginary component, or if the factor is zero. The summation of signed zeros
+     * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
+     * from the equivalent call to multiply by a real-only number.
      *
-     * <p>Warning: This method will generate a different result from
-     * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))} if the divisor is zero.
-     * In this case the divide method using a zero-valued Complex will produce the same result
-     * as dividing by a real-only zero. The output from dividing by imaginary zero will create
-     * infinite and NaN values in the same component parts as the output from
-     * {@code this.divide(Complex.ZERO).multiplyImaginary(1)}, however the sign
-     * of some infinite values may be negated.
-     *
-     * @param  divisor Value by which this complex number is to be divided.
-     * @return {@code this / divisor}.
-     * @see #divide(Complex)
-     * @see #divide(double)
-     */
-    public Complex divideImaginary(double divisor) {
-        return new Complex(imaginary / divisor, -real / divisor);
-    }
-
-    /**
-     * Test for equality with another object. If the other object is a {@code Complex} then a
-     * comparison is made of the real and imaginary parts; otherwise {@code false} is returned.
-     *
-     * <p>If both the real and imaginary parts of two complex numbers
-     * are exactly the same the two {@code Complex} objects are considered to be equal.
-     * For this purpose, two {@code double} values are considered to be
-     * the same if and only if the method {@link Double #doubleToLongBits(double)}
-     * returns the identical {@code long} value when applied to each.
-     *
-     * <p>Note that in most cases, for two instances of class
-     * {@code Complex}, {@code c1} and {@code c2}, the
-     * value of {@code c1.equals(c2)} is {@code true} if and only if
-     *
-     * <pre>
-     *  {@code c1.getReal() == c2.getReal() && c1.getImaginary() == c2.getImaginary()}</pre>
-     *
-     * <p>also has the value {@code true}. However, there are exceptions:
-     *
-     * <ul>
-     *  <li>
-     *   Instances that contain {@code NaN} values in the same part
-     *   are considered to be equal for that part, even though {@code Double.NaN == Double.NaN}
-     *   has the value {@code false}.
-     *  </li>
-     *  <li>
-     *   Instances that share a {@code NaN} value in one part
-     *   but have different values in the other part are <em>not</em> considered equal.
-     *  </li>
-     *  <li>
-     *   Instances that contain different representations of zero in the same part
-     *   are <em>not</em> considered to be equal for that part, even though {@code -0.0 == 0.0}
-     *   has the value {@code true}.
-     *  </li>
-     * </ul>
-     *
-     * <p>The behavior is the same as if the components of the two complex numbers were passed
-     * to {@link java.util.Arrays#equals(double[], double[]) Arrays.equals(double[], double[])}:
-     *
-     * <pre>
-     *  Arrays.equals(new double[]{c1.getReal(), c1.getImaginary()},
-     *                new double[]{c2.getReal(), c2.getImaginary()}); </pre>
-     *
-     * @param other Object to test for equality with this instance.
-     * @return {@code true} if the objects are equal, {@code false} if object
-     * is {@code null}, not an instance of {@code Complex}, or not equal to
-     * this instance.
-     * @see java.lang.Double#doubleToLongBits(double)
-     * @see java.util.Arrays#equals(double[], double[])
-     */
-    @Override
-    public boolean equals(Object other) {
-        if (this == other) {
-            return true;
-        }
-        if (other instanceof Complex) {
-            final Complex c = (Complex) other;
-            return equals(real, c.real) &&
-                equals(imaginary, c.imaginary);
-        }
-        return false;
-    }
-
-    /**
-     * Gets a hash code for the complex number.
-     *
-     * <p>The behavior is the same as if the components of the complex number were passed
-     * to {@link java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])}:
-     *
-     * <pre>
-     *  {@code Arrays.hashCode(new double[] {getReal(), getImaginary()})}</pre>
-     *
-     * @return A hash code value for this object.
-     * @see java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])
+     * @param  factor Value to be multiplied by this complex number.
+     * @return {@code this * factor}.
+     * @see #multiply(Complex)
      */
-    @Override
-    public int hashCode() {
-        return 31 * (31 + Double.hashCode(real)) + Double.hashCode(imaginary);
+    public Complex multiply(double factor) {
+        return new Complex(real * factor, imaginary * factor);
     }
 
     /**
-     * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
+     * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
+     * interpreted as an imaginary number.
+     * Implements the formula:
      *
-     * @return The imaginary part.
-     */
-    public double getImaginary() {
-        return imaginary;
-    }
-
-    /**
-     * Gets the imaginary part \( b \) of this complex number \( (a + i b) \).
+     * <p>\[ (a + i b) id = (-bd) + i (ad) \]
      *
-     * <p>This method is the equivalent of the C++ method {@code std::complex::imag}.
+     * <p>This method can be used to compute the multiplication of this complex number \( z \)
+     * by \( i \) using a factor with magnitude 1.0. This should be used in preference to
+     * {@link #multiply(Complex) multiply(Complex.I)} with or without {@link #negate() negation}:</p>
      *
-     * @return The imaginary part.
-     * @see #getImaginary()
-     */
-    public double imag() {
-        return getImaginary();
-    }
-
-    /**
-     * Gets the real part \( a \) of this complex number \( (a + i b) \).
+     * \[ iz = (-b + i a) \\
+     *   -iz = (b - i a) \]
      *
-     * @return The real part.
-     */
-    public double getReal() {
-        return real;
-    }
-
-    /**
-     * Gets the real part \( a \) of this complex number \( (a + i b) \).
+     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
+     * imaginary-only and complex numbers.</p>
      *
-     * <p>This method is the equivalent of the C++ method {@code std::complex::real}.
+     * <p>Note: This method should be preferred over using
+     * {@link #multiply(Complex) multiply(Complex.ofCartesian(0, factor))}. Multiplication
+     * can generate signed zeros if either {@code this} complex has zeros for the real
+     * and/or imaginary component, or if the factor is zero. The summation of signed zeros
+     * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
+     * from the equivalent call to multiply by an imaginary-only number.
      *
-     * @return The real part.
-     * @see #getReal()
+     * @param  factor Value to be multiplied by this complex number.
+     * @return {@code this * factor}.
+     * @see #multiply(Complex)
      */
-    public double real() {
-        return getReal();
+    public Complex multiplyImaginary(double factor) {
+        return new Complex(-imaginary * factor, real * factor);
     }
 
     /**
-     * Returns a {@code Complex} whose value is {@code this * factor}.
+     * Returns a {@code Complex} whose value is {@code (this / divisor)}.
      * Implements the formula:
      *
-     * <p>\[ (a + i b)(c + i d) = (ac - bd) + i (ad + bc) \]
+     * <p>\[ \frac{a + i b}{c + i d} = \frac{(ac + bd) + i (bc - ad)}{c^2+d^2} \]
      *
-     * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
+     * <p>Re-calculates NaN result values to recover infinities as specified in C99 standard G.5.1.
      *
-     * @param  factor Value to be multiplied by this complex number.
-     * @return {@code this * factor}.
-     * @see <a href="http://mathworld.wolfram.com/ComplexMultiplication.html">Complex Muliplication</a>
+     * @param divisor Value by which this complex number is to be divided.
+     * @return {@code this / divisor}.
+     * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
      */
-    public Complex multiply(Complex factor) {
-        return multiply(real, imaginary, factor.real, factor.imaginary);
+    public Complex divide(Complex divisor) {
+        return divide(real, imaginary, divisor.real, divisor.imaginary);
     }
 
     /**
      * Returns a {@code Complex} whose value is:
      * <pre>
-     *  (a + i b)(c + i d) = (ac - bd) + i (ad + bc)</pre>
+     * <code>
+     *   a + i b     (ac + bd) + i (bc - ad)
+     *   -------  =  -----------------------
+     *   c + i d            c<sup>2</sup> + d<sup>2</sup>
+     * </code>
+     * </pre>
      *
-     * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
+     * <p>Recalculates to recover infinities as specified in C99
+     * standard G.5.1. Method is fully in accordance with
+     * C++11 standards for complex numbers.</p>
+     *
+     * <p>Note: In the event of divide by zero this method produces the same result
+     * as dividing by a real-only zero using {@link #divide(double)}.
      *
      * @param re1 Real component of first number.
      * @param im1 Imaginary component of first number.
      * @param re2 Real component of second number.
      * @param im2 Imaginary component of second number.
-     * @return (a + b i)(c + d i).
+     * @return (a + i b) / (c + i d).
+     * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
+     * @see #divide(double)
      */
-    private static Complex multiply(double re1, double im1, double re2, double im2) {
+    private static Complex divide(double re1, double im1, double re2, double im2) {
         double a = re1;
         double b = im1;
         double c = re2;
         double d = im2;
-        final double ac = a * c;
-        final double bd = b * d;
-        final double ad = a * d;
-        final double bc = b * c;
-        double x = ac - bd;
-        double y = ad + bc;
+        int ilogbw = 0;
+        // Get the exponent to scale the divisor parts to the range [1, 2).
+        final int exponent = getScale(c, d);
+        if (exponent <= Double.MAX_EXPONENT) {
+            ilogbw = exponent;
+            c = Math.scalb(c, -ilogbw);
+            d = Math.scalb(d, -ilogbw);
+        }
+        final double denom = c * c + d * d;
 
-        // --------------
-        // NaN can occur if:
-        // - any of (a,b,c,d) are NaN (for NaN or Infinite complex numbers)
-        // - a multiplication of infinity by zero (ac,bd,ad,bc).
-        // - a subtraction of infinity from infinity (e.g. ac - bd)
-        //   Note that (ac,bd,ad,bc) can be infinite due to overflow.
-        //
-        // Detect a NaN result and perform correction.
-        //
-        // Modification from the listing in ISO C99 G.5.1 (6)
-        // Do not correct infinity multiplied by zero. This is left as NaN.
-        // --------------
+        // Note: Modification from the listing in ISO C99 G.5.1 (8):
+        // Avoid overflow if a or b are very big.
+        // Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
+        // when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
+        // (bc - ad) with large negative values.
+        // Use the maximum exponent as an approximation to the magnitude.
+        if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
+            ilogbw -= 2;
+            a /= 4;
+            b /= 4;
+        }
 
+        double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
+        double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
+        // Recover infinities and zeros that computed as NaN+iNaN
+        // the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
         if (Double.isNaN(x) && Double.isNaN(y)) {
-            // Recover infinities that computed as NaN+iNaN ...
-            boolean recalc = false;
-            if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
-                isNotZero(c, d)) {
-                // This complex is infinite.
-                // "Box" the infinity and change NaNs in the other factor to 0.
+            if ((denom == 0.0) &&
+                    (!Double.isNaN(a) || !Double.isNaN(b))) {
+                // nonzero/zero
+                // This case produces the same result as divide by a real-only zero
+                // using Complex.divide(+/-0.0)
+                x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
+                y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
+            } else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
+                    Double.isFinite(c) && Double.isFinite(d)) {
+                // infinite/finite
                 a = boxInfinity(a);
                 b = boxInfinity(b);
-                c = changeNaNtoZero(c);
-                d = changeNaNtoZero(d);
-                recalc = true;
-            }
-            if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
-                isNotZero(a, b)) {
-                // The other complex is infinite.
-                // "Box" the infinity and change NaNs in the other factor to 0.
+                x = Double.POSITIVE_INFINITY * (a * c + b * d);
+                y = Double.POSITIVE_INFINITY * (b * c - a * d);
+            } else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
+                    Double.isFinite(a) && Double.isFinite(b)) {
+                // finite/infinite
                 c = boxInfinity(c);
                 d = boxInfinity(d);
-                a = changeNaNtoZero(a);
-                b = changeNaNtoZero(b);
-                recalc = true;
-            }
-            if (!recalc && (Double.isInfinite(ac) || Double.isInfinite(bd) ||
-                            Double.isInfinite(ad) || Double.isInfinite(bc))) {
-                // The result overflowed to infinity.
-                // Recover infinities from overflow by changing NaNs to 0 ...
-                a = changeNaNtoZero(a);
-                b = changeNaNtoZero(b);
-                c = changeNaNtoZero(c);
-                d = changeNaNtoZero(d);
-                recalc = true;
-            }
-            if (recalc) {
-                x = Double.POSITIVE_INFINITY * (a * c - b * d);
-                y = Double.POSITIVE_INFINITY * (a * d + b * c);
+                x = 0.0 * (a * c + b * d);
+                y = 0.0 * (b * c - a * d);
             }
         }
         return new Complex(x, y);
     }
 
     /**
-     * Box values for the real or imaginary component of an infinite complex number.
-     * Any infinite value will be returned as one. Non-infinite values will be returned as zero.
-     * The sign is maintained.
-     *
-     * <pre>
-     *  inf  =  1
-     * -inf  = -1
-     *  x    =  0
-     * -x    = -0
-     * </pre>
-     *
-     * @param component the component
-     * @return The boxed value
-     */
-    private static double boxInfinity(double component) {
-        return Math.copySign(Double.isInfinite(component) ? 1.0 : 0.0, component);
-    }
-
-    /**
-     * Checks if the complex number is not zero.
-     *
-     * @param real the real component
-     * @param imaginary the imaginary component
-     * @return true if the complex is not zero
-     */
-    private static boolean isNotZero(double real, double imaginary) {
-        // The use of equals is deliberate.
-        // This method must distinguish NaN from zero thus ruling out:
-        // (real != 0.0 || imaginary != 0.0)
-        return !(real == 0.0 && imaginary == 0.0);
-    }
-
-    /**
-     * Change NaN to zero preserving the sign; otherwise return the value.
-     *
-     * @param value the value
-     * @return The new value
-     */
-    private static double changeNaNtoZero(double value) {
-        return Double.isNaN(value) ? Math.copySign(0.0, value) : value;
-    }
-
-    /**
-     * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
-     * interpreted as a real number.
+     * Returns a {@code Complex} whose value is {@code (this / divisor)},
+     * with {@code divisor} interpreted as a real number.
      * Implements the formula:
      *
-     * <p>\[ (a + i b) c =  (ac) + i (bc) \]
+     * <p>\[ \frac{a + i b}{c} = \frac{a}{c} + i \frac{b}{c} \]
      *
      * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
      * real-only and complex numbers.</p>
      *
      * <p>Note: This method should be preferred over using
-     * {@link #multiply(Complex) multiply(Complex.ofCartesian(factor, 0))}. Multiplication
-     * can generate signed zeros if either {@code this} complex has zeros for the real
-     * and/or imaginary component, or if the factor is zero. The summation of signed zeros
-     * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
-     * from the equivalent call to multiply by a real-only number.
+     * {@link #divide(Complex) divide(Complex.ofCartesian(divisor, 0))}. Division
+     * can generate signed zeros if {@code this} complex has zeros for the real
+     * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
+     * in {@link #divide(Complex)} may create zeros in the result that differ in sign
+     * from the equivalent call to divide by a real-only number.
      *
-     * @param  factor Value to be multiplied by this complex number.
-     * @return {@code this * factor}.
-     * @see #multiply(Complex)
+     * @param  divisor Value by which this complex number is to be divided.
+     * @return {@code this / divisor}.
+     * @see #divide(Complex)
      */
-    public Complex multiply(double factor) {
-        return new Complex(real * factor, imaginary * factor);
+    public Complex divide(double divisor) {
+        return new Complex(real / divisor, imaginary / divisor);
     }
 
     /**
-     * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor}
-     * interpreted as an imaginary number.
+     * Returns a {@code Complex} whose value is {@code (this / divisor)},
+     * with {@code divisor} interpreted as an imaginary number.
      * Implements the formula:
      *
-     * <p>\[ (a + i b) id = (-bd) + i (ad) \]
-     *
-     * <p>This method can be used to compute the multiplication of this complex number \( z \)
-     * by \( i \) using a factor with magnitude 1.0. This should be used in preference to
-     * {@link #multiply(Complex) multiply(Complex.I)} with or without {@link #negate() negation}:</p>
-     *
-     * \[ iz = (-b + i a) \\
-     *   -iz = (b - i a) \]
+     * <p>\[ \frac{a + i b}{id} = \frac{b}{d} - i \frac{a}{d} \]
      *
      * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
      * imaginary-only and complex numbers.</p>
      *
      * <p>Note: This method should be preferred over using
-     * {@link #multiply(Complex) multiply(Complex.ofCartesian(0, factor))}. Multiplication
-     * can generate signed zeros if either {@code this} complex has zeros for the real
-     * and/or imaginary component, or if the factor is zero. The summation of signed zeros
-     * in {@link #multiply(Complex)} may create zeros in the result that differ in sign
-     * from the equivalent call to multiply by an imaginary-only number.
-     *
-     * @param  factor Value to be multiplied by this complex number.
-     * @return {@code this * factor}.
-     * @see #multiply(Complex)
-     */
-    public Complex multiplyImaginary(double factor) {
-        return new Complex(-imaginary * factor, real * factor);
-    }
-
-    /**
-     * Returns a {@code Complex} whose value is the negation of both the real and imaginary parts
-     * of complex number \( z \).
+     * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))}. Division
+     * can generate signed zeros if {@code this} complex has zeros for the real
+     * and/or imaginary component, or the divisor is infinite. The summation of signed zeros
+     * in {@link #divide(Complex)} may create zeros in the result that differ in sign
+     * from the equivalent call to divide by an imaginary-only number.
      *
-     * <p>\[ z  =  a + i b \\
-     *      -z  = -a - i b \]
+     * <p>Warning: This method will generate a different result from
+     * {@link #divide(Complex) divide(Complex.ofCartesian(0, divisor))} if the divisor is zero.
+     * In this case the divide method using a zero-valued Complex will produce the same result
+     * as dividing by a real-only zero. The output from dividing by imaginary zero will create
+     * infinite and NaN values in the same component parts as the output from
+     * {@code this.divide(Complex.ZERO).multiplyImaginary(1)}, however the sign
+     * of some infinite values may be negated.
      *
-     * @return \( -z \).
+     * @param  divisor Value by which this complex number is to be divided.
+     * @return {@code this / divisor}.
+     * @see #divide(Complex)
+     * @see #divide(double)
      */
-    public Complex negate() {
-        return new Complex(-real, -imaginary);
+    public Complex divideImaginary(double divisor) {
+        return new Complex(imaginary / divisor, -real / divisor);
     }
 
     /**
-     * Returns a {@code Complex} whose value is {@code (this - subtrahend)}.
-     * Implements the formula:
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
+     * exponential function</a> of this complex number.
      *
-     * <p>\[ (a + i b) - (c + i d) = (a - c) + i (b - d) \]
+     * <p>\[ \exp(z) = e^z \]
      *
-     * @param  subtrahend Value to be subtracted from this complex number.
-     * @return {@code this - subtrahend}.
-     * @see <a href="http://mathworld.wolfram.com/ComplexSubtraction.html">Complex Subtraction</a>
-     */
-    public Complex subtract(Complex subtrahend) {
-        return new Complex(real - subtrahend.real,
-                           imaginary - subtrahend.imaginary);
-    }
-
-    /**
-     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
-     * with {@code subtrahend} interpreted as a real number.
-     * Implements the formula:
+     * <p>The exponential function of \( z \) is an entire function in the complex plane.
+     * Special cases:
      *
-     * <p>\[ (a + i b) - c = (a - c) + i b \]
+     * <ul>
+     * <li>{@code z.conj().exp() == z.exp().conj()}.
+     * <li>If {@code z} is ±0 + i0, returns 1 + i0.
+     * <li>If {@code z} is x + i∞ for finite x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
+     * <li>If {@code z} is −∞ + iy for finite y, returns +0 cis(y) (see {@link #ofCis(double)}).
+     * <li>If {@code z} is +∞ + iy for finite nonzero y, returns +∞ cis(y).
+     * <li>If {@code z} is −∞ + i∞, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
+     * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
+     * <li>If {@code z} is −∞ + iNaN, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
+     * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
+     * <li>If {@code z} is NaN + i0, returns NaN + i0.
+     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
+     * </ul>
      *
-     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
-     * real-only and complex numbers.</p>
+     * <p>Implements the formula:
      *
-     * @param  subtrahend Value to be subtracted from this complex number.
-     * @return {@code this - subtrahend}.
-     * @see #subtract(Complex)
+     * <p>\[ \exp(x + iy) = e^x (\cos(y) + i \sin(y)) \]
+     *
+     * @return The exponential of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</a>
      */
-    public Complex subtract(double subtrahend) {
-        return new Complex(real - subtrahend, imaginary);
+    public Complex exp() {
+        if (Double.isInfinite(real)) {
+            // Set the scale factor applied to cis(y)
+            double zeroOrInf;
+            if (real < 0) {
+                if (!Double.isFinite(imaginary)) {
+                    // (−∞ + i∞) or (−∞ + iNaN) returns (±0 ± i0) (where the signs of the
+                    // real and imaginary parts of the result are unspecified).
+                    // Here we preserve the conjugate equality.
+                    return new Complex(0, Math.copySign(0, imaginary));
+                }
+                // (−∞ + iy) returns +0 cis(y), for finite y
+                zeroOrInf = 0;
+            } else {
+                // (+∞ + i0) returns +∞ + i0.
+                if (imaginary == 0) {
+                    return this;
+                }
+                // (+∞ + i∞) or (+∞ + iNaN) returns (±∞ + iNaN) and raises the invalid
+                // floating-point exception (where the sign of the real part of the
+                // result is unspecified).
+                if (!Double.isFinite(imaginary)) {
+                    return new Complex(real, Double.NaN);
+                }
+                // (+∞ + iy) returns (+∞ cis(y)), for finite nonzero y.
+                zeroOrInf = real;
+            }
+            return new Complex(zeroOrInf * Math.cos(imaginary),
+                               zeroOrInf * Math.sin(imaginary));
+        } else if (Double.isNaN(real)) {
+            // (NaN + i0) returns (NaN + i0)
+            // (NaN + iy) returns (NaN + iNaN) and optionally raises the invalid floating-point exception
+            // (NaN + iNaN) returns (NaN + iNaN)
+            return imaginary == 0 ? this : NAN;
+        } else if (!Double.isFinite(imaginary)) {
+            // (x + i∞) or (x + iNaN) returns (NaN + iNaN) and raises the invalid
+            // floating-point exception, for finite x.
+            return NAN;
+        }
+        // real and imaginary are finite.
+        // Compute e^a * (cos(b) + i sin(b)).
+
+        // Special case:
+        // (±0 + i0) returns (1 + i0)
+        final double exp = Math.exp(real);
+        if (imaginary == 0) {
+            return new Complex(exp, imaginary);
+        }
+        return new Complex(exp * Math.cos(imaginary),
+                           exp * Math.sin(imaginary));
     }
 
     /**
-     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
-     * with {@code subtrahend} interpreted as an imaginary number.
-     * Implements the formula:
-     *
-     * <p>\[ (a + i b) - i d = a + i (b - d) \]
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html">
+     * natural logarithm</a> of this complex number.
      *
-     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
-     * imaginary-only and complex numbers.</p>
+     * <p>The natural logarithm of \( z \) is unbounded along the real axis and
+     * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
+     * natural logarithm has a branch cut along the negative real axis \( (-infty,0] \).
+     * Special cases:
      *
-     * @param  subtrahend Value to be subtracted from this complex number.
-     * @return {@code this - subtrahend}.
-     * @see #subtract(Complex)
-     */
-    public Complex subtractImaginary(double subtrahend) {
-        return new Complex(real, imaginary - subtrahend);
-    }
-
-    /**
-     * Returns a {@code Complex} whose value is {@code (minuend - this)},
-     * with {@code minuend} interpreted as a real number.
-     * Implements the formula:
-     * \[ c - (a + i b) = (c - a) - i b \]
-     *
-     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
-     * real-only and complex numbers.</p>
+     * <ul>
+     * <li>{@code z.conj().log() == z.log().conj()}.
+     * <li>If {@code z} is −0 + i0, returns −∞ + iπ ("divide-by-zero" floating-point operation).
+     * <li>If {@code z} is +0 + i0, returns −∞ + i0 ("divide-by-zero" floating-point operation).
+     * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
+     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +∞ + iπ.
+     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
+     * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
+     * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
+     * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
+     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
+     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
+     * </ul>
      *
-     * <p>Note: This method inverts the sign of the imaginary component \( b \) if it is {@code 0.0}.
-     * The sign would not be inverted if subtracting from \( c + i 0 \) using
-     * {@link #subtract(Complex) Complex.ofCartesian(minuend, 0).subtract(this)} since
-     * {@code 0.0 - 0.0 = 0.0}.
+     * <p>Implements the formula:
      *
-     * @param  minuend Value this complex number is to be subtracted from.
-     * @return {@code minuend - this}.
-     * @see #subtract(Complex)
-     * @see #ofCartesian(double, double)
-     */
-    public Complex subtractFrom(double minuend) {
-        return new Complex(minuend - real, -imaginary);
-    }
-
-    /**
-     * Returns a {@code Complex} whose value is {@code (this - subtrahend)},
-     * with {@code minuend} interpreted as an imaginary number.
-     * Implements the formula:
-     * \[ i d - (a + i b) = -a + i (d - b) \]
+     * <p>\[ \ln(z) = \ln |z| + i \arg(z) \]
      *
-     * <p>This method is included for compatibility with ISO C99 which defines arithmetic between
-     * imaginary-only and complex numbers.</p>
+     * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
      *
-     * <p>Note: This method inverts the sign of the real component \( a \) if it is {@code 0.0}.
-     * The sign would not be inverted if subtracting from \( 0 + i d \) using
-     * {@link #subtract(Complex) Complex.ofCartesian(0, minuend).subtract(this)} since
-     * {@code 0.0 - 0.0 = 0.0}.
+     * <p>The implementation is based on the method described in:</p>
+     * <blockquote>
+     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
+     * Implementing complex elementary functions using exception handling.
+     * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
+     * </blockquote>
      *
-     * @param  minuend Value this complex number is to be subtracted from.
-     * @return {@code this - subtrahend}.
-     * @see #subtract(Complex)
-     * @see #ofCartesian(double, double)
+     * @return The natural logarithm of this complex number.
+     * @see Math#log(double)
+     * @see #abs()
+     * @see #arg()
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
      */
-    public Complex subtractFromImaginary(double minuend) {
-        return new Complex(-real, minuend - imaginary);
+    public Complex log() {
+        return log(Math::log, HALF, LN_2, Complex::ofCartesian);
     }
 
     /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/InverseCosine.html">
-     * inverse cosine</a> of this complex number.
-     *
-     * <p>\[ \cos^{-1}(z) = \frac{\pi}{2} + i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
-     *
-     * <p>The inverse cosine of \( z \) is in the range \( [0, \pi) \) along the real axis and
-     * unbounded along the imaginary axis. Special cases:
-     *
-     * <ul>
-     * <li>{@code z.conj().acos() == z.acos().conj()}.
-     * <li>If {@code z} is ±0 + i0, returns π/2 − i0.
-     * <li>If {@code z} is ±0 + iNaN, returns π/2 + iNaN.
-     * <li>If {@code z} is x + i∞ for finite x, returns π/2 − i∞.
-     * <li>If {@code z} is x + iNaN, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns π − i∞.
-     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +0 − i∞.
-     * <li>If {@code z} is −∞ + i∞, returns 3π/4 − i∞.
-     * <li>If {@code z} is +∞ + i∞, returns π/4 − i∞.
-     * <li>If {@code z} is ±∞ + iNaN, returns NaN ± i∞ where the sign of the imaginary part of the result is unspecified.
-     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is NaN + i∞, returns NaN − i∞.
-     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
-     * </ul>
-     *
-     * <p>The inverse cosine is a multivalued function and requires a branch cut in
-     * the complex plane; the cut is conventionally placed at the line segments
-     * \( (-\infty,-1) \) and \( (1,\infty) \) of the real axis.
-     *
-     * <p>This function is implemented using real \( x \) and imaginary \( y \) parts:
+     * Returns the base 10
+     * <a href="http://mathworld.wolfram.com/CommonLogarithm.html">
+     * common logarithm</a> of this complex number.
      *
-     * <p>\[ \cos^{-1}(z) = \cos^{-1}(B) - i\ \text{sgn}(y) \ln\left(A + \sqrt{A^2-1}\right) \\
-     *   A = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
-     *   B = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \]
+     * <p>The common logarithm of \( z \) is unbounded along the real axis and
+     * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
+     * common logarithm has a branch cut along the negative real axis \( (-infty,0] \).
+     * Special cases are as defined in the {@link #log() natural logarithm}:
      *
-     * <p>where \( \text{sgn}(y) \) is the sign function implemented using
-     * {@link Math#copySign(double,double) copySign(1.0, y)}.
+     * <p>Implements the formula:
      *
-     * <p>The implementation is based on the method described in:</p>
-     * <blockquote>
-     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
-     * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
-     * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
-     * </blockquote>
+     * <p>\[ \log_{10}(z) = \log_{10} |z| + i \arg(z) \]
      *
-     * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
-     * {@code c++} implementation {@code <boost/math/complex/acos.hpp>}.
+     * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
      *
-     * @return The inverse cosine of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCos/">ArcCos</a>
+     * @return The base 10 logarithm of this complex number.
+     * @see Math#log10(double)
+     * @see #abs()
+     * @see #arg()
      */
-    public Complex acos() {
-        return acos(real, imaginary, Complex::ofCartesian);
+    public Complex log10() {
+        return log(Math::log10, LOG_10E_O_2, LOG10_2, Complex::ofCartesian);
     }
 
     /**
-     * Returns the inverse cosine of the complex number.
-     *
-     * <p>This function exists to allow implementation of the identity
-     * {@code acosh(z) = +-i acos(z)}.
+     * Returns the logarithm of this complex number using the provided function.
+     * Implements the formula:
      *
-     * <p>Adapted from {@code <boost/math/complex/acos.hpp>}.
-     * The original notice is shown below and the licence is shown in full in LICENSE.txt:
      * <pre>
-     * (C) Copyright John Maddock 2005.
-     * Distributed under the Boost Software License, Version 1.0. (See accompanying
-     * file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
-     * </pre>
+     *   log(x + i y) = log(|x + i y|) + i arg(x + i y)</pre>
      *
-     * @param real Real part.
-     * @param imaginary Imaginary part.
-     * @param constructor Constructor.
-     * @return The inverse cosine of the complex number.
+     * <p>Warning: The argument {@code logOf2} must be equal to {@code log(2)} using the
+     * provided log function otherwise scaling using powers of 2 in the case of overflow
+     * will be incorrect. This is provided as an internal optimisation.
+     *
+     * @param log Log function.
+     * @param logOfeOver2 The log function applied to e, then divided by 2.
+     * @param logOf2 The log function applied to 2.
+     * @param constructor Constructor for the returned complex.
+     * @return The logarithm of this complex number.
+     * @see #abs()
+     * @see #arg()
      */
-    private static Complex acos(final double real, final double imaginary,
-                                final ComplexConstructor constructor) {
-        // Compute with positive values and determine sign at the end
-        final double x = Math.abs(real);
-        final double y = Math.abs(imaginary);
-        // The result (without sign correction)
-        double re;
-        double im;
-
-        // Handle C99 special cases
-        if (isPosInfinite(x)) {
-            if (isPosInfinite(y)) {
-                re = PI_OVER_4;
-                im = y;
-            } else if (Double.isNaN(y)) {
-                // sign of the imaginary part of the result is unspecified
-                return constructor.create(imaginary, real);
-            } else {
-                re = 0;
-                im = Double.POSITIVE_INFINITY;
-            }
-        } else if (Double.isNaN(x)) {
-            if (isPosInfinite(y)) {
-                return constructor.create(x, -imaginary);
+    private Complex log(DoubleUnaryOperator log, double logOfeOver2, double logOf2, ComplexConstructor constructor) {
+        // Handle NaN
+        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+            // Return NaN unless infinite
+            if (isInfinite()) {
+                return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
             }
             // No-use of the input constructor
             return NAN;
-        } else if (isPosInfinite(y)) {
-            re = PI_OVER_2;
-            im = y;
-        } else if (Double.isNaN(y)) {
-            return constructor.create(x == 0 ? PI_OVER_2 : y, y);
-        } else {
-            // Special case for real numbers:
-            if (y == 0 && x <= 1) {
-                return constructor.create(x == 0 ? PI_OVER_2 : Math.acos(real), -imaginary);
-            }
+        }
 
-            final double xp1 = x + 1;
-            final double xm1 = x - 1;
+        // Returns the real part:
+        // log(sqrt(x^2 + y^2))
+        // log(x^2 + y^2) / 2
 
-            if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
-                final double yy = y * y;
-                final double r = Math.sqrt(xp1 * xp1 + yy);
-                final double s = Math.sqrt(xm1 * xm1 + yy);
-                final double a = 0.5 * (r + s);
-                final double b = x / a;
+        // Compute with positive values
+        double x = Math.abs(real);
+        double y = Math.abs(imaginary);
 
-                if (b <= B_CROSSOVER) {
-                    re = Math.acos(b);
-                } else {
-                    final double apx = a + x;
-                    if (x <= 1) {
-                        re = Math.atan(Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))) / x);
-                    } else {
-                        re = Math.atan((y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))) / x);
-                    }
-                }
+        // Find the larger magnitude.
+        if (x < y) {
+            final double tmp = x;
+            x = y;
+            y = tmp;
+        }
 
-                if (a <= A_CROSSOVER) {
-                    double am1;
-                    if (x < 1) {
-                        am1 = 0.5 * (yy / (r + xp1) + yy / (s - xm1));
-                    } else {
-                        am1 = 0.5 * (yy / (r + xp1) + (s + xm1));
-                    }
-                    im = Math.log1p(am1 + Math.sqrt(am1 * (a + 1)));
-                } else {
-                    im = Math.log(a + Math.sqrt(a * a - 1));
-                }
-            } else {
-                // Hull et al: Exception handling code from figure 6
-                if (y <= (EPSILON * Math.abs(xm1))) {
-                    if (x < 1) {
-                        re = Math.acos(x);
-                        im = y / Math.sqrt(xp1 * (1 - x));
-                    } else {
-                        // This deviates from Hull et al's paper as per
-                        // https://svn.boost.org/trac/boost/ticket/7290
-                        if ((Double.MAX_VALUE / xp1) > xm1) {
-                            // xp1 * xm1 won't overflow:
-                            re = y / Math.sqrt(xm1 * xp1);
-                            im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
-                        } else {
-                            re = y / x;
-                            im = LN_2 + Math.log(x);
-                        }
-                    }
-                } else if (y <= SAFE_MIN) {
-                    // Hull et al: Assume x == 1.
-                    // True if:
-                    // E^2 > 8*sqrt(u)
-                    //
-                    // E = Machine epsilon: (1 + epsilon) = 1
-                    // u = Double.MIN_NORMAL
-                    re = Math.sqrt(y);
-                    im = Math.sqrt(y);
-                } else if (EPSILON * y - 1 >= x) {
-                    re = PI_OVER_2;
-                    im = LN_2 + Math.log(y);
-                } else if (x > 1) {
-                    re = Math.atan(y / x);
-                    final double xoy = x / y;
-                    im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
-                } else {
-                    re = PI_OVER_2;
-                    final double a = Math.sqrt(1 + y * y);
-                    im = 0.5 * Math.log1p(2 * y * (y + a));
+        if (x == 0) {
+            // Handle zero: raises the ‘‘divide-by-zero’’ floating-point exception.
+            return constructor.create(Double.NEGATIVE_INFINITY,
+                                      negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary);
+        }
+
+        double re;
+
+        // This alters the implementation of Hull et al (1994) which used a standard
+        // precision representation of |z|: sqrt(x*x + y*y).
+        // This formula should use the same definition of the magnitude returned
+        // by Complex.abs() which is a high precision computation with scaling.
+        // The checks for overflow thus only require ensuring the output of |z|
+        // will not overflow or underflow.
+
+        if (x > HALF && x < ROOT2) {
+            // x^2+y^2 close to 1. Use log1p(x^2+y^2 - 1) / 2.
+            re = Math.log1p(x2y2m1(x, y)) * logOfeOver2;
+        } else {
+            // Check for over/underflow in |z|
+            // When scaling:
+            // log(a / b) = log(a) - log(b)
+            // So initialise the result with the log of the scale factor.
+            re = 0;
+            if (x > Double.MAX_VALUE / 2) {
+                // Potential overflow.
+                if (isPosInfinite(x)) {
+                    // Handle infinity
+                    return constructor.create(x, arg());
+                }
+                // Scale down.
+                x /= 2;
+                y /= 2;
+                // log(2)
+                re = logOf2;
+            } else if (y < Double.MIN_NORMAL) {
+                // Potential underflow.
+                if (y == 0) {
+                    // Handle real only number
+                    return constructor.create(log.applyAsDouble(x), arg());
                 }
+                // Scale up sub-normal numbers to make them normal by scaling by 2^54,
+                // i.e. more than the mantissa digits.
+                x *= 0x1.0p54;
+                y *= 0x1.0p54;
+                // log(2^-54) = -54 * log(2)
+                re = -54 * logOf2;
             }
+            re += log.applyAsDouble(abs(x, y));
         }
 
-        return constructor.create(negative(real) ? Math.PI - re : re,
-                                  negative(imaginary) ? im : -im);
+        // All ISO C99 edge cases for the imaginary are satisfied by the Math library.
+        return constructor.create(re, arg());
     }
 
     /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/InverseSine.html">
-     * inverse sine</a> of this complex number.
+     * Returns the complex power of this complex number raised to the power of {@code x}.
+     * Implements the formula:
      *
-     * <p>\[ \sin^{-1}(z) = - i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
+     * <p>\[ z^x = e^{x \ln(z)} \]
      *
-     * <p>The inverse sine of \( z \) is unbounded along the imaginary axis and
-     * in the range \( [-\pi, \pi] \) along the real axis. Special cases are handled
-     * as if the operation is implemented using \( \sin^{-1}(z) = -i \sinh^{-1}(iz) \).
+     * <p>If this complex number is zero then this method returns zero if {@code x} is positive
+     * in the real component and zero in the imaginary component;
+     * otherwise it returns NaN + iNaN.
      *
-     * <p>The inverse sine is a multivalued function and requires a branch cut in
-     * the complex plane; the cut is conventionally placed at the line segments
-     * \( (\infty,-1) \) and \( (1,\infty) \) of the real axis.
+     * @param  x The exponent to which this complex number is to be raised.
+     * @return This complex number raised to the power of {@code x}.
+     * @see #log()
+     * @see #multiply(Complex)
+     * @see #exp()
+     * @see <a href="http://mathworld.wolfram.com/ComplexExponentiation.html">Complex exponentiation</a>
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
+     */
+    public Complex pow(Complex x) {
+        if (real == 0 &&
+            imaginary == 0) {
+            // This value is zero. Test the other.
+            if (x.real > 0 &&
+                x.imaginary == 0) {
+                // 0 raised to positive number is 0
+                return ZERO;
+            }
+            // 0 raised to anything else is NaN
+            return NAN;
+        }
+        return log().multiply(x).exp();
+    }
+
+    /**
+     * Returns the complex power of this complex number raised to the power of {@code x},
+     * with {@code x} interpreted as a real number.
+     * Implements the formula:
      *
-     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+     * <p>\[ z^x = e^{x \ln(z)} \]
      *
-     * <p>\[ \sin^{-1}(z) = \sin^{-1}(B) + i\ \text{sgn}(y)\ln \left(A + \sqrt{A^2-1} \right) \\
-     *   A = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
-     *   B = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \]
+     * <p>If this complex number is zero then this method returns zero if {@code x} is positive;
+     * otherwise it returns NaN + iNaN.
      *
-     * <p>where \( \text{sgn}(y) \) is the sign function implemented using
-     * {@link Math#copySign(double,double) copySign(1.0, y)}.
+     * @param  x The exponent to which this complex number is to be raised.
+     * @return This complex number raised to the power of {@code x}.
+     * @see #log()
+     * @see #multiply(double)
+     * @see #exp()
+     * @see #pow(Complex)
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
+     */
+    public Complex pow(double x) {
+        if (real == 0 &&
+            imaginary == 0) {
+            // This value is zero. Test the other.
+            if (x > 0) {
+                // 0 raised to positive number is 0
+                return ZERO;
+            }
+            // 0 raised to anything else is NaN
+            return NAN;
+        }
+        return log().multiply(x).exp();
+    }
+
+    /**
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/SquareRoot.html">
+     * square root</a> of this complex number.
      *
-     * <p>The implementation is based on the method described in:</p>
+     * <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) \]
+     *
+     * <p>The square root of \( z \) is in the range \( [0, +\infty) \) along the real axis and
+     * is unbounded along the imaginary axis. The imaginary part of the square root has a
+     * branch cut along the negative real axis \( (-infty,0) \). Special cases:
+     *
+     * <ul>
+     * <li>{@code z.conj().sqrt() == z.sqrt().conj()}.
+     * <li>If {@code z} is ±0 + i0, returns +0 + i0.
+     * <li>If {@code z} is x + i∞ for all x (including NaN), returns +∞ + i∞.
+     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +0 + i∞.
+     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
+     * <li>If {@code z} is −∞ + iNaN, returns NaN ± i∞ (where the sign of the imaginary part of the result is unspecified).
+     * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
+     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
+     * </ul>
+     *
+     * <p>Implements the following algorithm to compute \( \sqrt{x + iy} \):
+     * <ol>
+     * <li>Let \( t = \sqrt{2 (|x| + |x + iy|)} \)
+     * <li>if \( x \geq 0 \) return \( \frac{t}{2} + i \frac{y}{t} \)
+     * <li>else return \( \frac{|y|}{t} + i\ \text{sgn}(y) \frac{t}{2} \)
+     * </ol>
+     * where:
+     * <ul>
+     * <li>\( |x| =\ \){@link Math#abs(double) abs}(x)
+     * <li>\( |x + y i| =\ \){@link Complex#abs}
+     * <li>\( \text{sgn}(y) =\ \){@link Math#copySign(double,double) copySign}(1.0, y)
+     * </ul>
+     *
+     * <p>The implementation is overflow and underflow safe based on the method described in:</p>
      * <blockquote>
-     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
-     * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
-     * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
+     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
+     * Implementing complex elementary functions using exception handling.
+     * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
      * </blockquote>
      *
-     * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
-     * {@code c++} implementation {@code <boost/math/complex/asin.hpp>}.
-     *
-     * @return The inverse sine of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSin/">ArcSin</a>
+     * @return The square root of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sqrt/">Sqrt</a>
      */
-    public Complex asin() {
-        return asin(real, imaginary, Complex::ofCartesian);
+    public Complex sqrt() {
+        return sqrt(real, imaginary);
     }
 
     /**
-     * Returns the inverse sine of the complex number.
-     *
-     * <p>This function exists to allow implementation of the identity
-     * {@code asinh(z) = -i asin(iz)}.
-     *
-     * <p>Adapted from {@code <boost/math/complex/asin.hpp>}.
-     * The original notice is shown below and the licence is shown in full in LICENSE.txt:
-     * <pre>
-     * (C) Copyright John Maddock 2005.
-     * Distributed under the Boost Software License, Version 1.0. (See accompanying
-     * file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
-     * </pre>
+     * Returns the square root of the complex number {@code sqrt(x + i y)}.
      *
-     * @param real Real part.
-     * @param imaginary Imaginary part.
-     * @param constructor Constructor.
-     * @return The inverse sine of this complex number.
+     * @param real Real component.
+     * @param imaginary Imaginary component.
+     * @return The square root of the complex number.
      */
-    private static Complex asin(final double real, final double imaginary,
-                                final ComplexConstructor constructor) {
+    private static Complex sqrt(double real, double imaginary) {
+        // Handle NaN
+        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+            // Check for infinite
+            if (Double.isInfinite(imaginary)) {
+                return new Complex(Double.POSITIVE_INFINITY, imaginary);
+            }
+            if (Double.isInfinite(real)) {
+                if (real == Double.NEGATIVE_INFINITY) {
+                    return new Complex(Double.NaN, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
+                }
+                return new Complex(Double.POSITIVE_INFINITY, Double.NaN);
+            }
+            return NAN;
+        }
+
         // Compute with positive values and determine sign at the end
         final double x = Math.abs(real);
         final double y = Math.abs(imaginary);
-        // The result (without sign correction)
-        double re;
-        double im;
 
-        // Handle C99 special cases
-        if (Double.isNaN(x)) {
-            if (isPosInfinite(y)) {
-                re = x;
-                im = y;
+        // Compute
+        double t;
+
+        // This alters the implementation of Hull et al (1994) which used a standard
+        // precision representation of |z|: sqrt(x*x + y*y).
+        // This formula should use the same definition of the magnitude returned
+        // by Complex.abs() which is a high precision computation with scaling.
+        // Worry about overflow if 2 * (|z| + |x|) will overflow.
+        // Worry about underflow if |z| or |x| are sub-normal components.
+
+        if (inRegion(x, y, Double.MIN_NORMAL, SQRT_SAFE_UPPER)) {
+            // No over/underflow
+            t = Math.sqrt(2 * (abs(x, y) + x));
+        } else {
+            // Potential over/underflow. First check infinites and real/imaginary only.
+
+            // Check for infinite
+            if (isPosInfinite(y)) {
+                return new Complex(Double.POSITIVE_INFINITY, imaginary);
+            } else if (isPosInfinite(x)) {
+                if (real == Double.NEGATIVE_INFINITY) {
+                    return new Complex(0, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
+                }
+                return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0, imaginary));
+            } else if (y == 0) {
+                // Real only
+                final double sqrtAbs = Math.sqrt(x);
+                if (real < 0) {
+                    return new Complex(0, Math.copySign(sqrtAbs, imaginary));
+                }
+                return new Complex(sqrtAbs, imaginary);
+            } else if (x == 0) {
+                // Imaginary only. This sets the two components to the same magnitude.
+                // Note: In polar coordinates this does not happen:
+                // real = sqrt(abs()) * Math.cos(arg() / 2)
+                // imag = sqrt(abs()) * Math.sin(arg() / 2)
+                // arg() / 2 = pi/4 and cos and sin should both return sqrt(2)/2 but
+                // are different by 1 ULP.
+                final double sqrtAbs = Math.sqrt(y) * ONE_OVER_ROOT2;
+                return new Complex(sqrtAbs, Math.copySign(sqrtAbs, imaginary));
+            } else {
+                // Over/underflow.
+                // Full scaling is not required as this is done in the hypotenuse function.
+                // Keep the number as big as possible for maximum precision in the second sqrt.
+                // Note if we scale by an even power of 2, we can re-scale by sqrt of the number.
+                // a = sqrt(b)
+                // a = sqrt(b/4) * sqrt(4)
+
+                double rescale;
+                double sx;
+                double sy;
+                if (Math.max(x, y) > SQRT_SAFE_UPPER) {
+                    // Overflow. Scale down by 16 and rescale by sqrt(16).
+                    sx = x / 16;
+                    sy = y / 16;
+                    rescale = 4;
+                } else {
+                    // Sub-normal numbers. Make them normal by scaling by 2^54,
+                    // i.e. more than the mantissa digits, and rescale by sqrt(2^54) = 2^27.
+                    sx = x * 0x1.0p54;
+                    sy = y * 0x1.0p54;
+                    rescale = 0x1.0p-27;
+                }
+                t = rescale * Math.sqrt(2 * (abs(sx, sy) + sx));
+            }
+        }
+
+        if (real >= 0) {
+            return new Complex(t / 2, imaginary / t);
+        }
+        return new Complex(y / t, Math.copySign(t / 2, imaginary));
+    }
+
+    /**
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/Sine.html">
+     * sine</a> of this complex number.
+     *
+     * <p>\[ \sin(z) = \frac{1}{2} i \left( e^{-iz} - e^{iz} \right) \]
+     *
+     * <p>This is an odd function: \( \sin(z) = -\sin(-z) \).
+     * The sine is an entire function and requires no branch cuts.
+     *
+     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+     *
+     * <p>\[ \sin(x + iy) = \sin(x)\cosh(y) + i \cos(x)\sinh(y) \]
+     *
+     * <p>As per the C99 standard this function is computed using the trigonomic identity:
+     *
+     * <p>\[ \sin(z) = -i \sinh(iz) \]
+     *
+     * @return The sine of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sin/">Sin</a>
+     */
+    public Complex sin() {
+        // Define in terms of sinh
+        // sin(z) = -i sinh(iz)
+        // Multiply this number by I, compute sinh, then multiply by back
+        return sinh(-imaginary, real, Complex::multiplyNegativeI);
+    }
+
+    /**
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/Cosine.html">
+     * cosine</a> of this complex number.
+     *
+     * <p>\[ \cos(z) = \frac{1}{2} \left( e^{iz} + e^{-iz} \right) \]
+     *
+     * <p>This is an even function: \( \cos(z) = \cos(-z) \).
+     * The cosine is an entire function and requires no branch cuts.
+     *
+     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+     *
+     * <p>\[ \cos(x + iy) = \cos(x)\cosh(y) - i \sin(x)\sinh(y) \]
+     *
+     * <p>As per the C99 standard this function is computed using the trigonomic identity:
+     *
+     * <p>\[ cos(z) = cosh(iz) \]
+     *
+     * @return The cosine of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cos/">Cos</a>
+     */
+    public Complex cos() {
+        // Define in terms of cosh
+        // cos(z) = cosh(iz)
+        // Multiply this number by I and compute cosh.
+        return cosh(-imaginary, real, Complex::ofCartesian);
+    }
+
+    /**
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/Tangent.html">
+     * tangent</a> of this complex number.
+     *
+     * <p>\[ \tan(z) = \frac{i(e^{-iz} - e^{iz})}{e^{-iz} + e^{iz}} \]
+     *
+     * <p>This is an odd function: \( \tan(z) = -\tan(-z) \).
+     * The tangent is an entire function and requires no branch cuts.
+     *
+     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:</p>
+     * \[ \tan(x + iy) = \frac{\sin(2x)}{\cos(2x)+\cosh(2y)} + i \frac{\sinh(2y)}{\cos(2x)+\cosh(2y)} \]
+     *
+     * <p>As per the C99 standard this function is computed using the trigonomic identity:</p>
+     * \[ \tan(z) = -i \tanh(iz) \]
+     *
+     * @return The tangent of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tan/">Tangent</a>
+     */
+    public Complex tan() {
+        // Define in terms of tanh
+        // tan(z) = -i tanh(iz)
+        // Multiply this number by I, compute tanh, then multiply by back
+        return tanh(-imaginary, real, Complex::multiplyNegativeI);
+    }
+
+    /**
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/InverseSine.html">
+     * inverse sine</a> of this complex number.
+     *
+     * <p>\[ \sin^{-1}(z) = - i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
+     *
+     * <p>The inverse sine of \( z \) is unbounded along the imaginary axis and
+     * in the range \( [-\pi, \pi] \) along the real axis. Special cases are handled
+     * as if the operation is implemented using \( \sin^{-1}(z) = -i \sinh^{-1}(iz) \).
+     *
+     * <p>The inverse sine is a multivalued function and requires a branch cut in
+     * the complex plane; the cut is conventionally placed at the line segments
+     * \( (\infty,-1) \) and \( (1,\infty) \) of the real axis.
+     *
+     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+     *
+     * <p>\[ \sin^{-1}(z) = \sin^{-1}(B) + i\ \text{sgn}(y)\ln \left(A + \sqrt{A^2-1} \right) \\
+     *   A = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
+     *   B = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \]
+     *
+     * <p>where \( \text{sgn}(y) \) is the sign function implemented using
+     * {@link Math#copySign(double,double) copySign(1.0, y)}.
+     *
+     * <p>The implementation is based on the method described in:</p>
+     * <blockquote>
+     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
+     * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
+     * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
+     * </blockquote>
+     *
+     * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+     * {@code c++} implementation {@code <boost/math/complex/asin.hpp>}.
+     *
+     * @return The inverse sine of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSin/">ArcSin</a>
+     */
+    public Complex asin() {
+        return asin(real, imaginary, Complex::ofCartesian);
+    }
+
+    /**
+     * Returns the inverse sine of the complex number.
+     *
+     * <p>This function exists to allow implementation of the identity
+     * {@code asinh(z) = -i asin(iz)}.
+     *
+     * <p>Adapted from {@code <boost/math/complex/asin.hpp>}.
+     * The original notice is shown below and the licence is shown in full in LICENSE.txt:
+     * <pre>
+     * (C) Copyright John Maddock 2005.
+     * Distributed under the Boost Software License, Version 1.0. (See accompanying
+     * file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+     * </pre>
+     *
+     * @param real Real part.
+     * @param imaginary Imaginary part.
+     * @param constructor Constructor.
+     * @return The inverse sine of this complex number.
+     */
+    private static Complex asin(final double real, final double imaginary,
+                                final ComplexConstructor constructor) {
+        // Compute with positive values and determine sign at the end
+        final double x = Math.abs(real);
+        final double y = Math.abs(imaginary);
+        // The result (without sign correction)
+        double re;
+        double im;
+
+        // Handle C99 special cases
+        if (Double.isNaN(x)) {
+            if (isPosInfinite(y)) {
+                re = x;
+                im = y;
             } else {
                 // No-use of the input constructor
                 return NAN;
@@ -1602,130 +1939,67 @@ public final class Complex implements Serializable  {
 
     /**
      * Returns the
-     * <a href="http://mathworld.wolfram.com/InverseTangent.html">
-     * inverse tangent</a> of this complex number.
-     *
-     * <p>\[ \tan^{-1}(z) = \frac{i}{2} \ln \left( \frac{i + z}{i - z} \right) \]
-     *
-     * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the imaginary axis and
-     * in the range \( [-\pi/2, \pi/2] \) along the real axis.
+     * <a href="http://mathworld.wolfram.com/InverseCosine.html">
+     * inverse cosine</a> of this complex number.
      *
-     * <p>The inverse tangent is a multivalued function and requires a branch cut in
-     * the complex plane; the cut is conventionally placed at the line segments
-     * \( (i \infty,-i] \) and \( [i,i \infty) \) of the imaginary axis.
+     * <p>\[ \cos^{-1}(z) = \frac{\pi}{2} + i \left(\ln{iz + \sqrt{1 - z^2}}\right) \]
      *
-     * <p>As per the C99 standard this function is computed using the trigonomic identity:
-     * \[ \tan^{-1}(z) = -i \tanh^{-1}(iz) \]
-     *
-     * @return The inverse tangent of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTan/">ArcTan</a>
-     */
-    public Complex atan() {
-        // Define in terms of atanh
-        // atan(z) = -i atanh(iz)
-        // Multiply this number by I, compute atanh, then multiply by back
-        return atanh(-imaginary, real, Complex::multiplyNegativeI);
-    }
-
-    /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/InverseHyperbolicSine.html">
-     * inverse hyperbolic sine</a> of this complex number.
-     *
-     * <p>\[ \sinh^{-1}(z) = \ln \left(z + \sqrt{1 + z^2} \right) \]
-     *
-     * <p>The inverse hyperbolic sine of \( z \) is unbounded along the real axis and
-     * in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
+     * <p>The inverse cosine of \( z \) is in the range \( [0, \pi) \) along the real axis and
+     * unbounded along the imaginary axis. Special cases:
      *
      * <ul>
-     * <li>{@code z.conj().asinh() == z.asinh().conj()}.
-     * <li>This is an odd function: \( \sinh^{-1}(z) = -\sinh^{-1}(-z) \).
-     * <li>If {@code z} is +0 + i0, returns 0 + i0.
-     * <li>If {@code z} is x + i∞ for positive-signed finite x, returns +∞ + iπ/2.
-     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
-     * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
-     * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
-     * <li>If {@code z} is NaN + i0, returns NaN + i0.
-     * <li>If {@code z} is NaN + iy for finite nonzero y, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is NaN + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
+     * <li>{@code z.conj().acos() == z.acos().conj()}.
+     * <li>If {@code z} is ±0 + i0, returns π/2 − i0.
+     * <li>If {@code z} is ±0 + iNaN, returns π/2 + iNaN.
+     * <li>If {@code z} is x + i∞ for finite x, returns π/2 − i∞.
+     * <li>If {@code z} is x + iNaN, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns π − i∞.
+     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +0 − i∞.
+     * <li>If {@code z} is −∞ + i∞, returns 3π/4 − i∞.
+     * <li>If {@code z} is +∞ + i∞, returns π/4 − i∞.
+     * <li>If {@code z} is ±∞ + iNaN, returns NaN ± i∞ where the sign of the imaginary part of the result is unspecified.
+     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + i∞, returns NaN − i∞.
      * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
      * </ul>
      *
-     * <p>The inverse hyperbolic sine is a multivalued function and requires a branch cut in
+     * <p>The inverse cosine is a multivalued function and requires a branch cut in
      * the complex plane; the cut is conventionally placed at the line segments
-     * \( (-i \infty,-i) \) and \( (i,i \infty) \) of the imaginary axis.
-     *
-     * <p>This function is computed using the trigonomic identity:
-     *
-     * <p>\[ \sinh^{-1}(z) = -i \sin^{-1}(iz) \]
-     *
-     * @return The inverse hyperbolic sine of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSinh/">ArcSinh</a>
-     */
-    public Complex asinh() {
-        // Define in terms of asin
-        // asinh(z) = -i asin(iz)
-        // Note: This is the opposite to the identity defined in the C99 standard:
-        // asin(z) = -i asinh(iz)
-        // Multiply this number by I, compute asin, then multiply by back
-        return asin(-imaginary, real, Complex::multiplyNegativeI);
-    }
-
-    /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/InverseHyperbolicTangent.html">
-     * inverse hyperbolic tangent</a> of this complex number.
-     *
-     * <p>\[ \tanh^{-1}(z) = \frac{1}{2} \ln \left( \frac{1 + z}{1 - z} \right) \]
-     *
-     * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the real axis and
-     * in the range \( [-\pi/2, \pi/2] \) along the imaginary axis. Special cases:
-     *
-     * <ul>
-     * <li>{@code z.conj().atanh() == z.atanh().conj()}.
-     * <li>This is an odd function: \( \tanh^{-1}(z) = -\tanh^{-1}(-z) \).
-     * <li>If {@code z} is +0 + i0, returns +0 + i0.
-     * <li>If {@code z} is +0 + iNaN, returns +0 + iNaN.
-     * <li>If {@code z} is +1 + i0, returns +∞ + i0 ("divide-by-zero" floating-point operation).
-     * <li>If {@code z} is x + i∞ for finite positive-signed x, returns +0 + iπ/2.
-     * <li>If {@code z} is x+iNaN for nonzero finite x, returns NaN+iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +0 + iπ/2.
-     * <li>If {@code z} is +∞ + i∞, returns +0 + iπ/2.
-     * <li>If {@code z} is +∞ + iNaN, returns +0 + iNaN.
-     * <li>If {@code z} is NaN+iy for finite y, returns NaN+iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is NaN + i∞, returns ±0 + iπ/2 (where the sign of the real part of the result is unspecified).
-     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
-     * </ul>
+     * \( (-\infty,-1) \) and \( (1,\infty) \) of the real axis.
      *
-     * <p>The inverse hyperbolic tangent is a multivalued function and requires a branch cut in
-     * the complex plane; the cut is conventionally placed at the line segments
-     * \( (\infty,-1] \) and \( [1,\infty) \) of the real axis.
+     * <p>This function is implemented using real \( x \) and imaginary \( y \) parts:
      *
-     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+     * <p>\[ \cos^{-1}(z) = \cos^{-1}(B) - i\ \text{sgn}(y) \ln\left(A + \sqrt{A^2-1}\right) \\
+     *   A = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} + \sqrt{(x-1)^2+y^2} \right] \\
+     *   B = \frac{1}{2} \left[ \sqrt{(x+1)^2+y^2} - \sqrt{(x-1)^2+y^2} \right] \]
      *
-     * <p>\[ \tanh^{-1}(z) = \frac{1}{4} \ln \left(1 + \frac{4x}{(1-x)^2+y^2} \right) + \\
-     *                     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) \]
+     * <p>where \( \text{sgn}(y) \) is the sign function implemented using
+     * {@link Math#copySign(double,double) copySign(1.0, y)}.
      *
-     * <p>The imaginary part is computed using {@link Math#atan2(double, double)} to ensure the
-     * correct quadrant is returned from \( \tan^{-1} \left(\frac{2y}{1-x^2-y^2} \right) \).
+     * <p>The implementation is based on the method described in:</p>
+     * <blockquote>
+     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1997)
+     * Implementing the complex Arcsine and Arccosine Functions using Exception Handling.
+     * ACM Transactions on Mathematical Software, Vol 23, No 3, pp 299-335.
+     * </blockquote>
      *
      * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
-     * {@code c++} implementation {@code <boost/math/complex/atanh.hpp>}.
+     * {@code c++} implementation {@code <boost/math/complex/acos.hpp>}.
      *
-     * @return The inverse hyperbolic tangent of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTanh/">ArcTanh</a>
+     * @return The inverse cosine of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCos/">ArcCos</a>
      */
-    public Complex atanh() {
-        return atanh(real, imaginary, Complex::ofCartesian);
+    public Complex acos() {
+        return acos(real, imaginary, Complex::ofCartesian);
     }
 
     /**
-     * Returns the inverse hyperbolic tangent of this complex number.
+     * Returns the inverse cosine of the complex number.
      *
      * <p>This function exists to allow implementation of the identity
-     * {@code atan(z) = -i atanh(iz)}.
+     * {@code acosh(z) = +-i acos(z)}.
      *
+     * <p>Adapted from {@code <boost/math/complex/acos.hpp>}.
      * The original notice is shown below and the licence is shown in full in LICENSE.txt:
      * <pre>
      * (C) Copyright John Maddock 2005.
@@ -1736,264 +2010,230 @@ public final class Complex implements Serializable  {
      * @param real Real part.
      * @param imaginary Imaginary part.
      * @param constructor Constructor.
-     * @return The inverse hyperbolic tangent of the complex number.
+     * @return The inverse cosine of the complex number.
      */
-    private static Complex atanh(final double real, final double imaginary,
-                                 final ComplexConstructor constructor) {
+    private static Complex acos(final double real, final double imaginary,
+                                final ComplexConstructor constructor) {
         // Compute with positive values and determine sign at the end
-        double x = Math.abs(real);
-        double y = Math.abs(imaginary);
+        final double x = Math.abs(real);
+        final double y = Math.abs(imaginary);
         // The result (without sign correction)
         double re;
         double im;
 
         // Handle C99 special cases
-        if (Double.isNaN(x)) {
+        if (isPosInfinite(x)) {
             if (isPosInfinite(y)) {
-                // The sign of the real part of the result is unspecified
-                return constructor.create(0, Math.copySign(PI_OVER_2, imaginary));
-            }
-            // No-use of the input constructor.
-            // Optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
-            return NAN;
-        } else if (Double.isNaN(y)) {
-            if (isPosInfinite(x)) {
-                return constructor.create(Math.copySign(0, real), Double.NaN);
+                re = PI_OVER_4;
+                im = y;
+            } else if (Double.isNaN(y)) {
+                // sign of the imaginary part of the result is unspecified
+                return constructor.create(imaginary, real);
+            } else {
+                re = 0;
+                im = Double.POSITIVE_INFINITY;
             }
-            if (x == 0) {
-                return constructor.create(real, Double.NaN);
+        } else if (Double.isNaN(x)) {
+            if (isPosInfinite(y)) {
+                return constructor.create(x, -imaginary);
             }
             // No-use of the input constructor
             return NAN;
+        } else if (isPosInfinite(y)) {
+            re = PI_OVER_2;
+            im = y;
+        } else if (Double.isNaN(y)) {
+            return constructor.create(x == 0 ? PI_OVER_2 : y, y);
         } else {
-            // x && y are finite or infinite.
-
-            // Check the safe region.
-            // The lower and upper bounds have been copied from boost::math::atanh.
-            // They are different from the safe region for asin and acos.
-            // x >= SAFE_UPPER: (1-x) == -x
-            // x <= SAFE_LOWER: 1 - x^2 = 1
+            // Special case for real numbers:
+            if (y == 0 && x <= 1) {
+                return constructor.create(x == 0 ? PI_OVER_2 : Math.acos(real), -imaginary);
+            }
 
-            if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
-                // Normal computation within a safe region.
+            final double xp1 = x + 1;
+            final double xm1 = x - 1;
 
-                // minus x plus 1: (-x+1)
-                final double mxp1 = 1 - x;
+            if (inRegion(x, y, SAFE_MIN, SAFE_MAX)) {
                 final double yy = y * y;
-                // The definition of real component is:
-                // real = log( ((x+1)^2+y^2) / ((1-x)^2+y^2) ) / 4
-                // This simplifies by adding 1 and subtracting 1 as a fraction:
-                //      = log(1 + ((x+1)^2+y^2) / ((1-x)^2+y^2) - ((1-x)^2+y^2)/((1-x)^2+y^2) ) / 4
-                //
-                // real(atanh(z)) == log(1 + 4*x / ((1-x)^2+y^2)) / 4
-                // imag(atanh(z)) == tan^-1 (2y, (1-x)(1+x) - y^2) / 2
-                // imag(atanh(z)) == tan^-1 (2y, (1 - x^2 - y^2) / 2
-                // The division is done at the end of the function.
-                re = Math.log1p(4 * x / (mxp1 * mxp1 + yy));
-                // Modified from boost which does not switch the magnitude of x and y.
-                // The denominator for atan2 is 1 - x^2 - y^2.
-                // This can be made more precise if |x| > |y|.
-                final double numerator = 2 * y;
-                double denominator;
-                if (x < y) {
-                    final double tmp = x;
-                    x = y;
-                    y = tmp;
-                }
-                // 1 - x is precise if |x| >= 1
-                if (x >= 1) {
-                    denominator = (1 - x) * (1 + x) - y * y;
-                } else {
-                    // |x| < 1: Use high precision if possible:
-                    // 1 - x^2 - y^2 = -(x^2 + y^2 - 1)
-                    denominator = -x2y2m1(x, y);
-                }
-                im = Math.atan2(numerator, denominator);
-            } else {
-                // This section handles exception cases that would normally cause
-                // underflow or overflow in the main formulas.
+                final double r = Math.sqrt(xp1 * xp1 + yy);
+                final double s = Math.sqrt(xm1 * xm1 + yy);
+                final double a = 0.5 * (r + s);
+                final double b = x / a;
 
-                // C99. G.7: Special case for imaginary only numbers
-                if (x == 0) {
-                    if (imaginary == 0) {
-                        return constructor.create(real, imaginary);
+                if (b <= B_CROSSOVER) {
+                    re = Math.acos(b);
+                } else {
+                    final double apx = a + x;
+                    if (x <= 1) {
+                        re = Math.atan(Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))) / x);
+                    } else {
+                        re = Math.atan((y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))) / x);
                     }
-                    // atanh(iy) = i atan(y)
-                    return constructor.create(real, Math.atan(imaginary));
                 }
 
-                // Real part:
-                // real = Math.log1p(4x / ((1-x)^2 + y^2))
-                // real = Math.log1p(4x / (1 - 2x + x^2 + y^2))
-                // real = Math.log1p(4x / (1 + x(x-2) + y^2))
-                // without either overflow or underflow in the squared terms.
-                if (x >= SAFE_UPPER) {
-                    // (1-x) = -x to machine precision:
-                    // log1p(4x / (x^2 + y^2))
-                    if (isPosInfinite(x) || isPosInfinite(y)) {
-                        re = 0;
-                    } else if (y >= SAFE_UPPER) {
-                        // Big x and y: divide by x*y
-                        re = Math.log1p((4 / y) / (x / y + y / x));
-                    } else if (y > 1) {
-                        // Big x: divide through by x:
-                        re = Math.log1p(4 / (x + y * y / x));
-                    } else {
-                        // Big x small y, as above but neglect y^2/x:
-                        re = Math.log1p(4 / x);
-                    }
-                } else if (y >= SAFE_UPPER) {
-                    if (x > 1) {
-                        // Big y, medium x, divide through by y:
-                        final double mxp1 = 1 - x;
-                        re = Math.log1p((4 * x / y) / (mxp1 * mxp1 / y + y));
+                if (a <= A_CROSSOVER) {
+                    double am1;
+                    if (x < 1) {
+                        am1 = 0.5 * (yy / (r + xp1) + yy / (s - xm1));
                     } else {
-                        // Big y, small x, as above but neglect (1-x)^2/y:
-                        // Note: log1p(v) == v - v^2/2 + v^3/3 ... Taylor series when v is small.
-                        // Here v is so small only the first term matters.
-                        re = 4 * x / y / y;
+                        am1 = 0.5 * (yy / (r + xp1) + (s + xm1));
                     }
-                } else if (x == 1) {
-                    // x = 1, small y:
-                    // Special case when x == 1 as (1-x) is invalid.
-                    // Simplify the following formula:
-                    // real = log( sqrt((x+1)^2+y^2) ) / 2 - log( sqrt((1-x)^2+y^2) ) / 2
-                    //      = log( sqrt(4+y^2) ) / 2 - log(y) / 2
-                    // if: 4+y^2 -> 4
-                    //      = log( 2 ) / 2 - log(y) / 2
-                    //      = (log(2) - log(y)) / 2
-                    // Multiply by 2 as it will be divided by 4 at the end.
-                    // C99: if y=0 raises the ‘‘divide-by-zero’’ floating-point exception.
-                    re = 2 * (LN_2 - Math.log(y));
+                    im = Math.log1p(am1 + Math.sqrt(am1 * (a + 1)));
                 } else {
-                    // Modified from boost which checks y > SAFE_LOWER.
-                    // if y*y -> 0 it will be ignored so always include it.
-                    final double mxp1 = 1 - x;
-                    re = Math.log1p((4 * x) / (mxp1 * mxp1 + y * y));
+                    im = Math.log(a + Math.sqrt(a * a - 1));
                 }
-
-                // Imaginary part:
-                // imag = atan2(2y, (1-x)(1+x) - y^2)
-                // if x or y are large, then the formula:
-                //   atan2(2y, (1-x)(1+x) - y^2)
-                // evaluates to +(PI - theta) where theta is negligible compared to PI.
-                if ((x >= SAFE_UPPER) || (y >= SAFE_UPPER)) {
-                    im = Math.PI;
-                } else if (x <= SAFE_LOWER) {
-                    // (1-x)^2 -> 1
-                    if (y <= SAFE_LOWER) {
-                        // 1 - y^2 -> 1
-                        im = Math.atan2(2 * y, 1);
+            } else {
+                // Hull et al: Exception handling code from figure 6
+                if (y <= (EPSILON * Math.abs(xm1))) {
+                    if (x < 1) {
+                        re = Math.acos(x);
+                        im = y / Math.sqrt(xp1 * (1 - x));
                     } else {
-                        im = Math.atan2(2 * y, 1 - y * y);
+                        // This deviates from Hull et al's paper as per
+                        // https://svn.boost.org/trac/boost/ticket/7290
+                        if ((Double.MAX_VALUE / xp1) > xm1) {
+                            // xp1 * xm1 won't overflow:
+                            re = y / Math.sqrt(xm1 * xp1);
+                            im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
+                        } else {
+                            re = y / x;
+                            im = LN_2 + Math.log(x);
+                        }
                     }
+                } else if (y <= SAFE_MIN) {
+                    // Hull et al: Assume x == 1.
+                    // True if:
+                    // E^2 > 8*sqrt(u)
+                    //
+                    // E = Machine epsilon: (1 + epsilon) = 1
+                    // u = Double.MIN_NORMAL
+                    re = Math.sqrt(y);
+                    im = Math.sqrt(y);
+                } else if (EPSILON * y - 1 >= x) {
+                    re = PI_OVER_2;
+                    im = LN_2 + Math.log(y);
+                } else if (x > 1) {
+                    re = Math.atan(y / x);
+                    final double xoy = x / y;
+                    im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
                 } else {
-                    // Medium x, small y.
-                    // Modified from boost which checks (y == 0) && (x == 1) and sets re = 0.
-                    // This is same as the result from calling atan2(0, 0) so exclude this case.
-                    // 1 - y^2 = 1 so ignore subtracting y^2
-                    im = Math.atan2(2 * y, (1 - x) * (1 + x));
+                    re = PI_OVER_2;
+                    final double a = Math.sqrt(1 + y * y);
+                    im = 0.5 * Math.log1p(2 * y * (y + a));
                 }
             }
         }
 
-        re /= 4;
-        im /= 2;
-        return constructor.create(changeSign(re, real),
-                                  changeSign(im, imaginary));
+        return constructor.create(negative(real) ? Math.PI - re : re,
+                                  negative(imaginary) ? im : -im);
     }
 
     /**
      * Returns the
-     * <a href="http://mathworld.wolfram.com/InverseHyperbolicCosine.html">
-     * inverse hyperbolic cosine</a> of this complex number.
-     *
-     * <p>\[ \cosh^{-1}(z) = \ln \left(z + \sqrt{z + 1} \sqrt{z - 1} \right) \]
-     *
-     * <p>The inverse hyperbolic cosine of \( z \) is in the range \( [0, \infty) \) along the
-     * real axis and in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
-     *
-     * <ul>
-     * <li>{@code z.conj().acosh() == z.acosh().conj()}.
-     * <li>If {@code z} is ±0 + i0, returns +0 + iπ/2.
-     * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
-     * <li>If {@code z} is 0 + iNaN, returns NaN + iπ/2 <sup>[1]</sup>.
-     * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns +∞ + iπ.
-     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
-     * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
-     * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
-     * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
-     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
-     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
-     * </ul>
-     *
-     * <p>Special cases include the technical corrigendum
-     * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
-     * DR 471: Complex math functions cacosh and ctanh</a>.
+     * <a href="http://mathworld.wolfram.com/InverseTangent.html">
+     * inverse tangent</a> of this complex number.
      *
-     * <p>The inverse hyperbolic cosine is a multivalued function and requires a branch cut in
-     * the complex plane; the cut is conventionally placed at the line segment
-     * \( (-\infty,-1) \) of the real axis.
+     * <p>\[ \tan^{-1}(z) = \frac{i}{2} \ln \left( \frac{i + z}{i - z} \right) \]
      *
-     * <p>This function is computed using the trigonomic identity:
+     * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the imaginary axis and
+     * in the range \( [-\pi/2, \pi/2] \) along the real axis.
      *
-     * <p>\[ \cosh^{-1}(z) = \pm i \cos^{-1}(z) \]
+     * <p>The inverse tangent is a multivalued function and requires a branch cut in
+     * the complex plane; the cut is conventionally placed at the line segments
+     * \( (i \infty,-i] \) and \( [i,i \infty) \) of the imaginary axis.
      *
-     * <p>The sign of the multiplier is chosen to give {@code z.acosh().real() >= 0}
-     * and compatibility with the C99 standard.
+     * <p>As per the C99 standard this function is computed using the trigonomic identity:
+     * \[ \tan^{-1}(z) = -i \tanh^{-1}(iz) \]
      *
-     * @return The inverse hyperbolic cosine of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCosh/">ArcCosh</a>
+     * @return The inverse tangent of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTan/">ArcTan</a>
      */
-    public Complex acosh() {
-        // Define in terms of acos
-        // acosh(z) = +-i acos(z)
-        // Note the special case:
-        // acos(+-0 + iNaN) = π/2 + iNaN
-        // acosh(0 + iNaN) = NaN + iπ/2
-        // will not appropriately multiply by I to maintain positive imaginary if
-        // acos() imaginary computes as NaN. So do this explicitly.
-        if (Double.isNaN(imaginary) && real == 0) {
-            return new Complex(Double.NaN, PI_OVER_2);
-        }
-        return acos(real, imaginary, (re, im) ->
-            // Set the sign appropriately for real >= 0
-            (negative(im)) ?
-                // Multiply by I
-                new Complex(-im, re) :
-                // Multiply by -I
-                new Complex(im, -re)
-        );
+    public Complex atan() {
+        // Define in terms of atanh
+        // atan(z) = -i atanh(iz)
+        // Multiply this number by I, compute atanh, then multiply by back
+        return atanh(-imaginary, real, Complex::multiplyNegativeI);
     }
 
     /**
      * Returns the
-     * <a href="http://mathworld.wolfram.com/Cosine.html">
-     * cosine</a> of this complex number.
+     * <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
+     * hyperbolic sine</a> of this complex number.
      *
-     * <p>\[ \cos(z) = \frac{1}{2} \left( e^{iz} + e^{-iz} \right) \]
+     * <p>\[ \sinh(z) = \frac{1}{2} \left( e^{z} - e^{-z} \right) \]
      *
-     * <p>This is an even function: \( \cos(z) = \cos(-z) \).
-     * The cosine is an entire function and requires no branch cuts.
+     * <p>The hyperbolic sine of \( z \) is an entire function in the complex plane
+     * and is periodic with respect to the imaginary component with period \( 2\pi i \).
+     * Special cases:
+     *
+     * <ul>
+     * <li>{@code z.conj().sinh() == z.sinh().conj()}.
+     * <li>This is an odd function: \( \sinh(z) = -\sinh(-z) \).
+     * <li>If {@code z} is +0 + i0, returns +0 + i0.
+     * <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).
+     * <li>If {@code z} is +0 + iNaN, returns ±0 + iNaN (where the sign of the real part of the result is unspecified).
+     * <li>If {@code z} is x + i∞ for positive finite x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
+     * <li>If {@code z} is +∞ + iy for positive finite y, returns +∞ cis(y) (see {@link #ofCis(double)}.
+     * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
+     * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
+     * <li>If {@code z} is NaN + i0, returns NaN + i0.
+     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
+     * </ul>
      *
      * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
      *
-     * <p>\[ \cos(x + iy) = \cos(x)\cosh(y) - i \sin(x)\sinh(y) \]
+     * <p>\[ \sinh(x + iy) = \sinh(x)\cos(y) + i \cosh(x)\sin(y) \]
      *
-     * <p>As per the C99 standard this function is computed using the trigonomic identity:
+     * @return The hyperbolic sine of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sinh/">Sinh</a>
+     */
+    public Complex sinh() {
+        return sinh(real, imaginary, Complex::ofCartesian);
+    }
+
+    /**
+     * Returns the hyperbolic sine of the complex number.
      *
-     * <p>\[ cos(z) = cosh(iz) \]
+     * <p>This function exists to allow implementation of the identity
+     * {@code sin(z) = -i sinh(iz)}.<p>
      *
-     * @return The cosine of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cos/">Cos</a>
+     * @param real Real part.
+     * @param imaginary Imaginary part.
+     * @param constructor Constructor.
+     * @return The hyperbolic sine of the complex number.
      */
-    public Complex cos() {
-        // Define in terms of cosh
-        // cos(z) = cosh(iz)
-        // Multiply this number by I and compute cosh.
-        return cosh(-imaginary, real, Complex::ofCartesian);
+    private static Complex sinh(double real, double imaginary, ComplexConstructor constructor) {
+        if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
+            return constructor.create(real, Double.NaN);
+        }
+        if (real == 0) {
+            // Imaginary-only sinh(iy) = i sin(y).
+            if (Double.isFinite(imaginary)) {
+                // Maintain periodic property with respect to the imaginary component.
+                // sinh(+/-0.0) * cos(+/-x) = +/-0 * cos(x)
+                return constructor.create(changeSign(real, Math.cos(imaginary)),
+                                          Math.sin(imaginary));
+            }
+            // If imaginary is inf/NaN the sign of the real part is unspecified.
+            // Returning the same real value maintains the conjugate equality.
+            // It is not possible to also maintain the odd function (hence the unspecified sign).
+            return constructor.create(real, Double.NaN);
+        }
+        if (imaginary == 0) {
+            // Real-only sinh(x).
+            return constructor.create(Math.sinh(real), imaginary);
+        }
+        final double x = Math.abs(real);
+        if (x > SAFE_EXP) {
+            // Approximate sinh/cosh(x) using exp^|x| / 2
+            return coshsinh(x, real, imaginary, true, constructor);
+        }
+        // No overflow of sinh/cosh
+        return constructor.create(Math.sinh(real) * Math.cos(imaginary),
+                                  Math.cosh(real) * Math.sin(imaginary));
     }
 
     /**
@@ -2145,1025 +2385,727 @@ public final class Complex implements Serializable  {
 
     /**
      * Returns the
-     * <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
-     * exponential function</a> of this complex number.
+     * <a href="http://mathworld.wolfram.com/HyperbolicTangent.html">
+     * hyperbolic tangent</a> of this complex number.
      *
-     * <p>\[ \exp(z) = e^z \]
+     * <p>\[ \tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}} \]
      *
-     * <p>The exponential function of \( z \) is an entire function in the complex plane.
-     * Special cases:
+     * <p>The hyperbolic tangent of \( z \) is an entire function in the complex plane
+     * and is periodic with respect to the imaginary component with period \( \pi i \)
+     * and has poles of the first order along the imaginary line, at coordinates
+     * \( (0, \pi(\frac{1}{2} + n)) \).
+     * Note that the {@code double} floating-point representation is unable to exactly represent
+     * \( \pi/2 \) and there is no value for which a pole error occurs. Special cases:
      *
      * <ul>
-     * <li>{@code z.conj().exp() == z.exp().conj()}.
-     * <li>If {@code z} is ±0 + i0, returns 1 + i0.
-     * <li>If {@code z} is x + i∞ for finite x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
-     * <li>If {@code z} is −∞ + iy for finite y, returns +0 cis(y) (see {@link #ofCis(double)}).
-     * <li>If {@code z} is +∞ + iy for finite nonzero y, returns +∞ cis(y).
-     * <li>If {@code z} is −∞ + i∞, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
-     * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
-     * <li>If {@code z} is −∞ + iNaN, returns ±0 ± i0 (where the signs of the real and imaginary parts of the result are unspecified).
-     * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
+     * <li>{@code z.conj().tanh() == z.tanh().conj()}.
+     * <li>This is an odd function: \( \tanh(z) = -\tanh(-z) \).
+     * <li>If {@code z} is +0 + i0, returns +0 + i0.
+     * <li>If {@code z} is 0 + i∞, returns 0 + iNaN.
+     * <li>If {@code z} is x + i∞ for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is 0 + iNaN, returns 0 + iNAN.
+     * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns 1 + i0 sin(2y).
+     * <li>If {@code z} is +∞ + i∞, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
+     * <li>If {@code z} is +∞ + iNaN, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
      * <li>If {@code z} is NaN + i0, returns NaN + i0.
      * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
      * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
      * </ul>
      *
-     * <p>Implements the formula:
+     * <p>Special cases include the technical corrigendum
+     * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
+     * DR 471: Complex math functions cacosh and ctanh</a>.
      *
-     * <p>\[ \exp(x + iy) = e^x (\cos(y) + i \sin(y)) \]
+     * <p>This is defined using real \( x \) and imaginary \( y \) parts:
      *
-     * @return The exponential of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</a>
+     * <p>\[ \tan(x + iy) = \frac{\sinh(2x)}{\cosh(2x)+\cos(2y)} + i \frac{\sin(2y)}{\cosh(2x)+\cos(2y)} \]
+     *
+     * <p>The implementation uses double-angle identities to avoid overflow of {@code 2x}
+     * and {@code 2y}.
+     *
+     * @return The hyperbolic tangent of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
      */
-    public Complex exp() {
-        if (Double.isInfinite(real)) {
-            // Set the scale factor applied to cis(y)
-            double zeroOrInf;
-            if (real < 0) {
-                if (!Double.isFinite(imaginary)) {
-                    // (−∞ + i∞) or (−∞ + iNaN) returns (±0 ± i0) (where the signs of the
-                    // real and imaginary parts of the result are unspecified).
-                    // Here we preserve the conjugate equality.
-                    return new Complex(0, Math.copySign(0, imaginary));
-                }
-                // (−∞ + iy) returns +0 cis(y), for finite y
-                zeroOrInf = 0;
-            } else {
-                // (+∞ + i0) returns +∞ + i0.
-                if (imaginary == 0) {
-                    return this;
-                }
-                // (+∞ + i∞) or (+∞ + iNaN) returns (±∞ + iNaN) and raises the invalid
-                // floating-point exception (where the sign of the real part of the
-                // result is unspecified).
-                if (!Double.isFinite(imaginary)) {
-                    return new Complex(real, Double.NaN);
+    public Complex tanh() {
+        return tanh(real, imaginary, Complex::ofCartesian);
+    }
+
+    /**
+     * Returns the hyperbolic tangent of this complex number.
+     *
+     * <p>This function exists to allow implementation of the identity
+     * {@code tan(z) = -i tanh(iz)}.<p>
+     *
+     * @param real Real part.
+     * @param imaginary Imaginary part.
+     * @param constructor Constructor.
+     * @return The hyperbolic tangent of the complex number.
+     */
+    private static Complex tanh(double real, double imaginary, ComplexConstructor constructor) {
+        // Cache the absolute real value
+        final double x = Math.abs(real);
+
+        // Handle inf or nan.
+        if (!isPosFinite(x) || !Double.isFinite(imaginary)) {
+            if (isPosInfinite(x)) {
+                if (Double.isFinite(imaginary)) {
+                    // The sign is copied from sin(2y)
+                    // The identity sin(2a) = 2 sin(a) cos(a) is used for consistency
+                    // with the computation below. Only the magnitude is important
+                    // so drop the 2. When |y| is small sign(sin(2y)) = sign(y).
+                    final double sign = Math.abs(imaginary) < PI_OVER_2 ?
+                                        imaginary :
+                                        Math.sin(imaginary) * Math.cos(imaginary);
+                    return constructor.create(Math.copySign(1, real),
+                                              Math.copySign(0, sign));
                 }
-                // (+∞ + iy) returns (+∞ cis(y)), for finite nonzero y.
-                zeroOrInf = real;
+                // imaginary is infinite or NaN
+                return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
             }
-            return new Complex(zeroOrInf * Math.cos(imaginary),
-                               zeroOrInf * Math.sin(imaginary));
-        } else if (Double.isNaN(real)) {
-            // (NaN + i0) returns (NaN + i0)
-            // (NaN + iy) returns (NaN + iNaN) and optionally raises the invalid floating-point exception
-            // (NaN + iNaN) returns (NaN + iNaN)
-            return imaginary == 0 ? this : NAN;
-        } else if (!Double.isFinite(imaginary)) {
-            // (x + i∞) or (x + iNaN) returns (NaN + iNaN) and raises the invalid
-            // floating-point exception, for finite x.
-            return NAN;
+            // Remaining cases:
+            // (0 + i inf), returns (0 + i NaN)
+            // (0 + i NaN), returns (0 + i NaN)
+            // (x + i inf), returns (NaN + i NaN) for non-zero x (including infinite)
+            // (x + i NaN), returns (NaN + i NaN) for non-zero x (including infinite)
+            // (NaN + i 0), returns (NaN + i 0)
+            // (NaN + i y), returns (NaN + i NaN) for non-zero y (including infinite)
+            // (NaN + i NaN), returns (NaN + i NaN)
+            return constructor.create(real == 0 ? real : Double.NaN,
+                                      imaginary == 0 ? imaginary : Double.NaN);
         }
-        // real and imaginary are finite.
-        // Compute e^a * (cos(b) + i sin(b)).
 
-        // Special case:
-        // (±0 + i0) returns (1 + i0)
-        final double exp = Math.exp(real);
+        // Finite components
+        // tanh(x+iy) = (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y))
+
+        if (real == 0) {
+            // Imaginary-only tanh(iy) = i tan(y)
+            // Identity: sin 2y / (1 + cos 2y) = tan(y)
+            return constructor.create(real, Math.tan(imaginary));
+        }
         if (imaginary == 0) {
-            return new Complex(exp, imaginary);
+            // Identity: sinh 2x / (1 + cosh 2x) = tanh(x)
+            return constructor.create(Math.tanh(real), imaginary);
         }
-        return new Complex(exp * Math.cos(imaginary),
-                           exp * Math.sin(imaginary));
+
+        // The double angles can be avoided using the identities:
+        // sinh(2x) = 2 sinh(x) cosh(x)
+        // sin(2y) = 2 sin(y) cos(y)
+        // cosh(2x) = 2 sinh^2(x) + 1
+        // cos(2y) = 2 cos^2(y) - 1
+        // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
+        // To avoid a junction when swapping between the double angles and the identities
+        // the identities are used in all cases.
+
+        if (x > SAFE_EXP / 2) {
+            // Potential overflow in sinh/cosh(2x).
+            // Approximate sinh/cosh using exp^x.
+            // Ignore cos^2(y) in the divisor as it is insignificant.
+            // real = sinh(x)cosh(x) / sinh^2(x) = +/-1
+            final double re = Math.copySign(1, real);
+            // imag = sin(2y) / 2 sinh^2(x)
+            // sinh(x) -> sign(x) * e^|x| / 2 when x is large.
+            // sinh^2(x) -> e^2|x| / 4 when x is large.
+            // imag = sin(2y) / 2 (e^2|x| / 4) = 2 sin(2y) / e^2|x|
+            //      = 4 * sin(y) cos(y) / e^2|x|
+            // Underflow safe divide as e^2|x| may overflow:
+            // imag = 4 * sin(y) cos(y) / e^m / e^(2|x| - m)
+            // (|im| is a maximum of 2)
+            double im = Math.sin(imaginary) * Math.cos(imaginary);
+            if (x > SAFE_EXP) {
+                // e^2|x| > e^m * e^m
+                // This will underflow 2.0 / e^m / e^m
+                im = Math.copySign(0.0, im);
+            } else {
+                // e^2|x| = e^m * e^(2|x| - m)
+                im = 4 * im / EXP_M / Math.exp(2 * x - SAFE_EXP);
+            }
+            return constructor.create(re, im);
+        }
+
+        // No overflow of sinh(2x) and cosh(2x)
+
+        // Note: This does not use the definitional formula but uses the identity:
+        // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
+
+        final double sinhx = Math.sinh(real);
+        final double coshx = Math.cosh(real);
+        final double siny = Math.sin(imaginary);
+        final double cosy = Math.cos(imaginary);
+        final double divisor = sinhx * sinhx + cosy * cosy;
+        return constructor.create(sinhx * coshx / divisor,
+                                  siny * cosy / divisor);
     }
 
     /**
      * Returns the
-     * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html">
-     * natural logarithm</a> of this complex number.
+     * <a href="http://mathworld.wolfram.com/InverseHyperbolicSine.html">
+     * inverse hyperbolic sine</a> of this complex number.
      *
-     * <p>The natural logarithm of \( z \) is unbounded along the real axis and
-     * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
-     * natural logarithm has a branch cut along the negative real axis \( (-infty,0] \).
-     * Special cases:
+     * <p>\[ \sinh^{-1}(z) = \ln \left(z + \sqrt{1 + z^2} \right) \]
+     *
+     * <p>The inverse hyperbolic sine of \( z \) is unbounded along the real axis and
+     * in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
      *
      * <ul>
-     * <li>{@code z.conj().log() == z.log().conj()}.
-     * <li>If {@code z} is −0 + i0, returns −∞ + iπ ("divide-by-zero" floating-point operation).
-     * <li>If {@code z} is +0 + i0, returns −∞ + i0 ("divide-by-zero" floating-point operation).
-     * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
+     * <li>{@code z.conj().asinh() == z.asinh().conj()}.
+     * <li>This is an odd function: \( \sinh^{-1}(z) = -\sinh^{-1}(-z) \).
+     * <li>If {@code z} is +0 + i0, returns 0 + i0.
+     * <li>If {@code z} is x + i∞ for positive-signed finite x, returns +∞ + iπ/2.
      * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +∞ + iπ.
-     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
-     * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
+     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
      * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
-     * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
-     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
+     * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
+     * <li>If {@code z} is NaN + i0, returns NaN + i0.
+     * <li>If {@code z} is NaN + iy for finite nonzero y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
      * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
      * </ul>
      *
-     * <p>Implements the formula:
-     *
-     * <p>\[ \ln(z) = \ln |z| + i \arg(z) \]
+     * <p>The inverse hyperbolic sine is a multivalued function and requires a branch cut in
+     * the complex plane; the cut is conventionally placed at the line segments
+     * \( (-i \infty,-i) \) and \( (i,i \infty) \) of the imaginary axis.
      *
-     * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
+     * <p>This function is computed using the trigonomic identity:
      *
-     * <p>The implementation is based on the method described in:</p>
-     * <blockquote>
-     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
-     * Implementing complex elementary functions using exception handling.
-     * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
-     * </blockquote>
+     * <p>\[ \sinh^{-1}(z) = -i \sin^{-1}(iz) \]
      *
-     * @return The natural logarithm of this complex number.
-     * @see Math#log(double)
-     * @see #abs()
-     * @see #arg()
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
+     * @return The inverse hyperbolic sine of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSinh/">ArcSinh</a>
      */
-    public Complex log() {
-        return log(Math::log, HALF, LN_2, Complex::ofCartesian);
+    public Complex asinh() {
+        // Define in terms of asin
+        // asinh(z) = -i asin(iz)
+        // Note: This is the opposite to the identity defined in the C99 standard:
+        // asin(z) = -i asinh(iz)
+        // Multiply this number by I, compute asin, then multiply by back
+        return asin(-imaginary, real, Complex::multiplyNegativeI);
     }
 
     /**
-     * Returns the base 10
-     * <a href="http://mathworld.wolfram.com/CommonLogarithm.html">
-     * common logarithm</a> of this complex number.
-     *
-     * <p>The common logarithm of \( z \) is unbounded along the real axis and
-     * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the
-     * common logarithm has a branch cut along the negative real axis \( (-infty,0] \).
-     * Special cases are as defined in the {@link #log() natural logarithm}:
-     *
-     * <p>Implements the formula:
+     * Returns the
+     * <a href="http://mathworld.wolfram.com/InverseHyperbolicCosine.html">
+     * inverse hyperbolic cosine</a> of this complex number.
      *
-     * <p>\[ \log_{10}(z) = \log_{10} |z| + i \arg(z) \]
+     * <p>\[ \cosh^{-1}(z) = \ln \left(z + \sqrt{z + 1} \sqrt{z - 1} \right) \]
      *
-     * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
-     *
-     * @return The base 10 logarithm of this complex number.
-     * @see Math#log10(double)
-     * @see #abs()
-     * @see #arg()
-     */
-    public Complex log10() {
-        return log(Math::log10, LOG_10E_O_2, LOG10_2, Complex::ofCartesian);
-    }
-
-    /**
-     * Returns the logarithm of this complex number using the provided function.
-     * Implements the formula:
-     *
-     * <pre>
-     *   log(x + i y) = log(|x + i y|) + i arg(x + i y)</pre>
-     *
-     * <p>Warning: The argument {@code logOf2} must be equal to {@code log(2)} using the
-     * provided log function otherwise scaling using powers of 2 in the case of overflow
-     * will be incorrect. This is provided as an internal optimisation.
-     *
-     * @param log Log function.
-     * @param logOfeOver2 The log function applied to e, then divided by 2.
-     * @param logOf2 The log function applied to 2.
-     * @param constructor Constructor for the returned complex.
-     * @return The logarithm of this complex number.
-     * @see #abs()
-     * @see #arg()
-     */
-    private Complex log(DoubleUnaryOperator log, double logOfeOver2, double logOf2, ComplexConstructor constructor) {
-        // Handle NaN
-        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
-            // Return NaN unless infinite
-            if (isInfinite()) {
-                return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
-            }
-            // No-use of the input constructor
-            return NAN;
-        }
-
-        // Returns the real part:
-        // log(sqrt(x^2 + y^2))
-        // log(x^2 + y^2) / 2
-
-        // Compute with positive values
-        double x = Math.abs(real);
-        double y = Math.abs(imaginary);
-
-        // Find the larger magnitude.
-        if (x < y) {
-            final double tmp = x;
-            x = y;
-            y = tmp;
-        }
-
-        if (x == 0) {
-            // Handle zero: raises the ‘‘divide-by-zero’’ floating-point exception.
-            return constructor.create(Double.NEGATIVE_INFINITY,
-                                      negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary);
-        }
-
-        double re;
-
-        // This alters the implementation of Hull et al (1994) which used a standard
-        // precision representation of |z|: sqrt(x*x + y*y).
-        // This formula should use the same definition of the magnitude returned
-        // by Complex.abs() which is a high precision computation with scaling.
-        // The checks for overflow thus only require ensuring the output of |z|
-        // will not overflow or underflow.
-
-        if (x > HALF && x < ROOT2) {
-            // x^2+y^2 close to 1. Use log1p(x^2+y^2 - 1) / 2.
-            re = Math.log1p(x2y2m1(x, y)) * logOfeOver2;
-        } else {
-            // Check for over/underflow in |z|
-            // When scaling:
-            // log(a / b) = log(a) - log(b)
-            // So initialise the result with the log of the scale factor.
-            re = 0;
-            if (x > Double.MAX_VALUE / 2) {
-                // Potential overflow.
-                if (isPosInfinite(x)) {
-                    // Handle infinity
-                    return constructor.create(x, arg());
-                }
-                // Scale down.
-                x /= 2;
-                y /= 2;
-                // log(2)
-                re = logOf2;
-            } else if (y < Double.MIN_NORMAL) {
-                // Potential underflow.
-                if (y == 0) {
-                    // Handle real only number
-                    return constructor.create(log.applyAsDouble(x), arg());
-                }
-                // Scale up sub-normal numbers to make them normal by scaling by 2^54,
-                // i.e. more than the mantissa digits.
-                x *= 0x1.0p54;
-                y *= 0x1.0p54;
-                // log(2^-54) = -54 * log(2)
-                re = -54 * logOf2;
-            }
-            re += log.applyAsDouble(abs(x, y));
-        }
-
-        // All ISO C99 edge cases for the imaginary are satisfied by the Math library.
-        return constructor.create(re, arg());
-    }
-
-    /**
-     * Compute {@code x^2 + y^2 - 1} in high precision.
-     * Assumes that the values x and y can be multiplied without overflow; that
-     * {@code x >= y}; and both values are positive.
-     *
-     * @param x the x value
-     * @param y the y value
-     * @return {@code x^2 + y^2 - 1}.
-     */
-    private static double x2y2m1(double x, double y) {
-        // Hull et al used (x-1)*(x+1)+y*y.
-        // From the paper on page 236:
-
-        // If x == 1 there is no cancellation.
-
-        // If x > 1, there is also no cancellation, but the argument is now accurate
-        // only to within a factor of 1 + 3 EPSILSON (note that x – 1 is exact),
-        // so that error = 3 EPSILON.
-
-        // If x < 1, there can be serious cancellation:
-
-        // If 4 y^2 < |x^2 – 1| the cancellation is not serious ... the argument is accurate
-        // only to within a factor of 1 + 4 EPSILSON so that error = 4 EPSILON.
-
-        // Otherwise there can be serious cancellation and the relative error in the real part
-        // could be enormous.
-
-        final double xx = x * x;
-        final double yy = y * y;
-        // Modify to use high precision before the threshold set by Hull et al.
-        // This is to preserve the monotonic output of the computation at the switch.
-        // Set the threshold when x^2 + y^2 is above 0.5 thus subtracting 1 results in a number
-        // that can be expressed with a higher precision than any number in the range 0.5-1.0
-        // due to the variable exponent used below 0.5.
-        if (x < 1 && xx + yy > 0.5) {
-            // Large relative error.
-            // This does not use o.a.c.numbers.LinearCombination.value(x, x, y, y, 1, -1).
-            // It is optimised knowing that:
-            // - the products are squares
-            // - the final term is -1 (which does not require split multiplication and addition)
-            // - The answer will not be NaN as the terms are not NaN components
-            // - The order is known to be 1 > |x| >= |y|
-            // The squares are computed using a split multiply algorithm and
-            // the summation using an extended precision summation algorithm.
-
-            // Split x and y as one 26 bits number and one 27 bits number
-            final double xHigh = splitHigh(x);
-            final double xLow  = x - xHigh;
-            final double yHigh = splitHigh(y);
-            final double yLow  = y - yHigh;
-
-            // Accurate split multiplication x * x and y * y
-            final double x2Low = squareLow(xLow, xHigh, xx);
-            final double y2Low = squareLow(yLow, yHigh, yy);
-
-            return sumx2y2m1(xx, x2Low, yy, y2Low);
-        }
-        return (x - 1) * (x + 1) + yy;
-    }
-
-    /**
-     * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) create
-     * a big value from which to derive the two split parts.
-     * <pre>
-     * c = (2^s + 1) * a
-     * a_big = c - a
-     * a_hi = c - a_big
-     * a_lo = a - a_hi
-     * a = a_hi + a_lo
-     * </pre>
-     *
-     * <p>The multiplicand must be odd allowing a p-bit value to be split into
-     * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.
-     * Combined they have (p􏰔-1) bits of significand but the sign bit of {@code a_lo}
-     * contains a bit of information.
-     *
-     * @param a Value.
-     * @return the high part of the value.
-     * @see <a href="https://doi.org/10.1007/BF01397083">
-     * Dekker (1971) A floating-point technique for extending the available precision</a>
-     */
-    private static double splitHigh(double a) {
-        final double c = MULTIPLIER * a;
-        return c - (c - a);
-    }
-
-    /**
-     * Compute the round-off from the square of a split number with {@code low} and {@code high}
-     * components. Uses Dekker's algorithm for split multiplication modified for a square product.
-     *
-     * <p>Note: This is candidate to be replaced with {@code Math.fma(x, x, -x * x)} to compute
-     * the round-off from the square product {@code x * x}. This would remove the requirement
-     * to compute the split number and make this method redundant. {@code Math.fma} requires
-     * JDK 9 and FMA hardware support.
-     *
-     * @param low Low part of number.
-     * @param high High part of number.
-     * @param square Square of the number.
-     * @return <code>low * low - (((product - high * high) - low * high) - high * low)</code>
-     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
-     * Shewchuk (1997) Theorum 18</a>
-     */
-    private static double squareLow(double low, double high, double square) {
-        final double lh = low * high;
-        return low * low - (((square - high * high) - lh) - lh);
-    }
-
-    /**
-     * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
-     * Dekker's two-sum algorithm. The values are required to be ordered by magnitude:
-     * {@code |a| >= |b|}.
-     *
-     * @param a First part of sum.
-     * @param b Second part of sum.
-     * @param x Sum.
-     * @return <code>b - (x - a)</code>
-     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
-     * Shewchuk (1997) Theorum 6</a>
-     */
-    private static double fastSumLow(double a, double b, double x) {
-        // x = a + b
-        // bVirtual = x - a
-        // y = b - bVirtual
-        return b - (x - a);
-    }
-
-    /**
-     * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
-     * Knuth's two-sum algorithm. The values are not required to be ordered by magnitude.
-     *
-     * @param a First part of sum.
-     * @param b Second part of sum.
-     * @param x Sum.
-     * @return <code>(a - (x - (x - a))) + (b - (x - a))</code>
-     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
-     * Shewchuk (1997) Theorum 7</a>
-     */
-    private static double sumLow(double a, double b, double x) {
-        // x = a + b
-        // bVirtual = x - a
-        // aVirtual = x - bVirtual
-        // bRoundoff = b - bVirtual
-        // aRoundoff = a - aVirtual
-        // y = aRoundoff + bRoundoff
-        final double bVirtual = x - a;
-        return (a - (x - bVirtual)) + (b - bVirtual);
-    }
-
-    /**
-     * Sum x^2 + y^2 - 1. It is assumed that {@code y <= x < 1}.
-     *
-     * <p>Implement Shewchuk's expansion-sum algorithm: [x2Low, x2High] + [-1] + [y2Low, y2High].
-     *
-     * @param x2High High part of x^2.
-     * @param x2Low Low part of x^2.
-     * @param y2High High part of y^2.
-     * @param y2Low Low part of y^2.
-     * @return x^2 + y^2 - 1
-     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
-     * Shewchuk (1997) Theorum 12</a>
-     */
-    private static double sumx2y2m1(double x2High, double x2Low, double y2High, double y2Low) {
-        // Let e and f be non-overlapping expansions of components of length m and n.
-        // The following algorithm will produce a non-overlapping expansion h where the
-        // sum h_i = e + f and components of h are in increasing order of magnitude.
-
-        // Expansion-sum proceeds by a grow-expansion of the first part from one expansion
-        // into the other, extending its length by 1. The process repeats for the next part
-        // but the grow-expansion starts at the previous merge position + 1.
-        // Thus expansion-sum requires mn two-sum operations to merge length m into length n
-        // resulting in length m+n-1.
-
-        // Variables numbered from 1 as per Figure 7 (p.12). The output expansion h is placed
-        // into e increasing its length for each grow expansion.
-
-        // We have two expansions for x^2 and y^2 and the whole number -1.
-        // Expecting (x^2 + y^2) close to 1 we generate first the intermediate expansion
-        // (x^2 - 1) moving the result away from 1 where there are sparse floating point
-        // representations. This is then added to a similar magnitude y^2. Leaving the -1
-        // until last suffers from 1 ulp rounding errors more often and the requirement
-        // for a distillation sum to reduce rounding error frequency.
-
-        // Note: Do not use the alternative fast-expansion-sum of the parts sorted by magnitude.
-        // The parts can be ordered with a single comparison into:
-        // [y2Low, (y2High|x2Low), x2High, -1]
-        // The fast-two-sum saves 1 fast-two-sum and 3 two-sum operations (21 additions) and
-        // adds a penalty of a single branch condition.
-        // However the order in not "strongly non-overlapping" and the fast-expansion-sum
-        // output will not be strongly non-overlapping. The sum of the output has 1 ulp error
-        // on random cis numbers approximately 1 in 160 events. This can be removed by a
-        // distillation two-sum pass over the final expansion as a cost of 1 fast-two-sum and
-        // 3 two-sum operations! So we use the expansion sum with the same operations and
-        // no branches.
-
-        // q=running sum
-        double q = x2Low - 1;
-        double e1 = fastSumLow(-1, x2Low, q);
-        double e3 = q + x2High;
-        double e2 = sumLow(q, x2High, e3);
-
-        final double f1 = y2Low;
-        final double f2 = y2High;
-
-        // Grow expansion of f1 into e
-        q = f1 + e1;
-        e1 = sumLow(f1, e1, q);
-        double p = q + e2;
-        e2 = sumLow(q, e2, p);
-        double e4 = p + e3;
-        e3 = sumLow(p, e3, e4);
-
-        // Grow expansion of f2 into e (only required to start at e2)
-        q = f2 + e2;
-        e2 = sumLow(f2, e2, q);
-        p = q + e3;
-        e3 = sumLow(q, e3, p);
-        final double e5 = p + e4;
-        e4 = sumLow(p, e4, e5);
-
-        // Final summation:
-        // The sum of the parts is within 1 ulp of the true expansion value e:
-        // |e - sum| < ulp(sum).
-        // To achieve the exact result requires iteration of a distillation two-sum through
-        // the expansion until convergence, i.e. no smaller term changes higher terms.
-        // This requires (n-1) iterations for length n. Here we neglect this as
-        // although the method is not ensured to be exact is it robust on random
-        // cis numbers.
-        return e1 + e2 + e3 + e4 + e5;
-    }
-
-    /**
-     * Returns the complex power of this complex number raised to the power of {@code x}.
-     * Implements the formula:
-     *
-     * <p>\[ z^x = e^{x \ln(z)} \]
-     *
-     * <p>If this complex number is zero then this method returns zero if {@code x} is positive
-     * in the real component and zero in the imaginary component;
-     * otherwise it returns NaN + iNaN.
-     *
-     * @param  x The exponent to which this complex number is to be raised.
-     * @return This complex number raised to the power of {@code x}.
-     * @see #log()
-     * @see #multiply(Complex)
-     * @see #exp()
-     * @see <a href="http://mathworld.wolfram.com/ComplexExponentiation.html">Complex exponentiation</a>
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
-     */
-    public Complex pow(Complex x) {
-        if (real == 0 &&
-            imaginary == 0) {
-            // This value is zero. Test the other.
-            if (x.real > 0 &&
-                x.imaginary == 0) {
-                // 0 raised to positive number is 0
-                return ZERO;
-            }
-            // 0 raised to anything else is NaN
-            return NAN;
-        }
-        return log().multiply(x).exp();
-    }
-
-    /**
-     * Returns the complex power of this complex number raised to the power of {@code x},
-     * with {@code x} interpreted as a real number.
-     * Implements the formula:
-     *
-     * <p>\[ z^x = e^{x \ln(z)} \]
-     *
-     * <p>If this complex number is zero then this method returns zero if {@code x} is positive;
-     * otherwise it returns NaN + iNaN.
-     *
-     * @param  x The exponent to which this complex number is to be raised.
-     * @return This complex number raised to the power of {@code x}.
-     * @see #log()
-     * @see #multiply(double)
-     * @see #exp()
-     * @see #pow(Complex)
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
-     */
-    public Complex pow(double x) {
-        if (real == 0 &&
-            imaginary == 0) {
-            // This value is zero. Test the other.
-            if (x > 0) {
-                // 0 raised to positive number is 0
-                return ZERO;
-            }
-            // 0 raised to anything else is NaN
-            return NAN;
-        }
-        return log().multiply(x).exp();
-    }
-
-    /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/Sine.html">
-     * sine</a> of this complex number.
-     *
-     * <p>\[ \sin(z) = \frac{1}{2} i \left( e^{-iz} - e^{iz} \right) \]
-     *
-     * <p>This is an odd function: \( \sin(z) = -\sin(-z) \).
-     * The sine is an entire function and requires no branch cuts.
-     *
-     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
-     *
-     * <p>\[ \sin(x + iy) = \sin(x)\cosh(y) + i \cos(x)\sinh(y) \]
-     *
-     * <p>As per the C99 standard this function is computed using the trigonomic identity:
-     *
-     * <p>\[ \sin(z) = -i \sinh(iz) \]
-     *
-     * @return The sine of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sin/">Sin</a>
-     */
-    public Complex sin() {
-        // Define in terms of sinh
-        // sin(z) = -i sinh(iz)
-        // Multiply this number by I, compute sinh, then multiply by back
-        return sinh(-imaginary, real, Complex::multiplyNegativeI);
-    }
-
-    /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
-     * hyperbolic sine</a> of this complex number.
-     *
-     * <p>\[ \sinh(z) = \frac{1}{2} \left( e^{z} - e^{-z} \right) \]
-     *
-     * <p>The hyperbolic sine of \( z \) is an entire function in the complex plane
-     * and is periodic with respect to the imaginary component with period \( 2\pi i \).
-     * Special cases:
+     * <p>The inverse hyperbolic cosine of \( z \) is in the range \( [0, \infty) \) along the
+     * real axis and in the range \( [-\pi, \pi] \) along the imaginary axis. Special cases:
      *
      * <ul>
-     * <li>{@code z.conj().sinh() == z.sinh().conj()}.
-     * <li>This is an odd function: \( \sinh(z) = -\sinh(-z) \).
-     * <li>If {@code z} is +0 + i0, returns +0 + i0.
-     * <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).
-     * <li>If {@code z} is +0 + iNaN, returns ±0 + iNaN (where the sign of the real part of the result is unspecified).
-     * <li>If {@code z} is x + i∞ for positive finite x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is x + iNaN for finite nonzero x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is +∞ + i0, returns +∞ + i0.
-     * <li>If {@code z} is +∞ + iy for positive finite y, returns +∞ cis(y) (see {@link #ofCis(double)}.
-     * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified; "invalid" floating-point operation).
-     * <li>If {@code z} is +∞ + iNaN, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
-     * <li>If {@code z} is NaN + i0, returns NaN + i0.
-     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>{@code z.conj().acosh() == z.acosh().conj()}.
+     * <li>If {@code z} is ±0 + i0, returns +0 + iπ/2.
+     * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2.
+     * <li>If {@code z} is 0 + iNaN, returns NaN + iπ/2 <sup>[1]</sup>.
+     * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is −∞ + iy for positive-signed finite y, returns +∞ + iπ.
+     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns +∞ + i0.
+     * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4.
+     * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4.
+     * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN.
+     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN.
      * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
      * </ul>
      *
-     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+     * <p>Special cases include the technical corrigendum
+     * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
+     * DR 471: Complex math functions cacosh and ctanh</a>.
      *
-     * <p>\[ \sinh(x + iy) = \sinh(x)\cos(y) + i \cosh(x)\sin(y) \]
+     * <p>The inverse hyperbolic cosine is a multivalued function and requires a branch cut in
+     * the complex plane; the cut is conventionally placed at the line segment
+     * \( (-\infty,-1) \) of the real axis.
      *
-     * @return The hyperbolic sine of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sinh/">Sinh</a>
-     */
-    public Complex sinh() {
-        return sinh(real, imaginary, Complex::ofCartesian);
-    }
-
-    /**
-     * Returns the hyperbolic sine of the complex number.
+     * <p>This function is computed using the trigonomic identity:
      *
-     * <p>This function exists to allow implementation of the identity
-     * {@code sin(z) = -i sinh(iz)}.<p>
+     * <p>\[ \cosh^{-1}(z) = \pm i \cos^{-1}(z) \]
      *
-     * @param real Real part.
-     * @param imaginary Imaginary part.
-     * @param constructor Constructor.
-     * @return The hyperbolic sine of the complex number.
+     * <p>The sign of the multiplier is chosen to give {@code z.acosh().real() >= 0}
+     * and compatibility with the C99 standard.
+     *
+     * @return The inverse hyperbolic cosine of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCosh/">ArcCosh</a>
      */
-    private static Complex sinh(double real, double imaginary, ComplexConstructor constructor) {
-        if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
-            return constructor.create(real, Double.NaN);
-        }
-        if (real == 0) {
-            // Imaginary-only sinh(iy) = i sin(y).
-            if (Double.isFinite(imaginary)) {
-                // Maintain periodic property with respect to the imaginary component.
-                // sinh(+/-0.0) * cos(+/-x) = +/-0 * cos(x)
-                return constructor.create(changeSign(real, Math.cos(imaginary)),
-                                          Math.sin(imaginary));
-            }
-            // If imaginary is inf/NaN the sign of the real part is unspecified.
-            // Returning the same real value maintains the conjugate equality.
-            // It is not possible to also maintain the odd function (hence the unspecified sign).
-            return constructor.create(real, Double.NaN);
-        }
-        if (imaginary == 0) {
-            // Real-only sinh(x).
-            return constructor.create(Math.sinh(real), imaginary);
-        }
-        final double x = Math.abs(real);
-        if (x > SAFE_EXP) {
-            // Approximate sinh/cosh(x) using exp^|x| / 2
-            return coshsinh(x, real, imaginary, true, constructor);
+    public Complex acosh() {
+        // Define in terms of acos
+        // acosh(z) = +-i acos(z)
+        // Note the special case:
+        // acos(+-0 + iNaN) = π/2 + iNaN
+        // acosh(0 + iNaN) = NaN + iπ/2
+        // will not appropriately multiply by I to maintain positive imaginary if
+        // acos() imaginary computes as NaN. So do this explicitly.
+        if (Double.isNaN(imaginary) && real == 0) {
+            return new Complex(Double.NaN, PI_OVER_2);
         }
-        // No overflow of sinh/cosh
-        return constructor.create(Math.sinh(real) * Math.cos(imaginary),
-                                  Math.cosh(real) * Math.sin(imaginary));
+        return acos(real, imaginary, (re, im) ->
+            // Set the sign appropriately for real >= 0
+            (negative(im)) ?
+                // Multiply by I
+                new Complex(-im, re) :
+                // Multiply by -I
+                new Complex(im, -re)
+        );
     }
 
     /**
      * Returns the
-     * <a href="http://mathworld.wolfram.com/SquareRoot.html">
-     * square root</a> of this complex number.
-     *
-     * <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) \]
+     * <a href="http://mathworld.wolfram.com/InverseHyperbolicTangent.html">
+     * inverse hyperbolic tangent</a> of this complex number.
      *
-     * <p>The square root of \( z \) is in the range \( [0, +\infty) \) along the real axis and
-     * is unbounded along the imaginary axis. The imaginary part of the square root has a
-     * branch cut along the negative real axis \( (-infty,0) \). Special cases:
+     * <p>\[ \tanh^{-1}(z) = \frac{1}{2} \ln \left( \frac{1 + z}{1 - z} \right) \]
      *
-     * <ul>
-     * <li>{@code z.conj().sqrt() == z.sqrt().conj()}.
-     * <li>If {@code z} is ±0 + i0, returns +0 + i0.
-     * <li>If {@code z} is x + i∞ for all x (including NaN), returns +∞ + i∞.
-     * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +0 + i∞.
-     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0.
-     * <li>If {@code z} is −∞ + iNaN, returns NaN ± i∞ (where the sign of the imaginary part of the result is unspecified).
-     * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
-     * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
-     * </ul>
+     * <p>The inverse hyperbolic tangent of \( z \) is unbounded along the real axis and
+     * in the range \( [-\pi/2, \pi/2] \) along the imaginary axis. Special cases:
      *
-     * <p>Implements the following algorithm to compute \( \sqrt{x + iy} \):
-     * <ol>
-     * <li>Let \( t = \sqrt{2 (|x| + |x + iy|)} \)
-     * <li>if \( x \geq 0 \) return \( \frac{t}{2} + i \frac{y}{t} \)
-     * <li>else return \( \frac{|y|}{t} + i\ \text{sgn}(y) \frac{t}{2} \)
-     * </ol>
-     * where:
      * <ul>
-     * <li>\( |x| =\ \){@link Math#abs(double) abs}(x)
-     * <li>\( |x + y i| =\ \){@link Complex#abs}
-     * <li>\( \text{sgn}(y) =\ \){@link Math#copySign(double,double) copySign}(1.0, y)
+     * <li>{@code z.conj().atanh() == z.atanh().conj()}.
+     * <li>This is an odd function: \( \tanh^{-1}(z) = -\tanh^{-1}(-z) \).
+     * <li>If {@code z} is +0 + i0, returns +0 + i0.
+     * <li>If {@code z} is +0 + iNaN, returns +0 + iNaN.
+     * <li>If {@code z} is +1 + i0, returns +∞ + i0 ("divide-by-zero" floating-point operation).
+     * <li>If {@code z} is x + i∞ for finite positive-signed x, returns +0 + iπ/2.
+     * <li>If {@code z} is x+iNaN for nonzero finite x, returns NaN+iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +0 + iπ/2.
+     * <li>If {@code z} is +∞ + i∞, returns +0 + iπ/2.
+     * <li>If {@code z} is +∞ + iNaN, returns +0 + iNaN.
+     * <li>If {@code z} is NaN+iy for finite y, returns NaN+iNaN ("invalid" floating-point operation).
+     * <li>If {@code z} is NaN + i∞, returns ±0 + iπ/2 (where the sign of the real part of the result is unspecified).
+     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
      * </ul>
      *
-     * <p>The implementation is overflow and underflow safe based on the method described in:</p>
-     * <blockquote>
-     * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994)
-     * Implementing complex elementary functions using exception handling.
-     * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
-     * </blockquote>
+     * <p>The inverse hyperbolic tangent is a multivalued function and requires a branch cut in
+     * the complex plane; the cut is conventionally placed at the line segments
+     * \( (\infty,-1] \) and \( [1,\infty) \) of the real axis.
      *
-     * @return The square root of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sqrt/">Sqrt</a>
+     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:
+     *
+     * <p>\[ \tanh^{-1}(z) = \frac{1}{4} \ln \left(1 + \frac{4x}{(1-x)^2+y^2} \right) + \\
+     *                     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) \]
+     *
+     * <p>The imaginary part is computed using {@link Math#atan2(double, double)} to ensure the
+     * correct quadrant is returned from \( \tan^{-1} \left(\frac{2y}{1-x^2-y^2} \right) \).
+     *
+     * <p>The code has been adapted from the <a href="https://www.boost.org/">Boost</a>
+     * {@code c++} implementation {@code <boost/math/complex/atanh.hpp>}.
+     *
+     * @return The inverse hyperbolic tangent of this complex number.
+     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTanh/">ArcTanh</a>
      */
-    public Complex sqrt() {
-        return sqrt(real, imaginary);
+    public Complex atanh() {
+        return atanh(real, imaginary, Complex::ofCartesian);
     }
 
     /**
-     * Returns the square root of the complex number {@code sqrt(x + i y)}.
+     * Returns the inverse hyperbolic tangent of this complex number.
      *
-     * @param real Real component.
-     * @param imaginary Imaginary component.
-     * @return The square root of the complex number.
+     * <p>This function exists to allow implementation of the identity
+     * {@code atan(z) = -i atanh(iz)}.
+     *
+     * The original notice is shown below and the licence is shown in full in LICENSE.txt:
+     * <pre>
+     * (C) Copyright John Maddock 2005.
+     * Distributed under the Boost Software License, Version 1.0. (See accompanying
+     * file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+     * </pre>
+     *
+     * @param real Real part.
+     * @param imaginary Imaginary part.
+     * @param constructor Constructor.
+     * @return The inverse hyperbolic tangent of the complex number.
      */
-    private static Complex sqrt(double real, double imaginary) {
-        // Handle NaN
-        if (Double.isNaN(real) || Double.isNaN(imaginary)) {
-            // Check for infinite
-            if (Double.isInfinite(imaginary)) {
-                return new Complex(Double.POSITIVE_INFINITY, imaginary);
+    private static Complex atanh(final double real, final double imaginary,
+                                 final ComplexConstructor constructor) {
+        // Compute with positive values and determine sign at the end
+        double x = Math.abs(real);
+        double y = Math.abs(imaginary);
+        // The result (without sign correction)
+        double re;
+        double im;
+
+        // Handle C99 special cases
+        if (Double.isNaN(x)) {
+            if (isPosInfinite(y)) {
+                // The sign of the real part of the result is unspecified
+                return constructor.create(0, Math.copySign(PI_OVER_2, imaginary));
             }
-            if (Double.isInfinite(real)) {
-                if (real == Double.NEGATIVE_INFINITY) {
-                    return new Complex(Double.NaN, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
-                }
-                return new Complex(Double.POSITIVE_INFINITY, Double.NaN);
+            // No-use of the input constructor.
+            // Optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
+            return NAN;
+        } else if (Double.isNaN(y)) {
+            if (isPosInfinite(x)) {
+                return constructor.create(Math.copySign(0, real), Double.NaN);
+            }
+            if (x == 0) {
+                return constructor.create(real, Double.NaN);
             }
+            // No-use of the input constructor
             return NAN;
-        }
-
-        // Compute with positive values and determine sign at the end
-        final double x = Math.abs(real);
-        final double y = Math.abs(imaginary);
-
-        // Compute
-        double t;
+        } else {
+            // x && y are finite or infinite.
 
-        // This alters the implementation of Hull et al (1994) which used a standard
-        // precision representation of |z|: sqrt(x*x + y*y).
-        // This formula should use the same definition of the magnitude returned
-        // by Complex.abs() which is a high precision computation with scaling.
-        // Worry about overflow if 2 * (|z| + |x|) will overflow.
-        // Worry about underflow if |z| or |x| are sub-normal components.
+            // Check the safe region.
+            // The lower and upper bounds have been copied from boost::math::atanh.
+            // They are different from the safe region for asin and acos.
+            // x >= SAFE_UPPER: (1-x) == -x
+            // x <= SAFE_LOWER: 1 - x^2 = 1
 
-        if (inRegion(x, y, Double.MIN_NORMAL, SQRT_SAFE_UPPER)) {
-            // No over/underflow
-            t = Math.sqrt(2 * (abs(x, y) + x));
-        } else {
-            // Potential over/underflow. First check infinites and real/imaginary only.
+            if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
+                // Normal computation within a safe region.
 
-            // Check for infinite
-            if (isPosInfinite(y)) {
-                return new Complex(Double.POSITIVE_INFINITY, imaginary);
-            } else if (isPosInfinite(x)) {
-                if (real == Double.NEGATIVE_INFINITY) {
-                    return new Complex(0, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
+                // minus x plus 1: (-x+1)
+                final double mxp1 = 1 - x;
+                final double yy = y * y;
+                // The definition of real component is:
+                // real = log( ((x+1)^2+y^2) / ((1-x)^2+y^2) ) / 4
+                // This simplifies by adding 1 and subtracting 1 as a fraction:
+                //      = log(1 + ((x+1)^2+y^2) / ((1-x)^2+y^2) - ((1-x)^2+y^2)/((1-x)^2+y^2) ) / 4
+                //
+                // real(atanh(z)) == log(1 + 4*x / ((1-x)^2+y^2)) / 4
+                // imag(atanh(z)) == tan^-1 (2y, (1-x)(1+x) - y^2) / 2
+                // imag(atanh(z)) == tan^-1 (2y, (1 - x^2 - y^2) / 2
+                // The division is done at the end of the function.
+                re = Math.log1p(4 * x / (mxp1 * mxp1 + yy));
+                // Modified from boost which does not switch the magnitude of x and y.
+                // The denominator for atan2 is 1 - x^2 - y^2.
+                // This can be made more precise if |x| > |y|.
+                final double numerator = 2 * y;
+                double denominator;
+                if (x < y) {
+                    final double tmp = x;
+                    x = y;
+                    y = tmp;
                 }
-                return new Complex(Double.POSITIVE_INFINITY, Math.copySign(0, imaginary));
-            } else if (y == 0) {
-                // Real only
-                final double sqrtAbs = Math.sqrt(x);
-                if (real < 0) {
-                    return new Complex(0, Math.copySign(sqrtAbs, imaginary));
+                // 1 - x is precise if |x| >= 1
+                if (x >= 1) {
+                    denominator = (1 - x) * (1 + x) - y * y;
+                } else {
+                    // |x| < 1: Use high precision if possible:
+                    // 1 - x^2 - y^2 = -(x^2 + y^2 - 1)
+                    denominator = -x2y2m1(x, y);
                 }
-                return new Complex(sqrtAbs, imaginary);
-            } else if (x == 0) {
-                // Imaginary only. This sets the two components to the same magnitude.
-                // Note: In polar coordinates this does not happen:
-                // real = sqrt(abs()) * Math.cos(arg() / 2)
-                // imag = sqrt(abs()) * Math.sin(arg() / 2)
-                // arg() / 2 = pi/4 and cos and sin should both return sqrt(2)/2 but
-                // are different by 1 ULP.
-                final double sqrtAbs = Math.sqrt(y) * ONE_OVER_ROOT2;
-                return new Complex(sqrtAbs, Math.copySign(sqrtAbs, imaginary));
+                im = Math.atan2(numerator, denominator);
             } else {
-                // Over/underflow.
-                // Full scaling is not required as this is done in the hypotenuse function.
-                // Keep the number as big as possible for maximum precision in the second sqrt.
-                // Note if we scale by an even power of 2, we can re-scale by sqrt of the number.
-                // a = sqrt(b)
-                // a = sqrt(b/4) * sqrt(4)
+                // This section handles exception cases that would normally cause
+                // underflow or overflow in the main formulas.
 
-                double rescale;
-                double sx;
-                double sy;
-                if (Math.max(x, y) > SQRT_SAFE_UPPER) {
-                    // Overflow. Scale down by 16 and rescale by sqrt(16).
-                    sx = x / 16;
-                    sy = y / 16;
-                    rescale = 4;
+                // C99. G.7: Special case for imaginary only numbers
+                if (x == 0) {
+                    if (imaginary == 0) {
+                        return constructor.create(real, imaginary);
+                    }
+                    // atanh(iy) = i atan(y)
+                    return constructor.create(real, Math.atan(imaginary));
+                }
+
+                // Real part:
+                // real = Math.log1p(4x / ((1-x)^2 + y^2))
+                // real = Math.log1p(4x / (1 - 2x + x^2 + y^2))
+                // real = Math.log1p(4x / (1 + x(x-2) + y^2))
+                // without either overflow or underflow in the squared terms.
+                if (x >= SAFE_UPPER) {
+                    // (1-x) = -x to machine precision:
+                    // log1p(4x / (x^2 + y^2))
+                    if (isPosInfinite(x) || isPosInfinite(y)) {
+                        re = 0;
+                    } else if (y >= SAFE_UPPER) {
+                        // Big x and y: divide by x*y
+                        re = Math.log1p((4 / y) / (x / y + y / x));
+                    } else if (y > 1) {
+                        // Big x: divide through by x:
+                        re = Math.log1p(4 / (x + y * y / x));
+                    } else {
+                        // Big x small y, as above but neglect y^2/x:
+                        re = Math.log1p(4 / x);
+                    }
+                } else if (y >= SAFE_UPPER) {
+                    if (x > 1) {
+                        // Big y, medium x, divide through by y:
+                        final double mxp1 = 1 - x;
+                        re = Math.log1p((4 * x / y) / (mxp1 * mxp1 / y + y));
+                    } else {
+                        // Big y, small x, as above but neglect (1-x)^2/y:
+                        // Note: log1p(v) == v - v^2/2 + v^3/3 ... Taylor series when v is small.
+                        // Here v is so small only the first term matters.
+                        re = 4 * x / y / y;
+                    }
+                } else if (x == 1) {
+                    // x = 1, small y:
+                    // Special case when x == 1 as (1-x) is invalid.
+                    // Simplify the following formula:
+                    // real = log( sqrt((x+1)^2+y^2) ) / 2 - log( sqrt((1-x)^2+y^2) ) / 2
+                    //      = log( sqrt(4+y^2) ) / 2 - log(y) / 2
+                    // if: 4+y^2 -> 4
+                    //      = log( 2 ) / 2 - log(y) / 2
+                    //      = (log(2) - log(y)) / 2
+                    // Multiply by 2 as it will be divided by 4 at the end.
+                    // C99: if y=0 raises the ‘‘divide-by-zero’’ floating-point exception.
+                    re = 2 * (LN_2 - Math.log(y));
                 } else {
-                    // Sub-normal numbers. Make them normal by scaling by 2^54,
-                    // i.e. more than the mantissa digits, and rescale by sqrt(2^54) = 2^27.
-                    sx = x * 0x1.0p54;
-                    sy = y * 0x1.0p54;
-                    rescale = 0x1.0p-27;
+                    // Modified from boost which checks y > SAFE_LOWER.
+                    // if y*y -> 0 it will be ignored so always include it.
+                    final double mxp1 = 1 - x;
+                    re = Math.log1p((4 * x) / (mxp1 * mxp1 + y * y));
+                }
+
+                // Imaginary part:
+                // imag = atan2(2y, (1-x)(1+x) - y^2)
+                // if x or y are large, then the formula:
+                //   atan2(2y, (1-x)(1+x) - y^2)
+                // evaluates to +(PI - theta) where theta is negligible compared to PI.
+                if ((x >= SAFE_UPPER) || (y >= SAFE_UPPER)) {
+                    im = Math.PI;
+                } else if (x <= SAFE_LOWER) {
+                    // (1-x)^2 -> 1
+                    if (y <= SAFE_LOWER) {
+                        // 1 - y^2 -> 1
+                        im = Math.atan2(2 * y, 1);
+                    } else {
+                        im = Math.atan2(2 * y, 1 - y * y);
+                    }
+                } else {
+                    // Medium x, small y.
+                    // Modified from boost which checks (y == 0) && (x == 1) and sets re = 0.
+                    // This is same as the result from calling atan2(0, 0) so exclude this case.
+                    // 1 - y^2 = 1 so ignore subtracting y^2
+                    im = Math.atan2(2 * y, (1 - x) * (1 + x));
                 }
-                t = rescale * Math.sqrt(2 * (abs(sx, sy) + sx));
             }
         }
 
-        if (real >= 0) {
-            return new Complex(t / 2, imaginary / t);
-        }
-        return new Complex(y / t, Math.copySign(t / 2, imaginary));
+        re /= 4;
+        im /= 2;
+        return constructor.create(changeSign(re, real),
+                                  changeSign(im, imaginary));
     }
 
     /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/Tangent.html">
-     * tangent</a> of this complex number.
-     *
-     * <p>\[ \tan(z) = \frac{i(e^{-iz} - e^{iz})}{e^{-iz} + e^{iz}} \]
-     *
-     * <p>This is an odd function: \( \tan(z) = -\tan(-z) \).
-     * The tangent is an entire function and requires no branch cuts.
-     *
-     * <p>This is implemented using real \( x \) and imaginary \( y \) parts:</p>
-     * \[ \tan(x + iy) = \frac{\sin(2x)}{\cos(2x)+\cosh(2y)} + i \frac{\sinh(2y)}{\cos(2x)+\cosh(2y)} \]
-     *
-     * <p>As per the C99 standard this function is computed using the trigonomic identity:</p>
-     * \[ \tan(z) = -i \tanh(iz) \]
+     * Compute {@code x^2 + y^2 - 1} in high precision.
+     * Assumes that the values x and y can be multiplied without overflow; that
+     * {@code x >= y}; and both values are positive.
      *
-     * @return The tangent of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tan/">Tangent</a>
+     * @param x the x value
+     * @param y the y value
+     * @return {@code x^2 + y^2 - 1}.
      */
-    public Complex tan() {
-        // Define in terms of tanh
-        // tan(z) = -i tanh(iz)
-        // Multiply this number by I, compute tanh, then multiply by back
-        return tanh(-imaginary, real, Complex::multiplyNegativeI);
-    }
+    private static double x2y2m1(double x, double y) {
+        // Hull et al used (x-1)*(x+1)+y*y.
+        // From the paper on page 236:
 
-    /**
-     * Returns the
-     * <a href="http://mathworld.wolfram.com/HyperbolicTangent.html">
-     * hyperbolic tangent</a> of this complex number.
-     *
-     * <p>\[ \tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}} \]
-     *
-     * <p>The hyperbolic tangent of \( z \) is an entire function in the complex plane
-     * and is periodic with respect to the imaginary component with period \( \pi i \)
-     * and has poles of the first order along the imaginary line, at coordinates
-     * \( (0, \pi(\frac{1}{2} + n)) \).
-     * Note that the {@code double} floating-point representation is unable to exactly represent
-     * \( \pi/2 \) and there is no value for which a pole error occurs. Special cases:
-     *
-     * <ul>
-     * <li>{@code z.conj().tanh() == z.tanh().conj()}.
-     * <li>This is an odd function: \( \tanh(z) = -\tanh(-z) \).
-     * <li>If {@code z} is +0 + i0, returns +0 + i0.
-     * <li>If {@code z} is 0 + i∞, returns 0 + iNaN.
-     * <li>If {@code z} is x + i∞ for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is 0 + iNaN, returns 0 + iNAN.
-     * <li>If {@code z} is x + iNaN for finite non-zero x, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is +∞ + iy for positive-signed finite y, returns 1 + i0 sin(2y).
-     * <li>If {@code z} is +∞ + i∞, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
-     * <li>If {@code z} is +∞ + iNaN, returns 1 ± i0 (where the sign of the imaginary part of the result is unspecified).
-     * <li>If {@code z} is NaN + i0, returns NaN + i0.
-     * <li>If {@code z} is NaN + iy for all nonzero numbers y, returns NaN + iNaN ("invalid" floating-point operation).
-     * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN.
-     * </ul>
-     *
-     * <p>Special cases include the technical corrigendum
-     * <a href="http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471">
-     * DR 471: Complex math functions cacosh and ctanh</a>.
-     *
-     * <p>This is defined using real \( x \) and imaginary \( y \) parts:
-     *
-     * <p>\[ \tan(x + iy) = \frac{\sinh(2x)}{\cosh(2x)+\cos(2y)} + i \frac{\sin(2y)}{\cosh(2x)+\cos(2y)} \]
-     *
-     * <p>The implementation uses double-angle identities to avoid overflow of {@code 2x}
-     * and {@code 2y}.
-     *
-     * @return The hyperbolic tangent of this complex number.
-     * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
-     */
-    public Complex tanh() {
-        return tanh(real, imaginary, Complex::ofCartesian);
-    }
+        // If x == 1 there is no cancellation.
 
-    /**
-     * Returns the hyperbolic tangent of this complex number.
-     *
-     * <p>This function exists to allow implementation of the identity
-     * {@code tan(z) = -i tanh(iz)}.<p>
-     *
-     * @param real Real part.
-     * @param imaginary Imaginary part.
-     * @param constructor Constructor.
-     * @return The hyperbolic tangent of the complex number.
-     */
-    private static Complex tanh(double real, double imaginary, ComplexConstructor constructor) {
-        // Cache the absolute real value
-        final double x = Math.abs(real);
+        // If x > 1, there is also no cancellation, but the argument is now accurate
+        // only to within a factor of 1 + 3 EPSILSON (note that x – 1 is exact),
+        // so that error = 3 EPSILON.
 
-        // Handle inf or nan.
-        if (!isPosFinite(x) || !Double.isFinite(imaginary)) {
-            if (isPosInfinite(x)) {
-                if (Double.isFinite(imaginary)) {
-                    // The sign is copied from sin(2y)
-                    // The identity sin(2a) = 2 sin(a) cos(a) is used for consistency
-                    // with the computation below. Only the magnitude is important
-                    // so drop the 2. When |y| is small sign(sin(2y)) = sign(y).
-                    final double sign = Math.abs(imaginary) < PI_OVER_2 ?
-                                        imaginary :
-                                        Math.sin(imaginary) * Math.cos(imaginary);
-                    return constructor.create(Math.copySign(1, real),
-                                              Math.copySign(0, sign));
-                }
-                // imaginary is infinite or NaN
-                return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
-            }
-            // Remaining cases:
-            // (0 + i inf), returns (0 + i NaN)
-            // (0 + i NaN), returns (0 + i NaN)
-            // (x + i inf), returns (NaN + i NaN) for non-zero x (including infinite)
-            // (x + i NaN), returns (NaN + i NaN) for non-zero x (including infinite)
-            // (NaN + i 0), returns (NaN + i 0)
-            // (NaN + i y), returns (NaN + i NaN) for non-zero y (including infinite)
-            // (NaN + i NaN), returns (NaN + i NaN)
-            return constructor.create(real == 0 ? real : Double.NaN,
-                                      imaginary == 0 ? imaginary : Double.NaN);
-        }
+        // If x < 1, there can be serious cancellation:
 
-        // Finite components
-        // tanh(x+iy) = (sinh(2x) + i sin(2y)) / (cosh(2x) + cos(2y))
+        // If 4 y^2 < |x^2 – 1| the cancellation is not serious ... the argument is accurate
+        // only to within a factor of 1 + 4 EPSILSON so that error = 4 EPSILON.
 
-        if (real == 0) {
-            // Imaginary-only tanh(iy) = i tan(y)
-            // Identity: sin 2y / (1 + cos 2y) = tan(y)
-            return constructor.create(real, Math.tan(imaginary));
-        }
-        if (imaginary == 0) {
-            // Identity: sinh 2x / (1 + cosh 2x) = tanh(x)
-            return constructor.create(Math.tanh(real), imaginary);
-        }
+        // Otherwise there can be serious cancellation and the relative error in the real part
+        // could be enormous.
 
-        // The double angles can be avoided using the identities:
-        // sinh(2x) = 2 sinh(x) cosh(x)
-        // sin(2y) = 2 sin(y) cos(y)
-        // cosh(2x) = 2 sinh^2(x) + 1
-        // cos(2y) = 2 cos^2(y) - 1
-        // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
-        // To avoid a junction when swapping between the double angles and the identities
-        // the identities are used in all cases.
+        final double xx = x * x;
+        final double yy = y * y;
+        // Modify to use high precision before the threshold set by Hull et al.
+        // This is to preserve the monotonic output of the computation at the switch.
+        // Set the threshold when x^2 + y^2 is above 0.5 thus subtracting 1 results in a number
+        // that can be expressed with a higher precision than any number in the range 0.5-1.0
+        // due to the variable exponent used below 0.5.
+        if (x < 1 && xx + yy > 0.5) {
+            // Large relative error.
+            // This does not use o.a.c.numbers.LinearCombination.value(x, x, y, y, 1, -1).
+            // It is optimised knowing that:
+            // - the products are squares
+            // - the final term is -1 (which does not require split multiplication and addition)
+            // - The answer will not be NaN as the terms are not NaN components
+            // - The order is known to be 1 > |x| >= |y|
+            // The squares are computed using a split multiply algorithm and
+            // the summation using an extended precision summation algorithm.
 
-        if (x > SAFE_EXP / 2) {
-            // Potential overflow in sinh/cosh(2x).
-            // Approximate sinh/cosh using exp^x.
-            // Ignore cos^2(y) in the divisor as it is insignificant.
-            // real = sinh(x)cosh(x) / sinh^2(x) = +/-1
-            final double re = Math.copySign(1, real);
-            // imag = sin(2y) / 2 sinh^2(x)
-            // sinh(x) -> sign(x) * e^|x| / 2 when x is large.
-            // sinh^2(x) -> e^2|x| / 4 when x is large.
-            // imag = sin(2y) / 2 (e^2|x| / 4) = 2 sin(2y) / e^2|x|
-            //      = 4 * sin(y) cos(y) / e^2|x|
-            // Underflow safe divide as e^2|x| may overflow:
-            // imag = 4 * sin(y) cos(y) / e^m / e^(2|x| - m)
-            // (|im| is a maximum of 2)
-            double im = Math.sin(imaginary) * Math.cos(imaginary);
-            if (x > SAFE_EXP) {
-                // e^2|x| > e^m * e^m
-                // This will underflow 2.0 / e^m / e^m
-                im = Math.copySign(0.0, im);
-            } else {
-                // e^2|x| = e^m * e^(2|x| - m)
-                im = 4 * im / EXP_M / Math.exp(2 * x - SAFE_EXP);
-            }
-            return constructor.create(re, im);
-        }
+            // Split x and y as one 26 bits number and one 27 bits number
+            final double xHigh = splitHigh(x);
+            final double xLow  = x - xHigh;
+            final double yHigh = splitHigh(y);
+            final double yLow  = y - yHigh;
 
-        // No overflow of sinh(2x) and cosh(2x)
+            // Accurate split multiplication x * x and y * y
+            final double x2Low = squareLow(xLow, xHigh, xx);
+            final double y2Low = squareLow(yLow, yHigh, yy);
 
-        // Note: This does not use the definitional formula but uses the identity:
-        // tanh(x+iy) = (sinh(x)cosh(x) + i sin(y)cos(y)) / (sinh^2(x) + cos^2(y))
+            return sumx2y2m1(xx, x2Low, yy, y2Low);
+        }
+        return (x - 1) * (x + 1) + yy;
+    }
 
-        final double sinhx = Math.sinh(real);
-        final double coshx = Math.cosh(real);
-        final double siny = Math.sin(imaginary);
-        final double cosy = Math.cos(imaginary);
-        final double divisor = sinhx * sinhx + cosy * cosy;
-        return constructor.create(sinhx * coshx / divisor,
-                                  siny * cosy / divisor);
+    /**
+     * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) create
+     * a big value from which to derive the two split parts.
+     * <pre>
+     * c = (2^s + 1) * a
+     * a_big = c - a
+     * a_hi = c - a_big
+     * a_lo = a - a_hi
+     * a = a_hi + a_lo
+     * </pre>
+     *
+     * <p>The multiplicand must be odd allowing a p-bit value to be split into
+     * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.
+     * Combined they have (p􏰔-1) bits of significand but the sign bit of {@code a_lo}
+     * contains a bit of information.
+     *
+     * @param a Value.
+     * @return the high part of the value.
+     * @see <a href="https://doi.org/10.1007/BF01397083">
+     * Dekker (1971) A floating-point technique for extending the available precision</a>
+     */
+    private static double splitHigh(double a) {
+        final double c = MULTIPLIER * a;
+        return c - (c - a);
+    }
+
+    /**
+     * Compute the round-off from the square of a split number with {@code low} and {@code high}
+     * components. Uses Dekker's algorithm for split multiplication modified for a square product.
+     *
+     * <p>Note: This is candidate to be replaced with {@code Math.fma(x, x, -x * x)} to compute
+     * the round-off from the square product {@code x * x}. This would remove the requirement
+     * to compute the split number and make this method redundant. {@code Math.fma} requires
+     * JDK 9 and FMA hardware support.
+     *
+     * @param low Low part of number.
+     * @param high High part of number.
+     * @param square Square of the number.
+     * @return <code>low * low - (((product - high * high) - low * high) - high * low)</code>
+     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+     * Shewchuk (1997) Theorum 18</a>
+     */
+    private static double squareLow(double low, double high, double square) {
+        final double lh = low * high;
+        return low * low - (((square - high * high) - lh) - lh);
     }
 
     /**
-     * Returns the argument of this complex number.
+     * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
+     * Dekker's two-sum algorithm. The values are required to be ordered by magnitude:
+     * {@code |a| >= |b|}.
      *
-     * <p>The argument is the angle phi between the positive real axis and
-     * the point representing this number in the complex plane.
-     * The value returned is between \( -\pi \) (not inclusive)
-     * and \( \pi \) (inclusive), with negative values returned for numbers with
-     * negative imaginary parts.
+     * @param a First part of sum.
+     * @param b Second part of sum.
+     * @param x Sum.
+     * @return <code>b - (x - a)</code>
+     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+     * Shewchuk (1997) Theorum 6</a>
+     */
+    private static double fastSumLow(double a, double b, double x) {
+        // x = a + b
+        // bVirtual = x - a
+        // y = b - bVirtual
+        return b - (x - a);
+    }
+
+    /**
+     * Compute the round-off from the sum of two numbers {@code a} and {@code b} using
+     * Knuth's two-sum algorithm. The values are not required to be ordered by magnitude.
      *
-     * <p>If either real or imaginary part (or both) is NaN, then the result is NaN.
-     * Infinite parts are handled as {@linkplain Math#atan2} handles them,
-     * essentially treating finite parts as zero in the presence of an
-     * infinite coordinate and returning a multiple of \( \frac{\pi}{4} \) depending on
-     * the signs of the infinite parts.
+     * @param a First part of sum.
+     * @param b Second part of sum.
+     * @param x Sum.
+     * @return <code>(a - (x - (x - a))) + (b - (x - a))</code>
+     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+     * Shewchuk (1997) Theorum 7</a>
+     */
+    private static double sumLow(double a, double b, double x) {
+        // x = a + b
+        // bVirtual = x - a
+        // aVirtual = x - bVirtual
+        // bRoundoff = b - bVirtual
+        // aRoundoff = a - aVirtual
+        // y = aRoundoff + bRoundoff
+        final double bVirtual = x - a;
+        return (a - (x - bVirtual)) + (b - bVirtual);
+    }
+
+    /**
+     * Sum x^2 + y^2 - 1. It is assumed that {@code y <= x < 1}.
      *
-     * <p>This code follows the
-     * <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
-     * in calculating the returned value using the {@code atan2(y, x)} method for complex
-     * \( x + iy \).
+     * <p>Implement Shewchuk's expansion-sum algorithm: [x2Low, x2High] + [-1] + [y2Low, y2High].
      *
-     * @return The argument of this complex number.
-     * @see Math#atan2(double, double)
+     * @param x2High High part of x^2.
+     * @param x2Low Low part of x^2.
+     * @param y2High High part of y^2.
+     * @param y2Low Low part of y^2.
+     * @return x^2 + y^2 - 1
+     * @see <a href="http://www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps">
+     * Shewchuk (1997) Theorum 12</a>
      */
-    public double arg() {
-        // Delegate
-        return Math.atan2(imaginary, real);
+    private static double sumx2y2m1(double x2High, double x2Low, double y2High, double y2Low) {
+        // Let e and f be non-overlapping expansions of components of length m and n.
+        // The following algorithm will produce a non-overlapping expansion h where the
+        // sum h_i = e + f and components of h are in increasing order of magnitude.
+
+        // Expansion-sum proceeds by a grow-expansion of the first part from one expansion
+        // into the other, extending its length by 1. The process repeats for the next part
+        // but the grow-expansion starts at the previous merge position + 1.
+        // Thus expansion-sum requires mn two-sum operations to merge length m into length n
+        // resulting in length m+n-1.
+
+        // Variables numbered from 1 as per Figure 7 (p.12). The output expansion h is placed
+        // into e increasing its length for each grow expansion.
+
+        // We have two expansions for x^2 and y^2 and the whole number -1.
+        // Expecting (x^2 + y^2) close to 1 we generate first the intermediate expansion
+        // (x^2 - 1) moving the result away from 1 where there are sparse floating point
+        // representations. This is then added to a similar magnitude y^2. Leaving the -1
+        // until last suffers from 1 ulp rounding errors more often and the requirement
+        // for a distillation sum to reduce rounding error frequency.
+
+        // Note: Do not use the alternative fast-expansion-sum of the parts sorted by magnitude.
+        // The parts can be ordered with a single comparison into:
+        // [y2Low, (y2High|x2Low), x2High, -1]
+        // The fast-two-sum saves 1 fast-two-sum and 3 two-sum operations (21 additions) and
+        // adds a penalty of a single branch condition.
+        // However the order in not "strongly non-overlapping" and the fast-expansion-sum
+        // output will not be strongly non-overlapping. The sum of the output has 1 ulp error
+        // on random cis numbers approximately 1 in 160 events. This can be removed by a
+        // distillation two-sum pass over the final expansion as a cost of 1 fast-two-sum and
+        // 3 two-sum operations! So we use the expansion sum with the same operations and
+        // no branches.
+
+        // q=running sum
+        double q = x2Low - 1;
+        double e1 = fastSumLow(-1, x2Low, q);
+        double e3 = q + x2High;
+        double e2 = sumLow(q, x2High, e3);
+
+        final double f1 = y2Low;
+        final double f2 = y2High;
+
+        // Grow expansion of f1 into e
+        q = f1 + e1;
+        e1 = sumLow(f1, e1, q);
+        double p = q + e2;
+        e2 = sumLow(q, e2, p);
+        double e4 = p + e3;
+        e3 = sumLow(p, e3, e4);
+
+        // Grow expansion of f2 into e (only required to start at e2)
+        q = f2 + e2;
+        e2 = sumLow(f2, e2, q);
+        p = q + e3;
+        e3 = sumLow(q, e3, p);
+        final double e5 = p + e4;
+        e4 = sumLow(p, e4, e5);
+
+        // Final summation:
+        // The sum of the parts is within 1 ulp of the true expansion value e:
+        // |e - sum| < ulp(sum).
+        // To achieve the exact result requires iteration of a distillation two-sum through
+        // the expansion until convergence, i.e. no smaller term changes higher terms.
+        // This requires (n-1) iterations for length n. Here we neglect this as
+        // although the method is not ensured to be exact is it robust on random
+        // cis numbers.
+        return e1 + e2 + e3 + e4 + e5;
     }
 
     /**
@@ -3210,6 +3152,86 @@ public final class Complex implements Serializable  {
     }
 
     /**
+     * Test for equality with another object. If the other object is a {@code Complex} then a
+     * comparison is made of the real and imaginary parts; otherwise {@code false} is returned.
+     *
+     * <p>If both the real and imaginary parts of two complex numbers
+     * are exactly the same the two {@code Complex} objects are considered to be equal.
+     * For this purpose, two {@code double} values are considered to be
+     * the same if and only if the method {@link Double #doubleToLongBits(double)}
+     * returns the identical {@code long} value when applied to each.
+     *
+     * <p>Note that in most cases, for two instances of class
+     * {@code Complex}, {@code c1} and {@code c2}, the
+     * value of {@code c1.equals(c2)} is {@code true} if and only if
+     *
+     * <pre>
+     *  {@code c1.getReal() == c2.getReal() && c1.getImaginary() == c2.getImaginary()}</pre>
+     *
+     * <p>also has the value {@code true}. However, there are exceptions:
+     *
+     * <ul>
+     *  <li>
+     *   Instances that contain {@code NaN} values in the same part
+     *   are considered to be equal for that part, even though {@code Double.NaN == Double.NaN}
+     *   has the value {@code false}.
+     *  </li>
+     *  <li>
+     *   Instances that share a {@code NaN} value in one part
+     *   but have different values in the other part are <em>not</em> considered equal.
+     *  </li>
+     *  <li>
+     *   Instances that contain different representations of zero in the same part
+     *   are <em>not</em> considered to be equal for that part, even though {@code -0.0 == 0.0}
+     *   has the value {@code true}.
+     *  </li>
+     * </ul>
+     *
+     * <p>The behavior is the same as if the components of the two complex numbers were passed
+     * to {@link java.util.Arrays#equals(double[], double[]) Arrays.equals(double[], double[])}:
+     *
+     * <pre>
+     *  Arrays.equals(new double[]{c1.getReal(), c1.getImaginary()},
+     *                new double[]{c2.getReal(), c2.getImaginary()}); </pre>
+     *
+     * @param other Object to test for equality with this instance.
+     * @return {@code true} if the objects are equal, {@code false} if object
+     * is {@code null}, not an instance of {@code Complex}, or not equal to
+     * this instance.
+     * @see java.lang.Double#doubleToLongBits(double)
+     * @see java.util.Arrays#equals(double[], double[])
+     */
+    @Override
+    public boolean equals(Object other) {
+        if (this == other) {
+            return true;
+        }
+        if (other instanceof Complex) {
+            final Complex c = (Complex) other;
+            return equals(real, c.real) &&
+                equals(imaginary, c.imaginary);
+        }
+        return false;
+    }
+
+    /**
+     * Gets a hash code for the complex number.
+     *
+     * <p>The behavior is the same as if the components of the complex number were passed
+     * to {@link java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])}:
+     *
+     * <pre>
+     *  {@code Arrays.hashCode(new double[] {getReal(), getImaginary()})}</pre>
+     *
+     * @return A hash code value for this object.
+     * @see java.util.Arrays#hashCode(double[]) Arrays.hashCode(double[])
+     */
+    @Override
+    public int hashCode() {
+        return 31 * (31 + Double.hashCode(real)) + Double.hashCode(imaginary);
+    }
+
+    /**
      * Returns a string representation of the complex number.
      *
      * <p>The string will represent the numeric values of the real and imaginary parts.
@@ -3425,28 +3447,6 @@ public final class Complex implements Serializable  {
     }
 
     /**
-     * Creates an exception.
-     *
-     * @param message Message prefix.
-     * @param error Input that caused the error.
-     * @param cause Underlying exception (if any).
-     * @return A new instance.
-     */
-    private static NumberFormatException parsingException(String message,
-                                                          Object error,
-                                                          Throwable cause) {
-        // Not called with a null message or error
-        final StringBuilder sb = new StringBuilder(100)
-            .append(message)
-            .append(" '").append(error).append('\'');
-        if (cause != null) {
-            sb.append(": ").append(cause.getMessage());
-        }
-
-        return new NumberFormatException(sb.toString());
-    }
-
-    /**
      * Returns {@code sqrt(x^2 + y^2)} without intermediate overflow or underflow.
      *
      * <p>Special cases:


Mime
View raw message