You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2020/04/01 10:58:59 UTC
[commons-numbers] 05/08: Rearrange Complex source code.
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 <ah...@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: