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 2022/07/25 10:07:50 UTC
[commons-numbers] 01/03: NUMBERS-188: Refactor complex unuary functions to static methods
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch complex-gsoc-2022
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit ab614de1b4a8d6f1da2624ba52d3b910a1ef5dd6
Author: Sumanth Rajkumar <ra...@gmail.com>
AuthorDate: Mon Jul 25 10:31:46 2022 +0100
NUMBERS-188: Refactor complex unuary functions to static methods
Refactor complex-to-double and complex-to-boolean functions as static methods.
---
.../apache/commons/numbers/complex/Complex.java | 1159 +----------
.../commons/numbers/complex/ComplexFunctions.java | 2183 ++++++++++++++++++--
.../commons/numbers/complex/CReferenceTest.java | 51 +-
.../commons/numbers/complex/CStandardTest.java | 691 +++----
.../numbers/complex/ComplexEdgeCaseTest.java | 276 +--
.../commons/numbers/complex/ComplexTest.java | 125 +-
.../apache/commons/numbers/complex/TestUtils.java | 63 +
.../checkstyle/checkstyle-suppressions.xml | 1 +
8 files changed, 2634 insertions(+), 1915 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 5bb38cfd..4dc47117 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
@@ -82,100 +82,13 @@ public final class Complex implements Serializable {
/** A complex number representing {@code NaN + i NaN}. */
private static final Complex NAN = new Complex(Double.NaN, Double.NaN);
- /** π/2. */
- private static final double PI_OVER_2 = 0.5 * Math.PI;
- /** π/4. */
- private static final double PI_OVER_4 = 0.25 * Math.PI;
- /** Natural logarithm of 2 (ln(2)). */
- private static final double LN_2 = Math.log(2);
- /** {@code 1.0 / sqrt(2)}.
- * This is pre-computed to the closest double from the exact result.
- * It is 1 ULP different from 1.0 / Math.sqrt(2) but equal to Math.sqrt(2) / 2.
- */
- private static final double ONE_OVER_ROOT2 = 0.7071067811865476;
/** Exponent offset in IEEE754 representation. */
private static final int EXPONENT_OFFSET = 1023;
- /**
- * Largest double-precision floating-point number such that
- * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
- * bound on the relative error due to rounding real numbers to double
- * precision floating-point numbers.
- *
- * <p>In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
- * Copied from o.a.c.numbers.Precision.
- *
- * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
- */
- private static final double EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
/** Mask to remove the sign bit from a long. */
private static final long UNSIGN_MASK = 0x7fff_ffff_ffff_ffffL;
/** Mask to extract the 52-bit mantissa from a long representation of a double. */
private static final long MANTISSA_MASK = 0x000f_ffff_ffff_ffffL;
- /**
- * Crossover point to switch computation for asin/acos factor A.
- * This has been updated from the 1.5 value used by Hull et al to 10
- * as used in boost::math::complex.
- * @see <a href="https://svn.boost.org/trac/boost/ticket/7290">Boost ticket 7290</a>
- */
- private static final double A_CROSSOVER = 10.0;
- /** Crossover point to switch computation for asin/acos factor B. */
- private static final double B_CROSSOVER = 0.6471;
- /**
- * The safe maximum double value {@code x} to avoid loss of precision in asin/acos.
- * Equal to sqrt(M) / 8 in Hull, et al (1997) with M the largest normalised floating-point value.
- */
- private static final double SAFE_MAX = Math.sqrt(Double.MAX_VALUE) / 8;
- /**
- * The safe minimum double value {@code x} to avoid loss of precision/underflow in asin/acos.
- * Equal to sqrt(u) * 4 in Hull, et al (1997) with u the smallest normalised floating-point value.
- */
- private static final double SAFE_MIN = Math.sqrt(Double.MIN_NORMAL) * 4;
- /**
- * The safe maximum double value {@code x} to avoid loss of precision in atanh.
- * Equal to sqrt(M) / 2 with M the largest normalised floating-point value.
- */
- private static final double SAFE_UPPER = Math.sqrt(Double.MAX_VALUE) / 2;
- /**
- * The safe minimum double value {@code x} to avoid loss of precision/underflow in atanh.
- * Equal to sqrt(u) * 2 with u the smallest normalised floating-point value.
- */
- private static final double SAFE_LOWER = Math.sqrt(Double.MIN_NORMAL) * 2;
- /** The safe maximum double value {@code x} to avoid overflow in sqrt. */
- private static final double SQRT_SAFE_UPPER = Double.MAX_VALUE / 8;
- /**
- * A safe maximum double value {@code m} where {@code e^m} is not infinite.
- * This can be used when functions require approximations of sinh(x) or cosh(x)
- * when x is large using exp(x):
- * <pre>
- * sinh(x) = (e^x - e^-x) / 2 = sign(x) * e^|x| / 2
- * cosh(x) = (e^x + e^-x) / 2 = e^|x| / 2 </pre>
- *
- * <p>This value can be used to approximate e^x using a product:
- *
- * <pre>
- * e^x = product_n (e^m) * e^(x-nm)
- * n = (int) x/m
- * e.g. e^2000 = e^m * e^m * e^(2000 - 2m) </pre>
- *
- * <p>The value should be below ln(max_value) ~ 709.783.
- * The value m is set to an integer for less error when subtracting m and chosen as
- * even (m=708) as it is used as a threshold in tanh with m/2.
- *
- * <p>The value is used to compute e^x multiplied by a small number avoiding
- * overflow (sinh/cosh) or a small number divided by e^x without underflow due to
- * infinite e^x (tanh). The following conditions are used:
- * <pre>
- * 0.5 * e^m * Double.MIN_VALUE * e^m * e^m = Infinity
- * 2.0 / e^m / e^m = 0.0 </pre>
- */
- private static final double SAFE_EXP = 708;
- /**
- * The value of Math.exp(SAFE_EXP): e^708.
- * To be used in overflow/underflow safe products of e^m to approximate e^x where x > m.
- */
- private static final double EXP_M = Math.exp(SAFE_EXP);
-
/** Serializable version identifier. */
private static final long serialVersionUID = 20180201L;
@@ -204,22 +117,6 @@ public final class Complex implements Serializable {
/** The real part. */
private final double real;
- /**
- * Define a constructor for a Complex.
- * This is used in functions that implement trigonomic identities.
- */
- @FunctionalInterface
- private interface ComplexConstructor {
- /**
- * Create a complex number given the real and imaginary parts.
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @return {@code Complex} object.
- */
- Complex create(double real, double imaginary);
- }
-
/**
* Private default constructor.
*
@@ -484,33 +381,6 @@ public final class Complex implements Serializable {
* @see <a href="http://mathworld.wolfram.com/ComplexModulus.html">Complex modulus</a>
*/
public double abs() {
- return abs(real, imaginary);
- }
-
- /**
- * Returns the absolute value of the complex number.
- * <pre>abs(x + i y) = sqrt(x^2 + y^2)</pre>
- *
- * <p>This should satisfy the special cases of the hypot function in ISO C99 F.9.4.3:
- * "The hypot functions compute the square root of the sum of the squares of x and y,
- * without undue overflow or underflow."
- *
- * <ul>
- * <li>hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent.
- * <li>hypot(x, ±0) is equivalent to |x|.
- * <li>hypot(±∞, y) returns +∞, even if y is a NaN.
- * </ul>
- *
- * <p>This method is called by all methods that require the absolute value of the complex
- * number, e.g. abs(), sqrt() and log().
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @return The absolute value.
- */
- private static double abs(double real, double imaginary) {
- // Specialised implementation of hypot.
- // See NUMBERS-143
return ComplexFunctions.abs(real, imaginary);
}
@@ -565,10 +435,7 @@ public final class Complex implements Serializable {
* @see <a href="http://mathworld.wolfram.com/AbsoluteSquare.html">Absolute square</a>
*/
public double norm() {
- if (isInfinite()) {
- return Double.POSITIVE_INFINITY;
- }
- return real * real + imaginary * imaginary;
+ return ComplexFunctions.norm(real, imaginary);
}
/**
@@ -588,10 +455,7 @@ public final class Complex implements Serializable {
* @see #equals(Object) Complex.equals(Object)
*/
public boolean isNaN() {
- if (Double.isNaN(real) || Double.isNaN(imaginary)) {
- return !isInfinite();
- }
- return false;
+ return ComplexFunctions.isNaN(real, imaginary);
}
/**
@@ -614,7 +478,7 @@ public final class Complex implements Serializable {
* @see Double#isFinite(double)
*/
public boolean isFinite() {
- return Double.isFinite(real) && Double.isFinite(imaginary);
+ return ComplexFunctions.isFinite(real, imaginary);
}
/**
@@ -652,7 +516,7 @@ public final class Complex implements Serializable {
* <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:
+ * <p>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>
*
@@ -1217,55 +1081,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</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);
- }
- // (+∞ + 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));
+ return ComplexFunctions.exp(real, imaginary, Complex::ofCartesian);
}
/**
@@ -1455,107 +1271,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sqrt/">Sqrt</a>
*/
public Complex sqrt() {
- return sqrt(real, imaginary);
- }
-
- /**
- * Returns the square root of the complex number {@code sqrt(x + i y)}.
- *
- * @param real Real component.
- * @param imaginary Imaginary component.
- * @return The square root 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);
- }
- 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);
-
- // 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));
+ return ComplexFunctions.sqrt(real, imaginary, Complex::ofCartesian);
}
/**
@@ -1580,10 +1296,7 @@ public final class Complex implements Serializable {
* @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);
+ return ComplexFunctions.sin(real, imaginary, Complex::ofCartesian);
}
/**
@@ -1608,10 +1321,7 @@ public final class Complex implements Serializable {
* @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);
+ return ComplexFunctions.cos(real, imaginary, Complex::ofCartesian);
}
/**
@@ -1634,10 +1344,7 @@ public final class Complex implements Serializable {
* @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);
+ return ComplexFunctions.tan(real, imaginary, Complex::ofCartesian);
}
/**
@@ -1679,145 +1386,7 @@ public final class Complex implements Serializable {
* @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>}. This method only (and not
- * invoked methods within) is distributed under the Boost Software License V1.0.
- * The original notice is shown below and the licence is shown in full in LICENSE:
- * <pre>
- * (C) Copyright John Maddock 2005.
- * Distributed under the Boost Software License, Version 1.0. (See accompanying
- * file LICENSE or copy at https://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;
- }
- } else if (Double.isNaN(y)) {
- if (x == 0) {
- re = 0;
- im = y;
- } else if (isPosInfinite(x)) {
- re = y;
- im = x;
- } else {
- // No-use of the input constructor
- return NAN;
- }
- } else if (isPosInfinite(x)) {
- re = isPosInfinite(y) ? PI_OVER_4 : PI_OVER_2;
- im = x;
- } else if (isPosInfinite(y)) {
- re = 0;
- im = y;
- } else {
- // Special case for real numbers:
- if (y == 0 && x <= 1) {
- return constructor.create(Math.asin(real), imaginary);
- }
-
- final double xp1 = x + 1;
- final double xm1 = x - 1;
-
- 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;
-
- if (b <= B_CROSSOVER) {
- re = Math.asin(b);
- } else {
- final double apx = a + x;
- if (x <= 1) {
- re = Math.atan(x / Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))));
- } else {
- re = Math.atan(x / (y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))));
- }
- }
-
- 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 4
- if (y <= (EPSILON * Math.abs(xm1))) {
- if (x < 1) {
- re = Math.asin(x);
- im = y / Math.sqrt(xp1 * (1 - x));
- } else {
- re = PI_OVER_2;
- if ((Double.MAX_VALUE / xp1) > xm1) {
- // xp1 * xm1 won't overflow:
- im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
- } else {
- 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 = PI_OVER_2 - Math.sqrt(y);
- im = Math.sqrt(y);
- } else if (EPSILON * y - 1 >= x) {
- // Possible underflow:
- re = x / y;
- im = LN_2 + Math.log(y);
- } else if (x > 1) {
- re = Math.atan(x / y);
- final double xoy = x / y;
- im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
- } else {
- final double a = Math.sqrt(1 + y * y);
- // Possible underflow:
- re = x / a;
- im = 0.5 * Math.log1p(2 * y * (y + a));
- }
- }
- }
-
- return constructor.create(changeSign(re, real),
- changeSign(im, imaginary));
+ return ComplexFunctions.asin(real, imaginary, Complex::ofCartesian);
}
/**
@@ -1874,143 +1443,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCos/">ArcCos</a>
*/
public Complex acos() {
- return acos(real, imaginary, 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)}.
- *
- * <p>Adapted from {@code <boost/math/complex/acos.hpp>}. This method only (and not
- * invoked methods within) is distributed under the Boost Software License V1.0.
- * The original notice is shown below and the licence is shown in full in LICENSE:
- * <pre>
- * (C) Copyright John Maddock 2005.
- * Distributed under the Boost Software License, Version 1.0. (See accompanying
- * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
- * </pre>
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @param constructor Constructor.
- * @return The inverse cosine of the complex number.
- */
- 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);
- }
- // 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;
-
- 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;
-
- 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);
- }
- }
-
- 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));
- }
- }
- }
-
- return constructor.create(negative(real) ? Math.PI - re : re,
- negative(imaginary) ? im : -im);
+ return ComplexFunctions.acos(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2034,10 +1467,7 @@ public final class Complex implements Serializable {
* @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);
+ return ComplexFunctions.atan(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2076,49 +1506,7 @@ public final class Complex implements Serializable {
* @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 exists to allow implementation of the identity
- * {@code sin(z) = -i sinh(iz)}.<p>
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @param constructor Constructor.
- * @return The hyperbolic sine of the complex number.
- */
- 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));
+ return ComplexFunctions.sinh(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2157,115 +1545,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cosh/">Cosh</a>
*/
public Complex cosh() {
- return cosh(real, imaginary, Complex::ofCartesian);
- }
-
- /**
- * Returns the hyperbolic cosine of the complex number.
- *
- * <p>This function exists to allow implementation of the identity
- * {@code cos(z) = cosh(iz)}.<p>
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @param constructor Constructor.
- * @return The hyperbolic cosine of the complex number.
- */
- private static Complex cosh(double real, double imaginary, ComplexConstructor constructor) {
- // ISO C99: Preserve the even function by mapping to positive
- // f(z) = f(-z)
- if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
- return constructor.create(Math.abs(real), Double.NaN);
- }
- if (real == 0) {
- // Imaginary-only cosh(iy) = cos(y).
- if (Double.isFinite(imaginary)) {
- // Maintain periodic property with respect to the imaginary component.
- // sinh(+/-0.0) * sin(+/-x) = +/-0 * sin(x)
- return constructor.create(Math.cos(imaginary),
- changeSign(real, Math.sin(imaginary)));
- }
- // If imaginary is inf/NaN the sign of the imaginary part is unspecified.
- // Although not required by C99 changing the sign maintains the conjugate equality.
- // It is not possible to also maintain the even function (hence the unspecified sign).
- return constructor.create(Double.NaN, changeSign(real, imaginary));
- }
- if (imaginary == 0) {
- // Real-only cosh(x).
- // Change sign to preserve conjugate equality and even function.
- // sin(+/-0) * sinh(+/-x) = +/-0 * +/-a (sinh is monotonic and same sign)
- // => change the sign of imaginary using real. Handles special case of infinite real.
- // If real is NaN the sign of the imaginary part is unspecified.
- return constructor.create(Math.cosh(real), changeSign(imaginary, real));
- }
- final double x = Math.abs(real);
- if (x > SAFE_EXP) {
- // Approximate sinh/cosh(x) using exp^|x| / 2
- return coshsinh(x, real, imaginary, false, constructor);
- }
- // No overflow of sinh/cosh
- return constructor.create(Math.cosh(real) * Math.cos(imaginary),
- Math.sinh(real) * Math.sin(imaginary));
- }
-
- /**
- * Compute cosh or sinh when the absolute real component |x| is large. In this case
- * cosh(x) and sinh(x) can be approximated by exp(|x|) / 2:
- *
- * <pre>
- * cosh(x+iy) real = (e^|x| / 2) * cos(y)
- * cosh(x+iy) imag = (e^|x| / 2) * sin(y) * sign(x)
- * sinh(x+iy) real = (e^|x| / 2) * cos(y) * sign(x)
- * sinh(x+iy) imag = (e^|x| / 2) * sin(y)
- * </pre>
- *
- * @param x Absolute real component |x|.
- * @param real Real part (x).
- * @param imaginary Imaginary part (y).
- * @param sinh Set to true to compute sinh, otherwise cosh.
- * @param constructor Constructor.
- * @return The hyperbolic sine/cosine of the complex number.
- */
- private static Complex coshsinh(double x, double real, double imaginary, boolean sinh,
- ComplexConstructor constructor) {
- // Always require the cos and sin.
- double re = Math.cos(imaginary);
- double im = Math.sin(imaginary);
- // Compute the correct function
- if (sinh) {
- re = changeSign(re, real);
- } else {
- im = changeSign(im, real);
- }
- // Multiply by (e^|x| / 2).
- // Overflow safe computation since sin/cos can be very small allowing a result
- // when e^x overflows: e^x / 2 = (e^m / 2) * e^m * e^(x-2m)
- if (x > SAFE_EXP * 3) {
- // e^x > e^m * e^m * e^m
- // y * (e^m / 2) * e^m * e^m will overflow when starting with Double.MIN_VALUE.
- // Note: Do not multiply by +inf to safeguard against sin(y)=0.0 which
- // will create 0 * inf = nan.
- re *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
- im *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
- } else {
- // Initial part of (e^x / 2) using (e^m / 2)
- re *= EXP_M / 2;
- im *= EXP_M / 2;
- double xm;
- if (x > SAFE_EXP * 2) {
- // e^x = e^m * e^m * e^(x-2m)
- re *= EXP_M;
- im *= EXP_M;
- xm = x - SAFE_EXP * 2;
- } else {
- // e^x = e^m * e^(x-m)
- xm = x - SAFE_EXP;
- }
- final double exp = Math.exp(xm);
- re *= exp;
- im *= exp;
- }
- return constructor.create(re, im);
+ return ComplexFunctions.cosh(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2313,113 +1593,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
*/
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));
- }
- // 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);
- }
-
- // 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) {
- // Identity: sinh 2x / (1 + cosh 2x) = tanh(x)
- return constructor.create(Math.tanh(real), 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);
+ return ComplexFunctions.tanh(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2459,12 +1633,7 @@ public final class Complex implements Serializable {
* @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);
+ return ComplexFunctions.asinh(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2512,24 +1681,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCosh/">ArcCosh</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)
- );
+ return ComplexFunctions.acosh(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2577,207 +1729,7 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTanh/">ArcTanh</a>
*/
public Complex atanh() {
- return atanh(real, imaginary, Complex::ofCartesian);
- }
-
- /**
- * Returns the inverse hyperbolic tangent of this complex number.
- *
- * <p>This function exists to allow implementation of the identity
- * {@code atan(z) = -i atanh(iz)}.
- *
- * <p>Adapted from {@code <boost/math/complex/atanh.hpp>}. This method only (and not
- * invoked methods within) is distributed under the Boost Software License V1.0.
- * The original notice is shown below and the licence is shown in full in LICENSE:
- * <pre>
- * (C) Copyright John Maddock 2005.
- * Distributed under the Boost Software License, Version 1.0. (See accompanying
- * file LICENSE or copy at https://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 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));
- }
- // 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;
- } 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
-
- if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
- // Normal computation within a safe region.
-
- // 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;
- }
- // 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)
- // Modified from boost to use the custom high precision method.
- 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.
-
- // 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 {
- // 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));
- }
- }
- }
-
- re /= 4;
- im /= 2;
- return constructor.create(changeSign(re, real),
- changeSign(im, imaginary));
- }
-
- /**
- * 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}.
- */
- //TODO - make it private in ComplexFunctions in future
- private static double x2y2m1(double x, double y) {
- return ComplexFunctions.x2y2m1(x, y);
+ return ComplexFunctions.atanh(real, imaginary, Complex::ofCartesian);
}
/**
@@ -2956,68 +1908,6 @@ public final class Complex implements Serializable {
return ComplexFunctions.negative(d);
}
- /**
- * Check that a value is positive infinity. Used to replace {@link Double#isInfinite()}
- * when the input value is known to be positive (i.e. in the case where it has been
- * set using {@link Math#abs(double)}).
- *
- * @param d Value.
- * @return {@code true} if {@code d} is +inf.
- */
- //TODO - make private in ComplexFunctions in future
- private static boolean isPosInfinite(double d) {
- return ComplexFunctions.isPosInfinite(d);
- }
-
- /**
- * Check that an absolute value is finite. Used to replace {@link Double#isFinite(double)}
- * when the input value is known to be positive (i.e. in the case where it has been
- * set using {@link Math#abs(double)}).
- *
- * @param d Value.
- * @return {@code true} if {@code d} is +finite.
- */
- private static boolean isPosFinite(double d) {
- return d <= Double.MAX_VALUE;
- }
-
- /**
- * Create a complex number given the real and imaginary parts, then multiply by {@code -i}.
- * This is used in functions that implement trigonomic identities. It is the functional
- * equivalent of:
- *
- * <pre>
- * z = new Complex(real, imaginary).multiplyImaginary(-1);</pre>
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @return {@code Complex} object.
- */
- private static Complex multiplyNegativeI(double real, double imaginary) {
- return new Complex(imaginary, -real);
- }
-
- /**
- * Change the sign of the magnitude based on the signed value.
- *
- * <p>If the signed value is negative then the result is {@code -magnitude}; otherwise
- * return {@code magnitude}.
- *
- * <p>A signed value of {@code -0.0} is treated as negative. A signed value of {@code NaN}
- * is treated as positive.
- *
- * <p>This is not the same as {@link Math#copySign(double, double)} as this method
- * will change the sign based on the signed value rather than copy the sign.
- *
- * @param magnitude the magnitude
- * @param signedValue the signed value
- * @return magnitude or -magnitude.
- * @see #negative(double)
- */
- private static double changeSign(double magnitude, double signedValue) {
- return negative(signedValue) ? -magnitude : magnitude;
- }
-
/**
* Returns a scale suitable for use with {@link Math#scalb(double, int)} to normalise
* the number to the interval {@code [1, 2)}.
@@ -3106,17 +1996,4 @@ public final class Complex implements Serializable {
// A speed test is required to determine performance.
return Math.max(Math.getExponent(a), Math.getExponent(b));
}
-
- /**
- * Checks if both x and y are in the region defined by the minimum and maximum.
- *
- * @param x x value.
- * @param y y value.
- * @param min the minimum (exclusive).
- * @param max the maximum (exclusive).
- * @return true if inside the region.
- */
- private static boolean inRegion(double x, double y, double min, double max) {
- return x < max && x > min && y < max && y > min;
- }
}
diff --git a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java
index 83c0e5e5..e4d783d1 100644
--- a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java
+++ b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java
@@ -58,21 +58,112 @@ import java.util.function.DoubleUnaryOperator;
*/
public final class ComplexFunctions {
- /** Mask to remove the sign bit from a long. */
- static final long UNSIGN_MASK = 0x7fff_ffff_ffff_ffffL;
-
+ /** π/2. */
+ private static final double PI_OVER_2 = 0.5 * Math.PI;
+ /** π/4. */
+ private static final double PI_OVER_4 = 0.25 * Math.PI;
/** Natural logarithm of 2 (ln(2)). */
private static final double LN_2 = Math.log(2);
- /** {@code 1/2}. */
- private static final double HALF = 0.5;
/** Base 10 logarithm of e divided by 2 (log10(e)/2). */
private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2;
/** Base 10 logarithm of 2 (log10(2)). */
private static final double LOG10_2 = Math.log10(2);
+ /** {@code 1/2}. */
+ private static final double HALF = 0.5;
/** {@code sqrt(2)}. */
private static final double ROOT2 = 1.4142135623730951;
+ /** {@code 1.0 / sqrt(2)}.
+ * This is pre-computed to the closest double from the exact result.
+ * It is 1 ULP different from 1.0 / Math.sqrt(2) but equal to Math.sqrt(2) / 2.
+ */
+ private static final double ONE_OVER_ROOT2 = 0.7071067811865476;
/** The bit representation of {@code -0.0}. */
private static final long NEGATIVE_ZERO_LONG_BITS = Double.doubleToLongBits(-0.0);
+ /** Exponent offset in IEEE754 representation. */
+ private static final int EXPONENT_OFFSET = 1023;
+ /**
+ * Largest double-precision floating-point number such that
+ * {@code 1 + EPSILON} is numerically equal to 1. This value is an upper
+ * bound on the relative error due to rounding real numbers to double
+ * precision floating-point numbers.
+ *
+ * <p>In IEEE 754 arithmetic, this is 2<sup>-53</sup>.
+ * Copied from o.a.c.numbers.Precision.
+ *
+ * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
+ */
+ private static final double EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
+ /** Mask to remove the sign bit from a long. */
+ private static final long UNSIGN_MASK = 0x7fff_ffff_ffff_ffffL;
+ /** The multiplier used to split the double value into hi and low parts. This must be odd
+ * and a value of 2^s + 1 in the range {@code p/2 <= s <= p-1} where p is the number of
+ * bits of precision of the floating point number. Here {@code s = 27}.*/
+ private static final double MULTIPLIER = 1.34217729E8;
+
+ /**
+ * Crossover point to switch computation for asin/acos factor A.
+ * This has been updated from the 1.5 value used by Hull et al to 10
+ * as used in boost::math::complex.
+ * @see <a href="https://svn.boost.org/trac/boost/ticket/7290">Boost ticket 7290</a>
+ */
+ private static final double A_CROSSOVER = 10.0;
+ /** Crossover point to switch computation for asin/acos factor B. */
+ private static final double B_CROSSOVER = 0.6471;
+ /**
+ * The safe maximum double value {@code x} to avoid loss of precision in asin/acos.
+ * Equal to sqrt(M) / 8 in Hull, et al (1997) with M the largest normalised floating-point value.
+ */
+ private static final double SAFE_MAX = Math.sqrt(Double.MAX_VALUE) / 8;
+ /**
+ * The safe minimum double value {@code x} to avoid loss of precision/underflow in asin/acos.
+ * Equal to sqrt(u) * 4 in Hull, et al (1997) with u the smallest normalised floating-point value.
+ */
+ private static final double SAFE_MIN = Math.sqrt(Double.MIN_NORMAL) * 4;
+ /**
+ * The safe maximum double value {@code x} to avoid loss of precision in atanh.
+ * Equal to sqrt(M) / 2 with M the largest normalised floating-point value.
+ */
+ private static final double SAFE_UPPER = Math.sqrt(Double.MAX_VALUE) / 2;
+ /**
+ * The safe minimum double value {@code x} to avoid loss of precision/underflow in atanh.
+ * Equal to sqrt(u) * 2 with u the smallest normalised floating-point value.
+ */
+ private static final double SAFE_LOWER = Math.sqrt(Double.MIN_NORMAL) * 2;
+ /** The safe maximum double value {@code x} to avoid overflow in sqrt. */
+ private static final double SQRT_SAFE_UPPER = Double.MAX_VALUE / 8;
+ /**
+ * A safe maximum double value {@code m} where {@code e^m} is not infinite.
+ * This can be used when functions require approximations of sinh(x) or cosh(x)
+ * when x is large using exp(x):
+ * <pre>
+ * sinh(x) = (e^x - e^-x) / 2 = sign(x) * e^|x| / 2
+ * cosh(x) = (e^x + e^-x) / 2 = e^|x| / 2 </pre>
+ *
+ * <p>This value can be used to approximate e^x using a product:
+ *
+ * <pre>
+ * e^x = product_n (e^m) * e^(x-nm)
+ * n = (int) x/m
+ * e.g. e^2000 = e^m * e^m * e^(2000 - 2m) </pre>
+ *
+ * <p>The value should be below ln(max_value) ~ 709.783.
+ * The value m is set to an integer for less error when subtracting m and chosen as
+ * even (m=708) as it is used as a threshold in tanh with m/2.
+ *
+ * <p>The value is used to compute e^x multiplied by a small number avoiding
+ * overflow (sinh/cosh) or a small number divided by e^x without underflow due to
+ * infinite e^x (tanh). The following conditions are used:
+ * <pre>
+ * 0.5 * e^m * Double.MIN_VALUE * e^m * e^m = Infinity
+ * 2.0 / e^m / e^m = 0.0 </pre>
+ */
+ private static final double SAFE_EXP = 708;
+ /**
+ * The value of Math.exp(SAFE_EXP): e^708.
+ * To be used in overflow/underflow safe products of e^m to approximate e^x where x > m.
+ */
+ private static final double EXP_M = Math.exp(SAFE_EXP);
+
/** 54 shifted 20-bits to align with the exponent of the upper 32-bits of a double. */
private static final int EXP_54 = 0x36_00000;
/** Represents an exponent of 500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */
@@ -87,17 +178,48 @@ public final class ComplexFunctions {
/** 2^-600. */
private static final double TWO_POW_NEG_600 = 0x1.0p-600;
- /** The multiplier used to split the double value into hi and low parts. This must be odd
- * and a value of 2^s + 1 in the range {@code p/2 <= s <= p-1} where p is the number of
- * bits of precision of the floating point number. Here {@code s = 27}.*/
- private static final double MULTIPLIER = 1.34217729E8;
-
/**
* Private constructor for utility class.
*/
private ComplexFunctions() {
}
+ /**
+ * Returns the absolute value of the complex number. This is also called complex norm, modulus,
+ * or magnitude.
+ *
+ * <p>\[ \text{abs}(x + i y) = \sqrt{(x^2 + y^2)} \]
+ *
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>{@code abs(x + iy) == abs(y + ix) == abs(x - iy)}.
+ * <li>If {@code z} is ±∞ + iy for any y, returns +∞.
+ * <li>If {@code z} is x + iNaN for non-infinite x, returns NaN.
+ * <li>If {@code z} is x + i0, returns |x|.
+ * </ul>
+ *
+ * <p>The cases ensure that if either component is infinite then the result is positive
+ * infinity. If either component is NaN and this is not {@link #isInfinite(double, double) infinite} then
+ * the result is NaN.
+ *
+ * <p>This method follows the
+ * <a href="http://www.iso-9899.info/wiki/The_Standard">ISO C Standard</a>, Annex G,
+ * in calculating the returned value without intermediate overflow or underflow.
+ *
+ * <p>The computed result will be within 1 ulp of the exact result.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @return The absolute value.
+ * @see #isInfinite(double, double)
+ * @see #isNaN(double, double)
+ * @see <a href="http://mathworld.wolfram.com/ComplexModulus.html">Complex modulus</a>
+ */
+ public static double abs(double real, double imaginary) {
+ return computeAbs(real, imaginary);
+ }
+
/**
* Returns the absolute value of the complex number.
* <pre>abs(x + i y) = sqrt(x^2 + y^2)</pre>
@@ -115,78 +237,1601 @@ public final class ComplexFunctions {
* <p>This method is called by all methods that require the absolute value of the complex
* number, e.g. abs(), sqrt() and log().
*
- * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @return The absolute value.
+ */
+ private static double computeAbs(double real, double imaginary) {
+ // Specialised implementation of hypot.
+ // See NUMBERS-143
+ return hypot(real, imaginary);
+ }
+
+ /**
+ * Returns the argument of the 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 \).
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \). part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @return The argument of the complex number.
+ * @see Math#atan2(double, double)
+ */
+ public static double arg(double real, double imaginary) {
+ // Delegate
+ return Math.atan2(imaginary, real);
+ }
+
+ /**
+ * Returns the squared norm value of the complex number. This is also called the absolute
+ * square.
+ *
+ * <p>\[ \text{norm}(x + i y) = x^2 + y^2 \]
+ *
+ * <p>If either component is infinite then the result is positive infinity. If either
+ * component is NaN and this is not {@link #isInfinite(double, double) infinite} then the result is NaN.
+ *
+ * <p>Note: This method may not return the same value as the square of {@link #abs(double, double)} as
+ * that method uses an extended precision computation.
+ *
+ * <p>{@code norm()} can be used as a faster alternative than {@code abs()} for ranking by
+ * magnitude. If used for ranking any overflow to infinity will create an equal ranking for
+ * values that may be still distinguished by {@code abs()}.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @return The square norm value.
+ * @see #isInfinite(double, double)
+ * @see #isNaN(double, double)
+ * @see #abs(double, double)
+ * @see <a href="http://mathworld.wolfram.com/AbsoluteSquare.html">Absolute square</a>
+ */
+ public static double norm(double real, double imaginary) {
+ if (isInfinite(real, imaginary)) {
+ return Double.POSITIVE_INFINITY;
+ }
+ return real * real + imaginary * imaginary;
+ }
+
+ /**
+ * 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>
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @return {@code true} if the complex number contains NaN and no infinite parts.
+ * @see Double#isNaN(double)
+ * @see #isInfinite(double, double)
+ * @see #equals(Object) Complex.equals(Object)
+ */
+ public static boolean isNaN(double real, double imaginary) {
+ if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+ return !isInfinite(real, imaginary);
+ }
+ 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).
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @return {@code true} if the complex number contains an infinite value.
+ * @see Double#isInfinite(double)
+ */
+ public static boolean isInfinite(double real, double imaginary) {
+ return Double.isInfinite(real) || Double.isInfinite(imaginary);
+ }
+
+ /**
+ * Returns {@code true} if both real and imaginary component of the complex number are finite.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @return {@code true} if the complex number instance contains finite values.
+ * @see Double#isFinite(double)
+ */
+ public static boolean isFinite(double real, double imaginary) {
+ return Double.isFinite(real) && Double.isFinite(imaginary);
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/ComplexConjugate.html">conjugate</a>
+ * \( \overline{z} \) of the complex number \( z \).
+ *
+ * <p>\[ \begin{aligned}
+ * z &= a + i b \\
+ * \overline{z} &= a - i b \end{aligned}\]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the exponential of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ public static <R> R conj(double real, double imaginary, ComplexSink<R> action) {
+ return action.apply(real, -imaginary);
+ }
+
+ /**
+ * Returns the negation of both the real and imaginary parts of the complex number \( z \).
+ *
+ * <p>\[ \begin{aligned}
+ * z &= a + i b \\
+ * -z &= -a - i b \end{aligned} \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the exponential of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ public static <R> R negate(double real, double imaginary, ComplexSink<R> action) {
+ return action.apply(-real, -imaginary);
+ }
+
+ /**
+ * Returns the projection of the 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.
+ *
+ * <p>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>
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the exponential of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see #isInfinite(double, double)
+ * @see <a href="http://pubs.opengroup.org/onlinepubs/9699919799/functions/cproj.html">
+ * IEEE and ISO C standards: cproj</a>
+ */
+ public static <R> R proj(double real, double imaginary, ComplexSink<R> action) {
+ if (isInfinite(real, imaginary)) {
+ return action.apply(Double.POSITIVE_INFINITY, Math.copySign(0.0, imaginary));
+ }
+ return action.apply(real, imaginary);
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
+ * exponential function</a> of complex number using its real and
+ * imaginary part.
+ *
+ * <p>\[ \exp(z) = e^z \]
+ *
+ * <p>The exponential function of \( z \) is an entire function in the complex plane.
+ * 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)
+ * <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>Implements the formula:
+ *
+ * <p>\[ \exp(x + iy) = e^x (\cos(y) + i \sin(y)) \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \). part \( a \) of the complex number \( (a +ib) \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the exponential of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</a>
+ */
+ public static <R> R exp(double real, double imaginary, ComplexSink<R> action) {
+ 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 action.apply(0, Math.copySign(0, imaginary));
+ }
+ // (−∞ + iy) returns +0 cis(y), for finite y
+ zeroOrInf = 0;
+ } else {
+ // (+∞ + i0) returns +∞ + i0.
+ if (imaginary == 0) {
+ return action.apply(real, imaginary);
+ }
+ // (+∞ + 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 action.apply(real, Double.NaN);
+ }
+ // (+∞ + iy) returns (+∞ cis(y)), for finite nonzero y.
+ zeroOrInf = real;
+ }
+ return action.apply(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)
+ if (imaginary == 0) {
+ return action.apply(real, imaginary);
+ } else {
+ return action.apply(Double.NaN, Double.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 action.apply(Double.NaN, 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);
+ if (imaginary == 0) {
+ return action.apply(exp, imaginary);
+ }
+ return action.apply(exp * Math.cos(imaginary),
+ exp * Math.sin(imaginary));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html">
+ * natural logarithm</a> of the complex number using its real and imaginary parts.
+ *
+ * <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:
+ *
+ * <ul>
+ * <li>{@code log(conj(z)) == conj(log(z))}, where {@code conj} is the conjugate function.
+ * <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>Implements the formula:
+ *
+ * <p>\[ \ln(z) = \ln |z| + i \arg(z) \]
+ *
+ * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
+ *
+ * <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 real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the natural logarithm of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see Math#log(double)
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
+ */
+ public static <R> R log(double real, double imaginary, ComplexSink<R> action) {
+ return log(Math::log, HALF, LN_2, real, imaginary, action);
+ }
+
+ /**
+ * Returns the base 10
+ * <a href="http://mathworld.wolfram.com/CommonLogarithm.html">
+ * common logarithm</a> of the complex number using its real and imaginary parts.
+ *
+ * <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 log:
+ *
+ * <p>Implements the formula:
+ *
+ * <p>\[ \log_{10}(z) = \log_{10} |z| + i \arg(z) \]
+ *
+ * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the base 10 common logarithm of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see Math#log10(double)
+ */
+ public static <R> R log10(double real, double imaginary, ComplexSink<R> action) {
+ return log(Math::log10, LOG_10E_O_2, LOG10_2, real, imaginary, action);
+ }
+
+ /**
+ * Returns the logarithm of complex number using its real and
+ * imaginary parts and 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 real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the natural logarithm of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ private static <R> R log(DoubleUnaryOperator log,
+ double logOfeOver2,
+ double logOf2,
+ double real,
+ double imaginary,
+ ComplexSink<R> action) {
+
+ // Handle NaN
+ if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+ // Return NaN unless infinite
+ if (isInfinite(real, imaginary)) {
+ return action.apply(Double.POSITIVE_INFINITY, Double.NaN);
+ }
+ return action.apply(Double.NaN, Double.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 action.apply(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 initialize 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 action.apply(x, arg(real, imaginary));
+ }
+ // 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 action.apply(log.applyAsDouble(x), arg(real, imaginary));
+ }
+ // 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 action.apply(re, arg(real, imaginary));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/SquareRoot.html">
+ * square root</a> of the 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) \]
+ *
+ * <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 (1994)
+ * Implementing complex elementary functions using exception handling.
+ * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244.
+ * </blockquote>
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the square root of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sqrt/">Sqrt</a>
+ */
+ public static <R> R sqrt(double real, double imaginary, ComplexSink<R> action) {
+ return computeSqrt(real, imaginary, action);
+ }
+
+ /**
+ * Returns the square root of the complex number using its real
+ * and imaginary parts {@code sqrt(x + i y)}.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the square root of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ private static <R> R computeSqrt(double real, double imaginary, ComplexSink<R> action) {
+ // Handle NaN
+ if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+ // Check for infinite
+ if (Double.isInfinite(imaginary)) {
+ return action.apply(Double.POSITIVE_INFINITY, imaginary);
+ }
+ if (Double.isInfinite(real)) {
+ if (real == Double.NEGATIVE_INFINITY) {
+ return action.apply(Double.NaN, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
+ }
+ return action.apply(Double.POSITIVE_INFINITY, Double.NaN);
+ }
+ return action.apply(Double.NaN, Double.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;
+
+ // 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 action.apply(Double.POSITIVE_INFINITY, imaginary);
+ } else if (isPosInfinite(x)) {
+ if (real == Double.NEGATIVE_INFINITY) {
+ return action.apply(0, Math.copySign(Double.POSITIVE_INFINITY, imaginary));
+ }
+ return action.apply(Double.POSITIVE_INFINITY, Math.copySign(0, imaginary));
+ } else if (y == 0) {
+ // Real only
+ final double sqrtAbs = Math.sqrt(x);
+ if (real < 0) {
+ return action.apply(0, Math.copySign(sqrtAbs, imaginary));
+ }
+ return action.apply(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 action.apply(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 action.apply(t / 2, imaginary / t);
+ }
+ return action.apply(y / t, Math.copySign(t / 2, imaginary));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/Sine.html">
+ * sine</a> of the 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) \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the sine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sin/">Sin</a>
+ */
+ public static <R> R sin(double real, double imaginary, ComplexSink<R> action) {
+ // 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, (r, i) -> multiplyNegativeI(r, i, action));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/Cosine.html">
+ * cosine</a> of the 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) \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the cosine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cos/">Cos</a>
+ */
+ public static <R> R cos(double real, double imaginary, ComplexSink<R> action) {
+ // Define in terms of cosh
+ // cos(z) = cosh(iz)
+ // Multiply this number by I and compute cosh.
+ return cosh(-imaginary, real, action);
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/Tangent.html">
+ * tangent</a> of the 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) \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the tangent of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tan/">Tangent</a>
+ */
+ public static <R> R tan(double real, double imaginary, ComplexSink<R> action) {
+ // 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, (r, i) -> multiplyNegativeI(r, i, action));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/InverseSine.html">
+ * inverse sine</a> of the 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>\[ \begin{aligned}
+ * \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] \end{aligned} \]
+ *
+ * <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>}.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse sine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSin/">ArcSin</a>
+ */
+ public static <R> R asin(final double real, final double imaginary,
+ ComplexSink<R> action) {
+ return computeAsin(real, imaginary, action);
+ }
+
+ /**
+ * 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>}. This method only (and not
+ * invoked methods within) is distributed under the Boost Software License V1.0.
+ * The original notice is shown below and the licence is shown in full in LICENSE:
+ * <pre>
+ * (C) Copyright John Maddock 2005.
+ * Distributed under the Boost Software License, Version 1.0. (See accompanying
+ * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
+ * </pre>
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse sine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ private static <R> R computeAsin(final double real, final double imaginary,
+ ComplexSink<R> action) {
+ // 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 {
+ return action.apply(Double.NaN, Double.NaN);
+ }
+ } else if (Double.isNaN(y)) {
+ if (x == 0) {
+ re = 0;
+ im = y;
+ } else if (isPosInfinite(x)) {
+ re = y;
+ im = x;
+ } else {
+ return action.apply(Double.NaN, Double.NaN);
+ }
+ } else if (isPosInfinite(x)) {
+ re = isPosInfinite(y) ? PI_OVER_4 : PI_OVER_2;
+ im = x;
+ } else if (isPosInfinite(y)) {
+ re = 0;
+ im = y;
+ } else {
+ // Special case for real numbers:
+ if (y == 0 && x <= 1) {
+ return action.apply(Math.asin(real), imaginary);
+ }
+
+ final double xp1 = x + 1;
+ final double xm1 = x - 1;
+
+ 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;
+
+ if (b <= B_CROSSOVER) {
+ re = Math.asin(b);
+ } else {
+ final double apx = a + x;
+ if (x <= 1) {
+ re = Math.atan(x / Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))));
+ } else {
+ re = Math.atan(x / (y * Math.sqrt(0.5 * (apx / (r + xp1) + apx / (s + xm1)))));
+ }
+ }
+
+ 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 4
+ if (y <= (EPSILON * Math.abs(xm1))) {
+ if (x < 1) {
+ re = Math.asin(x);
+ im = y / Math.sqrt(xp1 * (1 - x));
+ } else {
+ re = PI_OVER_2;
+ if ((Double.MAX_VALUE / xp1) > xm1) {
+ // xp1 * xm1 won't overflow:
+ im = Math.log1p(xm1 + Math.sqrt(xp1 * xm1));
+ } else {
+ 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 = PI_OVER_2 - Math.sqrt(y);
+ im = Math.sqrt(y);
+ } else if (EPSILON * y - 1 >= x) {
+ // Possible underflow:
+ re = x / y;
+ im = LN_2 + Math.log(y);
+ } else if (x > 1) {
+ re = Math.atan(x / y);
+ final double xoy = x / y;
+ im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
+ } else {
+ final double a = Math.sqrt(1 + y * y);
+ // Possible underflow:
+ re = x / a;
+ im = 0.5 * Math.log1p(2 * y * (y + a));
+ }
+ }
+ }
+ return action.apply(changeSign(re, real),
+ changeSign(im, imaginary));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/InverseCosine.html">
+ * inverse cosine</a> of the 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:
+ *
+ * <p>\[ \begin{aligned}
+ * \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] \end{aligned} \]
+ *
+ * <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/acos.hpp>}.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse cosine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCos/">ArcCos</a>
+ */
+ public static <R> R acos(final double real, final double imaginary,
+ final ComplexSink<R> action) {
+ return computeAcos(real, imaginary, action);
+ }
+
+ /**
+ * Returns the inverse cosine of the complex number.
+ *
+ * <p>This function exists to allow implementation of the identity
+ * {@code acosh(z) = +-i acos(z)}.
+ *
+ * <p>Adapted from {@code <boost/math/complex/acos.hpp>}. This method only (and not
+ * invoked methods within) is distributed under the Boost Software License V1.0.
+ * The original notice is shown below and the licence is shown in full in LICENSE:
+ * <pre>
+ * (C) Copyright John Maddock 2005.
+ * Distributed under the Boost Software License, Version 1.0. (See accompanying
+ * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
+ * </pre>
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse cosine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ private static <R> R computeAcos(final double real, final double imaginary,
+ final ComplexSink<R> action) {
+ // 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 action.apply(imaginary, real);
+ } else {
+ re = 0;
+ im = Double.POSITIVE_INFINITY;
+ }
+ } else if (Double.isNaN(x)) {
+ if (isPosInfinite(y)) {
+ return action.apply(x, -imaginary);
+ }
+ return action.apply(Double.NaN, Double.NaN);
+ } else if (isPosInfinite(y)) {
+ re = PI_OVER_2;
+ im = y;
+ } else if (Double.isNaN(y)) {
+ return action.apply(x == 0 ? PI_OVER_2 : y, y);
+ } else {
+ // Special case for real numbers:
+ if (y == 0 && x <= 1) {
+ return action.apply(x == 0 ? PI_OVER_2 : Math.acos(real), -imaginary);
+ }
+
+ final double xp1 = x + 1;
+ final double xm1 = x - 1;
+
+ 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;
+
+ 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);
+ }
+ }
+
+ 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));
+ }
+ }
+ }
+ return action.apply(negative(real) ? Math.PI - re : re,
+ negative(imaginary) ? im : -im);
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/InverseTangent.html">
+ * inverse tangent</a> of the 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.
+ *
+ * <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>As per the C99 standard this function is computed using the trigonomic identity:
+ * \[ \tan^{-1}(z) = -i \tanh^{-1}(iz) \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse tangent of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTan/">ArcTan</a>
+ */
+ public static <R> R atan(double real, double imaginary, ComplexSink<R> action) {
+ // 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, (r, i) -> multiplyNegativeI(r, i, action));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
+ * hyperbolic sine</a> of the complexnumber.
+ *
+ * <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:
+ *
+ * <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 ofCis(double) in Complex class).
+ * <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>\[ \sinh(x + iy) = \sinh(x)\cos(y) + i \cosh(x)\sin(y) \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the hyperbolic sine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action..
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sinh/">Sinh</a>
+ */
+ public static <R> R sinh(double real, double imaginary, ComplexSink<R> action) {
+ return computeSinh(real, imaginary, action);
+ }
+
+ /**
+ * Returns the hyperbolic sine of the complex number using its real
+ * and imaginary parts.
+ *
+ * <p>This function exists to allow implementation of the identity
+ * {@code sin(z) = -i sinh(iz)}.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the hyperbolic sine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ public static <R> R computeSinh(double real, double imaginary, ComplexSink<R> action) {
+ if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
+ return action.apply(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 action.apply(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 action.apply(real, Double.NaN);
+ }
+ if (imaginary == 0) {
+ // Real-only sinh(x).
+ return action.apply(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, action);
+ }
+ // No overflow of sinh/cosh
+ return action.apply(Math.sinh(real) * Math.cos(imaginary),
+ Math.cosh(real) * Math.sin(imaginary));
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
+ * hyperbolic cosine</a> of the complexnumber.
+ *
+ * <p>\[ \cosh(z) = \frac{1}{2} \left( e^{z} + e^{-z} \right) \]
+ *
+ * <p>The hyperbolic cosine 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().cosh() == z.cosh().conj()}.
+ * <li>This is an even function: \( \cosh(z) = \cosh(-z) \).
+ * <li>If {@code z} is +0 + i0, returns 1 + i0.
+ * <li>If {@code z} is +0 + i∞, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified; "invalid" floating-point operation).
+ * <li>If {@code z} is +0 + iNaN, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified).
+ * <li>If {@code z} is x + i∞ for finite nonzero 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 finite nonzero y, returns +∞ cis(y) (see ofCis(double) in Complex class).
+ * <li>If {@code z} is +∞ + i∞, returns ±∞ + iNaN (where the sign of the real part of the result is unspecified).
+ * <li>If {@code z} is +∞ + iNaN, returns +∞ + iNaN.
+ * <li>If {@code z} is NaN + i0, returns NaN ± i0 (where the sign of the imaginary part of the result is unspecified).
+ * <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>\[ \cosh(x + iy) = \cosh(x)\cos(y) + i \sinh(x)\sin(y) \]
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the hyperbolic tangent of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cosh/">Cosh</a>
+ */
+ public static <R> R cosh(double real, double imaginary, ComplexSink<R> action) {
+ return computeCosh(real, imaginary, action);
+ }
+
+ /**
+ * Returns the hyperbolic cosine of the complex number using its real and
+ * imaginary parts.
+ *
+ * <p>This function exists to allow implementation of the identity
+ * {@code cos(z) = cosh(iz)}.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the hyperbolic cosine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ public static <R> R computeCosh(double real, double imaginary, ComplexSink<R> action) {
+ // ISO C99: Preserve the even function by mapping to positive
+ // f(z) = f(-z)
+ if (Double.isInfinite(real) && !Double.isFinite(imaginary)) {
+ return action.apply(Math.abs(real), Double.NaN);
+ }
+ if (real == 0) {
+ // Imaginary-only cosh(iy) = cos(y).
+ if (Double.isFinite(imaginary)) {
+ // Maintain periodic property with respect to the imaginary component.
+ // sinh(+/-0.0) * sin(+/-x) = +/-0 * sin(x)
+ return action.apply(Math.cos(imaginary),
+ changeSign(real, Math.sin(imaginary)));
+ }
+ // If imaginary is inf/NaN the sign of the imaginary part is unspecified.
+ // Although not required by C99 changing the sign maintains the conjugate equality.
+ // It is not possible to also maintain the even function (hence the unspecified sign).
+ return action.apply(Double.NaN, changeSign(real, imaginary));
+ }
+ if (imaginary == 0) {
+ // Real-only cosh(x).
+ // Change sign to preserve conjugate equality and even function.
+ // sin(+/-0) * sinh(+/-x) = +/-0 * +/-a (sinh is monotonic and same sign)
+ // => change the sign of imaginary using real. Handles special case of infinite real.
+ // If real is NaN the sign of the imaginary part is unspecified.
+ return action.apply(Math.cosh(real), changeSign(imaginary, real));
+ }
+ final double x = Math.abs(real);
+ if (x > SAFE_EXP) {
+ // Approximate sinh/cosh(x) using exp^|x| / 2
+ return coshsinh(x, real, imaginary, false, action);
+ }
+ // No overflow of sinh/cosh
+ return action.apply(Math.cosh(real) * Math.cos(imaginary),
+ Math.sinh(real) * Math.sin(imaginary));
+ }
+
+ /**
+ * Compute cosh or sinh when the absolute real component |x| is large. In this case
+ * cosh(x) and sinh(x) can be approximated by exp(|x|) / 2:
+ *
+ * <pre>
+ * cosh(x+iy) real = (e^|x| / 2) * cos(y)
+ * cosh(x+iy) imag = (e^|x| / 2) * sin(y) * sign(x)
+ * sinh(x+iy) real = (e^|x| / 2) * cos(y) * sign(x)
+ * sinh(x+iy) imag = (e^|x| / 2) * sin(y)
+ * </pre>
+ *
+ * @param x Absolute real component |x|.
+ * @param real Real part (x).
+ * @param imaginary Imaginary part (y).
+ * @param sinh Set to true to compute sinh, otherwise cosh.
+ * @param action Consumer for the cosh or sinh of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ private static <R> R coshsinh(double x, double real, double imaginary, boolean sinh,
+ ComplexSink<R> action) {
+ // Always require the cos and sin.
+ double re = Math.cos(imaginary);
+ double im = Math.sin(imaginary);
+ // Compute the correct function
+ if (sinh) {
+ re = changeSign(re, real);
+ } else {
+ im = changeSign(im, real);
+ }
+ // Multiply by (e^|x| / 2).
+ // Overflow safe computation since sin/cos can be very small allowing a result
+ // when e^x overflows: e^x / 2 = (e^m / 2) * e^m * e^(x-2m)
+ if (x > SAFE_EXP * 3) {
+ // e^x > e^m * e^m * e^m
+ // y * (e^m / 2) * e^m * e^m will overflow when starting with Double.MIN_VALUE.
+ // Note: Do not multiply by +inf to safeguard against sin(y)=0.0 which
+ // will create 0 * inf = nan.
+ re *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
+ im *= Double.MAX_VALUE * Double.MAX_VALUE * Double.MAX_VALUE;
+ } else {
+ // Initial part of (e^x / 2) using (e^m / 2)
+ re *= EXP_M / 2;
+ im *= EXP_M / 2;
+ double xm;
+ if (x > SAFE_EXP * 2) {
+ // e^x = e^m * e^m * e^(x-2m)
+ re *= EXP_M;
+ im *= EXP_M;
+ xm = x - SAFE_EXP * 2;
+ } else {
+ // e^x = e^m * e^(x-m)
+ xm = x - SAFE_EXP;
+ }
+ final double exp = Math.exp(xm);
+ re *= exp;
+ im *= exp;
+ }
+ return action.apply(re, im);
+ }
+
+ /**
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/HyperbolicTangent.html">
+ * hyperbolic tangent</a> of the complexnumber.
+ *
+ * <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}.
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
* @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
- * @return The absolute value.
+ * @param action Consumer for the hyperbolic tangent of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
*/
- public static double abs(double real, double imaginary) {
- // Specialised implementation of hypot.
- // See NUMBERS-143
- return hypot(real, imaginary);
+ public static <R> R tanh(double real, double imaginary, ComplexSink<R> action) {
+ return computeTanh(real, imaginary, action);
}
/**
- * Returns the argument of the 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.
+ * Returns the hyperbolic tangent of the complex number using its real
+ * and imaginary 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 \).
+ * <p>This function exists to allow implementation of the identity
+ * {@code tan(z) = -i tanh(iz)}.
*
- * @param r Real part \( a \) of the complex number \( (a +ib) \).
- * @param i Imaginary part \( b \) of the complex number \( (a +ib) \).
- * @return The argument of the complex number.
- * @see Math#atan2(double, double)
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the hyperbolic tangent of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
*/
- public static double arg(double r, double i) {
- // Delegate
- return Math.atan2(i, r);
+ public static <R> R computeTanh(double real, double imaginary, ComplexSink<R> action) {
+ // 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 action.apply(Math.copySign(1, real),
+ Math.copySign(0, sign));
+ }
+ // imaginary is infinite or NaN
+ return action.apply(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 action.apply(real == 0 ? real : Double.NaN,
+ imaginary == 0 ? imaginary : Double.NaN);
+ }
+
+ // 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 action.apply(real, Math.tan(imaginary));
+ }
+ if (imaginary == 0) {
+ // Identity: sinh 2x / (1 + cosh 2x) = tanh(x)
+ return action.apply(Math.tanh(real), 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 action.apply(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 action.apply(sinhx * coshx / divisor,
+ siny * cosy / divisor);
}
/**
- * Returns {@code true} if either real or imaginary component of the complex number is infinite.
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/InverseHyperbolicSine.html">
+ * inverse hyperbolic sine</a> of the complex number.
*
- * <p>Note: A complex number with at least one infinite part is regarded
- * as an infinity (even if its other part is a NaN).
- * @param real Real part \( a \) of the complex number \( (a +ib) \).
- * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
- * @return {@code true} if the complex number contains an infinite value.
- * @see Double#isInfinite(double)
+ * <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().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>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
+ * 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) \]
+ *
+ * @param real part of Complex number
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
+ * @param action Consumer for the inverse hyperbolic sine of the complex number.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSinh/">ArcSinh</a>
*/
- public static boolean isInfinite(double real, double imaginary) {
- return Double.isInfinite(real) || Double.isInfinite(imaginary);
+ public static <R> R asinh(double real, double imaginary, ComplexSink<R> action) {
+ // 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, (r, i) -> multiplyNegativeI(r, i, action));
}
/**
* Returns the
- * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html">
- * natural logarithm</a> of the complex number using its real and imaginary parts.
+ * <a href="http://mathworld.wolfram.com/InverseHyperbolicCosine.html">
+ * inverse hyperbolic cosine</a> of complex number using its real and imaginary parts.
*
- * <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>\[ \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 log(conj(z)) == conj(log(z))}, where {@code conj} is the conjugate function.
- * <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>{@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 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 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.
@@ -195,163 +1840,271 @@ public final class ComplexFunctions {
* <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>\[ \ln(z) = \ln |z| + i \arg(z) \]
+ * <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>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>\[ \cosh^{-1}(z) = \pm i \cos^{-1}(z) \]
+ *
+ * <p>The sign of the multiplier is chosen to give {@code z.acosh().real() >= 0}
+ * and compatibility with the C99 standard.
*
* @param real Real part \( a \) of the complex number \( (a +ib \).
- * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
- * @param action Consumer for the natural logarithm of the complex number.
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse hyperbolic cosine of the complex number.
* @param <R> the return type of the supplied action.
* @return the object returned by the supplied action.
- * @see Math#log(double)
- * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCosh/">ArcCosh</a>
*/
- public static <R> R log(double real, double imaginary, ComplexSink<R> action) {
- return log(Math::log, HALF, LN_2, real, imaginary, action);
+ public static <R> R acosh(double real, double imaginary, ComplexSink<R> action) {
+ // 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 action.apply(Double.NaN, PI_OVER_2);
+ }
+ return acos(real, imaginary, (re, im) ->
+ // Set the sign appropriately for real >= 0
+ (negative(im)) ?
+ // Multiply by I
+ action.apply(-im, re) :
+ // Multiply by -I
+ action.apply(im, -re)
+ );
}
/**
- * Returns the base 10
- * <a href="http://mathworld.wolfram.com/CommonLogarithm.html">
- * common logarithm</a> of the complex number using its real and imaginary parts.
+ * Returns the
+ * <a href="http://mathworld.wolfram.com/InverseHyperbolicTangent.html">
+ * inverse hyperbolic tangent</a> of the 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 log:
+ * <p>\[ \tanh^{-1}(z) = \frac{1}{2} \ln \left( \frac{1 + z}{1 - z} \right) \]
*
- * <p>Implements the formula:
+ * <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>\[ \log_{10}(z) = \log_{10} |z| + i \arg(z) \]
+ * <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>
*
- * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument.
+ * <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 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>}.
*
* @param real Real part \( a \) of the complex number \( (a +ib \).
- * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
- * @param action Consumer for the base 10 common logarithm of the complex number.
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse hyperbolic tangent of the complex number.
* @param <R> the return type of the supplied action.
* @return the object returned by the supplied action.
- * @see Math#log10(double)
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTanh/">ArcTanh</a>
*/
- public static <R> R log10(double real, double imaginary, ComplexSink<R> action) {
- return log(Math::log10, LOG_10E_O_2, LOG10_2, real, imaginary, action);
+ public static <R> R atanh(final double real, final double imaginary,
+ final ComplexSink<R> action) {
+ return computeAtanh(real, imaginary, action);
}
/**
- * Returns the logarithm of complex number using its real and
- * imaginary parts and the provided function.
- * Implements the formula:
+ * Returns the inverse hyperbolic tangent of the complex number using its
+ * real and imaginary parts.
*
- * <pre>
- * log(x + i y) = log(|x + i y|) + i arg(x + i y)</pre>
+ * <p>This function exists to allow implementation of the identity
+ * {@code atan(z) = -i atanh(iz)}.
*
- * <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.
+ * <p>Adapted from {@code <boost/math/complex/atanh.hpp>}. This method only (and not
+ * invoked methods within) is distributed under the Boost Software License V1.0.
+ * The original notice is shown below and the licence is shown in full in LICENSE:
+ * <pre>
+ * (C) Copyright John Maddock 2005.
+ * Distributed under the Boost Software License, Version 1.0. (See accompanying
+ * file LICENSE or copy at https://www.boost.org/LICENSE_1_0.txt)
+ * </pre>
*
- * @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 real Real part \( a \) of the complex number \( (a +ib \).
- * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib \).
- * @param action Consumer for the natural logarithm of the complex number.
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the inverse hyperbolic tangent of the complex number.
* @param <R> the return type of the supplied action.
* @return the object returned by the supplied action.
*/
- private static <R> R log(DoubleUnaryOperator log,
- double logOfeOver2,
- double logOf2,
- double real,
- double imaginary,
- ComplexSink<R> action) {
-
- // Handle NaN
- if (Double.isNaN(real) || Double.isNaN(imaginary)) {
- // Return NaN unless infinite
- if (isInfinite(real, imaginary)) {
- return action.apply(Double.POSITIVE_INFINITY, Double.NaN);
- }
- return action.apply(Double.NaN, Double.NaN);
- }
-
- // Returns the real part:
- // log(sqrt(x^2 + y^2))
- // log(x^2 + y^2) / 2
-
- // Compute with positive values
+ private static <R> R computeAtanh(final double real, final double imaginary,
+ final ComplexSink<R> action) {
+ // Compute with positive values and determine sign at the end
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 action.apply(Double.NEGATIVE_INFINITY,
- negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary);
- }
-
+ // The result (without sign correction)
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;
+ 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 action.apply(0, Math.copySign(PI_OVER_2, imaginary));
+ }
+ return action.apply(Double.NaN, Double.NaN);
+ } else if (Double.isNaN(y)) {
+ if (isPosInfinite(x)) {
+ return action.apply(Math.copySign(0, real), Double.NaN);
+ }
+ if (x == 0) {
+ return action.apply(real, Double.NaN);
+ }
+ return action.apply(Double.NaN, Double.NaN);
} else {
- // Check for over/underflow in |z|
- // When scaling:
- // log(a / b) = log(a) - log(b)
- // So initialize 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 action.apply(x, arg(real, imaginary));
+ // 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
+ if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
+ // Normal computation within a safe region.
+ // 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;
}
- // 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 action.apply(log.applyAsDouble(x), arg(real, 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)
+ // Modified from boost to use the custom high precision method.
+ 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.
+ // C99. G.7: Special case for imaginary only numbers
+ if (x == 0) {
+ if (imaginary == 0) {
+ return action.apply(real, imaginary);
+ }
+ // atanh(iy) = i atan(y)
+ return action.apply(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 {
+ // 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));
}
- // 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 action.apply(re, arg(real, imaginary));
+ re /= 4;
+ im /= 2;
+ return action.apply(changeSign(re, real),
+ changeSign(im, imaginary));
}
/**
@@ -595,8 +2348,7 @@ public final class ComplexFunctions {
* @param y the y value
* @return {@code x^2 + y^2 - 1}.
*/
- //TODO - make it private in future
- static double x2y2m1(double x, double y) {
+ private static double x2y2m1(double x, double y) {
// Hull et al used (x-1)*(x+1)+y*y.
// From the paper on page 236:
@@ -840,8 +2592,71 @@ public final class ComplexFunctions {
* @param d Value.
* @return {@code true} if {@code d} is +inf.
*/
- //TODO - make private in future
- static boolean isPosInfinite(double d) {
+ private static boolean isPosInfinite(double d) {
return d == Double.POSITIVE_INFINITY;
}
+
+ /**
+ * Check that an absolute value is finite. Used to replace {@link Double#isFinite(double)}
+ * when the input value is known to be positive (i.e. in the case where it has been
+ * set using {@link Math#abs(double)}).
+ *
+ * @param d Value.
+ * @return {@code true} if {@code d} is +finite.
+ */
+ private static boolean isPosFinite(double d) {
+ return d <= Double.MAX_VALUE;
+ }
+
+ /**
+ * Create a complex number given the real and imaginary parts, then multiply by {@code -i}.
+ * This is used in functions that implement trigonomic identities. It is the functional
+ * equivalent of:
+ *
+ * <pre>
+ * z = new Complex(real, imaginary).multiplyImaginary(-1);</pre>
+ *
+ * @param real Real part \( a \) of the complex number \( (a +ib \).
+ * @param imaginary Imaginary part \( b \) of the complex number \( (a +ib) \).
+ * @param action Consumer for the complex number multiplied by {@code -i}.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ */
+ private static <R> R multiplyNegativeI(double real, double imaginary, ComplexSink<R> action) {
+ return action.apply(imaginary, -real);
+ }
+
+ /**
+ * Change the sign of the magnitude based on the signed value.
+ *
+ * <p>If the signed value is negative then the result is {@code -magnitude}; otherwise
+ * return {@code magnitude}.
+ *
+ * <p>A signed value of {@code -0.0} is treated as negative. A signed value of {@code NaN}
+ * is treated as positive.
+ *
+ * <p>This is not the same as {@link Math#copySign(double, double)} as this method
+ * will change the sign based on the signed value rather than copy the sign.
+ *
+ * @param magnitude the magnitude
+ * @param signedValue the signed value
+ * @return magnitude or -magnitude.
+ * @see #negative(double)
+ */
+ private static double changeSign(double magnitude, double signedValue) {
+ return negative(signedValue) ? -magnitude : magnitude;
+ }
+
+ /**
+ * Checks if both x and y are in the region defined by the minimum and maximum.
+ *
+ * @param x x value.
+ * @param y y value.
+ * @param min the minimum (exclusive).
+ * @param max the maximum (exclusive).
+ * @return true if inside the region.
+ */
+ private static boolean inRegion(double x, double y, double min, double max) {
+ return x < max && x > min && y < max && y > min;
+ }
}
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java
index 1a139053..62dc190e 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java
@@ -223,19 +223,21 @@ class CReferenceTest {
/**
* Assert the operation using the data loaded from test resources.
*
- * @param name the operation name
- * @param operation the operation
- * @param maxUlps the maximum units of least precision between the two values
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex object.
+ * @param operation2 Operation on the complex real and imaginary parts.
+ * @param maxUlps Maximum units of least precision between the two values.
*/
private static void assertOperation(String name,
- UnaryOperator<Complex> operation, long maxUlps) {
+ UnaryOperator<Complex> operation1,
+ ComplexUnaryOperator<ComplexNumber> operation2,
+ long maxUlps) {
final List<Complex[]> data = loadTestData(name);
final long ulps = getTestUlps(maxUlps);
for (final Complex[] pair : data) {
- assertComplex(pair[0], name, operation, pair[1], ulps);
+ assertComplex(pair[0], name, operation1, operation2, pair[1], ulps);
}
}
-
/**
* Assert the operation using the data loaded from test resources.
*
@@ -252,25 +254,6 @@ class CReferenceTest {
}
}
- /**
- * Assert the operation using the data loaded from test resources.
- *
- * @param name Operation name.
- * @param operation1 Operation on the Complex object.
- * @param operation2 Operation on the complex real and imaginary parts.
- * @param maxUlps Maximum units of least precision between the two values.
- */
- private static void assertOperation(String name,
- UnaryOperator<Complex> operation1,
- ComplexUnaryOperator<ComplexNumber> operation2,
- long maxUlps) {
- final List<Complex[]> data = loadTestData(name);
- final long ulps = getTestUlps(maxUlps);
- for (final Complex[] pair : data) {
- assertComplex(pair[0], name, operation1, operation2, pair[1], ulps);
- }
- }
-
/**
* Assert the operation using the data loaded from test resources.
*
@@ -314,47 +297,47 @@ class CReferenceTest {
@Test
void testAcos() {
- assertOperation("acos", Complex::acos, 2);
+ assertOperation("acos", Complex::acos, ComplexFunctions::acos, 2);
}
@Test
void testAcosh() {
- assertOperation("acosh", Complex::acosh, 2);
+ assertOperation("acosh", Complex::acosh, ComplexFunctions::acosh, 2);
}
@Test
void testAsinh() {
// Odd function: negative real cases defined by positive real cases
- assertOperation("asinh", Complex::asinh, 3);
+ assertOperation("asinh", Complex::asinh, ComplexFunctions::asinh, 3);
}
@Test
void testAtanh() {
// Odd function: negative real cases defined by positive real cases
- assertOperation("atanh", Complex::atanh, 1);
+ assertOperation("atanh", Complex::atanh, ComplexFunctions::atanh, 1);
}
@Test
void testCosh() {
// Even function: negative real cases defined by positive real cases
- assertOperation("cosh", Complex::cosh, 2);
+ assertOperation("cosh", Complex::cosh, ComplexFunctions::cosh, 2);
}
@Test
void testSinh() {
// Odd function: negative real cases defined by positive real cases
- assertOperation("sinh", Complex::sinh, 2);
+ assertOperation("sinh", Complex::sinh, ComplexFunctions::sinh, 2);
}
@Test
void testTanh() {
// Odd function: negative real cases defined by positive real cases
- assertOperation("tanh", Complex::tanh, 2);
+ assertOperation("tanh", Complex::tanh, ComplexFunctions::tanh, 2);
}
@Test
void testExp() {
- assertOperation("exp", Complex::exp, 2);
+ assertOperation("exp", Complex::exp, ComplexFunctions::exp, 2);
}
@Test
@@ -364,7 +347,7 @@ class CReferenceTest {
@Test
void testSqrt() {
- assertOperation("sqrt", Complex::sqrt, 1);
+ assertOperation("sqrt", Complex::sqrt, ComplexFunctions::sqrt, 1);
}
@Test
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java
index 5d02360e..2ac487b0 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java
@@ -218,88 +218,6 @@ class CStandardTest {
() -> String.format("%s expected: %s %s %s = %s", expectedName, c1, operationName, c2, z));
}
- /**
- * Assert the operation on the complex number satisfies the conjugate equality.
- *
- * <pre>
- * op(conj(z)) = conj(op(z))
- * </pre>
- *
- * <p>The results must be binary equal. This includes the sign of zero.
- *
- * <h2>ISO C99 equalities</h2>
- *
- * <p>Note that this method currently enforces the conjugate equalities for some cases
- * where the sign of the real/imaginary parts are unspecified in ISO C99. This is
- * allowed (since they are unspecified). The sign specification is appropriately
- * handled during testing of odd/even functions. There are some functions where it
- * is not possible to satisfy the conjugate equality and also the odd/even rule.
- * The compromise made here is to satisfy only one and the other is allowed to fail
- * only on the sign of the output result. Known functions where this applies are:
- *
- * <ul>
- * <li>asinh(NaN, inf)
- * <li>atanh(NaN, inf)
- * <li>cosh(NaN, 0.0)
- * <li>sinh(inf, inf)
- * <li>sinh(inf, nan)
- * </ul>
- *
- * @param operation the operation
- */
- private static void assertConjugateEquality(UnaryOperator<Complex> operation) {
- // Edge cases. Inf/NaN are specifically handled in the C99 test cases
- // but are repeated here to enforce the conjugate equality even when the C99
- // standard does not specify a sign. This may be revised in the future.
- final double[] parts = {Double.NEGATIVE_INFINITY, -1, -0.0, 0.0, 1,
- Double.POSITIVE_INFINITY, Double.NaN};
- for (final double x : parts) {
- for (final double y : parts) {
- // No conjugate for imaginary NaN
- if (!Double.isNaN(y)) {
- assertConjugateEquality(complex(x, y), operation, UnspecifiedSign.NONE);
- }
- }
- }
- // Random numbers
- final UniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create();
- for (int i = 0; i < 100; i++) {
- final double x = next(rng);
- final double y = next(rng);
- assertConjugateEquality(complex(x, y), operation, UnspecifiedSign.NONE);
- }
- }
-
- /**
- * Assert the operation on the complex number satisfies the conjugate equality.
- *
- * <pre>
- * op(conj(z)) = conj(op(z))
- * </pre>
- *
- * <p>The results must be binary equal; the sign of the complex number is first processed
- * using the provided sign specification.
- *
- * @param z the complex number
- * @param operation the operation
- * @param sign the sign specification
- */
- private static void assertConjugateEquality(Complex z,
- UnaryOperator<Complex> operation, UnspecifiedSign sign) {
- final Complex c1 = operation.apply(z.conj());
- final Complex c2 = operation.apply(z).conj();
- final Complex t1 = sign.removeSign(c1);
- final Complex t2 = sign.removeSign(c2);
-
- // Test for binary equality
- if (!equals(t1.getReal(), t2.getReal()) ||
- !equals(t1.getImaginary(), t2.getImaginary())) {
- Assertions.fail(
- String.format("Conjugate equality failed (z=%s). Expected: %s but was: %s (Unspecified sign = %s)",
- z, c1, c2, sign));
- }
- }
-
/**
* Assert the operation on the complex number satisfies the conjugate equality.
*
@@ -383,9 +301,10 @@ class CStandardTest {
ComplexUnaryOperator<ComplexNumber> operation2,
UnspecifiedSign sign) {
- final Complex zConj = z.conj();
+ final Complex zConj = TestUtils.assertSame(z, "conj", Complex::conj, ComplexFunctions::conj);
final Complex c1 = TestUtils.assertSame(zConj, name, operation1, operation2);
- final Complex c2 = TestUtils.assertSame(z, name, operation1, operation2).conj();
+ final Complex result = TestUtils.assertSame(z, name, operation1, operation2);
+ final Complex c2 = TestUtils.assertSame(result, "conj", Complex::conj, ComplexFunctions::conj);
final Complex t1 = sign.removeSign(c1);
final Complex t2 = sign.removeSign(c2);
@@ -409,10 +328,18 @@ class CStandardTest {
*
* <p>The results must be binary equal. This includes the sign of zero.
*
- * @param operation the operation
+ * <p>Assert the operation on the complex number is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts.
+ *
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex object.
+ * @param operation2 Operation on the complex real and imaginary parts.
* @param type the type
*/
- private static void assertFunctionType(UnaryOperator<Complex> operation, FunctionType type) {
+ private static void assertFunctionType(String name,
+ UnaryOperator<Complex> operation1,
+ ComplexUnaryOperator<ComplexNumber> operation2,
+ FunctionType type) {
// Note: It may not be possible to satisfy the conjugate equality
// and be an odd/even function with regard to zero.
// The C99 standard allows for these cases to have unspecified sign.
@@ -426,7 +353,7 @@ class CStandardTest {
final double[] parts = {-2, -1, -0.0, 0.0, 1, 2};
for (final double x : parts) {
for (final double y : parts) {
- assertFunctionType(complex(x, y), operation, type, UnspecifiedSign.NONE);
+ assertFunctionType(complex(x, y), name, operation1, operation2, type, UnspecifiedSign.NONE);
}
}
// Random numbers
@@ -434,7 +361,7 @@ class CStandardTest {
for (int i = 0; i < 100; i++) {
final double x = next(rng);
final double y = next(rng);
- assertFunctionType(complex(x, y), operation, type, UnspecifiedSign.NONE);
+ assertFunctionType(complex(x, y), name, operation1, operation2, type, UnspecifiedSign.NONE);
}
}
@@ -449,15 +376,23 @@ class CStandardTest {
* <p>The results must be binary equal; the sign of the complex number is first processed
* using the provided sign specification.
*
+ * <p>Assert the operation on the complex number is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts.
+ *
* @param z the complex number
- * @param operation the operation
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex object.
+ * @param operation2 Operation on the complex real and imaginary parts.
* @param type the type (assumed to be ODD/EVEN)
* @param sign the sign specification
*/
private static void assertFunctionType(Complex z,
- UnaryOperator<Complex> operation, FunctionType type, UnspecifiedSign sign) {
- final Complex c1 = operation.apply(z);
- Complex c2 = operation.apply(z.negate());
+ String name,
+ UnaryOperator<Complex> operation1,
+ ComplexUnaryOperator<ComplexNumber> operation2,
+ FunctionType type, UnspecifiedSign sign) {
+ final Complex c1 = TestUtils.assertSame(z, name, operation1, operation2);
+ Complex c2 = TestUtils.assertSame(z.negate(), name, operation1, operation2);
if (type == FunctionType.ODD) {
c2 = c2.negate();
}
@@ -484,127 +419,20 @@ class CStandardTest {
*
* <p>The results must be binary equal. This includes the sign of zero.
*
- * @param z the complex
- * @param operation the operation
- * @param expected the expected
- */
- private static void assertComplex(Complex z,
- UnaryOperator<Complex> operation, Complex expected) {
- assertComplex(z, operation, expected, FunctionType.NONE, UnspecifiedSign.NONE);
- }
-
- /**
- * Assert the operation on the complex number is equal to the expected value. If the
- * imaginary part is not NaN the operation must also satisfy the conjugate equality.
- *
- * <pre>
- * op(conj(z)) = conj(op(z))
- * </pre>
- *
- * <p>The results must be binary equal. This includes the sign of zero.
- *
- * @param z the complex
- * @param operation the operation
- * @param expected the expected
- * @param sign the sign specification
- */
- private static void assertComplex(Complex z,
- UnaryOperator<Complex> operation, Complex expected, UnspecifiedSign sign) {
- assertComplex(z, operation, expected, FunctionType.NONE, sign);
- }
-
- /**
- * Assert the operation on the complex number is equal to the expected value. If the
- * imaginary part is not NaN the operation must also satisfy the conjugate equality.
- *
- * <pre>
- * op(conj(z)) = conj(op(z))
- * </pre>
- *
- * <p>If the function type is ODD/EVEN the operation must satisfy the function type equality.
- *
- * <pre>
- * Odd : f(z) = -f(-z)
- * Even: f(z) = f(-z)
- * </pre>
- *
- * <p>The results must be binary equal. This includes the sign of zero.
- *
- * @param z the complex
- * @param operation the operation
- * @param expected the expected
- * @param type the type
- */
- private static void assertComplex(Complex z,
- UnaryOperator<Complex> operation, Complex expected,
- FunctionType type) {
- assertComplex(z, operation, expected, type, UnspecifiedSign.NONE);
- }
-
- /**
- * Assert the operation on the complex number is equal to the expected value. If the
- * imaginary part is not NaN the operation must also satisfy the conjugate equality.
- *
- * <pre>
- * op(conj(z)) = conj(op(z))
- * </pre>
- *
- * <p>If the function type is ODD/EVEN the operation must satisfy the function type equality.
- *
- * <pre>
- * Odd : f(z) = -f(-z)
- * Even: f(z) = f(-z)
- * </pre>
- *
- * <p>An ODD/EVEN function is also tested that the conjugate equalities hold with {@code -z}.
- * This effectively enumerates testing: (re, im); (re, -im); (-re, -im); and (-re, im).
- *
- * <p>The results must be binary equal; the sign of the complex number is first processed
- * using the provided sign specification.
+ * <p>Assert the operation on the complex number is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts.
*
- * @param z the complex
- * @param operation the operation
- * @param expected the expected
- * @param type the type
- * @param sign the sign specification
+ * @param z Input complex number.
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex object.
+ * @param operation2 Operation on the complex real and imaginary parts.
+ * @param expected Expected complex number.
*/
private static void assertComplex(Complex z,
- UnaryOperator<Complex> operation, Complex expected,
- FunctionType type, UnspecifiedSign sign) {
- // Developer note: Set the sign specification to UnspecifiedSign.NONE
- // to see which equalities fail. They should be for input defined
- // in ISO C99 with an unspecified output sign, e.g.
- // sign = UnspecifiedSign.NONE
-
- // Test the operation
- final Complex c = operation.apply(z);
- final Complex t1 = sign.removeSign(c);
- final Complex t2 = sign.removeSign(expected);
- if (!equals(t1.getReal(), t2.getReal()) ||
- !equals(t1.getImaginary(), t2.getImaginary())) {
- Assertions.fail(
- String.format("Operation failed (z=%s). Expected: %s but was: %s (Unspecified sign = %s)",
- z, expected, c, sign));
- }
-
- if (!Double.isNaN(z.getImaginary())) {
- assertConjugateEquality(z, operation, sign);
- }
-
- if (type != FunctionType.NONE) {
- assertFunctionType(z, operation, type, sign);
-
- // An odd/even function should satisfy the conjugate equality
- // on the negated complex. This ensures testing the equalities
- // hold for:
- // (re, im) = (re, -im)
- // (re, im) = (-re, -im) (even)
- // = -(-re, -im) (odd)
- // (-re, -im) = (-re, im)
- if (!Double.isNaN(z.getImaginary())) {
- assertConjugateEquality(z.negate(), operation, sign);
- }
- }
+ String name, UnaryOperator<Complex> operation1,
+ ComplexUnaryOperator<ComplexNumber> operation2,
+ Complex expected) {
+ assertComplex(z, name, operation1, operation2, expected, FunctionType.NONE, UnspecifiedSign.NONE);
}
/**
@@ -625,12 +453,13 @@ class CStandardTest {
* @param operation1 Operation on the Complex object.
* @param operation2 Operation on the complex real and imaginary parts.
* @param expected Expected complex number.
+ * @param sign the sign specification
*/
private static void assertComplex(Complex z,
String name, UnaryOperator<Complex> operation1,
ComplexUnaryOperator<ComplexNumber> operation2,
- Complex expected) {
- assertComplex(z, name, operation1, operation2, expected, FunctionType.NONE, UnspecifiedSign.NONE);
+ Complex expected, UnspecifiedSign sign) {
+ assertComplex(z, name, operation1, operation2, expected, FunctionType.NONE, sign);
}
/**
@@ -691,7 +520,7 @@ class CStandardTest {
}
if (type != FunctionType.NONE) {
- assertFunctionType(z, operation1, type, sign);
+ assertFunctionType(z, name, operation1, operation2, type, sign);
// An odd/even function should satisfy the conjugate equality
// on the negated complex. This ensures testing the equalities
@@ -706,6 +535,42 @@ class CStandardTest {
}
}
+ /**
+ * Assert the operation on the complex number is equal to the expected value. If the
+ * imaginary part is not NaN the operation must also satisfy the conjugate equality.
+ *
+ * <pre>
+ * op(conj(z)) = conj(op(z))
+ * </pre>
+ *
+ * <p>If the function type is ODD/EVEN the operation must satisfy the function type equality.
+ *
+ * <pre>
+ * Odd : f(z) = -f(-z)
+ * Even: f(z) = f(-z)
+ * </pre>
+ *
+ * <p>The results must be binary equal. This includes the sign of zero.
+ *
+ * <p>Assert the operation on the complex number is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts.
+ *
+ * @param z Input complex number
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex object.
+ * @param operation2 Operation on the complex real and imaginary parts.
+ * @param expected Expected complex number.
+ * @param type Function type.
+ */
+ private static void assertComplex(Complex z,
+ String name,
+ UnaryOperator<Complex> operation1,
+ ComplexUnaryOperator<ComplexNumber> operation2,
+ Complex expected,
+ FunctionType type) {
+ assertComplex(z, name, operation1, operation2, expected, type, UnspecifiedSign.NONE);
+ }
+
/**
* Assert {@link Complex#abs()} is functionally equivalent to using
* {@link Math#hypot(double, double)}. If the results differ the true result
@@ -735,7 +600,7 @@ class CStandardTest {
double y = z.getImaginary();
// For speed use Math.hypot as the reference, not BigDecimal computation.
final double expected = Math.hypot(x, y);
- final double observed = z.abs();
+ final double observed = TestUtils.assertSame(z, "abs", Complex::abs, ComplexFunctions::abs);
if (expected == observed) {
// This condition will occur in the majority of cases.
return;
@@ -799,6 +664,18 @@ class CStandardTest {
}
}
+ /**
+ * Verifies that the expected {@code abs} value is <em>exactly</em> equal to the
+ * actual {@code abs} value which is calculated by the given input argument.
+ *
+ * @param expected Expected double value.
+ * @param z Input complex number.
+ */
+ private static void assertAbs(double expected, Complex z) {
+ final double actual = TestUtils.assertSame(z, "abs", Complex::abs, ComplexFunctions::abs);
+ Assertions.assertEquals(expected, actual);
+ }
+
/**
* Creates a sub-normal number with up to 52-bits in the mantissa. The number of bits
* to drop must be in the range [0, 51].
@@ -853,7 +730,7 @@ class CStandardTest {
*/
private static ArrayList<Complex> createInfinites() {
final double[] values = {0, 1, inf, negInf, nan};
- return createCombinations(values, Complex::isInfinite);
+ return createCombinations(values, "isInfinite", Complex::isInfinite, ComplexFunctions::isInfinite);
}
/**
@@ -883,7 +760,7 @@ class CStandardTest {
*/
private static ArrayList<Complex> createNaNs() {
final double[] values = {0, 1, nan};
- return createCombinations(values, Complex::isNaN);
+ return createCombinations(values, "isNaN", Complex::isNaN, ComplexFunctions::isNaN);
}
/**
@@ -907,6 +784,35 @@ class CStandardTest {
return list;
}
+ /**
+ * Creates a list of Complex numbers as an all-vs-all combinations that pass the
+ * condition.
+ *
+ * <p>Assert the condition operation on the complex number is <em>exactly</em> equal to the condition
+ * operation on complex real and imaginary parts.
+ *
+ * @param values the values
+ * @param name Operation name.
+ * @param operation1 Condition operation on the Complex object.
+ * @param operation2 Condition operation on the complex real and imaginary parts.
+ * @return the list
+ */
+ private static ArrayList<Complex> createCombinations(final double[] values,
+ String name,
+ Predicate<Complex> operation1,
+ TestUtils.DoubleBinaryPredicate operation2) {
+ final ArrayList<Complex> list = new ArrayList<>();
+ for (final double re : values) {
+ for (final double im : values) {
+ final Complex z = complex(re, im);
+ if (TestUtils.assertSame(z, name, operation1, operation2)) {
+ list.add(z);
+ }
+ }
+ }
+ return list;
+ }
+
/**
* Checks if the complex is zero. This method uses the {@code ==} operator and allows
* equality between signed zeros: {@code -0.0 == 0.0}.
@@ -1049,8 +955,13 @@ class CStandardTest {
*/
@Test
void testSqrt1() {
- assertComplex(complex(-2.0, 0.0).sqrt(), complex(0.0, Math.sqrt(2)));
- assertComplex(complex(-2.0, -0.0).sqrt(), complex(0.0, -Math.sqrt(2)));
+ Complex z = complex(-2.0, 0.0);
+ Complex actual = TestUtils.assertSame(z, "sqrt", Complex::sqrt, ComplexFunctions::sqrt);
+ assertComplex(actual, complex(0.0, Math.sqrt(2)));
+
+ z = complex(-2.0, -0.0);
+ actual = TestUtils.assertSame(z, "sqrt", Complex::sqrt, ComplexFunctions::sqrt);
+ assertComplex(actual, complex(0.0, -Math.sqrt(2)));
}
/**
@@ -1064,11 +975,26 @@ class CStandardTest {
final double im = next(rng);
final Complex z = complex(re, im);
final Complex iz = z.multiplyImaginary(1);
- assertComplex(z.asin(), iz.asinh().multiplyImaginary(-1));
- assertComplex(z.atan(), iz.atanh().multiplyImaginary(-1));
- assertComplex(z.cos(), iz.cosh());
- assertComplex(z.sin(), iz.sinh().multiplyImaginary(-1));
- assertComplex(z.tan(), iz.tanh().multiplyImaginary(-1));
+
+ Complex actual = TestUtils.assertSame(z, "asin", Complex::asin, ComplexFunctions::asin);
+ Complex expected = TestUtils.assertSame(iz, "asinh", Complex::asinh, ComplexFunctions::asinh);
+ assertComplex(actual, expected.multiplyImaginary(-1));
+
+ actual = TestUtils.assertSame(z, "atan", Complex::atan, ComplexFunctions::atan);
+ expected = TestUtils.assertSame(iz, "atanh", Complex::atanh, ComplexFunctions::atanh);
+ assertComplex(actual, expected.multiplyImaginary(-1));
+
+ actual = TestUtils.assertSame(z, "cos", Complex::cos, ComplexFunctions::cos);
+ expected = TestUtils.assertSame(iz, "cosh", Complex::cosh, ComplexFunctions::cosh);
+ assertComplex(actual, expected);
+
+ actual = TestUtils.assertSame(z, "sin", Complex::sin, ComplexFunctions::sin);
+ expected = TestUtils.assertSame(iz, "sinh", Complex::sinh, ComplexFunctions::sinh);
+ assertComplex(actual, expected.multiplyImaginary(-1));
+
+ actual = TestUtils.assertSame(z, "tan", Complex::tan, ComplexFunctions::tan);
+ expected = TestUtils.assertSame(iz, "tanh", Complex::tanh, ComplexFunctions::tanh);
+ assertComplex(actual, expected.multiplyImaginary(-1));
}
}
@@ -1089,18 +1015,18 @@ class CStandardTest {
*/
@Test
void testAbs() {
- Assertions.assertEquals(inf, complex(inf, nan).abs());
- Assertions.assertEquals(inf, complex(negInf, nan).abs());
+ assertAbs(inf, complex(inf, nan));
+ assertAbs(inf, complex(negInf, nan));
final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create();
for (int i = 0; i < 10; i++) {
final double x = next(rng);
final double y = next(rng);
- Assertions.assertEquals(complex(x, y).abs(), complex(y, x).abs());
- Assertions.assertEquals(complex(x, y).abs(), complex(x, -y).abs());
- Assertions.assertEquals(Math.abs(x), complex(x, 0.0).abs());
- Assertions.assertEquals(Math.abs(x), complex(x, -0.0).abs());
- Assertions.assertEquals(inf, complex(inf, y).abs());
- Assertions.assertEquals(inf, complex(negInf, y).abs());
+ assertAbs(complex(x, y).abs(), complex(y, x));
+ assertAbs(complex(x, y).abs(), complex(x, -y));
+ assertAbs(Math.abs(x), complex(x, 0.0));
+ assertAbs(Math.abs(x), complex(x, -0.0));
+ assertAbs(inf, complex(inf, y));
+ assertAbs(inf, complex(negInf, y));
}
// Test verses Math.hypot due to the use of a custom implementation.
@@ -1109,7 +1035,7 @@ class CStandardTest {
Double.NEGATIVE_INFINITY, Double.NaN};
for (final double x : parts) {
for (final double y : parts) {
- Assertions.assertEquals(Math.hypot(x, y), complex(x, y).abs());
+ assertAbs(Math.hypot(x, y), complex(x, y));
}
}
@@ -1124,6 +1050,7 @@ class CStandardTest {
{1.40905821964671, 1.4090583434236112},
{1.912164268932753, 1.9121638616231227}}) {
final Complex z = complex(pair[0], pair[1]);
+ assertAbs(z.abs(), z.multiplyImaginary(1));
Assertions.assertEquals(z.abs(), z.multiplyImaginary(1).abs(), "Expected |z| == |iz|");
}
@@ -1174,33 +1101,35 @@ class CStandardTest {
*/
@Test
void testAcos() {
- final UnaryOperator<Complex> operation = Complex::acos;
- assertConjugateEquality(operation);
- assertComplex(Complex.ZERO, operation, piTwoNegZero);
- assertComplex(negZeroZero, operation, piTwoNegZero);
- assertComplex(zeroNaN, operation, piTwoNaN);
- assertComplex(negZeroNaN, operation, piTwoNaN);
+ final String name = "acos";
+ final UnaryOperator<Complex> operation1 = Complex::acos;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::acos;
+ assertConjugateEquality(name, operation1, operation2);
+ assertComplex(Complex.ZERO, name, operation1, operation2, piTwoNegZero);
+ assertComplex(negZeroZero, name, operation1, operation2, piTwoNegZero);
+ assertComplex(zeroNaN, name, operation1, operation2, piTwoNaN);
+ assertComplex(negZeroNaN, name, operation1, operation2, piTwoNaN);
for (double x : finite) {
- assertComplex(complex(x, inf), operation, piTwoNegInf);
+ assertComplex(complex(x, inf), name, operation1, operation2, piTwoNegInf);
}
for (double x : nonZeroFinite) {
- assertComplex(complex(x, nan), operation, NAN);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
for (double y : positiveFinite) {
- assertComplex(complex(-inf, y), operation, piNegInf);
+ assertComplex(complex(-inf, y), name, operation1, operation2, piNegInf);
}
for (double y : positiveFinite) {
- assertComplex(complex(inf, y), operation, zeroNegInf);
+ assertComplex(complex(inf, y), name, operation1, operation2, zeroNegInf);
}
- assertComplex(negInfInf, operation, threePiFourNegInf);
- assertComplex(infInf, operation, piFourNegInf);
- assertComplex(infNaN, operation, nanInf, UnspecifiedSign.IMAGINARY);
- assertComplex(negInfNaN, operation, nanNegInf, UnspecifiedSign.IMAGINARY);
+ assertComplex(negInfInf, name, operation1, operation2, threePiFourNegInf);
+ assertComplex(infInf, name, operation1, operation2, piFourNegInf);
+ assertComplex(infNaN, name, operation1, operation2, nanInf, UnspecifiedSign.IMAGINARY);
+ assertComplex(negInfNaN, name, operation1, operation2, nanNegInf, UnspecifiedSign.IMAGINARY);
for (double y : finite) {
- assertComplex(complex(nan, y), operation, NAN);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
- assertComplex(nanInf, operation, nanNegInf);
- assertComplex(NAN, operation, NAN);
+ assertComplex(nanInf, name, operation1, operation2, nanNegInf);
+ assertComplex(NAN, name, operation1, operation2, NAN);
}
/**
@@ -1211,33 +1140,36 @@ class CStandardTest {
*/
@Test
void testAcosh() {
- final UnaryOperator<Complex> operation = Complex::acosh;
- assertConjugateEquality(operation);
- assertComplex(Complex.ZERO, operation, zeroPiTwo);
- assertComplex(negZeroZero, operation, zeroPiTwo);
+ final String name = "acosh";
+ final UnaryOperator<Complex> operation1 = Complex::acosh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::acosh;
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertComplex(Complex.ZERO, name, operation1, operation2, zeroPiTwo);
+ assertComplex(negZeroZero, name, operation1, operation2, zeroPiTwo);
for (double x : finite) {
- assertComplex(complex(x, inf), operation, infPiTwo);
+ assertComplex(complex(x, inf), name, operation1, operation2, infPiTwo);
}
- assertComplex(zeroNaN, operation, nanPiTwo);
- assertComplex(negZeroNaN, operation, nanPiTwo);
+ assertComplex(zeroNaN, name, operation1, operation2, nanPiTwo);
+ assertComplex(negZeroNaN, name, operation1, operation2, nanPiTwo);
for (double x : nonZeroFinite) {
- assertComplex(complex(x, nan), operation, NAN);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
for (double y : positiveFinite) {
- assertComplex(complex(-inf, y), operation, infPi);
+ assertComplex(complex(-inf, y), name, operation1, operation2, infPi);
}
for (double y : positiveFinite) {
- assertComplex(complex(inf, y), operation, infZero);
+ assertComplex(complex(inf, y), name, operation1, operation2, infZero);
}
- assertComplex(negInfInf, operation, infThreePiFour);
- assertComplex(infInf, operation, infPiFour);
- assertComplex(infNaN, operation, infNaN);
- assertComplex(negInfNaN, operation, infNaN);
+ assertComplex(negInfInf, name, operation1, operation2, infThreePiFour);
+ assertComplex(infInf, name, operation1, operation2, infPiFour);
+ assertComplex(infNaN, name, operation1, operation2, infNaN);
+ assertComplex(negInfNaN, name, operation1, operation2, infNaN);
for (double y : finite) {
- assertComplex(complex(nan, y), operation, NAN);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
- assertComplex(nanInf, operation, infNaN);
- assertComplex(NAN, operation, NAN);
+ assertComplex(nanInf, name, operation1, operation2, infNaN);
+ assertComplex(NAN, name, operation1, operation2, NAN);
}
/**
@@ -1245,28 +1177,31 @@ class CStandardTest {
*/
@Test
void testAsinh() {
- final UnaryOperator<Complex> operation = Complex::asinh;
+ final String name = "asinh";
+ final UnaryOperator<Complex> operation1 = Complex::asinh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::asinh;
final FunctionType type = FunctionType.ODD;
- assertConjugateEquality(operation);
- assertFunctionType(operation, type);
- assertComplex(Complex.ZERO, operation, Complex.ZERO, type);
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertFunctionType(name, operation1, operation2, type);
+ assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ZERO, type);
for (double x : positiveFinite) {
- assertComplex(complex(x, inf), operation, infPiTwo, type);
+ assertComplex(complex(x, inf), name, operation1, operation2, infPiTwo, type);
}
for (double x : finite) {
- assertComplex(complex(x, nan), operation, NAN, type);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
for (double y : positiveFinite) {
- assertComplex(complex(inf, y), operation, infZero, type);
+ assertComplex(complex(inf, y), name, operation1, operation2, infZero, type);
}
- assertComplex(infInf, operation, infPiFour, type);
- assertComplex(infNaN, operation, infNaN, type);
- assertComplex(nanZero, operation, nanZero, type);
+ assertComplex(infInf, name, operation1, operation2, infPiFour, type);
+ assertComplex(infNaN, name, operation1, operation2, infNaN, type);
+ assertComplex(nanZero, name, operation1, operation2, nanZero, type);
for (double y : nonZeroFinite) {
- assertComplex(complex(nan, y), operation, NAN, type);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
- assertComplex(nanInf, operation, infNaN, type, UnspecifiedSign.REAL);
- assertComplex(NAN, operation, NAN, type);
+ assertComplex(nanInf, name, operation1, operation2, infNaN, type, UnspecifiedSign.REAL);
+ assertComplex(NAN, name, operation1, operation2, NAN, type);
}
/**
@@ -1274,30 +1209,33 @@ class CStandardTest {
*/
@Test
void testAtanh() {
- final UnaryOperator<Complex> operation = Complex::atanh;
+ final String name = "atanh";
+ final UnaryOperator<Complex> operation1 = Complex::atanh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::atanh;
final FunctionType type = FunctionType.ODD;
- assertConjugateEquality(operation);
- assertFunctionType(operation, type);
- assertComplex(Complex.ZERO, operation, Complex.ZERO, type);
- assertComplex(zeroNaN, operation, zeroNaN, type);
- assertComplex(oneZero, operation, infZero, type);
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertFunctionType(name, operation1, operation2, type);
+ assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ZERO, type);
+ assertComplex(zeroNaN, name, operation1, operation2, zeroNaN, type);
+ assertComplex(oneZero, name, operation1, operation2, infZero, type);
for (double x : positiveFinite) {
- assertComplex(complex(x, inf), operation, zeroPiTwo, type);
+ assertComplex(complex(x, inf), name, operation1, operation2, zeroPiTwo, type);
}
for (double x : nonZeroFinite) {
- assertComplex(complex(x, nan), operation, NAN, type);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
for (double y : positiveFinite) {
- assertComplex(complex(inf, y), operation, zeroPiTwo, type);
+ assertComplex(complex(inf, y), name, operation1, operation2, zeroPiTwo, type);
}
- assertComplex(infInf, operation, zeroPiTwo, type);
- assertComplex(infNaN, operation, zeroNaN, type);
- assertComplex(nanZero, operation, NAN, type);
+ assertComplex(infInf, name, operation1, operation2, zeroPiTwo, type);
+ assertComplex(infNaN, name, operation1, operation2, zeroNaN, type);
+ assertComplex(nanZero, name, operation1, operation2, NAN, type);
for (double y : finite) {
- assertComplex(complex(nan, y), operation, NAN, type);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
- assertComplex(nanInf, operation, zeroPiTwo, type, UnspecifiedSign.REAL);
- assertComplex(NAN, operation, NAN, type);
+ assertComplex(nanInf, name, operation1, operation2, zeroPiTwo, type, UnspecifiedSign.REAL);
+ assertComplex(NAN, name, operation1, operation2, NAN, type);
}
/**
@@ -1305,30 +1243,33 @@ class CStandardTest {
*/
@Test
void testCosh() {
- final UnaryOperator<Complex> operation = Complex::cosh;
+ final String name = "cosh";
+ final UnaryOperator<Complex> operation1 = Complex::cosh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::cosh;
final FunctionType type = FunctionType.EVEN;
- assertConjugateEquality(operation);
- assertFunctionType(operation, type);
- assertComplex(Complex.ZERO, operation, Complex.ONE, type);
- assertComplex(zeroInf, operation, nanZero, type, UnspecifiedSign.IMAGINARY);
- assertComplex(zeroNaN, operation, nanZero, type, UnspecifiedSign.IMAGINARY);
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertFunctionType(name, operation1, operation2, type);
+ assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ONE, type);
+ assertComplex(zeroInf, name, operation1, operation2, nanZero, type, UnspecifiedSign.IMAGINARY);
+ assertComplex(zeroNaN, name, operation1, operation2, nanZero, type, UnspecifiedSign.IMAGINARY);
for (double x : nonZeroFinite) {
- assertComplex(complex(x, inf), operation, NAN, type);
+ assertComplex(complex(x, inf), name, operation1, operation2, NAN, type);
}
for (double x : nonZeroFinite) {
- assertComplex(complex(x, nan), operation, NAN, type);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
- assertComplex(infZero, operation, infZero, type);
+ assertComplex(infZero, name, operation1, operation2, infZero, type);
for (double y : nonZeroFinite) {
- assertComplex(complex(inf, y), operation, Complex.ofCis(y).multiply(inf), type);
+ assertComplex(complex(inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(inf), type);
}
- assertComplex(infInf, operation, infNaN, type, UnspecifiedSign.REAL);
- assertComplex(infNaN, operation, infNaN, type);
- assertComplex(nanZero, operation, nanZero, type, UnspecifiedSign.IMAGINARY);
+ assertComplex(infInf, name, operation1, operation2, infNaN, type, UnspecifiedSign.REAL);
+ assertComplex(infNaN, name, operation1, operation2, infNaN, type);
+ assertComplex(nanZero, name, operation1, operation2, nanZero, type, UnspecifiedSign.IMAGINARY);
for (double y : nonZero) {
- assertComplex(complex(nan, y), operation, NAN, type);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
- assertComplex(NAN, operation, NAN, type);
+ assertComplex(NAN, name, operation1, operation2, NAN, type);
}
/**
@@ -1336,31 +1277,34 @@ class CStandardTest {
*/
@Test
void testSinh() {
- final UnaryOperator<Complex> operation = Complex::sinh;
+ final String name = "sinh";
+ final UnaryOperator<Complex> operation1 = Complex::sinh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::sinh;
final FunctionType type = FunctionType.ODD;
- assertConjugateEquality(operation);
- assertFunctionType(operation, type);
- assertComplex(Complex.ZERO, operation, Complex.ZERO, type);
- assertComplex(zeroInf, operation, zeroNaN, type, UnspecifiedSign.REAL);
- assertComplex(zeroNaN, operation, zeroNaN, type, UnspecifiedSign.REAL);
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertFunctionType(name, operation1, operation2, type);
+ assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ZERO, type);
+ assertComplex(zeroInf, name, operation1, operation2, zeroNaN, type, UnspecifiedSign.REAL);
+ assertComplex(zeroNaN, name, operation1, operation2, zeroNaN, type, UnspecifiedSign.REAL);
for (double x : nonZeroPositiveFinite) {
- assertComplex(complex(x, inf), operation, NAN, type);
+ assertComplex(complex(x, inf), name, operation1, operation2, NAN, type);
}
for (double x : nonZeroFinite) {
- assertComplex(complex(x, nan), operation, NAN, type);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
- assertComplex(infZero, operation, infZero, type);
+ assertComplex(infZero, name, operation1, operation2, infZero, type);
// Note: Error in the ISO C99 reference to use positive finite y but the zero case is different
for (double y : nonZeroFinite) {
- assertComplex(complex(inf, y), operation, Complex.ofCis(y).multiply(inf), type);
+ assertComplex(complex(inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(inf), type);
}
- assertComplex(infInf, operation, infNaN, type, UnspecifiedSign.REAL);
- assertComplex(infNaN, operation, infNaN, type, UnspecifiedSign.REAL);
- assertComplex(nanZero, operation, nanZero, type);
+ assertComplex(infInf, name, operation1, operation2, infNaN, type, UnspecifiedSign.REAL);
+ assertComplex(infNaN, name, operation1, operation2, infNaN, type, UnspecifiedSign.REAL);
+ assertComplex(nanZero, name, operation1, operation2, nanZero, type);
for (double y : nonZero) {
- assertComplex(complex(nan, y), operation, NAN, type);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
- assertComplex(NAN, operation, NAN, type);
+ assertComplex(NAN, name, operation1, operation2, NAN, type);
}
/**
@@ -1371,29 +1315,32 @@ class CStandardTest {
*/
@Test
void testTanh() {
- final UnaryOperator<Complex> operation = Complex::tanh;
+ final String name = "tanh";
+ final UnaryOperator<Complex> operation1 = Complex::tanh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::tanh;
final FunctionType type = FunctionType.ODD;
- assertConjugateEquality(operation);
- assertFunctionType(operation, type);
- assertComplex(Complex.ZERO, operation, Complex.ZERO, type);
- assertComplex(zeroInf, operation, zeroNaN, type);
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertFunctionType(name, operation1, operation2, type);
+ assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ZERO, type);
+ assertComplex(zeroInf, name, operation1, operation2, zeroNaN, type);
for (double x : nonZeroFinite) {
- assertComplex(complex(x, inf), operation, NAN, type);
+ assertComplex(complex(x, inf), name, operation1, operation2, NAN, type);
}
- assertComplex(zeroNaN, operation, zeroNaN, type);
+ assertComplex(zeroNaN, name, operation1, operation2, zeroNaN, type);
for (double x : nonZeroFinite) {
- assertComplex(complex(x, nan), operation, NAN, type);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
for (double y : positiveFinite) {
- assertComplex(complex(inf, y), operation, complex(1.0, Math.copySign(0, Math.sin(2 * y))), type);
+ assertComplex(complex(inf, y), name, operation1, operation2, complex(1.0, Math.copySign(0, Math.sin(2 * y))), type);
}
- assertComplex(infInf, operation, oneZero, type, UnspecifiedSign.IMAGINARY);
- assertComplex(infNaN, operation, oneZero, type, UnspecifiedSign.IMAGINARY);
- assertComplex(nanZero, operation, nanZero, type);
+ assertComplex(infInf, name, operation1, operation2, oneZero, type, UnspecifiedSign.IMAGINARY);
+ assertComplex(infNaN, name, operation1, operation2, oneZero, type, UnspecifiedSign.IMAGINARY);
+ assertComplex(nanZero, name, operation1, operation2, nanZero, type);
for (double y : nonZero) {
- assertComplex(complex(nan, y), operation, NAN, type);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
- assertComplex(NAN, operation, NAN, type);
+ assertComplex(NAN, name, operation1, operation2, NAN, type);
}
/**
@@ -1401,32 +1348,35 @@ class CStandardTest {
*/
@Test
void testExp() {
- final UnaryOperator<Complex> operation = Complex::exp;
- assertConjugateEquality(operation);
- assertComplex(Complex.ZERO, operation, oneZero);
- assertComplex(negZeroZero, operation, oneZero);
+ final String name = "exp";
+ final UnaryOperator<Complex> operation1 = Complex::exp;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::exp;
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertComplex(Complex.ZERO, name, operation1, operation2, oneZero);
+ assertComplex(negZeroZero, name, operation1, operation2, oneZero);
for (double x : finite) {
- assertComplex(complex(x, inf), operation, NAN);
+ assertComplex(complex(x, inf), name, operation1, operation2, NAN);
}
for (double x : finite) {
- assertComplex(complex(x, nan), operation, NAN);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
- assertComplex(infZero, operation, infZero);
+ assertComplex(infZero, name, operation1, operation2, infZero);
for (double y : finite) {
- assertComplex(complex(-inf, y), operation, Complex.ofCis(y).multiply(0.0));
+ assertComplex(complex(-inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(0.0));
}
for (double y : nonZeroFinite) {
- assertComplex(complex(inf, y), operation, Complex.ofCis(y).multiply(inf));
+ assertComplex(complex(inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(inf));
}
- assertComplex(negInfInf, operation, Complex.ZERO, UnspecifiedSign.REAL_IMAGINARY);
- assertComplex(infInf, operation, infNaN, UnspecifiedSign.REAL);
- assertComplex(negInfNaN, operation, Complex.ZERO, UnspecifiedSign.REAL_IMAGINARY);
- assertComplex(infNaN, operation, infNaN, UnspecifiedSign.REAL);
- assertComplex(nanZero, operation, nanZero);
+ assertComplex(negInfInf, name, operation1, operation2, Complex.ZERO, UnspecifiedSign.REAL_IMAGINARY);
+ assertComplex(infInf, name, operation1, operation2, infNaN, UnspecifiedSign.REAL);
+ assertComplex(negInfNaN, name, operation1, operation2, Complex.ZERO, UnspecifiedSign.REAL_IMAGINARY);
+ assertComplex(infNaN, name, operation1, operation2, infNaN, UnspecifiedSign.REAL);
+ assertComplex(nanZero, name, operation1, operation2, nanZero);
for (double y : nonZero) {
- assertComplex(complex(nan, y), operation, NAN);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
- assertComplex(NAN, operation, NAN);
+ assertComplex(NAN, name, operation1, operation2, NAN);
}
/**
@@ -1505,31 +1455,34 @@ class CStandardTest {
*/
@Test
void testSqrt() {
- final UnaryOperator<Complex> operation = Complex::sqrt;
- assertConjugateEquality(operation);
- assertComplex(negZeroZero, operation, Complex.ZERO);
- assertComplex(Complex.ZERO, operation, Complex.ZERO);
+ final String name = "sqrt";
+ final UnaryOperator<Complex> operation1 = Complex::sqrt;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::sqrt;
+
+ assertConjugateEquality(name, operation1, operation2);
+ assertComplex(negZeroZero, name, operation1, operation2, Complex.ZERO);
+ assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ZERO);
for (double x : finite) {
- assertComplex(complex(x, inf), operation, infInf);
+ assertComplex(complex(x, inf), name, operation1, operation2, infInf);
}
// Include infinity and nan for (x, inf).
- assertComplex(infInf, operation, infInf);
- assertComplex(negInfInf, operation, infInf);
- assertComplex(nanInf, operation, infInf);
+ assertComplex(infInf, name, operation1, operation2, infInf);
+ assertComplex(negInfInf, name, operation1, operation2, infInf);
+ assertComplex(nanInf, name, operation1, operation2, infInf);
for (double x : finite) {
- assertComplex(complex(x, nan), operation, NAN);
+ assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
for (double y : positiveFinite) {
- assertComplex(complex(-inf, y), operation, zeroInf);
+ assertComplex(complex(-inf, y), name, operation1, operation2, zeroInf);
}
for (double y : positiveFinite) {
- assertComplex(complex(inf, y), operation, infZero);
+ assertComplex(complex(inf, y), name, operation1, operation2, infZero);
}
- assertComplex(negInfNaN, operation, nanInf, UnspecifiedSign.IMAGINARY);
- assertComplex(infNaN, operation, infNaN);
+ assertComplex(negInfNaN, name, operation1, operation2, nanInf, UnspecifiedSign.IMAGINARY);
+ assertComplex(infNaN, name, operation1, operation2, infNaN);
for (double y : finite) {
- assertComplex(complex(nan, y), operation, NAN);
+ assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
- assertComplex(NAN, operation, NAN);
+ assertComplex(NAN, name, operation1, operation2, NAN);
}
}
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
index 3fd467ff..dfd53f45 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
@@ -182,7 +182,8 @@ class ComplexEdgeCaseTest {
void testAcos() {
// acos(z) = (pi / 2) + i ln(iz + sqrt(1 - z^2))
final String name = "acos";
- final UnaryOperator<Complex> operation = Complex::acos;
+ final UnaryOperator<Complex> operation1 = Complex::acos;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::acos;
// Edge cases are when values are big but not infinite and small but not zero.
// Big and small are set using the limits in atanh.
@@ -193,23 +194,23 @@ class ComplexEdgeCaseTest {
final double big = Math.sqrt(Double.MAX_VALUE) / 8;
final double medium = 100;
final double small = Math.sqrt(Double.MIN_NORMAL) * 4;
- assertComplex(huge, big, name, operation, 0.06241880999595735, -356.27960012801969);
- assertComplex(huge, medium, name, operation, 3.7291703656001039e-153, -356.27765080781188);
- assertComplex(huge, small, name, operation, 2.2250738585072019e-308, -356.27765080781188);
- assertComplex(big, big, name, operation, 0.78539816339744828, -353.85163567585209);
- assertComplex(big, medium, name, operation, 5.9666725849601662e-152, -353.50506208557209);
- assertComplex(big, small, name, operation, 3.560118173611523e-307, -353.50506208557209);
- assertComplex(medium, big, name, operation, 1.5707963267948966, -353.50506208557209);
- assertComplex(medium, medium, name, operation, 0.78541066339744181, -5.6448909570623842);
- assertComplex(medium, small, name, operation, 5.9669709409662999e-156, -5.298292365610485);
- assertComplex(small, big, name, operation, 1.5707963267948966, -353.50506208557209);
- assertComplex(small, medium, name, operation, 1.5707963267948966, -5.2983423656105888);
- assertComplex(small, small, name, operation, 1.5707963267948966, -5.9666725849601654e-154);
+ assertComplex(huge, big, name, operation1, operation2, 0.06241880999595735, -356.27960012801969);
+ assertComplex(huge, medium, name, operation1, operation2, 3.7291703656001039e-153, -356.27765080781188);
+ assertComplex(huge, small, name, operation1, operation2, 2.2250738585072019e-308, -356.27765080781188);
+ assertComplex(big, big, name, operation1, operation2, 0.78539816339744828, -353.85163567585209);
+ assertComplex(big, medium, name, operation1, operation2, 5.9666725849601662e-152, -353.50506208557209);
+ assertComplex(big, small, name, operation1, operation2, 3.560118173611523e-307, -353.50506208557209);
+ assertComplex(medium, big, name, operation1, operation2, 1.5707963267948966, -353.50506208557209);
+ assertComplex(medium, medium, name, operation1, operation2, 0.78541066339744181, -5.6448909570623842);
+ assertComplex(medium, small, name, operation1, operation2, 5.9669709409662999e-156, -5.298292365610485);
+ assertComplex(small, big, name, operation1, operation2, 1.5707963267948966, -353.50506208557209);
+ assertComplex(small, medium, name, operation1, operation2, 1.5707963267948966, -5.2983423656105888);
+ assertComplex(small, small, name, operation1, operation2, 1.5707963267948966, -5.9666725849601654e-154);
// Additional cases to achieve full coverage
// xm1 = 0
- assertComplex(1, small, name, operation, 2.4426773395109241e-77, -2.4426773395109241e-77);
+ assertComplex(1, small, name, operation1, operation2, 2.4426773395109241e-77, -2.4426773395109241e-77);
// https://svn.boost.org/trac10/ticket/7290
- assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation, 2.4252018043912224e-196, -0.00023604834149293664);
+ assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation1, operation2, 2.4252018043912224e-196, -0.00023604834149293664);
}
// acosh is defined by acos so is not tested
@@ -218,7 +219,8 @@ class ComplexEdgeCaseTest {
void testAsin() {
// asin(z) = -i (ln(iz + sqrt(1 - z^2)))
final String name = "asin";
- final UnaryOperator<Complex> operation = Complex::asin;
+ final UnaryOperator<Complex> operation1 = Complex::asin;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::asin;
// This method is essentially the same as acos and the edge cases are the same.
// The results have been generated using g++ -std=c++11 asin.
@@ -226,23 +228,23 @@ class ComplexEdgeCaseTest {
final double big = Math.sqrt(Double.MAX_VALUE) / 8;
final double medium = 100;
final double small = Math.sqrt(Double.MIN_NORMAL) * 4;
- assertComplex(huge, big, name, operation, 1.5083775167989393, 356.27960012801969);
- assertComplex(huge, medium, name, operation, 1.5707963267948966, 356.27765080781188);
- assertComplex(huge, small, name, operation, 1.5707963267948966, 356.27765080781188);
- assertComplex(big, big, name, operation, 0.78539816339744828, 353.85163567585209);
- assertComplex(big, medium, name, operation, 1.5707963267948966, 353.50506208557209);
- assertComplex(big, small, name, operation, 1.5707963267948966, 353.50506208557209);
- assertComplex(medium, big, name, operation, 5.9666725849601662e-152, 353.50506208557209);
- assertComplex(medium, medium, name, operation, 0.78538566339745486, 5.6448909570623842);
- assertComplex(medium, small, name, operation, 1.5707963267948966, 5.298292365610485);
- assertComplex(small, big, name, operation, 3.560118173611523e-307, 353.50506208557209);
- assertComplex(small, medium, name, operation, 5.9663742737040751e-156, 5.2983423656105888);
- assertComplex(small, small, name, operation, 5.9666725849601654e-154, 5.9666725849601654e-154);
+ assertComplex(huge, big, name, operation1, operation2, 1.5083775167989393, 356.27960012801969);
+ assertComplex(huge, medium, name, operation1, operation2, 1.5707963267948966, 356.27765080781188);
+ assertComplex(huge, small, name, operation1, operation2, 1.5707963267948966, 356.27765080781188);
+ assertComplex(big, big, name, operation1, operation2, 0.78539816339744828, 353.85163567585209);
+ assertComplex(big, medium, name, operation1, operation2, 1.5707963267948966, 353.50506208557209);
+ assertComplex(big, small, name, operation1, operation2, 1.5707963267948966, 353.50506208557209);
+ assertComplex(medium, big, name, operation1, operation2, 5.9666725849601662e-152, 353.50506208557209);
+ assertComplex(medium, medium, name, operation1, operation2, 0.78538566339745486, 5.6448909570623842);
+ assertComplex(medium, small, name, operation1, operation2, 1.5707963267948966, 5.298292365610485);
+ assertComplex(small, big, name, operation1, operation2, 3.560118173611523e-307, 353.50506208557209);
+ assertComplex(small, medium, name, operation1, operation2, 5.9663742737040751e-156, 5.2983423656105888);
+ assertComplex(small, small, name, operation1, operation2, 5.9666725849601654e-154, 5.9666725849601654e-154);
// Additional cases to achieve full coverage
// xm1 = 0
- assertComplex(1, small, name, operation, 1.5707963267948966, 2.4426773395109241e-77);
+ assertComplex(1, small, name, operation1, operation2, 1.5707963267948966, 2.4426773395109241e-77);
// https://svn.boost.org/trac10/ticket/7290
- assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation, 1.5707963267948966, 0.00023604834149293664);
+ assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation1, operation2, 1.5707963267948966, 0.00023604834149293664);
}
// asinh is defined by asin so is not tested
@@ -252,7 +254,8 @@ class ComplexEdgeCaseTest {
// atanh(z) = (1/2) ln((1 + z) / (1 - z))
// Odd function: negative real cases defined by positive real cases
final String name = "atanh";
- final UnaryOperator<Complex> operation = Complex::atanh;
+ final UnaryOperator<Complex> operation1 = Complex::atanh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::atanh;
// Edge cases are when values are big but not infinite and small but not zero.
// Big and small are set using the limits in atanh.
@@ -262,18 +265,18 @@ class ComplexEdgeCaseTest {
final double big = Math.sqrt(Double.MAX_VALUE) / 2;
final double medium = 100;
final double small = Math.sqrt(Double.MIN_NORMAL) * 2;
- assertComplex(big, big, name, operation, 7.4583407312002067e-155, 1.5707963267948966);
- assertComplex(big, medium, name, operation, 1.4916681462400417e-154, 1.5707963267948966);
- assertComplex(big, small, name, operation, 1.4916681462400417e-154, 1.5707963267948966);
- assertComplex(medium, big, name, operation, 2.225073858507202e-306, 1.5707963267948966);
- assertComplex(medium, medium, name, operation, 0.0049999166641667555, 1.5657962434640633);
- assertComplex(medium, small, name, operation, 0.010000333353334761, 1.5707963267948966);
- assertComplex(small, big, name, operation, 0, 1.5707963267948966);
- assertComplex(small, medium, name, operation, 2.9830379886812147e-158, 1.5607966601082315);
- assertComplex(small, small, name, operation, 2.9833362924800827e-154, 2.9833362924800827e-154);
+ assertComplex(big, big, name, operation1, operation2, 7.4583407312002067e-155, 1.5707963267948966);
+ assertComplex(big, medium, name, operation1, operation2, 1.4916681462400417e-154, 1.5707963267948966);
+ assertComplex(big, small, name, operation1, operation2, 1.4916681462400417e-154, 1.5707963267948966);
+ assertComplex(medium, big, name, operation1, operation2, 2.225073858507202e-306, 1.5707963267948966);
+ assertComplex(medium, medium, name, operation1, operation2, 0.0049999166641667555, 1.5657962434640633);
+ assertComplex(medium, small, name, operation1, operation2, 0.010000333353334761, 1.5707963267948966);
+ assertComplex(small, big, name, operation1, operation2, 0, 1.5707963267948966);
+ assertComplex(small, medium, name, operation1, operation2, 2.9830379886812147e-158, 1.5607966601082315);
+ assertComplex(small, small, name, operation1, operation2, 2.9833362924800827e-154, 2.9833362924800827e-154);
// Additional cases to achieve full coverage
- assertComplex(inf, big, name, operation, 0, 1.5707963267948966);
- assertComplex(big, inf, name, operation, 0, 1.5707963267948966);
+ assertComplex(inf, big, name, operation1, operation2, 0, 1.5707963267948966);
+ assertComplex(big, inf, name, operation1, operation2, 0, 1.5707963267948966);
}
@Test
@@ -281,22 +284,23 @@ class ComplexEdgeCaseTest {
// cosh(a + b i) = cosh(a)cos(b) + i sinh(a)sin(b)
// Even function: negative real cases defined by positive real cases
final String name = "cosh";
- final UnaryOperator<Complex> operation = Complex::cosh;
+ final UnaryOperator<Complex> operation1 = Complex::cosh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::cosh;
// Implementation defers to java.util.Math.
// Hit edge cases for extreme values.
final double big = Double.MAX_VALUE;
final double medium = 2;
final double small = Double.MIN_NORMAL;
- assertComplex(big, big, name, operation, -inf, inf);
- assertComplex(big, medium, name, operation, -inf, inf);
- assertComplex(big, small, name, operation, inf, inf);
- assertComplex(medium, big, name, operation, -3.7621493762972804, 0.017996317370418576);
- assertComplex(medium, medium, name, operation, -1.5656258353157435, 3.297894836311237);
- assertComplex(medium, small, name, operation, 3.7621956910836314, 8.0700322819551687e-308);
- assertComplex(small, big, name, operation, -0.99998768942655991, 1.1040715888508271e-310);
- assertComplex(small, medium, name, operation, -0.41614683654714241, 2.0232539340376892e-308);
- assertComplex(small, small, name, operation, 1, 0);
+ assertComplex(big, big, name, operation1, operation2, -inf, inf);
+ assertComplex(big, medium, name, operation1, operation2, -inf, inf);
+ assertComplex(big, small, name, operation1, operation2, inf, inf);
+ assertComplex(medium, big, name, operation1, operation2, -3.7621493762972804, 0.017996317370418576);
+ assertComplex(medium, medium, name, operation1, operation2, -1.5656258353157435, 3.297894836311237);
+ assertComplex(medium, small, name, operation1, operation2, 3.7621956910836314, 8.0700322819551687e-308);
+ assertComplex(small, big, name, operation1, operation2, -0.99998768942655991, 1.1040715888508271e-310);
+ assertComplex(small, medium, name, operation1, operation2, -0.41614683654714241, 2.0232539340376892e-308);
+ assertComplex(small, small, name, operation1, operation2, 1, 0);
// Overflow test.
// Based on MATH-901 discussion of FastMath functionality.
@@ -311,34 +315,34 @@ class ComplexEdgeCaseTest {
final double x = 709.783;
Assertions.assertEquals(inf, Math.exp(x));
// As computed by GNU g++
- assertComplex(x, 0, name, operation, 8.9910466927705402e+307, 0.0);
- assertComplex(-x, 0, name, operation, 8.9910466927705402e+307, -0.0);
+ assertComplex(x, 0, name, operation1, operation2, 8.9910466927705402e+307, 0.0);
+ assertComplex(-x, 0, name, operation1, operation2, 8.9910466927705402e+307, -0.0);
// sub-normal number x:
// cos(x) = 1 => real = (e^|x| / 2)
// sin(x) = x => imaginary = x * (e^|x| / 2)
- assertComplex(x, small, name, operation, 8.9910466927705402e+307, 2.0005742956701358);
- assertComplex(-x, small, name, operation, 8.9910466927705402e+307, -2.0005742956701358);
- assertComplex(x, tiny, name, operation, 8.9910466927705402e+307, 4.4421672910524807e-16);
- assertComplex(-x, tiny, name, operation, 8.9910466927705402e+307, -4.4421672910524807e-16);
+ assertComplex(x, small, name, operation1, operation2, 8.9910466927705402e+307, 2.0005742956701358);
+ assertComplex(-x, small, name, operation1, operation2, 8.9910466927705402e+307, -2.0005742956701358);
+ assertComplex(x, tiny, name, operation1, operation2, 8.9910466927705402e+307, 4.4421672910524807e-16);
+ assertComplex(-x, tiny, name, operation1, operation2, 8.9910466927705402e+307, -4.4421672910524807e-16);
// Should not overflow imaginary.
- assertComplex(2 * x, tiny, name, operation, inf, 7.9879467061901743e+292);
- assertComplex(-2 * x, tiny, name, operation, inf, -7.9879467061901743e+292);
+ assertComplex(2 * x, tiny, name, operation1, operation2, inf, 7.9879467061901743e+292);
+ assertComplex(-2 * x, tiny, name, operation1, operation2, inf, -7.9879467061901743e+292);
// Test when large enough to overflow any non-zero value to infinity. Result should be
// as if x was infinite and y was finite.
- assertComplex(3 * x, tiny, name, operation, inf, inf);
- assertComplex(-3 * x, tiny, name, operation, inf, -inf);
+ assertComplex(3 * x, tiny, name, operation1, operation2, inf, inf);
+ assertComplex(-3 * x, tiny, name, operation1, operation2, inf, -inf);
// pi / 2 x:
// cos(x) = ~0 => real = x * (e^|x| / 2)
// sin(x) = ~1 => imaginary = (e^|x| / 2)
final double pi2 = Math.PI / 2;
- assertComplex(x, pi2, name, operation, 5.5054282766429199e+291, 8.9910466927705402e+307);
- assertComplex(-x, pi2, name, operation, 5.5054282766429199e+291, -8.9910466927705402e+307);
- assertComplex(2 * x, pi2, name, operation, inf, inf);
- assertComplex(-2 * x, pi2, name, operation, inf, -inf);
+ assertComplex(x, pi2, name, operation1, operation2, 5.5054282766429199e+291, 8.9910466927705402e+307);
+ assertComplex(-x, pi2, name, operation1, operation2, 5.5054282766429199e+291, -8.9910466927705402e+307);
+ assertComplex(2 * x, pi2, name, operation1, operation2, inf, inf);
+ assertComplex(-2 * x, pi2, name, operation1, operation2, inf, -inf);
// Test when large enough to overflow any non-zero value to infinity. Result should be
// as if x was infinite and y was finite.
- assertComplex(3 * x, pi2, name, operation, inf, inf);
- assertComplex(-3 * x, pi2, name, operation, inf, -inf);
+ assertComplex(3 * x, pi2, name, operation1, operation2, inf, inf);
+ assertComplex(-3 * x, pi2, name, operation1, operation2, inf, -inf);
}
@Test
@@ -346,22 +350,23 @@ class ComplexEdgeCaseTest {
// sinh(a + b i) = sinh(a)cos(b) + i cosh(a)sin(b)
// Odd function: negative real cases defined by positive real cases
final String name = "sinh";
- final UnaryOperator<Complex> operation = Complex::sinh;
+ final UnaryOperator<Complex> operation1 = Complex::sinh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::sinh;
// Implementation defers to java.util.Math.
// Hit edge cases for extreme values.
final double big = Double.MAX_VALUE;
final double medium = 2;
final double small = Double.MIN_NORMAL;
- assertComplex(big, big, name, operation, -inf, inf);
- assertComplex(big, medium, name, operation, -inf, inf);
- assertComplex(big, small, name, operation, inf, inf);
- assertComplex(medium, big, name, operation, -3.6268157591156114, 0.018667844927220067);
- assertComplex(medium, medium, name, operation, -1.5093064853236158, 3.4209548611170133);
- assertComplex(medium, small, name, operation, 3.626860407847019, 8.3711632828186228e-308);
- assertComplex(small, big, name, operation, -2.2250464665720564e-308, 0.004961954789184062);
- assertComplex(small, medium, name, operation, -9.2595744730151568e-309, 0.90929742682568171);
- assertComplex(small, small, name, operation, 2.2250738585072014e-308, 2.2250738585072014e-308);
+ assertComplex(big, big, name, operation1, operation2, -inf, inf);
+ assertComplex(big, medium, name, operation1, operation2, -inf, inf);
+ assertComplex(big, small, name, operation1, operation2, inf, inf);
+ assertComplex(medium, big, name, operation1, operation2, -3.6268157591156114, 0.018667844927220067);
+ assertComplex(medium, medium, name, operation1, operation2, -1.5093064853236158, 3.4209548611170133);
+ assertComplex(medium, small, name, operation1, operation2, 3.626860407847019, 8.3711632828186228e-308);
+ assertComplex(small, big, name, operation1, operation2, -2.2250464665720564e-308, 0.004961954789184062);
+ assertComplex(small, medium, name, operation1, operation2, -9.2595744730151568e-309, 0.90929742682568171);
+ assertComplex(small, small, name, operation1, operation2, 2.2250738585072014e-308, 2.2250738585072014e-308);
// Overflow test.
// As per cosh with sign changes to real and imaginary
@@ -374,34 +379,34 @@ class ComplexEdgeCaseTest {
final double x = 709.783;
Assertions.assertEquals(inf, Math.exp(x));
// As computed by GNU g++
- assertComplex(x, 0, name, operation, 8.9910466927705402e+307, 0.0);
- assertComplex(-x, 0, name, operation, -8.9910466927705402e+307, 0.0);
+ assertComplex(x, 0, name, operation1, operation2, 8.9910466927705402e+307, 0.0);
+ assertComplex(-x, 0, name, operation1, operation2, -8.9910466927705402e+307, 0.0);
// sub-normal number:
// cos(x) = 1 => real = (e^|x| / 2)
// sin(x) = x => imaginary = x * (e^|x| / 2)
- assertComplex(x, small, name, operation, 8.9910466927705402e+307, 2.0005742956701358);
- assertComplex(-x, small, name, operation, -8.9910466927705402e+307, 2.0005742956701358);
- assertComplex(x, tiny, name, operation, 8.9910466927705402e+307, 4.4421672910524807e-16);
- assertComplex(-x, tiny, name, operation, -8.9910466927705402e+307, 4.4421672910524807e-16);
+ assertComplex(x, small, name, operation1, operation2, 8.9910466927705402e+307, 2.0005742956701358);
+ assertComplex(-x, small, name, operation1, operation2, -8.9910466927705402e+307, 2.0005742956701358);
+ assertComplex(x, tiny, name, operation1, operation2, 8.9910466927705402e+307, 4.4421672910524807e-16);
+ assertComplex(-x, tiny, name, operation1, operation2, -8.9910466927705402e+307, 4.4421672910524807e-16);
// Should not overflow imaginary.
- assertComplex(2 * x, tiny, name, operation, inf, 7.9879467061901743e+292);
- assertComplex(-2 * x, tiny, name, operation, -inf, 7.9879467061901743e+292);
+ assertComplex(2 * x, tiny, name, operation1, operation2, inf, 7.9879467061901743e+292);
+ assertComplex(-2 * x, tiny, name, operation1, operation2, -inf, 7.9879467061901743e+292);
// Test when large enough to overflow any non-zero value to infinity. Result should be
// as if x was infinite and y was finite.
- assertComplex(3 * x, tiny, name, operation, inf, inf);
- assertComplex(-3 * x, tiny, name, operation, -inf, inf);
+ assertComplex(3 * x, tiny, name, operation1, operation2, inf, inf);
+ assertComplex(-3 * x, tiny, name, operation1, operation2, -inf, inf);
// pi / 2 x:
// cos(x) = ~0 => real = x * (e^|x| / 2)
// sin(x) = ~1 => imaginary = (e^|x| / 2)
final double pi2 = Math.PI / 2;
- assertComplex(x, pi2, name, operation, 5.5054282766429199e+291, 8.9910466927705402e+307);
- assertComplex(-x, pi2, name, operation, -5.5054282766429199e+291, 8.9910466927705402e+307);
- assertComplex(2 * x, pi2, name, operation, inf, inf);
- assertComplex(-2 * x, pi2, name, operation, -inf, inf);
+ assertComplex(x, pi2, name, operation1, operation2, 5.5054282766429199e+291, 8.9910466927705402e+307);
+ assertComplex(-x, pi2, name, operation1, operation2, -5.5054282766429199e+291, 8.9910466927705402e+307);
+ assertComplex(2 * x, pi2, name, operation1, operation2, inf, inf);
+ assertComplex(-2 * x, pi2, name, operation1, operation2, -inf, inf);
// Test when large enough to overflow any non-zero value to infinity. Result should be
// as if x was infinite and y was finite.
- assertComplex(3 * x, pi2, name, operation, inf, inf);
- assertComplex(-3 * x, pi2, name, operation, -inf, inf);
+ assertComplex(3 * x, pi2, name, operation1, operation2, inf, inf);
+ assertComplex(-3 * x, pi2, name, operation1, operation2, -inf, inf);
}
@Test
@@ -409,34 +414,34 @@ class ComplexEdgeCaseTest {
// tan(a + b i) = sinh(2a)/(cosh(2a)+cos(2b)) + i [sin(2b)/(cosh(2a)+cos(2b))]
// Odd function: negative real cases defined by positive real cases
final String name = "tanh";
- final UnaryOperator<Complex> operation = Complex::tanh;
-
+ final UnaryOperator<Complex> operation1 = Complex::tanh;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::tanh;
// Overflow on 2b:
// cos(2b) = cos(inf) = NaN
// sin(2b) = sin(inf) = NaN
- assertComplex(1, Double.MAX_VALUE, name, operation, 0.76160203106265523, -0.0020838895895863505);
+ assertComplex(1, Double.MAX_VALUE, name, operation1, operation2, 0.76160203106265523, -0.0020838895895863505);
// Underflow on 2b:
// cos(2b) -> 1
// sin(2b) -> 0
- assertComplex(1, Double.MIN_NORMAL, name, operation, 0.76159415595576485, 9.344739287691424e-309);
- assertComplex(1, Double.MIN_VALUE, name, operation, 0.76159415595576485, 0);
+ assertComplex(1, Double.MIN_NORMAL, name, operation1, operation2, 0.76159415595576485, 9.344739287691424e-309);
+ assertComplex(1, Double.MIN_VALUE, name, operation1, operation2, 0.76159415595576485, 0);
// Overflow on 2a:
// sinh(2a) = sinh(inf) = inf
// cosh(2a) = cosh(inf) = inf
// Test all sign variants as this execution path to treat real as infinite
// is not tested else where.
- assertComplex(Double.MAX_VALUE, 1, name, operation, 1, 0.0);
- assertComplex(Double.MAX_VALUE, -1, name, operation, 1, -0.0);
- assertComplex(-Double.MAX_VALUE, 1, name, operation, -1, 0.0);
- assertComplex(-Double.MAX_VALUE, -1, name, operation, -1, -0.0);
+ assertComplex(Double.MAX_VALUE, 1, name, operation1, operation2, 1, 0.0);
+ assertComplex(Double.MAX_VALUE, -1, name, operation1, operation2, 1, -0.0);
+ assertComplex(-Double.MAX_VALUE, 1, name, operation1, operation2, -1, 0.0);
+ assertComplex(-Double.MAX_VALUE, -1, name, operation1, operation2, -1, -0.0);
// Underflow on 2a:
// sinh(2a) -> 0
// cosh(2a) -> 0
- assertComplex(Double.MIN_NORMAL, 1, name, operation, 7.6220323800193346e-308, 1.5574077246549021);
- assertComplex(Double.MIN_VALUE, 1, name, operation, 1.4821969375237396e-323, 1.5574077246549021);
+ assertComplex(Double.MIN_NORMAL, 1, name, operation1, operation2, 7.6220323800193346e-308, 1.5574077246549021);
+ assertComplex(Double.MIN_VALUE, 1, name, operation1, operation2, 1.4821969375237396e-323, 1.5574077246549021);
// Underflow test.
// sinh(x) can be approximated by exp(x) but must be overflow safe.
@@ -448,29 +453,29 @@ class ComplexEdgeCaseTest {
Assertions.assertEquals(1.0, Math.sin(2 * y), 1e-16);
Assertions.assertEquals(Double.POSITIVE_INFINITY, Math.exp(2 * x));
// As computed by GNU g++
- assertComplex(x, y, name, operation, 1, 1.1122175583895849e-308);
+ assertComplex(x, y, name, operation1, operation2, 1, 1.1122175583895849e-308);
}
@Test
void testExp() {
final String name = "exp";
- final UnaryOperator<Complex> operation = Complex::exp;
-
+ final UnaryOperator<Complex> operation1 = Complex::exp;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::exp;
// exp(a + b i) = exp(a) (cos(b) + i sin(b))
// Overflow if exp(a) == inf
- assertComplex(1000, 0, name, operation, inf, 0.0);
- assertComplex(1000, 1, name, operation, inf, inf);
- assertComplex(1000, 2, name, operation, -inf, inf);
- assertComplex(1000, 3, name, operation, -inf, inf);
- assertComplex(1000, 4, name, operation, -inf, -inf);
+ assertComplex(1000, 0, name, operation1, operation2, inf, 0.0);
+ assertComplex(1000, 1, name, operation1, operation2, inf, inf);
+ assertComplex(1000, 2, name, operation1, operation2, -inf, inf);
+ assertComplex(1000, 3, name, operation1, operation2, -inf, inf);
+ assertComplex(1000, 4, name, operation1, operation2, -inf, -inf);
// Underflow if exp(a) == 0
- assertComplex(-1000, 0, name, operation, 0.0, 0.0);
- assertComplex(-1000, 1, name, operation, 0.0, 0.0);
- assertComplex(-1000, 2, name, operation, -0.0, 0.0);
- assertComplex(-1000, 3, name, operation, -0.0, 0.0);
- assertComplex(-1000, 4, name, operation, -0.0, -0.0);
+ assertComplex(-1000, 0, name, operation1, operation2, 0.0, 0.0);
+ assertComplex(-1000, 1, name, operation1, operation2, 0.0, 0.0);
+ assertComplex(-1000, 2, name, operation1, operation2, -0.0, 0.0);
+ assertComplex(-1000, 3, name, operation1, operation2, -0.0, 0.0);
+ assertComplex(-1000, 4, name, operation1, operation2, -0.0, -0.0);
}
@Test
@@ -610,7 +615,8 @@ class ComplexEdgeCaseTest {
@Test
void testSqrt() {
final String name = "sqrt";
- final UnaryOperator<Complex> operation = Complex::sqrt;
+ final UnaryOperator<Complex> operation1 = Complex::sqrt;
+ final ComplexUnaryOperator<ComplexNumber> operation2 = ComplexFunctions::sqrt;
// Test real/imaginary only numbers satisfy the definition using polar coordinates:
// real = sqrt(abs()) * Math.cos(arg() / 2)
@@ -645,14 +651,14 @@ class ComplexEdgeCaseTest {
Assertions.assertEquals(root2over2, sinArgIm, ulp, "Expected sin(pi/4) to be 1 ulp from sqrt(2) / 2");
for (final double a : new double[] {0.5, 1.0, 1.2322, 345345.234523}) {
final double rootA = Math.sqrt(a);
- assertComplex(a, 0, name, operation, rootA * cosArgRe, rootA * sinArgRe, 0);
+ assertComplex(a, 0, name, operation1, operation2, rootA * cosArgRe, rootA * sinArgRe, 0);
// This should be exact. It will fail if using the polar computation
// real = sqrt(abs()) * Math.cos(arg() / 2) as cos(pi/2) is not 0.0 but 6.123233995736766e-17
- assertComplex(-a, 0, name, operation, rootA * sinArgRe, rootA * cosArgRe, 0);
+ assertComplex(-a, 0, name, operation1, operation2, rootA * sinArgRe, rootA * cosArgRe, 0);
// This should be exact. It won't be if Complex is using polar computation
// with sin/cos which does not output the same result for angle pi/4.
- assertComplex(0, a, name, operation, rootA * root2over2, rootA * root2over2, 0);
- assertComplex(0, -a, name, operation, rootA * root2over2, -rootA * root2over2, 0);
+ assertComplex(0, a, name, operation1, operation2, rootA * root2over2, rootA * root2over2, 0);
+ assertComplex(0, -a, name, operation1, operation2, rootA * root2over2, -rootA * root2over2, 0);
}
// Check overflow safe.
@@ -664,26 +670,26 @@ class ComplexEdgeCaseTest {
// new BigDecimal(b).multiply(new BigDecimal(b)))
// .sqrt(MathContext.DECIMAL128).sqrt(MathContext.DECIMAL128).doubleValue()
final double newAbs = 1.3612566508088272E154;
- assertComplex(a, b, name, operation, newAbs * Math.cos(0.5 * Math.atan2(b, a)),
- newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3);
- assertComplex(b, a, name, operation, newAbs * Math.cos(0.5 * Math.atan2(a, b)),
- newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2);
+ assertComplex(a, b, name, operation1, operation2, newAbs * Math.cos(0.5 * Math.atan2(b, a)),
+ newAbs * Math.sin(0.5 * Math.atan2(b, a)), 3);
+ assertComplex(b, a, name, operation1, operation2, newAbs * Math.cos(0.5 * Math.atan2(a, b)),
+ newAbs * Math.sin(0.5 * Math.atan2(a, b)), 2);
// Note that the computation is possible in polar coords if abs() does not overflow.
a = Double.MAX_VALUE / 2;
- assertComplex(-a, a, name, operation, 4.3145940638864765e+153, 1.0416351505169177e+154, 2);
- assertComplex(a, a, name, operation, 1.0416351505169177e+154, 4.3145940638864758e+153);
- assertComplex(-a, -a, name, operation, 4.3145940638864765e+153, -1.0416351505169177e+154, 2);
- assertComplex(a, -a, name, operation, 1.0416351505169177e+154, -4.3145940638864758e+153);
+ assertComplex(-a, a, name, operation1, operation2, 4.3145940638864765e+153, 1.0416351505169177e+154, 2);
+ assertComplex(a, a, name, operation1, operation2, 1.0416351505169177e+154, 4.3145940638864758e+153);
+ assertComplex(-a, -a, name, operation1, operation2, 4.3145940638864765e+153, -1.0416351505169177e+154, 2);
+ assertComplex(a, -a, name, operation1, operation2, 1.0416351505169177e+154, -4.3145940638864758e+153);
// Check minimum normal value conditions
// Computing in polar coords produces a very different result with
// MIN_VALUE so use MIN_NORMAL
a = Double.MIN_NORMAL;
- assertComplex(-a, a, name, operation, 6.7884304867749663e-155, 1.6388720948399111e-154);
- assertComplex(a, a, name, operation, 1.6388720948399111e-154, 6.7884304867749655e-155);
- assertComplex(-a, -a, name, operation, 6.7884304867749663e-155, -1.6388720948399111e-154);
- assertComplex(a, -a, name, operation, 1.6388720948399111e-154, -6.7884304867749655e-155);
+ assertComplex(-a, a, name, operation1, operation2, 6.7884304867749663e-155, 1.6388720948399111e-154);
+ assertComplex(a, a, name, operation1, operation2, 1.6388720948399111e-154, 6.7884304867749655e-155);
+ assertComplex(-a, -a, name, operation1, operation2, 6.7884304867749663e-155, -1.6388720948399111e-154);
+ assertComplex(a, -a, name, operation1, operation2, 1.6388720948399111e-154, -6.7884304867749655e-155);
}
// Note: inf/nan edge cases for
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java
index 657d2d75..35ebd50b 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java
@@ -294,25 +294,31 @@ class ComplexTest {
@Test
void testAbs() {
- final Complex z = Complex.ofCartesian(3.0, 4.0);
- Assertions.assertEquals(5.0, z.abs());
+ final Complex input = Complex.ofCartesian(3.0, 4.0);
+ double z = TestUtils.assertSame(input, "abs", Complex::abs, ComplexFunctions::abs);
+ Assertions.assertEquals(5.0, z);
}
@Test
void testAbsNaN() {
// The result is NaN if either argument is NaN and the other is not infinite
- Assertions.assertEquals(nan, NAN.abs());
- Assertions.assertEquals(nan, Complex.ofCartesian(3.0, nan).abs());
- Assertions.assertEquals(nan, Complex.ofCartesian(nan, 3.0).abs());
+ assertAbs(nan, NAN);
+ assertAbs(nan, Complex.ofCartesian(3.0, nan));
+ assertAbs(nan, Complex.ofCartesian(nan, 3.0));
// The result is positive infinite if either argument is infinite
- Assertions.assertEquals(inf, Complex.ofCartesian(inf, nan).abs());
- Assertions.assertEquals(inf, Complex.ofCartesian(-inf, nan).abs());
- Assertions.assertEquals(inf, Complex.ofCartesian(nan, inf).abs());
- Assertions.assertEquals(inf, Complex.ofCartesian(nan, -inf).abs());
- Assertions.assertEquals(inf, Complex.ofCartesian(inf, 3.0).abs());
- Assertions.assertEquals(inf, Complex.ofCartesian(-inf, 3.0).abs());
- Assertions.assertEquals(inf, Complex.ofCartesian(3.0, inf).abs());
- Assertions.assertEquals(inf, Complex.ofCartesian(3.0, -inf).abs());
+ assertAbs(inf, Complex.ofCartesian(inf, nan));
+ assertAbs(inf, Complex.ofCartesian(-inf, nan));
+ assertAbs(inf, Complex.ofCartesian(nan, inf));
+ assertAbs(inf, Complex.ofCartesian(nan, -inf));
+ assertAbs(inf, Complex.ofCartesian(inf, 3.0));
+ assertAbs(inf, Complex.ofCartesian(-inf, 3.0));
+ assertAbs(inf, Complex.ofCartesian(3.0, inf));
+ assertAbs(inf, Complex.ofCartesian(3.0, -inf));
+ }
+
+ private static void assertAbs(double expected, Complex input) {
+ final double actual = TestUtils.assertSame(input, "abs", Complex::abs, ComplexFunctions::abs);
+ Assertions.assertEquals(expected, actual);
}
/**
@@ -371,28 +377,33 @@ class ComplexTest {
}
private static void assertArgument(double expected, Complex complex, double delta) {
- final double actual = complex.arg();
+ final double actual = TestUtils.assertSame(complex, "arg", Complex::arg, ComplexFunctions::arg);
Assertions.assertEquals(expected, actual, delta);
- Assertions.assertEquals(actual, complex.arg(), delta);
}
@Test
void testNorm() {
final Complex z = Complex.ofCartesian(3.0, 4.0);
- Assertions.assertEquals(25.0, z.norm());
+ double norm = TestUtils.assertSame(z, "norm", Complex::norm, ComplexFunctions::norm);
+ Assertions.assertEquals(25.0, norm);
}
@Test
void testNormNaN() {
// The result is NaN if either argument is NaN and the other is not infinite
- Assertions.assertEquals(nan, NAN.norm());
- Assertions.assertEquals(nan, Complex.ofCartesian(3.0, nan).norm());
- Assertions.assertEquals(nan, Complex.ofCartesian(nan, 3.0).norm());
+ assertNorm(nan, NAN);
+ assertNorm(nan, Complex.ofCartesian(3.0, nan));
+ assertNorm(nan, Complex.ofCartesian(nan, 3.0));
// The result is positive infinite if either argument is infinite
- Assertions.assertEquals(inf, Complex.ofCartesian(inf, nan).norm());
- Assertions.assertEquals(inf, Complex.ofCartesian(-inf, nan).norm());
- Assertions.assertEquals(inf, Complex.ofCartesian(nan, inf).norm());
- Assertions.assertEquals(inf, Complex.ofCartesian(nan, -inf).norm());
+ assertNorm(inf, Complex.ofCartesian(inf, nan));
+ assertNorm(inf, Complex.ofCartesian(-inf, nan));
+ assertNorm(inf, Complex.ofCartesian(nan, inf));
+ assertNorm(inf, Complex.ofCartesian(nan, -inf));
+ }
+
+ private static void assertNorm(double expected, Complex input) {
+ final double actual = TestUtils.assertSame(input, "norm", Complex::norm, ComplexFunctions::norm);
+ Assertions.assertEquals(expected, actual);
}
/**
@@ -431,9 +442,9 @@ class ComplexTest {
*/
private static void assertNumberType(double real, double imaginary, NumberType type) {
final Complex z = Complex.ofCartesian(real, imaginary);
- final boolean isNaN = z.isNaN();
- final boolean isInfinite = z.isInfinite();
- final boolean isFinite = z.isFinite();
+ final boolean isNaN = TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN);
+ final boolean isInfinite = TestUtils.assertSame(z, "isInfinite", Complex::isInfinite, ComplexFunctions::isInfinite);
+ final boolean isFinite = TestUtils.assertSame(z, "isFinite", Complex::isFinite, ComplexFunctions::isFinite);
// A number can be only one
int count = isNaN ? 1 : 0;
count += isInfinite ? 1 : 0;
@@ -459,53 +470,61 @@ class ComplexTest {
@Test
void testConjugate() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
- final Complex z = x.conj();
+ final Complex z = TestUtils.assertSame(x, "conj", Complex::conj, ComplexFunctions::conj);
Assertions.assertEquals(3.0, z.getReal());
Assertions.assertEquals(-4.0, z.getImaginary());
}
@Test
void testConjugateNaN() {
- final Complex z = NAN.conj();
- Assertions.assertTrue(z.isNaN());
+ final Complex z = TestUtils.assertSame(NAN, "conj", Complex::conj, ComplexFunctions::conj);
+ Assertions.assertTrue(TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN));
}
@Test
void testConjugateInfinite() {
Complex z = Complex.ofCartesian(0, inf);
- Assertions.assertEquals(neginf, z.conj().getImaginary());
+ Complex zConj = TestUtils.assertSame(z, "conj", Complex::conj, ComplexFunctions::conj);
+ Assertions.assertEquals(neginf, zConj.getImaginary());
+
z = Complex.ofCartesian(0, neginf);
- Assertions.assertEquals(inf, z.conj().getImaginary());
+ zConj = TestUtils.assertSame(z, "conj", Complex::conj, ComplexFunctions::conj);
+ Assertions.assertEquals(inf, zConj.getImaginary());
}
@Test
void testNegate() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
- final Complex z = x.negate();
+ final Complex z = TestUtils.assertSame(x, "negate", Complex::negate, ComplexFunctions::negate);
Assertions.assertEquals(-3.0, z.getReal());
Assertions.assertEquals(-4.0, z.getImaginary());
}
@Test
void testNegateNaN() {
- final Complex z = NAN.negate();
- Assertions.assertTrue(z.isNaN());
+ final Complex z = TestUtils.assertSame(NAN, "negate", Complex::negate, ComplexFunctions::negate);
+ Assertions.assertTrue(TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN));
}
@Test
void testProj() {
final Complex z = Complex.ofCartesian(3.0, 4.0);
- Assertions.assertSame(z, z.proj());
+ assertProj(z, z);
// Sign must be the same for projection
- TestUtils.assertSame(infZero, Complex.ofCartesian(inf, 4.0).proj());
- TestUtils.assertSame(infZero, Complex.ofCartesian(inf, inf).proj());
- TestUtils.assertSame(infZero, Complex.ofCartesian(inf, nan).proj());
- TestUtils.assertSame(infZero, Complex.ofCartesian(3.0, inf).proj());
- TestUtils.assertSame(infZero, Complex.ofCartesian(nan, inf).proj());
- TestUtils.assertSame(infNegZero, Complex.ofCartesian(inf, -4.0).proj());
- TestUtils.assertSame(infNegZero, Complex.ofCartesian(inf, -inf).proj());
- TestUtils.assertSame(infNegZero, Complex.ofCartesian(3.0, -inf).proj());
- TestUtils.assertSame(infNegZero, Complex.ofCartesian(nan, -inf).proj());
+ assertProj(infZero, Complex.ofCartesian(inf, 4.0));
+ assertProj(infZero, Complex.ofCartesian(inf, inf));
+ assertProj(infZero, Complex.ofCartesian(inf, nan));
+ assertProj(infZero, Complex.ofCartesian(3.0, inf));
+ assertProj(infZero, Complex.ofCartesian(nan, inf));
+ assertProj(infNegZero, Complex.ofCartesian(inf, -4.0));
+ assertProj(infNegZero, Complex.ofCartesian(inf, -inf));
+ assertProj(infNegZero, Complex.ofCartesian(3.0, -inf));
+ assertProj(infNegZero, Complex.ofCartesian(nan, -inf));
+ }
+
+ private static void assertProj(Complex expected, Complex input) {
+ final Complex actual = TestUtils.assertSame(input, "proj", Complex::proj, ComplexFunctions::proj);
+ TestUtils.assertSame(expected, actual);
}
@Test
@@ -834,7 +853,7 @@ class ComplexTest {
void testMultiplyInfInf() {
final Complex z = infInf.multiply(infInf);
// Assert.assertTrue(z.isNaN()); // MATH-620
- Assertions.assertTrue(z.isInfinite());
+ Assertions.assertTrue(TestUtils.assertSame(z, "isInfinite", Complex::isInfinite, ComplexFunctions::isInfinite));
// Expected results from g++:
Assertions.assertEquals(Complex.ofCartesian(nan, inf), infInf.multiply(infInf));
@@ -1098,7 +1117,7 @@ class ComplexTest {
void testDivideNaN() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
final Complex z = x.divide(NAN);
- Assertions.assertTrue(z.isNaN());
+ Assertions.assertTrue(TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN));
}
@Test
@@ -1510,7 +1529,8 @@ class ComplexTest {
theta += pi / 12;
final Complex z = Complex.ofPolar(r, theta);
final Complex sqrtz = Complex.ofPolar(Math.sqrt(r), theta / 2);
- TestUtils.assertEquals(sqrtz, z.sqrt(), tol);
+ Complex actual = TestUtils.assertSame(z, "sqrt", Complex::sqrt, ComplexFunctions::sqrt);
+ TestUtils.assertEquals(sqrtz, actual, tol);
}
}
}
@@ -1952,7 +1972,8 @@ class ComplexTest {
@Test
void testAtanhEdgeConditions() {
// Hits the edge case when imaginary == 0 but real != 0 or 1
- final Complex c = Complex.ofCartesian(2, 0).atanh();
+ final Complex z = Complex.ofCartesian(2, 0);
+ final Complex c = TestUtils.assertSame(z, "atanh", Complex::atanh, ComplexFunctions::atanh);
// Answer from g++
Assertions.assertEquals(0.54930614433405489, c.getReal());
Assertions.assertEquals(1.5707963267948966, c.getImaginary());
@@ -2092,7 +2113,7 @@ class ComplexTest {
// is consistent.
for (int i = 0; i < samples; i++) {
final Complex z = supplier.get();
- final double abs = z.abs();
+ final double abs = TestUtils.assertSame(z, "abs", Complex::abs, ComplexFunctions::abs);
final double x = Math.abs(z.getReal());
final double y = Math.abs(z.getImaginary());
@@ -2105,7 +2126,7 @@ class ComplexTest {
// argument:
// real = sqrt(|z|) * cos(0.5 * arg(z))
// imag = sqrt(|z|) * sin(0.5 * arg(z))
- final Complex c = z.sqrt();
+ final Complex c = TestUtils.assertSame(z, "sqrt", Complex::sqrt, ComplexFunctions::sqrt);
final double t = Math.sqrt(2 * (x + abs));
if (z.getReal() >= 0) {
Assertions.assertEquals(t / 2, c.getReal());
@@ -2158,13 +2179,13 @@ class ComplexTest {
// is consistent.
for (int i = 0; i < samples; i++) {
final Complex z = supplier.get();
- final double abs = z.abs();
+ final double abs = TestUtils.assertSame(z, "abs", Complex::abs, ComplexFunctions::abs);
final double x = Math.abs(z.getReal());
final double y = Math.abs(z.getImaginary());
// log(x + iy) = log(|x + i y|) + i arg(x + i y)
// Only test the real component
- final Complex c = z.log();
+ final Complex c = TestUtils.assertSame(z, "log", Complex::log, ComplexFunctions::log);
Assertions.assertEquals(Math.log(abs), c.getReal());
}
}
diff --git a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/TestUtils.java b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/TestUtils.java
index e3d86a60..a601df6c 100644
--- a/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/TestUtils.java
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/TestUtils.java
@@ -27,6 +27,9 @@ import java.io.ObjectOutputStream;
import java.util.ArrayList;
import java.util.List;
import java.util.function.Consumer;
+import java.util.function.DoubleBinaryOperator;
+import java.util.function.Predicate;
+import java.util.function.ToDoubleFunction;
import java.util.function.UnaryOperator;
import org.apache.commons.numbers.core.Precision;
@@ -38,6 +41,24 @@ import org.junit.jupiter.api.Assertions;
*/
public final class TestUtils {
+ /**
+ * Represents a predicate (boolean-valued function) of two {@code double}-valued
+ * arguments. This is the {@code double}-consuming primitive type specialization
+ * of {@link java.util.function.BiPredicate BiPredicate}.
+ */
+ @FunctionalInterface
+ public interface DoubleBinaryPredicate {
+ /**
+ * Evaluates this predicate on the given argument.
+ *
+ * @param a the first input argument
+ * @param b the second input argument
+ * @return {@code true} if the input arguments match the predicate,
+ * otherwise {@code false}
+ */
+ boolean test(double a, double b);
+ }
+
/**
* The option for how to process test data lines flagged (prefixed)
* with the {@code ;} character
@@ -410,4 +431,46 @@ public final class TestUtils {
Assertions.assertEquals(z.imag(), z2.getImaginary(), () -> "Unary operator mismatch: " + name + " imaginary");
return z;
}
+
+ /**
+ * Assert the operation on the complex number is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts.
+ *
+ * @param c Input complex number.
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex object.
+ * @param operation2 Operation on the complex real and imaginary parts.
+ * @return Result double from the given operation.
+ */
+ public static double assertSame(Complex c,
+ String name,
+ ToDoubleFunction<Complex> operation1,
+ DoubleBinaryOperator operation2) {
+ final double r1 = operation1.applyAsDouble(c);
+ // Test operation2 produces the exact same result
+ final double r2 = operation2.applyAsDouble(c.real(), c.imag());
+ Assertions.assertEquals(r1, r2, () -> "Double operator mismatch: " + name);
+ return r1;
+ }
+
+ /**
+ * Assert the condition operation on the complex number is <em>exactly</em> equal to the condition
+ * operation on complex real and imaginary parts.
+ *
+ * @param c Input complex number.
+ * @param name Operation name.
+ * @param operation1 Condition operation on the Complex object.
+ * @param operation2 Condition peration on the complex real and imaginary parts.
+ * @return Result boolean from the given operation.
+ */
+ public static boolean assertSame(Complex c,
+ String name,
+ Predicate<Complex> operation1,
+ DoubleBinaryPredicate operation2) {
+ final boolean b1 = operation1.test(c);
+ // Test operation2 produces the exact same result
+ final boolean b2 = operation2.test(c.real(), c.imag());
+ Assertions.assertEquals(b1, b2, () -> "Predicate operator mismatch: " + name);
+ return b1;
+ }
}
diff --git a/src/main/resources/checkstyle/checkstyle-suppressions.xml b/src/main/resources/checkstyle/checkstyle-suppressions.xml
index 96e82e5f..f0c7931e 100644
--- a/src/main/resources/checkstyle/checkstyle-suppressions.xml
+++ b/src/main/resources/checkstyle/checkstyle-suppressions.xml
@@ -22,6 +22,7 @@
<suppress checks="Indentation" files=".*[/\\]combinatorics[/\\]Factorial\.java" />
<suppress checks="ParameterNumber" files=".*[/\\]core[/\\]LinearCombination[s]?\.java" />
<suppress checks="FileLengthCheck" files=".*[/\\]Complex(Test)?\.java" />
+ <suppress checks="FileLengthCheck" files=".*[/\\]ComplexFunctions?\.java" />
<suppress checks="FileLengthCheck" files=".*jmh[/\\]core[/\\]LinearCombinations.java" />
<suppress checks="MethodLength" files=".*[/\\]Complex\.java" />
<!-- Naming and method logic preserved from the original Boost C++ code -->