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:49 UTC
[commons-numbers] branch complex-gsoc-2022 updated (5939b5ea -> 1b0bfa60)
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a change to branch complex-gsoc-2022
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
from 5939b5ea Remove redundant method
new ab614de1 NUMBERS-188: Refactor complex unuary functions to static methods
new e4c7642f Remove redundant javadoc in ComplexFunctions
new 1b0bfa60 Test code clean-up
The 3 revisions listed above as "new" are entirely new to this
repository and will be described in separate emails. The revisions
listed as "add" were already present in the repository and have only
been added to this reference.
Summary of changes:
.../apache/commons/numbers/complex/Complex.java | 1159 +----------
.../commons/numbers/complex/ComplexFunctions.java | 2111 ++++++++++++++++++--
.../commons/numbers/complex/CReferenceTest.java | 51 +-
.../commons/numbers/complex/CStandardTest.java | 863 ++++----
.../numbers/complex/ComplexEdgeCaseTest.java | 321 ++-
.../commons/numbers/complex/ComplexTest.java | 125 +-
.../apache/commons/numbers/complex/TestUtils.java | 65 +-
.../checkstyle/checkstyle-suppressions.xml | 1 +
8 files changed, 2639 insertions(+), 2057 deletions(-)
[commons-numbers] 02/03: Remove redundant javadoc in ComplexFunctions
Posted by ah...@apache.org.
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 e4c7642f87183005b4ac80f41f0b1a0cf55821a0
Author: aherbert <ah...@apache.org>
AuthorDate: Mon Jul 25 10:51:40 2022 +0100
Remove redundant javadoc in ComplexFunctions
Drop some internal compute methods where the additional javadoc is
redundant.
Minor formatting changes.
---
.../commons/numbers/complex/ComplexFunctions.java | 120 +++------------------
1 file changed, 12 insertions(+), 108 deletions(-)
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 e4d783d1..460e1f39 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
@@ -217,31 +217,6 @@ public final class ComplexFunctions {
* @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>
- *
- * <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 \( 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);
@@ -438,7 +413,7 @@ public final class ComplexFunctions {
* <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 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).
@@ -459,6 +434,7 @@ public final class ComplexFunctions {
* @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>
+ * @see Complex#ofCis(double)
*/
public static <R> R exp(double real, double imaginary, ComplexSink<R> action) {
if (Double.isInfinite(real)) {
@@ -753,20 +729,6 @@ public final class ComplexFunctions {
* @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
@@ -1003,9 +965,6 @@ public final class ComplexFunctions {
/**
* 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:
@@ -1022,7 +981,7 @@ public final class ComplexFunctions {
* @return the object returned by the supplied action.
*/
private static <R> R computeAsin(final double real, final double imaginary,
- ComplexSink<R> action) {
+ 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);
@@ -1201,9 +1160,6 @@ public final class ComplexFunctions {
/**
* 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:
@@ -1220,7 +1176,7 @@ public final class ComplexFunctions {
* @return the object returned by the supplied action.
*/
private static <R> R computeAcos(final double real, final double imaginary,
- final ComplexSink<R> action) {
+ 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);
@@ -1367,7 +1323,7 @@ public final class ComplexFunctions {
/**
* Returns the
* <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
- * hyperbolic sine</a> of the complexnumber.
+ * hyperbolic sine</a> of the complex number.
*
* <p>\[ \sinh(z) = \frac{1}{2} \left( e^{z} - e^{-z} \right) \]
*
@@ -1384,7 +1340,7 @@ public final class ComplexFunctions {
* <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 +∞ + iy for positive finite y, returns +∞ cis(y).
* <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.
@@ -1402,25 +1358,9 @@ public final class ComplexFunctions {
* @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>
+ * @see Complex#ofCis(double)
*/
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);
}
@@ -1454,7 +1394,7 @@ public final class ComplexFunctions {
/**
* Returns the
* <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
- * hyperbolic cosine</a> of the complexnumber.
+ * hyperbolic cosine</a> of the complex number.
*
* <p>\[ \cosh(z) = \frac{1}{2} \left( e^{z} + e^{-z} \right) \]
*
@@ -1471,7 +1411,7 @@ public final class ComplexFunctions {
* <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 +∞ + iy for finite nonzero y, returns +∞ cis(y).
* <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).
@@ -1489,25 +1429,9 @@ public final class ComplexFunctions {
* @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>
+ * @see Complex#ofCis(double)
*/
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)) {
@@ -1608,7 +1532,7 @@ public final class ComplexFunctions {
/**
* Returns the
* <a href="http://mathworld.wolfram.com/HyperbolicTangent.html">
- * hyperbolic tangent</a> of the complexnumber.
+ * hyperbolic tangent</a> of the complex number.
*
* <p>\[ \tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}} \]
*
@@ -1654,23 +1578,6 @@ public final class ComplexFunctions {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
*/
public static <R> R tanh(double real, double imaginary, ComplexSink<R> action) {
- return computeTanh(real, imaginary, action);
- }
-
- /**
- * Returns the hyperbolic tangent of the complex number using its real
- * and imaginary parts.
- *
- * <p>This function exists to allow implementation of the identity
- * {@code 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 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 <R> R computeTanh(double real, double imaginary, ComplexSink<R> action) {
// Cache the absolute real value
final double x = Math.abs(real);
@@ -1940,9 +1847,6 @@ public final class ComplexFunctions {
* Returns the inverse hyperbolic tangent of the complex number using its
* real and imaginary parts.
*
- * <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:
@@ -1959,7 +1863,7 @@ public final class ComplexFunctions {
* @return the object returned by the supplied action.
*/
private static <R> R computeAtanh(final double real, final double imaginary,
- final ComplexSink<R> action) {
+ final ComplexSink<R> action) {
// Compute with positive values and determine sign at the end
double x = Math.abs(real);
double y = Math.abs(imaginary);
[commons-numbers] 01/03: NUMBERS-188: Refactor complex unuary functions to static methods
Posted by ah...@apache.org.
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 -->
[commons-numbers] 03/03: Test code clean-up
Posted by ah...@apache.org.
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 1b0bfa60cd8ae023fc75eac7f91e97986a887fa6
Author: aherbert <ah...@apache.org>
AuthorDate: Mon Jul 25 10:59:35 2022 +0100
Test code clean-up
Remove unused assert methods (replace by assertions with
UnaryOperator<Complex> operation1 and
ComplexUnaryOperator<ComplexNumber> operation2).
Use final.
Correct typos.
---
.../commons/numbers/complex/CStandardTest.java | 106 ++++++++++-----------
.../numbers/complex/ComplexEdgeCaseTest.java | 45 +--------
.../commons/numbers/complex/ComplexTest.java | 6 +-
.../apache/commons/numbers/complex/TestUtils.java | 8 +-
4 files changed, 62 insertions(+), 103 deletions(-)
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 2ac487b0..85bc5590 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
@@ -658,8 +658,8 @@ class CStandardTest {
ToDoubleFunction<UniformRandomProvider> fy,
int samples) {
for (int i = 0; i < samples; i++) {
- double x = fx.applyAsDouble(rng);
- double y = fy.applyAsDouble(rng);
+ final double x = fx.applyAsDouble(rng);
+ final double y = fy.applyAsDouble(rng);
assertAbs(Complex.ofCartesian(x, y), 1);
}
}
@@ -1084,7 +1084,7 @@ class CStandardTest {
@Override
public double applyAsDouble(UniformRandomProvider rng) {
if (Double.isNaN(tmp)) {
- double u = rng.nextDouble() * Math.PI;
+ final double u = rng.nextDouble() * Math.PI;
tmp = Math.cos(u);
return Math.sin(u);
}
@@ -1109,23 +1109,23 @@ class CStandardTest {
assertComplex(negZeroZero, name, operation1, operation2, piTwoNegZero);
assertComplex(zeroNaN, name, operation1, operation2, piTwoNaN);
assertComplex(negZeroNaN, name, operation1, operation2, piTwoNaN);
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, inf), name, operation1, operation2, piTwoNegInf);
}
- for (double x : nonZeroFinite) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(-inf, y), name, operation1, operation2, piNegInf);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, zeroNegInf);
}
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) {
+ for (final double y : finite) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
assertComplex(nanInf, name, operation1, operation2, nanNegInf);
@@ -1147,25 +1147,25 @@ class CStandardTest {
assertConjugateEquality(name, operation1, operation2);
assertComplex(Complex.ZERO, name, operation1, operation2, zeroPiTwo);
assertComplex(negZeroZero, name, operation1, operation2, zeroPiTwo);
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, inf), name, operation1, operation2, infPiTwo);
}
assertComplex(zeroNaN, name, operation1, operation2, nanPiTwo);
assertComplex(negZeroNaN, name, operation1, operation2, nanPiTwo);
- for (double x : nonZeroFinite) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(-inf, y), name, operation1, operation2, infPi);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, infZero);
}
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) {
+ for (final double y : finite) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
assertComplex(nanInf, name, operation1, operation2, infNaN);
@@ -1185,19 +1185,19 @@ class CStandardTest {
assertConjugateEquality(name, operation1, operation2);
assertFunctionType(name, operation1, operation2, type);
assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ZERO, type);
- for (double x : positiveFinite) {
+ for (final double x : positiveFinite) {
assertComplex(complex(x, inf), name, operation1, operation2, infPiTwo, type);
}
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, infZero, 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) {
+ for (final double y : nonZeroFinite) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
assertComplex(nanInf, name, operation1, operation2, infNaN, type, UnspecifiedSign.REAL);
@@ -1219,19 +1219,19 @@ class CStandardTest {
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) {
+ for (final double x : positiveFinite) {
assertComplex(complex(x, inf), name, operation1, operation2, zeroPiTwo, type);
}
- for (double x : nonZeroFinite) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, zeroPiTwo, 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) {
+ for (final double y : finite) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
assertComplex(nanInf, name, operation1, operation2, zeroPiTwo, type, UnspecifiedSign.REAL);
@@ -1253,20 +1253,20 @@ class CStandardTest {
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) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, inf), name, operation1, operation2, NAN, type);
}
- for (double x : nonZeroFinite) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
assertComplex(infZero, name, operation1, operation2, infZero, type);
- for (double y : nonZeroFinite) {
+ for (final double y : nonZeroFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(inf), type);
}
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) {
+ for (final double y : nonZero) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
assertComplex(NAN, name, operation1, operation2, NAN, type);
@@ -1287,21 +1287,21 @@ class CStandardTest {
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) {
+ for (final double x : nonZeroPositiveFinite) {
assertComplex(complex(x, inf), name, operation1, operation2, NAN, type);
}
- for (double x : nonZeroFinite) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN, 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) {
+ for (final double y : nonZeroFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(inf), 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) {
+ for (final double y : nonZero) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
assertComplex(NAN, name, operation1, operation2, NAN, type);
@@ -1324,20 +1324,20 @@ class CStandardTest {
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) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, inf), name, operation1, operation2, NAN, type);
}
assertComplex(zeroNaN, name, operation1, operation2, zeroNaN, type);
- for (double x : nonZeroFinite) {
+ for (final double x : nonZeroFinite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN, type);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, complex(1.0, Math.copySign(0, Math.sin(2 * y))), 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) {
+ for (final double y : nonZero) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN, type);
}
assertComplex(NAN, name, operation1, operation2, NAN, type);
@@ -1355,17 +1355,17 @@ class CStandardTest {
assertConjugateEquality(name, operation1, operation2);
assertComplex(Complex.ZERO, name, operation1, operation2, oneZero);
assertComplex(negZeroZero, name, operation1, operation2, oneZero);
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, inf), name, operation1, operation2, NAN);
}
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
assertComplex(infZero, name, operation1, operation2, infZero);
- for (double y : finite) {
+ for (final double y : finite) {
assertComplex(complex(-inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(0.0));
}
- for (double y : nonZeroFinite) {
+ for (final double y : nonZeroFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, Complex.ofCis(y).multiply(inf));
}
assertComplex(negInfInf, name, operation1, operation2, Complex.ZERO, UnspecifiedSign.REAL_IMAGINARY);
@@ -1373,7 +1373,7 @@ class CStandardTest {
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) {
+ for (final double y : nonZero) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
assertComplex(NAN, name, operation1, operation2, NAN);
@@ -1391,23 +1391,23 @@ class CStandardTest {
assertConjugateEquality(name, operation1, operation2);
assertComplex(negZeroZero, name, operation1, operation2, negInfPi);
assertComplex(Complex.ZERO, name, operation1, operation2, negInfZero);
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, inf), name, operation1, operation2, infPiTwo);
}
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(-inf, y), name, operation1, operation2, infPi);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, infZero);
}
assertComplex(negInfInf, name, operation1, operation2, infThreePiFour);
assertComplex(infInf, name, operation1, operation2, infPiFour);
assertComplex(negInfNaN, name, operation1, operation2, infNaN);
assertComplex(infNaN, name, operation1, operation2, infNaN);
- for (double y : finite) {
+ for (final double y : finite) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
assertComplex(nanInf, name, operation1, operation2, infNaN);
@@ -1427,23 +1427,23 @@ class CStandardTest {
assertConjugateEquality(name, operation1, operation2);
assertComplex(negZeroZero, name, operation1, operation2, negInfPi);
assertComplex(Complex.ZERO, name, operation1, operation2, negInfZero);
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, inf), name, operation1, operation2, infPiTwo);
}
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(-inf, y), name, operation1, operation2, infPi);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, infZero);
}
assertComplex(negInfInf, name, operation1, operation2, infThreePiFour);
assertComplex(infInf, name, operation1, operation2, infPiFour);
assertComplex(negInfNaN, name, operation1, operation2, infNaN);
assertComplex(infNaN, name, operation1, operation2, infNaN);
- for (double y : finite) {
+ for (final double y : finite) {
assertComplex(complex(nan, y), name, operation1, operation2, NAN);
}
assertComplex(nanInf, name, operation1, operation2, infNaN);
@@ -1462,25 +1462,25 @@ class CStandardTest {
assertConjugateEquality(name, operation1, operation2);
assertComplex(negZeroZero, name, operation1, operation2, Complex.ZERO);
assertComplex(Complex.ZERO, name, operation1, operation2, Complex.ZERO);
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, inf), name, operation1, operation2, infInf);
}
// Include infinity and nan for (x, inf).
assertComplex(infInf, name, operation1, operation2, infInf);
assertComplex(negInfInf, name, operation1, operation2, infInf);
assertComplex(nanInf, name, operation1, operation2, infInf);
- for (double x : finite) {
+ for (final double x : finite) {
assertComplex(complex(x, nan), name, operation1, operation2, NAN);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(-inf, y), name, operation1, operation2, zeroInf);
}
- for (double y : positiveFinite) {
+ for (final double y : positiveFinite) {
assertComplex(complex(inf, y), name, operation1, operation2, infZero);
}
assertComplex(negInfNaN, name, operation1, operation2, nanInf, UnspecifiedSign.IMAGINARY);
assertComplex(infNaN, name, operation1, operation2, infNaN);
- for (double y : finite) {
+ for (final double y : finite) {
assertComplex(complex(nan, y), name, operation1, operation2, 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 dfd53f45..72cc88d3 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
@@ -40,47 +40,6 @@ class ComplexEdgeCaseTest {
private static final double inf = Double.POSITIVE_INFINITY;
private static final double nan = Double.NaN;
- /**
- * Assert the operation on the complex number is equal to the expected value.
- *
- * <p>The results are are considered equal if there are no floating-point values between them.
- *
- * @param a Real part.
- * @param b Imaginary part.
- * @param name The operation name.
- * @param operation The operation.
- * @param x Expected real part.
- * @param y Expected imaginary part.
- */
- private static void assertComplex(double a, double b,
- String name, UnaryOperator<Complex> operation,
- double x, double y) {
- assertComplex(a, b, name, operation, x, y, 1);
- }
-
- /**
- * Assert the operation on the complex number is equal to the expected value.
- *
- * <p>The results are considered equal within the provided units of least
- * precision. The maximum count of numbers allowed between the two values is
- * {@code maxUlps - 1}.
- *
- * @param a Real part.
- * @param b Imaginary part.
- * @param name The operation name.
- * @param operation The operation.
- * @param x Expected real part.
- * @param y Expected imaginary part.
- * @param maxUlps the maximum units of least precision between the two values
- */
- private static void assertComplex(double a, double b,
- String name, UnaryOperator<Complex> operation,
- double x, double y, long maxUlps) {
- final Complex c = Complex.ofCartesian(a, b);
- final Complex e = Complex.ofCartesian(x, y);
- CReferenceTest.assertComplex(c, name, operation, e, maxUlps);
- }
-
/**
* Assert the operation on the complex number is equal to the expected value.
*
@@ -448,8 +407,8 @@ class ComplexEdgeCaseTest {
// im = 2 sin(2y) / e^2|x|
// This can be computed when e^2|x| only just overflows.
// Set a case where e^2|x| overflows but the imaginary can be computed
- double x = 709.783 / 2;
- double y = Math.PI / 4;
+ final double x = 709.783 / 2;
+ final double y = Math.PI / 4;
Assertions.assertEquals(1.0, Math.sin(2 * y), 1e-16);
Assertions.assertEquals(Double.POSITIVE_INFINITY, Math.exp(2 * x));
// As computed by GNU g++
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 35ebd50b..23665b71 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
@@ -295,7 +295,7 @@ class ComplexTest {
@Test
void testAbs() {
final Complex input = Complex.ofCartesian(3.0, 4.0);
- double z = TestUtils.assertSame(input, "abs", Complex::abs, ComplexFunctions::abs);
+ final double z = TestUtils.assertSame(input, "abs", Complex::abs, ComplexFunctions::abs);
Assertions.assertEquals(5.0, z);
}
@@ -384,7 +384,7 @@ class ComplexTest {
@Test
void testNorm() {
final Complex z = Complex.ofCartesian(3.0, 4.0);
- double norm = TestUtils.assertSame(z, "norm", Complex::norm, ComplexFunctions::norm);
+ final double norm = TestUtils.assertSame(z, "norm", Complex::norm, ComplexFunctions::norm);
Assertions.assertEquals(25.0, norm);
}
@@ -1529,7 +1529,7 @@ class ComplexTest {
theta += pi / 12;
final Complex z = Complex.ofPolar(r, theta);
final Complex sqrtz = Complex.ofPolar(Math.sqrt(r), theta / 2);
- Complex actual = TestUtils.assertSame(z, "sqrt", Complex::sqrt, ComplexFunctions::sqrt);
+ final Complex actual = TestUtils.assertSame(z, "sqrt", Complex::sqrt, ComplexFunctions::sqrt);
TestUtils.assertEquals(sqrtz, actual, tol);
}
}
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 a601df6c..63556769 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
@@ -61,7 +61,7 @@ public final class TestUtils {
/**
* The option for how to process test data lines flagged (prefixed)
- * with the {@code ;} character
+ * with the {@code ;} character.
*/
public enum TestDataFlagOption {
/** Ignore the line. */
@@ -449,7 +449,7 @@ public final class TestUtils {
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);
+ Assertions.assertEquals(r1, r2, () -> "ToDouble function mismatch: " + name);
return r1;
}
@@ -460,7 +460,7 @@ public final class TestUtils {
* @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.
+ * @param operation2 Condition operation on the complex real and imaginary parts.
* @return Result boolean from the given operation.
*/
public static boolean assertSame(Complex c,
@@ -470,7 +470,7 @@ public final class TestUtils {
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);
+ Assertions.assertEquals(b1, b2, () -> "Predicate mismatch: " + name);
return b1;
}
}