aherbert commented on code in PR #117: URL: https://github.com/apache/commons-numbers/pull/117#discussion_r922998982
########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexResult.java: ########## @@ -0,0 +1,40 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.numbers.complex; + +/** + * Represents a Terminating consumer for the complex result, used to construct a result of type {@code R}. + * The operation may return an object representation of the complex result. + * + * <p>This is a functional interface whose functional method is + * {@link #apply(double, double)}. + * + * @param <R> The type of the complex result + * @since 1.1 + */ +@FunctionalInterface +public interface ComplexResult<R> { Review Comment: Rename to `ComplexSink`. For consistency the use of `a + ib` should be in brackets: `\( (a + i b) \)` There are a few place where the brackets are not used in Complex (these should be changed in future); most cases seem to have the enclosing brackets. ########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexResult.java: ########## @@ -0,0 +1,40 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.numbers.complex; + +/** + * Represents a Terminating consumer for the complex result, used to construct a result of type {@code R}. Review Comment: terminating in lower case. Since we are replacing the name `result` with `sink` the description should avoid the word result. The interface may be used for example as the operation to add numbers to an extendable list: ```Java // list returns itself list.apply(x, y).apply(r, i). ``` I would keep the doc simple (which can be expanded later): ``` Represents a data sink for a complex number \( (a + i b) \); operations return a result of type {@code R}. ``` ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java: ########## @@ -145,6 +145,39 @@ static void assertEquals(Supplier<String> msg, double expected, double actual, l } } + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on + * complex real and imaginary parts. + * + * <p>The results are considered equal within the provided units of least + * precision. The maximum count of numbers allowed between the two values is + * {@code maxUlps - 1}. + * + * <p>Numbers must have the same sign. Thus -0.0 and 0.0 are never equal. + * + * @param c Input number. + * @param name the operation name + * @param operation1 the operation on the Complex object + * @param operation2 the operation on the complex real and imaginary parts + * @param expected Expected result. + */ + static void assertComplex(Complex c, + String name, + UnaryOperator<Complex> operation1, + ComplexUnaryOperator<ComplexNumber> operation2, + Complex expected, long maxUlps) { + + final Complex z = operation1.apply(c); + final ComplexNumber y = operation2.apply(c.getReal(), c.getImaginary(), ComplexNumber::new); + + assertEquals(() -> "UnaryOperator " + name + "(" + c + "): real", expected.real(), z.real(), maxUlps); Review Comment: These messages on the object are written in OO form (as per the other messages in the class): ``` c + "." + name + "(): e.g. c = (3 + 4i), name = sin (3, 4).sin(): ``` The procedural form should reflect that: ``` name + c: e.g. c = (3 + 4i), name = sin sin(3, 4): ``` ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java: ########## @@ -227,7 +260,27 @@ private static void assertBiOperation(String name, /** * Assert the operation using the data loaded from test resources. * - * @param testData Test data resource name. + * @param name the operation name + * @param operation1 the operation on the Complex object + * @param operation2 the operation on the complex real and imaginary parts + * @param maxUlps the maximum units of least precision between the two values + */ + private static void assertOperation(String name, + UnaryOperator<Complex> operation1, + ComplexUnaryOperator<ComplexNumber> operation2, + long maxUlps) { + Review Comment: Remove extra line ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexNumber.java: ########## @@ -0,0 +1,75 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.numbers.complex; + +/** + * Cartesian representation of a complex number. The complex number is expressed + * in the form \( a + ib \) where \( a \) and \( b \) are real numbers and \( i \) + * is the imaginary unit which satisfies the equation \( i^2 = -1 \). For the + * complex number \( a + ib \), \( a \) is called the <em>real part</em> and + * \( b \) is called the <em>imaginary part</em>. + */ +class ComplexNumber { + + /** The real part. */ + private final double real; + /** The imaginary part. */ + private final double imaginary; + + /** + * Constructor representing a complex number by its real and imaginary parts. + * Takes in real and imaginary and sets it to this complex number's real and imaginary + * + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * + */ + ComplexNumber(double real, double imaginary) { + this.real = real; + this.imaginary = imaginary; + } + + /** + * Creates a conjugated complex number given the real and imaginary parts. + * Review Comment: ``` , that is for the argument (a + ib), returns (a - ib). ``` ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java: ########## @@ -81,6 +81,55 @@ private static void assertComplex(double a, double b, CReferenceTest.assertComplex(c, name, operation, e, maxUlps); } + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on + * complex real and imaginary parts. + * + * <p>The results are considered equal if there are no floating-point values between them. + * + * @param a Real part. + * @param b Imaginary part. + * @param name The operation name. + * @param operation1 the operation on the Complex object. + * @param operation2 the operation on the complex real and imaginary parts + * @param x Expected real part. + * @param y Expected imaginary part. + */ + private static void assertComplex(double a, double b, + String name, UnaryOperator<Complex> operation1, + ComplexUnaryOperator<ComplexNumber> operation2, + double x, double y) { + assertComplex(a, b, name, operation1, operation2, x, y, 1); + } + + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on Review Comment: Move the statement about the second operation to the end. ########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java: ########## @@ -0,0 +1,584 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.numbers.complex; + +import java.util.function.DoubleUnaryOperator; + +/** + * Cartesian representation of a complex number. The complex number is expressed + * in the form \( a + ib \) where \( a \) and \( b \) are real numbers and \( i \) + * is the imaginary unit which satisfies the equation \( i^2 = -1 \). For the + * complex number \( a + ib \), \( a \) is called the <em>real part</em> and + * \( b \) is called the <em>imaginary part</em>. + * + * <p>This class is immutable, it is to change from an OO representation of a Complex number + * to a representation as its real and imaginary parts. All arithmetic will create a new instance for the + * result.</p> + * + * <p>Arithmetic in this class conforms to the C99 standard for complex numbers + * defined in ISO/IEC 9899, Annex G. Methods have been named using the equivalent + * method in ISO C99. The behavior for special cases is listed as defined in C99.</p> + * + * <p>For functions \( f \) which obey the conjugate equality \( conj(f(z)) = f(conj(z)) \), + * the specifications for the upper half-plane imply the specifications for the lower + * half-plane.</p> + * + * <p>For functions that are either odd, \( f(z) = -f(-z) \), or even, \( f(z) = f(-z) \), + * the specifications for the first quadrant imply the specifications for the other three + * quadrants.</p> + * + * <p>Special cases of <a href="http://mathworld.wolfram.com/BranchCut.html";>branch cuts</a> + * for multivalued functions adopt the principle value convention from C99. Specials cases + * from C99 that raise the "invalid" or "divide-by-zero" + * <a href="https://en.cppreference.com/w/c/numeric/fenv/FE_exceptions";>floating-point + * exceptions</a> return the documented value without an explicit mechanism to notify + * of the exception case, that is no exceptions are thrown during computations in-line with + * the convention of the corresponding single-valued functions in + * {@link java.lang.Math java.lang.Math}. + * These cases are documented in the method special cases as "invalid" or "divide-by-zero" + * floating-point operation. + * Note: Invalid floating-point exception cases will result in a complex number where the + * cardinality of NaN component parts has increased as a real or imaginary part could + * not be computed and is set to NaN. + * + * @see <a href="http://www.open-std.org/JTC1/SC22/WG14/www/standards";> + * ISO/IEC 9899 - Programming languages - C</a> + */ +public final class ComplexFunctions { + + /** Mask to remove the sign bit from a long. */ + static final long UNSIGN_MASK = 0x7fff_ffff_ffff_ffffL; + + /** Natural logarithm of 2 (ln(2)). */ + private static final double LN_2 = Math.log(2); + /** {@code 1/2}. */ + private static final double HALF = 0.5; + /** Base 10 logarithm of 10 divided by 2 (log10(e)/2). */ + private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2; + /** Base 10 logarithm of 2 (log10(2)). */ + private static final double LOG10_2 = Math.log10(2); + /** {@code sqrt(2)}. */ + private static final double ROOT2 = 1.4142135623730951; + + /** 54 shifted 20-bits to align with the exponent of the upper 32-bits of a double. */ + private static final int EXP_54 = 0x36_00000; + /** Represents an exponent of 500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */ + private static final int EXP_500 = 0x5f3_00000; + /** Represents an exponent of 1024 in unbiased form (infinite or nan) + * shifted 20-bits to align with the upper 32-bits of a double. */ + private static final int EXP_1024 = 0x7ff_00000; + /** Represents an exponent of -500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */ + private static final int EXP_NEG_500 = 0x20b_00000; + /** 2^600. */ + private static final double TWO_POW_600 = 0x1.0p+600; + /** 2^-600. */ + private static final double TWO_POW_NEG_600 = 0x1.0p-600; + + /** + * Private constructor for utility class. + */ + private ComplexFunctions() { + } + + /** + * Returns the absolute value of the complex number. + * <pre>abs(x + i y) = sqrt(x^2 + y^2)</pre> + * + * <p>This should satisfy the special cases of the hypot function in ISO C99 F.9.4.3: + * "The hypot functions compute the square root of the sum of the squares of x and y, + * without undue overflow or underflow." + * + * <ul> + * <li>hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent. + * <li>hypot(x, ±0) is equivalent to |x|. + * <li>hypot(±∞, y) returns +∞, even if y is a NaN. + * </ul> + * + * <p>This method is called by all methods that require the absolute value of the complex + * number, e.g. abs(), sqrt() and log(). + * + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @return The absolute value. + */ + public static double abs(double real, double imaginary) { + // Specialised implementation of hypot. + // See NUMBERS-143 + return hypot(real, imaginary); + } + + /** + * Returns the argument of the complex number. + * + * <p>The argument is the angle phi between the positive real axis and + * the point representing this number in the complex plane. + * The value returned is between \( -\pi \) (not inclusive) + * and \( \pi \) (inclusive), with negative values returned for numbers with + * negative imaginary parts. + * + * <p>If either real or imaginary part (or both) is NaN, then the result is NaN. + * Infinite parts are handled as {@linkplain Math#atan2} handles them, + * essentially treating finite parts as zero in the presence of an + * infinite coordinate and returning a multiple of \( \frac{\pi}{4} \) depending on + * the signs of the infinite parts. + * + * <p>This code follows the + * <a href="http://www.iso-9899.info/wiki/The_Standard";>ISO C Standard</a>, Annex G, + * in calculating the returned value using the {@code atan2(y, x)} method for complex + * \( x + iy \). + * + * @param r Real part \( a \) of the complex number \(a +ib \). + * @param i Imaginary part \( b \) of the complex number \(a +ib \). + * @return The argument of the complex number. + * @see Math#atan2(double, double) + */ + public static double arg(double r, double i) { + // Delegate + return Math.atan2(i, r); + } + + /** + * Returns {@code true} if either real or imaginary component of the complex number is infinite. + * + * <p>Note: A complex number with at least one infinite part is regarded + * as an infinity (even if its other part is a NaN). + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @return {@code true} if the complex number contains an infinite value. + * @see Double#isInfinite(double) + */ + public static boolean isInfinite(double real, double imaginary) { + return Double.isInfinite(real) || Double.isInfinite(imaginary); + } + + /** + * Returns the + * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html";> + * natural logarithm</a> of the complex number using its real and imaginary parts. + * + * <p>The natural logarithm of \( z \) is unbounded along the real axis and + * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the + * natural logarithm has a branch cut along the negative real axis \( (-infty,0] \). + * Special cases: + * + * <ul> + * <li>{@code z.conj().log() == z.log().conj()}. + * <li>If {@code z} is −0 + i0, returns −∞ + iπ ("divide-by-zero" floating-point operation). + * <li>If {@code z} is +0 + i0, returns −∞ + i0 ("divide-by-zero" floating-point operation). + * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2. + * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation). + * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +∞ + iπ. + * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0. + * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4. + * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4. + * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN. + * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation). + * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN. + * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN. + * </ul> + * + * <p>Implements the formula: + * + * <p>\[ \ln(z) = \ln |z| + i \arg(z) \] + * + * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument. + * + * <p>The implementation is based on the method described in:</p> + * <blockquote> + * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994) + * Implementing complex elementary functions using exception handling. + * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244. + * </blockquote> + * + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @param constructor Terminating consumer for the complex result, used to construct a result of type {@code R}. + * @param <R> the object type produced by the supplied constructor. + * @return The natural logarithm of the complex number. + * @see Math#log(double) + * @see <a href="http://functions.wolfram.com/ElementaryFunctions/Log/";>Log</a> + */ + public static <R> R log(double real, double imaginary, ComplexResult<R> constructor) { + return log(Math::log, HALF, LN_2, real, imaginary, constructor); + } + + /** + * Returns the base 10 + * <a href="http://mathworld.wolfram.com/CommonLogarithm.html";> + * common logarithm</a> of the complex number using its real and imaginary parts. + * + * <p>The common logarithm of \( z \) is unbounded along the real axis and + * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the + * common logarithm has a branch cut along the negative real axis \( (-infty,0] \). + * Special cases are as defined in the log: + * + * <p>Implements the formula: + * + * <p>\[ \log_{10}(z) = \log_{10} |z| + i \arg(z) \] + * + * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument. + * + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @param constructor Terminating consumer for the complex result, used to construct a result of type {@code R}. + * @param <R> the object type produced by the supplied constructor. + * @return The base 10 logarithm of the complex number. + * @see Math#log10(double) + */ + public static <R> R log10(double real, double imaginary, ComplexResult<R> constructor) { + return log(Math::log10, LOG_10E_O_2, LOG10_2, real, imaginary, constructor); + } + + /** + * Returns the logarithm of complex number using its real and + * imaginary parts and the provided function. + * Implements the formula: + * + * <pre> + * log(x + i y) = log(|x + i y|) + i arg(x + i y)</pre> + * + * <p>Warning: The argument {@code logOf2} must be equal to {@code log(2)} using the + * provided log function otherwise scaling using powers of 2 in the case of overflow + * will be incorrect. This is provided as an internal optimisation. + * + * @param log Log function. + * @param logOfeOver2 The log function applied to e, then divided by 2. + * @param logOf2 The log function applied to 2. + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @param constructor Terminating consumer for the complex result, used to construct a result of type {@code R}. + * @param <R> the object type produced by the supplied constructor. + * @return The logarithm of the complex number. + */ + private static <R> R log(DoubleUnaryOperator log, + double logOfeOver2, + double logOf2, + double real, + double imaginary, + ComplexResult<R> constructor) { + + // Handle NaN + if (Double.isNaN(real) || Double.isNaN(imaginary)) { + // Return NaN unless infinite + if (isInfinite(real, imaginary)) { + return constructor.apply(Double.POSITIVE_INFINITY, Double.NaN); + } + return constructor.apply(Double.NaN, Double.NaN); + } + + // Returns the real part: + // log(sqrt(x^2 + y^2)) + // log(x^2 + y^2) / 2 + + // Compute with positive values + double x = Math.abs(real); + double y = Math.abs(imaginary); + + // Find the larger magnitude. + if (x < y) { + final double tmp = x; + x = y; + y = tmp; + } + + if (x == 0) { + // Handle zero: raises the ‘‘divide-by-zero’’ floating-point exception. + return constructor.apply(Double.NEGATIVE_INFINITY, + Complex.negative(real) ? Math.copySign(Math.PI, imaginary) : imaginary); + } + + double re; + + // This alters the implementation of Hull et al (1994) which used a standard + // precision representation of |z|: sqrt(x*x + y*y). + // This formula should use the same definition of the magnitude returned + // by Complex.abs() which is a high precision computation with scaling. + // The checks for overflow thus only require ensuring the output of |z| + // will not overflow or underflow. + + if (x > HALF && x < ROOT2) { + // x^2+y^2 close to 1. Use log1p(x^2+y^2 - 1) / 2. + re = Math.log1p(Complex.x2y2m1(x, y)) * logOfeOver2; + } else { + // Check for over/underflow in |z| + // When scaling: + // log(a / b) = log(a) - log(b) + // So initialize the result with the log of the scale factor. + re = 0; + if (x > Double.MAX_VALUE / 2) { + // Potential overflow. + if (Complex.isPosInfinite(x)) { + // Handle infinity + return constructor.apply(x, arg(real, imaginary)); + } + // Scale down. + x /= 2; + y /= 2; + // log(2) + re = logOf2; + } else if (y < Double.MIN_NORMAL) { + // Potential underflow. + if (y == 0) { + // Handle real only number + return constructor.apply(log.applyAsDouble(x), arg(real, imaginary)); + } + // Scale up sub-normal numbers to make them normal by scaling by 2^54, + // i.e. more than the mantissa digits. + x *= 0x1.0p54; + y *= 0x1.0p54; + // log(2^-54) = -54 * log(2) + re = -54 * logOf2; + } + re += log.applyAsDouble(abs(x, y)); + } + + // All ISO C99 edge cases for the imaginary are satisfied by the Math library. + return constructor.apply(re, arg(real, imaginary)); + } + + /** + * Returns {@code sqrt(x^2 + y^2)} without intermediate overflow or underflow. + * + * <p>Special cases: + * <ul> + * <li>If either argument is infinite, then the result is positive infinity. + * <li>If either argument is NaN and neither argument is infinite, then the result is NaN. + * </ul> + * + * <p>The computed result is expected to be within 1 ulp of the exact result. + * + * <p>This method is a replacement for {@link Math#hypot(double, double)}. There + * will be differences between this method and {@code Math.hypot(double, double)} due + * to the use of a different algorithm to compute the high precision sum of + * {@code x^2 + y^2}. This method has been tested to have a lower maximum error from + * the exact result; any differences are expected to be 1 ULP indicating a rounding + * change in the sum. + * + * <p>JDK9 ported the hypot function to Java for bug JDK-7130085 due to the slow performance + * of the method as a native function. Benchmarks of the Complex class for functions that + * use hypot confirm this is slow pre-Java 9. This implementation outperforms the new faster + * {@code Math.hypot(double, double)} on JDK 11 (LTS). See the Commons numbers examples JMH + * module for benchmarks. Comparisons with alternative implementations indicate + * performance gains are related to edge case handling and elimination of an unpredictable + * branch in the computation of {@code x^2 + y^2}. + * + * <p>This port was adapted from the "Freely Distributable Math Library" hypot function. + * This method only (and not invoked methods within) is distributed under the terms of the + * original notice as shown below: + * <pre> + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * </pre> + * + * <p>Note: The fdlibm c code makes use of the language ability to read and write directly + * to the upper and lower 32-bits of the 64-double. The function performs + * checking on the upper 32-bits for the magnitude of the two numbers by accessing + * the exponent and 20 most significant bits of the mantissa. These upper bits + * are manipulated during scaling and then used to perform extended precision + * computation of the sum {@code x^2 + y^2} where the high part of the number has 20-bit + * precision. Manipulation of direct bits has no equivalent in Java + * other than use of {@link Double#doubleToLongBits(double)} and + * {@link Double#longBitsToDouble(long)}. To avoid conversion to and from long and double + * representations this implementation only scales the double representation. The high + * and low parts of a double for the extended precision computation are extracted + * using the method of Dekker (1971) to create two 26-bit numbers. This works for sub-normal + * numbers and reduces the maximum error in comparison to fdlibm hypot which does not + * use a split number algorithm for sub-normal numbers. + * + * @param x Value x + * @param y Value y + * @return sqrt(x^2 + y^2) + * @see Math#hypot(double, double) + * @see <a href="https://www.netlib.org/fdlibm/e_hypot.c";>fdlibm e_hypot.c</a> + * @see <a href="https://bugs.java.com/bugdatabase/view_bug.do?bug_id=7130085";>JDK-7130085 : Port fdlibm hypot to Java</a> + */ + private static double hypot(double x, double y) { + // Differences to the fdlibm reference: + // + // 1. fdlibm orders the two parts using the magnitude of the upper 32-bits. + // This incorrectly orders numbers which differ only in the lower 32-bits. + // This invalidates hypot(x, y) = hypot(y, x) for small sub-normal numbers and a minority + // of cases of normal numbers. This implementation forces the |x| >= |y| order + // using the entire 63-bits of the unsigned doubles to ensure the function + // is commutative. + // + // 2. fdlibm computed scaling by directly writing changes to the exponent bits + // and maintained the high part (ha) during scaling for use in the high + // precision sum x^2 + y^2. Since exponent scaling cannot be applied to sub-normals + // the original version dropped the split number representation for sub-normals + // and can produce maximum errors above 1 ULP for sub-normal numbers. + // This version uses Dekker's method to split the number. This can be applied to + // sub-normals and allows dropping the condition to check for sub-normal numbers + // since all small numbers are handled with a single scaling factor. + // The effect is increased precision for the majority of sub-normal cases where + // the implementations compute a different result. + // + // 3. An alteration is done here to add an 'else if' instead of a second + // 'if' statement. Thus you cannot scale down and up at the same time. + // + // 4. There is no use of the absolute double value. The magnitude comparison is + // performed using the long bit representation. The computation x^2+y^2 is + // insensitive to the sign bit. Thus use of Math.abs(double) is only in edge-case + // branches. + // + // 5. The exponent different to ignore the smaller component has changed from 60 to 54. Review Comment: Correct original typo to `difference` ########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexUnaryOperator.java: ########## @@ -0,0 +1,44 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.numbers.complex; + +/** + * Represents a unary operation on a Cartesian form of a complex number \( a + ib \) + * where \( a \) and \( b \) are real numbers represented as two {@code double} + * parts. The operation creates a complex number result; the result is supplied + * to a terminating consumer function which may return an object representation + * of the complex result. + * + * <p>This is a functional interface whose functional method is + * {@link #apply(double, double, ComplexResult)}. + * + * @param <R> The type of the complex result + * @since 1.1 + */ +@FunctionalInterface +public interface ComplexUnaryOperator<R> { + + /** + * Represents an operator that accepts real and imaginary parts of a complex number and a complex constructor to produce and return the result. Review Comment: `and supplies the complex result to the provided consumer` ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java: ########## @@ -145,6 +145,39 @@ static void assertEquals(Supplier<String> msg, double expected, double actual, l } } + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on Review Comment: I would the describe the equality to the expected result (which has further conditions), then state the two operations should be equal. ``` Assert the result is equal to the expected value. ... <p>Assert the operation on the complex number is <em>exactly</em> equal to the operation on the complex real and imaginary parts. ``` ########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexUnaryOperator.java: ########## @@ -0,0 +1,44 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.numbers.complex; + +/** + * Represents a unary operation on a Cartesian form of a complex number \( a + ib \) + * where \( a \) and \( b \) are real numbers represented as two {@code double} + * parts. The operation creates a complex number result; the result is supplied + * to a terminating consumer function which may return an object representation + * of the complex result. + * + * <p>This is a functional interface whose functional method is + * {@link #apply(double, double, ComplexResult)}. + * + * @param <R> The type of the complex result + * @since 1.1 + */ +@FunctionalInterface +public interface ComplexUnaryOperator<R> { + + /** + * Represents an operator that accepts real and imaginary parts of a complex number and a complex constructor to produce and return the result. + * @param r Real part \( a \) of the complex number \(a +ib \). + * @param i Imaginary part \( b \) of the complex number \(a +ib \). + * @param out Terminating consumer for the complex result, used to construct a result of type {@code R}. + * @return the object created by the supplied constructor. Review Comment: ``` @param out Consumer for the complex result @return the object returned by the provided consumer ``` ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java: ########## @@ -300,6 +300,103 @@ private static void assertConjugateEquality(Complex z, } } + /** + * Assert the operation on the complex number satisfies the conjugate equality. + * Assert the operation on the complex number is exactly equal to the operation on Review Comment: Move this before the `<h2>` section, i.e. describe the conjugate equality test then describe the requirement for exact equality between the two operations. `<em>exactly</em>` ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java: ########## @@ -300,6 +300,103 @@ private static void assertConjugateEquality(Complex z, } } + /** + * Assert the operation on the complex number satisfies the conjugate equality. + * Assert the operation on the complex number is exactly equal to the operation on + * complex real and imaginary parts. + * <pre> + * op(conj(z)) = conj(op(z)) + * </pre> + * + * <p>The results must be binary equal. This includes the sign of zero. + * + * <h2>ISO C99 equalities</h2> + * + * <p>Note that this method currently enforces the conjugate equalities for some cases + * where the sign of the real/imaginary parts are unspecified in ISO C99. This is + * allowed (since they are unspecified). The sign specification is appropriately + * handled during testing of odd/even functions. There are some functions where it + * is not possible to satisfy the conjugate equality and also the odd/even rule. + * The compromise made here is to satisfy only one and the other is allowed to fail + * only on the sign of the output result. Known functions where this applies are: + * + * <ul> + * <li>asinh(NaN, inf) + * <li>atanh(NaN, inf) + * <li>cosh(NaN, 0.0) + * <li>sinh(inf, inf) + * <li>sinh(inf, nan) + * </ul> + * + * @param operation the operation + */ + private static void assertConjugateEquality(UnaryOperator<Complex> operation, + ComplexUnaryOperator<ComplexNumber> operation2) { + // Edge cases. Inf/NaN are specifically handled in the C99 test cases + // but are repeated here to enforce the conjugate equality even when the C99 + // standard does not specify a sign. This may be revised in the future. + final double[] parts = {Double.NEGATIVE_INFINITY, -1, -0.0, 0.0, 1, + Double.POSITIVE_INFINITY, Double.NaN}; + for (final double x : parts) { + for (final double y : parts) { + // No conjugate for imaginary NaN + if (!Double.isNaN(y)) { + assertConjugateEquality(complex(x, y), operation, operation2, UnspecifiedSign.NONE); + } + } + } + // Random numbers + final UniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create(); + for (int i = 0; i < 100; i++) { + final double x = next(rng); + final double y = next(rng); + assertConjugateEquality(complex(x, y), operation, operation2, UnspecifiedSign.NONE); + } + } + + /** + * Assert the operation on the complex number satisfies the conjugate equality. + * Assert the operation on the complex number is exactly equal to the operation on Review Comment: Again move this to after the description of conjugate equality. ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java: ########## @@ -300,6 +300,103 @@ private static void assertConjugateEquality(Complex z, } } + /** + * Assert the operation on the complex number satisfies the conjugate equality. + * Assert the operation on the complex number is exactly equal to the operation on + * complex real and imaginary parts. + * <pre> + * op(conj(z)) = conj(op(z)) + * </pre> + * + * <p>The results must be binary equal. This includes the sign of zero. + * + * <h2>ISO C99 equalities</h2> + * + * <p>Note that this method currently enforces the conjugate equalities for some cases + * where the sign of the real/imaginary parts are unspecified in ISO C99. This is + * allowed (since they are unspecified). The sign specification is appropriately + * handled during testing of odd/even functions. There are some functions where it + * is not possible to satisfy the conjugate equality and also the odd/even rule. + * The compromise made here is to satisfy only one and the other is allowed to fail + * only on the sign of the output result. Known functions where this applies are: + * + * <ul> + * <li>asinh(NaN, inf) + * <li>atanh(NaN, inf) + * <li>cosh(NaN, 0.0) + * <li>sinh(inf, inf) + * <li>sinh(inf, nan) + * </ul> + * + * @param operation the operation + */ + private static void assertConjugateEquality(UnaryOperator<Complex> operation, + ComplexUnaryOperator<ComplexNumber> operation2) { + // Edge cases. Inf/NaN are specifically handled in the C99 test cases + // but are repeated here to enforce the conjugate equality even when the C99 + // standard does not specify a sign. This may be revised in the future. + final double[] parts = {Double.NEGATIVE_INFINITY, -1, -0.0, 0.0, 1, + Double.POSITIVE_INFINITY, Double.NaN}; + for (final double x : parts) { + for (final double y : parts) { + // No conjugate for imaginary NaN + if (!Double.isNaN(y)) { + assertConjugateEquality(complex(x, y), operation, operation2, UnspecifiedSign.NONE); + } + } + } + // Random numbers + final UniformRandomProvider rng = RandomSource.SPLIT_MIX_64.create(); + for (int i = 0; i < 100; i++) { + final double x = next(rng); + final double y = next(rng); + assertConjugateEquality(complex(x, y), operation, operation2, UnspecifiedSign.NONE); + } + } + + /** + * Assert the operation on the complex number satisfies the conjugate equality. + * Assert the operation on the complex number is exactly equal to the operation on + * complex real and imaginary parts. + * <pre> + * op(conj(z)) = conj(op(z)) + * </pre> + * + * <p>The results must be binary equal; the sign of the complex number is first processed + * using the provided sign specification. + * + * @param z the complex number + * @param operation1 the operation on the Complex object + * @param operation2 the operation on the complex real and imaginary parts + * @param sign the sign specification + */ + private static void assertConjugateEquality(Complex z, + UnaryOperator<Complex> operation1, + ComplexUnaryOperator<ComplexNumber> operation2, + UnspecifiedSign sign) { + + final Complex zConj = z.conj(); + final Complex c1 = operation1.apply(zConj); + final Complex c2 = operation1.apply(z).conj(); + + final ComplexNumber cn1 = operation2.apply(zConj.getReal(), zConj.getImaginary(), ComplexNumber::new); Review Comment: I would move the assertions using operation2 after the check for conjugate equality. Currently it is confusing to do the conjugate operations, then do them again and test they are the same, then test the output of the first conjugate operations. ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java: ########## @@ -508,6 +605,107 @@ private static void assertComplex(Complex z, } } + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on Review Comment: Move the statement about the second operation to the end of the javadoc. ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CStandardTest.java: ########## @@ -508,6 +605,107 @@ private static void assertComplex(Complex z, } } + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on + * complex real and imaginary parts. If the imaginary part is not NaN the operation + * must also satisfy the conjugate equality. + * + * <pre> + * op(conj(z)) = conj(op(z)) + * </pre> + * + * <p>The results must be binary equal. This includes the sign of zero. + * @param z the complex + * @param operation1 the operation on the Complex object + * @param operation2 the operation on the complex real and imaginary parts + * @param expected the expected + */ + private static void assertComplex(Complex z, + UnaryOperator<Complex> operation1, + ComplexUnaryOperator<ComplexNumber> operation2, + Complex expected) { + + assertComplex(z, operation1, operation2, expected, FunctionType.NONE, UnspecifiedSign.NONE); + } + + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on Review Comment: Again, move the statement about the second operation to the end. ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/CReferenceTest.java: ########## @@ -145,6 +145,39 @@ static void assertEquals(Supplier<String> msg, double expected, double actual, l } } + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on + * complex real and imaginary parts. + * + * <p>The results are considered equal within the provided units of least + * precision. The maximum count of numbers allowed between the two values is + * {@code maxUlps - 1}. + * + * <p>Numbers must have the same sign. Thus -0.0 and 0.0 are never equal. + * + * @param c Input number. + * @param name the operation name + * @param operation1 the operation on the Complex object + * @param operation2 the operation on the complex real and imaginary parts + * @param expected Expected result. + */ + static void assertComplex(Complex c, + String name, + UnaryOperator<Complex> operation1, + ComplexUnaryOperator<ComplexNumber> operation2, + Complex expected, long maxUlps) { + + final Complex z = operation1.apply(c); + final ComplexNumber y = operation2.apply(c.getReal(), c.getImaginary(), ComplexNumber::new); Review Comment: Do not create y until you need it. So either test the exact equality first, or test the equals to the expected and then generate y. The later is a fail fast approach so I prefer that. ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexTest.java: ########## @@ -1419,6 +1419,9 @@ void testLog10() { final Complex z = Complex.ofCartesian(rng.nextDouble() * 2, rng.nextDouble() * 2); final Complex lnz = z.log(); final Complex log10z = z.log10(); + final ComplexNumber log10z2 = ComplexFunctions.log10(z.getReal(), z.getImaginary(), ComplexNumber::new); Review Comment: There is no test for ComplexFunctions.log. In this case it may be convenient to add a method to TestUtils to assert two operations are the same, then return the result. You can then do further assertions on the result. ```Java static Complex assertSame(Complex c, UnaryOperator<Complex> operation1, ComplexUnaryOperator<ComplexNumber> operation2, String name) { final Complex z = operation1.apply(c); // Test operation2 produces the exact same result operation2.apply(c.real(), c.imag(), (x, y) -> { Assertions.assertEquals(z.real(), x, () -> name + " real"); Assertions.assertEquals(z.imag(), y, () -> name + " imaginary"); return null; }); return z; } Complex lnz = TestUtils.assertSame(z, Complex::log, ComplexFunctions::log, "log"); Complex ln10z = TestUtils.assertSame(z, Complex::log10, ComplexFunctions::log10, "log10"); // Test lnz ... ``` ########## commons-numbers-complex/src/test/java/org/apache/commons/numbers/complex/ComplexEdgeCaseTest.java: ########## @@ -81,6 +81,55 @@ private static void assertComplex(double a, double b, CReferenceTest.assertComplex(c, name, operation, e, maxUlps); } + /** + * Assert the operation on the complex number is equal to the expected value. + * Assert the operation on the complex number is exactly equal to the operation on Review Comment: Move statement about operation2 to the end. ########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java: ########## @@ -0,0 +1,584 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.numbers.complex; + +import java.util.function.DoubleUnaryOperator; + +/** + * Cartesian representation of a complex number. The complex number is expressed + * in the form \( a + ib \) where \( a \) and \( b \) are real numbers and \( i \) + * is the imaginary unit which satisfies the equation \( i^2 = -1 \). For the + * complex number \( a + ib \), \( a \) is called the <em>real part</em> and + * \( b \) is called the <em>imaginary part</em>. + * + * <p>This class is immutable, it is to change from an OO representation of a Complex number Review Comment: The class is a utility class with only static methods so immutability is not relevant. I think just remove this paragraph. ########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/ComplexFunctions.java: ########## @@ -0,0 +1,584 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.numbers.complex; + +import java.util.function.DoubleUnaryOperator; + +/** + * Cartesian representation of a complex number. The complex number is expressed + * in the form \( a + ib \) where \( a \) and \( b \) are real numbers and \( i \) + * is the imaginary unit which satisfies the equation \( i^2 = -1 \). For the + * complex number \( a + ib \), \( a \) is called the <em>real part</em> and + * \( b \) is called the <em>imaginary part</em>. + * + * <p>This class is immutable, it is to change from an OO representation of a Complex number + * to a representation as its real and imaginary parts. All arithmetic will create a new instance for the + * result.</p> + * + * <p>Arithmetic in this class conforms to the C99 standard for complex numbers + * defined in ISO/IEC 9899, Annex G. Methods have been named using the equivalent + * method in ISO C99. The behavior for special cases is listed as defined in C99.</p> + * + * <p>For functions \( f \) which obey the conjugate equality \( conj(f(z)) = f(conj(z)) \), + * the specifications for the upper half-plane imply the specifications for the lower + * half-plane.</p> + * + * <p>For functions that are either odd, \( f(z) = -f(-z) \), or even, \( f(z) = f(-z) \), + * the specifications for the first quadrant imply the specifications for the other three + * quadrants.</p> + * + * <p>Special cases of <a href="http://mathworld.wolfram.com/BranchCut.html";>branch cuts</a> + * for multivalued functions adopt the principle value convention from C99. Specials cases + * from C99 that raise the "invalid" or "divide-by-zero" + * <a href="https://en.cppreference.com/w/c/numeric/fenv/FE_exceptions";>floating-point + * exceptions</a> return the documented value without an explicit mechanism to notify + * of the exception case, that is no exceptions are thrown during computations in-line with + * the convention of the corresponding single-valued functions in + * {@link java.lang.Math java.lang.Math}. + * These cases are documented in the method special cases as "invalid" or "divide-by-zero" + * floating-point operation. + * Note: Invalid floating-point exception cases will result in a complex number where the + * cardinality of NaN component parts has increased as a real or imaginary part could + * not be computed and is set to NaN. + * + * @see <a href="http://www.open-std.org/JTC1/SC22/WG14/www/standards";> + * ISO/IEC 9899 - Programming languages - C</a> + */ +public final class ComplexFunctions { + + /** Mask to remove the sign bit from a long. */ + static final long UNSIGN_MASK = 0x7fff_ffff_ffff_ffffL; + + /** Natural logarithm of 2 (ln(2)). */ + private static final double LN_2 = Math.log(2); + /** {@code 1/2}. */ + private static final double HALF = 0.5; + /** Base 10 logarithm of 10 divided by 2 (log10(e)/2). */ + private static final double LOG_10E_O_2 = Math.log10(Math.E) / 2; + /** Base 10 logarithm of 2 (log10(2)). */ + private static final double LOG10_2 = Math.log10(2); + /** {@code sqrt(2)}. */ + private static final double ROOT2 = 1.4142135623730951; + + /** 54 shifted 20-bits to align with the exponent of the upper 32-bits of a double. */ + private static final int EXP_54 = 0x36_00000; + /** Represents an exponent of 500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */ + private static final int EXP_500 = 0x5f3_00000; + /** Represents an exponent of 1024 in unbiased form (infinite or nan) + * shifted 20-bits to align with the upper 32-bits of a double. */ + private static final int EXP_1024 = 0x7ff_00000; + /** Represents an exponent of -500 in unbiased form shifted 20-bits to align with the upper 32-bits of a double. */ + private static final int EXP_NEG_500 = 0x20b_00000; + /** 2^600. */ + private static final double TWO_POW_600 = 0x1.0p+600; + /** 2^-600. */ + private static final double TWO_POW_NEG_600 = 0x1.0p-600; + + /** + * Private constructor for utility class. + */ + private ComplexFunctions() { + } + + /** + * Returns the absolute value of the complex number. + * <pre>abs(x + i y) = sqrt(x^2 + y^2)</pre> + * + * <p>This should satisfy the special cases of the hypot function in ISO C99 F.9.4.3: + * "The hypot functions compute the square root of the sum of the squares of x and y, + * without undue overflow or underflow." + * + * <ul> + * <li>hypot(x, y), hypot(y, x), and hypot(x, −y) are equivalent. + * <li>hypot(x, ±0) is equivalent to |x|. + * <li>hypot(±∞, y) returns +∞, even if y is a NaN. + * </ul> + * + * <p>This method is called by all methods that require the absolute value of the complex + * number, e.g. abs(), sqrt() and log(). + * + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @return The absolute value. + */ + public static double abs(double real, double imaginary) { + // Specialised implementation of hypot. + // See NUMBERS-143 + return hypot(real, imaginary); + } + + /** + * Returns the argument of the complex number. + * + * <p>The argument is the angle phi between the positive real axis and + * the point representing this number in the complex plane. + * The value returned is between \( -\pi \) (not inclusive) + * and \( \pi \) (inclusive), with negative values returned for numbers with + * negative imaginary parts. + * + * <p>If either real or imaginary part (or both) is NaN, then the result is NaN. + * Infinite parts are handled as {@linkplain Math#atan2} handles them, + * essentially treating finite parts as zero in the presence of an + * infinite coordinate and returning a multiple of \( \frac{\pi}{4} \) depending on + * the signs of the infinite parts. + * + * <p>This code follows the + * <a href="http://www.iso-9899.info/wiki/The_Standard";>ISO C Standard</a>, Annex G, + * in calculating the returned value using the {@code atan2(y, x)} method for complex + * \( x + iy \). + * + * @param r Real part \( a \) of the complex number \(a +ib \). + * @param i Imaginary part \( b \) of the complex number \(a +ib \). + * @return The argument of the complex number. + * @see Math#atan2(double, double) + */ + public static double arg(double r, double i) { + // Delegate + return Math.atan2(i, r); + } + + /** + * Returns {@code true} if either real or imaginary component of the complex number is infinite. + * + * <p>Note: A complex number with at least one infinite part is regarded + * as an infinity (even if its other part is a NaN). + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @return {@code true} if the complex number contains an infinite value. + * @see Double#isInfinite(double) + */ + public static boolean isInfinite(double real, double imaginary) { + return Double.isInfinite(real) || Double.isInfinite(imaginary); + } + + /** + * Returns the + * <a href="http://mathworld.wolfram.com/NaturalLogarithm.html";> + * natural logarithm</a> of the complex number using its real and imaginary parts. + * + * <p>The natural logarithm of \( z \) is unbounded along the real axis and + * in the range \( [-\pi, \pi] \) along the imaginary axis. The imaginary part of the + * natural logarithm has a branch cut along the negative real axis \( (-infty,0] \). + * Special cases: + * + * <ul> + * <li>{@code z.conj().log() == z.log().conj()}. + * <li>If {@code z} is −0 + i0, returns −∞ + iπ ("divide-by-zero" floating-point operation). + * <li>If {@code z} is +0 + i0, returns −∞ + i0 ("divide-by-zero" floating-point operation). + * <li>If {@code z} is x + i∞ for finite x, returns +∞ + iπ/2. + * <li>If {@code z} is x + iNaN for finite x, returns NaN + iNaN ("invalid" floating-point operation). + * <li>If {@code z} is −∞ + iy for finite positive-signed y, returns +∞ + iπ. + * <li>If {@code z} is +∞ + iy for finite positive-signed y, returns +∞ + i0. + * <li>If {@code z} is −∞ + i∞, returns +∞ + i3π/4. + * <li>If {@code z} is +∞ + i∞, returns +∞ + iπ/4. + * <li>If {@code z} is ±∞ + iNaN, returns +∞ + iNaN. + * <li>If {@code z} is NaN + iy for finite y, returns NaN + iNaN ("invalid" floating-point operation). + * <li>If {@code z} is NaN + i∞, returns +∞ + iNaN. + * <li>If {@code z} is NaN + iNaN, returns NaN + iNaN. + * </ul> + * + * <p>Implements the formula: + * + * <p>\[ \ln(z) = \ln |z| + i \arg(z) \] + * + * <p>where \( |z| \) is the absolute and \( \arg(z) \) is the argument. + * + * <p>The implementation is based on the method described in:</p> + * <blockquote> + * T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang (1994) + * Implementing complex elementary functions using exception handling. + * ACM Transactions on Mathematical Software, Vol 20, No 2, pp 215-244. + * </blockquote> + * + * @param real Real part \( a \) of the complex number \(a +ib \). + * @param imaginary Imaginary part \( b \) of the complex number \(a +ib \). + * @param constructor Terminating consumer for the complex result, used to construct a result of type {@code R}. Review Comment: Rename `constructor` to `action`. Move the javadoc description from the return tag to this tag. The return tag can be made more generic as you do not know what the object will represent. ``` @param action Consumer for the natural logarithm of the complex number @return the object returned by the supplied action ... ``` ########## commons-numbers-complex/src/main/java/org/apache/commons/numbers/complex/Complex.java: ########## @@ -2904,7 +2782,8 @@ private static Complex atanh(final double real, final double imaginary, * @param y the y value * @return {@code x^2 + y^2 - 1}. */ - private static double x2y2m1(double x, double y) { + //TODO - make it private in future Review Comment: Since the intent is to eventually move all these private functions to ComplexFunctions, then please start the process now. Make any function that is currently also required in Complex package-private and add the TODO note to make private in CompexFunctions in the future. It may be convenient to leave the functions locally as pass through functions: ```java private static double x2y2m1(double x, double y) { return ComplexFunctions.x2y2m1(x, y); } ``` This will prevent having to update all the call sites in Complex and minimise the diff. -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. To unsubscribe, e-mail: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
