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 594ffaddb28b553558138c2606924b2182b84d91 Author: aherbert <[email protected]> AuthorDate: Thu Nov 7 18:21:07 2019 +0000 [NUMBERS-78] Increase test coverage of divide/multiply. This involves fixing the edge case correction algorithm. Multiply has been updated to explicitly detect infinities multiplied by zero and not correct the NaN,NaN result. Divide algorithm has been updated to handle overflow in infinite/finite and finite/infinite correction. --- .../apache/commons/numbers/complex/Complex.java | 210 +++++++++++++----- .../commons/numbers/complex/CStandardTest.java | 234 +++++++++++++++++++++ 2 files changed, 391 insertions(+), 53 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 083f774..3a86fd3 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 @@ -354,7 +354,7 @@ public final class Complex implements Serializable { double d = divisor.getImaginary(); int ilogbw = 0; final double logbw = Math.log(Math.max(Math.abs(c), Math.abs(d))) / Math.log(2); - if (!Double.isInfinite(logbw)) { + if (Double.isFinite(logbw)) { ilogbw = (int)logbw; c = Math.scalb(c, -ilogbw); d = Math.scalb(d, -ilogbw); @@ -362,29 +362,83 @@ public final class Complex implements Serializable { final double denom = c * c + d * d; double x = Math.scalb((a * c + b * d) / denom, -ilogbw); double y = Math.scalb((b * c - a * d) / denom, -ilogbw); + // Recover infinities and zeros that computed as NaN+iNaN + // the only cases are nonzero/zero, infinite/finite, and finite/infinite, ... + // -------------- + // Modification from the listing in ISO C.99 G.5.1 (8): + // Prevent overflow in (a * c + b * d) and (b * c - a * d). + // It is only the sign that is important. not the magnitude. + // -------------- if (Double.isNaN(x) && Double.isNaN(y)) { if ((denom == 0.0) && (!Double.isNaN(a) || !Double.isNaN(b))) { + // nonzero/zero x = Math.copySign(Double.POSITIVE_INFINITY, c) * a; y = Math.copySign(Double.POSITIVE_INFINITY, c) * b; - } else if ((Double.isInfinite(a) && Double.isInfinite(b)) && - !Double.isInfinite(c) && !Double.isInfinite(d)) { - a = Math.copySign(Double.isInfinite(a) ? 1.0 : 0.0, a); - b = Math.copySign(Double.isInfinite(b) ? 1.0 : 0.0, b); - x = Double.POSITIVE_INFINITY * (a * c + b * d); - y = Double.POSITIVE_INFINITY * (b * c - a * d); - } else if (Double.isInfinite(logbw) && - !Double.isInfinite(a) && !Double.isInfinite(b)) { - c = Math.copySign(Double.isInfinite(c) ? 1.0 : 0.0, c); - d = Math.copySign(Double.isInfinite(d) ? 1.0 : 0.0, d); - x = 0.0 * (a * c + b * d); - y = 0.0 * (b * c - a * d); + } else if ((Double.isInfinite(a) || Double.isInfinite(b)) && + Double.isFinite(c) && Double.isFinite(d)) { + // infinite/finite + a = boxInfinity(a); + b = boxInfinity(b); + x = Double.POSITIVE_INFINITY * computeACplusBD(a, b, c, d); + y = Double.POSITIVE_INFINITY * computeBCminusAD(a, b, c, d); + } else if (divisor.isInfinite() && + Double.isFinite(a) && Double.isFinite(b)) { + // finite/infinite + c = boxInfinity(c); + d = boxInfinity(d); + x = 0.0 * computeACplusBD(a, b, c, d); + y = 0.0 * computeBCminusAD(a, b, c, d); } } return new Complex(x, y); } /** + * Compute {@code a*c + b*d} without overflow. + * It is assumed: either {@code a} and {@code b} or {@code c} and {@code d} are + * either zero or one (i.e. a boxed infinity); and the sign of the result is important, + * not the value. + * + * @param a the a + * @param b the b + * @param c the c + * @param d the d + * @return the result + */ + private static double computeACplusBD(double a, double b, double c, double d) { + final double ac = a * c; + final double bd = b * d; + double result = ac + bd; + return Double.isFinite(result) ? + result : + // Overflow. Just divide by 2 as it is the sign of the result that matters. + ac * 0.5 + bd * 0.5; + } + + /** + * Compute {@code b*c - a*d} without overflow. + * It is assumed: either {@code a} and {@code b} or {@code c} and {@code d} are + * either zero or one (i.e. a boxed infinity); and the sign of the result is important, + * not the value. + * + * @param a the a + * @param b the b + * @param c the c + * @param d the d + * @return the result + */ + private static double computeBCminusAD(double a, double b, double c, double d) { + final double bc = b * c; + final double ad = a * d; + final double result = bc - ad; + return Double.isFinite(result) ? + result : + // Overflow. Just divide by 2 as it is the sign of the result that matters. + bc * 0.5 + ad * 0.5; + } + + /** * Returns a {@code Complex} whose value is {@code (this / divisor)}, * with {@code divisor} interpreted as a real number. * @@ -612,44 +666,51 @@ public final class Complex implements Serializable { final double bc = b * c; double x = ac - bd; double y = ad + bc; - if (Double.isNaN(a) && Double.isNaN(b)) { + + // -------------- + // NaN can occur if: + // - any of (a,b,c,d) are NaN (for NaN or Infinite complex numbers) + // - a multiplication of infinity by zero (ac,bd,ad,bc). + // - a subtraction of infinity from infinity (e.g. ac - bd) + // Note that (ac,bd,ad,bc) can be infinite due to overflow. + // + // Detect a NaN result and perform correction. + // + // Modification from the listing in ISO C.99 G.5.1 (6) + // Do not correct infinity multiplied by zero. This is left as NaN. + // -------------- + + if (Double.isNaN(x) && Double.isNaN(y)) { + // Recover infinities that computed as NaN+iNaN ... boolean recalc = false; - if (Double.isInfinite(a) || Double.isInfinite(b)) { - a = Math.copySign(Double.isInfinite(a) ? 1.0 : 0.0, a); - b = Math.copySign(Double.isInfinite(a) ? 1.0 : 0.0, a); - if (Double.isNaN(c)) { - c = Math.copySign(0.0, c); - } - if (Double.isNaN(d)) { - d = Math.copySign(0.0, d); - } + if (isInfinite() && isNotZero(c, d)) { + // This complex is infinite. + // "Box" the infinity and change NaNs in the other factor to 0. + a = boxInfinity(a); + b = boxInfinity(b); + c = changeNaNtoZero(c); + d = changeNaNtoZero(d); recalc = true; } - if (Double.isInfinite(c) || Double.isInfinite(d)) { - c = Math.copySign(Double.isInfinite(c) ? 1.0 : 0.0, c); - d = Math.copySign(Double.isInfinite(d) ? 1.0 : 0.0, d); - if (Double.isNaN(a)) { - a = Math.copySign(0.0, a); - } - if (Double.isNaN(b)) { - b = Math.copySign(0.0, b); - } + // (c, d) may have been corrected so do not use factor.isInfinite(). + if ((Double.isInfinite(c) || Double.isInfinite(d)) && + isNotZero(a, b)) { + // This other complex is infinite. + // "Box" the infinity and change NaNs in the other factor to 0. + c = boxInfinity(c); + d = boxInfinity(d); + a = changeNaNtoZero(a); + b = changeNaNtoZero(b); recalc = true; } if (!recalc && (Double.isInfinite(ac) || Double.isInfinite(bd) || - Double.isInfinite(ad) || Double.isInfinite(bc))) { - if (Double.isNaN(a)) { - a = Math.copySign(0.0, a); - } - if (Double.isNaN(b)) { - b = Math.copySign(0.0, b); - } - if (Double.isNaN(c)) { - c = Math.copySign(0.0, c); - } - if (Double.isNaN(d)) { - d = Math.copySign(0.0, d); - } + Double.isInfinite(ad) || Double.isInfinite(bc))) { + // The result overflowed to infinity. + // Recover infinities from overflow by changing NaNs to 0 ... + a = changeNaNtoZero(a); + b = changeNaNtoZero(b); + c = changeNaNtoZero(c); + d = changeNaNtoZero(d); recalc = true; } if (recalc) { @@ -661,6 +722,49 @@ public final class Complex implements Serializable { } /** + * Box values for the real or imaginary component of an infinite complex number. + * Any infinite value will be returned as one. Non-infinite values will be returned as zero. + * The sign is maintained. + * + * <pre> + * inf = 1 + * -inf = -1 + * x = 0 + * -x = -0 + * </pre> + * + * @param component the component + * @return the boxed value + */ + private static double boxInfinity(double component) { + return Math.copySign(Double.isInfinite(component) ? 1.0 : 0.0, component); + } + + /** + * Checks if the complex number is not zero. + * + * @param real the real component + * @param imaginary the imaginary component + * @return true if the complex is zero + */ + private static boolean isNotZero(double real, double imaginary) { + // The use of equals is deliberate. + // This method must distinguish NaN from zero thus ruling out: + // (real != 0.0 || imaginary != 0.0) + return !(real == 0.0 && imaginary == 0.0); + } + + /** + * Change NaN to zero preserving the sign; otherwise return the value. + * + * @param value the value + * @return the new value + */ + private static double changeNaNtoZero(double value) { + return Double.isNaN(value) ? Math.copySign(0.0, value) : value; + } + + /** * Returns a {@code Complex} whose value is {@code this * factor}, with {@code factor} * interpreted as a integer number. * @@ -830,7 +934,7 @@ public final class Complex implements Serializable { imaginary == Double.POSITIVE_INFINITY) { return new Complex(Double.POSITIVE_INFINITY, PI_OVER_2); } else if (real == Double.POSITIVE_INFINITY && - !Double.isInfinite(imaginary) && !Double.isNaN(imaginary)) { + Double.isFinite(imaginary)) { return new Complex(Double.POSITIVE_INFINITY, 0); } else if (real == Double.POSITIVE_INFINITY && imaginary == Double.POSITIVE_INFINITY) { @@ -1073,8 +1177,8 @@ public final class Complex implements Serializable { return new Complex(Double.POSITIVE_INFINITY, PI_OVER_4); } else if ((real == Double.POSITIVE_INFINITY && Double.isNaN(imaginary)) || - ((Double.isNaN(real) && - imaginary == Double.POSITIVE_INFINITY))) { + (Double.isNaN(real) && + imaginary == Double.POSITIVE_INFINITY)) { return new Complex(Double.POSITIVE_INFINITY, Double.NaN); } return new Complex(Math.log(abs()), @@ -1095,8 +1199,8 @@ public final class Complex implements Serializable { return new Complex(Double.POSITIVE_INFINITY, PI_OVER_4); } else if ((real == Double.POSITIVE_INFINITY && Double.isNaN(imaginary)) || - ((Double.isNaN(real) && - imaginary == Double.POSITIVE_INFINITY))) { + (Double.isNaN(real) && + imaginary == Double.POSITIVE_INFINITY)) { return new Complex(Double.POSITIVE_INFINITY, Double.NaN); } return new Complex(Math.log(abs()) / Math.log(10), @@ -1479,15 +1583,15 @@ public final class Complex implements Serializable { * <li>it is not infinite,</li> * <li>it is positive signed,</li> * </ul> - * - * <p>Note: This is true for negative zero.</p> + * + * <p>Note: This is false for negative zero.</p> * * @param d Value. * @return {@code true} if {@code d} meets all the conditions and * {@code false} otherwise. */ private static boolean positiveSignedFinite(double d) { - return Double.isFinite(d) && d >= 0; + return (Double.isFinite(d) && d > 0) || Double.doubleToLongBits(d) == 0; } /** See {@link #parse(String)}. */ 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 e336b0c..4a8fe97 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 @@ -23,6 +23,10 @@ import org.apache.commons.rng.simple.RandomSource; import org.junit.jupiter.api.Assertions; import org.junit.jupiter.api.Test; +import java.util.ArrayList; +import java.util.function.BiFunction; +import java.util.function.Predicate; + /** * Tests the standards defined by the C.99 standard for complex numbers * defined in ISO/IEC 9899, Annex G. @@ -36,6 +40,7 @@ public class CStandardTest { private static final double inf = Double.POSITIVE_INFINITY; private static final double negInf = Double.NEGATIVE_INFINITY; private static final double nan = Double.NaN; + private static final double max = Double.MAX_VALUE; private static final double piOverFour = Math.PI / 4.0; private static final double piOverTwo = Math.PI / 2.0; private static final double threePiOverFour = 3.0 * Math.PI / 4.0; @@ -74,6 +79,9 @@ public class CStandardTest { private static final Complex threePiFourNegInf = complex(threePiOverFour, negInf); private static final Complex piFourNegInf = complex(piOverFour, negInf); private static final Complex NAN = complex(nan, nan); + private static final Complex maxMax = complex(max, max); + private static final Complex maxNan = complex(max, nan); + private static final Complex nanMax = complex(nan, max); // CHECKSTYLE: resume ConstantName /** @@ -107,6 +115,24 @@ public class CStandardTest { } /** + * Assert the operation on the two complex numbers. + * + * @param c1 the first complex + * @param c2 the second complex + * @param operation the operation + * @param operationName the operation name + * @param expected the expected + * @param expectedName the expected name + */ + private static void assertOperation(Complex c1, Complex c2, + BiFunction<Complex, Complex, Complex> operation, String operationName, + Predicate<Complex> expected, String expectedName) { + final Complex z = operation.apply(c1, c2); + Assertions.assertTrue(expected.test(z), + () -> String.format("%s expected: %s %s %s = %s", expectedName, c1, operationName, c2, z)); + } + + /** * Utility to create a Complex. * * @param real the real @@ -118,6 +144,203 @@ public class CStandardTest { } /** + * Creates a list of Complex infinites. + * + * @return the list + */ + private static ArrayList<Complex> createInfinites() { + final double[] values = {0, 1, inf, negInf, nan}; + return createCombinations(values, Complex::isInfinite); + } + + /** + * Creates a list of Complex finites that are not zero. + * + * @return the list + */ + private static ArrayList<Complex> createNonZeroFinites() { + final double[] values = {-1, -0, 0, 1, Double.MAX_VALUE}; + return createCombinations(values, c -> !CStandardTest.isZero(c)); + } + + /** + * Creates a list of Complex finites that are zero: [0,0], [-0,0], [0,-0], [-0,-0]. + * + * @return the list + */ + private static ArrayList<Complex> createZeroFinites() { + final double[] values = {-0, 0}; + return createCombinations(values, c -> true); + } + + /** + * Creates a list of Complex NaNs. + * + * @return the list + */ + private static ArrayList<Complex> createNaNs() { + final double[] values = {0, 1, nan}; + return createCombinations(values, Complex::isNaN); + } + + /** + * Creates a list of Complex numbers as an all-vs-all combinations that pass the + * condition. + * + * @param values the values + * @param condition the condition + * @return the list + */ + private static ArrayList<Complex> createCombinations(final double[] values, Predicate<Complex> condition) { + final ArrayList<Complex> list = new ArrayList<>(); + for (double re : values) { + for (double im : values) { + final Complex z = complex(re, im); + if (condition.test(z)) { + list.add(z); + } + } + } + return list; + } + + /** + * Checks if the complex is zero. + * + * @param c the complex + * @return true if zero + */ + private static boolean isZero(Complex c) { + return Complex.equals(c, Complex.ZERO, 0); + } + + /** + * ISO C Standard G.5 (4). + */ + @Test + public void testMultiply() { + final ArrayList<Complex> infinites = createInfinites(); + final ArrayList<Complex> nonZeroFinites = createNonZeroFinites(); + final ArrayList<Complex> zeroFinites = createZeroFinites(); + + // C.99 refers to non-zero finites. + // Standard multiplication of zero with infinites is not defined. + Assertions.assertEquals(nan, 0.0 * inf, "0 * inf"); + Assertions.assertEquals(nan, 0.0 * negInf, "0 * -inf"); + Assertions.assertEquals(nan, -0.0 * inf, "-0 * inf"); + Assertions.assertEquals(nan, -0.0 * negInf, "-0 * -inf"); + + // "if one operand is an infinity and the other operand is a nonzero finite number or an + // infinity, then the result of the * operator is an infinity;" + for (Complex z : infinites) { + for (Complex w : infinites) { + assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, "Inf"); + } + for (Complex w : nonZeroFinites) { + assertOperation(z, w, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(w, z, Complex::multiply, "*", Complex::isInfinite, "Inf"); + } + // C.99 refers to non-zero finites. + // Infer that Complex multiplication of zero with infinites is not defined. + for (Complex w : zeroFinites) { + assertOperation(z, w, Complex::multiply, "*", Complex::isNaN, "NaN"); + assertOperation(w, z, Complex::multiply, "*", Complex::isNaN, "NaN"); + } + } + + // ISO C Standard in Annex G is missing an explicit definition of how to handle NaNs. + // We will assume multiplication by (nan,nan) is not allowed. + // It is undefined how to multiply when a complex has only one NaN component. + // The reference implementation allows it. + + // The GNU g++ compiler computes: + // (1e300 + i 1e300) * (1e30 + i NAN) = inf + i inf + // Thus this is allowing some computations with NaN. + + // Check multiply with (NaN,NaN) is not corrected + final double[] values = {0, 1, inf, negInf, nan}; + for (Complex z : createCombinations(values, c -> true)) { + assertOperation(z, NAN, Complex::multiply, "*", Complex::isNaN, "NaN"); + assertOperation(NAN, z, Complex::multiply, "*", Complex::isNaN, "NaN"); + } + + // Test multiply cases which result in overflow are corrected to infinity + assertOperation(maxMax, maxMax, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(maxNan, maxNan, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(nanMax, maxNan, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(maxNan, nanMax, Complex::multiply, "*", Complex::isInfinite, "Inf"); + assertOperation(nanMax, nanMax, Complex::multiply, "*", Complex::isInfinite, "Inf"); + } + + /** + * ISO C Standard G.5 (4). + */ + @Test + public void testDivide() { + final ArrayList<Complex> infinites = createInfinites(); + final ArrayList<Complex> nonZeroFinites = createNonZeroFinites(); + final ArrayList<Complex> zeroFinites = createZeroFinites(); + final ArrayList<Complex> nans = createNaNs(); + final ArrayList<Complex> finites = new ArrayList<>(nonZeroFinites); + finites.addAll(zeroFinites); + + // "if the first operand is an infinity and the second operand is a finite number, then the + // result of the / operator is an infinity;" + for (Complex z : infinites) { + for (Complex w : nonZeroFinites) { + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); + } + for (Complex w : zeroFinites) { + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); + } + // Check inf/inf cannot be done + for (Complex w : infinites) { + assertOperation(z, w, Complex::divide, "/", Complex::isNaN, "NaN"); + } + } + + // "if the first operand is a finite number and the second operand is an infinity, then the + // result of the / operator is a zero;" + for (Complex z : finites) { + for (Complex w : infinites) { + assertOperation(z, w, Complex::divide, "/", CStandardTest::isZero, "Zero"); + } + } + + // "if the first operand is a nonzero finite number or an infinity and the second operand is + // a zero, then the result of the / operator is an infinity." + for (Complex w : zeroFinites) { + for (Complex z : nonZeroFinites) { + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); + } + for (Complex z : infinites) { + assertOperation(z, w, Complex::divide, "/", Complex::isInfinite, "Inf"); + } + } + + // ISO C Standard in Annex G is missing an explicit definition of how to handle NaNs. + // The reference implementation does not correct for divide by NaN components unless + // infinite. + for (Complex w : nans) { + for (Complex z : finites) { + assertOperation(z, w, Complex::divide, "/", c -> NAN.equals(c), "(NaN,NaN)"); + } + for (Complex z : infinites) { + assertOperation(z, w, Complex::divide, "/", c -> NAN.equals(c), "(NaN,NaN)"); + } + } + + // Check (NaN,NaN) divide is not corrected for the edge case of divide by zero or infinite + for (Complex w : zeroFinites) { + assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, "NaN"); + } + for (Complex w : infinites) { + assertOperation(NAN, w, Complex::divide, "/", Complex::isNaN, "NaN"); + } + } + + /** * ISO C Standard G.6 (3). */ @Test @@ -210,7 +433,9 @@ public class CStandardTest { assertComplex(zeroNaN.acosh(), NAN); assertComplex(oneNaN.acosh(), NAN); assertComplex(negInfOne.acosh(), infPi); + assertComplex(negInfZero.acosh(), infPi); assertComplex(infOne.acosh(), infZero); + assertComplex(infZero.acosh(), infZero); assertComplex(negInfPosInf.acosh(), infThreePiFour); assertComplex(infInf.acosh(), infPiFour); assertComplex(infNaN.acosh(), infNaN); @@ -218,6 +443,15 @@ public class CStandardTest { assertComplex(nanOne.acosh(), NAN); assertComplex(nanInf.acosh(), infNaN); assertComplex(NAN.acosh(), NAN); + // The standard mentions positive-signed y + Assertions.assertThrows(AssertionError.class, () -> { + // Not −∞ + iy, positive signed y + assertComplex(complex(negInf, -0.0).acosh(), infPi); + }); + Assertions.assertThrows(AssertionError.class, () -> { + // Not ∞ + iy, positive signed y + assertComplex(complex(inf, -0.0).acosh(), infPi); + }); } /**
