diff --git a/Sources/RealModule/Angle.swift b/Sources/RealModule/Angle.swift new file mode 100644 index 00000000..b951b09a --- /dev/null +++ b/Sources/RealModule/Angle.swift @@ -0,0 +1,313 @@ +//===--- Angle.swift ------------------------------------------*- swift -*-===// +// +// This source file is part of the Swift Numerics open source project +// +// Copyright (c) 2020 Apple Inc. and the Swift Numerics project authors +// Licensed under Apache License v2.0 with Runtime Library Exception +// +// See https://swift.org/LICENSE.txt for license information +// See https://swift.org/CONTRIBUTORS.txt for the list of Swift project authors +// +//===----------------------------------------------------------------------===// + +/// A wrapper type for angle operations and functions +/// +/// All trigonometric functions expect the argument to be passed as radians (Real), but this is not enforced by the type system. +/// This type serves exactly this purpose, and can be seen as an alternative to the underlying Real implementation. +public struct Angle: Equatable { + fileprivate var degreesPart: T = 0 + fileprivate var radiansPart: T = 0 + + fileprivate init() {} + + fileprivate init(degreesPart: T, radiansPart: T) { + self.degreesPart = degreesPart + self.radiansPart = radiansPart + } + + public init(radians: T) { + radiansPart = radians + } + + public static func radians(_ value: T) -> Angle { + .init(radians: value) + } + + public init(degrees: T) { + degreesPart = degrees + } + + public static func degrees(_ value: T) -> Angle { + .init(degrees: value) + } + + public var radians: T { + radiansPart + degreesPart * .pi / 180 + } + + public var degrees: T { + radiansPart * 180 / .pi + degreesPart + } +} + +extension Angle: AdditiveArithmetic { + public static var zero: Angle { .init() } + + public static func + (lhs: Angle, rhs: Angle) -> Angle { + .init(degreesPart: lhs.degreesPart + rhs.degreesPart, + radiansPart: lhs.radiansPart + rhs.radiansPart) + } + + public static func += (lhs: inout Angle, rhs: Angle) { + lhs.degreesPart += rhs.degreesPart + lhs.radiansPart += rhs.radiansPart + } + + public static func - (lhs: Angle, rhs: Angle) -> Angle { + .init(degreesPart: lhs.degreesPart - rhs.degreesPart, + radiansPart: lhs.radiansPart - rhs.radiansPart) + } + + public static func -= (lhs: inout Angle, rhs: Angle) { + lhs.degreesPart -= rhs.degreesPart + lhs.radiansPart -= rhs.radiansPart + } +} + +extension Angle { + public static func * (lhs: Angle, rhs: T) -> Angle { + .init(degreesPart: lhs.degreesPart * rhs, + radiansPart: lhs.radiansPart * rhs) + } + + public static func *= (lhs: inout Angle, rhs: T) { + lhs.degreesPart *= rhs + lhs.radiansPart *= rhs + } + + public static func * (lhs: T, rhs: Angle) -> Angle { + return rhs * lhs + } + + public static func / (lhs: Angle, rhs: T) -> Angle { + assert(rhs != 0) + return .init(degreesPart: lhs.degreesPart / rhs, + radiansPart: lhs.radiansPart / rhs) + } + + public static func /= (lhs: inout Angle, rhs: T) { + assert(rhs != 0) + lhs.degreesPart /= rhs + lhs.radiansPart /= rhs + } +} + +private extension Angle { + var piTimesFromDegrees: T { + if degreesPart.magnitude <= 180 { + return degreesPart / 180 + } + let remainder = degreesPart.remainder(dividingBy: 180) + return remainder / 180 + } + + var piTimesFromRadians: T { + if radiansPart.magnitude <= .pi { + return radiansPart / .pi + } + let remainder = radiansPart.remainder(dividingBy: .pi) + return remainder / .pi + } +} + +extension ElementaryFunctions +where Self: Real { + /// The cos of the angle. + /// + /// The degrees and radians parts are treated separately and then combined together + /// using standard trigonometric [identities]. For each part, the corresponding remainder + /// by pi or 180° is found, and the higher precision `cos(piTimes:)` function is used + /// + /// See also: + /// - + /// `ElementaryFunctions.cos()` + /// [identities]: https://en.wikipedia.org/wiki/List_of_trigonometric_identities#Angle_sum_and_difference_identities + public static func cos(_ angle: Angle) -> Self { + let piTimesDegrees = angle.piTimesFromDegrees + let piTimesRadians = angle.piTimesFromRadians + let cosa = cos(piTimes: piTimesDegrees) + let cosb = cos(piTimes: piTimesRadians) + let sina = sin(piTimes: piTimesDegrees) + let sinb = sin(piTimes: piTimesRadians) + return cosa * cosb - sina * sinb + } +} + +extension ElementaryFunctions +where Self: Real { + /// The sine of the angle. + /// + /// + /// The degrees and radians parts are treated separately and then combined together + /// using standard trigonometric [identities]. For each part, the corresponding remainder + /// by pi or 180° is found, and the higher precision `sin(piTimes:)` function is used + /// + /// See also: + /// - + /// `ElementaryFunctions.sin()` + /// [identities]: https://en.wikipedia.org/wiki/List_of_trigonometric_identities#Angle_sum_and_difference_identities + public static func sin(_ angle: Angle) -> Self { + let piTimesDegrees = angle.piTimesFromDegrees + let piTimesRadians = angle.piTimesFromRadians + let cosa = cos(piTimes: piTimesDegrees) + let cosb = cos(piTimes: piTimesRadians) + let sina = sin(piTimes: piTimesDegrees) + let sinb = sin(piTimes: piTimesRadians) + return sina * cosb + cosa * sinb + } +} + +extension ElementaryFunctions +where Self: Real { + /// The tangent of the angle. + /// + /// The degrees and radians parts are treated separately and then combined together + /// using standard trigonometric [identities]. For each part, the corresponding remainder + /// by pi or 180° is found, and the higher precision `tan(piTimes:)` function is used + /// + /// See also: + /// - + /// `ElementaryFunctions.tan()` + /// [identities]: https://en.wikipedia.org/wiki/List_of_trigonometric_identities#Angle_sum_and_difference_identities + public static func tan(_ angle: Angle) -> Self { + let piTimesDegrees = angle.piTimesFromDegrees + let piTimesRadians = angle.piTimesFromRadians + let tana = tan(piTimes: piTimesDegrees) + let tanb = tan(piTimes: piTimesRadians) + switch (tana.isFinite, tanb) { + case (false, 0): + return tana + case (false, _): + return -1 / tanb + default: + return (tana + tanb) / (1 - tana * tanb) + } + } +} + +extension Angle { + /// See also: + /// - + /// `ElementaryFunctions.acos()` + public static func acos(_ x: T) -> Self { + .radians(T.acos(x)) + } + + /// See also: + /// - + /// `ElementaryFunctions.asin()` + public static func asin(_ x: T) -> Self { + .radians(T.asin(x)) + } + + /// See also: + /// - + /// `ElementaryFunctions.atan()` + public static func atan(_ x: T) -> Self { + .radians(T.atan(x)) + } + + /// The 2-argument atan function. + /// + ///- Precondition: `x` and `y` cannot be both 0 at the same time + /// + /// See also: + /// - + /// `RealFunctions.atan2()` + public static func atan2(y: T, x: T) -> Self { + .radians(T.atan2(y: y, x: x)) + } +} + +extension Angle { + /// Checks whether the current angle is contained within a range, defined from a start and end angle. + /// + /// The comparison is performed based on the equivalent normalized angles in [-pi, pi]. + /// + /// Examples: + /// + /// ```swift + /// let angle = Angle(degrees: 175) + /// + /// // returns true + /// angle.isInRange(start: Angle(degrees: 170), end:Angle(degrees: -170)) + /// + /// // returns false + /// angle.isInRange(start: Angle(degrees: -170), end:Angle(degrees: 170)) + /// + /// // returns true + /// angle.isInRange(start: Angle(degrees: 170), end:Angle(degrees: 180)) + /// + /// // returns false + /// angle.isInRange(start: Angle(degrees: 30), end:Angle(degrees: 60)) + /// ``` + /// + /// - Parameters: + /// + /// - start: The start of the range, within which containment is checked. + /// + /// - end: The end of the range, within which containment is checked. + public func isInRange(start: Angle, end: Angle) -> Bool { + let fullNormalized = normalize(value: degrees, limit: 180) + let normalizedStart = normalize(value: start.degrees, limit: 180) + var normalizedEnd = normalize(value: end.degrees, limit: 180) + if normalizedEnd < normalizedStart { + normalizedEnd += 360 + } + return (normalizedStart <= fullNormalized && fullNormalized <= normalizedEnd) + || (normalizedStart <= fullNormalized + 360 && fullNormalized + 360 <= normalizedEnd) + } +} + +extension Angle { + /// Checks whether the current angle is close to another angle within a given tolerance + /// + /// - Precondition: `tolerance` must positive, otherwise the return value is always false + /// + /// - Parameters: + /// + /// - other: the angle from which the distance is controlled. + /// + /// - tolerance: the tolerance around `other` for which the result will be true + /// + /// - Returns: `true` if the current angle falls within the range ```[self - tolerance, self + tolerance]```, otherwise false + public func isClose(to other: Angle, within tolerance: Angle) -> Bool { + precondition(tolerance.degrees >= 0) + return isInRange(start: other - tolerance, end: other + tolerance) + } +} + +extension Angle: Comparable { + public static func < (lhs: Angle, rhs: Angle) -> Bool { + guard lhs != rhs else { + return false + } + return lhs.radians < rhs.radians + } +} + +extension Angle { + fileprivate func normalize(value: T, limit: T) -> T { + var normalized = value + + while normalized > limit { + normalized -= 2 * limit + } + + while normalized < -limit { + normalized += 2 * limit + } + + return normalized + } +} diff --git a/Sources/RealModule/Double+Real.swift b/Sources/RealModule/Double+Real.swift index 3e060f3f..e7678b87 100644 --- a/Sources/RealModule/Double+Real.swift +++ b/Sources/RealModule/Double+Real.swift @@ -108,6 +108,21 @@ extension Double: Real { } #if os(macOS) || os(iOS) || os(tvOS) || os(watchOS) + @_transparent + public static func cos(piTimes x: Double) -> Double { + libm_cospi(x) + } + + @_transparent + public static func sin(piTimes x: Double) -> Double { + libm_sinpi(x) + } + + @_transparent + public static func tan(piTimes x: Double) -> Double { + libm_tanpi(x) + } + @_transparent public static func exp10(_ x: Double) -> Double { libm_exp10(x) diff --git a/Sources/RealModule/Float+Real.swift b/Sources/RealModule/Float+Real.swift index 4147e8b8..e24282fd 100644 --- a/Sources/RealModule/Float+Real.swift +++ b/Sources/RealModule/Float+Real.swift @@ -108,6 +108,21 @@ extension Float: Real { } #if os(macOS) || os(iOS) || os(tvOS) || os(watchOS) + @_transparent + public static func cos(piTimes x: Float) -> Float { + libm_cospif(x) + } + + @_transparent + public static func sin(piTimes x: Float) -> Float { + libm_sinpif(x) + } + + @_transparent + public static func tan(piTimes x: Float) -> Float { + libm_tanpif(x) + } + @_transparent public static func exp10(_ x: Float) -> Float { libm_exp10f(x) diff --git a/Sources/RealModule/Real.swift b/Sources/RealModule/Real.swift index 0ef492da..9c8ce969 100644 --- a/Sources/RealModule/Real.swift +++ b/Sources/RealModule/Real.swift @@ -35,8 +35,92 @@ public protocol Real: FloatingPoint, RealFunctions, AlgebraicField { // it does allow us to default the implementation of a few operations, // and also provides `signGamma`. extension Real { - // Most math libraries do not provide exp10, so we need a default - // implementation. + + @_transparent + public static func cos(piTimes x: Self) -> Self { + // Cosine is even, so all we need is the magnitude. + let x = x.magnitude + // If x is not finite, the result is nan. + guard x.isFinite else { return .nan } + // If x is finite and at least .radix / .ulpOfOne, it is an even + // integer, which means that cos(piTimes: x) is 1.0 + if x >= Self(Self.radix) / .ulpOfOne { return 1 } + // Break x up as x = n/2 + f where n is an integer. In binary, the + // following computation is always exact, and trivially gives the + // correct result. + // TODO: analyze and fixup for decimal types + let n = (2*x).rounded(.toNearestOrEven) + let f = x.addingProduct(-1/2, n) + // Because cosine is 2π-periodic, we don't actually care about + // most of n; we only need the two least significant bits of n + // represented as an integer: + let quadrant = n._lowWord & 0x3 + switch quadrant { + case 0: return cos(.pi * f) + case 1: return -sin(.pi * f) + case 2: return -cos(.pi * f) + case 3: return sin(.pi * f) + default: fatalError() + } + } + + @_transparent + public static func sin(piTimes x: Self) -> Self { + // If x is not finite, the result is nan. + guard x.isFinite else { return .nan } + // If x.magnitude is finite and at least 1 / .ulpOfOne, it is an + // integer, which means that sin(piTimes: x) is ±0.0 + if x.magnitude >= 1 / .ulpOfOne { + return Self(signOf: x, magnitudeOf: 0) + } + // Break x up as x = n/2 + f where n is an integer. In binary, the + // following computation is always exact, and trivially gives the + // correct result. + // TODO: analyze and fixup for decimal types + let n = (2*x).rounded(.toNearestOrEven) + let f = x.addingProduct(-1/2, n) + // Because sine is 2π-periodic, we don't actually care about + // most of n; we only need the two least significant bits of n + // represented as an integer: + let quadrant = n._lowWord & 0x3 + switch quadrant { + case 0: return sin(.pi * f) + case 1: return cos(.pi * f) + case 2: return -sin(.pi * f) + case 3: return -cos(.pi * f) + default: fatalError() + } + } + + @_transparent + public static func tan(piTimes x: Self) -> Self { + // If x is not finite, the result is nan. + guard x.isFinite else { return .nan } + // TODO: choose policy for exact 0, 1, infinity cases and implement as + // appropriate. + // If x.magnitude is finite and at least .radix / .ulpOfOne, it is an + // even integer, which means that sin(piTimes: x) is ±0.0 and + // cos(piTimes: x) is 1.0. + if x.magnitude >= Self(Self.radix) / .ulpOfOne { + return Self(signOf: x, magnitudeOf: 0) + } + // Break x up as x = n/2 + f where n is an integer. In binary, the + // following computation is always exact, and trivially gives the + // correct result. + // TODO: analyze and fixup for decimal types + let n = (2*x).rounded(.toNearestOrEven) + let f = x.addingProduct(-1/2, n) + // Because tangent is π-periodic, we don't actually care about + // most of n; we only need the least significant bit of n represented + // as an integer: + let sector = n._lowWord & 0x1 + switch sector { + case 0: return tan(.pi * f) + case 1: return -1/tan(.pi * f) + default: fatalError() + } + } + @_transparent public static func exp10(_ x: Self) -> Self { return pow(10, x) @@ -153,3 +237,24 @@ extension Real { return nil } } + +// MARK: Implementation details +extension Real where Self: BinaryFloatingPoint { + @_transparent + public var _lowWord: UInt { + // If magnitude is small enough, we can simply convert to Int64 and then + // wrap to UInt. + if magnitude < 0x1.0p63 { + return UInt(truncatingIfNeeded: Int64(self.rounded(.down))) + } + precondition(isFinite) + // Clear any bits above bit 63; the result of this expression is + // strictly in the range [0, 0x1p64). (Note that if we had not eliminated + // small magnitudes already, the range would include tiny negative values + // which would then produce the wrong result; the branch above is not + // only for performance. + let cleared = self - 0x1p64*(self * 0x1p-64).rounded(.down) + // Now we can unconditionally convert to UInt64, and then wrap to UInt. + return UInt(truncatingIfNeeded: UInt64(cleared)) + } +} diff --git a/Sources/RealModule/RealFunctions.swift b/Sources/RealModule/RealFunctions.swift index 4701f1d2..7344f6be 100644 --- a/Sources/RealModule/RealFunctions.swift +++ b/Sources/RealModule/RealFunctions.swift @@ -2,7 +2,7 @@ // // This source file is part of the Swift Numerics open source project // -// Copyright (c) 2019 Apple Inc. and the Swift Numerics project authors +// Copyright (c) 2019-2020 Apple Inc. and the Swift Numerics project authors // Licensed under Apache License v2.0 with Runtime Library Exception // // See https://swift.org/LICENSE.txt for license information @@ -10,93 +10,212 @@ //===----------------------------------------------------------------------===// public protocol RealFunctions: ElementaryFunctions { - /// `atan(y/x)`, with sign selected according to the quadrant of `(x, y)`. + + /// `atan(y/x)`, with representative selected using the quadrant`(x,y)`. + /// + /// The [atan2 function][wiki] computes the angle (in radians, in the + /// range [-π, π]) formed between the positive real axis and the point + /// `(x,y)`. The sign of the result always matches the sign of y. + /// + /// - Warning: + /// Note the parameter ordering of this function; the `y` parameter + /// comes *before* the `x` parameter. This is a historical curiosity + /// going back to early FORTRAN math libraries. In order to minimize + /// opportunities for confusion and subtle bugs, we require explicit + /// parameter labels with this function. /// /// See also: /// - - /// - `atan()` + /// - `ElementaryFunctions.atan(_:)` + /// + /// [wiki]: https://en.wikipedia.org/wiki/Atan2 static func atan2(y: Self, x: Self) -> Self - /// The error function evaluated at `x`. + /// `cos(πx)` + /// + /// Computes the cosine of π times `x`. + /// + /// Because π is not representable in any `FloatingPoint` type, for large + /// `x`, `.cos(.pi * x)` can have arbitrarily large relative error; + /// `.cos(piTimes: x)` always provides a result with small relative error. + /// + /// This is observable even for modest arguments; consider `0.5`: + /// ```swift + /// Float.cos(.pi * 0.5) // 7.54979e-08 + /// Float.cos(piTimes: 0.5) // 0.0 + /// ``` + /// It's important to be clear that there is no bug in the example + /// given above. Every step of both computations is producing the most + /// accurate possible result. + /// + /// Symmetries: + /// - + /// - `.cos(piTimes: -x) = .cos(piTimes: x)`. + /// + /// See also: + /// - + /// - `sin(piTimes:)` + /// - `tan(piTimes:)` + /// - `ElementaryFunctions.cos(_:)` + static func cos(piTimes x: Self) -> Self + + /// `sin(πx)` + /// + /// Computes the sine of π times `x`. + /// + /// Because π is not representable in any `FloatingPoint` type, for large + /// `x`, `.sin(.pi * x)` can have arbitrarily large relative error; + /// `.sin(piTimes: x)` always provides a result with small relative error. + /// + /// This is observable even for modest arguments; consider `10`: + /// ```swift + /// Float.sin(.pi * 10) // -2.4636322e-06 + /// Float.sin(piTimes: 10) // 0.0 + /// ``` + /// It's important to be clear that there is no bug in the example + /// given above. Every step of both computations is producing the most + /// accurate possible result. + /// + /// Symmetry: + /// - + /// `.sin(piTimes: -x) = -.sin(piTimes: x)`. /// /// See also: /// - - /// - `erfc()` + /// - `cos(piTimes:)` + /// - `tan(piTimes:)` + /// - `ElementaryFunctions.sin(_:)` + static func sin(piTimes x: Self) -> Self + + /// `tan(πx)` + /// + /// Computes the tangent of π times `x`. + /// + /// Because π is not representable in any `FloatingPoint` type, for large + /// `x`, `.tan(.pi * x)` can have arbitrarily large relative error; + /// `.tan(piTimes: x)` always provides a result with small relative error. + /// + /// This is observable even for modest arguments; consider `0.5`: + /// ```swift + /// Float.tan(.pi * 0.5) // 13245402.0 + /// Float.tan(piTimes: 0.5) // infinity + /// ``` + /// It's important to be clear that there is no bug in the example + /// given above. Every step of both computations is producing the most + /// accurate possible result. + /// + /// Symmetry: + /// - + /// `.tan(piTimes: -x) = -.tan(piTimes: x)`. + /// + /// See also: + /// - + /// - `cos(piTimes:)` + /// - `sin(piTimes:)` + /// - `ElementaryFunctions.tan(_:)` + static func tan(piTimes x: Self) -> Self + + /// The [error function][wiki] evaluated at `x`. + /// + /// See also: + /// - + /// - `erfc(_:)` + /// + /// [wiki]: https://en.wikipedia.org/wiki/Error_function static func erf(_ x: Self) -> Self - /// The complimentary error function evaluated at `x`. + /// The complimentary [error function][wiki] evaluated at `x`. /// /// See also: /// - - /// - `erf()` + /// - `erf(_:)` + /// + /// [wiki]: https://en.wikipedia.org/wiki/Error_function static func erfc(_ x: Self) -> Self - /// 2^x + /// 2ˣ /// /// See also: /// - - /// - `exp()` - /// - `expMinusOne()` - /// - `exp10()` - /// - `log2()` - /// - `pow()` + /// - `ElementaryFunctions.exp(_:)` + /// - `ElementaryFunctions.expMinusOne(_:)` + /// - `exp10(_:)` + /// - `log2(_:)` + /// - `ElementaryFunctions.pow(_:)` static func exp2(_ x: Self) -> Self - /// 10^x + /// 10ˣ /// /// See also: /// - - /// - `exp()` - /// - `expMinusOne()` - /// - `exp2()` - /// - `log10()` - /// - `pow()` + /// - `ElementaryFunctions.exp(_:)` + /// - `ElementaryFunctions.expMinusOne(_:)` + /// - `exp2(_:)` + /// - `log10(_:)` + /// - `ElementaryFunctions.pow(_:)` static func exp10(_ x: Self) -> Self - /// `sqrt(x*x + y*y)`, computed in a manner that avoids spurious overflow or - /// underflow. +/// The square root of the sum of squares of `x` and `y`. +/// +/// The naive expression `.sqrt(x*x + y*y)` and overflow +/// or underflow if `x` or `y` is not well-scaled, producing zero or +/// infinity, even when the mathematical result is representable. +/// +/// The [hypot][wiki] takes care to avoid this, and always +/// produces an accurate result when one is available. +/// +/// See also: +/// - +/// - `ElementaryFunctions.sqrt(_:)` +/// +/// [wiki]: https://en.wikipedia.org/wiki/Hypot static func hypot(_ x: Self, _ y: Self) -> Self - /// The gamma function Γ(x). + /// The [gamma function][wiki] Γ(x). /// /// See also: /// - - /// - `logGamma()` - /// - `signGamma()` + /// - `logGamma(_:)` + /// - `signGamma(_:)` + /// + /// [wiki]: https://en.wikipedia.org/wiki/Gamma_function static func gamma(_ x: Self) -> Self /// The base-2 logarithm of `x`. /// /// See also: /// - - /// - `exp2()` - /// - `log()` - /// - `log(onePlus:)` - /// - `log10()` + /// - `exp2(_:)` + /// - `ElementaryFunctions.log(_:)` + /// - `ElementaryFunctions.log(onePlus:)` + /// - `log10(_:)` static func log2(_ x: Self) -> Self /// The base-10 logarithm of `x`. /// /// See also: /// - - /// - `exp10()` - /// - `log()` - /// - `log(onePlus:)` - /// - `log2()` + /// - `exp10(_:)` + /// - `ElementaryFunctions.log(_:)` + /// - `ElementaryFunctions.log(onePlus:)` + /// - `log2(_:)` static func log10(_ x: Self) -> Self #if !os(Windows) - /// The logarithm of the absolute value of the gamma function, log(|Γ(x)|). + /// The logarithm of the absolute value of the [gamma function][wiki], log(|Γ(x)|). /// - /// Not available on Windows targets. + /// - Warning: + /// Not available on Windows. /// /// See also: /// - - /// - `gamma()` - /// - `signGamma()` + /// - `gamma(_:)` + /// - `signGamma(_:)` + /// + /// [wiki]: https://en.wikipedia.org/wiki/Gamma_function static func logGamma(_ x: Self) -> Self - /// The sign of the gamma function, Γ(x). + /// The sign of the [gamma function][wiki], Γ(x). /// /// For `x >= 0`, `signGamma(x)` is `.plus`. For negative `x`, `signGamma(x)` /// is `.plus` when `x` is an integer, and otherwise it is `.minus` whenever @@ -105,12 +224,15 @@ public protocol RealFunctions: ElementaryFunctions { /// This function is used together with `logGamma`, which computes the /// logarithm of the absolute value of Γ(x), to recover the sign information. /// - /// Not available on Windows targets. + /// - Warning: + /// Not available on Windows. /// /// See also: /// - - /// - `gamma()` - /// - `logGamma()` + /// - `gamma(_:)` + /// - `logGamma(_:)` + /// + /// [wiki]: https://en.wikipedia.org/wiki/Gamma_function static func signGamma(_ x: Self) -> FloatingPointSign #endif @@ -118,4 +240,10 @@ public protocol RealFunctions: ElementaryFunctions { /// /// Whichever is faster should be chosen by the compiler statically. static func _mulAdd(_ a: Self, _ b: Self, _ c: Self) -> Self + + // MARK: Implementation details + + /// The low-word of the integer formed by truncating this value. + var _lowWord: UInt { get } } + diff --git a/Sources/_NumericsShims/include/_NumericsShims.h b/Sources/_NumericsShims/include/_NumericsShims.h index 7e25eab6..0927d025 100644 --- a/Sources/_NumericsShims/include/_NumericsShims.h +++ b/Sources/_NumericsShims/include/_NumericsShims.h @@ -116,6 +116,21 @@ HEADER_SHIM float libm_exp2f(float x) { } #if __APPLE__ +HEADER_SHIM float libm_cospif(float x) { + extern float __cospif(float); + return __cospif(x); +} + +HEADER_SHIM float libm_sinpif(float x) { + extern float __sinpif(float); + return __sinpif(x); +} + +HEADER_SHIM float libm_tanpif(float x) { + extern float __tanpif(float); + return __tanpif(x); +} + HEADER_SHIM float libm_exp10f(float x) { extern float __exp10f(float); return __exp10f(x); @@ -241,6 +256,21 @@ HEADER_SHIM double libm_exp2(double x) { } #if __APPLE__ +HEADER_SHIM double libm_cospi(double x) { + extern double __cospi(double); + return __cospi(x); +} + +HEADER_SHIM double libm_sinpi(double x) { + extern double __sinpi(double); + return __sinpi(x); +} + +HEADER_SHIM double libm_tanpi(double x) { + extern double __tanpi(double); + return __tanpi(x); +} + HEADER_SHIM double libm_exp10(double x) { extern double __exp10(double); return __exp10(x); diff --git a/Tests/RealTests/AngleTests.swift b/Tests/RealTests/AngleTests.swift new file mode 100644 index 00000000..d85f1e48 --- /dev/null +++ b/Tests/RealTests/AngleTests.swift @@ -0,0 +1,244 @@ +//===--- AngleTests.swift -------------------------------------*- swift -*-===// +// +// This source file is part of the Swift Numerics open source project +// +// Copyright (c) 2020 Apple Inc. and the Swift Numerics project authors +// Licensed under Apache License v2.0 with Runtime Library Exception +// +// See https://swift.org/LICENSE.txt for license information +// See https://swift.org/CONTRIBUTORS.txt for the list of Swift project authors +// +//===----------------------------------------------------------------------===// + +import RealModule +import XCTest +import _TestSupport + +internal extension Real +where Self: BinaryFloatingPoint { + static func conversionBetweenRadiansAndDegreesChecks() { + let angleFromRadians = Angle(radians: .pi / 3) + assertClose(60, angleFromRadians.degrees) + + let angleFromDegrees = Angle(degrees: 120) + assertClose(2 * .pi / 3, angleFromDegrees.radians) + } + + static func inverseTrigonometricFunctionChecks() { + assertClose(0, Angle.acos(1).degrees) + assertClose(30, Angle.acos(sqrt(3)/2).degrees) + assertClose(45, Angle.acos(sqrt(2)/2).degrees) + assertClose(60, Angle.acos(0.5).degrees) + assertClose(90, Angle.acos(0).degrees) + assertClose(120, Angle.acos(-0.5).degrees) + assertClose(135, Angle.acos(-sqrt(2)/2).degrees) + assertClose(150, Angle.acos(-sqrt(3)/2).degrees) + assertClose(180, Angle.acos(-1).degrees) + + assertClose(-90, Angle.asin(-1).degrees) + assertClose(-60, Angle.asin(-sqrt(3)/2).degrees) + assertClose(-45, Angle.asin(-sqrt(2)/2).degrees) + assertClose(-30, Angle.asin(-0.5).degrees) + assertClose(0, Angle.asin(0).degrees) + assertClose(30, Angle.asin(0.5).degrees) + assertClose(45, Angle.asin(sqrt(2)/2).degrees) + assertClose(60, Angle.asin(sqrt(3)/2).degrees) + assertClose(90, Angle.asin(1).degrees) + + assertClose(-90, Angle.atan(-.infinity).degrees) + assertClose(-60, Angle.atan(-sqrt(3)).degrees) + assertClose(-45, Angle.atan(-1).degrees) + assertClose(-30, Angle.atan(-sqrt(3)/3).degrees) + assertClose(0, Angle.atan(0).degrees) + assertClose(30, Angle.atan(sqrt(3)/3).degrees) + assertClose(45, Angle.atan(1).degrees) + assertClose(60, Angle.atan(sqrt(3)).degrees) + assertClose(90, Angle.atan(.infinity).degrees) + + assertClose(-150, Angle.atan2(y:-sqrt(3), x:-3).degrees) + assertClose(-135, Angle.atan2(y:-1, x:-1).degrees) + assertClose(-120, Angle.atan2(y:-sqrt(3), x:-1).degrees) + assertClose(-90, Angle.atan2(y:-1, x:0).degrees) + assertClose(-60, Angle.atan2(y:-sqrt(3), x:1).degrees) + assertClose(-45, Angle.atan2(y:-1, x:1).degrees) + assertClose(-30, Angle.atan2(y:-sqrt(3), x:3).degrees) + assertClose(0, Angle.atan2(y:0, x:1).degrees) + assertClose(30, Angle.atan2(y:sqrt(3), x:3).degrees) + assertClose(45, Angle.atan2(y:1, x:1).degrees) + assertClose(60, Angle.atan2(y:sqrt(3), x:1).degrees) + assertClose(90, Angle.atan2(y:1, x:0).degrees) + assertClose(120, Angle.atan2(y:sqrt(3), x:-1).degrees) + assertClose(135, Angle.atan2(y:1, x:-1).degrees) + assertClose(150, Angle.atan2(y:sqrt(3), x:-3).degrees) + assertClose(180, Angle.atan2(y:0, x:-1).degrees) + + assertClose(1.1863995522992575361931268186727044683, Angle.acos(0.375).radians) + assertClose(0.3843967744956390830381948729670469737, Angle.asin(0.375).radians) + assertClose(0.3587706702705722203959200639264604997, Angle.atan(0.375).radians) + assertClose(0.54041950027058415544357836460859991, Angle.atan2(y: 0.375, x: 0.625).radians) + } + + static func trigonometricFunctionChecks() { + XCTAssertEqual(1, cos(Angle(degrees: -360))) + XCTAssertEqual(0, cos(Angle(degrees: -270))) + XCTAssertEqual(-1, cos(Angle(degrees: -180))) + assertClose(-0.86602540378443864676372317075293618347, cos(Angle(degrees: -150))) + assertClose(-0.70710678118654752440084436210484903929, cos(Angle(degrees: -135))) + assertClose(-0.5, cos(Angle(degrees: -120))) + XCTAssertEqual(0, cos(Angle(degrees: -90))) + assertClose(0.5, cos(Angle(degrees: -60))) + assertClose(0.70710678118654752440084436210484903929, cos(Angle(degrees: -45))) + assertClose(0.86602540378443864676372317075293618347, cos(Angle(degrees: -30))) + XCTAssertEqual(1, cos(Angle(degrees: 0))) + assertClose(0.86602540378443864676372317075293618347, cos(Angle(degrees: 30))) + assertClose(0.70710678118654752440084436210484903929, cos(Angle(degrees: 45))) + assertClose(0.5, cos(Angle(degrees: 60))) + XCTAssertEqual(0, cos(Angle(degrees: 90))) + assertClose(-0.5, cos(Angle(degrees: 120))) + assertClose(-0.70710678118654752440084436210484903929, cos(Angle(degrees: 135))) + assertClose(-0.86602540378443864676372317075293618347, cos(Angle(degrees: 150))) + XCTAssertEqual(-1, cos(Angle(degrees: 180))) + XCTAssertEqual(0, cos(Angle(degrees: 270))) + XCTAssertEqual(1, cos(Angle(degrees: 360))) + + XCTAssertEqual(0, sin(Angle(degrees: -360))) + XCTAssertEqual(1, sin(Angle(degrees: -270))) + XCTAssertEqual(0, sin(Angle(degrees: -180))) + assertClose(-0.5, sin(Angle(degrees: -150))) + assertClose(-0.70710678118654752440084436210484903929, sin(Angle(degrees: -135))) + assertClose(-0.86602540378443864676372317075293618347, sin(Angle(degrees: -120))) + XCTAssertEqual(-1, sin(Angle(degrees: -90))) + assertClose(-0.86602540378443864676372317075293618347, sin(Angle(degrees: -60))) + assertClose(-0.70710678118654752440084436210484903929, sin(Angle(degrees: -45))) + assertClose(-0.5, sin(Angle(degrees: -30))) + XCTAssertEqual(0, sin(Angle(degrees: 0))) + assertClose(0.5, sin(Angle(degrees: 30))) + assertClose(0.70710678118654752440084436210484903929, sin(Angle(degrees: 45))) + assertClose(0.86602540378443864676372317075293618347, sin(Angle(degrees: 60))) + XCTAssertEqual(1, sin(Angle(degrees: 90))) + assertClose(0.86602540378443864676372317075293618347, sin(Angle(degrees: 120))) + assertClose(0.70710678118654752440084436210484903929, sin(Angle(degrees: 135))) + assertClose(0.5, sin(Angle(degrees: 150))) + XCTAssertEqual(0, sin(Angle(degrees: 180))) + XCTAssertEqual(-1, sin(Angle(degrees: 270))) + XCTAssertEqual(0, sin(Angle(degrees: 360))) + + XCTAssertEqual(0, tan(Angle(degrees: -360))) + XCTAssertEqual(.infinity, tan(Angle(degrees: -270))) + XCTAssertEqual(0, tan(Angle(degrees: -180))) + assertClose(0.57735026918962576450914878050195745565, tan(Angle(degrees: -150))) + XCTAssertEqual(1, tan(Angle(degrees: -135))) + assertClose(1.7320508075688772935274463415058723669, tan(Angle(degrees: -120))) + XCTAssertEqual(-.infinity, tan(Angle(degrees: -90))) + assertClose(-1.7320508075688772935274463415058723669, tan(Angle(degrees: -60))) + XCTAssertEqual(-1, tan(Angle(degrees: -45))) + assertClose(-0.57735026918962576450914878050195745565, tan(Angle(degrees: -30))) + XCTAssertEqual(0, tan(Angle(degrees: 0))) + assertClose(0.57735026918962576450914878050195745565, tan(Angle(degrees: 30))) + XCTAssertEqual(1, tan(Angle(degrees: 45))) + assertClose(1.7320508075688772935274463415058723669, tan(Angle(degrees: 60))) + XCTAssertEqual(.infinity, tan(Angle(degrees: 90))) + assertClose(-1.7320508075688772935274463415058723669, tan(Angle(degrees: 120))) + XCTAssertEqual(-1, tan(Angle(degrees: 135))) + assertClose(-0.57735026918962576450914878050195745565, tan(Angle(degrees: 150))) + XCTAssertEqual(0, tan(Angle(degrees: 180))) + XCTAssertEqual(-.infinity, tan(Angle(degrees: 270))) + XCTAssertEqual(0, tan(Angle(degrees: 360))) + + assertClose(0.9305076219123142911494767922295555080, cos(Angle(radians: 0.375))) + assertClose(0.3662725290860475613729093517162641571, sin(Angle(radians: 0.375))) + assertClose(0.3936265759256327582294137871012180981, tan(Angle(radians: 0.375))) + } + + static func additiveArithmeticTests() { + let angle1 = Angle(degrees: 90) + let angle2 = Angle(radians: .pi) + let sum = angle1 + angle2 + XCTAssertEqual(270, sum.degrees) + XCTAssertEqual(3 * .pi / 2, sum.radians) + XCTAssertEqual(360, (sum + angle1).degrees) + XCTAssertEqual(2 * .pi, (sum + angle1).radians) + var angle = Angle(degrees: 30) + assertClose(50, (angle + Angle(degrees: 20)).degrees) + assertClose(10, (angle - Angle(degrees: 20)).degrees) + XCTAssertEqual(Angle(degrees: 60), angle * 2) + XCTAssertEqual(Angle(degrees: 60), 2 * angle) + XCTAssertEqual(Angle(degrees: 15), angle / 2) + angle += Angle(degrees: 10) + XCTAssertEqual(Angle(degrees: 40), angle) + angle -= Angle(degrees: 20) + XCTAssertEqual(Angle(degrees: 20), angle) + angle *= 3 + XCTAssertEqual(Angle(degrees: 60), angle) + angle /= 6 + XCTAssertEqual(Angle(degrees: 10), angle) + } + + static func rangeContainmentTests() { + let angle175Deg = Angle(degrees: -5) + Angle.radians(Self.pi) + let angle170Deg = Angle(degrees: 350) + Angle.radians(-Self.pi) + let angleMinus170Deg = Angle(degrees: -350) + Angle.radians(Self.pi) + + XCTAssertTrue(angle175Deg.isInRange(start: angle170Deg, end: angleMinus170Deg)) + XCTAssertTrue(Angle(degrees: -175).isInRange(start: angle170Deg, end: angleMinus170Deg)) + + XCTAssertFalse(angle175Deg.isInRange(start: angleMinus170Deg, end: angle170Deg)) + XCTAssertFalse(Angle(degrees: -175).isInRange(start: angleMinus170Deg, end: angle170Deg)) + + XCTAssertTrue(Angle(degrees: 10).isInRange(start: Angle(degrees: 10), end: Angle(degrees: 60))) + XCTAssertTrue(Angle(degrees: 60).isInRange(start: Angle(degrees: 10), end: Angle(degrees: 60))) + XCTAssertTrue(Angle(degrees: 30).isInRange(start: Angle(degrees: 10), end: Angle(degrees: 60))) + + XCTAssertFalse(Angle(degrees: 0).isInRange(start: Angle(degrees: 10), end: Angle(degrees: 60))) + XCTAssertFalse(Angle(degrees: 70).isInRange(start: Angle(degrees: 10), end: Angle(degrees: 60))) + } + + static func distanceChecks() { + let angle175Deg = Angle(degrees: -5) + Angle.radians(Self.pi) + let angle170Deg = Angle(degrees: 350) + Angle.radians(-Self.pi) + let angleMinus170Deg = Angle(degrees: -350) + Angle.radians(Self.pi) + + XCTAssertFalse(angle170Deg.isClose(to: angle175Deg, within: Angle(degrees: 2))) + XCTAssertFalse(angle175Deg.isClose(to: angle170Deg, within: Angle(degrees: 2))) + + XCTAssertTrue(angle170Deg.isClose(to: angle175Deg, within: Angle(degrees: 10))) + XCTAssertTrue(angle175Deg.isClose(to: angle170Deg, within: Angle(degrees: 5))) + + XCTAssertTrue(angleMinus170Deg.isClose(to: angle175Deg, within: Angle(degrees: 20))) + XCTAssertFalse(angleMinus170Deg.isClose(to: angle175Deg, within: Angle(degrees: 10))) + } +} + +final class AngleTests: XCTestCase { + private func execute(`for` Type: T.Type) { + Type.conversionBetweenRadiansAndDegreesChecks() + Type.inverseTrigonometricFunctionChecks() + Type.trigonometricFunctionChecks() + Type.additiveArithmeticTests() + Type.rangeContainmentTests() + Type.distanceChecks() + } + + + #if swift(>=5.3) && !(os(macOS) || os(iOS) && targetEnvironment(macCatalyst)) + func testFloat16() { + if #available(iOS 14.0, watchOS 14.0, tvOS 7.0, *) { + execute(for: Float16.self) + } + } + #endif + + func testFloat() { + execute(for: Float.self) + } + + func testDouble() { + execute(for: Double.self) + } + + #if (arch(i386) || arch(x86_64)) && !os(Windows) && !os(Android) + func testFloat80() { + execute(for: Float80.self) + } + #endif +} diff --git a/Tests/RealTests/ElementaryFunctionChecks.swift b/Tests/RealTests/ElementaryFunctionChecks.swift index 7dfd7caf..2574c57c 100644 --- a/Tests/RealTests/ElementaryFunctionChecks.swift +++ b/Tests/RealTests/ElementaryFunctionChecks.swift @@ -1,4 +1,4 @@ -//===--- ElementaryFunctionChecks.swift ------------------------*- swift -*-===// +//===--- ElementaryFunctionChecks.swift -----------------------*- swift -*-===// // // This source file is part of the Swift Numerics open source project // @@ -27,6 +27,11 @@ internal func assertClose( file: StaticString = #file, line: UInt = #line ) -> T where T: BinaryFloatingPoint { + // we need to first check if the values are zero before we check the signs + // otherwise, "0" and "-0" compare as unequal (eg. sin(-180) == 0) + let expectedT = T(expected) + if observed.isZero && expectedT.isZero { return 0 } + // Shortcut relative-error check if we got the sign wrong; it's OK to // underflow to zero, but we do not want to allow going right through // zero and getting the sign wrong. @@ -38,8 +43,7 @@ internal func assertClose( if observed.isNaN && expected.isNaN { return 0 } // If T(expected) is zero or infinite, and matches observed, the error // is zero. - let expectedT = T(expected) - if observed.isZero && expectedT.isZero { return 0 } + if observed.isInfinite && expectedT.isInfinite { return 0 } // Special-case where only one of expectedT or observed is infinity. // Artificially knock everything down a binade, treat actual infinity as @@ -83,6 +87,11 @@ internal func assertClose( )) } +extension FloatingPoint { + var isPositiveZero: Bool { self == 0 && sign == .plus } + var isNegativeZero: Bool { self == 0 && sign == .minus } +} + internal extension ElementaryFunctions where Self: BinaryFloatingPoint { static func elementaryFunctionChecks() { assertClose(1.1863995522992575361931268186727044683, Self.acos(0.375)) @@ -153,6 +162,7 @@ final class ElementaryFunctionChecks: XCTestCase { func testFloat80() { Float80.elementaryFunctionChecks() Float80.realFunctionChecks() + Float80.cosPiTests() } #endif }