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 2019/12/10 23:57:22 UTC
[commons-numbers] branch master updated (9b3c706 -> 32318c6)
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a change to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git.
from 9b3c706 Use identities for tanh() with real / imaginary parts of zero.
new ff1b7d1 Ensure log() and log10() avoid overflow.
new 4429840 Added a ComplexEdgeCaseTest.
new 32318c6 Removed unused field and clarified multiply() comment regarding NaN.
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 | 54 +++-
.../commons/numbers/complex/CReferenceTest.java | 6 +-
.../commons/numbers/complex/CStandardTest.java | 3 +-
.../numbers/complex/ComplexEdgeCaseTest.java | 300 +++++++++++++++++++++
4 files changed, 352 insertions(+), 11 deletions(-)
create mode 100644 commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
[commons-numbers] 03/03: Removed unused field and clarified
multiply() comment regarding NaN.
Posted by ah...@apache.org.
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 32318c672ce084bd2a6f4a791b0f10f210e22f2a
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Tue Dec 10 23:57:18 2019 +0000
Removed unused field and clarified multiply() comment regarding NaN.
---
.../test/java/org/apache/commons/numbers/complex/CStandardTest.java | 3 +--
1 file changed, 1 insertion(+), 2 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 11a0bc2..1a6deb3 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
@@ -49,7 +49,6 @@ public class CStandardTest {
private static final Complex zeroNegInf = complex(0, negInf);
private static final Complex zeroNaN = complex(0, nan);
private static final Complex zeroPiTwo = complex(0.0, piOverTwo);
- private static final Complex negZeroPiTwo = complex(-0.0, piOverTwo);
private static final Complex negZeroZero = complex(-0.0, 0);
private static final Complex negZeroNaN = complex(-0.0, nan);
private static final Complex negI = complex(0.0, -1.0);
@@ -637,7 +636,7 @@ public class CStandardTest {
// ISO C Standard in Annex G is missing an explicit definition of how to handle NaNs.
// We will assume multiplication by (nan,nan) is not allowed.
// It is undefined how to multiply when a complex has only one NaN component.
- // The reference implementation allows it.
+ // The reference implementation in Annex G allows it.
// The GNU g++ compiler computes:
// (1e300 + i 1e300) * (1e30 + i NAN) = inf + i inf
[commons-numbers] 02/03: Added a ComplexEdgeCaseTest.
Posted by ah...@apache.org.
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit 44298409892b82031cef4540274f913e52d583f0
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Tue Dec 10 23:56:01 2019 +0000
Added a ComplexEdgeCaseTest.
Overflow and underflow are tested for log, exp and sqrt.
The test is incomplete. Placeholders are left for the other C99
functions to test boundary conditions.
---
.../commons/numbers/complex/CReferenceTest.java | 6 +-
.../numbers/complex/ComplexEdgeCaseTest.java | 300 +++++++++++++++++++++
2 files changed, 303 insertions(+), 3 deletions(-)
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 e11a4e3..57ad59f 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
@@ -96,7 +96,7 @@ public class CReferenceTest {
* @param actual the actual
* @param maxUlps the maximum units of least precision between the two values
*/
- private static void assertEquals(Supplier<String> msg, double expected, double actual, long maxUlps) {
+ static void assertEquals(Supplier<String> msg, double expected, double actual, long maxUlps) {
final long e = Double.doubleToLongBits(expected);
final long a = Double.doubleToLongBits(actual);
@@ -160,7 +160,7 @@ public class CReferenceTest {
* @param expected Expected result.
* @param maxUlps the maximum units of least precision between the two values
*/
- private static void assertComplex(Complex c,
+ static void assertComplex(Complex c,
String name, UnaryOperator<Complex> operation,
Complex expected, long maxUlps) {
final Complex z = operation.apply(c);
@@ -184,7 +184,7 @@ public class CReferenceTest {
* @param expected Expected real part.
* @param maxUlps the maximum units of least precision between the two values
*/
- private static void assertComplex(Complex c1, Complex c2,
+ static void assertComplex(Complex c1, Complex c2,
String name, BiFunction<Complex, Complex, Complex> operation,
Complex expected, long maxUlps) {
final Complex z = operation.apply(c1, c2);
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
new file mode 100644
index 0000000..84d856f
--- /dev/null
+++ b/commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java
@@ -0,0 +1,300 @@
+/*
+ * 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;
+
+import org.junit.jupiter.api.Assertions;
+import org.junit.jupiter.api.Test;
+
+import java.util.function.BiFunction;
+import java.util.function.UnaryOperator;
+
+/**
+ * Edge case tests for the functions defined by the C.99 standard for complex numbers
+ * defined in ISO/IEC 9899, Annex G.
+ *
+ * <p>The test contained here are specifically written to target edge cases of finite valued
+ * input values that cause overflow/underflow during the computation.
+ *
+ * <p>The test data is generated from a known implementation of the standard.
+ *
+ * @see <a href="http://www.open-std.org/JTC1/SC22/WG14/www/standards">
+ * ISO/IEC 9899 - Programming languages - C</a>
+ */
+public 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 z = operation.apply(c);
+ CReferenceTest.assertComplex(c, name, operation, z, maxUlps);
+ }
+
+ /**
+ * Assert the operation on the complex numbers is equal to the expected value.
+ *
+ * <p>The results are considered equal if there are no floating-point values between them.
+ *
+ * @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 x Expected real part.
+ * @param y Expected imaginary part.
+ */
+ // CHECKSTYLE: stop ParameterNumberCheck
+ private static void assertComplex(double a, double b, double c, double d,
+ String name, BiFunction<Complex, Complex, Complex> operation,
+ double x, double y) {
+ assertComplex(a, b, c, d, name, operation, x, y, 1);
+ }
+
+ /**
+ * Assert the operation on the complex numbers 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 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 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,
+ double x, double y, long maxUlps) {
+ final Complex c1 = Complex.ofCartesian(a, b);
+ final Complex c2 = Complex.ofCartesian(c, d);
+ final Complex z = operation.apply(c1, c2);
+ CReferenceTest.assertComplex(c1, c2, name, operation, z, maxUlps);
+ }
+
+ @Test
+ public void testAcos() {
+ // acos(z) = (pi / 2) + i ln(iz + sqrt(1 - z^2))
+ // TODO
+ }
+
+ @Test
+ public void testAcosh() {
+ // Defined by acos() so edge cases are the same
+ // TODO
+ }
+
+ @Test
+ public void testAsinh() {
+ // asinh(z) = ln(z + sqrt(1 + z^2))
+ // Odd function: negative real cases defined by positive real cases
+ // TODO
+ }
+
+ @Test
+ public void testAtanh() {
+ // atanh(z) = (1/2) ln((1 + z) / (1 - z))
+ // Odd function: negative real cases defined by positive real cases
+ // TODO
+ }
+
+ @Test
+ public void testCosh() {
+ // cosh(a + b i) = cosh(a)cos(b) + i sinh(a)sin(b)
+ // Even function: negative real cases defined by positive real cases
+ // TODO
+ }
+
+ @Test
+ public void testSinh() {
+ // sinh(a + b i) = sinh(a)cos(b)) + i cosh(a)sin(b)
+ // Odd function: negative real cases defined by positive real cases
+ // TODO
+ }
+
+ @Test
+ public void testTanh() {
+ // 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
+ // TODO
+ }
+
+ @Test
+ public void testExp() {
+ final String name = "exp";
+ final UnaryOperator<Complex> operation = Complex::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);
+
+ // 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);
+ }
+
+ @Test
+ public void testLog() {
+ final String name = "log";
+ final UnaryOperator<Complex> operation = Complex::log;
+
+ // ln(a + b i) = ln(|a + b i|) + i arg(a + b i)
+ // |a + b i| = sqrt(a^2 + b^2)
+ // arg(a + b i) = Math.atan2(imaginary, real)
+
+ // Overflow if sqrt(a^2 + b^2) == inf.
+ // Matlab computes this.
+ assertComplex(-Double.MAX_VALUE, Double.MAX_VALUE, name, operation, 7.101292864836639e2, Math.PI * 3 / 4);
+ assertComplex(Double.MAX_VALUE, Double.MAX_VALUE, name, operation, 7.101292864836639e2, Math.PI / 4);
+ assertComplex(-Double.MAX_VALUE, Double.MAX_VALUE / 4, name, operation, 7.098130252042921e2, 2.896613990462929);
+ assertComplex(Double.MAX_VALUE, Double.MAX_VALUE / 4, name, operation, 7.098130252042921e2, 2.449786631268641e-1, 2);
+
+ // Underflow if sqrt(a^2 + b^2) == 0
+ assertComplex(-Double.MIN_VALUE, Double.MIN_VALUE, name, operation, -744.44007192138122, 2.3561944901923448);
+ assertComplex(Double.MIN_VALUE, Double.MIN_VALUE, name, operation, -744.44007192138122, 0.78539816339744828);
+ }
+
+ @Test
+ public void testSqrt() {
+ final String name = "sqrt";
+ final UnaryOperator<Complex> operation = Complex::sqrt;
+
+ // Computed in polar coordinates:
+ // real = (x^2 + y^2)^0.25 * cos(0.5 * atan2(y, x))
+ // imag = (x^2 + y^2)^0.25 * sin(0.5 * atan2(y, x))
+ // ---
+ // Note:
+ // If x is positive and y is +/-0.0 atan2 returns +/-0.
+ // If x is negative and y is +/-0.0 atan2 returns +/-PI.
+ // This causes problems as
+ // cos(0.5 * PI) = 6.123233995736766e-17
+ // assert: Math.cos(Math.acos(0)) != 0.0
+ // Thus an unchecked polar computation will produce incorrect output when
+ // there is no imaginary component and real is negative.
+ // The computation should be done for real only numbers separately.
+ // This condition is tested in the reference test against known data
+ // so not repeated here.
+ // ---
+
+ // Check overflow safe.
+ double a = Double.MAX_VALUE;
+ final double b = a / 4;
+ Assertions.assertEquals(inf, Complex.ofCartesian(a, b).abs(), "Expected overflow");
+ // Compute the expected new magnitude by expressing b as a scale factor of a:
+ // (x^2 + y^2)^0.25
+ // = sqrt(sqrt(a^2 + (b/a)^2 * a^2))
+ // = sqrt(sqrt((1+(b/a)^2) * a^2))
+ // = sqrt(sqrt((1+(b/a)^2))) * sqrt(a)
+ final double newAbs = Math.sqrt(Math.sqrt(1 + (b / a) * (b / a))) * Math.sqrt(a);
+ 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);
+
+ // In polar coords:
+ // real = sqrt(abs()) * Math.cos(arg() / 2)
+ // imag = sqrt(abs()) * Math.sin(arg() / 2)
+ // This is possible 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);
+
+ // 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);
+ }
+
+ // Note:
+ // multiply is tested in CStandardTest
+ // divide is tested in CStandardTest
+
+ @Test
+ public void testPow() {
+ final String name = "pow";
+ final BiFunction<Complex, Complex, Complex> operation = Complex::pow;
+
+ // pow(Complex) is log().multiply(Complex).exp()
+ // All are overflow safe and handle infinities as defined in the C99 standard.
+ // TODO: Test edge cases with:
+ // Double.MAX_VALUE, Double.MIN_NORMAL, Inf
+ // 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);
+ }
+}
[commons-numbers] 01/03: Ensure log() and log10() avoid overflow.
Posted by ah...@apache.org.
This is an automated email from the ASF dual-hosted git repository.
aherbert pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-numbers.git
commit ff1b7d1d197fac7ce58864a8c89d032132ee53f0
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Tue Dec 10 23:23:08 2019 +0000
Ensure log() and log10() avoid overflow.
---
.../apache/commons/numbers/complex/Complex.java | 54 +++++++++++++++++++---
1 file changed, 48 insertions(+), 6 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 c9cd5ea..0181189 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
@@ -96,6 +96,21 @@ public final class Complex implements Serializable {
}
/**
+ * Define a unary operation on a double.
+ * This is used in the log() and log10() functions.
+ */
+ @FunctionalInterface
+ private interface UnaryOperation {
+ /**
+ * Apply an operation to a value.
+ *
+ * @param value The value.
+ * @return The result.
+ */
+ double apply(double value);
+ }
+
+ /**
* Private default constructor.
*
* @param real Real part.
@@ -1469,13 +1484,11 @@ public final class Complex implements Serializable {
* @see #arg()
*/
public Complex log() {
- // All edge cases satisfied by the Math library
- return new Complex(Math.log(abs()),
- arg());
+ return log(Math::log);
}
/**
- * Compute the base 10 or
+ * Compute the base 10
* <a href="http://mathworld.wolfram.com/CommonLogarithm.html">
* common logarithm</a> of this complex number.
* Implements the formula:
@@ -1489,8 +1502,37 @@ public final class Complex implements Serializable {
* @see #arg()
*/
public Complex log10() {
- // All edge cases satisfied by the Math library
- return new Complex(Math.log10(abs()),
+ return log(Math::log10);
+ }
+
+ /**
+ * Compute the logarithm of this complex number using the provided function.
+ * Implements the formula:
+ * <pre>
+ * log(a + bi) = log(|a + b i|) + i arg(a + b i)
+ * </pre>
+ *
+ * @param log Log function.
+ * @return the logarithm of {@code this}.
+ * @see #abs()
+ * @see #arg()
+ */
+ private Complex log(UnaryOperation log) {
+ // All ISO C99 edge cases satisfied by the Math library.
+ // Make computation overflow safe.
+ final double abs = abs();
+ if (abs == Double.POSITIVE_INFINITY && isFinite()) {
+ // Edge-case where the |a + b i| overflows.
+ // |a + b i| = sqrt(a^2 + b^2)
+ // This can be scaled linearly.
+ // Scale the absolute and exploit:
+ // ln(abs / s) = ln(abs) - ln(s)
+ // ln(abs) = ln(abs / s) + ln(s)
+ final double s = Math.max(Math.abs(real), Math.abs(imaginary));
+ final double absOs = Math.hypot(real / s, imaginary / s);
+ return new Complex(log.apply(absOs) + log.apply(s), arg());
+ }
+ return new Complex(log.apply(abs),
arg());
}