From 65aac0be46725deac5fddc8494f8c0e8de1a5c92 Mon Sep 17 00:00:00 2001 From: Stephen Canon Date: Fri, 29 Nov 2019 17:24:26 -0500 Subject: [PATCH 1/4] Add test cases to stress pow(_:Float/Double,:Int) These are the hardest cases to get right; most of the fail at present. Also some improvements to sanityCheck to handle negative values correctly and treat infinity as the next binade up (which isn't ideal for general use, but is good for this sort of rough testing.) --- Tests/RealTests/RealTests.swift | 121 +++++++++++++++++++++++++++----- 1 file changed, 102 insertions(+), 19 deletions(-) diff --git a/Tests/RealTests/RealTests.swift b/Tests/RealTests/RealTests.swift index 4a569639..917d031c 100644 --- a/Tests/RealTests/RealTests.swift +++ b/Tests/RealTests/RealTests.swift @@ -1,4 +1,4 @@ -//===--- ElementaryFunctionTests.swift ------------------------*- swift -*-===// +//===--- ElementaryFunctionChecks.swift ------------------------*- swift -*-===// // // This source file is part of the Swift Numerics open source project // @@ -22,11 +22,15 @@ func sanityCheck(_ expected: TestLiteralType, _ actual: T, ulps allowed: T = 16, file: StaticString = #file, line: UInt = #line) where T: BinaryFloatingPoint { + // Shortcut relative-error check if we got the sign wrong; it's OK to + // underflow to zero a little bit early, but we don't want to allow going + // right through zero to the other side. + XCTAssert(actual.sign == expected.sign, "\(actual) != \(expected) as \(T.self)", file: file, line: line) // Default tolerance is 16 ulps; It's OK to relax this as needed for new - // platforms, as these tests are *not* intended to validate the math + // platforms, as these Checks are *not* intended to validate the math // library--they are only intended to check that the Swift bindings are // calling the right functions in the math library. It's important, however - // not to relax the tolerance beyond a few hundred ulps, because these tests + // not to relax the tolerance beyond a few hundred ulps, because these Checks // need to detect errors where the *wrong function* is being called; e.g. // we need to flag an implentation that inadvertently called the C hypotf // function instead of hypot. This is especially important because the C @@ -34,16 +38,26 @@ func sanityCheck(_ expected: TestLiteralType, _ actual: T, if actual == T(expected) || actual.isNaN && expected.isNaN { return } - // Compute error in ulp, compare to tolerance. - let absoluteError = T(abs(TestLiteralType(actual) - expected)) - let ulpError = absoluteError / T(expected).ulp - XCTAssert(ulpError <= allowed, "\(actual) != \(expected) as \(T.self)" + - "\n \(ulpError)-ulp error exceeds \(allowed)-ulp tolerance.", - file: file, line: line) + // Special-case where expected or observed is infinity. + // Artificially knock everything down a binade, treat actual infinity as + // the base of the next binade up. + if actual.isInfinite || T(expected).isInfinite { + let scaledExpected = TestLiteralType(signOf: expected, + magnitudeOf: expected.isInfinite ? TestLiteralType.greatestFiniteMagnitude.binade : 0.5 * expected + ) + let scaledActual = T(signOf: actual, + magnitudeOf: actual.isInfinite ? T.greatestFiniteMagnitude.binade : 0.5 * actual + ) + return sanityCheck(scaledExpected, scaledActual, ulps: allowed, file: file, line: line) + } + // Compute error in ulp, compare to tolerance. + let absoluteError = T(abs(TestLiteralType(actual) - expected)).magnitude + let ulpError = absoluteError / max(T(expected).magnitude, T.leastNormalMagnitude).ulp + XCTAssert(ulpError <= allowed, "\(actual) != \(expected) as \(T.self)\n\(ulpError)-ulp error exceeds \(allowed)-ulp tolerance.", file: file, line: line) } internal extension ElementaryFunctions where Self: BinaryFloatingPoint { - static func elementaryFunctionTests() { + static func elementaryFunctionChecks() { sanityCheck(1.1863995522992575361931268186727044683, Self.acos(0.375)) sanityCheck(0.3843967744956390830381948729670469737, Self.asin(0.375)) sanityCheck(0.3587706702705722203959200639264604997, Self.atan(0.375)) @@ -69,7 +83,7 @@ internal extension ElementaryFunctions where Self: BinaryFloatingPoint { } internal extension Real where Self: BinaryFloatingPoint { - static func realFunctionTests() { + static func realFunctionChecks() { sanityCheck(1.2968395546510096659337541177924511598, Self.exp2(0.375)) sanityCheck(2.3713737056616552616517527574788898386, Self.exp10(0.375)) sanityCheck(-1.415037499278843818546261056052183491, Self.log2(0.375)) @@ -85,25 +99,94 @@ internal extension Real where Self: BinaryFloatingPoint { XCTAssertEqual(.minus, Self.signGamma(-2.375)) #endif } + + static func testPownCommon() { + // If x is -1, then the result is ±1 with sign chosen by parity of n. + // Simply converting n to Real will flip parity when n is large, so + // first check that we get those cases right. + XCTAssertEqual(Self.pow(-1, 0), 1) + XCTAssertEqual(Self.pow(-1, 1), -1) + XCTAssertEqual(Self.pow(-1, -1), -1) + XCTAssertEqual(Self.pow(-1, 2), 1) + XCTAssertEqual(Self.pow(-1, -2), 1) + XCTAssertEqual(Self.pow(-1, Int.max - 1), 1) + XCTAssertEqual(Self.pow(-1, -Int.max + 1), 1) + XCTAssertEqual(Self.pow(-1, Int.max), -1) + XCTAssertEqual(Self.pow(-1, -Int.max), -1) + XCTAssertEqual(Self.pow(-1, Int.min), 1) + } +} + +extension Float { + static func testPown() { + testPownCommon() + let u = Float(1).nextUp + let d = Float(1).nextDown + // Smallest exponents not exactly representable as Float. + sanityCheck(-7.3890560989306677280287919329569359, Float.pow(-u, 0x1000001)) + sanityCheck(-0.3678794082804575860056608283059288, Float.pow(-d, 0x1000001)) + // Exponents close to overflow boundary. + sanityCheck(-3.4028231352500001570898203463449749e38, Float.pow(-u, 744261161)) + sanityCheck( 3.4028235408981285772043562848249166e38, Float.pow(-u, 744261162)) + sanityCheck(-3.4028239465463053543440887892352174e38, Float.pow(-u, 744261163)) + sanityCheck( 3.4028233551634475284795244782720072e38, Float.pow(-d, -1488522190)) + sanityCheck(-3.4028235579875369356575053576685267e38, Float.pow(-d, -1488522191)) + sanityCheck( 3.4028237608116384320940078199368685e38, Float.pow(-d, -1488522192)) + // Exponents close to underflow boundary. + sanityCheck( 7.0064936491761438872280296737844625e-46, Float.pow(-u, -872181048)) + sanityCheck(-7.0064928139371132951305928725186420e-46, Float.pow(-u, -872181049)) + sanityCheck( 7.0064919786981822712727285793333389e-46, Float.pow(-u, -872181050)) + sanityCheck(-7.0064924138100205091278464932003585e-46, Float.pow(-d, 1744361943)) + sanityCheck( 7.0064919961905290625123586120258840e-46, Float.pow(-d, 1744361944)) + sanityCheck(-7.0064915785710625079583096856510544e-46, Float.pow(-d, 1744361945)) + } } -final class ElementaryFunctionTests: XCTestCase { +extension Double { + static func testPown() { + testPownCommon() + let u: Double = 1.nextUp + let d: Double = 1.nextDown + // Smallest exponent not exactly representable as Double. + sanityCheck(-7.3890560989306502272304274605750685, Double.pow(-u, 0x20000000000001)) + sanityCheck(-0.1353352832366126918939994949724833, Double.pow(-u, -0x20000000000001)) + sanityCheck(-0.3678794411714422603312898889458068, Double.pow(-d, 0x20000000000001)) + sanityCheck(-2.7182818284590456880451484776630468, Double.pow(-d, -0x20000000000001)) + // Exponents close to overflow boundary. + sanityCheck( 1.7976931348623151738531864721534215e308, Double.pow(-u, 3196577161300664268)) + sanityCheck(-1.7976931348623155730212483790972209e308, Double.pow(-u, 3196577161300664269)) + sanityCheck( 1.7976931348623159721893102860411089e308, Double.pow(-u, 3196577161300664270)) + sanityCheck( 1.7976931348623157075547244136070910e308, Double.pow(-d, -6393154322601327474)) + sanityCheck(-1.7976931348623159071387553670790721e308, Double.pow(-d, -6393154322601327475)) + sanityCheck( 1.7976931348623161067227863205510754e308, Double.pow(-d, -6393154322601327476)) + // Exponents close to underflow boundary. + sanityCheck( 2.4703282292062334560337346683707907e-324, Double.pow(-u, -3355781687888880946)) + sanityCheck(-2.4703282292062329075106789791206172e-324, Double.pow(-u, -3355781687888880947)) + sanityCheck( 2.4703282292062323589876232898705654e-324, Double.pow(-u, -3355781687888880948)) + sanityCheck(-2.4703282292062332640976590913373022e-324, Double.pow(-d, 6711563375777760775)) + sanityCheck( 2.4703282292062329898361312467121758e-324, Double.pow(-d, 6711563375777760776)) + sanityCheck(-2.4703282292062327155746034020870799e-324, Double.pow(-d, 6711563375777760777)) + } +} + +final class ElementaryFunctionChecks: XCTestCase { func testFloat() { - Float.elementaryFunctionTests() - Float.realFunctionTests() + Float.elementaryFunctionChecks() + Float.realFunctionChecks() + Float.testPown() } func testDouble() { - Double.elementaryFunctionTests() - Double.realFunctionTests() + Double.elementaryFunctionChecks() + Double.realFunctionChecks() + Double.testPown() } - #if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android) func testFloat80() { - Float80.elementaryFunctionTests() - Float80.realFunctionTests() + Float80.elementaryFunctionChecks() + Float80.realFunctionChecks() } #endif } From 5bf2f828ef17f83da84b6e99bc627e97c893fece Mon Sep 17 00:00:00 2001 From: Stephen Canon Date: Fri, 29 Nov 2019 22:56:14 -0500 Subject: [PATCH 2/4] Initial pass at implementing pow(Self,Int) with desired semantics. We don't need to do anything special for Float80, but Float and Double require some handholding because Int values are not always representable. --- Sources/Real/ScalarConformances.swift | 48 +++++++++++++++++++-------- 1 file changed, 34 insertions(+), 14 deletions(-) diff --git a/Sources/Real/ScalarConformances.swift b/Sources/Real/ScalarConformances.swift index 62fc21cd..1e85cfcf 100644 --- a/Sources/Real/ScalarConformances.swift +++ b/Sources/Real/ScalarConformances.swift @@ -50,11 +50,21 @@ extension Float: Real { } @_transparent public static func pow(_ x: Float, _ n: Int) -> Float { - // TODO: this implementation is not quite correct, because n may be - // rounded in conversion to Float. This only effects very extreme cases, - // so we'll leave it alone for now; however, it gets the sign wrong if - // it rounds an odd number to an even number, so we should fix it soon. - return libm_powf(x, Float(n)) + // If n is exactly representable as Float, we can just call powf: + if let y = Float(exactly: n) { + return libm_powf(x, y) + } + // Otherwise, n is too large to losslessly represent as Float. + // The range of "interesting" n is -1488522191 ... 1744361944; outside + // of this range, all x != 1 overflow or underflow, so only the parity + // of x matters. We don't really care about the specific range at all, + // only that the bounds fit exactly into two Floats. Mask the low 24 + // bits of n, get pow with that exponent (this contains the parity), + // then get pow with the rest (this may round, but if it does we've + // saturated anyway, so it doesn't matter). + let low = n & 0xffffff + let high = n - low + return libm_powf(x, Float(low)) * libm_powf(x, Float(high)) } @_transparent public static func root(_ x: Float, _ n: Int) -> Float { @@ -129,11 +139,22 @@ extension Double: Real { } @_transparent public static func pow(_ x: Double, _ n: Int) -> Double { - // TODO: this implementation is not quite correct, because n may be - // rounded in conversion to Double. This only effects very extreme cases, - // so we'll leave it alone for now; however, it gets the sign wrong if - // it rounds an odd number to an even number, so we should fix it soon. - return libm_pow(x, Double(n)) + // If n is exactly representable as Double, we can just call pow: + // Note that all calls on a 32b platform go down this path. + if let y = Double(exactly: n) { + return libm_pow(x, y) + } + // Otherwise, n is too large to losslessly represent as Double, so we + // just split it into two parts, high and low. This is always exact, + // so the only source of error is pow itself and the multiplication. + // + // mask constant is spelled in this funny way because if we just anded + // with the hex value, we'd get a compile error on 32b platforms, even + // though this whole branch is dead code on 32b. + let mask = Int(truncatingIfNeeded: 0x1fffffffffffff as UInt64) + let low = n & mask + let high = n - low + return libm_pow(x, Double(low)) * libm_pow(x, Double(high)) } @_transparent public static func root(_ x: Double, _ n: Int) -> Double { @@ -191,10 +212,9 @@ extension Float80: Real { } @_transparent public static func pow(_ x: Float80, _ n: Int) -> Float80 { - // TODO: this implementation is not quite correct, because n may be - // rounded in conversion to Float80. This only effects very extreme cases, - // so we'll leave it alone for now; however, it gets the sign wrong if - // it rounds an odd number to an even number, so we should fix it soon. + // Every Int value is exactly representable as Float80, so we don't need + // to do anything fancy--unlike Float and Double, we can just call the + // libm pow function. return libm_powl(x, Float80(n)) } From f7842a76dbdb8aebf54844726aaa22f98ba3d783 Mon Sep 17 00:00:00 2001 From: Stephen Canon Date: Fri, 29 Nov 2019 23:00:45 -0500 Subject: [PATCH 3/4] Capitalization search/replace-o --- Tests/RealTests/RealTests.swift | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Tests/RealTests/RealTests.swift b/Tests/RealTests/RealTests.swift index 917d031c..7ce30b69 100644 --- a/Tests/RealTests/RealTests.swift +++ b/Tests/RealTests/RealTests.swift @@ -27,10 +27,10 @@ func sanityCheck(_ expected: TestLiteralType, _ actual: T, // right through zero to the other side. XCTAssert(actual.sign == expected.sign, "\(actual) != \(expected) as \(T.self)", file: file, line: line) // Default tolerance is 16 ulps; It's OK to relax this as needed for new - // platforms, as these Checks are *not* intended to validate the math + // platforms, as these checks are *not* intended to validate the math // library--they are only intended to check that the Swift bindings are // calling the right functions in the math library. It's important, however - // not to relax the tolerance beyond a few hundred ulps, because these Checks + // not to relax the tolerance beyond a few hundred ulps, because these checks // need to detect errors where the *wrong function* is being called; e.g. // we need to flag an implentation that inadvertently called the C hypotf // function instead of hypot. This is especially important because the C From de5399615c7ba8cec63b82617db4c9e30c0dc41a Mon Sep 17 00:00:00 2001 From: Stephen Canon Date: Fri, 29 Nov 2019 23:11:15 -0500 Subject: [PATCH 4/4] Additional test coverage of saturating values. Also skip out of tests that won't compile and wouldn't make sense anyway on 32b systems. --- Tests/RealTests/RealTests.swift | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/Tests/RealTests/RealTests.swift b/Tests/RealTests/RealTests.swift index 7ce30b69..7ff931db 100644 --- a/Tests/RealTests/RealTests.swift +++ b/Tests/RealTests/RealTests.swift @@ -139,12 +139,26 @@ extension Float { sanityCheck(-7.0064924138100205091278464932003585e-46, Float.pow(-d, 1744361943)) sanityCheck( 7.0064919961905290625123586120258840e-46, Float.pow(-d, 1744361944)) sanityCheck(-7.0064915785710625079583096856510544e-46, Float.pow(-d, 1744361945)) + // Just hammer max/min exponents, these always saturate, but this will reveal + // errors in some implementations that one could try. + sanityCheck( .infinity, Self.pow(-u, Int.max - 1)) + sanityCheck( 0.0, Self.pow(-d, Int.max - 1)) + sanityCheck( 0.0, Self.pow(-u, -Int.max + 1)) + sanityCheck( .infinity, Self.pow(-d, -Int.max + 1)) + sanityCheck(-.infinity, Self.pow(-u, Int.max)) + sanityCheck(-0.0, Self.pow(-d, Int.max)) + sanityCheck(-0.0, Self.pow(-u, -Int.max)) + sanityCheck(-.infinity, Self.pow(-d, -Int.max)) + sanityCheck( 0.0, Self.pow(-u, Int.min)) + sanityCheck( .infinity, Self.pow(-d, Int.min)) } } extension Double { static func testPown() { testPownCommon() + // Following tests only make sense (and are only necessary) on 64b platforms. +#if arch(arm64) || arch(x86_64) let u: Double = 1.nextUp let d: Double = 1.nextDown // Smallest exponent not exactly representable as Double. @@ -166,6 +180,19 @@ extension Double { sanityCheck(-2.4703282292062332640976590913373022e-324, Double.pow(-d, 6711563375777760775)) sanityCheck( 2.4703282292062329898361312467121758e-324, Double.pow(-d, 6711563375777760776)) sanityCheck(-2.4703282292062327155746034020870799e-324, Double.pow(-d, 6711563375777760777)) + // Just hammer max/min exponents, these always saturate, but this will reveal + // errors in some implementations that one could try. + sanityCheck( .infinity, Self.pow(-u, Int.max - 1)) + sanityCheck( 0.0, Self.pow(-d, Int.max - 1)) + sanityCheck( 0.0, Self.pow(-u, -Int.max + 1)) + sanityCheck( .infinity, Self.pow(-d, -Int.max + 1)) + sanityCheck(-.infinity, Self.pow(-u, Int.max)) + sanityCheck(-0.0, Self.pow(-d, Int.max)) + sanityCheck(-0.0, Self.pow(-u, -Int.max)) + sanityCheck(-.infinity, Self.pow(-d, -Int.max)) + sanityCheck( 0.0, Self.pow(-u, Int.min)) + sanityCheck( .infinity, Self.pow(-d, Int.min)) +#endif } }