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/20 17:57:51 UTC
[commons-numbers] branch master updated (9c0cb71 -> 32d8fc1)
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 9c0cb71 Use multiplyImaginary instead of multiply by Complex.I
new 4dd3fc8 Updated sinh/asinh using an overflow/underflow safe method.
new 96fe378 Updated cosh/acosh using an overflow/underflow safe method.
new 2c77d00 Fix unspecified sign in C99 standard for acos.
new 93d6890 Update the max ULP for acos/acosh from 36 to 2 for the C reference test.
new a2cc4ec Removed unused static asinh method and constants.
new 5346128 Updated atanh using an overflow/underflow safe method.
new 26b9ae7 Update the max ULP for atanh from 17 to 1 for the C reference test.
new bdaecba Use final
new 67dc19d Compute log() within a safe region without using the sqrt.
new 4bbbef9 Update the max ULP for log from 3 to 1 for the C reference test.
new 1d20e9f Removed unused constants
new 5d9156b Added atanh assumptions test to check the safe upper and lower limits.
new 12c0b57 Drop exact log10 test. The method of computation cannot be assumed.
new 913b012 Trailing whitespace.
new 9aa23ab Switch atanh if statements to fix PMD warning.
new 244f6e7 Remove use of log1p when value is very small.
new e335125 Remove unused/redundant methods.
new bf2b32c Moved the check for the safe region to a method.
new 45d1a46 Added edge case tests for atanh.
new dfcc118 Added edge case tests for acos.
new f5ba494 Add atanh assumption on x close to 1.
new e1ad2a0 Added edge case tests for asin.
new 8375852 Fix javadoc comment on SAFE_MIN/MAX
new 3366f35 Fix overflow in tanh. Added edge case tests.
new 1f7af2c Added links to Wolfram elementary functions where applicable.
new ba7fe3a Edge cases for cosh and sinh.
new 3614d51 Simplify sqrt() edge cases.
new e2e7b8b Simplify exp() edge cases.
new 1eae0a4 Simplify sinh, cosh and tanh
new 32d8fc1 Overflow test for pow()
The 30 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 | 1158 ++++++++++++--------
.../commons/numbers/complex/CReferenceTest.java | 10 +-
.../commons/numbers/complex/CStandardTest.java | 2 +-
.../numbers/complex/ComplexEdgeCaseTest.java | 164 ++-
.../commons/numbers/complex/ComplexTest.java | 358 +++---
5 files changed, 1080 insertions(+), 612 deletions(-)
[commons-numbers] 02/30: Updated cosh/acosh using an
overflow/underflow safe method.
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 96fe378a62aa7c5151c35721d05676c27b5164d6
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 11:30:35 2019 +0000
Updated cosh/acosh using an overflow/underflow safe method.
---
.../apache/commons/numbers/complex/Complex.java | 209 +++++++++++++++------
1 file changed, 148 insertions(+), 61 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 feb777c..010b40b 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
@@ -78,9 +78,14 @@ public final class Complex implements Serializable {
/** Base 10 logarithm of 2 (log10(2)). */
private static final double LOG10_2 = Math.log10(2);
- /** Crossover point to switch computation for asin factor A. */
+ /**
+ * 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;
- /** Crossover point to switch computation for asin factor B. */
+ /** 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.
@@ -1261,6 +1266,27 @@ public final class Complex implements Serializable {
* </code>
* </pre>
*
+ * <p>This is implemented using real {@code x} and imaginary {@code y} parts:</p>
+ * <pre>
+ * <code>
+ * acos(z) = acos(B) - i ln(A + sqrt(A<sup>2</sup>-1))
+ * A = 0.5 [ sqrt((x+1)<sup>2</sup>+y<sup>2</sup>) + sqrt((x-1)<sup>2</sup>+y<sup>2</sup>) ]
+ * B = 0.5 [ sqrt((x+1)<sup>2</sup>+y<sup>2</sup>) - sqrt((x-1)<sup>2</sup>+y<sup>2</sup>) ]
+ * </code>
+ * </pre>
+ *
+ * <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>}. The function is well
+ * defined over the entire complex number range, and produces accurate values even at the
+ * extremes due to special handling of overflow and underflow conditions.</p>
+ *
* @return the inverse cosine of this complex number.
*/
public Complex acos() {
@@ -1278,58 +1304,120 @@ public final class Complex implements Serializable {
* @param constructor Constructor.
* @return the inverse cosine of the complex number.
*/
- private static Complex acos(double real, double imaginary, ComplexConstructor constructor) {
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- // Special case for real numbers
- if (imaginary == 0 && Math.abs(real) <= 1) {
- return constructor.create(real == 0.0 ? PI_OVER_2 : Math.acos(real),
- Math.copySign(0, -imaginary));
- }
- // ISO C99: Preserve the equality
- // acos(conj(z)) = conj(acos(z))
- // by always computing on a positive imaginary Complex number.
- final double a = real;
- final double b = Math.abs(imaginary);
- final Complex z2 = multiply(a, b, a, b);
- // sqrt(1 - z^2)
- final Complex sqrt1mz2 = sqrt(1 - z2.real, -z2.imaginary);
- // Compute the rest inline to avoid Complex object creation.
- // (x + y i) = iz + sqrt(1 - z^2)
- final double x = -b + sqrt1mz2.real;
- final double y = a + sqrt1mz2.imaginary;
- // (re + im i) = (pi / 2) + i ln(iz + sqrt(1 - z^2))
- final double re = PI_OVER_2 - getArgument(x, y);
- final double im = Math.log(getAbsolute(x, y));
- // Map back to the correct sign
- return constructor.create(re, changeSign(im, imaginary));
+ private static Complex acos(final double real, final double imaginary,
+ final ComplexConstructor constructor) {
+ // Compute with positive values and determine sign at the end
+ double x = Math.abs(real);
+ double y = Math.abs(imaginary);
+ // 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;
}
- if (Double.isInfinite(imaginary)) {
- return constructor.create(PI_OVER_2, Math.copySign(Double.POSITIVE_INFINITY, -imaginary));
+ } else if (Double.isNaN(x)) {
+ if (isPosInfinite(y)) {
+ return constructor.create(x, -imaginary);
}
- // imaginary is NaN
- // Special case for real == 0
- return real == 0 ? constructor.create(PI_OVER_2, Double.NaN) : NAN;
- }
- if (Double.isInfinite(real)) {
- if (Double.isFinite(imaginary)) {
- final double re = real == Double.NEGATIVE_INFINITY ? Math.PI : 0;
- return constructor.create(re, Math.copySign(Double.POSITIVE_INFINITY, -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);
}
- if (Double.isInfinite(imaginary)) {
- final double re = real == Double.NEGATIVE_INFINITY ? PI_3_OVER_4 : PI_OVER_4;
- return constructor.create(re, Math.copySign(Double.POSITIVE_INFINITY, -imaginary));
+
+ double xp1 = x + 1;
+ double xm1 = x - 1;
+
+ if ((x < SAFE_MAX) && (x > SAFE_MIN) && (y < SAFE_MAX) && (y > SAFE_MIN)) {
+ double yy = y * y;
+ double r = Math.sqrt(xp1 * xp1 + yy);
+ double s = Math.sqrt(xm1 * xm1 + yy);
+ double a = 0.5 * (r + s);
+ double b = x / a;
+
+ if (b <= B_CROSSOVER) {
+ re = Math.acos(b);
+ } else {
+ 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);
+ double xoy = x / y;
+ im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
+ } else {
+ re = PI_OVER_2;
+ double a = Math.sqrt(1 + y * y);
+ im = 0.5 * Math.log1p(2 * y * (y + a));
+ }
}
- // imaginary is NaN
- // Swap real and imaginary
- return constructor.create(Double.NaN, real);
- }
- // real is NaN
- if (Double.isInfinite(imaginary)) {
- return constructor.create(Double.NaN, -imaginary);
}
- // optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
- return NAN;
+
+ return constructor.create(negative(real) ? Math.PI - re : re,
+ negative(imaginary) ? im : -im);
}
/**
@@ -1345,7 +1433,7 @@ public final class Complex implements Serializable {
* <p>This is implemented using real {@code x} and imaginary {@code y} parts:</p>
* <pre>
* <code>
- * asin(z) = asin(B) + i sign(y)log(A + sqrt(A<sup>2</sup>-1))
+ * asin(z) = asin(B) + i sign(y)ln(A + sqrt(A<sup>2</sup>-1))
* A = 0.5 [ sqrt((x+1)<sup>2</sup>+y<sup>2</sup>) + sqrt((x-1)<sup>2</sup>+y<sup>2</sup>) ]
* B = 0.5 [ sqrt((x+1)<sup>2</sup>+y<sup>2</sup>) - sqrt((x-1)<sup>2</sup>+y<sup>2</sup>) ]
* sign(y) = {@link Math#copySign(double,double) copySign(1.0, y)}
@@ -1384,7 +1472,8 @@ public final class Complex implements Serializable {
* @param constructor Constructor.
* @return the inverse sine of this complex number
*/
- private static Complex asin(double real, double imaginary, ComplexConstructor constructor) {
+ private static Complex asin(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);
@@ -1497,7 +1586,8 @@ public final class Complex implements Serializable {
}
}
- return constructor.create(changeSign(re, real), changeSign(im, imaginary));
+ return constructor.create(changeSign(re, real),
+ changeSign(im, imaginary));
}
/**
@@ -1733,16 +1823,13 @@ public final class Complex implements Serializable {
if (Double.isNaN(imaginary) && Double.isFinite(real)) {
return NAN;
}
- // ISO C99: Preserve the equality
- // acos(conj(z)) = conj(acos(z))
- // by always computing on a positive imaginary Complex number.
- return acos(real, Math.abs(imaginary), (re, im) ->
- // Set the sign appropriately for C99 equalities.
- (im > 0) ?
- // Multiply by -I and map back to the correct sign
- new Complex(im, changeSign(-re, imaginary)) :
+ return acos(real, imaginary, (re, im) ->
+ // Set the sign appropriately for real >= 0
+ (negative(im)) ?
// Multiply by I
- new Complex(-im, changeSign(re, imaginary))
+ new Complex(-im, re) :
+ // Multiply by -I
+ new Complex(im, -re)
);
}
[commons-numbers] 19/30: Added edge case tests for atanh.
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 45d1a4618ed0022ae669dd04561f00e6a7823312
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Dec 19 22:05:10 2019 +0000
Added edge case tests for atanh.
---
.../numbers/complex/ComplexEdgeCaseTest.java | 23 +++++++++++++++++++++-
1 file changed, 22 insertions(+), 1 deletion(-)
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 be47bbb..f94eddb 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
@@ -150,7 +150,28 @@ public class ComplexEdgeCaseTest {
public void testAtanh() {
// atanh(z) = (1/2) ln((1 + z) / (1 - z))
// Odd function: negative real cases defined by positive real cases
- // TODO
+ final String name = "atanh";
+ final UnaryOperator<Complex> operation = Complex::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.
+ // A medium value is used to test outside the range of the CReferenceTest.
+ // It hits an edge case when x is big and y > 1.
+ // The results have been generated using g++ -std=c++11 atanh.
+ double big = Math.sqrt(Double.MAX_VALUE) / 2;
+ double medium = 100;
+ double small = Math.sqrt(Double.MIN_NORMAL) * 2;
+ assertComplex(inf, big, name, operation, 0, 1.5707963267948966);
+ assertComplex(big, inf, name, operation, 0, 1.5707963267948966);
+ 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);
}
@Test
[commons-numbers] 30/30: Overflow test for pow()
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 32d8fc1e5c1f21a336ea3b20712777ec3a0de1be
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 17:29:27 2019 +0000
Overflow test for pow()
---
.../java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java | 4 ++++
1 file changed, 4 insertions(+)
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 bee0aff..4bfbabd 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
@@ -428,5 +428,9 @@ public class ComplexEdgeCaseTest {
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);
+
+ // Test overflow
+ assertComplex(Double.MAX_VALUE, 1, 2, 0, name, operation, inf, inf);
+ assertComplex(1, Double.MAX_VALUE, 2, 0, name, operation, -inf, inf);
}
}
[commons-numbers] 20/30: Added edge case tests for acos.
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 dfcc1188e1457e6f385996e4ebe8dafb94950e09
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Dec 19 22:50:04 2019 +0000
Added edge case tests for acos.
---
.../numbers/complex/ComplexEdgeCaseTest.java | 35 ++++++++++++++++++++--
1 file changed, 32 insertions(+), 3 deletions(-)
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 f94eddb..42ca2cb 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
@@ -130,7 +130,35 @@ public class ComplexEdgeCaseTest {
@Test
public void testAcos() {
// acos(z) = (pi / 2) + i ln(iz + sqrt(1 - z^2))
- // TODO
+ final String name = "acos";
+ final UnaryOperator<Complex> operation = Complex::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.
+ // A medium value is used to test outside the range of the CReferenceTest.
+ // The results have been generated using g++ -std=c++11 acos.
+ // xp1 * xm1 will overflow:
+ double huge = Math.sqrt(Double.MAX_VALUE) * 2;
+ double big = Math.sqrt(Double.MAX_VALUE) / 8;
+ double medium = 100;
+ 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);
+ // Additional cases to achieve full coverage
+ // xm1 = 0
+ assertComplex(1, small, name, operation, 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);
}
@Test
@@ -161,8 +189,6 @@ public class ComplexEdgeCaseTest {
double big = Math.sqrt(Double.MAX_VALUE) / 2;
double medium = 100;
double small = Math.sqrt(Double.MIN_NORMAL) * 2;
- assertComplex(inf, big, name, operation, 0, 1.5707963267948966);
- assertComplex(big, inf, name, operation, 0, 1.5707963267948966);
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);
@@ -172,6 +198,9 @@ public class ComplexEdgeCaseTest {
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);
+ // Additional cases to achieve full coverage
+ assertComplex(inf, big, name, operation, 0, 1.5707963267948966);
+ assertComplex(big, inf, name, operation, 0, 1.5707963267948966);
}
@Test
[commons-numbers] 05/30: Removed unused static asinh method and
constants.
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 a2cc4ec1cc601dec340a6d2c514f2c79d57c6670
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 11:35:46 2019 +0000
Removed unused static asinh method and constants.
---
.../apache/commons/numbers/complex/Complex.java | 78 +---------------------
1 file changed, 2 insertions(+), 76 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 010b40b..33aaec2 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
@@ -65,8 +65,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);
- /** 3*π/4. */
- private static final double PI_3_OVER_4 = 0.75 * Math.PI;
/** π/2. */
private static final double PI_OVER_2 = 0.5 * Math.PI;
/** π/4. */
@@ -1304,7 +1302,7 @@ public final class Complex implements Serializable {
* @param constructor Constructor.
* @return the inverse cosine of the complex number.
*/
- private static Complex acos(final double real, final double imaginary,
+ private static Complex acos(final double real, final double imaginary,
final ComplexConstructor constructor) {
// Compute with positive values and determine sign at the end
double x = Math.abs(real);
@@ -1472,7 +1470,7 @@ public final class Complex implements Serializable {
* @param constructor Constructor.
* @return the inverse sine of this complex number
*/
- private static Complex asin(final double real, final double imaginary,
+ private static Complex asin(final double real, final double imaginary,
final ComplexConstructor constructor) {
// Compute with positive values and determine sign at the end
double x = Math.abs(real);
@@ -1642,78 +1640,6 @@ public final class Complex implements Serializable {
}
/**
- * Compute the inverse hyperbolic sine of the complex number.
- *
- * <p>This function exists to allow implementation of the identity
- * {@code asin(z) = -i asinh(iz)}.<p>
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @param constructor Constructor.
- * @return the inverse hyperbolic sine of the complex number
- */
- private static Complex asinh(double real, double imaginary, ComplexConstructor constructor) {
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- // ISO C99: Preserve the equality
- // asinh(conj(z)) = conj(asinh(z))
- // and the odd function: f(z) = -f(-z)
- // by always computing on a positive valued Complex number.
- final double a = Math.abs(real);
- final double b = Math.abs(imaginary);
- // C99. G.7: Special case for imaginary only numbers
- if (a == 0 && b <= 1.0) {
- if (imaginary == 0) {
- return constructor.create(real, imaginary);
- }
- // asinh(iy) = i asin(y)
- final double im = Math.asin(imaginary);
- return constructor.create(real, im);
- }
- // square() is implemented using multiply
- final Complex z2 = multiply(a, b, a, b);
- // sqrt(1 + z^2)
- final Complex sqrt1pz2 = sqrt(1 + z2.real, z2.imaginary);
- // Compute the rest inline to avoid Complex object creation.
- // (x + y i) = z + sqrt(1 + z^2)
- final double x = a + sqrt1pz2.real;
- final double y = b + sqrt1pz2.imaginary;
- // (re + im i) = ln(z + sqrt(1 + z^2))
- final double re = Math.log(getAbsolute(x, y));
- final double im = getArgument(x, y);
- // Map back to the correct sign
- return constructor.create(changeSign(re, real),
- changeSign(im, imaginary));
- }
- if (Double.isInfinite(imaginary)) {
- return constructor.create(Math.copySign(Double.POSITIVE_INFINITY, real),
- Math.copySign(PI_OVER_2, imaginary));
- }
- // imaginary is NaN
- return NAN;
- }
- if (Double.isInfinite(real)) {
- if (Double.isFinite(imaginary)) {
- return constructor.create(real, Math.copySign(0, imaginary));
- }
- if (Double.isInfinite(imaginary)) {
- return constructor.create(real, Math.copySign(PI_OVER_4, imaginary));
- }
- // imaginary is NaN
- return constructor.create(real, Double.NaN);
- }
- // real is NaN
- if (imaginary == 0) {
- return constructor.create(Double.NaN, Math.copySign(0, imaginary));
- }
- if (Double.isInfinite(imaginary)) {
- return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
- }
- // optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
- return NAN;
- }
-
- /**
* Compute the
* <a href="http://mathworld.wolfram.com/InverseHyperbolicTangent.html">
* inverse hyperbolic tangent</a> of this complex number.
[commons-numbers] 03/30: Fix unspecified sign in C99 standard for
acos.
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 2c77d008b7232290a127993b41be6bd44becf894
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 11:31:13 2019 +0000
Fix unspecified sign in C99 standard for acos.
---
.../src/test/java/org/apache/commons/numbers/complex/CStandardTest.java | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
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 f8f150d..293b478 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
@@ -807,7 +807,7 @@ public class CStandardTest {
}
assertComplex(negInfInf, operation, threePiFourNegInf);
assertComplex(infInf, operation, piFourNegInf);
- assertComplex(infNaN, operation, nanInf);
+ assertComplex(infNaN, operation, nanInf, UnspecifiedSign.IMAGINARY);
assertComplex(negInfNaN, operation, nanNegInf, UnspecifiedSign.IMAGINARY);
for (double y : finite) {
assertComplex(complex(nan, y), operation, NAN);
[commons-numbers] 17/30: Remove unused/redundant methods.
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 e3351253466aab6ffc940e2dcd64cf7841ab5d0c
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Dec 19 19:08:08 2019 +0000
Remove unused/redundant methods.
Overflow safe average method is only used when one argument is less than
the other so compute
inline.
The getArgument method is not used.
---
.../apache/commons/numbers/complex/Complex.java | 32 ++--------------------
1 file changed, 2 insertions(+), 30 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 33b6056..b0583ff 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
@@ -2359,8 +2359,8 @@ public final class Complex implements Serializable {
// sqrt(2^exponent) = (2^exponent)^0.5 = 2^(exponent*0.5)
t = Math.scalb(Math.sqrt((scaleA + absC) / 2), exponent / 2);
} else {
- // Over-flow safe average
- t = Math.sqrt(average(absA, absC));
+ // Over-flow safe average: absA < absC and abdC is finite.
+ t = Math.sqrt(absA + (absC - absA) / 2);
}
if (real >= 0) {
@@ -2386,19 +2386,6 @@ public final class Complex implements Serializable {
}
/**
- * Compute the average of two positive finite values in an overflow safe manner.
- *
- * @param a the first value
- * @param b the second value
- * @return the average
- */
- private static double average(double a, double b) {
- return (a < b) ?
- a + (b - a) / 2 :
- b + (a - b) / 2;
- }
-
- /**
* Compute the
* <a href="http://mathworld.wolfram.com/Tangent.html">
* tangent</a> of this complex number.
@@ -2532,21 +2519,6 @@ public final class Complex implements Serializable {
}
/**
- * Compute the argument of the complex number.
- *
- * <p>This function exists for use in trigonomic functions.
- *
- * @param real Real part.
- * @param imaginary Imaginary part.
- * @return the argument.
- * @see Math#atan2(double, double)
- */
- private static double getArgument(double real, double imaginary) {
- // Delegate
- return Math.atan2(imaginary, real);
- }
-
- /**
* Computes the n-th roots of this complex number.
* The nth roots are defined by the formula:
* <pre>
[commons-numbers] 12/30: Added atanh assumptions test to check the
safe upper and lower limits.
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 5d9156bfe79858e804d91b3a2da3e9353b2f2506
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 17:31:59 2019 +0000
Added atanh assumptions test to check the safe upper and lower limits.
Updated the atanh code to use the boost assumptions.
---
.../apache/commons/numbers/complex/Complex.java | 25 +-
.../commons/numbers/complex/ComplexTest.java | 343 ++++++++++++---------
2 files changed, 198 insertions(+), 170 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 11f468f..684faa7 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
@@ -1719,13 +1719,12 @@ public final class Complex implements Serializable {
return NAN;
} else {
// x && y are finite or infinite.
- // Cases for very large finite are handled as if 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^2 -> inf
- // x <= SAFE_LOWER: x^2 -> 0
+ // x >= SAFE_UPPER: (1-x) == x
+ // x <= SAFE_LOWER: 1 - x^2 = 1
if ((x > SAFE_LOWER) && (x < SAFE_UPPER) && (y > SAFE_LOWER) && (y < SAFE_UPPER)) {
// Normal computation within a safe region.
@@ -1762,35 +1761,29 @@ public final class Complex implements Serializable {
// 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
+ // (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
- // This has been modified from the boost version to
- // include 1/(x*y) and -2/y. These are harmless if
- // machine precision prevents their addition to have an effect:
- // 1/(x*y) -> 0
- // (x-2) -> x
- re = Math.log1p((4 / y) / (1 / (x * y) + (x - 2) / y + y / x));
+ re = Math.log1p((4 / y) / (x / y + y / x));
} else if (y > 1) {
// Big x: divide through by x:
- // This has been modified from the boost version to
- // include 1/x and -2:
- re = Math.log1p(4 / (1 / x + x - 2 + y * y / 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 / (1 / x + x - 2));
+ 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) / (y + mxp1 * mxp1 / y));
+ re = Math.log1p((4 * x / y) / (mxp1 * mxp1 / y + y));
} else {
// Big y, small x, as above but neglect (1-x)^2/y:
// Note: The boost version has no log1p here.
- // This will tend towards 0 and log1p(0) = 0.
+ // This will tend towards 0 and log1p(0) = 0 so it may not matter.
re = Math.log1p(4 * x / y / y);
}
} else if (x != 1) {
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 ec395e2..d5f54cd 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
@@ -58,9 +58,7 @@ public class ComplexTest {
* Used to test the number category of a Complex.
*/
private enum NumberType {
- NAN,
- INFINITE,
- FINITE
+ NAN, INFINITE, FINITE
}
/**
@@ -121,7 +119,8 @@ public class ComplexTest {
assertNumberType(-inf, 0, NumberType.INFINITE);
assertNumberType(0, inf, NumberType.INFINITE);
assertNumberType(0, -inf, NumberType.INFINITE);
- // A complex or imaginary value with at least one infinite part is regarded as an infinity
+ // A complex or imaginary value with at least one infinite part is regarded as an
+ // infinity
// (even if its other part is a NaN).
assertNumberType(inf, nan, NumberType.INFINITE);
assertNumberType(-inf, nan, NumberType.INFINITE);
@@ -151,8 +150,8 @@ public class ComplexTest {
count += isInfinite ? 1 : 0;
count += isFinite ? 1 : 0;
Assertions.assertEquals(1, count,
- () -> String.format("Complex can be only one type: isNaN=%s, isInfinite=%s, isFinite=%s: %s",
- isNaN, isInfinite, isFinite, z));
+ () -> String.format("Complex can be only one type: isNaN=%s, isInfinite=%s, isFinite=%s: %s", isNaN,
+ isInfinite, isFinite, z));
switch (type) {
case FINITE:
Assertions.assertTrue(isFinite, () -> "not finite: " + z);
@@ -168,7 +167,6 @@ public class ComplexTest {
}
}
-
@Test
public void testProj() {
final Complex z = Complex.ofCartesian(3.0, 4.0);
@@ -759,8 +757,8 @@ public class ComplexTest {
final Complex z = c.multiply(Complex.I);
// Does not work when imaginary part is +0.0.
if (Double.compare(b, 0.0) == 0) {
- // (-0.0, 0.0).multiply( (0,1) ) => (-0.0, 0.0) expected (-0.0,-0.0)
- // ( 0.0, 0.0).multiply( (0,1) ) => ( 0.0, 0.0) expected (-0.0, 0.0)
+ // (-0.0, 0.0).multiply( (0,1) ) => (-0.0, 0.0) expected (-0.0,-0.0)
+ // ( 0.0, 0.0).multiply( (0,1) ) => ( 0.0, 0.0) expected (-0.0, 0.0)
Assertions.assertEquals(0, z.getReal(), 0.0);
Assertions.assertEquals(0, z.getImaginary(), 0.0);
Assertions.assertNotEquals(x, z);
@@ -788,16 +786,20 @@ public class ComplexTest {
final Complex z2 = c.multiply(Complex.I).negate();
// Does not work when imaginary part is -0.0.
if (Double.compare(b, -0.0) == 0) {
- // (-0.0,-0.0).multiply( (-0.0,-1) ) => ( 0.0, 0.0) expected (-0.0, 0.0)
- // ( 0.0,-0.0).multiply( (-0.0,-1) ) => (-0.0, 0.0) expected (-0.0,-0.0)
+ // (-0.0,-0.0).multiply( (-0.0,-1) ) => ( 0.0, 0.0) expected (-0.0,
+ // 0.0)
+ // ( 0.0,-0.0).multiply( (-0.0,-1) ) => (-0.0, 0.0) expected
+ // (-0.0,-0.0)
Assertions.assertEquals(0, z.getReal(), 0.0);
Assertions.assertEquals(0, z.getImaginary(), 0.0);
Assertions.assertNotEquals(x, z);
- // When multiply by I.negate() fails multiply by I then negate() works!
+ // When multiply by I.negate() fails multiply by I then negate()
+ // works!
Assertions.assertEquals(x, z2);
} else {
Assertions.assertEquals(x, z);
- // When multiply by I.negate() works multiply by I then negate() fails!
+ // When multiply by I.negate() works multiply by I then negate()
+ // fails!
Assertions.assertNotEquals(x, z2);
}
}
@@ -805,23 +807,23 @@ public class ComplexTest {
}
/**
- * Arithmetic test using combinations of +/- x for real, imaginary and
- * and the double argument for add, subtract, subtractFrom, multiply and divide,
- * where x is zero or non-zero.
+ * Arithmetic test using combinations of +/- x for real, imaginary and and the double
+ * argument for add, subtract, subtractFrom, multiply and divide, where x is zero or
+ * non-zero.
*
- * <p>The differences to the same argument as a Complex are tested. The only differences
- * should be the sign of zero in certain cases.
+ * <p>The differences to the same argument as a Complex are tested. The only
+ * differences should be the sign of zero in certain cases.
*/
@Test
public void testSignedArithmetic() {
// The following lists the conditions for the double primitive operation where
// the Complex operation is different. Here the double argument can be:
- // x : any value
- // +x : positive
+ // x : any value
+ // +x : positive
// +0.0: positive zero
- // -x : negative
+ // -x : negative
// -0.0: negative zero
- // 0 : any zero
+ // 0 : any zero
// use y for any non-zero value
// Check the known fail cases using an integer as a bit set.
@@ -839,17 +841,22 @@ public class ComplexTest {
// and the javadoc in Complex does not break down the actual cases.
// 16: (x,-0.0) + x
- assertSignedZeroArithmetic("addReal", Complex::add, ComplexTest::ofReal, Complex::add, 0b1111000000000000111100000000000011110000000000001111L);
+ assertSignedZeroArithmetic("addReal", Complex::add, ComplexTest::ofReal, Complex::add,
+ 0b1111000000000000111100000000000011110000000000001111L);
// 16: (-0.0,x) + x
- assertSignedZeroArithmetic("addImaginary", Complex::addImaginary, ComplexTest::ofImaginary, Complex::add, 0b1111111111111111L);
+ assertSignedZeroArithmetic("addImaginary", Complex::addImaginary, ComplexTest::ofImaginary, Complex::add,
+ 0b1111111111111111L);
// 0:
assertSignedZeroArithmetic("subtractReal", Complex::subtract, ComplexTest::ofReal, Complex::subtract, 0);
// 0:
- assertSignedZeroArithmetic("subtractImaginary", Complex::subtractImaginary, ComplexTest::ofImaginary, Complex::subtract, 0);
+ assertSignedZeroArithmetic("subtractImaginary", Complex::subtractImaginary, ComplexTest::ofImaginary,
+ Complex::subtract, 0);
// 16: x - (x,+0.0)
- assertSignedZeroArithmetic("subtractFromReal", Complex::subtractFrom, ComplexTest::ofReal, (y, z) -> z.subtract(y), 0b11110000000000001111000000000000111100000000000011110000L);
+ assertSignedZeroArithmetic("subtractFromReal", Complex::subtractFrom, ComplexTest::ofReal,
+ (y, z) -> z.subtract(y), 0b11110000000000001111000000000000111100000000000011110000L);
// 16: x - (+0.0,x)
- assertSignedZeroArithmetic("subtractFromImaginary", Complex::subtractFromImaginary, ComplexTest::ofImaginary, (y, z) -> z.subtract(y), 0b11111111111111110000000000000000L);
+ assertSignedZeroArithmetic("subtractFromImaginary", Complex::subtractFromImaginary, ComplexTest::ofImaginary,
+ (y, z) -> z.subtract(y), 0b11111111111111110000000000000000L);
// 4: (-0.0,-x) * +x
// 4: (+0.0,-0.0) * x
// 4: (+0.0,x) * -x
@@ -859,7 +866,8 @@ public class ComplexTest {
// 2: (+y,-x) * -0.0
// 2: (+x,-y) * +0.0
// 2: (+x,+y) * -0.0
- assertSignedZeroArithmetic("multiplyReal", Complex::multiply, ComplexTest::ofReal, Complex::multiply, 0b1001101011011000000100000001000010111010111110000101000001010L);
+ assertSignedZeroArithmetic("multiplyReal", Complex::multiply, ComplexTest::ofReal, Complex::multiply,
+ 0b1001101011011000000100000001000010111010111110000101000001010L);
// 4: (-0.0,+x) * +x
// 2: (+0.0,-0.0) * -x
// 4: (+0.0,+0.0) * x
@@ -867,23 +875,25 @@ public class ComplexTest {
// 2: (-y,+x) * +0.0
// 4: (+y,x) * -0.0
// 2: (+0.0,+/-y) * -/+0
- // 2: (+y,+/-0.0) * +/-y (sign 0.0 matches sign y)
+ // 2: (+y,+/-0.0) * +/-y (sign 0.0 matches sign y)
// 2: (+y,+x) * +0.0
- assertSignedZeroArithmetic("multiplyImaginary", Complex::multiplyImaginary, ComplexTest::ofImaginary, Complex::multiply, 0b11000110110101001000000010000001110001111101011010000010100000L);
+ assertSignedZeroArithmetic("multiplyImaginary", Complex::multiplyImaginary, ComplexTest::ofImaginary,
+ Complex::multiply, 0b11000110110101001000000010000001110001111101011010000010100000L);
// 2: (-0.0,0) / +y
// 2: (+0.0,+x) / -y
// 2: (-x,0) / -y
// 1: (-0.0,+y) / +y
// 1: (-y,+0.0) / -y
- assertSignedZeroArithmetic("divideReal", Complex::divide, ComplexTest::ofReal, Complex::divide, 0b100100001000000010000001000000011001000L);
+ assertSignedZeroArithmetic("divideReal", Complex::divide, ComplexTest::ofReal, Complex::divide,
+ 0b100100001000000010000001000000011001000L);
- // DivideImaginary has its own test as the result is not always equal ignoring the sign.
+ // DivideImaginary has its own test as the result is not always equal ignoring the
+ // sign.
}
- private static void assertSignedZeroArithmetic(String name,
- BiFunction<Complex, Double, Complex> doubleOperation,
- DoubleFunction<Complex> doubleToComplex,
- BiFunction<Complex, Complex, Complex> complexOperation, long expectedFailures) {
+ private static void assertSignedZeroArithmetic(String name, BiFunction<Complex, Double, Complex> doubleOperation,
+ DoubleFunction<Complex> doubleToComplex, BiFunction<Complex, Complex, Complex> complexOperation,
+ long expectedFailures) {
// With an operation on zero or non-zero arguments
final double[] arguments = {-0.0, 0.0, -2, 3};
for (final double a : arguments) {
@@ -896,21 +906,22 @@ public class ComplexTest {
expectedFailures >>>= 1;
// Check the same answer. Sign is allowed to be different for zero.
Assertions.assertEquals(y.getReal(), z.getReal(), 0, () -> c + " " + name + " " + arg + ": real");
- Assertions.assertEquals(y.getImaginary(), z.getImaginary(), 0, () -> c + " " + name + " " + arg + ": imaginary");
- Assertions.assertEquals(expectedFailure, !y.equals(z), () -> c + " " + name + " " + arg + ": sign-difference");
+ Assertions.assertEquals(y.getImaginary(), z.getImaginary(), 0,
+ () -> c + " " + name + " " + arg + ": imaginary");
+ Assertions.assertEquals(expectedFailure, !y.equals(z),
+ () -> c + " " + name + " " + arg + ": sign-difference");
}
}
}
}
/**
- * Arithmetic test using combinations of +/- x for real, imaginary and
- * and the double argument for divideImaginary,
- * where x is zero or non-zero.
+ * Arithmetic test using combinations of +/- x for real, imaginary and and the double
+ * argument for divideImaginary, where x is zero or non-zero.
*
- * <p>The differences to the same argument as a Complex are tested. This checks for sign
- * differences of zero or, if divide by zero, that the result is equal
- * to divide by zero using a Complex then multiplied by I.
+ * <p>The differences to the same argument as a Complex are tested. This checks for
+ * sign differences of zero or, if divide by zero, that the result is equal to divide
+ * by zero using a Complex then multiplied by I.
*/
@Test
public void testDivideImaginaryArithmetic() {
@@ -923,7 +934,8 @@ public class ComplexTest {
// 2: (+0.0,+/-y) / +0.0
// 4: (-y,x) / +0.0
// 4: (y,x) / +0.0
- // If multiplied by -I all the divide by -0.0 cases have sign errors and / +0.0 is OK.
+ // If multiplied by -I all the divide by -0.0 cases have sign errors and / +0.0 is
+ // OK.
long expectedFailures = 0b11001101111011001100110011001110110011110010000111001101000000L;
// With an operation on zero or non-zero arguments
final double[] arguments = {-0.0, 0.0, -2, 3};
@@ -935,7 +947,8 @@ public class ComplexTest {
Complex z = c.divide(ofImaginary(arg));
final boolean expectedFailure = (expectedFailures & 0x1) == 1;
expectedFailures >>>= 1;
- // If divide by zero then the divide(Complex) method matches divide by real.
+ // If divide by zero then the divide(Complex) method matches divide by
+ // real.
// To match divide by imaginary requires multiplication by I.
if (arg == 0) {
// Same result if multiplied by I. The sign may not match so
@@ -948,10 +961,14 @@ public class ComplexTest {
Assertions.assertEquals(ya, za, () -> c + " divideImaginary " + arg + ": real");
Assertions.assertEquals(yb, zb, () -> c + " divideImaginary " + arg + ": imaginary");
} else {
- // Check the same answer. Sign is allowed to be different for zero.
- Assertions.assertEquals(y.getReal(), z.getReal(), 0, () -> c + " divideImaginary " + arg + ": real");
- Assertions.assertEquals(y.getImaginary(), z.getImaginary(), 0, () -> c + " divideImaginary " + arg + ": imaginary");
- Assertions.assertEquals(expectedFailure, !y.equals(z), () -> c + " divideImaginary " + arg + ": sign-difference");
+ // Check the same answer. Sign is allowed to be different for
+ // zero.
+ Assertions.assertEquals(y.getReal(), z.getReal(), 0,
+ () -> c + " divideImaginary " + arg + ": real");
+ Assertions.assertEquals(y.getImaginary(), z.getImaginary(), 0,
+ () -> c + " divideImaginary " + arg + ": imaginary");
+ Assertions.assertEquals(expectedFailure, !y.equals(z),
+ () -> c + " divideImaginary " + arg + ": sign-difference");
}
}
}
@@ -1175,15 +1192,13 @@ public class ComplexTest {
@Test
public void testFloatingPointEqualsPrecondition1() {
Assertions.assertThrows(NullPointerException.class,
- () -> Complex.equals(Complex.ofCartesian(3.0, 4.0), null, 3)
- );
+ () -> Complex.equals(Complex.ofCartesian(3.0, 4.0), null, 3));
}
@Test
public void testFloatingPointEqualsPrecondition2() {
Assertions.assertThrows(NullPointerException.class,
- () -> Complex.equals(null, Complex.ofCartesian(3.0, 4.0), 3)
- );
+ () -> Complex.equals(null, Complex.ofCartesian(3.0, 4.0), 3));
}
@Test
@@ -1304,28 +1319,23 @@ public class ComplexTest {
/**
* Test {@link Complex#equals(Object)}. It should be consistent with
- * {@link Arrays#equals(double[], double[])} called using the components of two complex numbers.
+ * {@link Arrays#equals(double[], double[])} called using the components of two
+ * complex numbers.
*/
@Test
public void testEqualsIsConsistentWithArraysEquals() {
// Explicit check of the cases documented in the Javadoc:
- assertEqualsIsConsistentWithArraysEquals(
- Complex.ofCartesian(Double.NaN, 0.0),
- Complex.ofCartesian(Double.NaN, 1.0), "NaN real and different non-NaN imaginary");
- assertEqualsIsConsistentWithArraysEquals(
- Complex.ofCartesian(0.0, Double.NaN),
- Complex.ofCartesian(1.0, Double.NaN), "Different non-NaN real and NaN imaginary");
- assertEqualsIsConsistentWithArraysEquals(
- Complex.ofCartesian(0.0, 0.0),
- Complex.ofCartesian(-0.0, 0.0), "Different real zeros");
- assertEqualsIsConsistentWithArraysEquals(
- Complex.ofCartesian(0.0, 0.0),
- Complex.ofCartesian(0.0, -0.0), "Different imaginary zeros");
+ assertEqualsIsConsistentWithArraysEquals(Complex.ofCartesian(Double.NaN, 0.0),
+ Complex.ofCartesian(Double.NaN, 1.0), "NaN real and different non-NaN imaginary");
+ assertEqualsIsConsistentWithArraysEquals(Complex.ofCartesian(0.0, Double.NaN),
+ Complex.ofCartesian(1.0, Double.NaN), "Different non-NaN real and NaN imaginary");
+ assertEqualsIsConsistentWithArraysEquals(Complex.ofCartesian(0.0, 0.0), Complex.ofCartesian(-0.0, 0.0),
+ "Different real zeros");
+ assertEqualsIsConsistentWithArraysEquals(Complex.ofCartesian(0.0, 0.0), Complex.ofCartesian(0.0, -0.0),
+ "Different imaginary zeros");
// Test some values of edge cases
- final double[] values = {
- Double.NaN, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, -1, 0, 1
- };
+ final double[] values = {Double.NaN, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, -1, 0, 1};
final ArrayList<Complex> list = createCombinations(values);
for (final Complex c : list) {
@@ -1353,9 +1363,7 @@ public class ComplexTest {
@Test
public void testEqualsWithDifferentNaNs() {
// Test some NaN combinations
- final double[] values = {
- Double.NaN, 0, 1
- };
+ final double[] values = {Double.NaN, 0, 1};
final ArrayList<Complex> list = createCombinations(values);
// Is the all-vs-all comparison only the exact same values should be equal, e.g.
@@ -1374,34 +1382,32 @@ public class ComplexTest {
}
/**
- * Test the two complex numbers with {@link Complex#equals(Object)} and check the result
- * is consistent with {@link Arrays#equals(double[], double[])}.
+ * Test the two complex numbers with {@link Complex#equals(Object)} and check the
+ * result is consistent with {@link Arrays#equals(double[], double[])}.
*
* @param c1 the first complex
* @param c2 the second complex
* @param msg the message to append to an assertion error
*/
private static void assertEqualsIsConsistentWithArraysEquals(Complex c1, Complex c2, String msg) {
- final boolean expected = Arrays.equals(new double[]{c1.getReal(), c1.getImaginary()},
- new double[]{c2.getReal(), c2.getImaginary()});
+ final boolean expected = Arrays.equals(new double[] {c1.getReal(), c1.getImaginary()},
+ new double[] {c2.getReal(), c2.getImaginary()});
final boolean actual = c1.equals(c2);
- Assertions.assertEquals(expected, actual, () -> String.format(
- "equals(Object) is not consistent with Arrays.equals: %s. %s vs %s", msg, c1, c2));
+ Assertions.assertEquals(expected, actual,
+ () -> String.format("equals(Object) is not consistent with Arrays.equals: %s. %s vs %s", msg, c1, c2));
}
/**
* Test {@link Complex#hashCode()}. It should be consistent with
* {@link Arrays#hashCode(double[])} called using the components of the complex number
- * and fulfil the contract of {@link Object#hashCode()},
- * i.e. objects with different hash codes are {@code false} for {@link Object#equals(Object)}.
+ * and fulfil the contract of {@link Object#hashCode()}, i.e. objects with different
+ * hash codes are {@code false} for {@link Object#equals(Object)}.
*/
@Test
public void testHashCode() {
// Test some values match Arrays.hashCode(double[])
- final double[] values = {
- Double.NaN, Double.NEGATIVE_INFINITY, -3.45, -1, -0.0, 0.0,
- Double.MIN_VALUE, 1, 3.45, Double.POSITIVE_INFINITY
- };
+ final double[] values = {Double.NaN, Double.NEGATIVE_INFINITY, -3.45, -1, -0.0, 0.0, Double.MIN_VALUE, 1, 3.45,
+ Double.POSITIVE_INFINITY};
final ArrayList<Complex> list = createCombinations(values);
final String msg = "'equals' not compatible with 'hashCode'";
@@ -1413,7 +1419,8 @@ public class ComplexTest {
final int hash = c.hashCode();
Assertions.assertEquals(expected, hash, "hashCode does not match Arrays.hashCode({re, im})");
- // Test a copy has the same hash code, i.e. is not System.identityHashCode(Object)
+ // Test a copy has the same hash code, i.e. is not
+ // System.identityHashCode(Object)
final Complex copy = Complex.ofCartesian(real, imag);
Assertions.assertEquals(hash, copy.hashCode(), "Copy hash code is not equal");
@@ -1441,11 +1448,13 @@ public class ComplexTest {
/**
* Specific test that different representations of zero satisfy the contract of
* {@link Object#hashCode()}: if two objects have different hash codes, "equals" must
- * return false. This is an issue with using {@link Double#hashCode(double)} to create hash
- * codes and {@code ==} for equality when using different representations of zero:
- * Double.hashCode(-0.0) != Double.hashCode(0.0) but -0.0 == 0.0 is {@code true}.
+ * return false. This is an issue with using {@link Double#hashCode(double)} to create
+ * hash codes and {@code ==} for equality when using different representations of
+ * zero: Double.hashCode(-0.0) != Double.hashCode(0.0) but -0.0 == 0.0 is
+ * {@code true}.
*
- * @see <a href="https://issues.apache.org/jira/projects/MATH/issues/MATH-1118">MATH-1118</a>
+ * @see <a
+ * href="https://issues.apache.org/jira/projects/MATH/issues/MATH-1118">MATH-1118</a>
*/
@Test
public void testHashCodeWithDifferentZeros() {
@@ -1483,9 +1492,10 @@ public class ComplexTest {
}
/**
- * Perform the smallest change to the value. This returns the next double value adjacent to
- * d in the direction of infinity. Edge cases: if already infinity then return the next closest
- * in the direction of negative infinity; if nan then return 0.
+ * Perform the smallest change to the value. This returns the next double value
+ * adjacent to d in the direction of infinity. Edge cases: if already infinity then
+ * return the next closest in the direction of negative infinity; if nan then return
+ * 0.
*
* @param x the x
* @return the new value
@@ -1494,9 +1504,7 @@ public class ComplexTest {
if (Double.isNaN(x)) {
return 0;
}
- return x == Double.POSITIVE_INFINITY ?
- Math.nextDown(x) :
- Math.nextUp(x);
+ return x == Double.POSITIVE_INFINITY ? Math.nextDown(x) : Math.nextUp(x);
}
@Test
@@ -1506,22 +1514,23 @@ public class ComplexTest {
// CHECKSTYLE: stop Regexp
System.out.println(">>testJava()");
// MathTest#testExpSpecialCases() checks the following:
- // Assert.assertEquals("exp of -infinity should be 0.0", 0.0, Math.exp(Double.NEGATIVE_INFINITY), Precision.EPSILON);
+ // Assert.assertEquals("exp of -infinity should be 0.0", 0.0,
+ // Math.exp(Double.NEGATIVE_INFINITY), Precision.EPSILON);
// Let's check how well Math works:
System.out.println("Math.exp=" + Math.exp(Double.NEGATIVE_INFINITY));
- final String[] props = {
- "java.version", // Java Runtime Environment version
+ final String[] props = {"java.version", // Java Runtime Environment version
"java.vendor", // Java Runtime Environment vendor
- "java.vm.specification.version", // Java Virtual Machine specification version
- "java.vm.specification.vendor", // Java Virtual Machine specification vendor
- "java.vm.specification.name", // Java Virtual Machine specification name
+ "java.vm.specification.version", // Java Virtual Machine specification version
+ "java.vm.specification.vendor", // Java Virtual Machine specification vendor
+ "java.vm.specification.name", // Java Virtual Machine specification name
"java.vm.version", // Java Virtual Machine implementation version
- "java.vm.vendor", // Java Virtual Machine implementation vendor
- "java.vm.name", // Java Virtual Machine implementation name
- "java.specification.version", // Java Runtime Environment specification version
- "java.specification.vendor", // Java Runtime Environment specification vendor
+ "java.vm.vendor", // Java Virtual Machine implementation vendor
+ "java.vm.name", // Java Virtual Machine implementation name
+ "java.specification.version", // Java Runtime Environment specification
+ // version
+ "java.specification.vendor", // Java Runtime Environment specification vendor
"java.specification.name", // Java Runtime Environment specification name
- "java.class.version", // Java class format version number
+ "java.class.version", // Java class format version number
};
for (final String t : props) {
System.out.println(t + "=" + System.getProperty(t));
@@ -1627,6 +1636,7 @@ public class ComplexTest {
/**
* Test: computing <b>third roots</b> of z.
+ *
* <pre>
* <code>
* <b>z = -2 + 2 * i</b>
@@ -1645,18 +1655,19 @@ public class ComplexTest {
// Returned Collection must not be empty!
Assertions.assertEquals(3, thirdRootsOfZ.length);
// test z_0
- Assertions.assertEquals(1.0, thirdRootsOfZ[0].getReal(), 1.0e-5);
- Assertions.assertEquals(1.0, thirdRootsOfZ[0].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(1.0, thirdRootsOfZ[0].getReal(), 1.0e-5);
+ Assertions.assertEquals(1.0, thirdRootsOfZ[0].getImaginary(), 1.0e-5);
// test z_1
- Assertions.assertEquals(-1.3660254037844386, thirdRootsOfZ[1].getReal(), 1.0e-5);
- Assertions.assertEquals(0.36602540378443843, thirdRootsOfZ[1].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(-1.3660254037844386, thirdRootsOfZ[1].getReal(), 1.0e-5);
+ Assertions.assertEquals(0.36602540378443843, thirdRootsOfZ[1].getImaginary(), 1.0e-5);
// test z_2
- Assertions.assertEquals(0.366025403784439, thirdRootsOfZ[2].getReal(), 1.0e-5);
- Assertions.assertEquals(-1.3660254037844384, thirdRootsOfZ[2].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(0.366025403784439, thirdRootsOfZ[2].getReal(), 1.0e-5);
+ Assertions.assertEquals(-1.3660254037844384, thirdRootsOfZ[2].getImaginary(), 1.0e-5);
}
/**
* Test: computing <b>fourth roots</b> of z.
+ *
* <pre>
* <code>
* <b>z = 5 - 2 * i</b>
@@ -1676,21 +1687,22 @@ public class ComplexTest {
// Returned Collection must not be empty!
Assertions.assertEquals(4, fourthRootsOfZ.length);
// test z_0
- Assertions.assertEquals(1.5164629308487783, fourthRootsOfZ[0].getReal(), 1.0e-5);
- Assertions.assertEquals(-0.14469266210702247, fourthRootsOfZ[0].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(1.5164629308487783, fourthRootsOfZ[0].getReal(), 1.0e-5);
+ Assertions.assertEquals(-0.14469266210702247, fourthRootsOfZ[0].getImaginary(), 1.0e-5);
// test z_1
- Assertions.assertEquals(0.14469266210702256, fourthRootsOfZ[1].getReal(), 1.0e-5);
- Assertions.assertEquals(1.5164629308487783, fourthRootsOfZ[1].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(0.14469266210702256, fourthRootsOfZ[1].getReal(), 1.0e-5);
+ Assertions.assertEquals(1.5164629308487783, fourthRootsOfZ[1].getImaginary(), 1.0e-5);
// test z_2
- Assertions.assertEquals(-1.5164629308487783, fourthRootsOfZ[2].getReal(), 1.0e-5);
- Assertions.assertEquals(0.14469266210702267, fourthRootsOfZ[2].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(-1.5164629308487783, fourthRootsOfZ[2].getReal(), 1.0e-5);
+ Assertions.assertEquals(0.14469266210702267, fourthRootsOfZ[2].getImaginary(), 1.0e-5);
// test z_3
- Assertions.assertEquals(-0.14469266210702275, fourthRootsOfZ[3].getReal(), 1.0e-5);
- Assertions.assertEquals(-1.5164629308487783, fourthRootsOfZ[3].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(-0.14469266210702275, fourthRootsOfZ[3].getReal(), 1.0e-5);
+ Assertions.assertEquals(-1.5164629308487783, fourthRootsOfZ[3].getImaginary(), 1.0e-5);
}
/**
* Test: computing <b>third roots</b> of z.
+ *
* <pre>
* <code>
* <b>z = 8</b>
@@ -1710,19 +1722,19 @@ public class ComplexTest {
// Returned Collection must not be empty!
Assertions.assertEquals(3, thirdRootsOfZ.length);
// test z_0
- Assertions.assertEquals(2.0, thirdRootsOfZ[0].getReal(), 1.0e-5);
- Assertions.assertEquals(0.0, thirdRootsOfZ[0].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(2.0, thirdRootsOfZ[0].getReal(), 1.0e-5);
+ Assertions.assertEquals(0.0, thirdRootsOfZ[0].getImaginary(), 1.0e-5);
// test z_1
- Assertions.assertEquals(-1.0, thirdRootsOfZ[1].getReal(), 1.0e-5);
+ Assertions.assertEquals(-1.0, thirdRootsOfZ[1].getReal(), 1.0e-5);
Assertions.assertEquals(1.7320508075688774, thirdRootsOfZ[1].getImaginary(), 1.0e-5);
// test z_2
- Assertions.assertEquals(-1.0, thirdRootsOfZ[2].getReal(), 1.0e-5);
+ Assertions.assertEquals(-1.0, thirdRootsOfZ[2].getReal(), 1.0e-5);
Assertions.assertEquals(-1.732050807568877, thirdRootsOfZ[2].getImaginary(), 1.0e-5);
}
-
/**
* Test: computing <b>third roots</b> of z with real part 0.
+ *
* <pre>
* <code>
* <b>z = 2 * i</b>
@@ -1741,21 +1753,21 @@ public class ComplexTest {
// Returned Collection must not be empty!
Assertions.assertEquals(3, thirdRootsOfZ.length);
// test z_0
- Assertions.assertEquals(1.0911236359717216, thirdRootsOfZ[0].getReal(), 1.0e-5);
- Assertions.assertEquals(0.6299605249474365, thirdRootsOfZ[0].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(1.0911236359717216, thirdRootsOfZ[0].getReal(), 1.0e-5);
+ Assertions.assertEquals(0.6299605249474365, thirdRootsOfZ[0].getImaginary(), 1.0e-5);
// test z_1
- Assertions.assertEquals(-1.0911236359717216, thirdRootsOfZ[1].getReal(), 1.0e-5);
- Assertions.assertEquals(0.6299605249474365, thirdRootsOfZ[1].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(-1.0911236359717216, thirdRootsOfZ[1].getReal(), 1.0e-5);
+ Assertions.assertEquals(0.6299605249474365, thirdRootsOfZ[1].getImaginary(), 1.0e-5);
// test z_2
- Assertions.assertEquals(-2.3144374213981936E-16, thirdRootsOfZ[2].getReal(), 1.0e-5);
- Assertions.assertEquals(-1.2599210498948732, thirdRootsOfZ[2].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(-2.3144374213981936E-16, thirdRootsOfZ[2].getReal(), 1.0e-5);
+ Assertions.assertEquals(-1.2599210498948732, thirdRootsOfZ[2].getImaginary(), 1.0e-5);
}
/**
- * Test: compute <b>third roots</b> using a negative argument
- * to go clockwise around the unit circle. Fourth roots of one
- * are taken in both directions around the circle using
- * positive and negative arguments.
+ * Test: compute <b>third roots</b> using a negative argument to go clockwise around
+ * the unit circle. Fourth roots of one are taken in both directions around the circle
+ * using positive and negative arguments.
+ *
* <pre>
* <code>
* <b>z = 1</b>
@@ -1773,31 +1785,31 @@ public class ComplexTest {
// The List holding all fourth roots
Complex[] fourthRootsOfZ = z.nthRoot(4).toArray(new Complex[0]);
// test z_0
- Assertions.assertEquals(1, fourthRootsOfZ[0].getReal(), 1.0e-5);
- Assertions.assertEquals(0, fourthRootsOfZ[0].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(1, fourthRootsOfZ[0].getReal(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[0].getImaginary(), 1.0e-5);
// test z_1
- Assertions.assertEquals(0, fourthRootsOfZ[1].getReal(), 1.0e-5);
- Assertions.assertEquals(1, fourthRootsOfZ[1].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[1].getReal(), 1.0e-5);
+ Assertions.assertEquals(1, fourthRootsOfZ[1].getImaginary(), 1.0e-5);
// test z_2
- Assertions.assertEquals(-1, fourthRootsOfZ[2].getReal(), 1.0e-5);
- Assertions.assertEquals(0, fourthRootsOfZ[2].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(-1, fourthRootsOfZ[2].getReal(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[2].getImaginary(), 1.0e-5);
// test z_3
- Assertions.assertEquals(0, fourthRootsOfZ[3].getReal(), 1.0e-5);
- Assertions.assertEquals(-1, fourthRootsOfZ[3].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[3].getReal(), 1.0e-5);
+ Assertions.assertEquals(-1, fourthRootsOfZ[3].getImaginary(), 1.0e-5);
// go clockwise around the unit circle using negative argument
fourthRootsOfZ = z.nthRoot(-4).toArray(new Complex[0]);
// test z_0
- Assertions.assertEquals(1, fourthRootsOfZ[0].getReal(), 1.0e-5);
- Assertions.assertEquals(0, fourthRootsOfZ[0].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(1, fourthRootsOfZ[0].getReal(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[0].getImaginary(), 1.0e-5);
// test z_1
- Assertions.assertEquals(0, fourthRootsOfZ[1].getReal(), 1.0e-5);
- Assertions.assertEquals(-1, fourthRootsOfZ[1].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[1].getReal(), 1.0e-5);
+ Assertions.assertEquals(-1, fourthRootsOfZ[1].getImaginary(), 1.0e-5);
// test z_2
- Assertions.assertEquals(-1, fourthRootsOfZ[2].getReal(), 1.0e-5);
- Assertions.assertEquals(0, fourthRootsOfZ[2].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(-1, fourthRootsOfZ[2].getReal(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[2].getImaginary(), 1.0e-5);
// test z_3
- Assertions.assertEquals(0, fourthRootsOfZ[3].getReal(), 1.0e-5);
- Assertions.assertEquals(1, fourthRootsOfZ[3].getImaginary(), 1.0e-5);
+ Assertions.assertEquals(0, fourthRootsOfZ[3].getReal(), 1.0e-5);
+ Assertions.assertEquals(1, fourthRootsOfZ[3].getImaginary(), 1.0e-5);
}
@Test
@@ -1811,6 +1823,7 @@ public class ComplexTest {
Assertions.assertTrue(Double.isNaN(c.getImaginary()));
}
}
+
@Test
public void testNthRootInf() {
final int n = 3;
@@ -1882,8 +1895,8 @@ public class ComplexTest {
@Test
public void testParse() {
- final double[] parts = {Double.NEGATIVE_INFINITY, -1, -0.0, 0.0, 1, Math.PI,
- Double.POSITIVE_INFINITY, Double.NaN};
+ final double[] parts = {Double.NEGATIVE_INFINITY, -1, -0.0, 0.0, 1, Math.PI, Double.POSITIVE_INFINITY,
+ Double.NaN};
for (final double x : parts) {
for (final double y : parts) {
final Complex z = Complex.ofCartesian(x, y);
@@ -2004,4 +2017,26 @@ public class ComplexTest {
Assertions.assertEquals(0.54930614433405489, c.getReal());
Assertions.assertEquals(1.5707963267948966, c.getImaginary());
}
+
+ @Test
+ public void testAtanhAssumptions() {
+ // Compute the same constants used by atanh
+ final double safeUpper = Math.sqrt(Double.MAX_VALUE) / 2;
+ final double safeLower = Math.sqrt(Double.MIN_NORMAL) * 2;
+
+ // Can we assume (1+x) = x when x is large
+ Assertions.assertEquals(safeUpper, 1 + safeUpper);
+ // Can we assume (1-x) = -x when x is large
+ Assertions.assertEquals(-safeUpper, 1 - safeUpper);
+ // Can we assume (y^2/x) = 0 when y is small and x is large
+ Assertions.assertEquals(0, safeLower * safeLower / safeUpper);
+ // Can we assume (1-x)^2/y + y = y when x <= 1. Try with x = 0.
+ Assertions.assertEquals(safeUpper, 1 / safeUpper + safeUpper);
+ // Can we assume (4+y^2) = 4 when y is small
+ Assertions.assertEquals(4, 4 + safeLower * safeLower);
+ // Can we assume (1-x)^2 = 1 when x is small
+ Assertions.assertEquals(1, (1 - safeLower) * (1 - safeLower));
+ // Can we assume 1 - y^2 = 1 when y is small
+ Assertions.assertEquals(1, 1 - safeLower * safeLower);
+ }
}
[commons-numbers] 07/30: Update the max ULP for atanh from 17 to 1
for the C reference test.
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 26b9ae7953d0c50faf3946c6aaf5bfe27227edf2
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 16:22:19 2019 +0000
Update the max ULP for atanh from 17 to 1 for the C reference test.
---
.../test/java/org/apache/commons/numbers/complex/CReferenceTest.java | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
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 f26611c..38090c1 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
@@ -284,7 +284,7 @@ public class CReferenceTest {
@Test
public void testAtanh() {
// Odd function: negative real cases defined by positive real cases
- assertOperation("atanh", Complex::atanh, 17);
+ assertOperation("atanh", Complex::atanh, 1);
}
@Test
[commons-numbers] 24/30: Fix overflow in tanh. Added edge case
tests.
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 3366f35d8120ddc5c52f991c1d120a7c7f95c72f
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 12:53:06 2019 +0000
Fix overflow in tanh. Added edge case tests.
Overflow is cause by computation of Math.sin(2a) or Math.cos(2a) when a
is finite and 2a is infinite. This has been fixed to use trigonomic
identities for sin(2a) and cos(2a) in this case.
---
.../apache/commons/numbers/complex/Complex.java | 69 ++++++++++++++++++----
.../numbers/complex/ComplexEdgeCaseTest.java | 30 +++++++++-
2 files changed, 86 insertions(+), 13 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 c9be094..4516be4 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
@@ -2400,6 +2400,7 @@ public final class Complex implements Serializable {
* </pre>
*
* @return the tangent of {@code this}.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tan/">Tangent</a>
*/
public Complex tan() {
// Define in terms of tanh
@@ -2420,6 +2421,7 @@ public final class Complex implements Serializable {
* <p>This is an odd function: {@code f(z) = -f(-z)}.
*
* @return the hyperbolic tangent of {@code this}.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Tanh/">Tanh</a>
*/
public Complex tanh() {
return tanh(real, imaginary, Complex::ofCartesian);
@@ -2441,10 +2443,8 @@ public final class Complex implements Serializable {
// Perform edge-condition checks on twice the imaginary value.
// This handles very big imaginary numbers as infinite.
- final double imaginary2 = 2 * imaginary;
-
if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary2)) {
+ if (Double.isFinite(imaginary)) {
if (real == 0) {
// Identity: sin x / (1 + cos x) = tan(x/2)
return constructor.create(real, Math.tan(imaginary));
@@ -2458,26 +2458,26 @@ public final class Complex implements Serializable {
// Math.cosh returns positive infinity for infinity.
// cosh -> inf
- final double d = Math.cosh(real2) + Math.cos(imaginary2);
+ final double divisor = Math.cosh(real2) + cos2(imaginary);
// Math.sinh returns the input infinity for infinity.
// sinh -> inf for positive x; else -inf
final double sinhRe2 = Math.sinh(real2);
// Avoid inf / inf
- if (Double.isInfinite(sinhRe2) && Double.isInfinite(d)) {
- // Fall-through to the result if infinite
- return constructor.create(Math.copySign(1, real), Math.copySign(0, Math.sin(imaginary2)));
+ if (Double.isInfinite(sinhRe2) && Double.isInfinite(divisor)) {
+ // Handle as if real was infinite
+ return constructor.create(Math.copySign(1, real), Math.copySign(0, sin2(imaginary)));
}
- return constructor.create(sinhRe2 / d,
- Math.sin(imaginary2) / d);
+ return constructor.create(sinhRe2 / divisor,
+ sin2(imaginary) / divisor);
}
// imaginary is infinite or NaN
return NAN;
}
if (Double.isInfinite(real)) {
- if (Double.isFinite(imaginary2)) {
- return constructor.create(Math.copySign(1, real), Math.copySign(0, Math.sin(imaginary2)));
+ if (Double.isFinite(imaginary)) {
+ return constructor.create(Math.copySign(1, real), Math.copySign(0, sin2(imaginary)));
}
// imaginary is infinite or NaN
return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
@@ -2490,7 +2490,52 @@ public final class Complex implements Serializable {
return NAN;
}
- /**
+ /**
+ * Safely compute {@code cos(2*a)} when {@code a} is finite.
+ * Note that {@link Math#cos(double)} returns NaN when the input is infinite.
+ * If {@code 2*a} is finite use {@code Math.cos(2*a)}; otherwise use the identity:
+ * <pre>
+ * <code>
+ * cos(2a) = 2 cos<sup>2</sup>(a) - 1
+ * </code>
+ * </pre>
+ *
+ * @param a Angle a.
+ * @return the cosine of 2a
+ * @see Math#cos(double)
+ */
+ private static double cos2(double a) {
+ final double twoA = 2 * a;
+ if (Double.isFinite(twoA)) {
+ return Math.cos(twoA);
+ }
+ final double cosA = Math.cos(a);
+ return 2 * cosA * cosA - 1;
+ }
+
+ /**
+ * Safely compute {@code sin(2*a)} when {@code a} is finite.
+ * Note that {@link Math#sin(double)} returns NaN when the input is infinite.
+ * If {@code 2*a} is finite use {@code Math.sin(2*a)}; otherwise use the identity:
+ * <pre>
+ * <code>
+ * sin(2a) = 2 sin(a) cos(a)
+ * </code>
+ * </pre>
+ *
+ * @param a Angle a.
+ * @return the sine of 2a
+ * @see Math#sin(double)
+ */
+ private static double sin2(double a) {
+ final double twoA = 2 * a;
+ if (Double.isFinite(twoA)) {
+ return Math.sin(twoA);
+ }
+ return 2 * Math.sin(a) * Math.cos(a);
+ }
+
+ /**
* Compute the argument of this complex number.
*
* <p>The argument is the angle phi between the positive real axis and
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 fd163e8..0fb94af 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
@@ -243,7 +243,35 @@ public class ComplexEdgeCaseTest {
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
+ final String name = "tanh";
+ final UnaryOperator<Complex> operation = Complex::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);
+
+ // 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);
+
+ // 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);
+
+ // 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);
}
@Test
[commons-numbers] 26/30: Edge cases for cosh and sinh.
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 ba7fe3aa925744d4fba5eacf3758b6f8625ee715
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 13:21:03 2019 +0000
Edge cases for cosh and sinh.
---
.../numbers/complex/ComplexEdgeCaseTest.java | 36 ++++++++++++++++++++--
1 file changed, 34 insertions(+), 2 deletions(-)
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 0fb94af..bee0aff 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
@@ -229,14 +229,46 @@ public class ComplexEdgeCaseTest {
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
+ final String name = "cosh";
+ final UnaryOperator<Complex> operation = Complex::cosh;
+
+ // Implementation defers to java.util.Math.
+ // Hit edge cases for extreme values.
+ double big = Double.MAX_VALUE;
+ double medium = 2;
+ 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);
}
@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
+ final String name = "sinh";
+ final UnaryOperator<Complex> operation = Complex::sinh;
+
+ // Implementation defers to java.util.Math.
+ // Hit edge cases for extreme values.
+ double big = Double.MAX_VALUE;
+ double medium = 2;
+ 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);
}
@Test
[commons-numbers] 14/30: Trailing whitespace.
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 913b01293b9aba4058dee953422c03c3c9689c30
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 17:34:14 2019 +0000
Trailing whitespace.
---
.../java/org/apache/commons/numbers/complex/ComplexTest.java | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
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 d4189b5..28601fb 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
@@ -1636,7 +1636,7 @@ public class ComplexTest {
/**
* Test: computing <b>third roots</b> of z.
- *
+ *
* <pre>
* <code>
* <b>z = -2 + 2 * i</b>
@@ -1667,7 +1667,7 @@ public class ComplexTest {
/**
* Test: computing <b>fourth roots</b> of z.
- *
+ *
* <pre>
* <code>
* <b>z = 5 - 2 * i</b>
@@ -1702,7 +1702,7 @@ public class ComplexTest {
/**
* Test: computing <b>third roots</b> of z.
- *
+ *
* <pre>
* <code>
* <b>z = 8</b>
@@ -1734,7 +1734,7 @@ public class ComplexTest {
/**
* Test: computing <b>third roots</b> of z with real part 0.
- *
+ *
* <pre>
* <code>
* <b>z = 2 * i</b>
@@ -1767,7 +1767,7 @@ public class ComplexTest {
* Test: compute <b>third roots</b> using a negative argument to go clockwise around
* the unit circle. Fourth roots of one are taken in both directions around the circle
* using positive and negative arguments.
- *
+ *
* <pre>
* <code>
* <b>z = 1</b>
[commons-numbers] 22/30: Added edge case tests for asin.
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 e1ad2a0e9c6f05b38d26d6484c3a3509b6a5935a
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Dec 19 23:09:58 2019 +0000
Added edge case tests for asin.
---
.../numbers/complex/ComplexEdgeCaseTest.java | 40 +++++++++++++++++-----
1 file changed, 31 insertions(+), 9 deletions(-)
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 42ca2cb..fd163e8 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
@@ -161,19 +161,41 @@ public class ComplexEdgeCaseTest {
assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation, 2.4252018043912224e-196, -0.00023604834149293664);
}
- @Test
- public void testAcosh() {
- // Defined by acos() so edge cases are the same
- // TODO
- }
+ // acosh is defined by acos so is not tested
@Test
- public void testAsinh() {
- // asinh(z) = ln(z + sqrt(1 + z^2))
- // Odd function: negative real cases defined by positive real cases
- // TODO
+ public void testAsin() {
+ // asin(z) = -i (ln(iz + sqrt(1 - z^2)))
+ final String name = "asin";
+ final UnaryOperator<Complex> operation = Complex::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.
+ double huge = Math.sqrt(Double.MAX_VALUE) * 2;
+ double big = Math.sqrt(Double.MAX_VALUE) / 8;
+ double medium = 100;
+ 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);
+ // Additional cases to achieve full coverage
+ // xm1 = 0
+ assertComplex(1, small, name, operation, 1.5707963267948966, 2.4426773395109241e-77);
+ // https://svn.boost.org/trac10/ticket/7290
+ assertComplex(1.00000002785941, 5.72464869028403e-200, name, operation, 1.5707963267948966, 0.00023604834149293664);
}
+ // asinh is defined by asin so is not tested
+
@Test
public void testAtanh() {
// atanh(z) = (1/2) ln((1 + z) / (1 - z))
[commons-numbers] 13/30: Drop exact log10 test. The method of
computation cannot be assumed.
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 12c0b57ccda1269ae020636051adc8ca9614adfc
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 17:32:54 2019 +0000
Drop exact log10 test. The method of computation cannot be assumed.
---
.../src/test/java/org/apache/commons/numbers/complex/ComplexTest.java | 2 --
1 file changed, 2 deletions(-)
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 d5f54cd..d4189b5 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
@@ -2003,8 +2003,6 @@ public class ComplexTest {
// This is prone to floating-point error so use a delta
Assertions.assertEquals(lnz.getReal() / ln10, log10z.getReal(), 1e-12, "real");
// This test should be exact
- final double abs = z.abs();
- Assertions.assertEquals(Math.log10(abs), log10z.getReal(), "real");
Assertions.assertEquals(lnz.getImaginary(), log10z.getImaginary(), "imag");
}
}
[commons-numbers] 04/30: Update the max ULP for acos/acosh from 36
to 2 for the C reference test.
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 93d68901bbe9bad7041011b8a455a1779a2e918d
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 11:34:22 2019 +0000
Update the max ULP for acos/acosh from 36 to 2 for the C reference test.
---
.../test/java/org/apache/commons/numbers/complex/CReferenceTest.java | 4 ++--
1 file changed, 2 insertions(+), 2 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 0f944b4..f26611c 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
@@ -267,12 +267,12 @@ public class CReferenceTest {
@Test
public void testAcos() {
- assertOperation("acos", Complex::acos, 36);
+ assertOperation("acos", Complex::acos, 2);
}
@Test
public void testAcosh() {
- assertOperation("acosh", Complex::acosh, 36);
+ assertOperation("acosh", Complex::acosh, 2);
}
@Test
[commons-numbers] 23/30: Fix javadoc comment on SAFE_MIN/MAX
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 8375852acc8b58f12b8f0acf51046753569e0e52
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 10:17:46 2019 +0000
Fix javadoc comment on SAFE_MIN/MAX
---
.../src/main/java/org/apache/commons/numbers/complex/Complex.java | 4 ++--
1 file changed, 2 insertions(+), 2 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 c390ccf..c9be094 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
@@ -87,12 +87,12 @@ public final class Complex implements Serializable {
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 the Hull, et al with M the largest normalised floating-point value.
+ * 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 the Hull, et al with u the smallest normalised floating-point value.
+ * 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;
/**
[commons-numbers] 21/30: Add atanh assumption on x close to 1.
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 f5ba494b10ea41cf4dad439ce882457e52d84ceb
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Dec 19 22:56:10 2019 +0000
Add atanh assumption on x close to 1.
---
.../src/test/java/org/apache/commons/numbers/complex/ComplexTest.java | 3 +++
1 file changed, 3 insertions(+)
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 b14166b..9d89490 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
@@ -2041,5 +2041,8 @@ public class ComplexTest {
Assertions.assertEquals(result, Math.log1p(result));
Assertions.assertEquals(result, result - result * result / 2,
"Expected log1p Taylor series to be redundant");
+ // Can we assume if x != 1 then (x-1) is valid for multiplications.
+ Assertions.assertNotEquals(0, 1 - Math.nextUp(1));
+ Assertions.assertNotEquals(0, 1 - Math.nextDown(1));
}
}
[commons-numbers] 29/30: Simplify sinh, cosh and tanh
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 1eae0a46843ba8e0ee2ec785d3f8003af54f3608
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 17:29:01 2019 +0000
Simplify sinh, cosh and tanh
---
.../apache/commons/numbers/complex/Complex.java | 173 +++++++--------------
1 file changed, 52 insertions(+), 121 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 bc91838..a126040 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
@@ -1932,56 +1932,22 @@ public final class Complex implements Serializable {
* @return the hyperbolic cosine of the complex number
*/
private static Complex cosh(double real, double imaginary, ComplexConstructor constructor) {
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- return constructor.create(Math.cosh(real) * Math.cos(imaginary),
- Math.sinh(real) * Math.sin(imaginary));
- }
- // ISO C99: Preserve the even function by mapping to positive
- // f(z) = f(-z)
- double re;
- double im;
- if (negative(real)) {
- re = -real;
- im = -imaginary;
- } else {
- re = real;
- im = imaginary;
- }
- // Special case for real == 0
- return constructor.create(Double.NaN,
- re == 0 ? Math.copySign(0, im) : Double.NaN);
+ // 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 (Double.isInfinite(real)) {
- if (Double.isFinite(imaginary)) {
- if (imaginary == 0) {
- // Determine sign
- final double im = real > 0 ? imaginary : -imaginary;
- return constructor.create(Double.POSITIVE_INFINITY, im);
- }
- // inf * cis(y)
- // ISO C99: Preserve the even function
- // f(z) = f(-z)
- double re;
- double im;
- if (real < 0) {
- re = -real;
- im = -imaginary;
- } else {
- re = real;
- im = imaginary;
- }
- return constructor.create(re * Math.cos(im), re * Math.sin(im));
- }
- // imaginary is infinite or NaN
- return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
+ if (real == 0 && !Double.isFinite(imaginary)) {
+ return constructor.create(Double.NaN, changeSign(real, imaginary));
}
- // real is NaN
- if (imaginary == 0) {
- return constructor.create(Double.NaN, imaginary);
+ if (real == 0 && imaginary == 0) {
+ return constructor.create(1, changeSign(real, imaginary));
}
- // optionally raises the ‘‘invalid’’ floating-point exception, for nonzero y.
- return NAN;
+ if (imaginary == 0 && !Double.isFinite(real)) {
+ return constructor.create(Math.abs(real), changeSign(imaginary, real));
+ }
+ return constructor.create(Math.cosh(real) * Math.cos(imaginary),
+ Math.sinh(real) * Math.sin(imaginary));
}
/**
@@ -2008,7 +1974,7 @@ public final class Complex implements Serializable {
im = Math.copySign(1, im);
}
} else {
- if (im == 0 || !Double.isFinite(im)){
+ if (im == 0 || !Double.isFinite(im)) {
return Double.isInfinite(im) ?
new Complex(real, Double.NaN) :
this;
@@ -2236,41 +2202,15 @@ public final class Complex implements Serializable {
* @return the hyperbolic sine of the complex number
*/
private static Complex sinh(double real, double imaginary, ComplexConstructor constructor) {
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- return constructor.create(Math.sinh(real) * Math.cos(imaginary),
- Math.cosh(real) * Math.sin(imaginary));
- }
- // Special case for real == 0
- final double re = real == 0 ? real : Double.NaN;
- return constructor.create(re, Double.NaN);
+ if ((Double.isInfinite(real) && !Double.isFinite(imaginary)) ||
+ (real == 0 && !Double.isFinite(imaginary))) {
+ return constructor.create(real, Double.NaN);
}
- if (Double.isInfinite(real)) {
- if (Double.isFinite(imaginary)) {
- if (imaginary == 0) {
- return constructor.create(real, imaginary);
- }
- // inf * cis(y)
- // ISO C99: Preserve the equality
- // sinh(conj(z)) = conj(sinh(z))
- // and the odd function: f(z) = -f(-z)
- // by always computing on a positive valued Complex number.
- // Math.cos(-x) == Math.cos(x) so ignore sign transform.
- final double signIm = imaginary < 0 ? -1 : 1;
- final double re = Double.POSITIVE_INFINITY * Math.cos(imaginary);
- final double im = Double.POSITIVE_INFINITY * Math.sin(imaginary * signIm);
- // Transform back
- return constructor.create(real < 0 ? -re : re, im * signIm);
- }
- // imaginary is infinite or NaN
- return constructor.create(Double.POSITIVE_INFINITY, Double.NaN);
+ if (imaginary == 0 && !Double.isFinite(real)) {
+ return constructor.create(real, imaginary);
}
- // real is NaN
- if (imaginary == 0) {
- return constructor.create(Double.NaN, Math.copySign(0, imaginary));
- }
- // optionally raises the ‘‘invalid’’ floating-point exception, for nonzero y.
- return NAN;
+ return constructor.create(Math.sinh(real) * Math.cos(imaginary),
+ Math.cosh(real) * Math.sin(imaginary));
}
/**
@@ -2438,42 +2378,6 @@ public final class Complex implements Serializable {
* @return the hyperbolic tangent of the complex number
*/
private static Complex tanh(double real, double imaginary, ComplexConstructor constructor) {
- // Math.cos and Math.sin return NaN for infinity.
- // Perform edge-condition checks on twice the imaginary value.
- // This handles very big imaginary numbers as infinite.
-
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- if (real == 0) {
- // Identity: sin x / (1 + cos x) = tan(x/2)
- return constructor.create(real, Math.tan(imaginary));
- }
- if (imaginary == 0) {
- // Identity: sinh x / (1 + cosh x) = tanh(x/2)
- return constructor.create(Math.tanh(real), imaginary);
- }
-
- final double real2 = 2 * real;
-
- // Math.cosh returns positive infinity for infinity.
- // cosh -> inf
- final double divisor = Math.cosh(real2) + cos2(imaginary);
-
- // Math.sinh returns the input infinity for infinity.
- // sinh -> inf for positive x; else -inf
- final double sinhRe2 = Math.sinh(real2);
-
- // Avoid inf / inf
- if (Double.isInfinite(sinhRe2) && Double.isInfinite(divisor)) {
- // Handle as if real was infinite
- return constructor.create(Math.copySign(1, real), Math.copySign(0, sin2(imaginary)));
- }
- return constructor.create(sinhRe2 / divisor,
- sin2(imaginary) / divisor);
- }
- // imaginary is infinite or NaN
- return NAN;
- }
if (Double.isInfinite(real)) {
if (Double.isFinite(imaginary)) {
return constructor.create(Math.copySign(1, real), Math.copySign(0, sin2(imaginary)));
@@ -2481,12 +2385,39 @@ public final class Complex implements Serializable {
// imaginary is infinite or NaN
return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
}
- // real is NaN
+
+ if (real == 0) {
+ if (Double.isFinite(imaginary)) {
+ // Identity: sin x / (1 + cos x) = tan(x/2)
+ return constructor.create(real, Math.tan(imaginary));
+ }
+ return constructor.create(Double.NaN, Double.NaN);
+ }
if (imaginary == 0) {
- return constructor.create(Double.NaN, Math.copySign(0, imaginary));
+ if (Double.isNaN(real)) {
+ return constructor.create(Double.NaN, imaginary);
+ }
+ // Identity: sinh x / (1 + cosh x) = tanh(x/2)
+ return constructor.create(Math.tanh(real), imaginary);
+ }
+
+ final double real2 = 2 * real;
+
+ // Math.cosh returns positive infinity for infinity.
+ // cosh -> inf
+ final double divisor = Math.cosh(real2) + cos2(imaginary);
+
+ // Math.sinh returns the input infinity for infinity.
+ // sinh -> inf for positive x; else -inf
+ final double sinhRe2 = Math.sinh(real2);
+
+ // Avoid inf / inf
+ if (Double.isInfinite(sinhRe2) && Double.isInfinite(divisor)) {
+ // Handle as if real was infinite
+ return constructor.create(Math.copySign(1, real), Math.copySign(0, imaginary));
}
- // optionally raises the ‘‘invalid’’ floating-point exception, for nonzero y.
- return NAN;
+ return constructor.create(sinhRe2 / divisor,
+ sin2(imaginary) / divisor);
}
/**
[commons-numbers] 01/30: Updated sinh/asinh using an
overflow/underflow safe method.
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 4dd3fc82b1017834c2947682e20f2a13e92082c4
Author: aherbert <ah...@apache.org>
AuthorDate: Wed Dec 18 21:34:21 2019 +0000
Updated sinh/asinh using an overflow/underflow safe method.
---
.../apache/commons/numbers/complex/Complex.java | 218 +++++++++++++++++++--
.../commons/numbers/complex/CReferenceTest.java | 2 +-
2 files changed, 208 insertions(+), 12 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 3d5f264..feb777c 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
@@ -78,6 +78,36 @@ public final class Complex implements Serializable {
/** Base 10 logarithm of 2 (log10(2)). */
private static final double LOG10_2 = Math.log10(2);
+ /** Crossover point to switch computation for asin factor A. */
+ private static final double A_CROSSOVER = 10;
+ /** Crossover point to switch computation for asin factor B. */
+ private static final double B_CROSSOVER = 0.6471;
+ /**
+ * The safe maximum double value {@code x} to avoid loss of precision in asin.
+ * Equal to sqrt(M) / 8 in the Hull, et al 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.
+ * Equal to sqrt(u) * 4 in the Hull, et al with u the smallest normalised floating-point value.
+ */
+ private static final double SAFE_MIN = Math.sqrt(Double.MIN_NORMAL) * 4;
+ /** Exponent offset in IEEE754 representation. */
+ private static final long EXPONENT_OFFSET = 1023L;
+ /**
+ * 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>.</p>
+ *
+ * <p>Copied from o.a.c.numbers.core.Precision
+ *
+ * @see <a href="http://en.wikipedia.org/wiki/Machine_epsilon">Machine epsilon</a>
+ */
+ public static final double EPSILON = Double.longBitsToDouble((EXPONENT_OFFSET - 53L) << 52);
+
/** Serializable version identifier. */
private static final long serialVersionUID = 20180201L;
@@ -1312,18 +1342,162 @@ public final class Complex implements Serializable {
* </code>
* </pre>
*
- * <p>As per the C.99 standard this function is computed using the trigonomic identity:</p>
+ * <p>This is implemented using real {@code x} and imaginary {@code y} parts:</p>
* <pre>
- * asin(z) = -i asinh(iz)
+ * <code>
+ * asin(z) = asin(B) + i sign(y)log(A + sqrt(A<sup>2</sup>-1))
+ * A = 0.5 [ sqrt((x+1)<sup>2</sup>+y<sup>2</sup>) + sqrt((x-1)<sup>2</sup>+y<sup>2</sup>) ]
+ * B = 0.5 [ sqrt((x+1)<sup>2</sup>+y<sup>2</sup>) - sqrt((x-1)<sup>2</sup>+y<sup>2</sup>) ]
+ * sign(y) = {@link Math#copySign(double,double) copySign(1.0, y)}
+ * </code>
* </pre>
*
+ * <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>}. The function is well
+ * defined over the entire complex number range, and produces accurate values even at the
+ * extremes due to special handling of overflow and underflow conditions.</p>
+ *
* @return the inverse sine of this complex number
*/
public Complex asin() {
- // Define in terms of asinh
- // asin(z) = -i asinh(iz)
- // Multiply this number by I, compute asinh, then multiply by back
- return asinh(-imaginary, real, Complex::multiplyNegativeI);
+ return asin(real, imaginary, Complex::ofCartesian);
+ }
+
+ /**
+ * Compute the inverse sine of the complex number.
+ *
+ * <p>This function exists to allow implementation of the identity
+ * {@code asinh(z) = -i asin(iz)}.<p>
+ *
+ * <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>}.</p>
+ *
+ * @param real Real part.
+ * @param imaginary Imaginary part.
+ * @param constructor Constructor.
+ * @return the inverse sine of this complex number
+ */
+ private static Complex asin(double real, double imaginary, 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)) {
+ 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);
+ }
+
+ double xp1 = x + 1;
+ double xm1 = x - 1;
+
+ if ((x < SAFE_MAX) && (x > SAFE_MIN) && (y < SAFE_MAX) && (y > SAFE_MIN)) {
+ double yy = y * y;
+ double r = Math.sqrt(xp1 * xp1 + yy);
+ double s = Math.sqrt(xm1 * xm1 + yy);
+ double a = 0.5 * (r + s);
+ double b = x / a;
+
+ if (b <= B_CROSSOVER) {
+ re = Math.asin(b);
+ } else {
+ 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 3
+ if (y <= (EPSILON * Math.abs(xm1))) {
+ if (x < 1) {
+ re = Math.asin(x);
+ im = y / Math.sqrt(-xp1 * xm1);
+ } 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);
+ double xoy = x / y;
+ im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
+ } else {
+ 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));
}
/**
@@ -1361,17 +1535,27 @@ public final class Complex implements Serializable {
*
* <p>This is an odd function: {@code f(z) = -f(-z)}.
*
+ * <p>This function is computed using the trigonomic identity:</p>
+ * <pre>
+ * asinh(z) = -i asin(iz)
+ * </pre>
+ *
* @return the inverse hyperbolic sine of this complex number
*/
public Complex asinh() {
- return asinh(real, imaginary, Complex::ofCartesian);
+ // Define in terms of asin
+ // asinh(z) = -i asin(iz)
+ // Note: This is the opposite the the identity defined in the C.99 standard:
+ // asin(z) = -i asinh(iz)
+ // Multiply this number by I, compute asin, then multiply by back
+ return asin(-imaginary, real, Complex::multiplyNegativeI);
}
/**
* Compute the inverse hyperbolic sine of the complex number.
*
* <p>This function exists to allow implementation of the identity
- * {@code sin(z) = -i sinh(iz)}.<p>
+ * {@code asin(z) = -i asinh(iz)}.<p>
*
* @param real Real part.
* @param imaginary Imaginary part.
@@ -1460,7 +1644,7 @@ public final class Complex implements Serializable {
* Compute the inverse hyperbolic tangent of this complex number.
*
* <p>This function exists to allow implementation of the identity
- * {@code sin(z) = -i sinh(iz)}.<p>
+ * {@code atan(z) = -i atanh(iz)}.<p>
*
* @param real Real part.
* @param imaginary Imaginary part.
@@ -1962,7 +2146,7 @@ public final class Complex implements Serializable {
* <ul>
* <li>{@code |a| = }{@link Math#abs}(a)
* <li>{@code |a + b i| = }{@link Complex#abs}(a + b i)
- * <li>{@code sign(b) = }{@link Math#copySign(double,double) copySign(1.0, b)}
+ * <li>{@code sign(b) = }{@link Math#copySign(double,double) copySign(1.0, b)}
* </ul>
*
* @return the square root of {@code this}.
@@ -2320,12 +2504,24 @@ public final class Complex implements Serializable {
}
/**
+ * 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 have been
+ * set using {@link Math#abs(double)}).
+ *
+ * @param d Value.
+ * @return {@code true} if {@code d} is +inf.
+ */
+ private static boolean isPosInfinite(double d) {
+ return d == Double.POSITIVE_INFINITY;
+ }
+
+ /**
* 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).multiply(new Complex(0, -1));
+ * z = new Complex(real, imaginary).multiplyImaginary(-1);
* </pre>
*
* @param real Real part.
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 57ad59f..0f944b4 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
@@ -278,7 +278,7 @@ public class CReferenceTest {
@Test
public void testAsinh() {
// Odd function: negative real cases defined by positive real cases
- assertOperation("asinh", Complex::asinh, 1);
+ assertOperation("asinh", Complex::asinh, 3);
}
@Test
[commons-numbers] 10/30: Update the max ULP for log from 3 to 1 for
the C reference test.
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 4bbbef917075a6670ae5289f7ea646a86d21f6f3
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 16:52:49 2019 +0000
Update the max ULP for log from 3 to 1 for the C reference test.
---
.../test/java/org/apache/commons/numbers/complex/CReferenceTest.java | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
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 38090c1..3b02907 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
@@ -312,7 +312,7 @@ public class CReferenceTest {
@Test
public void testLog() {
- assertOperation("log", Complex::log, 3);
+ assertOperation("log", Complex::log, 1);
}
@Test
[commons-numbers] 25/30: Added links to Wolfram elementary
functions where applicable.
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 1f7af2c2a634b70121de3e4d079cfd508cffe5af
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 13:10:06 2019 +0000
Added links to Wolfram elementary functions where applicable.
---
.../java/org/apache/commons/numbers/complex/Complex.java | 15 +++++++++++++++
1 file changed, 15 insertions(+)
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 4516be4..c00ac5e 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
@@ -1296,6 +1296,7 @@ public final class Complex implements Serializable {
* extremes due to special handling of overflow and underflow conditions.</p>
*
* @return the inverse cosine of this complex number.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCos/">ArcCos</a>
*/
public Complex acos() {
return acos(real, imaginary, Complex::ofCartesian);
@@ -1461,6 +1462,7 @@ public final class Complex implements Serializable {
* extremes due to special handling of overflow and underflow conditions.</p>
*
* @return the inverse sine of this complex number
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSin/">ArcSin</a>
*/
public Complex asin() {
return asin(real, imaginary, Complex::ofCartesian);
@@ -1612,6 +1614,7 @@ public final class Complex implements Serializable {
* </pre>
*
* @return the inverse tangent of this complex number
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTan/">ArcTan</a>
*/
public Complex atan() {
// Define in terms of atanh
@@ -1639,6 +1642,7 @@ public final class Complex implements Serializable {
* </pre>
*
* @return the inverse hyperbolic sine of this complex number
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcSinh/">ArcSinh</a>
*/
public Complex asinh() {
// Define in terms of asin
@@ -1853,6 +1857,7 @@ public final class Complex implements Serializable {
* and compatibility with the C.99 standard.</p>
*
* @return the inverse hyperbolic cosine of this complex number
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcCosh/">ArcCosh</a>
*/
public Complex acosh() {
// Define in terms of acos
@@ -1888,6 +1893,7 @@ public final class Complex implements Serializable {
* </pre>
*
* @return the cosine of this complex number.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cos/">Cos</a>
*/
public Complex cos() {
// Define in terms of cosh
@@ -1908,6 +1914,7 @@ public final class Complex implements Serializable {
* <p>This is an even function: {@code f(z) = f(-z)}.
*
* @return the hyperbolic cosine of this complex number.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Cosh/">Cosh</a>
*/
public Complex cosh() {
return cosh(real, imaginary, Complex::ofCartesian);
@@ -1987,6 +1994,7 @@ public final class Complex implements Serializable {
* </pre>
*
* @return <code><i>e</i><sup>this</sup></code>.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</a>
*/
public Complex exp() {
if (Double.isFinite(real)) {
@@ -2046,6 +2054,7 @@ public final class Complex implements Serializable {
* @see Math#log(double)
* @see #abs()
* @see #arg()
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/">Log</a>
*/
public Complex log() {
return log(Math::log, LN_2, Complex::ofCartesian);
@@ -2136,6 +2145,7 @@ public final class Complex implements Serializable {
*
* @param x exponent to which this {@code Complex} is to be raised.
* @return <code>this<sup>x</sup></code>.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
*/
public Complex pow(Complex x) {
if (real == 0 &&
@@ -2167,6 +2177,7 @@ public final class Complex implements Serializable {
* @param x exponent to which this {@code Complex} is to be raised.
* @return <code>this<sup>x</sup></code>.
* @see #pow(Complex)
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Power/">Power</a>
*/
public Complex pow(double x) {
if (real == 0 &&
@@ -2197,6 +2208,7 @@ public final class Complex implements Serializable {
* </pre>
*
* @return the sine of this complex number.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sin/">Sin</a>
*/
public Complex sin() {
// Define in terms of sinh
@@ -2217,6 +2229,7 @@ public final class Complex implements Serializable {
* <p>This is an odd function: {@code f(z) = -f(-z)}.
*
* @return the hyperbolic sine of {@code this}.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sinh/">Sinh</a>
*/
public Complex sinh() {
return sinh(real, imaginary, Complex::ofCartesian);
@@ -2289,6 +2302,7 @@ public final class Complex implements Serializable {
* </ul>
*
* @return the square root of {@code this}.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Sqrt/">Sqrt</a>
*/
public Complex sqrt() {
return sqrt(real, imaginary);
@@ -2581,6 +2595,7 @@ public final class Complex implements Serializable {
* @param n Degree of root.
* @return a List of all {@code n}-th roots of {@code this}.
* @throws IllegalArgumentException if {@code n} is zero.
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Root/">Root</a>
*/
public List<Complex> nthRoot(int n) {
if (n == 0) {
[commons-numbers] 18/30: Moved the check for the safe region to a
method.
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 bf2b32cbc2a5dc3785ebdb52b23c90f1008f15d7
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Dec 19 19:27:20 2019 +0000
Moved the check for the safe region to a method.
---
.../org/apache/commons/numbers/complex/Complex.java | 21 +++++++++++++++++----
1 file changed, 17 insertions(+), 4 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 b0583ff..c390ccf 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
@@ -1353,7 +1353,7 @@ public final class Complex implements Serializable {
final double xp1 = x + 1;
final double xm1 = x - 1;
- if ((x < SAFE_MAX) && (x > SAFE_MIN) && (y < SAFE_MAX) && (y > SAFE_MIN)) {
+ 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);
@@ -1524,7 +1524,7 @@ public final class Complex implements Serializable {
final double xp1 = x + 1;
final double xm1 = x - 1;
- if ((x < SAFE_MAX) && (x > SAFE_MIN) && (y < SAFE_MAX) && (y > SAFE_MIN)) {
+ 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);
@@ -1723,10 +1723,10 @@ public final class Complex implements Serializable {
// 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_UPPER: (1-x) == -x
// x <= SAFE_LOWER: 1 - x^2 = 1
- if ((x > SAFE_LOWER) && (x < SAFE_UPPER) && (y > SAFE_LOWER) && (y < SAFE_UPPER)) {
+ if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
// Normal computation within a safe region.
// minus x plus 1: (-x+1)
@@ -2689,6 +2689,19 @@ public final class Complex implements Serializable {
}
/**
+ * 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);
+ }
+
+ /**
* Creates an exception.
*
* @param message Message prefix.
[commons-numbers] 28/30: Simplify exp() edge cases.
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 e2e7b8b9fba25f181d1254f41db189a4543b480b
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 16:47:45 2019 +0000
Simplify exp() edge cases.
---
.../apache/commons/numbers/complex/Complex.java | 65 +++++++++-------------
1 file changed, 27 insertions(+), 38 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 7c9d962..bc91838 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
@@ -1997,48 +1997,37 @@ public final class Complex implements Serializable {
* @see <a href="http://functions.wolfram.com/ElementaryFunctions/Exp/">Exp</a>
*/
public Complex exp() {
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- final double expReal = Math.exp(real);
- if (imaginary == 0) {
- // Preserve sign for conjugate equality
- return new Complex(expReal, imaginary);
- }
- return new Complex(expReal * Math.cos(imaginary),
- expReal * Math.sin(imaginary));
- }
- // Imaginary is infinite or nan
- return NAN;
- }
+ // Set the values used to compute exp(real) * cis(im)
+ double expReal;
+ double im = imaginary;
if (Double.isInfinite(real)) {
- if (Double.isFinite(imaginary)) {
- if (real == Double.POSITIVE_INFINITY) {
- if (imaginary == 0) {
- return this;
- }
- // inf * cis(y)
- final double re = Double.POSITIVE_INFINITY * Math.cos(imaginary);
- final double im = Double.POSITIVE_INFINITY * Math.sin(imaginary);
- return new Complex(re, im);
+ if (real < 0) {
+ expReal = 0;
+ if (!Double.isFinite(im)) {
+ // Preserve conjugate equality
+ im = Math.copySign(1, im);
}
- // +0 * cis(y)
- final double re = 0.0 * Math.cos(imaginary);
- final double im = 0.0 * Math.sin(imaginary);
- return new Complex(re, im);
- }
- // imaginary is infinite or NaN
- if (real == Double.POSITIVE_INFINITY) {
- return new Complex(Double.POSITIVE_INFINITY, Double.NaN);
+ } else {
+ if (im == 0 || !Double.isFinite(im)){
+ return Double.isInfinite(im) ?
+ new Complex(real, Double.NaN) :
+ this;
+ }
+ expReal = real;
}
- // Preserve sign for conjugate equality
- return new Complex(0, Math.copySign(0, imaginary));
- }
- // real is NaN
- if (imaginary == 0) {
- return new Complex(Double.NaN, Math.copySign(0, imaginary));
+ } else if (imaginary == 0) {
+ // Real-only number
+ return Double.isNaN(real) ?
+ this :
+ new Complex(Math.exp(real), imaginary);
+ } else if (Double.isNaN(real)) {
+ return NAN;
+ } else {
+ // real is finite
+ expReal = Math.exp(real);
}
- // optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
- return NAN;
+ return new Complex(expReal * Math.cos(im),
+ expReal * Math.sin(im));
}
/**
[commons-numbers] 15/30: Switch atanh if statements to fix PMD
warning.
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 9aa23ab5738c67892e204f1c19598936460c1aa7
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 17:43:12 2019 +0000
Switch atanh if statements to fix PMD warning.
if (x != y) { 1; } else { 2; }
to:
if (x == y) { 2; } else { 1; }
---
.../java/org/apache/commons/numbers/complex/Complex.java | 12 ++++++------
1 file changed, 6 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 684faa7..5e6fab3 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
@@ -1786,12 +1786,7 @@ public final class Complex implements Serializable {
// This will tend towards 0 and log1p(0) = 0 so it may not matter.
re = Math.log1p(4 * x / y / y);
}
- } else if (x != 1) {
- // 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));
- } else {
+ } else if (x == 1) {
// x = 1, small y:
// Special case when x == 1 as (1-x) is invalid.
// Simplify the following formula:
@@ -1803,6 +1798,11 @@ public final class Complex implements Serializable {
// 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:
[commons-numbers] 16/30: Remove use of log1p when value is very
small.
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 244f6e7a0d93ed92d0748b062d29a59333217b5e
Author: Alex Herbert <ah...@apache.org>
AuthorDate: Thu Dec 19 19:00:43 2019 +0000
Remove use of log1p when value is very small.
---
.../src/main/java/org/apache/commons/numbers/complex/Complex.java | 6 +++---
.../test/java/org/apache/commons/numbers/complex/ComplexTest.java | 5 +++++
2 files changed, 8 insertions(+), 3 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 5e6fab3..33b6056 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
@@ -1782,9 +1782,9 @@ public final class Complex implements Serializable {
re = Math.log1p((4 * x / y) / (mxp1 * mxp1 / y + y));
} else {
// Big y, small x, as above but neglect (1-x)^2/y:
- // Note: The boost version has no log1p here.
- // This will tend towards 0 and log1p(0) = 0 so it may not matter.
- re = Math.log1p(4 * x / y / 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:
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 28601fb..b14166b 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
@@ -2036,5 +2036,10 @@ public class ComplexTest {
Assertions.assertEquals(1, (1 - safeLower) * (1 - safeLower));
// Can we assume 1 - y^2 = 1 when y is small
Assertions.assertEquals(1, 1 - safeLower * safeLower);
+ // Can we assume Math.log1p(4 * x / y / y) = (4 * x / y / y) when big y and small x
+ final double result = 4 * safeLower / safeUpper / safeUpper;
+ Assertions.assertEquals(result, Math.log1p(result));
+ Assertions.assertEquals(result, result - result * result / 2,
+ "Expected log1p Taylor series to be redundant");
}
}
[commons-numbers] 08/30: Use final
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 bdaecba1eb69fc44349362640230c93ce1a07366
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 16:23:34 2019 +0000
Use final
---
.../apache/commons/numbers/complex/Complex.java | 62 +++++++++++-----------
1 file changed, 31 insertions(+), 31 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 f1244c7..7577ddf 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
@@ -1315,8 +1315,8 @@ public final class Complex implements Serializable {
private static Complex acos(final double real, final double imaginary,
final ComplexConstructor constructor) {
// Compute with positive values and determine sign at the end
- double x = Math.abs(real);
- double y = Math.abs(imaginary);
+ final double x = Math.abs(real);
+ final double y = Math.abs(imaginary);
// The result (without sign correction)
double re;
double im;
@@ -1350,20 +1350,20 @@ public final class Complex implements Serializable {
return constructor.create(x == 0 ? PI_OVER_2 : Math.acos(real), -imaginary);
}
- double xp1 = x + 1;
- double xm1 = x - 1;
+ final double xp1 = x + 1;
+ final double xm1 = x - 1;
if ((x < SAFE_MAX) && (x > SAFE_MIN) && (y < SAFE_MAX) && (y > SAFE_MIN)) {
- double yy = y * y;
- double r = Math.sqrt(xp1 * xp1 + yy);
- double s = Math.sqrt(xm1 * xm1 + yy);
- double a = 0.5 * (r + s);
- double b = x / a;
+ 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 {
- double apx = a + x;
+ final double apx = a + x;
if (x <= 1) {
re = Math.atan(Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))) / x);
} else {
@@ -1414,11 +1414,11 @@ public final class Complex implements Serializable {
im = LN_2 + Math.log(y);
} else if (x > 1) {
re = Math.atan(y / x);
- double xoy = x / y;
+ final double xoy = x / y;
im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
} else {
re = PI_OVER_2;
- double a = Math.sqrt(1 + y * y);
+ final double a = Math.sqrt(1 + y * y);
im = 0.5 * Math.log1p(2 * y * (y + a));
}
}
@@ -1483,8 +1483,8 @@ public final class Complex implements Serializable {
private static Complex asin(final double real, final double imaginary,
final ComplexConstructor constructor) {
// Compute with positive values and determine sign at the end
- double x = Math.abs(real);
- double y = Math.abs(imaginary);
+ final double x = Math.abs(real);
+ final double y = Math.abs(imaginary);
// The result (without sign correction)
double re;
double im;
@@ -1521,20 +1521,20 @@ public final class Complex implements Serializable {
return constructor.create(Math.asin(real), imaginary);
}
- double xp1 = x + 1;
- double xm1 = x - 1;
+ final double xp1 = x + 1;
+ final double xm1 = x - 1;
if ((x < SAFE_MAX) && (x > SAFE_MIN) && (y < SAFE_MAX) && (y > SAFE_MIN)) {
- double yy = y * y;
- double r = Math.sqrt(xp1 * xp1 + yy);
- double s = Math.sqrt(xm1 * xm1 + yy);
- double a = 0.5 * (r + s);
- double b = x / a;
+ 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 {
- double apx = a + x;
+ final double apx = a + x;
if (x <= 1) {
re = Math.atan(x / Math.sqrt(0.5 * apx * (yy / (r + xp1) + (s - xm1))));
} else {
@@ -1583,10 +1583,10 @@ public final class Complex implements Serializable {
im = LN_2 + Math.log(y);
} else if (x > 1) {
re = Math.atan(x / y);
- double xoy = x / y;
+ final double xoy = x / y;
im = LN_2 + Math.log(y) + 0.5 * Math.log1p(xoy * xoy);
} else {
- double a = Math.sqrt(1 + y * y);
+ final double a = Math.sqrt(1 + y * y);
// Possible underflow:
re = x / a;
im = 0.5 * Math.log1p(2 * y * (y + a));
@@ -1693,8 +1693,8 @@ public final class Complex implements Serializable {
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);
+ final double x = Math.abs(real);
+ final double y = Math.abs(imaginary);
// The result (without sign correction)
double re;
double im;
@@ -1731,8 +1731,8 @@ public final class Complex implements Serializable {
// Normal computation within a safe region.
// minus x plus 1: (-x+1)
- double mxp1 = 1 - x;
- double yy = y * y;
+ 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:
@@ -1785,7 +1785,7 @@ public final class Complex implements Serializable {
} else if (y >= SAFE_UPPER) {
if (x > 1) {
// Big y, medium x, divide through by y:
- double mxp1 = 1 - x;
+ final double mxp1 = 1 - x;
re = Math.log1p((4 * x / y) / (y + mxp1 * mxp1 / y));
} else {
// Big y, small x, as above but neglect (1-x)^2/y:
@@ -1796,7 +1796,7 @@ public final class Complex implements Serializable {
} else if (x != 1) {
// Modified from boost which checks y > SAFE_LOWER.
// if y*y -> 0 it will be ignored so always include it.
- double mxp1 = 1 - x;
+ final double mxp1 = 1 - x;
re = Math.log1p((4 * x) / (mxp1 * mxp1 + y * y));
} else {
// x = 1, small y:
@@ -2324,7 +2324,7 @@ public final class Complex implements Serializable {
return new Complex(sqrtAbs, imaginary);
}
// Get the absolute of the real
- double absA = Math.abs(real);
+ final double absA = Math.abs(real);
// Compute |a + b i|
double absC = getAbsolute(real, imaginary);
[commons-numbers] 11/30: Removed unused constants
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 1d20e9fa8ad6cd6b416fba8e4091cd2f992d4445
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 16:57:58 2019 +0000
Removed unused constants
---
.../test/java/org/apache/commons/numbers/complex/ComplexTest.java | 5 -----
1 file changed, 5 deletions(-)
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 3a0adce..ec395e2 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
@@ -42,19 +42,14 @@ public class ComplexTest {
private static final Complex infOne = Complex.ofCartesian(inf, 1);
private static final Complex infZero = Complex.ofCartesian(inf, 0);
private static final Complex infNegZero = Complex.ofCartesian(inf, -0.0);
- private static final Complex infNan = Complex.ofCartesian(inf, nan);
private static final Complex infNegInf = Complex.ofCartesian(inf, neginf);
private static final Complex infInf = Complex.ofCartesian(inf, inf);
private static final Complex negInfInf = Complex.ofCartesian(neginf, inf);
- private static final Complex negInfZero = Complex.ofCartesian(neginf, 0);
private static final Complex negInfOne = Complex.ofCartesian(neginf, 1);
- private static final Complex negInfNan = Complex.ofCartesian(neginf, nan);
private static final Complex negInfNegInf = Complex.ofCartesian(neginf, neginf);
private static final Complex oneNan = Complex.ofCartesian(1, nan);
private static final Complex zeroInf = Complex.ofCartesian(0, inf);
private static final Complex zeroNan = Complex.ofCartesian(0, nan);
- private static final Complex nanInf = Complex.ofCartesian(nan, inf);
- private static final Complex nanNegInf = Complex.ofCartesian(nan, neginf);
private static final Complex nanZero = Complex.ofCartesian(nan, 0);
private static final Complex NAN = Complex.ofCartesian(nan, nan);
private static final Complex INF = Complex.ofCartesian(inf, inf);
[commons-numbers] 27/30: Simplify sqrt() edge cases.
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 3614d51ec5a53324f268ba5b96c975a6e5ac04ce
Author: aherbert <ah...@apache.org>
AuthorDate: Fri Dec 20 15:58:23 2019 +0000
Simplify sqrt() edge cases.
---
.../apache/commons/numbers/complex/Complex.java | 112 ++++++++++-----------
1 file changed, 54 insertions(+), 58 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 c00ac5e..7c9d962 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
@@ -2107,7 +2107,7 @@ public final class Complex implements Serializable {
final double y = Math.abs(imaginary);
// Use the safe region defined for atanh to avoid over/underflow for x^2
- if ((x > SAFE_LOWER) && (x < SAFE_UPPER) && (y > SAFE_LOWER) && (y < SAFE_UPPER)) {
+ if (inRegion(x, y, SAFE_LOWER, SAFE_UPPER)) {
return constructor.create(0.5 * log.apply(x * x + y * y), arg());
}
@@ -2332,60 +2332,6 @@ public final class Complex implements Serializable {
if (Double.isInfinite(imaginary)) {
return new Complex(Double.POSITIVE_INFINITY, imaginary);
}
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- // Edge case for real numbers
- if (imaginary == 0) {
- final double sqrtAbs = Math.sqrt(Math.abs(real));
- if (real < 0) {
- return new Complex(0, Math.copySign(sqrtAbs, imaginary));
- }
- return new Complex(sqrtAbs, imaginary);
- }
- // Get the absolute of the real
- final double absA = Math.abs(real);
- // Compute |a + b i|
- double absC = getAbsolute(real, imaginary);
-
- // t = sqrt((|a| + |a + b i|) / 2)
- // This is always representable as this complex is finite.
- double t;
-
- // Overflow safe
- if (absC == Double.POSITIVE_INFINITY) {
- // Complex is too large.
- // Divide by the largest absolute component,
- // compute the required sqrt and then scale back.
- // Use the equality: sqrt(n) = sqrt(scale) * sqrt(n/scale)
- // t = sqrt(max) * sqrt((|a|/max + |a + b i|/max) / 2)
- // Note: The function may be non-monotonic at the junction.
- // The alternative of returning infinity for a finite input is worse.
- // Use precise scaling with:
- // scale ~ 2^exponent
- // Make exponent even for fast rescaling using sqrt(2^exponent).
- final int exponent = getMaxExponent(absA, imaginary) & MASK_INT_TO_EVEN;
- // Implement scaling using 2^-exponent
- final double scaleA = Math.scalb(absA, -exponent);
- final double scaleB = Math.scalb(imaginary, -exponent);
- absC = getAbsolute(scaleA, scaleB);
- // t = Math.sqrt(2^exponent) * Math.sqrt((scaleA + absC) / 2)
- // This works if exponent is even:
- // sqrt(2^exponent) = (2^exponent)^0.5 = 2^(exponent*0.5)
- t = Math.scalb(Math.sqrt((scaleA + absC) / 2), exponent / 2);
- } else {
- // Over-flow safe average: absA < absC and abdC is finite.
- t = Math.sqrt(absA + (absC - absA) / 2);
- }
-
- if (real >= 0) {
- return new Complex(t, imaginary / (2 * t));
- }
- return new Complex(Math.abs(imaginary) / (2 * t),
- Math.copySign(1.0, imaginary) * t);
- }
- // Imaginary is nan
- return NAN;
- }
if (Double.isInfinite(real)) {
// imaginary is finite or NaN
final double part = Double.isNaN(imaginary) ? Double.NaN : 0;
@@ -2394,9 +2340,59 @@ public final class Complex implements Serializable {
}
return new Complex(Double.POSITIVE_INFINITY, Math.copySign(part, imaginary));
}
- // real is NaN
- // optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
- return NAN;
+ if (Double.isNaN(real) || Double.isNaN(imaginary)) {
+ return NAN;
+ }
+ // Finite real and imaginary
+
+ // Edge case for real numbers
+ if (imaginary == 0) {
+ final double sqrtAbs = Math.sqrt(Math.abs(real));
+ if (real < 0) {
+ return new Complex(0, Math.copySign(sqrtAbs, imaginary));
+ }
+ return new Complex(sqrtAbs, imaginary);
+ }
+ // Get the absolute of the real
+ final double absA = Math.abs(real);
+ // Compute |a + b i|
+ double absC = getAbsolute(real, imaginary);
+
+ // t = sqrt((|a| + |a + b i|) / 2)
+ // This is always representable as this complex is finite.
+ double t;
+
+ // Overflow safe
+ if (absC == Double.POSITIVE_INFINITY) {
+ // Complex is too large.
+ // Divide by the largest absolute component,
+ // compute the required sqrt and then scale back.
+ // Use the equality: sqrt(n) = sqrt(scale) * sqrt(n/scale)
+ // t = sqrt(max) * sqrt((|a|/max + |a + b i|/max) / 2)
+ // Note: The function may be non-monotonic at the junction.
+ // The alternative of returning infinity for a finite input is worse.
+ // Use precise scaling with:
+ // scale ~ 2^exponent
+ // Make exponent even for fast rescaling using sqrt(2^exponent).
+ final int exponent = getMaxExponent(absA, imaginary) & MASK_INT_TO_EVEN;
+ // Implement scaling using 2^-exponent
+ final double scaleA = Math.scalb(absA, -exponent);
+ final double scaleB = Math.scalb(imaginary, -exponent);
+ absC = getAbsolute(scaleA, scaleB);
+ // t = Math.sqrt(2^exponent) * Math.sqrt((scaleA + absC) / 2)
+ // This works if exponent is even:
+ // sqrt(2^exponent) = (2^exponent)^0.5 = 2^(exponent*0.5)
+ t = Math.scalb(Math.sqrt((scaleA + absC) / 2), exponent / 2);
+ } else {
+ // Over-flow safe average: absA < absC and abdC is finite.
+ t = Math.sqrt(absA + (absC - absA) / 2);
+ }
+
+ if (real >= 0) {
+ return new Complex(t, imaginary / (2 * t));
+ }
+ return new Complex(Math.abs(imaginary) / (2 * t),
+ Math.copySign(1.0, imaginary) * t);
}
/**
[commons-numbers] 09/30: Compute log() within a safe region without
using the sqrt.
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 67dc19de65579b1a3335919e6459b8946d1fae75
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 16:49:25 2019 +0000
Compute log() within a safe region without using the sqrt.
---
.../java/org/apache/commons/numbers/complex/Complex.java | 12 ++++++++++++
1 file changed, 12 insertions(+)
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 7577ddf..11f468f 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
@@ -2097,6 +2097,18 @@ public final class Complex implements Serializable {
private Complex log(UnaryOperation log, double logOf2, ComplexConstructor constructor) {
// All ISO C99 edge cases satisfied by the Math library.
// Make computation overflow safe.
+
+ // Note:
+ // log(|a + b i|) = log(sqrt(a^2 + b^2)) = 0.5 * log(a^2 + b^2)
+ // If real and imaginary are with a safe region then omit the sqrt().
+ final double x = Math.abs(real);
+ final double y = Math.abs(imaginary);
+
+ // Use the safe region defined for atanh to avoid over/underflow for x^2
+ if ((x > SAFE_LOWER) && (x < SAFE_UPPER) && (y > SAFE_LOWER) && (y < SAFE_UPPER)) {
+ return constructor.create(0.5 * log.apply(x * x + y * y), arg());
+ }
+
final double abs = abs();
if (abs == Double.POSITIVE_INFINITY && isFinite()) {
// Edge-case where the |a + b i| overflows.
[commons-numbers] 06/30: Updated atanh using an overflow/underflow
safe method.
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 5346128780ffd61c914cb76eab3301080fa46c3f
Author: aherbert <ah...@apache.org>
AuthorDate: Thu Dec 19 15:09:07 2019 +0000
Updated atanh using an overflow/underflow safe method.
---
.../apache/commons/numbers/complex/Complex.java | 215 ++++++++++++++++-----
1 file changed, 168 insertions(+), 47 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 33aaec2..f1244c7 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
@@ -86,15 +86,25 @@ public final class Complex implements Serializable {
/** 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.
+ * The safe maximum double value {@code x} to avoid loss of precision in asin/acos.
* Equal to sqrt(M) / 8 in the Hull, et al 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.
+ * The safe minimum double value {@code x} to avoid loss of precision/underflow in asin/acos.
* Equal to sqrt(u) * 4 in the Hull, et al 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;
/** Exponent offset in IEEE754 representation. */
private static final long EXPONENT_OFFSET = 1023L;
/**
@@ -1650,7 +1660,20 @@ public final class Complex implements Serializable {
*
* <p>This is an odd function: {@code f(z) = -f(-z)}.
*
+ * <p>This is implemented using real {@code x} and imaginary {@code y} parts:</p>
+ * <pre>
+ * <code>
+ * atanh(z) = 0.25 ln(1 + 4x/((1-x)<sup>2</sup>+y<sup>2</sup>) + i 0.5 tan<sup>-1</sup>(2y, 1-x<sup>2</sup>-y<sup>2</sup>)
+ * </code>
+ * </pre>
+ *
+ * <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>}. The function is well
+ * defined over the entire complex number range, and produces accurate values even at the
+ * extremes due to special handling of overflow and underflow conditions.</p>
+ *
* @return the inverse hyperbolic tangent of this complex number
+ * @see <a href="http://functions.wolfram.com/ElementaryFunctions/ArcTanh/">ArcTanh</a>
*/
public Complex atanh() {
return atanh(real, imaginary, Complex::ofCartesian);
@@ -1667,59 +1690,157 @@ public final class Complex implements Serializable {
* @param constructor Constructor.
* @return the inverse hyperbolic tangent of the complex number
*/
- private static Complex atanh(double real, double imaginary, ComplexConstructor constructor) {
- if (Double.isFinite(real)) {
- if (Double.isFinite(imaginary)) {
- // ISO C99: Preserve the equality
- // atanh(conj(z)) = conj(atanh(z))
- // and the odd function: f(z) = -f(-z)
- // by always computing on a positive valued Complex number.
- final double a = Math.abs(real);
- final double b = Math.abs(imaginary);
- // Special case for divide-by-zero
- if (a == 1 && b == 0) {
- // raises the ‘‘divide-by-zero’’ floating-point exception.
- return constructor.create(Math.copySign(Double.POSITIVE_INFINITY, real), imaginary);
- }
+ private static Complex atanh(final double real, final double imaginary,
+ final ComplexConstructor constructor) {
+ // Compute with positive values and determine sign at the end
+ double x = Math.abs(real);
+ double y = Math.abs(imaginary);
+ // The result (without sign correction)
+ double re;
+ double im;
+
+ // Handle C99 special cases
+ if (Double.isNaN(x)) {
+ if (isPosInfinite(y)) {
+ // The sign of the real part of the result is unspecified
+ return constructor.create(0, Math.copySign(PI_OVER_2, imaginary));
+ }
+ // 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.
+ // Cases for very large finite are handled as if 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^2 -> inf
+ // x <= SAFE_LOWER: x^2 -> 0
+
+ if ((x > SAFE_LOWER) && (x < SAFE_UPPER) && (y > SAFE_LOWER) && (y < SAFE_UPPER)) {
+ // Normal computation within a safe region.
+
+ // minus x plus 1: (-x+1)
+ double mxp1 = 1 - x;
+ 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
+ // The division is done at the end of the function.
+ re = Math.log1p(4 * x / (mxp1 * mxp1 + yy));
+ im = Math.atan2(2 * y, mxp1 * (1 + x) - yy);
+ } 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 (a == 0) {
+ if (x == 0) {
if (imaginary == 0) {
return constructor.create(real, imaginary);
}
// atanh(iy) = i atan(y)
- final double im = Math.atan(imaginary);
- return constructor.create(real, im);
+ 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
+ if (isPosInfinite(x) || isPosInfinite(y)) {
+ re = 0;
+ } else if (y >= SAFE_UPPER) {
+ // Big x and y: divide by x*y
+ // This has been modified from the boost version to
+ // include 1/(x*y) and -2/y. These are harmless if
+ // machine precision prevents their addition to have an effect:
+ // 1/(x*y) -> 0
+ // (x-2) -> x
+ re = Math.log1p((4 / y) / (1 / (x * y) + (x - 2) / y + y / x));
+ } else if (y > 1) {
+ // Big x: divide through by x:
+ // This has been modified from the boost version to
+ // include 1/x and -2:
+ re = Math.log1p(4 / (1 / x + x - 2 + y * y / x));
+ } else {
+ // Big x small y, as above but neglect y^2/x:
+ re = Math.log1p(4 / (1 / x + x - 2));
+ }
+ } else if (y >= SAFE_UPPER) {
+ if (x > 1) {
+ // Big y, medium x, divide through by y:
+ double mxp1 = 1 - x;
+ re = Math.log1p((4 * x / y) / (y + mxp1 * mxp1 / y));
+ } else {
+ // Big y, small x, as above but neglect (1-x)^2/y:
+ // Note: The boost version has no log1p here.
+ // This will tend towards 0 and log1p(0) = 0.
+ re = Math.log1p(4 * x / y / y);
+ }
+ } else if (x != 1) {
+ // Modified from boost which checks y > SAFE_LOWER.
+ // if y*y -> 0 it will be ignored so always include it.
+ double mxp1 = 1 - x;
+ re = Math.log1p((4 * x) / (mxp1 * mxp1 + y * y));
+ } else {
+ // 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));
+ }
+
+ // 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 just do that.
+ // 1 - y^2 = 1 so ignore subtracting y^2
+ im = Math.atan2(2 * y, (1 - x) * (1 + x));
}
- // (1 + (a + b i)) / (1 - (a + b i))
- final Complex result = divide(1 + a, b, 1 - a, -b);
- // Compute the log:
- // (re + im i) = (1/2) * ln((1 + z) / (1 - z))
- return result.log(Math::log, LN_2, (re, im) ->
- // Divide log() by 2 and map back to the correct sign
- constructor.create(0.5 * changeSign(re, real),
- 0.5 * changeSign(im, imaginary))
- );
- }
- if (Double.isInfinite(imaginary)) {
- return constructor.create(Math.copySign(0, real), Math.copySign(PI_OVER_2, imaginary));
- }
- // imaginary is NaN
- // Special case for real == 0
- return real == 0 ? constructor.create(real, Double.NaN) : NAN;
- }
- if (Double.isInfinite(real)) {
- if (Double.isNaN(imaginary)) {
- return constructor.create(Math.copySign(0, real), Double.NaN);
}
- // imaginary is finite or infinite
- return constructor.create(Math.copySign(0, real), Math.copySign(PI_OVER_2, imaginary));
}
- // real is NaN
- if (Double.isInfinite(imaginary)) {
- return constructor.create(0, Math.copySign(PI_OVER_2, imaginary));
- }
- // optionally raises the ‘‘invalid’’ floating-point exception, for finite y.
- return NAN;
+
+ re /= 4;
+ im /= 2;
+ return constructor.create(changeSign(re, real),
+ changeSign(im, imaginary));
}
/**