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/26 19:12:59 UTC
[commons-numbers] branch complex-gsoc-2022 updated: NUMBERS-188: Refactor complex binary functions to static methods
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch complex-gsoc-2022
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
The following commit(s) were added to refs/heads/complex-gsoc-2022 by this push:
new 7386b095 NUMBERS-188: Refactor complex binary functions to static methods
7386b095 is described below
commit 7386b095d767ef184bed9bbbf6b2d84f46a617f4
Author: sumanth-rajkumar <53...@users.noreply.github.com>
AuthorDate: Tue Jul 26 15:12:54 2022 -0400
NUMBERS-188: Refactor complex binary functions to static methods
Refactored binary methods accepting a Complex object as the second argument.
---
.../apache/commons/numbers/complex/Complex.java | 304 +---------------
.../numbers/complex/ComplexBinaryOperator.java | 47 +++
.../commons/numbers/complex/ComplexFunctions.java | 389 ++++++++++++++++++++-
.../commons/numbers/complex/CReferenceTest.java | 61 ++--
.../commons/numbers/complex/CStandardTest.java | 70 ++--
.../numbers/complex/ComplexEdgeCaseTest.java | 65 ++--
.../commons/numbers/complex/ComplexTest.java | 59 ++--
.../apache/commons/numbers/complex/TestUtils.java | 25 ++
8 files changed, 595 insertions(+), 425 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 4dc47117..9fc10012 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,12 +82,6 @@ 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);
- /** Exponent offset in IEEE754 representation. */
- private static final int EXPONENT_OFFSET = 1023;
- /** 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;
/** Serializable version identifier. */
private static final long serialVersionUID = 20180201L;
@@ -705,129 +699,7 @@ public final class Complex implements Serializable {
* @see <a href="http://mathworld.wolfram.com/ComplexMultiplication.html">Complex Muliplication</a>
*/
public Complex multiply(Complex factor) {
- return multiply(real, imaginary, factor.real, factor.imaginary);
- }
-
- /**
- * Returns a {@code Complex} whose value is:
- * <pre>
- * (a + i b)(c + i d) = (ac - bd) + i (ad + bc)</pre>
- *
- * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
- *
- * @param re1 Real component of first number.
- * @param im1 Imaginary component of first number.
- * @param re2 Real component of second number.
- * @param im2 Imaginary component of second number.
- * @return (a + b i)(c + d i).
- */
- private static Complex multiply(double re1, double im1, double re2, double im2) {
- double a = re1;
- double b = im1;
- double c = re2;
- double d = im2;
- final double ac = a * c;
- final double bd = b * d;
- final double ad = a * d;
- final double bc = b * c;
- double x = ac - bd;
- double y = ad + bc;
-
- // --------------
- // NaN can occur if:
- // - any of (a,b,c,d) are NaN (for NaN or Infinite complex numbers)
- // - a multiplication of infinity by zero (ac,bd,ad,bc).
- // - a subtraction of infinity from infinity (e.g. ac - bd)
- // Note that (ac,bd,ad,bc) can be infinite due to overflow.
- //
- // Detect a NaN result and perform correction.
- //
- // Modification from the listing in ISO C99 G.5.1 (6)
- // Do not correct infinity multiplied by zero. This is left as NaN.
- // --------------
-
- if (Double.isNaN(x) && Double.isNaN(y)) {
- // Recover infinities that computed as NaN+iNaN ...
- boolean recalc = false;
- if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
- isNotZero(c, d)) {
- // This complex is infinite.
- // "Box" the infinity and change NaNs in the other factor to 0.
- a = boxInfinity(a);
- b = boxInfinity(b);
- c = changeNaNtoZero(c);
- d = changeNaNtoZero(d);
- recalc = true;
- }
- if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
- isNotZero(a, b)) {
- // The other complex is infinite.
- // "Box" the infinity and change NaNs in the other factor to 0.
- c = boxInfinity(c);
- d = boxInfinity(d);
- a = changeNaNtoZero(a);
- b = changeNaNtoZero(b);
- recalc = true;
- }
- if (!recalc && (Double.isInfinite(ac) || Double.isInfinite(bd) ||
- Double.isInfinite(ad) || Double.isInfinite(bc))) {
- // The result overflowed to infinity.
- // Recover infinities from overflow by changing NaNs to 0 ...
- a = changeNaNtoZero(a);
- b = changeNaNtoZero(b);
- c = changeNaNtoZero(c);
- d = changeNaNtoZero(d);
- recalc = true;
- }
- if (recalc) {
- x = Double.POSITIVE_INFINITY * (a * c - b * d);
- y = Double.POSITIVE_INFINITY * (a * d + b * c);
- }
- }
- return new Complex(x, y);
- }
-
- /**
- * Box values for the real or imaginary component of an infinite complex number.
- * Any infinite value will be returned as one. Non-infinite values will be returned as zero.
- * The sign is maintained.
- *
- * <pre>
- * inf = 1
- * -inf = -1
- * x = 0
- * -x = -0
- * </pre>
- *
- * @param component the component
- * @return The boxed value
- */
- private static double boxInfinity(double component) {
- return Math.copySign(Double.isInfinite(component) ? 1.0 : 0.0, component);
- }
-
- /**
- * Checks if the complex number is not zero.
- *
- * @param real the real component
- * @param imaginary the imaginary component
- * @return true if the complex is not zero
- */
- private static boolean isNotZero(double real, double imaginary) {
- // The use of equals is deliberate.
- // This method must distinguish NaN from zero thus ruling out:
- // (real != 0.0 || imaginary != 0.0)
- return !(real == 0.0 && imaginary == 0.0);
- }
-
- /**
- * Change NaN to zero preserving the sign; otherwise return the value.
- *
- * @param value the value
- * @return The new value
- */
- private static double changeNaNtoZero(double value) {
- return Double.isNaN(value) ? Math.copySign(0.0, value) : value;
+ return ComplexFunctions.multiply(real, imaginary, factor.real, factor.imaginary, Complex::ofCartesian);
}
/**
@@ -901,90 +773,7 @@ public final class Complex implements Serializable {
* @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
*/
public Complex divide(Complex divisor) {
- return divide(real, imaginary, divisor.real, divisor.imaginary);
- }
-
- /**
- * Returns a {@code Complex} whose value is:
- * <pre>
- * <code>
- * a + i b (ac + bd) + i (bc - ad)
- * ------- = -----------------------
- * c + i d c<sup>2</sup> + d<sup>2</sup>
- * </code>
- * </pre>
- *
- * <p>Recalculates to recover infinities as specified in C99
- * standard G.5.1. Method is fully in accordance with
- * C++11 standards for complex numbers.</p>
- *
- * <p>Note: In the event of divide by zero this method produces the same result
- * as dividing by a real-only zero using {@link #divide(double)}.
- *
- * @param re1 Real component of first number.
- * @param im1 Imaginary component of first number.
- * @param re2 Real component of second number.
- * @param im2 Imaginary component of second number.
- * @return (a + i b) / (c + i d).
- * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
- * @see #divide(double)
- */
- private static Complex divide(double re1, double im1, double re2, double im2) {
- double a = re1;
- double b = im1;
- double c = re2;
- double d = im2;
- int ilogbw = 0;
- // Get the exponent to scale the divisor parts to the range [1, 2).
- final int exponent = getScale(c, d);
- if (exponent <= Double.MAX_EXPONENT) {
- ilogbw = exponent;
- c = Math.scalb(c, -ilogbw);
- d = Math.scalb(d, -ilogbw);
- }
- final double denom = c * c + d * d;
-
- // Note: Modification from the listing in ISO C99 G.5.1 (8):
- // Avoid overflow if a or b are very big.
- // Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
- // when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
- // (bc - ad) with large negative values.
- // Use the maximum exponent as an approximation to the magnitude.
- if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
- ilogbw -= 2;
- a /= 4;
- b /= 4;
- }
-
- double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
- double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
- // Recover infinities and zeros that computed as NaN+iNaN
- // the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
- if (Double.isNaN(x) && Double.isNaN(y)) {
- if (denom == 0.0 &&
- (!Double.isNaN(a) || !Double.isNaN(b))) {
- // nonzero/zero
- // This case produces the same result as divide by a real-only zero
- // using Complex.divide(+/-0.0)
- x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
- y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
- } else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
- Double.isFinite(c) && Double.isFinite(d)) {
- // infinite/finite
- a = boxInfinity(a);
- b = boxInfinity(b);
- x = Double.POSITIVE_INFINITY * (a * c + b * d);
- y = Double.POSITIVE_INFINITY * (b * c - a * d);
- } else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
- Double.isFinite(a) && Double.isFinite(b)) {
- // finite/infinite
- c = boxInfinity(c);
- d = boxInfinity(d);
- x = 0.0 * (a * c + b * d);
- y = 0.0 * (b * c - a * d);
- }
- }
- return new Complex(x, y);
+ return ComplexFunctions.divide(real, imaginary, divisor.real, divisor.imaginary, Complex::ofCartesian);
}
/**
@@ -1907,93 +1696,4 @@ public final class Complex implements Serializable {
private static boolean negative(double d) {
return ComplexFunctions.negative(d);
}
-
- /**
- * Returns a scale suitable for use with {@link Math#scalb(double, int)} to normalise
- * the number to the interval {@code [1, 2)}.
- *
- * <p>The scale is typically the largest unbiased exponent used in the representation of the
- * two numbers. In contrast to {@link Math#getExponent(double)} this handles
- * sub-normal numbers by computing the number of leading zeros in the mantissa
- * and shifting the unbiased exponent. The result is that for all finite, non-zero,
- * numbers {@code a, b}, the magnitude of {@code scalb(x, -getScale(a, b))} is
- * always in the range {@code [1, 2)}, where {@code x = max(|a|, |b|)}.
- *
- * <p>This method is a functional equivalent of the c function ilogb(double) adapted for
- * two input arguments.
- *
- * <p>The result is to be used to scale a complex number using {@link Math#scalb(double, int)}.
- * Hence the special case of both zero arguments is handled using the return value for NaN
- * as zero cannot be scaled. This is different from {@link Math#getExponent(double)}
- * or {@link #getMaxExponent(double, double)}.
- *
- * <p>Special cases:
- *
- * <ul>
- * <li>If either argument is NaN or infinite, then the result is
- * {@link Double#MAX_EXPONENT} + 1.
- * <li>If both arguments are zero, then the result is
- * {@link Double#MAX_EXPONENT} + 1.
- * </ul>
- *
- * @param a the first value
- * @param b the second value
- * @return The maximum unbiased exponent of the values to be used for scaling
- * @see Math#getExponent(double)
- * @see Math#scalb(double, int)
- * @see <a href="http://www.cplusplus.com/reference/cmath/ilogb/">ilogb</a>
- */
- private static int getScale(double a, double b) {
- // Only interested in the exponent and mantissa so remove the sign bit
- final long x = Double.doubleToRawLongBits(a) & UNSIGN_MASK;
- final long y = Double.doubleToRawLongBits(b) & UNSIGN_MASK;
- // Only interested in the maximum
- final long bits = Math.max(x, y);
- // Get the unbiased exponent
- int exp = ((int) (bits >>> 52)) - EXPONENT_OFFSET;
-
- // No case to distinguish nan/inf
- // Handle sub-normal numbers
- if (exp == Double.MIN_EXPONENT - 1) {
- // Special case for zero, return as nan/inf to indicate scaling is not possible
- if (bits == 0) {
- return Double.MAX_EXPONENT + 1;
- }
- // A sub-normal number has an exponent below -1022. The amount below
- // is defined by the number of shifts of the most significant bit in
- // the mantissa that is required to get a 1 at position 53 (i.e. as
- // if it were a normal number with assumed leading bit)
- final long mantissa = bits & MANTISSA_MASK;
- exp -= Long.numberOfLeadingZeros(mantissa << 12);
- }
- return exp;
- }
-
- /**
- * Returns the largest unbiased exponent used in the representation of the
- * two numbers. Special cases:
- *
- * <ul>
- * <li>If either argument is NaN or infinite, then the result is
- * {@link Double#MAX_EXPONENT} + 1.
- * <li>If both arguments are zero or subnormal, then the result is
- * {@link Double#MIN_EXPONENT} -1.
- * </ul>
- *
- * <p>This is used by {@link #divide(double, double, double, double)} as
- * a simple detection that a number may overflow if multiplied
- * by a value in the interval [1, 2).
- *
- * @param a the first value
- * @param b the second value
- * @return The maximum unbiased exponent of the values.
- * @see Math#getExponent(double)
- * @see #divide(double, double, double, double)
- */
- private static int getMaxExponent(double a, double b) {
- // This could return:
- // Math.getExponent(Math.max(Math.abs(a), Math.abs(b)))
- // A speed test is required to determine performance.
- return Math.max(Math.getExponent(a), Math.getExponent(b));
- }
}
diff --git a/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexBinaryOperator.java b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexBinaryOperator.java
new file mode 100644
index 00000000..6714e78e
--- /dev/null
+++ b/commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexBinaryOperator.java
@@ -0,0 +1,47 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.numbers.complex;
+
+/**
+ * Represents a binary operation on a Cartesian form of a complex number \( a + ib \)
+ * where \( a \) and \( b \) are real numbers represented as two {@code double}
+ * parts. The operation creates a complex number result; the result is supplied
+ * to a terminating consumer function which may return an object representation
+ * of the complex result.
+ *
+ * <p>This is a functional interface whose functional method is
+ * {@link #apply(double, double, double, double, ComplexSink)}.
+ *
+ * @param <R> The type of the complex result
+ * @since 1.1
+ */
+@FunctionalInterface
+public interface ComplexBinaryOperator<R> {
+
+ /**
+ * Represents an operator that accepts real and imaginary parts of two complex numbers and supplies the complex result to the provided consumer.
+ *
+ * @param real1 Real part \( a \) of the first complex number \( (a +ib) \).
+ * @param imaginary1 Imaginary part \( b \) of the first complex number \( (a +ib) \).
+ * @param real2 Real part \( a \) of the second complex number \( (a +ib) \).
+ * @param imaginary2 Imaginary part \( b \) of the second complex number \( (a +ib) \).
+ * @param out Consumer for the complex result.
+ * @return the object returned by the provided consumer.
+ */
+ R apply(double real1, double imaginary1, double real2, double imaginary2, ComplexSink<R> out);
+}
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 460e1f39..6dd614fb 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
@@ -95,6 +95,8 @@ public final class ComplexFunctions {
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;
/** 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}.*/
@@ -345,7 +347,7 @@ public final class ComplexFunctions {
*
* @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 action Consumer for the conjugate of the complex number.
* @param <R> the return type of the supplied action.
* @return the object returned by the supplied action.
*/
@@ -362,7 +364,7 @@ public final class ComplexFunctions {
*
* @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 action Consumer for the negation of the complex number.
* @param <R> the return type of the supplied action.
* @return the object returned by the supplied action.
*/
@@ -382,7 +384,7 @@ public final class ComplexFunctions {
*
* @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 action Consumer for the projection 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)
@@ -396,6 +398,259 @@ public final class ComplexFunctions {
return action.apply(real, imaginary);
}
+ /**
+ * Returns a {@code Object} whose value is {@code (real1 + real2, imaginary1 + imaginary2)}.
+ * Implements the formula:
+ *
+ * <p>\[ (a + i b) + (c + i d) = (a + c) + i (b + d) \]
+ *
+ * @param real1 Real part \( a \) of the first complex number \( (a +ib) \).
+ * @param imaginary1 Imaginary part \( b \) of the first complex number \( (a +ib) \).
+ * @param real2 Real part \( a \) of the second complex number \( (a +ib) \).
+ * @param imaginary2 Imaginary part \( b \) of the second complex number \( (a +ib) \).
+ * @param action Consumer for the addition result.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://mathworld.wolfram.com/ComplexAddition.html">Complex Addition</a>
+ */
+ public static <R> R add(double real1, double imaginary1,
+ double real2, double imaginary2,
+ ComplexSink<R> action) {
+ return action.apply(real1 + real2,
+ imaginary1 + imaginary2);
+ }
+
+ /**
+ * Returns a {@code Object} whose value is {@code (real1 - real2, imaginary1 - imaginary2)}.
+ * Implements the formula:
+ *
+ * <p>\[ (a + i b) - (c + i d) = (a - c) + i (b - d) \]
+ *
+ * @param real1 Real part \( a \) of the first complex number \( (a +ib) \).
+ * @param imaginary1 Imaginary part \( b \) of the first complex number \( (a +ib) \).
+ * @param real2 Real part \( a \) of the second complex number \( (a +ib) \).
+ * @param imaginary2 Imaginary part \( b \) of the second complex number \( (a +ib) \).
+ * @param action Consumer for the subtraction result.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://mathworld.wolfram.com/ComplexSubtraction.html">Complex Subtraction</a>
+ */
+ public static <R> R subtract(double real1, double imaginary1,
+ double real2, double imaginary2,
+ ComplexSink<R> action) {
+ return action.apply(real1 - real2,
+ imaginary1 - imaginary2);
+ }
+
+ /**
+ * Returns a {@code Object} whose value is {@code (real1 + i*imaginary1) * (real2 + i*imaginary2))}.
+ * Implements the formula:
+ *
+ * <p>\[ (a + i b)(c + i d) = (ac - bd) + i (ad + bc) \]
+ *
+ * <p>Recalculates to recover infinities as specified in C99 standard G.5.1.
+ *
+ * @param real1 Real part \( a \) of the first complex number \( (a +ib) \).
+ * @param imaginary1 Imaginary part \( b \) of the first complex number \( (a +ib) \).
+ * @param real2 Real part \( a \) of the second complex number \( (a +ib) \).
+ * @param imaginary2 Imaginary part \( b \) of the second complex number \( (a +ib) \).
+ * @param action Consumer for the multiplication result.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://mathworld.wolfram.com/ComplexMultiplication.html">Complex Muliplication</a>
+ */
+ public static <R> R multiply(double real1, double imaginary1,
+ double real2, double imaginary2,
+ ComplexSink<R> action) {
+ double a = real1;
+ double b = imaginary1;
+ double c = real2;
+ double d = imaginary2;
+ final double ac = a * c;
+ final double bd = b * d;
+ final double ad = a * d;
+ final double bc = b * c;
+ double x = ac - bd;
+ double y = ad + bc;
+
+ // --------------
+ // NaN can occur if:
+ // - any of (a,b,c,d) are NaN (for NaN or Infinite complex numbers)
+ // - a multiplication of infinity by zero (ac,bd,ad,bc).
+ // - a subtraction of infinity from infinity (e.g. ac - bd)
+ // Note that (ac,bd,ad,bc) can be infinite due to overflow.
+ //
+ // Detect a NaN result and perform correction.
+ //
+ // Modification from the listing in ISO C99 G.5.1 (6)
+ // Do not correct infinity multiplied by zero. This is left as NaN.
+ // --------------
+
+ if (Double.isNaN(x) && Double.isNaN(y)) {
+ // Recover infinities that computed as NaN+iNaN ...
+ boolean recalc = false;
+ if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
+ isNotZero(c, d)) {
+ // This complex is infinite.
+ // "Box" the infinity and change NaNs in the other factor to 0.
+ a = boxInfinity(a);
+ b = boxInfinity(b);
+ c = changeNaNtoZero(c);
+ d = changeNaNtoZero(d);
+ recalc = true;
+ }
+ if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
+ isNotZero(a, b)) {
+ // The other complex is infinite.
+ // "Box" the infinity and change NaNs in the other factor to 0.
+ c = boxInfinity(c);
+ d = boxInfinity(d);
+ a = changeNaNtoZero(a);
+ b = changeNaNtoZero(b);
+ recalc = true;
+ }
+ if (!recalc && (Double.isInfinite(ac) || Double.isInfinite(bd) ||
+ Double.isInfinite(ad) || Double.isInfinite(bc))) {
+ // The result overflowed to infinity.
+ // Recover infinities from overflow by changing NaNs to 0 ...
+ a = changeNaNtoZero(a);
+ b = changeNaNtoZero(b);
+ c = changeNaNtoZero(c);
+ d = changeNaNtoZero(d);
+ recalc = true;
+ }
+ if (recalc) {
+ x = Double.POSITIVE_INFINITY * (a * c - b * d);
+ y = Double.POSITIVE_INFINITY * (a * d + b * c);
+ }
+ }
+ return action.apply(x, y);
+ }
+
+ /**
+ * Box values for the real or imaginary component of an infinite complex number.
+ * Any infinite value will be returned as one. Non-infinite values will be returned as zero.
+ * The sign is maintained.
+ *
+ * <pre>
+ * inf = 1
+ * -inf = -1
+ * x = 0
+ * -x = -0
+ * </pre>
+ *
+ * @param component the component
+ * @return The boxed value
+ */
+ private static double boxInfinity(double component) {
+ return Math.copySign(Double.isInfinite(component) ? 1.0 : 0.0, component);
+ }
+
+ /**
+ * Checks if the complex number is not zero.
+ *
+ * @param real the real component
+ * @param imaginary the imaginary component
+ * @return true if the complex is not zero
+ */
+ private static boolean isNotZero(double real, double imaginary) {
+ // The use of equals is deliberate.
+ // This method must distinguish NaN from zero thus ruling out:
+ // (real != 0.0 || imaginary != 0.0)
+ return !(real == 0.0 && imaginary == 0.0);
+ }
+
+ /**
+ * Change NaN to zero preserving the sign; otherwise return the value.
+ *
+ * @param value the value
+ * @return The new value
+ */
+ private static double changeNaNtoZero(double value) {
+ return Double.isNaN(value) ? Math.copySign(0.0, value) : value;
+ }
+
+ /**
+ * Returns a {@code Object} whose value is {@code (real1 + i*imaginary1)/(real2 + i*imaginary2)}.
+ * Implements the formula:
+ *
+ * <p>\[ \frac{a + i b}{c + i d} = \frac{(ac + bd) + i (bc - ad)}{c^2+d^2} \]
+ *
+ * <p>Re-calculates NaN result values to recover infinities as specified in C99 standard G.5.1.
+ *
+ * <p>Note: In the event of divide by zero this method produces the same result
+ * as dividing by a real-only zero using (add divide reference)
+ *
+ * @param real1 Real part \( a \) of the first complex number \( (a +ib) \).
+ * @param imaginary1 Imaginary part \( b \) of the first complex number \( (a +ib) \).
+ * @param real2 Real part \( a \) of the second complex number \( (a +ib) \).
+ * @param imaginary2 Imaginary part \( b \) of the second complex number \( (a +ib) \).
+ * @param action Consumer for the division result.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see <a href="http://mathworld.wolfram.com/ComplexDivision.html">Complex Division</a>
+ */
+ //TODO - add divide reference once moved to ComplexFunctions
+ public static <R> R divide(double real1, double imaginary1,
+ double real2, double imaginary2,
+ ComplexSink<R> action) {
+ double a = real1;
+ double b = imaginary1;
+ double c = real2;
+ double d = imaginary2;
+ int ilogbw = 0;
+ // Get the exponent to scale the divisor parts to the range [1, 2).
+ final int exponent = getScale(c, d);
+ if (exponent <= Double.MAX_EXPONENT) {
+ ilogbw = exponent;
+ c = Math.scalb(c, -ilogbw);
+ d = Math.scalb(d, -ilogbw);
+ }
+ final double denom = c * c + d * d;
+
+ // Note: Modification from the listing in ISO C99 G.5.1 (8):
+ // Avoid overflow if a or b are very big.
+ // Since (c, d) in the range [1, 2) the sum (ac + bd) could overflow
+ // when (a, b) are both above (Double.MAX_VALUE / 4). The same applies to
+ // (bc - ad) with large negative values.
+ // Use the maximum exponent as an approximation to the magnitude.
+ if (getMaxExponent(a, b) > Double.MAX_EXPONENT - 2) {
+ ilogbw -= 2;
+ a /= 4;
+ b /= 4;
+ }
+
+ double x = Math.scalb((a * c + b * d) / denom, -ilogbw);
+ double y = Math.scalb((b * c - a * d) / denom, -ilogbw);
+ // Recover infinities and zeros that computed as NaN+iNaN
+ // the only cases are nonzero/zero, infinite/finite, and finite/infinite, ...
+ if (Double.isNaN(x) && Double.isNaN(y)) {
+ if (denom == 0.0 &&
+ (!Double.isNaN(a) || !Double.isNaN(b))) {
+ // nonzero/zero
+ // This case produces the same result as divide by a real-only zero
+ // using Complex.divide(+/-0.0)
+ x = Math.copySign(Double.POSITIVE_INFINITY, c) * a;
+ y = Math.copySign(Double.POSITIVE_INFINITY, c) * b;
+ } else if ((Double.isInfinite(a) || Double.isInfinite(b)) &&
+ Double.isFinite(c) && Double.isFinite(d)) {
+ // infinite/finite
+ a = boxInfinity(a);
+ b = boxInfinity(b);
+ x = Double.POSITIVE_INFINITY * (a * c + b * d);
+ y = Double.POSITIVE_INFINITY * (b * c - a * d);
+ } else if ((Double.isInfinite(c) || Double.isInfinite(d)) &&
+ Double.isFinite(a) && Double.isFinite(b)) {
+ // finite/infinite
+ c = boxInfinity(c);
+ d = boxInfinity(d);
+ x = 0.0 * (a * c + b * d);
+ y = 0.0 * (b * c - a * d);
+ }
+ }
+ return action.apply(x, y);
+ }
+
/**
* Returns the
* <a href="http://mathworld.wolfram.com/ExponentialFunction.html">
@@ -677,6 +932,45 @@ public final class ComplexFunctions {
return action.apply(re, arg(real, imaginary));
}
+ /**
+ * Returns the complex power of the first complex number raised to the power of
+ * second complex number.
+ * Implements the formula:
+ *
+ * <p>\[ z^x = e^{x \ln(z)} \]
+ *
+ * <p>If the complex number is zero then this method returns zero if the second complex number is positive
+ * in the real component and zero in the imaginary component;
+ * otherwise it returns NaN + iNaN.
+ *
+ * @param real1 Real part \( a \) of the first complex number \( (a +ib) \).
+ * @param imaginary1 Imaginary part \( b \) of the first complex number \( (a +ib) \).
+ * @param real2 Real part \( a \) of the second complex number \( (a +ib) \).
+ * @param imaginary2 Imaginary part \( b \) of the second complex number \( (a +ib) \).
+ * @param action Consumer for the power result.
+ * @param <R> the return type of the supplied action.
+ * @return the object returned by the supplied action.
+ * @see #log(double, double, ComplexSink)
+ * @see #multiply(double, double, double, double, ComplexSink)
+ * @see #exp(double, double, ComplexSink)
+ * @see <a href="http://mathworld.wolfram.com/ComplexExponentiation.html">Complex exponentiation</a>
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
+ */
+ public static <R> R pow(double real1, double imaginary1, double real2, double imaginary2, ComplexSink<R> action) {
+ if (real1 == 0 &&
+ imaginary1 == 0) {
+ // This value is zero. Test the other.
+ if (real2 > 0 &&
+ imaginary2 == 0) {
+ // 0 raised to positive number is 0
+ return action.apply(0, 0);
+ }
+ // 0 raised to anything else is NaN
+ return action.apply(Double.NaN, Double.NaN);
+ }
+ return log(real1, imaginary1, (x, y) -> multiply(x, y, real2, imaginary2, (a, b) -> exp(a, b, action)));
+ }
+
/**
* Returns the
* <a href="http://mathworld.wolfram.com/SquareRoot.html">
@@ -2551,6 +2845,95 @@ public final class ComplexFunctions {
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)}.
+ *
+ * <p>The scale is typically the largest unbiased exponent used in the representation of the
+ * two numbers. In contrast to {@link Math#getExponent(double)} this handles
+ * sub-normal numbers by computing the number of leading zeros in the mantissa
+ * and shifting the unbiased exponent. The result is that for all finite, non-zero,
+ * numbers {@code a, b}, the magnitude of {@code scalb(x, -getScale(a, b))} is
+ * always in the range {@code [1, 2)}, where {@code x = max(|a|, |b|)}.
+ *
+ * <p>This method is a functional equivalent of the c function ilogb(double) adapted for
+ * two input arguments.
+ *
+ * <p>The result is to be used to scale a complex number using {@link Math#scalb(double, int)}.
+ * Hence the special case of both zero arguments is handled using the return value for NaN
+ * as zero cannot be scaled. This is different from {@link Math#getExponent(double)}
+ * or {@link #getMaxExponent(double, double)}.
+ *
+ * <p>Special cases:
+ *
+ * <ul>
+ * <li>If either argument is NaN or infinite, then the result is
+ * {@link Double#MAX_EXPONENT} + 1.
+ * <li>If both arguments are zero, then the result is
+ * {@link Double#MAX_EXPONENT} + 1.
+ * </ul>
+ *
+ * @param a the first value
+ * @param b the second value
+ * @return The maximum unbiased exponent of the values to be used for scaling
+ * @see Math#getExponent(double)
+ * @see Math#scalb(double, int)
+ * @see <a href="http://www.cplusplus.com/reference/cmath/ilogb/">ilogb</a>
+ */
+ private static int getScale(double a, double b) {
+ // Only interested in the exponent and mantissa so remove the sign bit
+ final long x = Double.doubleToRawLongBits(a) & UNSIGN_MASK;
+ final long y = Double.doubleToRawLongBits(b) & UNSIGN_MASK;
+ // Only interested in the maximum
+ final long bits = Math.max(x, y);
+ // Get the unbiased exponent
+ int exp = ((int) (bits >>> 52)) - EXPONENT_OFFSET;
+
+ // No case to distinguish nan/inf
+ // Handle sub-normal numbers
+ if (exp == Double.MIN_EXPONENT - 1) {
+ // Special case for zero, return as nan/inf to indicate scaling is not possible
+ if (bits == 0) {
+ return Double.MAX_EXPONENT + 1;
+ }
+ // A sub-normal number has an exponent below -1022. The amount below
+ // is defined by the number of shifts of the most significant bit in
+ // the mantissa that is required to get a 1 at position 53 (i.e. as
+ // if it were a normal number with assumed leading bit)
+ final long mantissa = bits & MANTISSA_MASK;
+ exp -= Long.numberOfLeadingZeros(mantissa << 12);
+ }
+ return exp;
+ }
+
+ /**
+ * Returns the largest unbiased exponent used in the representation of the
+ * two numbers. Special cases:
+ *
+ * <ul>
+ * <li>If either argument is NaN or infinite, then the result is
+ * {@link Double#MAX_EXPONENT} + 1.
+ * <li>If both arguments are zero or subnormal, then the result is
+ * {@link Double#MIN_EXPONENT} -1.
+ * </ul>
+ *
+ * <p>This is used by {@link #divide(double, double, double, double, ComplexSink)} as
+ * a simple detection that a number may overflow if multiplied
+ * by a value in the interval [1, 2).
+ *
+ * @param a the first value
+ * @param b the second value
+ * @return The maximum unbiased exponent of the values.
+ * @see Math#getExponent(double)
+ * @see #divide(double, double, double, double, ComplexSink)
+ */
+ private static int getMaxExponent(double a, double b) {
+ // This could return:
+ // Math.getExponent(Math.max(Math.abs(a), Math.abs(b)))
+ // 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.
*
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 62dc190e..a56cce8f 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
@@ -173,29 +173,6 @@ class CReferenceTest {
assertEquals(() -> c + "." + name + "(): imaginary", expected.imag(), z.imag(), maxUlps);
}
- /**
- * 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}.
- *
- * <p>Numbers must have the same sign. Thus -0.0 and 0.0 are never equal.
- *
- * @param c Input number.
- * @param name the operation name
- * @param operation the operation
- * @param expected Expected result.
- * @param maxUlps the maximum units of least precision between the two values
- */
- static void assertComplex(Complex c,
- String name, UnaryOperator<Complex> operation,
- Complex expected, long maxUlps) {
- final Complex z = operation.apply(c);
- assertEquals(() -> c + "." + name + "(): real", expected.real(), z.real(), maxUlps);
- assertEquals(() -> c + "." + name + "(): imaginary", expected.imag(), z.imag(), maxUlps);
- }
-
/**
* Assert the operation on the complex numbers is equal to the expected value.
*
@@ -205,17 +182,22 @@ class CReferenceTest {
*
* <p>Numbers must have the same sign. Thus -0.0 and 0.0 are never equal.
*
- * @param c1 First number.
- * @param c2 Second number.
- * @param name the operation name
- * @param operation the operation
- * @param expected Expected real part.
- * @param maxUlps the maximum units of least precision between the two values
+ * <p>Assert the operation on the complex numbers is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts of two complex numbers.
+ *
+ * @param c1 First input number.
+ * @param c2 Second input number.
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex objects.
+ * @param operation2 Operation on the complex real and imaginary parts of two complex numbers.
+ * @param expected Expected result.
+ * @param maxUlps Maximum units of least precision between the two values.
*/
static void assertComplex(Complex c1, Complex c2,
- String name, BiFunction<Complex, Complex, Complex> operation,
+ String name, BiFunction<Complex, Complex, Complex> operation1,
+ ComplexBinaryOperator<ComplexNumber> operation2,
Complex expected, long maxUlps) {
- final Complex z = operation.apply(c1, c2);
+ final Complex z = TestUtils.assertSame(c1, c2, name, operation1, operation2);
assertEquals(() -> c1 + "." + name + c2 + ": real", expected.real(), z.real(), maxUlps);
assertEquals(() -> c1 + "." + name + c2 + ": imaginary", expected.imag(), z.imag(), maxUlps);
}
@@ -241,16 +223,19 @@ class CReferenceTest {
/**
* Assert the operation using the data loaded from test resources.
*
- * @param name the operation name
- * @param operation the operation
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex objects.
+ * @param operation2 Operation on the complex real and imaginary parts of two complex numbers.
* @param maxUlps the maximum units of least precision between the two values
*/
private static void assertBiOperation(String name,
- BiFunction<Complex, Complex, Complex> operation, long maxUlps) {
+ BiFunction<Complex, Complex, Complex> operation1,
+ ComplexBinaryOperator<ComplexNumber> operation2,
+ long maxUlps) {
final List<Complex[]> data = loadTestData(name);
final long ulps = getTestUlps(maxUlps);
for (final Complex[] triple : data) {
- assertComplex(triple[0], triple[1], name, operation, triple[2], ulps);
+ assertComplex(triple[0], triple[1], name, operation1, operation2, triple[2], ulps);
}
}
@@ -352,16 +337,16 @@ class CReferenceTest {
@Test
void testMultiply() {
- assertBiOperation("multiply", Complex::multiply, 0);
+ assertBiOperation("multiply", Complex::multiply, ComplexFunctions::multiply, 0);
}
@Test
void testDivide() {
- assertBiOperation("divide", Complex::divide, 7);
+ assertBiOperation("divide", Complex::divide, ComplexFunctions::divide, 7);
}
@Test
void testPowComplex() {
- assertBiOperation("pow", Complex::pow, 9);
+ assertBiOperation("pow", Complex::pow, ComplexFunctions::pow, 9);
}
}
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 85bc5590..620e3a92 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
@@ -203,17 +203,23 @@ class CStandardTest {
/**
* Assert the operation on the two complex numbers.
*
- * @param c1 the first complex
- * @param c2 the second complex
- * @param operation the operation
- * @param operationName the operation name
- * @param expected the expected
- * @param expectedName the expected name
+ * <p>Assert the operation on the complex numbers is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts of two complex numbers.
+ *
+ * @param c1 First input complex number.
+ * @param c2 Second input complex number.
+ * @param operation1 Operation on the Complex objects.
+ * @param operation2 Operation on the complex real and imaginary parts of two complex numbers.
+ * @param operationName Operation name.
+ * @param expected Expected complex number.
+ * @param expectedName Expected name.
*/
private static void assertOperation(Complex c1, Complex c2,
- BiFunction<Complex, Complex, Complex> operation, String operationName,
- Predicate<Complex> expected, String expectedName) {
- final Complex z = operation.apply(c1, c2);
+ BiFunction<Complex, Complex, Complex> operation1,
+ ComplexBinaryOperator<ComplexNumber> operation2,
+ String operationName, Predicate<Complex> expected,
+ String expectedName) {
+ final Complex z = TestUtils.assertSame(c1, c2, operationName, operation1, operation2);
Assertions.assertTrue(expected.test(z),
() -> String.format("%s expected: %s %s %s = %s", expectedName, c1, operationName, c2, z));
}
@@ -844,18 +850,18 @@ class CStandardTest {
// infinity, then the result of the * operator is an infinity;"
for (final Complex z : infinites) {
for (final Complex w : infinites) {
- assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, "Inf");
- assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(z, w, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(w, z, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
}
for (final Complex w : nonZeroFinites) {
- assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, "Inf");
- assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(z, w, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(w, z, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
}
// C.99 refers to non-zero finites.
// Infer that Complex multiplication of zero with infinites is not defined.
for (final Complex w : zeroFinites) {
- assertOperation(z, w, Complex::multiply, "*", Complex::isNaN, "NaN");
- assertOperation(w, z, Complex::multiply, "*", Complex::isNaN, "NaN");
+ assertOperation(z, w, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isNaN, "NaN");
+ assertOperation(w, z, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isNaN, "NaN");
}
}
@@ -871,16 +877,16 @@ class CStandardTest {
// Check multiply with (NaN,NaN) is not corrected
final double[] values = {0, 1, inf, negInf, nan};
for (final Complex z : createCombinations(values, c -> true)) {
- assertOperation(z, NAN, Complex::multiply, "*", Complex::isNaN, "NaN");
- assertOperation(NAN, z, Complex::multiply, "*", Complex::isNaN, "NaN");
+ assertOperation(z, NAN, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isNaN, "NaN");
+ assertOperation(NAN, z, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isNaN, "NaN");
}
// Test multiply cases which result in overflow are corrected to infinity
- assertOperation(maxMax, maxMax, Complex::multiply, "*", Complex::isInfinite, "Inf");
- assertOperation(maxNan, maxNan, Complex::multiply, "*", Complex::isInfinite, "Inf");
- assertOperation(nanMax, maxNan, Complex::multiply, "*", Complex::isInfinite, "Inf");
- assertOperation(maxNan, nanMax, Complex::multiply, "*", Complex::isInfinite, "Inf");
- assertOperation(nanMax, nanMax, Complex::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(maxMax, maxMax, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(maxNan, maxNan, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(nanMax, maxNan, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(maxNan, nanMax, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
+ assertOperation(nanMax, nanMax, Complex::multiply, ComplexFunctions::multiply, "*", Complex::isInfinite, "Inf");
}
/**
@@ -899,14 +905,14 @@ class CStandardTest {
// result of the / operator is an infinity;"
for (final Complex z : infinites) {
for (final Complex w : nonZeroFinites) {
- assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", Complex::isInfinite, "Inf");
}
for (final Complex w : zeroFinites) {
- assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", Complex::isInfinite, "Inf");
}
// Check inf/inf cannot be done
for (final Complex w : infinites) {
- assertOperation(z, w, Complex::divide, "/", Complex::isNaN, "NaN");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", Complex::isNaN, "NaN");
}
}
@@ -914,7 +920,7 @@ class CStandardTest {
// result of the / operator is a zero;"
for (final Complex z : finites) {
for (final Complex w : infinites) {
- assertOperation(z, w, Complex::divide, "/", CStandardTest::isZero, "Zero");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", CStandardTest::isZero, "Zero");
}
}
@@ -922,10 +928,10 @@ class CStandardTest {
// a zero, then the result of the / operator is an infinity."
for (final Complex w : zeroFinites) {
for (final Complex z : nonZeroFinites) {
- assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", Complex::isInfinite, "Inf");
}
for (final Complex z : infinites) {
- assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", Complex::isInfinite, "Inf");
}
}
@@ -934,19 +940,19 @@ class CStandardTest {
// infinite.
for (final Complex w : nans) {
for (final Complex z : finites) {
- assertOperation(z, w, Complex::divide, "/", c -> NAN.equals(c), "(NaN,NaN)");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", c -> NAN.equals(c), "(NaN,NaN)");
}
for (final Complex z : infinites) {
- assertOperation(z, w, Complex::divide, "/", c -> NAN.equals(c), "(NaN,NaN)");
+ assertOperation(z, w, Complex::divide, ComplexFunctions::divide, "/", c -> NAN.equals(c), "(NaN,NaN)");
}
}
// Check (NaN,NaN) divide is not corrected for the edge case of divide by zero or infinite
for (final Complex w : zeroFinites) {
- assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, "NaN");
+ assertOperation(NAN, w, Complex::divide, ComplexFunctions::divide, "/", Complex::isNaN, "NaN");
}
for (final Complex w : infinites) {
- assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, "NaN");
+ assertOperation(NAN, w, Complex::divide, ComplexFunctions::divide, "/", Complex::isNaN, "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 72cc88d3..6759d35e 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
@@ -96,19 +96,24 @@ class ComplexEdgeCaseTest {
*
* <p>The results are considered equal if there are no floating-point values between them.
*
+ * <p>Assert the operation on the complex numbers is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts.
+ *
* @param a Real part of first number.
* @param b Imaginary part of first number.
* @param c Real part of second number.
* @param d Imaginary part of second number.
* @param name The operation name.
- * @param operation The operation.
+ * @param operation1 Operation on the Complex object.
+ * @param operation2 Operation on the complex real and imaginary parts.
* @param x Expected real part.
* @param y Expected imaginary part.
*/
private static void assertComplex(double a, double b, double c, double d,
- String name, BiFunction<Complex, Complex, Complex> operation,
+ String name, BiFunction<Complex, Complex, Complex> operation1,
+ ComplexBinaryOperator<ComplexNumber> operation2,
double x, double y) {
- assertComplex(a, b, c, d, name, operation, x, y, 1);
+ assertComplex(a, b, c, d, name, operation1, operation2, x, y, 1);
}
/**
@@ -118,23 +123,28 @@ class ComplexEdgeCaseTest {
* precision. The maximum count of numbers allowed between the two values is
* {@code maxUlps - 1}.
*
+ * <p>Assert the operation on the complex numbers is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts.
+ *
* @param a Real part of first number.
* @param b Imaginary part of first number.
* @param c Real part of second number.
* @param d Imaginary part of second number.
* @param name The operation name
- * @param operation the operation
+ * @param operation1 Operation on the Complex objects.
+ * @param operation2 Operation on the complex real and imaginary parts.
* @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, double c, double d,
- String name, BiFunction<Complex, Complex, Complex> operation,
+ String name, BiFunction<Complex, Complex, Complex> operation1,
+ ComplexBinaryOperator<ComplexNumber> operation2,
double x, double y, long maxUlps) {
final Complex c1 = Complex.ofCartesian(a, b);
final Complex c2 = Complex.ofCartesian(c, d);
final Complex e = Complex.ofCartesian(x, y);
- CReferenceTest.assertComplex(c1, c2, name, operation, e, maxUlps);
+ CReferenceTest.assertComplex(c1, c2, name, operation1, operation2, e, maxUlps);
}
@Test
@@ -657,7 +667,8 @@ class ComplexEdgeCaseTest {
@Test
void testDivide() {
final String name = "divide";
- final BiFunction<Complex, Complex, Complex> operation = Complex::divide;
+ final BiFunction<Complex, Complex, Complex> operation1 = Complex::divide;
+ final ComplexBinaryOperator<ComplexNumber> operation2 = ComplexFunctions::divide;
// Should be able to divide by a complex whose absolute (c*c+d*d)
// overflows or underflows including all sub-normal numbers.
@@ -674,30 +685,34 @@ class ComplexEdgeCaseTest {
// In other words the result is (x+iy) / (x+iy) = (1+i0)
// The result is the same if imaginary is zero (i.e. a real only divide)
- assertComplex(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, name, operation, 1.0, 0.0);
- assertComplex(Double.MAX_VALUE, 0.0, Double.MAX_VALUE, 0.0, name, operation, 1.0, 0.0);
+ assertComplex(Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, Double.MAX_VALUE, name, operation1, operation2, 1.0,
+ 0.0);
+ assertComplex(Double.MAX_VALUE, 0.0, Double.MAX_VALUE, 0.0, name, operation1, operation2, 1.0, 0.0);
- assertComplex(1.0, 1.0, 1.0, 1.0, name, operation, 1.0, 0.0);
- assertComplex(1.0, 0.0, 1.0, 0.0, name, operation, 1.0, 0.0);
+ assertComplex(1.0, 1.0, 1.0, 1.0, name, operation1, operation2, 1.0, 0.0);
+ assertComplex(1.0, 0.0, 1.0, 0.0, name, operation1, operation2, 1.0, 0.0);
// Should work for all small values
x = Double.MIN_NORMAL;
while (x != 0) {
- assertComplex(x, x, x, x, name, operation, 1.0, 0.0);
- assertComplex(x, 0, x, 0, name, operation, 1.0, 0.0);
+ assertComplex(x, x, x, x, name, operation1, operation2, 1.0, 0.0);
+ assertComplex(x, 0, x, 0, name, operation1, operation2, 1.0, 0.0);
x /= 2;
}
// Some cases of not self-divide
- assertComplex(1, 1, Double.MIN_VALUE, Double.MIN_VALUE, name, operation, inf, 0);
+ assertComplex(1, 1, Double.MIN_VALUE, Double.MIN_VALUE, name, operation1, operation2, inf, 0);
// As computed by GNU g++
- assertComplex(Double.MIN_NORMAL, Double.MIN_NORMAL, Double.MIN_VALUE, Double.MIN_VALUE, name, operation, 4503599627370496L, 0);
- assertComplex(Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation, 2.2204460492503131e-16, 0);
+ assertComplex(Double.MIN_NORMAL, Double.MIN_NORMAL, Double.MIN_VALUE, Double.MIN_VALUE, name, operation1, operation2,
+ 4503599627370496L, 0);
+ assertComplex(Double.MIN_VALUE, Double.MIN_VALUE, Double.MIN_NORMAL, Double.MIN_NORMAL, name, operation1, operation2,
+ 2.2204460492503131e-16, 0);
}
@Test
void testPow() {
final String name = "pow";
- final BiFunction<Complex, Complex, Complex> operation = Complex::pow;
+ final BiFunction<Complex, Complex, Complex> operation1 = Complex::pow;
+ final ComplexBinaryOperator<ComplexNumber> operation2 = ComplexFunctions::pow;
// pow(Complex) is log().multiply(Complex).exp()
// All are overflow safe and handle infinities as defined in the C99 standard.
@@ -706,15 +721,15 @@ class ComplexEdgeCaseTest {
// using other library implementations.
// Test NaN
- assertComplex(1, 1, nan, nan, name, operation, nan, nan);
- assertComplex(nan, nan, 1, 1, name, operation, nan, nan);
- assertComplex(nan, 1, 1, 1, name, operation, nan, nan);
- assertComplex(1, nan, 1, 1, name, operation, nan, nan);
- assertComplex(1, 1, nan, 1, name, operation, nan, nan);
- assertComplex(1, 1, 1, nan, name, operation, nan, nan);
+ assertComplex(1, 1, nan, nan, name, operation1, operation2, nan, nan);
+ assertComplex(nan, nan, 1, 1, name, operation1, operation2, nan, nan);
+ assertComplex(nan, 1, 1, 1, name, operation1, operation2, nan, nan);
+ assertComplex(1, nan, 1, 1, name, operation1, operation2, nan, nan);
+ assertComplex(1, 1, nan, 1, name, operation1, operation2, nan, nan);
+ assertComplex(1, 1, 1, nan, name, operation1, operation2, nan, nan);
// Test overflow.
- assertComplex(Double.MAX_VALUE, 1, 2, 2, name, operation, inf, -inf);
- assertComplex(1, Double.MAX_VALUE, 2, 2, name, operation, -inf, inf);
+ assertComplex(Double.MAX_VALUE, 1, 2, 2, name, operation1, operation2, inf, -inf);
+ assertComplex(1, Double.MAX_VALUE, 2, 2, name, operation1, operation2, -inf, inf);
}
}
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 23665b71..02653eb6 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
@@ -531,7 +531,7 @@ class ComplexTest {
void testAdd() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
final Complex y = Complex.ofCartesian(5.0, 6.0);
- final Complex z = x.add(y);
+ final Complex z = TestUtils.assertSame(x, y, "add", Complex::add, ComplexFunctions::add);
Assertions.assertEquals(8.0, z.getReal());
Assertions.assertEquals(10.0, z.getImaginary());
}
@@ -540,12 +540,13 @@ class ComplexTest {
void testAddInf() {
Complex x = Complex.ofCartesian(1, 1);
final Complex z = Complex.ofCartesian(inf, 0);
- final Complex w = x.add(z);
+ final Complex w = TestUtils.assertSame(x, z, "add", Complex::add, ComplexFunctions::add);
Assertions.assertEquals(1, w.getImaginary());
Assertions.assertEquals(inf, w.getReal());
x = Complex.ofCartesian(neginf, 0);
- Assertions.assertTrue(Double.isNaN(x.add(z).getReal()));
+ final Complex result = TestUtils.assertSame(x, z, "add", Complex::add, ComplexFunctions::add);
+ Assertions.assertTrue(Double.isNaN(result.getReal()));
}
@Test
@@ -644,7 +645,7 @@ class ComplexTest {
void testSubtract() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
final Complex y = Complex.ofCartesian(5.0, 7.0);
- final Complex z = x.subtract(y);
+ final Complex z = TestUtils.assertSame(x, y, "subtract", Complex::subtract, ComplexFunctions::subtract);
Assertions.assertEquals(-2.0, z.getReal());
Assertions.assertEquals(-3.0, z.getImaginary());
}
@@ -653,11 +654,11 @@ class ComplexTest {
void testSubtractInf() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
final Complex y = Complex.ofCartesian(inf, 7.0);
- Complex z = x.subtract(y);
+ Complex z = TestUtils.assertSame(x, y, "subtract", Complex::subtract, ComplexFunctions::subtract);
Assertions.assertEquals(neginf, z.getReal());
Assertions.assertEquals(-3.0, z.getImaginary());
- z = y.subtract(y);
+ z = TestUtils.assertSame(y, y, "subtract", Complex::subtract, ComplexFunctions::subtract);
Assertions.assertEquals(nan, z.getReal());
Assertions.assertEquals(0.0, z.getImaginary());
}
@@ -844,22 +845,27 @@ class ComplexTest {
void testMultiply() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
final Complex y = Complex.ofCartesian(5.0, 6.0);
- final Complex z = x.multiply(y);
+ final Complex z = TestUtils.assertSame(x, y, "multiply", Complex::multiply, ComplexFunctions::multiply);
Assertions.assertEquals(-9.0, z.getReal());
Assertions.assertEquals(38.0, z.getImaginary());
}
@Test
void testMultiplyInfInf() {
- final Complex z = infInf.multiply(infInf);
+ final Complex z = TestUtils.assertSame(infInf, infInf, "multiply", Complex::multiply, ComplexFunctions::multiply);
// Assert.assertTrue(z.isNaN()); // MATH-620
Assertions.assertTrue(TestUtils.assertSame(z, "isInfinite", Complex::isInfinite, ComplexFunctions::isInfinite));
// Expected results from g++:
- Assertions.assertEquals(Complex.ofCartesian(nan, inf), infInf.multiply(infInf));
- Assertions.assertEquals(Complex.ofCartesian(inf, nan), infInf.multiply(infNegInf));
- Assertions.assertEquals(Complex.ofCartesian(-inf, nan), infInf.multiply(negInfInf));
- Assertions.assertEquals(Complex.ofCartesian(nan, -inf), infInf.multiply(negInfNegInf));
+ assertMultiply(Complex.ofCartesian(nan, inf), infInf, infInf);
+ assertMultiply(Complex.ofCartesian(inf, nan), infInf, infNegInf);
+ assertMultiply(Complex.ofCartesian(-inf, nan), infInf, negInfInf);
+ assertMultiply(Complex.ofCartesian(nan, -inf), infInf, negInfNegInf);
+ }
+
+ private static void assertMultiply(Complex expected, Complex input1, Complex input2) {
+ final Complex actual = TestUtils.assertSame(input1, input2, "multiply", Complex::multiply, ComplexFunctions::multiply);
+ TestUtils.assertSame(expected, actual);
}
@Test
@@ -1094,7 +1100,7 @@ class ComplexTest {
void testDivide() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
final Complex y = Complex.ofCartesian(5.0, 6.0);
- final Complex z = x.divide(y);
+ final Complex z = TestUtils.assertSame(x, y, "divide", Complex::divide, ComplexFunctions::divide);
Assertions.assertEquals(39.0 / 61.0, z.getReal());
Assertions.assertEquals(2.0 / 61.0, z.getImaginary());
}
@@ -1102,35 +1108,35 @@ class ComplexTest {
@Test
void testDivideZero() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
- final Complex z = x.divide(Complex.ZERO);
+ final Complex z = TestUtils.assertSame(x, Complex.ZERO, "divide", Complex::divide, ComplexFunctions::divide);
Assertions.assertEquals(INF, z);
}
@Test
void testDivideZeroZero() {
final Complex x = Complex.ofCartesian(0.0, 0.0);
- final Complex z = x.divide(Complex.ZERO);
+ final Complex z = TestUtils.assertSame(x, Complex.ZERO, "divide", Complex::divide, ComplexFunctions::divide);
Assertions.assertEquals(NAN, z);
}
@Test
void testDivideNaN() {
final Complex x = Complex.ofCartesian(3.0, 4.0);
- final Complex z = x.divide(NAN);
+ final Complex z = TestUtils.assertSame(x, NAN, "divide", Complex::divide, ComplexFunctions::divide);
Assertions.assertTrue(TestUtils.assertSame(z, "isNaN", Complex::isNaN, ComplexFunctions::isNaN));
}
@Test
void testDivideNanInf() {
- Complex z = oneInf.divide(Complex.ONE);
+ Complex z = TestUtils.assertSame(oneInf, Complex.ONE, "divide", Complex::divide, ComplexFunctions::divide);
Assertions.assertTrue(Double.isNaN(z.getReal()));
Assertions.assertEquals(inf, z.getImaginary());
- z = negInfNegInf.divide(oneNan);
+ z = TestUtils.assertSame(negInfNegInf, oneNan, "divide", Complex::divide, ComplexFunctions::divide);
Assertions.assertTrue(Double.isNaN(z.getReal()));
Assertions.assertTrue(Double.isNaN(z.getImaginary()));
- z = negInfInf.divide(Complex.ONE);
+ z = TestUtils.assertSame(negInfInf, Complex.ONE, "divide", Complex::divide, ComplexFunctions::divide);
Assertions.assertTrue(Double.isInfinite(z.getReal()));
Assertions.assertTrue(Double.isInfinite(z.getImaginary()));
}
@@ -1401,7 +1407,7 @@ class ComplexTest {
final Complex c = Complex.ofCartesian(a, b);
for (final double arg : arguments) {
final Complex y = c.divideImaginary(arg);
- Complex z = c.divide(ofImaginary(arg));
+ Complex z = TestUtils.assertSame(c, ofImaginary(arg), "divide", Complex::divide, ComplexFunctions::divide);
final boolean expectedFailure = (expectedFailures & 0x1) == 1;
expectedFailures >>>= 1;
// If divide by zero then the divide(Complex) method matches divide by real.
@@ -1451,7 +1457,8 @@ class ComplexTest {
final Complex x = Complex.ofCartesian(3, 4);
final double yDouble = 5.0;
final Complex yComplex = ofReal(yDouble);
- Assertions.assertEquals(x.pow(yComplex), x.pow(yDouble));
+ final Complex expected = TestUtils.assertSame(x, yComplex, "pow", Complex::pow, ComplexFunctions::pow);
+ Assertions.assertEquals(expected, x.pow(yDouble));
}
@Test
@@ -1459,7 +1466,7 @@ class ComplexTest {
// Hits the edge case when real == 0 but imaginary != 0
final Complex x = Complex.ofCartesian(0, 1);
final Complex z = Complex.ofCartesian(2, 3);
- final Complex c = x.pow(z);
+ final Complex c = TestUtils.assertSame(x, z, "pow", Complex::pow, ComplexFunctions::pow);
// Answer from g++
Assertions.assertEquals(-0.008983291021129429, c.getReal());
Assertions.assertEquals(1.1001358594835313e-18, c.getImaginary());
@@ -1476,7 +1483,7 @@ class ComplexTest {
private static void assertPowComplexZeroBase(double re, double im, Complex expected) {
final Complex z = Complex.ofCartesian(re, im);
- final Complex c = Complex.ZERO.pow(z);
+ final Complex c = TestUtils.assertSame(Complex.ZERO, z, "pow", Complex::pow, ComplexFunctions::pow);
Assertions.assertEquals(expected, c);
}
@@ -1507,7 +1514,8 @@ class ComplexTest {
final Complex x = NAN;
final double yDouble = 5.0;
final Complex yComplex = ofReal(yDouble);
- Assertions.assertEquals(x.pow(yComplex), x.pow(yDouble));
+ final Complex expected = TestUtils.assertSame(x, yComplex, "pow", Complex::pow, ComplexFunctions::pow);
+ Assertions.assertEquals(expected, x.pow(yDouble));
}
@Test
@@ -1515,7 +1523,8 @@ class ComplexTest {
final Complex x = Complex.ofCartesian(3, 4);
final double yDouble = Double.NaN;
final Complex yComplex = ofReal(yDouble);
- Assertions.assertEquals(x.pow(yComplex), x.pow(yDouble));
+ final Complex expected = TestUtils.assertSame(x, yComplex, "pow", Complex::pow, ComplexFunctions::pow);
+ Assertions.assertEquals(expected, x.pow(yDouble));
}
@Test
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 63556769..e67bec5f 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
@@ -26,6 +26,7 @@ import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.util.ArrayList;
import java.util.List;
+import java.util.function.BiFunction;
import java.util.function.Consumer;
import java.util.function.DoubleBinaryOperator;
import java.util.function.Predicate;
@@ -432,6 +433,30 @@ public final class TestUtils {
return z;
}
+ /**
+ * Assert the operation on the complex numbers is <em>exactly</em> equal to the operation on
+ * complex real and imaginary parts of two complex numbers.
+ *
+ * @param c1 First input complex number.
+ * @param c2 Second input complex number.
+ * @param name Operation name.
+ * @param operation1 Operation on the Complex objects.
+ * @param operation2 Operation on the complex real and imaginary parts of two complex numbers.
+ * @return Result complex number from the given operation.
+ */
+ public static Complex assertSame(Complex c1,
+ Complex c2,
+ String name,
+ BiFunction<Complex, Complex, Complex> operation1,
+ ComplexBinaryOperator<ComplexNumber> operation2) {
+ final Complex z = operation1.apply(c1, c2);
+ // Test operation2 produces the exact same result
+ final ComplexNumber z2 = operation2.apply(c1.real(), c1.imag(), c2.real(), c2.imag(), ComplexNumber::new);
+ Assertions.assertEquals(z.real(), z2.getReal(), () -> "Binary operator mismatch: " + name + " real");
+ Assertions.assertEquals(z.imag(), z2.getImaginary(), () -> "Binary 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.