Update README.md

This PR adds a PodSpec file and a `.swift-version` file to support CocoaPods.

.codecov.yml

@@ -0,0 +1,2 @@
+ignore:
+  - "Tests/.*"

.gitignore

@@ -1,4 +1,5 @@
 .DS_Store
 /.build
+/build
 /Packages
 /*.xcodeproj

.swift-version

@@ -0,0 +1 @@
+4.0

.travis.yml

@@ -0,0 +1,20 @@
+matrix:
+  include:
+    - os: linux
+      language: generic
+      sudo: required
+      dist: trusty
+      before_install:
+        - wget -q -O - https://swift.org/keys/all-keys.asc | gpg --import -
+        - wget https://swift.org/builds/swift-4.1-release/ubuntu1404/swift-4.1-RELEASE/swift-4.1-RELEASE-ubuntu14.04.tar.gz
+        - tar xzf swift-4.1-RELEASE-ubuntu14.04.tar.gz
+        - export PATH=${PWD}/swift-4.1-RELEASE-ubuntu14.04/usr/bin:"${PATH}"
+      script:
+        - swift test -Xcc -D_GNU_SOURCE
+    - os: osx
+      language: objective-c
+      osx_image: xcode9.3beta
+      script:
+        - swift package generate-xcodeproj && xcodebuild clean test -project NumericAnnex.xcodeproj -scheme NumericAnnex-Package -enableCodeCoverage YES
+      after_success:
+        - bash <(curl -s https://codecov.io/bash) -X xcodeplist

CHANGELOG.md

@@ -0,0 +1,95 @@
+## 0.1.19: Bug fixes
+
+Added tie-breaking exponentiation operator implementations for concrete integer types in the standard library.
+
+## 0.1.18: Exponentiation operators
+
+* Added the exponentiation operator `**` and the compound assignment operator `**=`.
+* Deprecated `BinaryInteger.pow(_:_:)` and `Math.pow(_:_:)`.
+
+## 0.1.17: CocoaPods
+
+Added a PodSpec file (and a `.swift-version` file) for CocoaPods support.
+
+## 0.1.16: Bug fixes
+
+Worked around a compiler segmentation fault by assuming memory bound instead of binding memory.
+
+## 0.1.15: Bug fixes and documentation improvements
+
+* Updated doc comments.
+* Updated generic algorithms for revised integer protocols.
+
+## 0.1.14: Bug fixes
+
+* Removed `Complex.isCanonical` and `Complex.polar`.
+* Fixed implementation of `Complex.isNaN` and `Complex.isSignalingNaN` to align with C/C++ standards: a complex value is NaN if at least one of its components is NaN *and* the other is not infinite.
+* Fixed implementation of `Complex` multiplication and division to align handling of special values with C/C++ standards.
+* Fixed implementation of `Complex.tanh(_:)` to align handling of special values with C/C++ standards.
+* Refined implementation of `Complex.description` for values with a negative imaginary component.
+* Fixed implementation of `Rational` initializers that convert from a `BinaryFloatingPoint` value.
+
+## 0.1.13: Bug fixes
+
+Fixed implementation of `Rational.isProper`.
+
+## 0.1.12: Integer square root and cube root
+
+* Added integer square root and cube root functions.
+* Removed `power(of:)` requirement from `Math`, substituting `pow(_:_:)` (formerly an extension method).
+
+## 0.1.11: Bug fixes
+
+* Expanded tests.
+* Fixed implementation of certain default implementations for protocol requirements of `Math` and `Real`.
+
+## 0.1.10: Codable
+
+Conformed `Complex` and `Rational` to `Codable`.
+
+## 0.1.9: Bug fixes
+
+* Renamed `FloatingPointMath` to `Real`.
+* Updated implementation of `Complex.hyperbolicTangent` to align handling of special values to revised C standard.
+* Updated `swift-tools-version` to `4.0`.
+
+## 0.1.8: Bug fixes and documentation improvements
+
+* Removed initializer requirements from `FloatingPointMath`, moving functionality to an extension to `BinaryFloatingPoint`.
+* Reorganized code and edited comments to improve documentation.
+
+## 0.1.7: Random numbers
+
+* Added initial implementation of a pseudo-random number generator protocol and conforming types.
+* Restored constraints on `FloatingPointMath` to refine `FloatingPoint` and not `BinaryFloatingPoint`.
+
+## 0.1.6: Bug fixes
+
+Modified `FloatingPointMath` to refine `BinaryFloatingPoint` instead of `FloatingPoint`.
+
+## 0.1.5: Bug fixes
+
+Fixed implementation of comparison in `Rational<T>`.
+
+## 0.1.4: Bug fixes
+
+Fixed several implementations in `Rational<T>` and added tests.
+
+## 0.1.3: Rational conversions
+
+Added initializers to convert from floating-point or integer values to `Rational<T>`, and vice versa.
+
+## 0.1.2: Initial pre-release
+
+Initial implementation of:
+
+* `Math`
+* `FloatingPointMath`
+* `Complex`
+* `Rational`
+
+...and extensions to:
+
+* `UnsignedInteger`
+* `Float`
+* `Double`

CONTRIBUTING.md

@@ -0,0 +1,33 @@
+Thanks for taking the time to consider how you might contribute to the project!
+
+
+## Reporting Bugs
+
+First, check that the bug has not [already been reported][1]. If no existing
+issue addresses the bug, [open a new issue][2] on GitHub. Be sure to include **a
+title and a description** that clearly explain the problem with as much detail
+as possible; consider attaching **code samples** that demonstrate the unexpected
+behavior.
+
+
+## Patching Code
+
+When ready to submit a patch, [open a new pull request][3] on GitHub. Include
+**a clear description** of the problem and your solution, and reference any
+relevant open issues.
+
+Code should adhere to the prevailing conventions of the project, which are
+similar to that of the Swift standard library. Use two spaces instead of tabs
+and end each file with a single newline character.
+
+By submitting a patch, you represent that: (a) you have the right to license the
+contribution to this project; (b) you agree to license the contribution under
+the project's current [license][4]; and (c) you give permission to re-license
+the contribution under any other [Open Source Initiative (OSI)-approved
+license][5] that the project may adopt in the future.
+
+[1]: https://github.com/xwu/NumericAnnex/issues
+[2]: https://github.com/xwu/NumericAnnex/issues/new
+[3]: https://github.com/xwu/NumericAnnex/compare
+[4]: https://github.com/xwu/NumericAnnex/blob/master/LICENSE
+[5]: https://opensource.org/licenses

NumericAnnex.podspec

@@ -0,0 +1,19 @@
+Pod::Spec.new do |s|
+
+  s.name                  = "NumericAnnex"
+  s.version               = "0.1.19"
+  s.author                = "Xiaodi Wu"
+  s.social_media_url      = "https://twitter.com/xwu"
+  s.license               = "MIT"
+
+  s.homepage              = "https://github.com/xwu/NumericAnnex"
+  s.source                = { :git => "https://github.com/xwu/NumericAnnex.git", :tag => "#{s.version}" }
+  s.summary               = "A supplement to the numeric facilities provided in the Swift standard library."
+  s.documentation_url     = "https://xwu.github.io/NumericAnnex"
+
+  s.ios.deployment_target = "8.0"
+  s.framework             = "Security"
+  s.source_files          = "Sources"
+  # s.exclude_files       = "Sources/Exclude"
+
+end

Package.swift

@@ -1,7 +1,25 @@
-// swift-tools-version:3.1
+// swift-tools-version:4.0
 
 import PackageDescription
 
 let package = Package(
-  name: "NumericAnnex"
+  name: "NumericAnnex",
+  products: [
+    .library(
+      name: "NumericAnnex",
+      targets: ["NumericAnnex"]
+    ),
+  ],
+  dependencies: [],
+  targets: [
+    .target(
+      name: "NumericAnnex",
+      dependencies: [],
+      path: "Sources"
+    ),
+    .testTarget(
+      name: "NumericAnnexTests",
+      dependencies: ["NumericAnnex"]
+    ),
+  ]
 )

README.md

@@ -1,23 +1,26 @@
-# NumericAnnex
+# <img src="https://xwu.github.io/NumericAnnex/img/NumericAnnex-2017-09-02.svg" alt="NumericAnnex" height="72"><br>NumericAnnex
 
 NumericAnnex supplements the numeric facilities provided in the Swift standard
 library.
 
+[![Build Status](https://travis-ci.org/xwu/NumericAnnex.svg?branch=master)](https://travis-ci.org/xwu/NumericAnnex)
+[![codecov](https://codecov.io/gh/xwu/NumericAnnex/branch/master/graph/badge.svg)](https://codecov.io/gh/xwu/NumericAnnex)
+
 
 ## Features
 
-- [x] `BinaryInteger` exponentiation, greatest common divisor, and least common
-      multiple functions.
-- [x] `Math`, a protocol for types providing square root, cube root, and
-      elementary transcendental functions.
-- [x] `FloatingPointMath`, a protocol for floating-point types providing certain
-      special functions.
-- [x] `PRNG`, a protocol for pseudo-random number generators (PRNGs).
-- [x] `Complex`, a value type to represent complex values in Cartesian form.
-- [x] `Rational`, a value type to represent rational values, which supports
-      division by zero.
-- [x] `Random` and `Random.Xoroshiro`, two reference types implementing
-      efficient non-cryptographically-secure PRNGs.
+- The exponentiation operator `**` and the compound assignment operator `**=`.
+- Extension methods for `BinaryInteger` exponentiation, square root, cube root,
+  greatest common divisor, and least common multiple.
+- `Math`, a protocol for signed numeric types that support elementary functions.
+- `Real`, a protocol for floating-point types that support elementary functions
+  and a selection of special functions.
+- `PRNG`, a protocol for pseudo-random number generators.
+- `Rational`, a value type to represent rational values which supports division
+  by zero.
+- `Complex`, a value type to represent complex values in Cartesian form.
+- `Random` and `Random.Xoroshiro`, two reference types implementing efficient
+  pseudo-random number generators.
 
 > Note: This project is in the early stages of development and is not
 > production-ready at this time.
@@ -25,22 +28,28 @@ library.
 
 ## Requirements
 
-NumericAnnex now requires a recent development snapshot of Swift 4.0 that
-includes the revised numeric protocols. It requires either `Darwin.C` or `Glibc`
-for transcendental functions provided by the C standard library and either
-`Security` (Apple platforms) or `Glibc` (Linux) for cryptographically secure
-random bytes.
+NumericAnnex requires Swift 4.1 (`swift-4.1-branch`) or Swift 4.2 (`master`). On
+Apple platforms, it also requires the Security framework for cryptographically
+secure random bytes.
 
 
 ## Installation
 
-After NumericAnnex has been cloned or downloaded locally, build the library by
-invoking `swift build`, or run tests with `swift test`. An Xcode project can be
-generated by invoking `swift package generate-xcodeproj`.
+After NumericAnnex has been cloned or downloaded locally, build the library
+using the command `swift build` (macOS) or `swift build -Xcc -D_GNU_SOURCE`
+(Linux). Run tests with the command `swift test` (macOS) or
+`swift test -Xcc -D_GNU_SOURCE` (Linux). An Xcode project can be generated with
+the command `swift package generate-xcodeproj`.
+
+To add the package as a dependency using [CocoaPods](https://cocoapods.org),
+insert the following line in your `Podfile`:
+
+```ruby
+pod 'NumericAnnex', '~> 0.1.19'
+```
 
 [Swift Package Manager](https://swift.org/package-manager/) can also be used to
-add the package as a dependency for your own project. See Swift documentation
-for details.
+add the package as a dependency. See Swift documentation for details.
 
 
 ## Basic Usage
@@ -48,6 +57,18 @@ for details.
 ```swift
 import NumericAnnex
 
+print(2 ** 3)
+// Prints "8".
+
+print(4.0 ** 5.0)
+// Prints "1024.0".
+
+print(Int.cbrt(8))
+// Prints "2".
+
+print(Double.cbrt(27.0))
+// Prints "3.0".
+
 var x: Ratio = 1 / 4
 // Ratio is a type alias for Rational<Int>.
 
@@ -76,8 +97,9 @@ print(Complex.cos(z).real)
 
 ## Documentation
 
-All public protocols, types, and functions have been carefully documented in
-comments.
+All public protocols, types, and functions have been carefully documented in the
+code. See the [formatted reference](https://xwu.github.io/NumericAnnex/) for
+details.
 
 The project adheres to many design patterns found in the Swift standard library.
 For example, `Math` types provide methods such as `cubeRoot()` and `tangent()`
@@ -93,9 +115,8 @@ counterpart (e.g., `Complex128.cbrt(-8) != -2`).
 
 ## Future Directions
 
-- [ ] Add more tests, including performance tests
-- [ ] Design and implement additional methods on `PRNG`
-- [ ] Design and implement `BigInt`
+- Add more tests, including performance tests
+- Design and implement additional methods on `PRNG`
 
 
 ## License

Sources/Complex+Math.swift

@@ -1,484 +0,0 @@
-//
-//  Complex+Math.swift
-//  NumericAnnex
-//
-//  Created by Xiaodi Wu on 4/2/17.
-//
-//  Note
-//  ====
-//
-//  For maximum consistency with corresponding functions in C/C++, checks for
-//  special values in `naturalExponential()`, `squareRoot()`, trigonometric
-//  functions, and hyperbolic functions are adapted from libc++.
-//
-//  Code in libc++ is dual-licensed under the MIT and UIUC/NCSA licenses.
-//  Copyright © 2009-2017 contributors to the LLVM/libc++ project.
-
-extension Complex : Numeric {
-  @_transparent // @_inlineable
-  public init?<U>(exactly source: U) where U : BinaryInteger {
-    guard let t = T(exactly: source) else { return nil }
-    self.real = t
-    self.imaginary = 0
-  }
-
-  @_transparent // @_inlineable
-  public static func + (lhs: Complex, rhs: Complex) -> Complex {
-    return Complex(
-      real: lhs.real + rhs.real, imaginary: lhs.imaginary + rhs.imaginary
-    )
-  }
-
-  @_transparent // @_inlineable
-  public static func += (lhs: inout Complex, rhs: Complex) {
-    lhs.real += rhs.real
-    lhs.imaginary += rhs.imaginary
-  }
-
-  @_transparent // @_inlineable
-  public static func - (lhs: Complex, rhs: Complex) -> Complex {
-    return Complex(
-      real: lhs.real - rhs.real, imaginary: lhs.imaginary - rhs.imaginary
-    )
-  }
-
-  @_transparent // @_inlineable
-  public static func -= (lhs: inout Complex, rhs: Complex) {
-    lhs.real -= rhs.real
-    lhs.imaginary -= rhs.imaginary
-  }
-
-  @_transparent // @_inlineable
-  public static func * (lhs: Complex, rhs: Complex) -> Complex {
-    return Complex(
-      real: lhs.real * rhs.real - lhs.imaginary * rhs.imaginary,
-      imaginary: lhs.real * rhs.imaginary + lhs.imaginary * rhs.real
-    )
-  }
-
-  @_transparent // @_inlineable
-  public static func *= (lhs: inout Complex, rhs: Complex) {
-    let t = lhs.real
-    lhs.real = lhs.real * rhs.real - lhs.imaginary * rhs.imaginary
-    lhs.imaginary = t * rhs.imaginary + lhs.imaginary * rhs.real
-  }
-}
-
-extension Complex : SignedNumeric {
-  @_transparent // @_inlineable
-  public static prefix func - (operand: Complex) -> Complex {
-    return Complex(real: -operand.real, imaginary: -operand.imaginary)
-  }
-
-  @_transparent // @_inlineable
-  public mutating func negate() {
-    real.negate()
-    imaginary.negate()
-  }
-}
-
-extension Complex : Math {
-  @_transparent // @_inlineable
-  public static var pi: Complex {
-    return Complex(real: .pi)
-  }
-
-  @_transparent // @_inlineable
-  public static var e: Complex {
-    return Complex(real: .e)
-  }
-
-  @_transparent // @_inlineable
-  public static var phi: Complex {
-    return Complex(real: .phi)
-  }
-
-  @_transparent // @_inlineable
-  public static func / (lhs: Complex, rhs: Complex) -> Complex {
-    // Prevent avoidable overflow; see Numerical Recipes.
-    if rhs.real.magnitude >= rhs.imaginary.magnitude {
-      let ratio = rhs.imaginary / rhs.real
-      let denominator = rhs.real + rhs.imaginary * ratio
-      return Complex(
-        real: (lhs.real + lhs.imaginary * ratio) / denominator,
-        imaginary: (lhs.imaginary - lhs.real * ratio) / denominator
-      )
-    }
-    let ratio = rhs.real / rhs.imaginary
-    let denominator = rhs.real * ratio + rhs.imaginary
-    return Complex(
-      real: (lhs.real * ratio + lhs.imaginary) / denominator,
-      imaginary: (lhs.imaginary * ratio - lhs.real) / denominator
-    )
-
-    /*
-    let denominator = rhs.squaredMagnitude
-    return Complex(
-      real:
-        (lhs.real * rhs.real + lhs.imaginary * rhs.imaginary) / denominator,
-      imaginary:
-        (lhs.imaginary * rhs.real - lhs.real * rhs.imaginary) / denominator
-    )
-    */
-  }
-
-  @_transparent // @_inlineable
-  public static func /= (lhs: inout Complex, rhs: Complex) {
-    // Prevent avoidable overflow; see Numerical Recipes.
-    let t = lhs.real
-    if rhs.real.magnitude >= rhs.imaginary.magnitude {
-      let ratio = rhs.imaginary / rhs.real
-      let denominator = rhs.real + rhs.imaginary * ratio
-      lhs.real = (lhs.real + lhs.imaginary * ratio) / denominator
-      lhs.imaginary = (lhs.imaginary - t * ratio) / denominator
-    } else {
-      let ratio = rhs.real / rhs.imaginary
-      let denominator = rhs.real * ratio + rhs.imaginary
-      lhs.real = (lhs.real * ratio + lhs.imaginary) / denominator
-      lhs.imaginary = (lhs.imaginary * ratio - t) / denominator
-    }
-  }
-
-  // @_transparent // @_inlineable
-  public func naturalExponential() -> Complex {
-    if real.isNaN && imaginary == 0 { return self }
-    var im = imaginary
-    if real.isInfinite {
-      if real < 0 && !im.isFinite {
-        im = 1
-      } else if im == 0 || !im.isFinite {
-        if im.isInfinite { im = .nan }
-        return Complex(real: real, imaginary: im)
-      }
-    }
-    return Complex(r: T.exp(real), theta: im)
-  }
-
-  @_transparent // @_inlineable
-  public func naturalLogarithm() -> Complex {
-    return Complex(real: T.log(magnitude), imaginary: argument)
-  }
-
-  @_transparent // @_inlineable
-  public func commonLogarithm() -> Complex {
-    return Complex.log(self) / Complex.log(10 as Complex)
-  }
-
-  // @_transparent // @_inlineable
-  public func squareRoot() -> Complex {
-    if imaginary.isInfinite {
-      return Complex(real: .infinity, imaginary: imaginary)
-    }
-    if real.isInfinite {
-      if real > 0 {
-        return Complex(
-          real: real,
-          imaginary: imaginary.isNaN ? imaginary :
-            T(signOf: imaginary, magnitudeOf: 0)
-        )
-      }
-      return Complex(
-        real: imaginary.isNaN ? imaginary : 0,
-        imaginary: T(signOf: imaginary, magnitudeOf: real)
-      )
-    }
-    // Guard intermediate results and enforce a branch cut; see Numerical
-    // Recipes.
-    if real == 0 && imaginary == 0 { return 0 }
-    let w: T
-    let x = abs(real), y = abs(imaginary)
-    if x >= y {
-      let r = y / x
-      w = T.sqrt(x) * T.sqrt((1 + T.sqrt(1 + r * r)) / 2)
-    } else {
-      let r = x / y
-      w = T.sqrt(y) * T.sqrt((r + T.sqrt(1 + r * r)) / 2)
-    }
-    if real >= 0 {
-      return Complex(real: w, imaginary: imaginary / (2 * w))
-    }
-    /* real < 0 */
-    return Complex(real: y / (2 * w), imaginary: imaginary >= 0 ? w : -w)
-
-    /*
-    return Complex(r: T.sqrt(magnitude), theta: argument / 2)
-    */
-  }
-
-  @_transparent // @_inlineable
-  public func cubeRoot() -> Complex {
-    return Complex.exp(Complex.log(self) / 3)
-
-    /*
-    return Complex(r: T.cbrt(magnitude), theta: argument / 3)
-    */
-  }
-
-  @_transparent // @_inlineable
-  public func power(of base: Complex) -> Complex {
-    return Complex.exp(self * Complex.log(base))
-  }
-
-  @_transparent // @_inlineable
-  public func sine() -> Complex {
-    let sinh = Complex.sinh(Complex(real: -imaginary, imaginary: real))
-    return Complex(real: sinh.imaginary, imaginary: -sinh.real)
-  }
-
-  @_transparent // @_inlineable
-  public func cosine() -> Complex {
-    return Complex.cosh(Complex(real: -imaginary, imaginary: real))
-  }
-
-  @_transparent // @_inlineable
-  public func tangent() -> Complex {
-    let tanh = Complex.tanh(Complex(real: -imaginary, imaginary: real))
-    return Complex(real: tanh.imaginary, imaginary: -tanh.real)
-  }
-
-  @_transparent // @_inlineable
-  public func inverseSine() -> Complex {
-    let asinh = Complex.asinh(Complex(real: -imaginary, imaginary: real))
-    return Complex(real: asinh.imaginary, imaginary: -asinh.real)
-  }
-
-  // @_transparent // @_inlineable
-  public func inverseCosine() -> Complex {
-    if real.isNaN {
-      if imaginary.isInfinite {
-        return Complex(real: real, imaginary: -imaginary)
-      }
-      return Complex(real: real, imaginary: real)
-    }
-    if real.isInfinite {
-      if imaginary.isNaN {
-        return Complex(real: imaginary, imaginary: real)
-      }
-      if imaginary.isInfinite {
-        return Complex(
-          real: real < 0 ? .pi * 3 / 4 : .pi / 4,
-          imaginary: -imaginary
-        )
-      }
-      return real < 0 ?
-        Complex(real: .pi, imaginary: imaginary.sign == .minus ? -real : real) :
-        Complex(real: 0, imaginary: imaginary.sign == .minus ? real : -real)
-    }
-    if real == 0 && (imaginary.isNaN || imaginary == 0) {
-      return Complex(real: .pi / 2, imaginary: -imaginary)
-    }
-    if imaginary.isInfinite {
-      return Complex(real: .pi / 2, imaginary: -imaginary)
-    }
-    let a = Complex.pow(self, 2) - 1
-    let b = Complex.log(self + Complex.sqrt(a))
-    return Complex(
-      real: abs(b.imaginary),
-      imaginary: imaginary.sign == .minus ? abs(b.real) : -abs(b.real)
-    )
-  }
-
-  @_transparent // @_inlineable
-  public func inverseTangent() -> Complex {
-    let atanh = Complex.atanh(Complex(real: -imaginary, imaginary: real))
-    return Complex(real: atanh.imaginary, imaginary: -atanh.real)
-  }
-
-  // @_transparent // @_inlineable
-  public func hyperbolicSine() -> Complex {
-    if (real.isInfinite || real == 0) && !imaginary.isFinite {
-      return Complex(real: real, imaginary: .nan)
-    }
-    if imaginary == 0 && !real.isFinite {
-      return self
-    }
-    return Complex(
-      real: T.sinh(real) * T.cos(imaginary),
-      imaginary: T.cosh(real) * T.sin(imaginary)
-    )
-  }
-
-  // @_transparent // @_inlineable
-  public func hyperbolicCosine() -> Complex {
-    if real.isInfinite && !imaginary.isFinite {
-      return Complex(real: abs(real), imaginary: .nan)
-    }
-    if real == 0 && imaginary == 0 {
-      return Complex(real: 1, imaginary: imaginary)
-    }
-    if real == 0 && !imaginary.isFinite {
-      return Complex(real: .nan, imaginary: real)
-    }
-    if imaginary == 0 && !real.isFinite {
-      return Complex(real: abs(real), imaginary: imaginary)
-    }
-    return Complex(
-      real: T.cosh(real) * T.cos(imaginary),
-      imaginary: T.sinh(real) * T.sin(imaginary)
-    )
-  }
-
-  // @_transparent // @_inlineable
-  public func hyperbolicTangent() -> Complex {
-    if real.isNaN && imaginary == 0 { return self }
-    if real.isInfinite {
-      if !imaginary.isFinite { return 1 }
-      return Complex(
-        real: 1, imaginary: T(signOf: T.sin(2 * imaginary), magnitudeOf: 0)
-      )
-    }
-    // See AMS55 4.5.51
-    let twiceReal = 2 * real, twiceImaginary = 2 * imaginary
-    let denominator = T.cosh(twiceReal) + T.cos(twiceImaginary)
-    let sinh = T.sinh(twiceReal)
-    if sinh.isInfinite && denominator.isInfinite {
-      return Complex(
-        real: sinh > 0 ? (1 as T) : -(1 as T),
-        imaginary: twiceImaginary > 0 ? (0 as T) : -(0 as T)
-      )
-    }
-    return Complex(
-      real: sinh / denominator,
-      imaginary: T.sin(twiceImaginary) / denominator
-    )
-  }
-
-  // @_transparent // @_inlineable
-  public func inverseHyperbolicSine() -> Complex {
-    if real.isNaN {
-      if imaginary.isInfinite {
-        return Complex(real: imaginary, imaginary: real)
-      }
-      if imaginary == 0 { return self }
-      return Complex(real: real, imaginary: real)
-    }
-    if real.isInfinite {
-      if imaginary.isNaN { return self }
-      if imaginary.isInfinite {
-        return Complex(
-          real: real, imaginary: T(signOf: imaginary, magnitudeOf: .pi / 4)
-        )
-      }
-      return Complex(
-        real: real, imaginary: T(signOf: imaginary, magnitudeOf: 0)
-      )
-    }
-    if imaginary.isInfinite {
-      return Complex(
-        real: T(signOf: real, magnitudeOf: imaginary),
-        imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
-      )
-    }
-    let a = Complex.pow(self, 2) + 1
-    let b = Complex.log(self + Complex.sqrt(a))
-    return Complex(
-      real: T(signOf: real, magnitudeOf: b.real),
-      imaginary: T(signOf: imaginary, magnitudeOf: b.imaginary)
-    )
-  }
-
-  // @_transparent // @_inlineable
-  public func inverseHyperbolicCosine() -> Complex {
-    if real.isNaN {
-      if imaginary.isInfinite {
-        return Complex(real: .infinity, imaginary: real)
-      }
-      return Complex(real: real, imaginary: real)
-    }
-    if real.isInfinite {
-      if imaginary.isNaN {
-        return Complex(real: .infinity, imaginary: imaginary)
-      }
-      if imaginary.isInfinite {
-        switch (real.sign, imaginary.sign) {
-        case (.plus, .plus):
-          return Complex(real: .infinity, imaginary: .pi / 4)
-        case (.plus, .minus):
-          return Complex(real: .infinity, imaginary: -.pi / 4)
-        case (.minus, .plus):
-          return Complex(real: .infinity, imaginary: .pi * 3 / 4)
-        case (.minus, .minus):
-          return Complex(real: .infinity, imaginary: -.pi * 3 / 4)
-        }
-      }
-      switch (real.sign, imaginary.sign) {
-      case (.plus, .plus):
-        return Complex(real: .infinity, imaginary: 0)
-      case (.plus, .minus):
-        return Complex(real: .infinity, imaginary: -(0 as T))
-      case (.minus, .plus):
-        return Complex(real: .infinity, imaginary: .pi)
-      case (.minus, .minus):
-        return Complex(real: .infinity, imaginary: -.pi)
-      }
-    }
-    if imaginary.isNaN {
-      // See C11 DR 471.
-      return Complex(
-        real: imaginary,
-        imaginary: real == 0 ? .pi / 2 : imaginary
-      )
-    }
-    if imaginary.isInfinite {
-      return Complex(
-        real: .infinity, imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
-      )
-    }
-    let a = Complex.pow(self, 2) - 1
-    let b = Complex.log(self + Complex.sqrt(a))
-    return Complex(
-      real: T(signOf: 0, magnitudeOf: b.real),
-      imaginary: T(signOf: imaginary, magnitudeOf: b.imaginary)
-    )
-  }
-
-  // @_transparent // @_inlineable
-  public func inverseHyperbolicTangent() -> Complex {
-    if imaginary.isNaN {
-      if real.isInfinite || real == 0 {
-        return Complex(
-          real: T(signOf: real, magnitudeOf: 0), imaginary: imaginary
-        )
-      }
-      return Complex(real: imaginary, imaginary: imaginary)
-    }
-    if imaginary.isInfinite {
-      return Complex(
-        real: T(signOf: real, magnitudeOf: 0),
-        imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
-      )
-    }
-    if imaginary == 0 && abs(real) == 1 {
-      return Complex(
-        real: T(signOf: real, magnitudeOf: .infinity),
-        imaginary: T(signOf: imaginary, magnitudeOf: 0)
-      )
-    }
-    if real.isNaN {
-      return Complex(real: real, imaginary: real)
-    }
-    if real.isInfinite {
-      return Complex(
-        real: T(signOf: real, magnitudeOf: 0),
-        imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
-      )
-    }
-    let a = Complex.log((1 + self) / (1 - self)) / 2
-    return Complex(
-      real: T(signOf: real, magnitudeOf: a.real),
-      imaginary: T(signOf: imaginary, magnitudeOf: a.imaginary)
-    )
-  }
-}
-
-/// Returns the absolute value (magnitude, modulus) of `z`.
-@_transparent
-public func abs<T>(_ z: Complex<T>) -> Complex<T> {
-  return Complex(real: z.magnitude)
-}
-
-/// Returns the square root of `z`.
-@_transparent
-public func sqrt<T>(_ z: Complex<T>) -> Complex<T> {
-  return z.squareRoot()
-}

Sources/Complex.swift

@@ -4,11 +4,23 @@
 //
 //  Created by Xiaodi Wu on 3/25/17.
 //
+//  Note
+//  ====
+//
+//  For maximum consistency with corresponding functions in C/C++, checks for
+//  special values in `naturalExponential()`, `squareRoot()`, trigonometric
+//  functions, and hyperbolic functions are adapted from libc++.
+//
+//  Code in libc++ is dual-licensed under the MIT and UIUC/NCSA licenses.
+//  Copyright © 2009-2017 contributors to the LLVM/libc++ project.
 
 /// A type to represent a complex value in Cartesian form.
 ///
+/// - Note: `Complex64` is a type alias for `Complex<Float>` and `Complex128` is
+///   a type alias for `Complex<Double>`.
+///
 /// Create new instances of `Complex<T>` using integer or floating-point
-/// literals and the imaginary unit `.i`. For example:
+/// literals and the imaginary unit `Complex<T>.i`. For example:
 ///
 /// ```swift
 /// let x = 2 + 4 * .i // `x` is of type `Complex<Double>`
@@ -18,7 +30,7 @@
 /// ```
 ///
 /// Additional Considerations
-/// =========================
+/// -------------------------
 ///
 /// Floating-point types have special values that represent infinity or NaN
 /// ("not a number"). Complex functions in different languages may return
@@ -40,14 +52,22 @@
 /// [dfn]: http://mathworld.wolfram.com/BranchCut.html
 /// [std]: http://www.open-std.org/JTC1/SC22/WG14/www/standards.html#9899
 @_fixed_layout
-public struct Complex<T : FloatingPointMath>
-where T : _ExpressibleByBuiltinFloatLiteral {
+public struct Complex<T : Real> : Codable
+where T : Codable & _ExpressibleByBuiltinFloatLiteral {
+  // ---------------------------------------------------------------------------
+  // MARK: Stored Properties
+  // ---------------------------------------------------------------------------
+
   /// The real component of the complex value.
   public var real: T
 
   /// The imaginary component of the complex value.
   public var imaginary: T
 
+  // ---------------------------------------------------------------------------
+  // MARK: Initializers
+  // ---------------------------------------------------------------------------
+
   /// Creates a new value from the given real and imaginary components.
   ///
   /// - Parameters:
@@ -59,14 +79,17 @@ where T : _ExpressibleByBuiltinFloatLiteral {
     self.imaginary = imaginary
   }
 
-  /// The imaginary unit i.
+  /// Creates a new value from the given polar coordinates `(r, theta)`.
+  ///
+  /// - Parameters:
+  ///   - r: The new value's radial coordinate.
+  ///   - theta: The new value's angular coordinate.
   @_transparent // @_inlineable
-  public static var i: Complex {
-    return Complex(real: 0, imaginary: 1)
+  public init(r: T, theta: T) {
+    self.real = r * T.cos(theta)
+    self.imaginary = r * T.sin(theta)
   }
-}
 
-extension Complex {
   /// Creates a new value from the given real component, rounded to the closest
   /// possible representation.
   ///
@@ -214,9 +237,16 @@ extension Complex {
 }
 
 extension Complex where T : BinaryFloatingPoint {
+  // ---------------------------------------------------------------------------
+  // MARK: Initializers (Constrained)
+  // ---------------------------------------------------------------------------
+
   /// Creates a new value from the given real component, rounded to the closest
   /// possible representation.
   ///
+  /// - Note: This initializer creates only instances of
+  ///   `Complex<T> where T : BinaryFloatingPoint`.
+  ///
   /// - Parameters:
   ///   - real: The value to convert to a real component of type `T`.
   @_transparent // @_inlineable
@@ -228,6 +258,9 @@ extension Complex where T : BinaryFloatingPoint {
   /// Creates a new value from the given real component, rounded to the closest
   /// possible representation.
   ///
+  /// - Note: This initializer creates only instances of
+  ///   `Complex<T> where T : BinaryFloatingPoint`.
+  ///
   /// - Parameters:
   ///   - real: The value to convert to a real component of type `T`.
   @_transparent // @_inlineable
@@ -244,17 +277,20 @@ extension Complex where T : BinaryFloatingPoint {
 }
 
 extension Complex {
-  /// Creates a new value from the given polar coordinates `(r, theta)`.
-  ///
-  /// - Parameters:
-  ///   - r: The new value's radial coordinate.
-  ///   - theta: The new value's angular coordinate.
+  // ---------------------------------------------------------------------------
+  // MARK: Static Properties
+  // ---------------------------------------------------------------------------
+
+  /// The imaginary unit _i_.
   @_transparent // @_inlineable
-  public init(r: T, theta: T) {
-    self.real = r * T.cos(theta)
-    self.imaginary = r * T.sin(theta)
+  public static var i: Complex {
+    return Complex(real: 0, imaginary: 1)
   }
 
+  // ---------------------------------------------------------------------------
+  // MARK: Computed Properties
+  // ---------------------------------------------------------------------------
+
   /// The principal argument (phase angle) of this value.
   ///
   /// Special cases are handled as if calling `T.atan2(imaginary, real)`. The
@@ -264,27 +300,35 @@ extension Complex {
     return T.atan2(imaginary, real)
   }
 
-  /// A Boolean value indicating whether the instance's real and imaginary
-  /// components are both in canonical form, as defined in the [IEEE 754
-  /// specification][spec]. Every `Float` or `Double` value is canonical.
+#if false
+  /// A Boolean value indicating whether the instance's representation is in
+  /// canonical form.
+  ///
+  /// A complex value is represented in canonical form if the its real and
+  /// imaginary components are both represented in canonical form, as defined in
+  /// the [IEEE 754 specification][spec]. Every `Float` or `Double` value is
+  /// canonical.
   ///
   /// [spec]: http://ieeexplore.ieee.org/servlet/opac?punumber=4610933
   @_transparent // @_inlineable
   public var isCanonical: Bool {
     return real.isCanonical && imaginary.isCanonical
   }
+#endif
 
-  /// A Boolean value indicating whether the instance's real and imaginary
-  /// components are both finite.
+  /// A Boolean value indicating whether the instance is finite.
   ///
-  /// All values other than NaN and infinity are considered finite.
+  /// A complex value is finite if its real and imaginary components are both
+  /// finite. A component is finite if it is not infinity or NaN.
   @_transparent // @_inlineable
   public var isFinite: Bool {
     return real.isFinite && imaginary.isFinite
   }
 
-  /// A Boolean value indicating whether the instance's real and/or imaginary
-  /// components are infinite.
+  /// A Boolean value indicating whether the instance is infinite.
+  ///
+  /// A complex value is infinite if at least one of its components (real or
+  /// imaginary) is infinite, even if the other component is NaN.
   ///
   /// Note that `isFinite` and `isInfinite` do not form a dichotomy because NaN
   /// is neither finite nor infinite.
@@ -293,8 +337,10 @@ extension Complex {
     return real.isInfinite || imaginary.isInfinite
   }
 
-  /// A Boolean value indicating whether the instance's real and/or imaginary
-  /// components are NaN ("not a number").
+  /// A Boolean value indicating whether the instance is NaN ("not a number").
+  ///
+  /// A complex value is NaN if at least one of its components (real or
+  /// imaginary) is NaN and the other component is not infinite.
   ///
   /// Because NaN is not equal to any value, including NaN, use this property
   /// instead of the equal-to operator (`==`) or not-equal-to operator (`!=`) to
@@ -303,23 +349,27 @@ extension Complex {
   /// This property is `true` for both quiet and signaling NaNs.
   @_transparent // @_inlineable
   public var isNaN: Bool {
-    return real.isNaN || imaginary.isNaN
+    return (real.isNaN && !imaginary.isInfinite) ||
+      (imaginary.isNaN && !real.isInfinite)
   }
 
-  /// A Boolean value indicating whether the instance's real and/or imaginary
-  /// components are signaling NaNs.
+  /// A Boolean value indicating whether the instance is a signaling NaN.
+  ///
+  /// A complex value is a signaling NaN if at least one of its components (real
+  /// or imaginary) is a signaling NaN and the other component is not infinite.
   ///
   /// Signaling NaNs typically raise the Invalid flag when used in general
   /// computing operations.
   @_transparent // @_inlineable
   public var isSignalingNaN: Bool {
-    return real.isSignalingNaN || imaginary.isSignalingNaN
+    return (real.isSignalingNaN && !imaginary.isInfinite) ||
+      (imaginary.isSignalingNaN && !real.isInfinite)
   }
 
   /// A Boolean value indicating whether the instance is equal to zero.
   ///
-  /// The `isZero` property of a value `z` is `true` when both real and
-  /// imaginary components represent either `-0.0` or `+0.0`.
+  /// A complex value is equal to zero if its real and imaginary components both
+  /// represent either `-0.0` or `+0.0`.
   @_transparent // @_inlineable
   public var isZero: Bool {
     return real.isZero && imaginary.isZero
@@ -333,13 +383,6 @@ extension Complex {
     return T.hypot(real, imaginary)
   }
 
-  /// The polar coordinates representing this value, equivalent to
-  /// `(r: magnitude, theta: argument)`.
-  @_transparent // @_inlineable
-  public var polar: (r: T, theta: T) {
-    return (r: magnitude, theta: argument)
-  }
-
   /// The squared magnitude (field norm, absolute square) of this value.
   ///
   /// This is less costly to compute than `magnitude` and, in some cases, can be
@@ -352,6 +395,10 @@ extension Complex {
     return real * real + imaginary * imaginary
   }
 
+  // ---------------------------------------------------------------------------
+  // MARK: Methods
+  // ---------------------------------------------------------------------------
+
   /// Returns the complex conjugate of this value, obtained by reversing the
   /// sign of the imaginary component.
   @_transparent // @_inlineable
@@ -375,7 +422,8 @@ extension Complex {
     return self
   }
 
-  /// Returns the reciprocal (multiplicative inverse) of this value.
+  /// Returns the reciprocal (multiplicative inverse) of this value, or NaN
+  /// ("not a number") if this value is infinite, NaN, or zero.
   @_transparent // @_inlineable
   public func reciprocal() -> Complex {
     let denominator = squaredMagnitude
@@ -386,6 +434,10 @@ extension Complex {
 }
 
 extension Complex : ExpressibleByFloatLiteral {
+  // ---------------------------------------------------------------------------
+  // MARK: ExpressibleByFloatLiteral
+  // ---------------------------------------------------------------------------
+
   @_transparent // @_inlineable
   public init(floatLiteral value: T) {
     self.real = value
@@ -394,6 +446,10 @@ extension Complex : ExpressibleByFloatLiteral {
 }
 
 extension Complex : ExpressibleByIntegerLiteral {
+  // ---------------------------------------------------------------------------
+  // MARK: ExpressibleByIntegerLiteral
+  // ---------------------------------------------------------------------------
+
   @_transparent // @_inlineable
   public init(integerLiteral value: Int) {
     self.real = T(value)
@@ -402,13 +458,30 @@ extension Complex : ExpressibleByIntegerLiteral {
 }
 
 extension Complex : CustomStringConvertible {
+  // ---------------------------------------------------------------------------
+  // MARK: CustomStringConvertible
+  // ---------------------------------------------------------------------------
+
   @_transparent // @_inlineable
   public var description: String {
-    return "\(real) + \(imaginary)i"
+    return real.isNaN && real.sign == .minus
+      // At present, -NaN is described as "nan", which is acceptable for real
+      // values. However, it is arguably misleading to describe -NaN - NaNi as
+      // "nan + nani" or "nan - nani". Therefore, handle this case separately.
+      ? imaginary.sign == .minus
+        ? "-\(-real) - \(-imaginary)i"
+        : "-\(-real) + \(imaginary)i"
+      : imaginary.sign == .minus
+        ? "\(real) - \(-imaginary)i"
+        : "\(real) + \(imaginary)i"
   }
 }
 
 extension Complex : Equatable {
+  // ---------------------------------------------------------------------------
+  // MARK: Equatable
+  // ---------------------------------------------------------------------------
+
   @_transparent // @_inlineable
   public static func == (lhs: Complex, rhs: Complex) -> Bool {
     return lhs.real == rhs.real && lhs.imaginary == rhs.imaginary
@@ -416,11 +489,564 @@ extension Complex : Equatable {
 }
 
 extension Complex : Hashable {
+  // ---------------------------------------------------------------------------
+  // MARK: Hashable
+  // ---------------------------------------------------------------------------
+
   // @_transparent // @_inlineable
   public var hashValue: Int {
-    return _fnv1a(real, imaginary)
+    return _Hash._combine(real, imaginary)
   }
 }
 
+extension Complex : Numeric {
+  // ---------------------------------------------------------------------------
+  // MARK: Numeric
+  // ---------------------------------------------------------------------------
+
+  // @_transparent // @_inlineable
+  public init?<U>(exactly source: U) where U : BinaryInteger {
+    guard let t = T(exactly: source) else { return nil }
+    self.real = t
+    self.imaginary = 0
+  }
+
+  @_transparent // @_inlineable
+  public static func + (lhs: Complex, rhs: Complex) -> Complex {
+    var lhs = lhs
+    lhs += rhs
+    return lhs
+  }
+
+  @_transparent // @_inlineable
+  public static func += (lhs: inout Complex, rhs: Complex) {
+    lhs.real += rhs.real
+    lhs.imaginary += rhs.imaginary
+  }
+
+  @_transparent // @_inlineable
+  public static func - (lhs: Complex, rhs: Complex) -> Complex {
+    var lhs = lhs
+    lhs -= rhs
+    return lhs
+  }
+
+  @_transparent // @_inlineable
+  public static func -= (lhs: inout Complex, rhs: Complex) {
+    lhs.real -= rhs.real
+    lhs.imaginary -= rhs.imaginary
+  }
+
+  @_transparent // @_inlineable
+  public static func * (lhs: Complex, rhs: Complex) -> Complex {
+    var a = lhs.real, b = lhs.imaginary, c = rhs.real, d = rhs.imaginary
+    let ac = a * c, bd = b * d, ad = a * d, bc = b * c
+    let x = ac - bd
+    let y = ad + bc
+    // Recover infinities that computed as NaN + iNaN.
+    // See C11 Annex G.
+    if x.isNaN && y.isNaN {
+      var recalculate = false
+      if a.isInfinite || b.isInfinite {
+        // "Box" the infinity and change NaNs in the other operand to 0.
+        a = T(signOf: a, magnitudeOf: a.isInfinite ? 1 : 0)
+        b = T(signOf: b, magnitudeOf: b.isInfinite ? 1 : 0)
+        if c.isNaN { c = T(signOf: c, magnitudeOf: 0) }
+        if d.isNaN { d = T(signOf: d, magnitudeOf: 0) }
+        recalculate = true
+      }
+      if c.isInfinite || d.isInfinite {
+        // "Box" the infinity and change NaNs in the other operand to 0.
+        if a.isNaN { a = T(signOf: a, magnitudeOf: 0) }
+        if b.isNaN { b = T(signOf: b, magnitudeOf: 0) }
+        c = T(signOf: c, magnitudeOf: c.isInfinite ? 1 : 0)
+        d = T(signOf: d, magnitudeOf: d.isInfinite ? 1 : 0)
+        recalculate = true
+      }
+      if !recalculate &&
+        (ac.isInfinite || bd.isInfinite || ad.isInfinite || bc.isInfinite) {
+        // Recover infinities from overflow by changing NaNs to 0.
+        if a.isNaN { a = T(signOf: a, magnitudeOf: 0) }
+        if b.isNaN { b = T(signOf: b, magnitudeOf: 0) }
+        if c.isNaN { c = T(signOf: c, magnitudeOf: 0) }
+        if d.isNaN { d = T(signOf: d, magnitudeOf: 0) }
+        recalculate = true
+      }
+      if recalculate {
+        return Complex(
+          real: .infinity * (a * c - b * d),
+          imaginary: .infinity * (a * d + b * c)
+        )
+      }
+    }
+    return Complex(real: x, imaginary: y)
+    /*
+    return Complex(
+      real: lhs.real * rhs.real - lhs.imaginary * rhs.imaginary,
+      imaginary: lhs.real * rhs.imaginary + lhs.imaginary * rhs.real
+    )
+    */
+  }
+
+  @_transparent // @_inlineable
+  public static func *= (lhs: inout Complex, rhs: Complex) {
+    lhs = lhs * rhs
+  }
+}
+
+extension Complex : SignedNumeric {
+  // ---------------------------------------------------------------------------
+  // MARK: SignedNumeric
+  // ---------------------------------------------------------------------------
+
+  @_transparent // @_inlineable
+  public static prefix func - (operand: Complex) -> Complex {
+    return Complex(real: -operand.real, imaginary: -operand.imaginary)
+  }
+
+  @_transparent // @_inlineable
+  public mutating func negate() {
+    real.negate()
+    imaginary.negate()
+  }
+}
+
+extension Complex : Math {
+  // ---------------------------------------------------------------------------
+  // MARK: Math
+  // ---------------------------------------------------------------------------
+
+  @_transparent // @_inlineable
+  public static var pi: Complex {
+    return Complex(real: .pi)
+  }
+
+  @_transparent // @_inlineable
+  public static var e: Complex {
+    return Complex(real: .e)
+  }
+
+  @_transparent // @_inlineable
+  public static var phi: Complex {
+    return Complex(real: .phi)
+  }
+
+  @_transparent // @_inlineable
+  public static func / (lhs: Complex, rhs: Complex) -> Complex {
+    var a = lhs.real, b = lhs.imaginary, c = rhs.real, d = rhs.imaginary
+    var x: T
+    var y: T
+    // Prevent avoidable overflow; see Numerical Recipes.
+    if c.magnitude >= d.magnitude {
+      let ratio = d / c
+      let denominator = c + d * ratio
+      x = (a + b * ratio) / denominator
+      y = (b - a * ratio) / denominator
+    } else {
+      let ratio = c / d
+      let denominator = c * ratio + d
+      x = (a * ratio + b) / denominator
+      y = (b * ratio - a) / denominator
+    }
+    // Recover infinities and zeros that computed as NaN + iNaN.
+    // See C11 Annex G.
+    if x.isNaN && y.isNaN {
+      if c == 0 && d == 0 && (!a.isNaN || !b.isNaN) {
+        x = T(signOf: c, magnitudeOf: .infinity) * a
+        y = T(signOf: c /* sic */, magnitudeOf: .infinity) * b
+      } else if (a.isInfinite || b.isInfinite) && c.isFinite && d.isFinite {
+        a = T(signOf: a, magnitudeOf: a.isInfinite ? 1 : 0)
+        b = T(signOf: b, magnitudeOf: b.isInfinite ? 1 : 0)
+        x = .infinity * (a * c + b * d)
+        y = .infinity * (b * c - a * d)
+      } else if (c.isInfinite || d.isInfinite) && a.isFinite && b.isFinite {
+        c = T(signOf: c, magnitudeOf: c.isInfinite ? 1 : 0)
+        d = T(signOf: d, magnitudeOf: d.isInfinite ? 1 : 0)
+        x = 0 * (a * c + b * d)
+        y = 0 * (b * c - a * d)
+      }
+    }
+    return Complex(real: x, imaginary: y)
+    /*
+    let denominator = rhs.squaredMagnitude
+    return Complex(
+      real:
+        (lhs.real * rhs.real + lhs.imaginary * rhs.imaginary) / denominator,
+      imaginary:
+        (lhs.imaginary * rhs.real - lhs.real * rhs.imaginary) / denominator
+    )
+    */
+  }
+
+  @_transparent // @_inlineable
+  public static func /= (lhs: inout Complex, rhs: Complex) {
+    lhs = lhs / rhs
+  }
+
+  @_transparent // @_inlineable
+  public static func ** (lhs: Complex, rhs: Complex) -> Complex {
+    return Complex.exp(rhs * Complex.log(lhs))
+  }
+
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout Complex, rhs: Complex) {
+    lhs = lhs ** rhs
+  }
+
+  @_transparent // @_inlineable
+  public func naturalExponential() -> Complex {
+    if real.isNaN && imaginary == 0 { return self }
+    var im = imaginary
+    if real.isInfinite {
+      if real < 0 && !im.isFinite {
+        im = 1
+      } else if im == 0 || !im.isFinite {
+        if im.isInfinite { im = .nan }
+        return Complex(real: real, imaginary: im)
+      }
+    }
+    return Complex(r: T.exp(real), theta: im)
+  }
+
+  @_transparent // @_inlineable
+  public func naturalLogarithm() -> Complex {
+    return Complex(real: T.log(magnitude), imaginary: argument)
+  }
+
+  @_transparent // @_inlineable
+  public func commonLogarithm() -> Complex {
+    return Complex.log(self) / Complex.log(10 as Complex)
+  }
+
+  @_transparent // @_inlineable
+  public func squareRoot() -> Complex {
+    if imaginary.isInfinite {
+      return Complex(real: .infinity, imaginary: imaginary)
+    }
+    if real.isInfinite {
+      if real > 0 {
+        return Complex(
+          real: real,
+          imaginary: imaginary.isNaN ? imaginary :
+            T(signOf: imaginary, magnitudeOf: 0)
+        )
+      }
+      return Complex(
+        real: imaginary.isNaN ? imaginary : 0,
+        imaginary: T(signOf: imaginary, magnitudeOf: real)
+      )
+    }
+    // Guard intermediate results and enforce a branch cut; see Numerical
+    // Recipes.
+    if real == 0 && imaginary == 0 { return 0 }
+    let w: T
+    let x = abs(real), y = abs(imaginary)
+    if x >= y {
+      let r = y / x
+      w = T.sqrt(x) * T.sqrt((1 + T.sqrt(1 + r * r)) / 2)
+    } else {
+      let r = x / y
+      w = T.sqrt(y) * T.sqrt((r + T.sqrt(1 + r * r)) / 2)
+    }
+    if real >= 0 {
+      return Complex(real: w, imaginary: imaginary / (2 * w))
+    }
+    /* real < 0 */
+    return Complex(real: y / (2 * w), imaginary: imaginary >= 0 ? w : -w)
+    /*
+    return Complex(r: T.sqrt(magnitude), theta: argument / 2)
+    */
+  }
+
+  @_transparent // @_inlineable
+  public func cubeRoot() -> Complex {
+    return Complex.exp(Complex.log(self) / 3)
+    /*
+    return Complex(r: T.cbrt(magnitude), theta: argument / 3)
+    */
+  }
+
+  @_transparent // @_inlineable
+  public func sine() -> Complex {
+    let sinh = Complex.sinh(Complex(real: -imaginary, imaginary: real))
+    return Complex(real: sinh.imaginary, imaginary: -sinh.real)
+  }
+
+  @_transparent // @_inlineable
+  public func cosine() -> Complex {
+    return Complex.cosh(Complex(real: -imaginary, imaginary: real))
+  }
+
+  @_transparent // @_inlineable
+  public func tangent() -> Complex {
+    let tanh = Complex.tanh(Complex(real: -imaginary, imaginary: real))
+    return Complex(real: tanh.imaginary, imaginary: -tanh.real)
+  }
+
+  @_transparent // @_inlineable
+  public func inverseSine() -> Complex {
+    let asinh = Complex.asinh(Complex(real: -imaginary, imaginary: real))
+    return Complex(real: asinh.imaginary, imaginary: -asinh.real)
+  }
+
+  @_transparent // @_inlineable
+  public func inverseCosine() -> Complex {
+    if real.isNaN {
+      if imaginary.isInfinite {
+        return Complex(real: real, imaginary: -imaginary)
+      }
+      return Complex(real: real, imaginary: real)
+    }
+    if real.isInfinite {
+      if imaginary.isNaN {
+        return Complex(real: imaginary, imaginary: real)
+      }
+      if imaginary.isInfinite {
+        return Complex(
+          real: real < 0 ? .pi * 3 / 4 : .pi / 4,
+          imaginary: -imaginary
+        )
+      }
+      return real < 0
+        ? Complex(real: .pi, imaginary: imaginary.sign == .minus ? -real : real)
+        : Complex(real: 0, imaginary: imaginary.sign == .minus ? real : -real)
+    }
+    if real == 0 && (imaginary.isNaN || imaginary == 0) {
+      return Complex(real: .pi / 2, imaginary: -imaginary)
+    }
+    if imaginary.isInfinite {
+      return Complex(real: .pi / 2, imaginary: -imaginary)
+    }
+    let a = self ** 2 - 1
+    let b = Complex.log(self + Complex.sqrt(a))
+    return Complex(
+      real: abs(b.imaginary),
+      imaginary: imaginary.sign == .minus ? abs(b.real) : -abs(b.real)
+    )
+  }
+
+  @_transparent // @_inlineable
+  public func inverseTangent() -> Complex {
+    let atanh = Complex.atanh(Complex(real: -imaginary, imaginary: real))
+    return Complex(real: atanh.imaginary, imaginary: -atanh.real)
+  }
+
+  @_transparent // @_inlineable
+  public func hyperbolicSine() -> Complex {
+    if (real.isInfinite || real == 0) && !imaginary.isFinite {
+      return Complex(real: real, imaginary: .nan)
+    }
+    if imaginary == 0 && !real.isFinite {
+      return self
+    }
+    return Complex(
+      real: T.sinh(real) * T.cos(imaginary),
+      imaginary: T.cosh(real) * T.sin(imaginary)
+    )
+  }
+
+  @_transparent // @_inlineable
+  public func hyperbolicCosine() -> Complex {
+    if real.isInfinite && !imaginary.isFinite {
+      return Complex(real: abs(real), imaginary: .nan)
+    }
+    if real == 0 && imaginary == 0 {
+      return Complex(real: 1, imaginary: imaginary)
+    }
+    if real == 0 && !imaginary.isFinite {
+      return Complex(real: .nan, imaginary: real)
+    }
+    if imaginary == 0 && !real.isFinite {
+      return Complex(real: abs(real), imaginary: imaginary)
+    }
+    return Complex(
+      real: T.cosh(real) * T.cos(imaginary),
+      imaginary: T.sinh(real) * T.sin(imaginary)
+    )
+  }
+
+  @_transparent // @_inlineable
+  public func hyperbolicTangent() -> Complex {
+    if real.isNaN && imaginary == 0 { return self }
+    if real.isInfinite {
+      if !imaginary.isFinite {
+        return Complex(real: T(signOf: real, magnitudeOf: 1), imaginary: 0)
+      }
+      return Complex(
+        real: T(signOf: real, magnitudeOf: 1),
+        imaginary: T(signOf: T.sin(2 * imaginary), magnitudeOf: 0)
+      )
+    }
+    if real == 0 {
+      // See C11 DR 471.
+      // http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471
+      return Complex(real: real, imaginary: T.tan(imaginary))
+    }
+    // See AMS55 4.5.51.
+    let twiceReal = 2 * real, twiceImaginary = 2 * imaginary
+    let sinhTwiceReal = T.sinh(twiceReal)
+    let denominator = T.cosh(twiceReal) + T.cos(twiceImaginary)
+    // If `sinh(2 * real)` is infinite, then `cosh(2 * real)` is infinite.
+    //
+    // Thus, `denominator` is either:
+    // - infinite, if `cos(2 * imaginary)` is not NaN; or
+    // - NaN, if `cos(2 * imaginary)` is NaN (i.e., if `2 * imaginary` is
+    //   infinite or NaN).
+    //
+    // In either case, if `imaginary` is finite, the complex hyperbolic tangent
+    // should be 1 or -1.
+    if sinhTwiceReal.isInfinite && imaginary.isFinite {
+      return Complex(
+        real: T(signOf: sinhTwiceReal, magnitudeOf: 1),
+        imaginary: T(signOf: T.sin(twiceImaginary), magnitudeOf: 0)
+      )
+    }
+    return Complex(
+      real: sinhTwiceReal / denominator,
+      imaginary: T.sin(twiceImaginary) / denominator
+    )
+  }
+
+  @_transparent // @_inlineable
+  public func inverseHyperbolicSine() -> Complex {
+    if real.isNaN {
+      if imaginary.isInfinite {
+        return Complex(real: imaginary, imaginary: real)
+      }
+      if imaginary == 0 { return self }
+      return Complex(real: real, imaginary: real)
+    }
+    if real.isInfinite {
+      if imaginary.isNaN { return self }
+      if imaginary.isInfinite {
+        return Complex(
+          real: real, imaginary: T(signOf: imaginary, magnitudeOf: .pi / 4)
+        )
+      }
+      return Complex(
+        real: real, imaginary: T(signOf: imaginary, magnitudeOf: 0)
+      )
+    }
+    if imaginary.isInfinite {
+      return Complex(
+        real: T(signOf: real, magnitudeOf: imaginary),
+        imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
+      )
+    }
+    let a = self ** 2 + 1
+    let b = Complex.log(self + Complex.sqrt(a))
+    return Complex(
+      real: T(signOf: real, magnitudeOf: b.real),
+      imaginary: T(signOf: imaginary, magnitudeOf: b.imaginary)
+    )
+  }
+
+  @_transparent // @_inlineable
+  public func inverseHyperbolicCosine() -> Complex {
+    if real.isNaN {
+      if imaginary.isInfinite {
+        return Complex(real: .infinity, imaginary: real)
+      }
+      return Complex(real: real, imaginary: real)
+    }
+    if real.isInfinite {
+      if imaginary.isNaN {
+        return Complex(real: .infinity, imaginary: imaginary)
+      }
+      if imaginary.isInfinite {
+        switch (real.sign, imaginary.sign) {
+        case (.plus, .plus):
+          return Complex(real: .infinity, imaginary: .pi / 4)
+        case (.plus, .minus):
+          return Complex(real: .infinity, imaginary: -.pi / 4)
+        case (.minus, .plus):
+          return Complex(real: .infinity, imaginary: .pi * 3 / 4)
+        case (.minus, .minus):
+          return Complex(real: .infinity, imaginary: -.pi * 3 / 4)
+        }
+      }
+      switch (real.sign, imaginary.sign) {
+      case (.plus, .plus):
+        return Complex(real: .infinity, imaginary: 0)
+      case (.plus, .minus):
+        return Complex(real: .infinity, imaginary: -(0 as T))
+      case (.minus, .plus):
+        return Complex(real: .infinity, imaginary: .pi)
+      case (.minus, .minus):
+        return Complex(real: .infinity, imaginary: -.pi)
+      }
+    }
+    if imaginary.isNaN {
+      if real == 0 {
+        // See C11 DR 471.
+        // http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1892.htm#dr_471
+        return Complex(real: imaginary, imaginary: .pi / 2)
+      }
+      return Complex(real: imaginary, imaginary: imaginary)
+    }
+    if imaginary.isInfinite {
+      return Complex(
+        real: .infinity, imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
+      )
+    }
+    let a = self ** 2 - 1
+    let b = Complex.log(self + Complex.sqrt(a))
+    return Complex(
+      real: T(signOf: 0, magnitudeOf: b.real),
+      imaginary: T(signOf: imaginary, magnitudeOf: b.imaginary)
+    )
+  }
+
+  @_transparent // @_inlineable
+  public func inverseHyperbolicTangent() -> Complex {
+    if imaginary.isNaN {
+      if real.isInfinite || real == 0 {
+        return Complex(
+          real: T(signOf: real, magnitudeOf: 0), imaginary: imaginary
+        )
+      }
+      return Complex(real: imaginary, imaginary: imaginary)
+    }
+    if imaginary.isInfinite {
+      return Complex(
+        real: T(signOf: real, magnitudeOf: 0),
+        imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
+      )
+    }
+    if imaginary == 0 && abs(real) == 1 {
+      return Complex(
+        real: T(signOf: real, magnitudeOf: .infinity),
+        imaginary: T(signOf: imaginary, magnitudeOf: 0)
+      )
+    }
+    if real.isNaN {
+      return Complex(real: real, imaginary: real)
+    }
+    if real.isInfinite {
+      return Complex(
+        real: T(signOf: real, magnitudeOf: 0),
+        imaginary: T(signOf: imaginary, magnitudeOf: .pi / 2)
+      )
+    }
+    let a = Complex.log((1 + self) / (1 - self)) / 2
+    return Complex(
+      real: T(signOf: real, magnitudeOf: a.real),
+      imaginary: T(signOf: imaginary, magnitudeOf: a.imaginary)
+    )
+  }
+}
+
+/// Returns the absolute value (magnitude, modulus) of `z`.
+@_transparent
+public func abs<T>(_ z: Complex<T>) -> Complex<T> {
+  return Complex(real: z.magnitude)
+}
+
+/// Returns the square root of `z`.
+@_transparent
+public func sqrt<T>(_ z: Complex<T>) -> Complex<T> {
+  return z.squareRoot()
+}
+
 public typealias Complex64 = Complex<Float>
 public typealias Complex128 = Complex<Double>

Sources/Exponentiation.swift

@@ -1,32 +0,0 @@
-//
-//  Exponentiation.swift
-//  NumericAnnex
-//
-//  Created by Xiaodi Wu on 4/24/17.
-//
-
-extension BinaryInteger {
-  /// Returns the result of raising `base` to the power of `exponent`, rounded
-  /// to a representable value.
-  // @_transparent // @_inlineable
-  public static func pow(_ base: Self, _ exponent: Self) -> Self {
-    var x = base, n = exponent
-    if Self.isSigned && n < 0 {
-      x = 1 / x
-      n = 0 - n
-    } else if n == 0 {
-      return 1
-    }
-    // Exponentiate by iterative squaring.
-    var y = 1 as Self
-    while n > 1 {
-      if n % 2 == 1 {
-        y *= x
-        n -= 1
-      }
-      x *= x
-      n /= 2
-    }
-    return x * y
-  }
-}

Sources/ExponentiationOperators.swift

@@ -0,0 +1,19 @@
+//
+//  ExponentiationOperators.swift
+//  NumericAnnex
+//
+//  Created by Xiaodi Wu on 8/26/17.
+//
+
+precedencegroup ExponentiationPrecedence {
+  associativity: right
+  higherThan: MultiplicationPrecedence
+}
+
+// "Exponentiative"
+
+infix operator ** : ExponentiationPrecedence
+
+// Compound
+
+infix operator **= : AssignmentPrecedence

Sources/Factoring.swift

@@ -1,115 +0,0 @@
-//
-//  Factoring.swift
-//  NumericAnnex
-//
-//  Created by Xiaodi Wu on 4/15/17.
-//
-
-extension UnsignedInteger {
-  /// Returns the greatest common divisor of `a` and `b`.
-  // @_transparent // @_inlineable
-  public static func gcd(_ a: Self, _ b: Self) -> Self {
-    // An iterative version of Stein's algorithm.
-    if a == 0 { return b } // gcd(0, b) == b
-    if b == 0 { return a } // gcd(a, 0) == a
-
-    var a = a, b = b, shift = 0 as Self
-    while ((a | b) & 1) == 0 {
-      a &>>= 1
-      b &>>= 1
-      shift += 1
-    }
-    // Now, shift is equal to log2(k), where k is the greatest power of 2
-    // dividing a and b.
-    while (a & 1) == 0 { a &>>= 1 } // Now, a is odd.
-    repeat {
-      while (b & 1) == 0 { b &>>= 1 } // Now, b is odd.
-      if a > b { swap(&a, &b) } // Now, a < b.
-      b -= a
-    } while b != 0
-    // Restore common factors of 2.
-    return a &<< shift
-  }
-
-  /// Returns the least common multiple of `a` and `b`.
-  @_transparent // @_inlineable
-  public static func lcm(_ a: Self, _ b: Self) -> Self {
-    if a == 0 || b == 0 { return 0 }
-    return a / .gcd(a, b) * b
-  }
-}
-
-extension UnsignedInteger where Self : FixedWidthInteger {
-  /// Returns the least common multiple of `a` and `b` and a flag to indicate
-  /// whether overflow occurred during the operation.
-  // @_transparent // @_inlineable
-  public static func lcmReportingOverflow(_ a: Self, _ b: Self)
-    -> (partialValue: Self, overflow: ArithmeticOverflow) {
-    if a == 0 || b == 0 { return (partialValue: 0, overflow: .none) }
-    return (a / .gcd(a, b)).multipliedReportingOverflow(by: b)
-  }
-
-  /// Returns the high and low parts of the least common multiple of `a` and `b`
-  /// computed using full-width arithmetic.
-  // @_transparent // @_inlineable
-  public static func lcmFullWidth(_ a: Self, _ b: Self)
-    -> (high: Self, low: Self.Magnitude) {
-    if a == 0 || b == 0 { return (0, 0) }
-    return (a / .gcd(a, b)).multipliedFullWidth(by: b)
-  }
-}
-
-extension BinaryInteger where Magnitude : UnsignedInteger {
-  /// Returns the greatest common divisor of `a` and `b`.
-  @_transparent // @_inlineable
-  public static func gcd(_ a: Self, _ b: Self) -> Self {
-    return Self(Magnitude.gcd(a.magnitude, b.magnitude))
-  }
-
-  /// Returns the least common multiple of `a` and `b`.
-  @_transparent // @_inlineable
-  public static func lcm(_ a: Self, _ b: Self) -> Self {
-    return Self(Magnitude.lcm(a.magnitude, b.magnitude))
-  }
-}
-
-// `BinaryInteger where Self : FixedWidthInteger` may seem superfluous, but it
-// is necessary to disambiguate calls to `Magnitude.lcmReportingOverflow(_:_:)`
-// and `Magnitude.lcmFullWidth(_:_:)`.
-extension BinaryInteger
-where Self : FixedWidthInteger, Magnitude : FixedWidthInteger & UnsignedInteger,
-  Magnitude.Magnitude == Magnitude {
-  /// Returns the greatest common divisor of `a` and `b` and a flag to indicate
-  /// whether overflow occurred during the operation.
-  // @_transparent // @_inlineable
-  public static func gcdReportingOverflow(_ a: Self, _ b: Self)
-    -> (partialValue: Self, overflow: ArithmeticOverflow) {
-    let t = Self(extendingOrTruncating: Magnitude.gcd(a.magnitude, b.magnitude))
-    return (
-      partialValue: t,
-      overflow: ArithmeticOverflow(t < 0)
-    )
-  }
-
-  /// Returns the least common multiple of `a` and `b` and a flag to indicate
-  /// whether overflow occurred during the operation.
-  // @_transparent // @_inlineable
-  public static func lcmReportingOverflow(_ a: Self, _ b: Self)
-    -> (partialValue: Self, overflow: ArithmeticOverflow) {
-    let (t, overflow) = Magnitude.lcmReportingOverflow(a.magnitude, b.magnitude)
-    let u = Self(extendingOrTruncating: t)
-    return (
-      partialValue: u,
-      overflow: ArithmeticOverflow(overflow == .overflow || u < 0)
-    )
-  }
-
-  /// Returns the high and low parts of the least common multiple of `a` and `b`
-  /// computed using full-width arithmetic.
-  // @_transparent // @_inlineable
-  public static func lcmFullWidth(_ a: Self, _ b: Self)
-    -> (high: Self, low: Self.Magnitude) {
-    let t = Magnitude.lcmFullWidth(a.magnitude, b.magnitude)
-    return (high: Self(t.high), low: t.low)
-  }
-}

Sources/Hash.swift

@@ -5,30 +5,30 @@
 //  Created by Xiaodi Wu on 4/1/17.
 //
 
-// @_versioned
-internal func _fnv1a(_ data: AnyHashable...) -> Int {
-  let basis, prime: UInt
-  switch MemoryLayout<UInt>.size {
-  case 4:
-    basis = 2166136261
-    prime = 16777619
-  case 8:
-    basis = 14695981039346656037
-    prime = 1099511628211
-  /*
-  case 16:
-    basis = 144066263297769815596495629667062367629
-    prime = 309485009821345068724781371
-  */
-  default:
-    fatalError("Unsupported UInt bit width for FNV hash")
-  }
+/// A namespace for hash utilities.
+internal enum _Hash { }
 
-  let result = data.reduce(basis) { partialResult, element in
-    var h = element.hashValue
-    return withUnsafeBytes(of: &h) { bytes in
-      bytes.reduce(partialResult) { ($0 ^ UInt($1)) &* prime }
+extension _Hash {
+  // @_versioned
+  /// Returns the result of combining `seed` with the `hashValue` of the given
+  /// `values`.
+  ///
+  /// Combining is performed using [a hash function][ref] described by T.C. Hoad
+  /// and J. Zobel that is also adopted in the Boost C++ libraries.
+  ///
+  /// [ref]: http://goanna.cs.rmit.edu.au/~jz/fulltext/jasist-tch.pdf
+  internal static func _combine(seed: Int = 0, _ values: AnyHashable...) -> Int {
+    // Use a magic number based on the golden ratio
+    // (0x1.9e3779b97f4a7c15f39cc0605cedc8341082276bf3a27251f86c6a11d0c18e95p0).
+#if arch(i386) || arch(arm)
+    let magic = 0x9e3779b9 as UInt
+#else
+    let magic = 0x9e3779b97f4a7c15 as UInt
+#endif
+    var x = UInt(bitPattern: seed)
+    for v in values {
+      x ^= UInt(bitPattern: v.hashValue) &+ magic &+ (x &<< 6) &+ (x &>> 2)
     }
+    return Int(bitPattern: x)
   }
-  return Int(bitPattern: result)
 }

Sources/IntegerAlgorithms.swift

@@ -0,0 +1,422 @@
+//
+//  IntegerAlgorithms.swift
+//  NumericAnnex
+//
+//  Created by Xiaodi Wu on 4/15/17.
+//
+
+extension BinaryInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Exponentiation
+  // ---------------------------------------------------------------------------
+
+  // We need to factor out the implementation of `**` so that the tie-breaking
+  // operators implemented in extensions to concrete integer types can call it.
+  @_versioned
+  internal static func _pow(_ lhs: Self, _ rhs: Self) -> Self {
+    var x = lhs, n = rhs
+    if Self.isSigned && n < 0 {
+      x = 1 / x
+      n = 0 - n
+    } else if n == 0 {
+      return 1
+    }
+    // Exponentiate by iterative squaring.
+    var y = 1 as Self
+    while n > 1 {
+      if n % 2 == 1 {
+        y *= x
+      }
+      x *= x
+      n /= 2
+    }
+    return x * y
+  }
+
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: Self, rhs: Self) -> Self {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout Self, rhs: Self) {
+    lhs = lhs ** rhs
+  }
+
+  /// Returns the result of raising `base` to the power of `exponent`, rounded
+  /// to a representable value (deprecated).
+  @available(*, deprecated, message: "Use operator instead")
+  public static func pow(_ base: Self, _ exponent: Self) -> Self {
+    return base ** exponent
+  }
+
+  // ---------------------------------------------------------------------------
+  // MARK: Square Root
+  // ---------------------------------------------------------------------------
+
+  /// Returns the square root of `x`, rounding toward zero. If `x` is negative,
+  /// a runtime error may occur.
+  @_transparent // @_inlineable
+  public static func sqrt(_ x: Self) -> Self {
+    precondition(!Self.isSigned || x >= 0)
+    var shift = x.bitWidth - 1
+    shift -= shift % 2
+
+    var x = x
+    var result = 0 as Self
+    while shift >= 0 {
+      result *= 2
+      let temporary = 2 * result + 1
+      if temporary <= x >> shift {
+        x -= temporary << shift
+        result += 1
+      }
+      shift -= 2
+    }
+    return result
+  }
+}
+
+extension UnsignedInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Cube Root
+  // ---------------------------------------------------------------------------
+
+  /// Returns the cube root of `x`, rounding toward zero.
+  @_transparent // @_inlineable
+  public static func cbrt(_ x: Self) -> Self {
+    var shift = x.bitWidth - 1
+    shift -= shift % 3
+
+    var x = x
+    var result = 0 as Self
+    while shift >= 0 {
+      result *= 2
+      let temporary = 3 * result * (result + 1) + 1
+      if temporary <= x >> shift {
+        x -= temporary << shift
+        result += 1
+      }
+      shift -= 3
+    }
+    return result
+  }
+
+  // ---------------------------------------------------------------------------
+  // MARK: Factoring
+  // ---------------------------------------------------------------------------
+
+  /// Returns the greatest common divisor of `a` and `b`.
+  @_transparent // @_inlineable
+  public static func gcd(_ a: Self, _ b: Self) -> Self {
+    // An iterative version of Stein's algorithm.
+    if a == 0 { return b } // gcd(0, b) == b
+    if b == 0 { return a } // gcd(a, 0) == a
+
+    var a = a, b = b, shift = 0 as Self
+    while ((a | b) & 1) == 0 {
+      a >>= 1
+      b >>= 1
+      shift += 1
+    }
+    // Now, shift is equal to log2(k), where k is the greatest power of 2
+    // dividing a and b.
+    while (a & 1) == 0 { a >>= 1 } // Now, a is odd.
+    repeat {
+      while (b & 1) == 0 { b >>= 1 } // Now, b is odd.
+      if a > b { swap(&a, &b) } // Now, a < b.
+      b -= a
+    } while b != 0
+    // Restore common factors of 2.
+    return a << shift
+  }
+
+  /// Returns the least common multiple of `a` and `b`.
+  @_transparent // @_inlineable
+  public static func lcm(_ a: Self, _ b: Self) -> Self {
+    if a == 0 || b == 0 { return 0 }
+    return a / .gcd(a, b) * b
+  }
+}
+
+extension UnsignedInteger where Self : FixedWidthInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Factoring (Fixed-Width)
+  // ---------------------------------------------------------------------------
+
+  /// Returns the least common multiple of `a` and `b` and a flag to indicate
+  /// whether overflow occurred during the operation.
+  @_transparent // @_inlineable
+  public static func lcmReportingOverflow(_ a: Self, _ b: Self)
+    -> (partialValue: Self, overflow: Bool) {
+    if a == 0 || b == 0 { return (0, false) }
+    return (a / .gcd(a, b)).multipliedReportingOverflow(by: b)
+  }
+
+  /// Returns the high and low parts of the least common multiple of `a` and `b`
+  /// computed using full-width arithmetic.
+  @_transparent // @_inlineable
+  public static func lcmFullWidth(_ a: Self, _ b: Self)
+    -> (high: Self, low: Self.Magnitude) {
+    if a == 0 || b == 0 { return (0, 0) }
+    return (a / .gcd(a, b)).multipliedFullWidth(by: b)
+  }
+}
+
+extension BinaryInteger where Magnitude : UnsignedInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Cube Root
+  // ---------------------------------------------------------------------------
+
+  /// Returns the cube root of `x`, rounding toward zero.
+  @_transparent // @_inlineable
+  public static func cbrt(_ x: Self) -> Self {
+    return x < 0
+      ? 0 - Self(Magnitude.cbrt(x.magnitude))
+      : Self(Magnitude.cbrt(x.magnitude))
+  }
+
+  // ---------------------------------------------------------------------------
+  // MARK: Factoring
+  // ---------------------------------------------------------------------------
+
+  /// Returns the greatest common divisor of `a` and `b`.
+  @_transparent // @_inlineable
+  public static func gcd(_ a: Self, _ b: Self) -> Self {
+    return Self(Magnitude.gcd(a.magnitude, b.magnitude))
+  }
+
+  /// Returns the least common multiple of `a` and `b`.
+  @_transparent // @_inlineable
+  public static func lcm(_ a: Self, _ b: Self) -> Self {
+    return Self(Magnitude.lcm(a.magnitude, b.magnitude))
+  }
+}
+
+// `BinaryInteger where Self : FixedWidthInteger` may seem superfluous, but it
+// is necessary to disambiguate calls to `Magnitude.lcmReportingOverflow(_:_:)`
+// and `Magnitude.lcmFullWidth(_:_:)`.
+extension BinaryInteger
+where Self : FixedWidthInteger, Magnitude : UnsignedInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Factoring (Fixed-Width)
+  // ---------------------------------------------------------------------------
+
+  /// Returns the greatest common divisor of `a` and `b` and a flag to indicate
+  /// whether overflow occurred during the operation.
+  @_transparent // @_inlineable
+  public static func gcdReportingOverflow(_ a: Self, _ b: Self)
+    -> (partialValue: Self, overflow: Bool) {
+    let t = Self(truncatingIfNeeded: Magnitude.gcd(a.magnitude, b.magnitude))
+    return (t, t < 0)
+  }
+
+  /// Returns the least common multiple of `a` and `b` and a flag to indicate
+  /// whether overflow occurred during the operation.
+  @_transparent // @_inlineable
+  public static func lcmReportingOverflow(_ a: Self, _ b: Self)
+    -> (partialValue: Self, overflow: Bool) {
+    let (t, overflow) = Magnitude.lcmReportingOverflow(a.magnitude, b.magnitude)
+    let u = Self(truncatingIfNeeded: t)
+    return (u, overflow || u < 0)
+  }
+
+  /// Returns the high and low parts of the least common multiple of `a` and `b`
+  /// computed using full-width arithmetic.
+  @_transparent // @_inlineable
+  public static func lcmFullWidth(_ a: Self, _ b: Self)
+    -> (high: Self, low: Self.Magnitude) {
+    let t = Magnitude.lcmFullWidth(a.magnitude, b.magnitude)
+    return (Self(t.high), t.low)
+  }
+}
+
+// =============================================================================
+// TIE-BREAKING OPERATORS
+//
+// The following extensions are required so that an expression such as `2 ** 3`
+// is unambiguous as long as `IntegerLiteralType` is one of the standard library
+// types below.
+// =============================================================================
+
+extension Int {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: Int, rhs: Int) -> Int {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout Int, rhs: Int) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension Int8 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: Int8, rhs: Int8) -> Int8 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout Int8, rhs: Int8) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension Int16 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: Int16, rhs: Int16) -> Int16 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout Int16, rhs: Int16) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension Int32 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: Int32, rhs: Int32) -> Int32 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout Int32, rhs: Int32) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension Int64 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: Int64, rhs: Int64) -> Int64 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout Int64, rhs: Int64) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension UInt {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: UInt, rhs: UInt) -> UInt {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout UInt, rhs: UInt) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension UInt8 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: UInt8, rhs: UInt8) -> UInt8 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout UInt8, rhs: UInt8) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension UInt16 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: UInt16, rhs: UInt16) -> UInt16 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout UInt16, rhs: UInt16) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension UInt32 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: UInt32, rhs: UInt32) -> UInt32 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout UInt32, rhs: UInt32) {
+    lhs = lhs ** rhs
+  }
+}
+
+extension UInt64 {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: UInt64, rhs: UInt64) -> UInt64 {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout UInt64, rhs: UInt64) {
+    lhs = lhs ** rhs
+  }
+}
+
+#if false
+extension DoubleWidth {
+  /// Returns the result of raising `lhs` to the power of `rhs`, rounded to a
+  /// representable value.
+  @_transparent // @_inlineable
+  public static func ** (lhs: DoubleWidth, rhs: DoubleWidth) -> DoubleWidth {
+    return _pow(lhs, rhs)
+  }
+
+  /// Raises `lhs` to the power of `rhs` and stores the result in `lhs`, rounded
+  /// to a representable value.
+  @_transparent // @_inlineable
+  public static func **= (lhs: inout DoubleWidth, rhs: DoubleWidth) {
+    lhs = lhs ** rhs
+  }
+}
+#endif

Sources/Math.swift

@@ -5,14 +5,13 @@
 //  Created by Xiaodi Wu on 3/31/17.
 //
 
-/// A signed numeric type that provides square root, cube root, and elementary
-/// transcendental functions.
+/// A signed numeric type that supports elementary functions.
 ///
 /// The `Math` protocol provides a suitable basis for writing functions that
-/// work on any real or complex floating-point type which provides the required
+/// work on any real or complex floating-point type that supports the required
 /// functions.
 public protocol Math : SignedNumeric {
-  /// The mathematical constant pi (_π_).
+  /// The mathematical constant pi (π).
   ///
   /// This value should be rounded toward zero to keep user computations with
   /// angles from inadvertently ending up in the wrong quadrant. A type that
@@ -23,7 +22,7 @@ public protocol Math : SignedNumeric {
   /// The mathematical constant _e_, or Euler's number.
   static var e: Self { get }
 
-  /// The mathematical constant phi (_φ_), or golden ratio.
+  /// The mathematical constant phi (φ), or golden ratio.
   static var phi: Self { get }
 
   /// Returns the quotient obtained by dividing the first value by the second,
@@ -42,6 +41,22 @@ public protocol Math : SignedNumeric {
   ///   - rhs: The value by which to divide `lhs`.
   static func /= (lhs: inout Self, rhs: Self)
 
+  /// Returns the result of raising the first value to the power of the second,
+  /// rounded to a representable value.
+  ///
+  /// - Parameters:
+  ///   - lhs: The value to be raised to the power of `rhs`.
+  ///   - rhs: The value by which to raise `lhs`.
+  static func ** (lhs: Self, rhs: Self) -> Self
+
+  /// Raises the left-hand side to the power of the right-hand side and stores
+  /// the result in the left-hand side, rounded to a representable value.
+  ///
+  /// - Parameters:
+  ///   - lhs: The value to be raised to the power of `rhs`.
+  ///   - rhs: The value by which to raise `lhs`.
+  static func **= (lhs: inout Self, rhs: Self)
+
   /// Returns the natural exponential of the value, rounded to a representable
   /// value.
   ///
@@ -113,17 +128,6 @@ public protocol Math : SignedNumeric {
   /// - SeeAlso: `cbrt(_:)`
   func cubeRoot() -> Self
 
-  /// Returns the result of raising `base` to the power of the value, rounded to
-  /// a representable value.
-  ///
-  /// - Note: The argument is the _base_ and the receiver is the _exponent_.
-  ///
-  /// - Parameters:
-  ///   - base: The base.
-  ///
-  /// - SeeAlso: `pow(_:_:)`
-  func power(of base: Self) -> Self
-
   /// Returns the sine of the value (given in radians), rounded to a
   /// representable value.
   ///
@@ -198,16 +202,45 @@ public protocol Math : SignedNumeric {
 }
 
 extension Math {
-  // @_transparent
+  /// The mathematical constant _e_, or Euler's number (default implementation).
   public static var e: Self {
     return Self.exp(1 as Self)
   }
 
-  // @_transparent
+  /// The mathematical constant phi (φ), or golden ratio (default
+  /// implementation).
   public static var phi: Self {
-    return Self.sqrt(((1 as Self) + Self.sqrt(5 as Self)) / (2 as Self))
+    return ((1 as Self) + Self.sqrt(5 as Self)) / (2 as Self)
   }
 
+#if false
+  /// Returns the quotient obtained by dividing the first value by the second,
+  /// rounded to a representable value (default implementation).
+  ///
+  /// - Parameters:
+  ///   - lhs: The value to divide.
+  ///   - rhs: The value by which to divide `lhs`.
+  @_transparent
+  public static func / (lhs: Self, rhs: Self) -> Self {
+    var lhs = lhs
+    lhs /= rhs
+    return lhs
+  }
+
+  /// Returns the result of raising the first value to the power of the second,
+  /// rounded to a representable value (default implementation).
+  ///
+  /// - Parameters:
+  ///   - lhs: The value to be raised to the power of `rhs`.
+  ///   - rhs: The value by which to raise `lhs`.
+  @_transparent
+  public static func ** (lhs: Self, rhs: Self) -> Self {
+    var lhs = lhs
+    lhs **= rhs
+    return lhs
+  }
+#endif
+
   public func binaryExponential() -> Self {
     return Self.exp(self * Self.log(2 as Self))
   }
@@ -232,10 +265,6 @@ extension Math {
     return Self.log(self + (1 as Self))
   }
 
-  public func cubeRoot() -> Self {
-    return Self.pow(self, 1 / 3 as Self)
-  }
-
   public func tangent() -> Self {
     return sine() / cosine()
   }
@@ -246,6 +275,17 @@ extension Math {
 }
 
 extension Math {
+  /// Returns the result of raising `base` to the power of `exponent`, rounded
+  /// to a representable value (deprecated).
+  ///
+  /// - Parameters:
+  ///   - base: The base to be raised to the power of `exponent`.
+  ///   - exponent: The exponent by which to raise `base`.
+  @available(*, deprecated, message: "Use operator instead")
+  public static func pow(_ base: Self, _ exponent: Self) -> Self {
+    return base ** exponent
+  }
+
   /// Returns the natural exponential of `x`, rounded to a representable value.
   ///
   /// The natural exponential of a value `x` is _e_ (2.7182818...) raised to the
@@ -322,7 +362,7 @@ extension Math {
   /// Returns the natural (base _e_) logarithm of the result of adding one to
   /// `x`, rounded to a representable value.
   ///
-  /// - SeeAlso: `naturalLogarithmOfOnePlus()`
+  /// - SeeAlso: `naturalLogarithmOnePlus()`
   @_transparent // @_inlineable
   public static func log1p(_ x: Self) -> Self {
     return x.naturalLogarithmOnePlus()
@@ -344,15 +384,6 @@ extension Math {
     return x.cubeRoot()
   }
 
-  /// Returns the result of raising `base` to the power of `exponent`, rounded
-  /// to a representable value.
-  ///
-  /// - SeeAlso: `power(of:)`
-  @_transparent // @_inlineable
-  public static func pow(_ base: Self, _ exponent: Self) -> Self {
-    return exponent.power(of: base)
-  }
-
   /// Returns the sine of `x` (given in radians), rounded to a representable
   /// value.
   ///

Sources/PRNG.swift

@@ -11,23 +11,51 @@ import Glibc
 import Security
 #endif
 
-/// A pseudo-random number generator.
+/// A pseudo-random number generator (PRNG).
+///
+/// Reference types that conform to `PRNG` are infinite sequences of
+/// pseudo-random elements. Protocol extension methods iterate over such a
+/// sequence as necessary to generate pseudo-random values from the desired
+/// distribution.
+///
+/// Considerations for Conforming Types
+/// -----------------------------------
+///
+/// For clarity to end users, custom PRNGs may be implemented in an extension to
+/// `Random`. For instance, the `xoroshiro128+` algorithm is implemented in a
+/// final class named `Random.Xoroshiro`.
+///
+/// The static methods `_entropy(_:)` and `_entropy(_:count:)` return
+/// cryptographically secure random bytes that may be useful for seeding your
+/// custom PRNG. However, these methods may return `nil` if the requested number
+/// of random bytes is not available, and they are not recommended as a routine
+/// source of random data.
+///
+/// Adding Other Probability Distributions
+/// --------------------------------------
+///
+/// Many built-in protocol extension methods make use of the primitive,
+/// overloaded method `_random(_:bitCount:)`. You may wish to use the same
+/// method in new protocol extension methods that return pseudo-random
+/// values from other probability distributions.
+///
+/// The method `_random(_:bitCount:)` generates uniformly distributed binary
+/// floating-point values in the half-open range [0, 1) with a precision of
+/// either `bitCount` or the significand bit count of the floating-point type,
+/// whichever is less. Additionally, this method generates uniformly distributed
+/// unsigned integers in the half-open range [0, 2 ** _x_), where ** is the
+/// exponentiation operator and _x_ is the lesser of `bitCount` and the bit
+/// width of the integer type.
+///
+/// For end users, however, the recommended spelling for a uniformly distributed
+/// numeric value is `uniform()`; that method is overloaded to permit custom
+/// minimum and maximum values for the uniform distribution.
 public protocol PRNG : class, IteratorProtocol, Sequence
-where Element : FixedWidthInteger & UnsignedInteger, SubSequence : Sequence,
-  /* Element == Iterator.Element, Element == SubSequence.Iterator.Element, */
-  Element == SubSequence.Element {
+where Element : FixedWidthInteger & UnsignedInteger {
   /// A type that can represent the internal state of the pseudo-random number
   /// generator.
   associatedtype State
 
-  /// The maximum value that may be generated by the pseudo-random number
-  /// generator.
-  static var max: Element { get }
-
-  /// The minimum value that may be generated by the pseudo-random number
-  /// generator.
-  static var min: Element { get }
-
   /// The internal state of the pseudo-random number generator.
   var state: State { get set }
 
@@ -43,11 +71,23 @@ where Element : FixedWidthInteger & UnsignedInteger, SubSequence : Sequence,
   /// If cryptographically secure random bytes are unavailable, the result is
   /// `nil`.
   init?()
+  
+  /// The maximum value that may be generated by the pseudo-random number
+  /// generator.
+  static var max: Element { get }
+
+  /// The minimum value that may be generated by the pseudo-random number
+  /// generator.
+  static var min: Element { get }
 }
 
 extension PRNG {
+  /// The maximum value that may be generated by the pseudo-random number
+  /// generator (default implementation: `Element.max`).
   public static var max: Element { return Element.max }
 
+  /// The minimum value that may be generated by the pseudo-random number
+  /// generator (default implementation: `Element.min`).
   public static var min: Element { return Element.min }
 
   /// The number of pseudo-random bits available from a value generated by the
@@ -76,8 +116,7 @@ extension PRNG {
 #else
     // Sandboxing can make `urandom` unavailable.
     let result = withUnsafeMutableBytes(of: &value) { ptr -> Int32 in
-      let bytes = ptr.baseAddress!.bindMemory(to: UInt8.self, capacity: size)
-      defer { ptr.baseAddress!.bindMemory(to: T.self, capacity: 1) }
+      let bytes = ptr.baseAddress!.assumingMemoryBound(to: UInt8.self)
       return SecRandomCopyBytes(nil, size, bytes)
     }
     guard result == errSecSuccess else { return nil }
@@ -100,10 +139,8 @@ extension PRNG {
     guard read == count else { return nil }
 #else
     let result = value.withUnsafeMutableBytes { ptr -> Int32 in
-      let n = stride * count
-      let bytes = ptr.baseAddress!.bindMemory(to: UInt8.self, capacity: n)
-      defer { ptr.baseAddress!.bindMemory(to: T.self, capacity: count) }
-      return SecRandomCopyBytes(nil, n, bytes)
+      let bytes = ptr.baseAddress!.assumingMemoryBound(to: UInt8.self)
+      return SecRandomCopyBytes(nil, stride * count, bytes)
     }
     guard result == errSecSuccess else { return nil }
 #endif
@@ -125,7 +162,7 @@ extension PRNG {
       bitCount == T.bitWidth {
       // It is an awkward way of spelling `next()`, but it is necessary.
       guard let next = first(where: { _ in true }) else { fatalError() }
-      return T(extendingOrTruncating: next)
+      return T(truncatingIfNeeded: next)
     }
 
     let (quotient, remainder) =
@@ -135,15 +172,14 @@ extension PRNG {
     // Call `next()` at least `quotient` times.
     for i in 0..<quotient {
       guard let next = first(where: { $0 <= max }) else { fatalError() }
-      temporary +=
-        T(extendingOrTruncating: next) &<< (randomBitWidth * i)
+      temporary += T(truncatingIfNeeded: next) &<< (randomBitWidth * i)
     }
     // If `remainder != 0`, call `next()` at least one more time.
     if remainder != 0 {
       guard let next = first(where: { $0 <= max }) else { fatalError() }
       let mask = Element.max &>> (Element.bitWidth - remainder)
       temporary +=
-        T(extendingOrTruncating: next & mask) &<< (randomBitWidth * quotient)
+        T(truncatingIfNeeded: next & mask) &<< (randomBitWidth * quotient)
     }
     return temporary
   }
@@ -209,7 +245,7 @@ extension PRNG {
   /// through `b` (inclusive) from the discrete uniform distribution.
   public func uniform<T : FixedWidthInteger & SignedInteger>(
     _: T.Type = T.self, a: T, b: T
-  ) -> T where T.Magnitude : FixedWidthInteger & UnsignedInteger {
+  ) -> T where T.Magnitude : UnsignedInteger {
     precondition(
       b >= a,
       "Discrete uniform distribution parameter b should not be less than a"
@@ -237,7 +273,7 @@ extension PRNG {
   @_transparent // @_inlineable
   public func uniform<T : FixedWidthInteger & SignedInteger>(
     _: T.Type = T.self
-  ) -> T where T.Magnitude : FixedWidthInteger & UnsignedInteger {
+  ) -> T where T.Magnitude : UnsignedInteger {
     return uniform(a: T.min, b: T.max)
   }
 
@@ -248,7 +284,7 @@ extension PRNG {
   public func uniform<T : FixedWidthInteger & SignedInteger>(
     _: T.Type = T.self, a: T, b: T, count: Int
   ) -> UnfoldSequence<T, Int>
-  where T.Magnitude : FixedWidthInteger & UnsignedInteger {
+  where T.Magnitude : UnsignedInteger {
     precondition(count >= 0, "Element count should be non-negative")
     return sequence(state: 0) { (state: inout Int) -> T? in
       state += 1
@@ -263,13 +299,14 @@ extension PRNG {
   public func uniform<T : FixedWidthInteger & SignedInteger>(
     _: T.Type = T.self, count: Int
   ) -> UnfoldSequence<T, Int>
-  where T.Magnitude : FixedWidthInteger & UnsignedInteger {
+  where T.Magnitude : UnsignedInteger {
     return uniform(a: T.min, b: T.max, count: count)
   }
 }
 
 // FIXME: If `FloatingPoint.init(_: FixedWidthInteger)` is added
 // then it becomes possible to remove the constraint `Element == UInt64`.
+
 extension PRNG where Element == UInt64 {
   /// Generates a pseudo-random binary floating-point value of type `T` in the
   /// range from 0 to 1 (exclusive) with `min(bitCount, T.significandBitCount)`
@@ -342,7 +379,7 @@ extension PRNG where Element == UInt64 {
 
 #if false
 extension PRNG where Element == UInt64 {
-  public func bernoulli<T : BinaryFloatingPoint>(
+  /* public */ func bernoulli<T : BinaryFloatingPoint>(
     _: Bool.Type = Bool.self, p: T
   ) -> Bool {
     precondition(
@@ -357,12 +394,12 @@ extension PRNG where Element == UInt64 {
   }
 
   @_transparent // @_inlineable
-  public func bernoulli(_: Bool.Type = Bool.self) -> Bool {
+  /* public */ func bernoulli(_: Bool.Type = Bool.self) -> Bool {
     return bernoulli(p: 0.5)
   }
 
   @_transparent // @_inlineable
-  public func bernoulli<T : BinaryFloatingPoint>(
+  /* public */ func bernoulli<T : BinaryFloatingPoint>(
     _: Bool.Type = Bool.self, p: T, count: Int
   ) -> UnfoldSequence<Bool, Int> {
     precondition(count >= 0, "Element count should be non-negative")
@@ -373,13 +410,13 @@ extension PRNG where Element == UInt64 {
   }
 
   @_transparent // @_inlineable
-  public func bernoulli(
+  /* public */ func bernoulli(
     _: Bool.Type = Bool.self, count: Int
   ) -> UnfoldSequence<Bool, Int> {
     return bernoulli(p: 0.5, count: count)
   }
 
-  public func exponential<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func exponential<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self, lambda: T
   ) -> T {
     precondition(
@@ -394,14 +431,14 @@ extension PRNG where Element == UInt64 {
   }
 
   @_transparent // @_inlineable
-  public func exponential<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func exponential<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self
   ) -> T {
     return exponential(lambda: 1)
   }
 
   @_transparent // @_inlineable
-  public func exponential<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func exponential<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self, lambda: T, count: Int
   ) -> UnfoldSequence<T, Int> {
     precondition(count >= 0, "Element count should be non-negative")
@@ -412,13 +449,13 @@ extension PRNG where Element == UInt64 {
   }
 
   @_transparent // @_inlineable
-  public func exponential<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func exponential<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self, count: Int
   ) -> UnfoldSequence<T, Int> {
     return exponential(lambda: 1, count: count)
   }
   
-  public func weibull<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func weibull<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self, lambda: T, kappa: T
   ) -> T {
     precondition(
@@ -433,14 +470,14 @@ extension PRNG where Element == UInt64 {
   }
 
   @_transparent // @_inlineable
-  public func weibull<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func weibull<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self
   ) -> T {
     return weibull(lambda: 1, kappa: 1)
   }
 
   @_transparent // @_inlineable
-  public func weibull<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func weibull<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self, lambda: T, kappa: T, count: Int
   ) -> UnfoldSequence<T, Int> {
     precondition(count >= 0, "Element count should be non-negative")
@@ -451,7 +488,7 @@ extension PRNG where Element == UInt64 {
   }
 
   @_transparent // @_inlineable
-  public func weibull<T : BinaryFloatingPoint & FloatingPointMath>(
+  /* public */ func weibull<T : BinaryFloatingPoint & Real>(
     _: T.Type = T.self, count: Int
   ) -> UnfoldSequence<T, Int> {
     return weibull(lambda: 1, kappa: 1, count: count)

Sources/Random.Xoroshiro.swift

@@ -6,14 +6,69 @@
 //
 
 extension Random {
-  /// A pseudo-random number generator that implements [`xoroshiro128+`], a
+  /// A pseudo-random number generator (PRNG) that implements `xoroshiro128+`, a
   /// successor to `xorshift128+` devised by S. Vigna and D. Blackman.
+  ///
+  /// - Warning: Once seeded from cryptographically secure random bytes, an
+  ///   instance of `Random.Xoroshiro` generates high-quality random numbers but
+  ///   is _not_ a cryptographically secure PRNG.
+  ///
+  /// To generate random numbers, create your own instance of `Random.Xoroshiro`
+  /// with an internal state seeded from cryptographically secure random bytes:
+  ///
+  /// ```swift
+  /// let random = Random.Xoroshiro()!
+  /// let x = random.uniform() as Int
+  ///
+  /// // You can also pass the desired result type as an argument.
+  /// let y = random.uniform(Int.self)
+  ///
+  /// if x > y {
+  ///   print("Here's a random value between 0 and 42 (inclusive):")
+  ///   print(random.uniform(a: 0, b: 42))
+  /// } else {
+  ///   print("Here's a random value between -42 and 0 (inclusive):")
+  ///   print(random.uniform(a: -42, b: 0))
+  /// }
+  /// ```
+  ///
+  /// - SeeAlso: `Random`, `PRNG`
   public final class Xoroshiro : PRNG {
+    /// The internal state of the pseudo-random number generator.
     public var state: (UInt64, UInt64)
 
-    // TODO: Document this function.
+    /// Creates a pseudo-random number generator with the given internal state.
+    ///
+    /// - Parameters:
+    ///   - state: The value to be used as the generator's internal state.
+    public init(state: (UInt64, UInt64)) {
+      self.state = state
+    }
+
+    /// Creates a pseudo-random number generator with an internal state seeded
+    /// using cryptographically secure random bytes.
+    ///
+    /// If cryptographically secure random bytes are unavailable, the result is
+    /// `nil`.
+    public convenience init?() {
+      self.init(_entropy: Xoroshiro._entropy(UInt64.self, count: 2))
+    }
+
+    /// Creates a pseudo-random number generator with an internal state seeded
+    /// using cryptographically secure random bytes.
+    ///
+    /// If cryptographically secure random bytes are unavailable, the result is
+    /// `nil`.
+    internal init?(_entropy: @autoclosure () -> [UInt64]?) {
+      repeat {
+        guard let entropy = _entropy() else { return nil }
+        self.state = (entropy[0], entropy[1])
+      } while self.state == (0, 0)
+    }
+
     // @_versioned
     internal static func _rotl(_ value: UInt64, _ count: Int) -> UInt64 {
+      // TODO: Document this function.
       return (value &<< count) | (value &>> (64 &- count))
     }
 
@@ -26,18 +81,5 @@ extension Random {
       state.1 = Xoroshiro._rotl(y, 36)
       return result
     }
-
-    public init(state: (UInt64, UInt64)) {
-      self.state = state
-    }
-
-    public init?() {
-      repeat {
-        guard let entropy = Xoroshiro._entropy(UInt64.self, count: 2) else {
-          return nil
-        }
-        self.state = (entropy[0], entropy[1])
-      } while self.state == (0, 0)
-    }
   }
 }

Sources/Random.swift

@@ -5,13 +5,68 @@
 //  Created by Xiaodi Wu on 5/13/17.
 //
 
-/// A pseudo-random number generator that implements [`xorshift128+`][ref], an
-/// efficient linear-feedback shift register.
+/// A pseudo-random number generator (PRNG) that implements
+/// [`xorshift128+`][ref], an efficient linear-feedback shift register.
+///
+/// - Warning: Once seeded from cryptographically secure random bytes, an
+///   instance of `Random` generates high-quality random numbers but is _not_ a
+///   cryptographically secure PRNG.
+///
+/// To generate random numbers, create your own instance of `Random` with an
+/// internal state seeded from cryptographically secure random bytes:
+///
+/// ```swift
+/// let random = Random()!
+/// let x = random.uniform() as Int
+///
+/// // You can also pass the desired result type as an argument.
+/// let y = random.uniform(Int.self)
+///
+/// if x > y {
+///   print("Here's a random value between 0 and 42 (inclusive):")
+///   print(random.uniform(a: 0, b: 42))
+/// } else {
+///   print("Here's a random value between -42 and 0 (inclusive):")
+///   print(random.uniform(a: -42, b: 0))
+/// }
+/// ```
+///
+/// - SeeAlso: `Random.Xoroshiro`, `PRNG`
 ///
 /// [ref]: http://vigna.di.unimi.it/ftp/papers/xorshiftplus.pdf
 public final class Random : PRNG {
+  /// The internal state of the pseudo-random number generator.
   public var state: (UInt64, UInt64)
 
+  /// Creates a pseudo-random number generator with the given internal state.
+  ///
+  /// - Parameters:
+  ///   - state: The value to be used as the generator's internal state.
+  public init(state: (UInt64, UInt64)) {
+    self.state = state
+  }
+
+  /// Creates a pseudo-random number generator with an internal state seeded
+  /// using cryptographically secure random bytes.
+  ///
+  /// If cryptographically secure random bytes are unavailable, the result is
+  /// `nil`.
+  public convenience init?() {
+    self.init(_entropy: Random._entropy(UInt64.self, count: 2))
+  }
+
+  /// Creates a pseudo-random number generator with an internal state seeded
+  /// using cryptographically secure random bytes.
+  ///
+  /// If cryptographically secure random bytes are unavailable, the result is
+  /// `nil`.
+  internal init?(_entropy: @autoclosure () -> [UInt64]?) {
+    repeat {
+      guard let entropy = _entropy() else { return nil }
+      self.state = (entropy[0], entropy[1])
+    } while self.state == (0, 0)
+  }
+  
   public func next() -> UInt64? {
 #if true
     // An updated version of xorshift128+.
@@ -35,17 +90,4 @@ public final class Random : PRNG {
     return state.1 &+ y
 #endif
   }
-
-  public init(state: (UInt64, UInt64)) {
-    self.state = state
-  }
-
-  public init?() {
-    repeat {
-      guard let entropy = Random._entropy(UInt64.self, count: 2) else {
-        return nil
-      }
-      self.state = (entropy[0], entropy[1])
-    } while self.state == (0, 0)
-  }
 }

Sources/Rational+SignedNumeric.swift

@@ -1,268 +0,0 @@
-//
-//  Rational+SignedNumeric.swift
-//  NumericAnnex
-//
-//  Created by Xiaodi Wu on 4/15/17.
-//
-
-extension Rational : Numeric {
-  // @_transparent // @_inlineable
-  public init?<U>(exactly source: U) where U : BinaryInteger {
-    guard let t = T(exactly: source) else { return nil }
-    // Ensure that `t.magnitude` is representable as a `T`.
-    guard let _ = T(exactly: t.magnitude) else { return nil }
-    self.numerator = t
-    self.denominator = 1
-  }
-
-  // @_transparent // @_inlineable
-  public static func + (lhs: Rational, rhs: Rational) -> Rational {
-    if lhs.denominator == 0 {
-      if rhs.denominator != 0 || lhs.numerator == 0 { return lhs }
-      if lhs.numerator > 0 { return rhs.numerator < 0 ? .nan : rhs }
-      return rhs.numerator > 0 ? .nan : rhs
-    }
-    if rhs.denominator == 0 { return rhs }
-
-    let ldm = lhs.denominator.magnitude
-    let rdm = rhs.denominator.magnitude
-    let gcd = T.Magnitude.gcd(ldm, rdm)
-    let a = T(rdm / gcd * lhs.numerator.magnitude)
-    let b = T(ldm / gcd * rhs.numerator.magnitude)
-    let n = lhs.sign == .plus
-      ? (rhs.sign == .plus ? a + b : a - b)
-      : (rhs.sign == .plus ? b - a : -a - b)
-    let d = T(ldm / gcd * rdm)
-    return Rational(numerator: n, denominator: d).canonical
-  }
-
-  @_transparent // @_inlineable
-  public static func += (lhs: inout Rational, rhs: Rational) {
-    lhs = lhs + rhs
-  }
-
-  @_transparent // @_inlineable
-  public static func - (lhs: Rational, rhs: Rational) -> Rational {
-    return lhs + (-rhs)
-  }
-
-  @_transparent // @_inlineable
-  public static func -= (lhs: inout Rational, rhs: Rational) {
-    lhs = lhs + (-rhs)
-  }
-  
-  // @_transparent // @_inlineable
-  public static func * (lhs: Rational, rhs: Rational) -> Rational {
-    if lhs.denominator == 0 {
-      if rhs.numerator == 0 { return .nan }
-      return rhs.sign == .plus ? lhs : -lhs
-    }
-    if rhs.denominator == 0 {
-      if lhs.numerator == 0 { return .nan }
-      return lhs.sign == .plus ? rhs : -rhs
-    }
-    
-    let lnm = lhs.numerator.magnitude, ldm = lhs.denominator.magnitude
-    let rnm = rhs.numerator.magnitude, rdm = rhs.denominator.magnitude
-    // Note that if `T` is a signed fixed-width integer type, `gcd(lnm, rdm)` or
-    // `gcd(rnm, ldm)` could be equal to `-T.min`, which is not representable as
-    // a `T`. This is why the following arithmetic is performed with values of
-    // type `T.Magnitude`.
-    let a = T.Magnitude.gcd(lnm, rdm)
-    let b = T.Magnitude.gcd(rnm, ldm)
-    let n = lhs.sign == rhs.sign
-      ? T(lnm / a * (rnm / b))
-      : -T(lnm / a * (rnm / b))
-    let d = T(ldm / b * (rdm / a))
-    return Rational(numerator: n, denominator: d)
-  }
-
-  @_transparent // @_inlineable
-  public static func *= (lhs: inout Rational, rhs: Rational) {
-    lhs = lhs * rhs
-  }
-}
-
-extension BinaryInteger {
-  /// Creates a new binary integer from the given rational value, if it can be
-  /// represented exactly.
-  ///
-  /// If `source` is not representable exactly, the result is `nil`.
-  ///
-  /// - Parameters:
-  ///   - source: A rational value to convert to a binary integer.
-  @_transparent // @_inlineable
-  public init?<U>(exactly source: Rational<U>) {
-    let (whole, fraction) = source.mixed
-    guard fraction.isZero, let exact = Self(exactly: whole) else { return nil }
-    self = exact
-  }
-
-  /// Creates a new binary integer from the given rational value, truncating any
-  /// fractional part.
-  ///
-  /// If `source` is outside the bounds of this type after truncation, a runtime
-  /// error may occur.
-  ///
-  /// - Parameters:
-  ///   - source: A rational value to convert to a binary integer.
-  @_transparent // @_inlineable
-  public init<U>(_ source: Rational<U>) {
-    self = Self(source.mixed.whole)
-  }
-}
-
-extension Rational : SignedNumeric {
-  @_transparent // @_inlineable
-  public static prefix func - (operand: Rational) -> Rational {
-    return Rational(
-      numerator: -operand.numerator, denominator: operand.denominator
-    )
-  }
-
-  @_transparent // @_inlineable
-  public mutating func negate() {
-    numerator.negate()
-  }
-}
-
-extension Rational {
-  /// Returns the quotient obtained by dividing the first value by the second,
-  /// trapping in case of arithmetic overflow.
-  ///
-  /// - Parameters:
-  ///   - lhs: The value to divide.
-  ///   - rhs: The value by which to divide `lhs`.
-  @_transparent // @_inlineable
-  public static func / (lhs: Rational, rhs: Rational) -> Rational {
-    return lhs * rhs.reciprocal()
-  }
-
-  /// Divides the left-hand side by the right-hand side and stores the quotient
-  /// in the left-hand side, trapping in case of arithmetic overflow.
-  ///
-  /// - Parameters:
-  ///   - lhs: The value to divide.
-  ///   - rhs: The value by which to divide `lhs`.
-  @_transparent // @_inlineable
-  public static func /= (lhs: inout Rational, rhs: Rational) {
-    lhs = lhs * rhs.reciprocal()
-  }
-
-  /// Returns this value rounded to an integral value using the specified
-  /// rounding rule.
-  ///
-  /// ```swift
-  /// let x = 7 / 2 as Rational<Int>
-  /// print(x.rounded()) // Prints "4"
-  /// print(x.rounded(.towardZero)) // Prints "3"
-  /// print(x.rounded(.up)) // Prints "4"
-  /// print(x.rounded(.down)) // Prints "3"
-  /// ```
-  ///
-  /// See the `FloatingPointRoundingRule` enumeration for more information about
-  /// the available rounding rules.
-  ///
-  /// - Parameters:
-  ///   - rule: The rounding rule to use.
-  ///
-  /// - SeeAlso: `round(_:)`, `FloatingPointRoundingRule`
-  @_transparent // @_inlineable
-  public func rounded(
-    _ rule: RoundingRule = .toNearestOrAwayFromZero
-  ) -> Rational {
-    var t = self
-    t.round(rule)
-    return t
-  }
-
-  /// Rounds the value to an integral value using the specified rounding rule.
-  ///
-  /// ```swift
-  /// var x = 7 / 2 as Rational<Int>
-  /// x.round() // x == 4
-  ///
-  /// var x = 7 / 2 as Rational<Int>
-  /// x.round(.towardZero) // x == 3
-  ///
-  /// var x = 7 / 2 as Rational<Int>
-  /// x.round(.up) // x == 4
-  ///
-  /// var x = 7 / 2 as Rational<Int>
-  /// x.round(.down) // x == 3
-  /// ```
-  ///
-  /// See the `FloatingPointRoundingRule` enumeration for more information about
-  /// the available rounding rules.
-  ///
-  /// - Parameters:
-  ///   - rule: The rounding rule to use.
-  ///
-  /// - SeeAlso: `round(_:)`, `FloatingPointRoundingRule`
-  @_transparent // @_inlineable
-  public mutating func round(_ rule: RoundingRule = .toNearestOrAwayFromZero) {
-    if denominator == 0 { return }
-
-    let f: T
-    (numerator, f) = numerator.quotientAndRemainder(dividingBy: denominator)
-    // Rounding rules only come into play if the fractional part is non-zero.
-    if f != 0 {
-      switch rule {
-      case .toNearestOrAwayFromZero:
-        fallthrough
-      case .toNearestOrEven:
-        switch denominator.magnitude.quotientAndRemainder(
-          dividingBy: f.magnitude
-        ) {
-        case (2, 0): // Tie.
-          if rule == .toNearestOrEven && numerator % 2 == 0 { break }
-          fallthrough
-        case (1, _): // Nearest is away from zero.
-          if f > 0 { numerator += 1 } else { numerator -= 1 }
-        default: // Nearest is toward zero.
-          break
-        }
-      case .up:
-        if f > 0 { numerator += 1 }
-      case .down:
-        if f < 0 { numerator -= 1 }
-      case .towardZero:
-        break
-      case .awayFromZero:
-        if f > 0 { numerator += 1 } else { numerator -= 1 }
-      }
-    }
-    denominator = 1
-  }
-}
-
-/// Returns the absolute value (magnitude) of `x`.
-@_transparent
-public func abs<T>(_ x: Rational<T>) -> Rational<T> {
-  return x.magnitude
-}
-
-/// Returns the closest integral value greater than or equal to `x`.
-@_transparent
-public func ceil<T>(_ x: Rational<T>) -> Rational<T> {
-  return x.rounded(.up)
-}
-
-/// Returns the closest integral value less than or equal to `x`.
-@_transparent
-public func floor<T>(_ x: Rational<T>) -> Rational<T> {
-  return x.rounded(.down)
-}
-
-/// Returns the closest integral value; if two values are equally close, returns
-/// the one with greater magnitude.
-@_transparent
-public func round<T>(_ x: Rational<T>) -> Rational<T> {
-  return x.rounded()
-}
-
-/// Returns the closest integral value with magnitude less than or equal to `x`.
-@_transparent
-public func trunc<T>(_ x: Rational<T>) -> Rational<T> {
-  return x.rounded(.towardZero)
-}

Sources/Rational.swift

@@ -28,17 +28,16 @@
 /// print(a) // Prints "3/3"
 /// ```
 ///
-/// All arithmetic operations with values in canonical form (i.e. reduced to
-/// lowest terms) return results in canonical form. However, operations with
+/// All arithmetic operations on values in canonical form (i.e. reduced to
+/// lowest terms) return results in canonical form. However, operations on
 /// values not in canonical form may or may not return results that are
 /// themselves in canonical form. The property `canonicalized` is the canonical
 /// form of any value.
 ///
 /// Additional Considerations
-/// =========================
+/// -------------------------
 ///
-/// Special Values
-/// --------------
+/// ### Special Values
 ///
 /// `Rational<T>` does not prohibit zero as a denominator. Any instance with a
 /// positive numerator and zero denominator represents (positive) infinity; any
@@ -56,8 +55,7 @@
 /// a value is NaN. `Rational<T>` arithmetic operations are intended to
 /// propagate NaN in the same manner as analogous floating-point operations.
 ///
-/// Fixed-Width Binary Parts
-/// ------------------------
+/// ### Numerical Limits
 ///
 /// When a value of type `Rational<T>` is in canonical form, the sign of the
 /// numerator is the sign of the value; that is, in canonical form, the sign of
@@ -70,15 +68,23 @@
 /// occurs when the division (`/`) operator is used to create a value of type
 /// `Rational<T>` with numerator `T.min`.
 @_fixed_layout
-public struct Rational<T : SignedInteger>
-where T : _ExpressibleByBuiltinIntegerLiteral, T.Magnitude : UnsignedInteger,
-  T.Magnitude.Magnitude == T.Magnitude {
+public struct Rational<T : SignedInteger> : Codable
+where T : Codable & _ExpressibleByBuiltinIntegerLiteral,
+  T.Magnitude : UnsignedInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Stored Properties
+  // ---------------------------------------------------------------------------
+
   /// The numerator of the rational value.
   public var numerator: T
 
   /// The denominator of the rational value.
   public var denominator: T
 
+  // ---------------------------------------------------------------------------
+  // MARK: Initializers
+  // ---------------------------------------------------------------------------
+
   /// Creates a new value from the given numerator and denominator without
   /// computing its canonical form (i.e., without reducing to lowest terms).
   ///
@@ -99,6 +105,118 @@ where T : _ExpressibleByBuiltinIntegerLiteral, T.Magnitude : UnsignedInteger,
     self.denominator = denominator
   }
 
+  /// Creates a new rational value from the given binary integer.
+  ///
+  /// If `source` or its magnitude is not representable as a numerator of type
+  /// `T`, a runtime error may occur.
+  ///
+  /// - Parameters:
+  ///   - source: A binary integer to convert to a rational value.
+  @_transparent // @_inlineable
+  public init<Source : BinaryInteger>(_ source: Source) {
+    let t = T(source)
+    // Ensure that `t.magnitude` is representable as a `T`.
+    _ = T(t.magnitude)
+    self.numerator = t
+    self.denominator = 1
+  }
+
+  /// Creates a new rational value from the given binary floating-point value.
+  ///
+  /// If `source` or its magnitude is not representable exactly as a ratio of
+  /// two signed integers of type `T`, a runtime error may occur.
+  ///
+  /// - Parameters:
+  ///   - source: A binary floating-point value to convert to a rational value.
+  @_transparent // @_inlineable
+  public init<Source : BinaryFloatingPoint>(_ source: Source) {
+    if source.isNaN { self = .nan; return }
+    if source == .infinity { self = .infinity; return }
+    if source == -.infinity { self = -.infinity; return }
+    if source.isZero { self = 0; return }
+
+    let exponent = source.exponent
+    let significandWidth = source.significandWidth
+    let shift = Source.Exponent(significandWidth) - exponent
+
+    if shift <= 0 {
+      self.numerator = T(source)
+      self.denominator = 1
+      return
+    }
+
+    let numerator = T(
+      Source(
+        sign: source.sign,
+        exponent: exponent + shift,
+        significand: source.significand
+      )
+    )
+    let denominator = T(Source(sign: .plus, exponent: shift, significand: 1))
+    // Ensure that `numerator.magnitude` and `denominator.magnitude` are each
+    // representable as a `T`.
+    _ = T(numerator.magnitude)
+    _ = T(denominator.magnitude)
+    self.numerator = numerator
+    self.denominator = denominator
+  }
+}
+
+extension Rational where T : FixedWidthInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Initializers (Constrained)
+  // ---------------------------------------------------------------------------
+
+  /// Creates a new rational value from the given binary floating-point value,
+  /// if it can be represented exactly.
+  ///
+  /// If `source` or its magnitude is not representable exactly as a ratio of
+  /// two signed integers of type `T`, the result is `nil`.
+  ///
+  /// - Note: This initializer creates only instances of
+  ///   `Rational<T> where T : FixedWidthInteger`.
+  ///
+  /// - Parameters:
+  ///   - source: A floating-point value to convert to a rational value.
+  @_transparent // @_inlineable
+  public init?<Source : BinaryFloatingPoint>(exactly source: Source) {
+    // TODO: Document this initializer.
+    if source.isNaN { self = .nan; return }
+    if source == .infinity { self = .infinity; return }
+    if source == -.infinity { self = -.infinity; return }
+    if source.isZero { self = 0; return } // Consider -0.0 to be exactly 0.
+
+    let exponent = source.exponent
+    let significandWidth = source.significandWidth
+    let shift = Source.Exponent(significandWidth) - exponent
+    let bitWidth = T.bitWidth
+
+    if shift <= 0 {
+      guard exponent + 1 < bitWidth else { return nil }
+      self.numerator = T(source)
+      self.denominator = 1
+      return
+    }
+
+    guard significandWidth + 1 < bitWidth && shift + 1 < bitWidth else {
+      return nil
+    }
+    self.numerator = T(
+      Source(
+        sign: source.sign,
+        exponent: exponent + shift,
+        significand: source.significand
+      )
+    )
+    self.denominator = T(Source(sign: .plus, exponent: shift, significand: 1))
+  }
+}
+
+extension Rational {
+  // ---------------------------------------------------------------------------
+  // MARK: Static Properties
+  // ---------------------------------------------------------------------------
+
   /// Positive infinity.
   ///
   /// Infinity compares greater than all finite numbers and equal to other
@@ -116,52 +234,16 @@ where T : _ExpressibleByBuiltinIntegerLiteral, T.Magnitude : UnsignedInteger,
   public static var nan: Rational {
     return Rational(numerator: 0, denominator: 0)
   }
-}
 
-extension Rational {
-  // TODO: Document this initializer.
-  // @_transparent // @_inlineable
-  public init<Source : BinaryInteger>(_ source: Source) {
-    let t = T(source)
-    // Ensure that `t.magnitude` is representable as a `T`.
-    _ = T(t.magnitude)
-    self.numerator = t
-    self.denominator = 1
-  }
+  // ---------------------------------------------------------------------------
+  // MARK: Static Methods
+  // ---------------------------------------------------------------------------
 
-  // TODO: Document this initializer.
-  // @_transparent // @_inlineable
-  public init<Source : BinaryFloatingPoint>(_ source: Source) {
-    if source.isNaN { self = .nan; return }
-    if source == .infinity { self = .infinity; return }
-    if source == -.infinity { self = -.infinity; return }
-    if source.isZero { self = 0; return }
-
-    let exponent = source.exponent
-    let significandWidth = source.significandWidth
-    if significandWidth <= exponent {
-      self.numerator = T(source)
-      self.denominator = 1
-      return
-    }
-    let shift = significandWidth - Int(exponent)
-    let numerator = T(source * Source(1 &<< shift))
-    // Ensure that `numerator.magnitude` is representable as a `T`.
-    _ = T(numerator.magnitude)
-    let denominator = T(1 &<< shift)
-    // Ensure that `denominator.magnitude` is representable as a `T`.
-    _ = T(denominator.magnitude)
-    self.numerator = numerator
-    self.denominator = denominator
-  }
-}
-
-extension Rational {
   /// Compares the (finite) magnitude of two finite values, returning -1 if
   /// `lhs.magnitude` is less than `rhs.magnitude`, 0 if `lhs.magnitude` is
   /// equal to `rhs.magnitude`, or 1 if `lhs.magnitude` is greater than
   /// `rhs.magnitude`.
-  // @_versioned
+  @_versioned
   internal static func _compareFiniteMagnitude(
     _ lhs: Rational, _ rhs: Rational
   ) -> Int {
@@ -182,8 +264,34 @@ extension Rational {
     */
   }
 
+  /// Returns the quotient obtained by dividing the first value by the second,
+  /// trapping in case of arithmetic overflow.
+  ///
+  /// - Parameters:
+  ///   - lhs: The value to divide.
+  ///   - rhs: The value by which to divide `lhs`.
+  @_transparent // @_inlineable
+  public static func / (lhs: Rational, rhs: Rational) -> Rational {
+    return lhs * rhs.reciprocal()
+  }
+
+  /// Divides the left-hand side by the right-hand side and stores the quotient
+  /// in the left-hand side, trapping in case of arithmetic overflow.
+  ///
+  /// - Parameters:
+  ///   - lhs: The value to divide.
+  ///   - rhs: The value by which to divide `lhs`.
+  @_transparent // @_inlineable
+  public static func /= (lhs: inout Rational, rhs: Rational) {
+    lhs = lhs * rhs.reciprocal()
+  }
+
+  // ---------------------------------------------------------------------------
+  // MARK: Computed Properties
+  // ---------------------------------------------------------------------------
+
   /// The canonical representation of this value.
-  // @_transparent // @_inlineable
+  @_transparent // @_inlineable
   public var canonical: Rational {
     let nm = numerator.magnitude, dm = denominator.magnitude
     // Note that if `T` is a signed fixed-width integer type, `gcd(nm, dm)`
@@ -236,10 +344,11 @@ extension Rational {
 
   /// A Boolean value indicating whether the instance is a proper fraction.
   ///
-  /// A fraction `p / q` is proper iff `p > 0`, `q > 0`, and `p < q`.
+  /// A fraction is proper if and only if the absolute value of the fraction is
+  /// less than 1.
   @_transparent // @_inlineable
   public var isProper: Bool {
-    return numerator > 0 && denominator > 0 && numerator < denominator
+    return denominator != 0 && numerator / denominator == 0
   }
 
   /// A Boolean value indicating whether the instance is equal to zero.
@@ -276,6 +385,10 @@ extension Rational {
       : .minus
   }
 
+  // ---------------------------------------------------------------------------
+  // MARK: Methods
+  // ---------------------------------------------------------------------------
+
   /// Returns the reciprocal (multiplicative inverse) of this value.
   @_transparent // @_inlineable
   public func reciprocal() -> Rational {
@@ -283,9 +396,99 @@ extension Rational {
       ? Rational(numerator: -denominator, denominator: -numerator)
       : Rational(numerator: denominator, denominator: numerator)
   }
+
+  /// Returns this value rounded to an integral value using the specified
+  /// rounding rule.
+  ///
+  /// ```swift
+  /// let x = 7 / 2 as Rational<Int>
+  /// print(x.rounded()) // Prints "4"
+  /// print(x.rounded(.towardZero)) // Prints "3"
+  /// print(x.rounded(.up)) // Prints "4"
+  /// print(x.rounded(.down)) // Prints "3"
+  /// ```
+  ///
+  /// See the `FloatingPointRoundingRule` enumeration for more information about
+  /// the available rounding rules.
+  ///
+  /// - Parameters:
+  ///   - rule: The rounding rule to use.
+  ///
+  /// - SeeAlso: `round(_:)`, `RoundingRule`
+  @_transparent // @_inlineable
+  public func rounded(
+    _ rule: RoundingRule = .toNearestOrAwayFromZero
+  ) -> Rational {
+    var t = self
+    t.round(rule)
+    return t
+  }
+
+  /// Rounds the value to an integral value using the specified rounding rule.
+  ///
+  /// ```swift
+  /// var x = 7 / 2 as Rational<Int>
+  /// x.round() // x == 4
+  ///
+  /// var x = 7 / 2 as Rational<Int>
+  /// x.round(.towardZero) // x == 3
+  ///
+  /// var x = 7 / 2 as Rational<Int>
+  /// x.round(.up) // x == 4
+  ///
+  /// var x = 7 / 2 as Rational<Int>
+  /// x.round(.down) // x == 3
+  /// ```
+  ///
+  /// See the `FloatingPointRoundingRule` enumeration for more information about
+  /// the available rounding rules.
+  ///
+  /// - Parameters:
+  ///   - rule: The rounding rule to use.
+  ///
+  /// - SeeAlso: `round(_:)`, `RoundingRule`
+  @_transparent // @_inlineable
+  public mutating func round(_ rule: RoundingRule = .toNearestOrAwayFromZero) {
+    if denominator == 0 { return }
+
+    let f: T
+    (numerator, f) = numerator.quotientAndRemainder(dividingBy: denominator)
+    // Rounding rules only come into play if the fractional part is non-zero.
+    if f != 0 {
+      switch rule {
+      case .toNearestOrAwayFromZero:
+        fallthrough
+      case .toNearestOrEven:
+        switch denominator.magnitude.quotientAndRemainder(
+          dividingBy: f.magnitude
+        ) {
+        case (2, 0): // Tie.
+          if rule == .toNearestOrEven && numerator % 2 == 0 { break }
+          fallthrough
+        case (1, _): // Nearest is away from zero.
+          if f > 0 { numerator += 1 } else { numerator -= 1 }
+        default: // Nearest is toward zero.
+          break
+        }
+      case .up:
+        if f > 0 { numerator += 1 }
+      case .down:
+        if f < 0 { numerator -= 1 }
+      case .towardZero:
+        break
+      case .awayFromZero:
+        if f > 0 { numerator += 1 } else { numerator -= 1 }
+      }
+    }
+    denominator = 1
+  }
 }
 
 extension Rational : ExpressibleByIntegerLiteral {
+  // ---------------------------------------------------------------------------
+  // MARK: ExpressibleByIntegerLiteral
+  // ---------------------------------------------------------------------------
+
   @_transparent // @_inlineable
   public init(integerLiteral value: T) {
     self.numerator = value
@@ -294,6 +497,10 @@ extension Rational : ExpressibleByIntegerLiteral {
 }
 
 extension Rational : CustomStringConvertible {
+  // ---------------------------------------------------------------------------
+  // MARK: CustomStringConvertible
+  // ---------------------------------------------------------------------------
+
   @_transparent // @_inlineable
   public var description: String {
     if numerator == 0 { return denominator == 0 ? "nan" : "0" }
@@ -303,7 +510,11 @@ extension Rational : CustomStringConvertible {
 }
 
 extension Rational : Equatable {
-  // @_transparent // @_inlineable
+  // ---------------------------------------------------------------------------
+  // MARK: Equatable
+  // ---------------------------------------------------------------------------
+
+  @_transparent // @_inlineable
   public static func == (lhs: Rational, rhs: Rational) -> Bool {
     if lhs.denominator == 0 {
       if lhs.numerator == 0 { return false }
@@ -317,15 +528,23 @@ extension Rational : Equatable {
 }
 
 extension Rational : Hashable {
+  // ---------------------------------------------------------------------------
+  // MARK: Hashable
+  // ---------------------------------------------------------------------------
+
   // @_transparent // @_inlineable
   public var hashValue: Int {
     let t = canonical
-    return _fnv1a(t.numerator, t.denominator)
+    return _Hash._combine(t.numerator, t.denominator)
   }
 }
 
 extension Rational : Comparable {
-  // @_transparent // @_inlineable
+  // ---------------------------------------------------------------------------
+  // MARK: Comparable
+  // ---------------------------------------------------------------------------
+
+  @_transparent // @_inlineable
   public static func < (lhs: Rational, rhs: Rational) -> Bool {
     if lhs.denominator == 0 {
       if lhs.numerator >= 0 { return false }
@@ -350,7 +569,7 @@ extension Rational : Comparable {
     return rhs < lhs
   }
 
-  // @_transparent // @_inlineable
+  @_transparent // @_inlineable
   public static func <= (lhs: Rational, rhs: Rational) -> Bool {
     if lhs.denominator == 0 {
       if lhs.numerator == 0 { return false }
@@ -377,7 +596,11 @@ extension Rational : Comparable {
   }
 }
 
-extension Rational : Strideable, _Strideable {
+extension Rational : Strideable {
+  // ---------------------------------------------------------------------------
+  // MARK: Strideable
+  // ---------------------------------------------------------------------------
+
   @_transparent // @_inlineable
   public func distance(to other: Rational) -> Rational {
     return other - self
@@ -389,38 +612,252 @@ extension Rational : Strideable, _Strideable {
   }
 }
 
-extension Rational
-where T : FixedWidthInteger, T.Magnitude : FixedWidthInteger {
-  // TODO: Document this initializer.
-  // @_transparent // @_inlineable
-  public init?<Source : BinaryFloatingPoint>(exactly source: Source) {
-    if source.isNaN { self = .nan; return }
-    if source == .infinity { self = .infinity; return }
-    if source == -.infinity { self = -.infinity; return }
-    if source.isZero { self = 0; return } // Consider -0.0 to be exactly 0.
+extension Rational : Numeric {
+  // ---------------------------------------------------------------------------
+  // MARK: Numeric
+  // ---------------------------------------------------------------------------
 
-    let exponent = source.exponent
-    let significandWidth = source.significandWidth
-    let bitWidth = T.bitWidth
-    if significandWidth <= exponent {
-      guard exponent + 1 < bitWidth else { return nil }
-      self.numerator = T(source)
-      self.denominator = 1
-      return
+  @_transparent // @_inlineable
+  public init?<U>(exactly source: U) where U : BinaryInteger {
+    guard let t = T(exactly: source) else { return nil }
+    // Ensure that `t.magnitude` is representable as a `T`.
+    guard let _ = T(exactly: t.magnitude) else { return nil }
+    self.numerator = t
+    self.denominator = 1
+  }
+
+  @_transparent // @_inlineable
+  public static func + (lhs: Rational, rhs: Rational) -> Rational {
+    if lhs.denominator == 0 {
+      if rhs.denominator != 0 || lhs.numerator == 0 { return lhs }
+      if lhs.numerator > 0 { return rhs.numerator < 0 ? .nan : rhs }
+      return rhs.numerator > 0 ? .nan : rhs
     }
-    let shift = significandWidth - Int(exponent)
-    guard significandWidth + 1 < bitWidth && shift < bitWidth else {
-      return nil
+    if rhs.denominator == 0 { return rhs }
+
+    let ldm = lhs.denominator.magnitude
+    let rdm = rhs.denominator.magnitude
+    let gcd = T.Magnitude.gcd(ldm, rdm)
+    let a = T(rdm / gcd * lhs.numerator.magnitude)
+    let b = T(ldm / gcd * rhs.numerator.magnitude)
+    let n = lhs.sign == .plus
+      ? (rhs.sign == .plus ? a + b : a - b)
+      : (rhs.sign == .plus ? b - a : -a - b)
+    let d = T(ldm / gcd * rdm)
+    return Rational(numerator: n, denominator: d).canonical
+  }
+
+  @_transparent // @_inlineable
+  public static func += (lhs: inout Rational, rhs: Rational) {
+    lhs = lhs + rhs
+  }
+
+  @_transparent // @_inlineable
+  public static func - (lhs: Rational, rhs: Rational) -> Rational {
+    return lhs + (-rhs)
+  }
+
+  @_transparent // @_inlineable
+  public static func -= (lhs: inout Rational, rhs: Rational) {
+    lhs = lhs + (-rhs)
+  }
+  
+  @_transparent // @_inlineable
+  public static func * (lhs: Rational, rhs: Rational) -> Rational {
+    if lhs.denominator == 0 {
+      if rhs.numerator == 0 { return .nan }
+      return rhs.sign == .plus ? lhs : -lhs
     }
-    let numerator = T(source * Source(1 &<< shift))
-    // Ensure that `numerator.magnitude` is representable as a `T`.
-    guard let _ = T(exactly: numerator.magnitude) else { return nil }
-    let denominator = T(1 &<< shift)
-    // Ensure that `denominator.magnitude` is representable as a `T`.
-    guard let _ = T(exactly: denominator.magnitude) else { return nil }
-    self.numerator = numerator
-    self.denominator = denominator
+    if rhs.denominator == 0 {
+      if lhs.numerator == 0 { return .nan }
+      return lhs.sign == .plus ? rhs : -rhs
+    }
+    
+    let lnm = lhs.numerator.magnitude, ldm = lhs.denominator.magnitude
+    let rnm = rhs.numerator.magnitude, rdm = rhs.denominator.magnitude
+    // Note that if `T` is a signed fixed-width integer type, `gcd(lnm, rdm)` or
+    // `gcd(rnm, ldm)` could be equal to `-T.min`, which is not representable as
+    // a `T`. This is why the following arithmetic is performed with values of
+    // type `T.Magnitude`.
+    let a = T.Magnitude.gcd(lnm, rdm)
+    let b = T.Magnitude.gcd(rnm, ldm)
+    let n = lhs.sign == rhs.sign
+      ? T(lnm / a * (rnm / b))
+      : -T(lnm / a * (rnm / b))
+    let d = T(ldm / b * (rdm / a))
+    return Rational(numerator: n, denominator: d)
+  }
+
+  @_transparent // @_inlineable
+  public static func *= (lhs: inout Rational, rhs: Rational) {
+    lhs = lhs * rhs
   }
 }
 
+extension Rational : SignedNumeric {
+  // ---------------------------------------------------------------------------
+  // MARK: SignedNumeric
+  // ---------------------------------------------------------------------------
+
+  @_transparent // @_inlineable
+  public static prefix func - (operand: Rational) -> Rational {
+    return Rational(
+      numerator: -operand.numerator, denominator: operand.denominator
+    )
+  }
+
+  @_transparent // @_inlineable
+  public mutating func negate() {
+    numerator.negate()
+  }
+}
+
+/// Returns the absolute value (magnitude) of `x`.
+@_transparent
+public func abs<T>(_ x: Rational<T>) -> Rational<T> {
+  return x.magnitude
+}
+
+/// Returns the closest integral value greater than or equal to `x`.
+@_transparent
+public func ceil<T>(_ x: Rational<T>) -> Rational<T> {
+  return x.rounded(.up)
+}
+
+/// Returns the closest integral value less than or equal to `x`.
+@_transparent
+public func floor<T>(_ x: Rational<T>) -> Rational<T> {
+  return x.rounded(.down)
+}
+
+/// Returns the closest integral value; if two values are equally close, returns
+/// the one with greater magnitude.
+@_transparent
+public func round<T>(_ x: Rational<T>) -> Rational<T> {
+  return x.rounded()
+}
+
+/// Returns the closest integral value with magnitude less than or equal to that
+/// of `x`.
+@_transparent
+public func trunc<T>(_ x: Rational<T>) -> Rational<T> {
+  return x.rounded(.towardZero)
+}
+
 public typealias Ratio = Rational<Int>
+
+// MARK: -
+
+extension BinaryInteger {
+  // ---------------------------------------------------------------------------
+  // MARK: Initializers
+  // ---------------------------------------------------------------------------
+
+  /// Creates a new binary integer from the given rational value, if it can be
+  /// represented exactly.
+  ///
+  /// If `source` is not representable exactly, the result is `nil`.
+  ///
+  /// - Parameters:
+  ///   - source: A rational value to convert to a binary integer.
+  @_transparent // @_inlineable
+  public init?<U>(exactly source: Rational<U>) {
+    let (whole, fraction) = source.mixed
+    guard fraction.isZero, let exact = Self(exactly: whole) else { return nil }
+    self = exact
+  }
+
+  /// Creates a new binary integer from the given rational value, rounding
+  /// toward zero.
+  ///
+  /// If `source` is outside the bounds of this type after rounding toward zero,
+  /// a runtime error may occur.
+  ///
+  /// - Parameters:
+  ///   - source: A rational value to convert to a binary integer.
+  @_transparent // @_inlineable
+  public init<U>(_ source: Rational<U>) {
+    self = Self(source.mixed.whole)
+  }
+}
+
+extension FloatingPoint {
+  // ---------------------------------------------------------------------------
+  // MARK: Initializers
+  // ---------------------------------------------------------------------------
+
+  /// Creates a new floating-point value from the given rational value, after
+  /// rounding the whole part, the numerator of the fractional part, and the
+  /// denominator of the fractional part each to the closest possible
+  /// representation.
+  ///
+  /// If two representable values are equally close, the result of rounding is
+  /// the value with more trailing zeros in its significand bit pattern.
+  ///
+  /// - Parameters:
+  ///   - source: The rational value to convert to a floating-point value.
+  public init(_ source: Rational<Int>) {
+    let (whole, fraction) = source.mixed
+    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
+  }
+
+  /// Creates a new floating-point value from the given rational value, after
+  /// rounding the whole part, the numerator of the fractional part, and the
+  /// denominator of the fractional part each to the closest possible
+  /// representation.
+  ///
+  /// If two representable values are equally close, the result of rounding is
+  /// the value with more trailing zeros in its significand bit pattern.
+  ///
+  /// - Parameters:
+  ///   - source: The rational value to convert to a floating-point value.
+  public init(_ source: Rational<Int8>) {
+    let (whole, fraction) = source.mixed
+    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
+  }
+
+  /// Creates a new floating-point value from the given rational value, after
+  /// rounding the whole part, the numerator of the fractional part, and the
+  /// denominator of the fractional part each to the closest possible
+  /// representation.
+  ///
+  /// If two representable values are equally close, the result of rounding is
+  /// the value with more trailing zeros in its significand bit pattern.
+  ///
+  /// - Parameters:
+  ///   - source: The rational value to convert to a floating-point value.
+  public init(_ source: Rational<Int16>) {
+    let (whole, fraction) = source.mixed
+    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
+  }
+
+  /// Creates a new floating-point value from the given rational value, after
+  /// rounding the whole part, the numerator of the fractional part, and the
+  /// denominator of the fractional part each to the closest possible
+  /// representation.
+  ///
+  /// If two representable values are equally close, the result of rounding is
+  /// the value with more trailing zeros in its significand bit pattern.
+  ///
+  /// - Parameters:
+  ///   - source: The rational value to convert to a floating-point value.
+  public init(_ source: Rational<Int32>) {
+    let (whole, fraction) = source.mixed
+    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
+  }
+
+  /// Creates a new floating-point value from the given rational value, after
+  /// rounding the whole part, the numerator of the fractional part, and the
+  /// denominator of the fractional part each to the closest possible
+  /// representation.
+  ///
+  /// If two representable values are equally close, the result of rounding is
+  /// the value with more trailing zeros in its significand bit pattern.
+  ///
+  /// - Parameters:
+  ///   - source: The rational value to convert to a floating-point value.
+  public init(_ source: Rational<Int64>) {
+    let (whole, fraction) = source.mixed
+    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
+  }
+}

Sources/Real.swift

@@ -1,5 +1,5 @@
 //
-//  FloatingPointMath.swift
+//  Real.swift
 //  NumericAnnex
 //
 //  Created by Xiaodi Wu on 4/1/17.
@@ -11,71 +11,12 @@ import Glibc
 import Darwin.C
 #endif
 
-/// A floating-point type that provides a selection of special functions.
+/// A floating-point type that supports elementary functions and a selection of
+/// special functions.
 ///
-/// The `FloatingPointMath` protocol provides a suitable basis for writing
-/// functions that work on any floating-point type that provides the required
-/// functions.
-public protocol FloatingPointMath : Math, FloatingPoint /*, Hashable */ {
-  /// Creates a new value from the given rational value, after rounding the
-  /// whole part and the numerator and denominator of the fractional part each
-  /// to the closest possible representation.
-  ///
-  /// If two representable values are equally close, the result is the value
-  /// with more trailing zeros in its significand bit pattern.
-  ///
-  /// - Parameters:
-  ///   - value: The rational value to convert to a floating-point value.
-  init(_ value: Rational<Int>)
-
-  /// Creates a new value from the given rational value, after rounding the
-  /// whole part and the numerator and denominator of the fractional part each
-  /// to the closest possible representation.
-  ///
-  /// If two representable values are equally close, the result is the value
-  /// with more trailing zeros in its significand bit pattern.
-  ///
-  /// - Parameters:
-  ///   - value: The rational value to convert to a floating-point value.
-  init(_ value: Rational<Int8>)
-
-  /// Creates a new value from the given rational value, after rounding the
-  /// whole part and the numerator and denominator of the fractional part each
-  /// to the closest possible representation.
-  ///
-  /// If two representable values are equally close, the result is the value
-  /// with more trailing zeros in its significand bit pattern.
-  ///
-  /// - Parameters:
-  ///   - value: The rational value to convert to a floating-point value.
-  init(_ value: Rational<Int16>)
-
-  /// Creates a new value from the given rational value, after rounding the
-  /// whole part and the numerator and denominator of the fractional part each
-  /// to the closest possible representation.
-  ///
-  /// If two representable values are equally close, the result is the value
-  /// with more trailing zeros in its significand bit pattern.
-  ///
-  /// - Parameters:
-  ///   - value: The rational value to convert to a floating-point value.
-  init(_ value: Rational<Int32>)
-
-  /// Creates a new value from the given rational value, after rounding the
-  /// whole part and the numerator and denominator of the fractional part each
-  /// to the closest possible representation.
-  ///
-  /// If two representable values are equally close, the result is the value
-  /// with more trailing zeros in its significand bit pattern.
-  ///
-  /// - Parameters:
-  ///   - value: The rational value to convert to a floating-point value.
-  init(_ value: Rational<Int64>)
-
-  // FIXME: If corresponding requirements are added to FloatingPoint
-  // add init<U : SignedInteger>(_: Rational<U>) as well as
-  // init?<U : SignedInteger>(exactly: Rational<U>).
-
+/// The `Real` protocol provides a suitable basis for writing functions that
+/// work on any floating-point type that supports the required functions.
+public protocol Real : Math, FloatingPoint {
   /// Returns the hypotenuse of a right-angle triangle with legs (catheti) of
   /// length `x` and `y`, preventing avoidable arithmetic overflow and
   /// underflow. The return value is the square root of the sum of squares of
@@ -124,43 +65,36 @@ public protocol FloatingPointMath : Math, FloatingPoint /*, Hashable */ {
   ///
   /// - SeeAlso: `lgamma(_:)`
   func logarithmicGamma() -> Self
+
+  // FIXME: Design and implement `logarithmicGammaAndSignOfGamma()`.
+
+  // FIXME: If corresponding requirements are added to FloatingPoint
+  // add init<U : SignedInteger>(_: Rational<U>) as well as
+  // init?<U : SignedInteger>(exactly: Rational<U>).
 }
 
-extension FloatingPointMath {
-  public init(_ value: Rational<Int>) {
-    let (whole, fraction) = value.mixed
-    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
-  }
-
-  public init(_ value: Rational<Int8>) {
-    let (whole, fraction) = value.mixed
-    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
-  }
-
-  public init(_ value: Rational<Int16>) {
-    let (whole, fraction) = value.mixed
-    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
-  }
-
-  public init(_ value: Rational<Int32>) {
-    let (whole, fraction) = value.mixed
-    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
-  }
-
-  public init(_ value: Rational<Int64>) {
-    let (whole, fraction) = value.mixed
-    self = Self(whole) + Self(fraction.numerator) / Self(fraction.denominator)
-  }
-
+extension Real {
   public static func hypot(_ x: Self, _ y: Self) -> Self {
     var x = abs(x), y = abs(y)
-    if x < y {
-      swap(&x, &y)
-    }
+    if x.isInfinite { return x }
+    if y.isInfinite { return y }
+    if x == 0 { return y }
+    if x < y { swap(&x, &y) }
     let ratio = y / x
     return x * .sqrt(1 + ratio * ratio)
   }
 
+  /// Returns the cube root of the value, rounded to a representable value
+  /// (default implementation).
+  ///
+  /// - Note: This method is required by the protocol `Math`; it is not an
+  ///   extension method.
+  ///
+  /// - SeeAlso: `cbrt(_:)`
+  public func cubeRoot() -> Self {
+    return sign == .minus ? -.exp(.log(-self) / 3) : .exp(.log(self) / 3)
+  }
+
   public func inverseTangent(dividingBy other: Self) -> Self {
     let y = self, x = other
     if y.isInfinite {
@@ -188,7 +122,7 @@ extension FloatingPointMath {
   }
 }
 
-extension FloatingPointMath {
+extension Real {
   /// Returns the inverse tangent of `y / x`, using the signs of `y` and `x` to
   /// determine the quadrant of the computed angle.
   ///
@@ -249,7 +183,11 @@ extension FloatingPointMath {
   }
 }
 
-extension Float : FloatingPointMath {
+extension Float : Real {
+  // ---------------------------------------------------------------------------
+  // MARK: Real
+  // ---------------------------------------------------------------------------
+
   @_transparent
   public static var e: Float {
     return Float(0x1.5bf0a8p1)
@@ -260,14 +198,24 @@ extension Float : FloatingPointMath {
     return Float(0x1.9e377ap0)
   }
 
+  @_transparent
+  public static func ** (lhs: Float, rhs: Float) -> Float {
+    return powf(lhs, rhs)
+  }
+
+  @_transparent
+  public static func **= (lhs: inout Float, rhs: Float) {
+    lhs = lhs ** rhs
+  }
+
   @_transparent
   public func naturalExponential() -> Float {
-    return expf(self)
+    return _exp(self)
   }
 
   @_transparent
   public func binaryExponential() -> Float {
-    return exp2f(self)
+    return _exp2(self)
   }
 
   @_transparent
@@ -286,17 +234,17 @@ extension Float : FloatingPointMath {
 
   @_transparent
   public func naturalLogarithm() -> Float {
-    return logf(self)
+    return _log(self)
   }
 
   @_transparent
   public func binaryLogarithm() -> Float {
-    return log2f(self)
+    return _log2(self)
   }
 
   @_transparent
   public func commonLogarithm() -> Float {
-    return log10f(self)
+    return _log10(self)
   }
 
   @_transparent
@@ -309,19 +257,14 @@ extension Float : FloatingPointMath {
     return cbrtf(self)
   }
 
-  @_transparent
-  public func power(of base: Float) -> Float {
-    return powf(base, self)
-  }
-
   @_transparent
   public func sine() -> Float {
-    return sinf(self)
+    return _sin(self)
   }
 
   @_transparent
   public func cosine() -> Float {
-    return cosf(self)
+    return _cos(self)
   }
 
   @_transparent
@@ -405,7 +348,11 @@ extension Float : FloatingPointMath {
   }
 }
 
-extension Double : FloatingPointMath {
+extension Double : Real {
+  // ---------------------------------------------------------------------------
+  // MARK: Real
+  // ---------------------------------------------------------------------------
+
   @_transparent
   public static var e: Double {
     return Double(0x1.5bf0a8b145769p1)
@@ -417,21 +364,27 @@ extension Double : FloatingPointMath {
   }
 
   @_transparent
-  public func naturalExponential() -> Double {
+  public static func ** (lhs: Double, rhs: Double) -> Double {
 #if os(Linux)
-    return Glibc.exp(self)
+    return Glibc.pow(lhs, rhs)
 #else
-    return Darwin.exp(self)
+    return Darwin.pow(lhs, rhs)
 #endif
   }
 
+  @_transparent
+  public static func **= (lhs: inout Double, rhs: Double) {
+    lhs = lhs ** rhs
+  }
+
+  @_transparent
+  public func naturalExponential() -> Double {
+    return _exp(self)
+  }
+
   @_transparent
   public func binaryExponential() -> Double {
-#if os(Linux)
-    return Glibc.exp2(self)
-#else
-    return Darwin.exp2(self)
-#endif
+    return _exp2(self)
   }
 
   @_transparent
@@ -454,29 +407,17 @@ extension Double : FloatingPointMath {
 
   @_transparent
   public func naturalLogarithm() -> Double {
-#if os(Linux)
-    return Glibc.log(self)
-#else
-    return Darwin.log(self)
-#endif
+    return _log(self)
   }
 
   @_transparent
   public func binaryLogarithm() -> Double {
-#if os(Linux)
-    return Glibc.log2(self)
-#else
-    return Darwin.log2(self)
-#endif
+    return _log2(self)
   }
 
   @_transparent
   public func commonLogarithm() -> Double {
-#if os(Linux)
-    return Glibc.log10(self)
-#else
-    return Darwin.log10(self)
-#endif
+    return _log10(self)
   }
 
   @_transparent
@@ -497,31 +438,14 @@ extension Double : FloatingPointMath {
 #endif
   }
 
-  @_transparent
-  public func power(of base: Double) -> Double {
-#if os(Linux)
-    return Glibc.pow(base, self)
-#else
-    return Darwin.pow(base, self)
-#endif
-  }
-
   @_transparent
   public func sine() -> Double {
-#if os(Linux)
-    return Glibc.sin(self)
-#else
-    return Darwin.sin(self)
-#endif
+    return _sin(self)
   }
 
   @_transparent
   public func cosine() -> Double {
-#if os(Linux)
-    return Glibc.cos(self)
-#else
-    return Darwin.cos(self)
-#endif
+    return _cos(self)
   }
 
   @_transparent

Tests/LinuxMain.swift

@@ -2,7 +2,11 @@ import XCTest
 @testable import NumericAnnexTests
 
 XCTMain([
+  testCase(DocumentationExampleTests.allTests),
+  testCase(ExponentiationTests.allTests),
+  testCase(RootExtractionTests.allTests),
   testCase(FactoringTests.allTests),
+  testCase(RealTests.allTests),
   testCase(RationalTests.allTests),
   testCase(ComplexTests.allTests),
   testCase(RandomXorshiftTests.allTests),

Tests/NumericAnnexTests/DocumentationExampleTests.swift

@@ -0,0 +1,58 @@
+import XCTest
+/* @testable */ import NumericAnnex
+
+class DocumentationExampleTests : XCTestCase {
+  // It reflects poorly on the project to provide examples that don't compile or
+  // give the wrong result, so test them here:
+  
+  func testReadmeExample() {
+    print(2 ** 3)
+    // Prints "8".
+    XCTAssertEqual("\(2 ** 3)", "8")
+
+    print(4.0 ** 5.0)
+    // Prints "1024.0".
+    XCTAssertEqual("\(4.0 ** 5.0)", "1024.0")
+
+    print(Int.cbrt(8))
+    // Prints "2".
+    XCTAssertEqual("\(Int.cbrt(8))", "2")
+
+    print(Double.cbrt(27.0))
+    // Prints "3.0".
+    XCTAssertEqual("\(Double.cbrt(27.0))", "3.0")
+
+    var x: Ratio = 1 / 4
+    // Ratio is a type alias for Rational<Int>.
+
+    print(x.reciprocal())
+    // Prints "4".
+    XCTAssertEqual("\(x.reciprocal())", "4")
+
+    x *= 8
+    print(x + x)
+    // Prints "4".
+    XCTAssertEqual("\(x + x)", "4")
+
+    x = Ratio(Float.phi) // Golden ratio.
+    print(x)
+    // Prints "13573053/8388608".
+    XCTAssertEqual("\(x)", "13573053/8388608")
+
+    var z: Complex64 = 42 * .i
+    // Complex64 is a type alias for Complex<Float>.
+
+    print(Complex.sqrt(z))
+    // Prints "4.58258 + 4.58258i".
+    XCTAssertEqual("\(Complex.sqrt(z))", "4.58258 + 4.58258i")
+
+    z = .pi + .i * .log(2 - .sqrt(3))
+    print(Complex.cos(z).real)
+    // Prints "-2.0".
+    XCTAssertEqual("\(Complex.cos(z).real)", "-2.0")
+  }
+
+  static var allTests = [
+    ("testReadmeExample", testReadmeExample),
+  ]
+}

Tests/NumericAnnexTests/ExponentiationTests.swift

@@ -0,0 +1,73 @@
+import XCTest
+@testable import NumericAnnex
+
+class ExponentiationTests : XCTestCase {
+  func testIntPow() {
+    for i in -21...21 {
+      for j in -7...7 {
+        if i == 0 && j < 0 {
+          // In this case, i ** j causes a division-by-zero error.
+          // Likewise, Double(i) ** Double(j) also causes the same error.
+          continue
+        }
+        let actual = i ** j
+        let expected = Int(Double(i) ** Double(j))
+        XCTAssertEqual(actual, expected)
+      }
+    }
+  }
+
+  func testUIntPow() {
+    for i in (0 as UInt)...21 {
+      for j in (0 as UInt)...7 {
+        let actual = i ** j
+        let expected = UInt(Double(i) ** Double(j))
+        XCTAssertEqual(actual, expected)
+      }
+    }
+  }
+
+  func testConcreteIntegerTypePow() {
+    XCTAssertEqual((2 as Int) ** 3, 8 as Int)
+    XCTAssertEqual((2 as Int8) ** 3, 8 as Int8)
+    XCTAssertEqual((2 as Int16) ** 3, 8 as Int16)
+    XCTAssertEqual((2 as Int32) ** 3, 8 as Int32)
+    XCTAssertEqual((2 as Int64) ** 3, 8 as Int64)
+
+    do { var x = 2 as Int; x **= 3; XCTAssertEqual(x, 8 as Int) }
+    do { var x = 2 as Int8; x **= 3; XCTAssertEqual(x, 8 as Int8) }
+    do { var x = 2 as Int16; x **= 3; XCTAssertEqual(x, 8 as Int16) }
+    do { var x = 2 as Int32; x **= 3; XCTAssertEqual(x, 8 as Int32) }
+    do { var x = 2 as Int64; x **= 3; XCTAssertEqual(x, 8 as Int64) }
+
+    XCTAssertEqual((2 as UInt) ** 3, 8 as UInt)
+    XCTAssertEqual((2 as UInt8) ** 3, 8 as UInt8)
+    XCTAssertEqual((2 as UInt16) ** 3, 8 as UInt16)
+    XCTAssertEqual((2 as UInt32) ** 3, 8 as UInt32)
+    XCTAssertEqual((2 as UInt64) ** 3, 8 as UInt64)
+
+    do { var x = 2 as UInt; x **= 3; XCTAssertEqual(x, 8 as UInt) }
+    do { var x = 2 as UInt8; x **= 3; XCTAssertEqual(x, 8 as UInt8) }
+    do { var x = 2 as UInt16; x **= 3; XCTAssertEqual(x, 8 as UInt16) }
+    do { var x = 2 as UInt32; x **= 3; XCTAssertEqual(x, 8 as UInt32) }
+    do { var x = 2 as UInt64; x **= 3; XCTAssertEqual(x, 8 as UInt64) }
+  }
+
+  func testGenericIntegerTypePow() {
+    func _testPow<
+      T : BinaryInteger
+    >(_: T.Type = T.self, lhs: T, rhs: T, expected: T) {
+      XCTAssertEqual(lhs ** rhs, expected)
+      do { var x = lhs; x **= rhs; XCTAssertEqual(x, expected) }
+    }
+    _testPow(lhs: 2 as Int, rhs: 3, expected: 8)
+    _testPow(lhs: 2 as UInt, rhs: 3, expected: 8)
+  }
+
+  static var allTests = [
+    ("testIntPow", testIntPow),
+    ("testUIntPow", testUIntPow),
+    ("testConcreteIntegerTypePow", testConcreteIntegerTypePow),
+    ("testGenericIntegerTypePow", testGenericIntegerTypePow),
+  ]
+}

Tests/NumericAnnexTests/FactoringTests.swift

@@ -1,28 +1,77 @@
 import XCTest
 @testable import NumericAnnex
 
-class FactoringTests: XCTestCase {
+class FactoringTests : XCTestCase {
   func testGCD() {
     XCTAssertEqual(UInt.gcd(18, 84), 6)
     XCTAssertEqual(UInt.gcd(24, 60), 12)
     XCTAssertEqual(UInt.gcd(42, 56), 14)
 
+    XCTAssertEqual(Int.gcd(18, 84), 6)
+    XCTAssertEqual(Int.gcd(24, 60), 12)
+    XCTAssertEqual(Int.gcd(42, 56), 14)
+
+    XCTAssertTrue(Int8.gcdReportingOverflow(-128, -128).overflow)
+
     // Test special values.
     XCTAssertEqual(UInt.gcd(0, 42), 42)
     XCTAssertEqual(UInt.gcd(42, 0), 42)
     XCTAssertEqual(UInt.gcd(0, 0), 0)
-
-    XCTAssertTrue(Int8.gcdReportingOverflow(-128, -128).overflow == .overflow)
   }
 
   func testLCM() {
     XCTAssertEqual(UInt.lcm(4, 6), 12)
     XCTAssertEqual(UInt.lcm(6, 21), 42)
 
+    XCTAssertEqual(Int.lcm(4, 6), 12)
+    XCTAssertEqual(Int.lcm(6, 21), 42)
+
+    XCTAssertEqual(Int8.lcmFullWidth(33, 48).high, 2)
+    XCTAssertEqual(Int8.lcmFullWidth(33, 48).low, 16)
+    XCTAssertEqual(Int8.lcmReportingOverflow(33, 48).partialValue, 16)
+    XCTAssertTrue(Int8.lcmReportingOverflow(33, 48).overflow)
+
+    XCTAssertEqual(UInt8.lcmFullWidth(33, 48).high, 2)
+    XCTAssertEqual(UInt8.lcmFullWidth(33, 48).low, 16)
+    XCTAssertEqual(UInt8.lcmReportingOverflow(33, 48).partialValue, 16)
+    XCTAssertTrue(UInt8.lcmReportingOverflow(33, 48).overflow)
+
     // Test special values.
     XCTAssertEqual(UInt.lcm(0, 42), 0)
     XCTAssertEqual(UInt.lcm(42, 0), 0)
     XCTAssertEqual(UInt.lcm(0, 0), 0)
+
+#if !swift(>=4.1)
+    XCTAssertEqual(Int8.lcmFullWidth(0, 42).high, 0)
+    XCTAssertEqual(Int8.lcmFullWidth(0, 42).low, 0)
+    XCTAssertEqual(Int8.lcmReportingOverflow(0, 42).partialValue, 0)
+    XCTAssertFalse(Int8.lcmReportingOverflow(0, 42).overflow)
+
+    XCTAssertEqual(Int8.lcmFullWidth(42, 0).high, 0)
+    XCTAssertEqual(Int8.lcmFullWidth(42, 0).low, 0)
+    XCTAssertEqual(Int8.lcmReportingOverflow(42, 0).partialValue, 0)
+    XCTAssertFalse(Int8.lcmReportingOverflow(42, 0).overflow)
+
+    XCTAssertEqual(Int8.lcmFullWidth(0, 0).high, 0)
+    XCTAssertEqual(Int8.lcmFullWidth(0, 0).low, 0)
+    XCTAssertEqual(Int8.lcmReportingOverflow(0, 0).partialValue, 0)
+    XCTAssertFalse(Int8.lcmReportingOverflow(0, 0).overflow)
+#endif
+
+    XCTAssertEqual(UInt8.lcmFullWidth(0, 42).high, 0)
+    XCTAssertEqual(UInt8.lcmFullWidth(0, 42).low, 0)
+    XCTAssertEqual(UInt8.lcmReportingOverflow(0, 42).partialValue, 0)
+    XCTAssertFalse(UInt8.lcmReportingOverflow(0, 42).overflow)
+
+    XCTAssertEqual(UInt8.lcmFullWidth(42, 0).high, 0)
+    XCTAssertEqual(UInt8.lcmFullWidth(42, 0).low, 0)
+    XCTAssertEqual(UInt8.lcmReportingOverflow(42, 0).partialValue, 0)
+    XCTAssertFalse(UInt8.lcmReportingOverflow(42, 0).overflow)
+
+    XCTAssertEqual(UInt8.lcmFullWidth(0, 0).high, 0)
+    XCTAssertEqual(UInt8.lcmFullWidth(0, 0).low, 0)
+    XCTAssertEqual(UInt8.lcmReportingOverflow(0, 0).partialValue, 0)
+    XCTAssertFalse(UInt8.lcmReportingOverflow(0, 0).overflow)
   }
 
   static var allTests = [

Tests/NumericAnnexTests/RandomXoroshiroTests.swift

@@ -1,7 +1,20 @@
 import XCTest
 @testable import NumericAnnex
 
-class RandomXoroshiroTests: XCTestCase {
+class RandomXoroshiroTests : XCTestCase {
+  func testEntropy() {
+    guard let rng = Random.Xoroshiro() else {
+      XCTFail()
+      return
+    }
+    XCTAssertTrue(rng.state != (0, 0))
+
+    let mock = {
+      return nil as [UInt64]?
+    }
+    XCTAssertNil(Random.Xoroshiro(_entropy: mock()))
+  }
+
   func test_1_2() {
     let rng = Random.Xoroshiro(state: (1, 2))
     // Values were generated using the reference implementation found at:
@@ -427,6 +440,7 @@ class RandomXoroshiroTests: XCTestCase {
   }
 
   static var allTests = [
+    ("testEntropy", testEntropy),
     ("test_1_2", test_1_2),
     ("test_42_42", test_42_42),
   ]

Tests/NumericAnnexTests/RandomXorshiftTests.swift

@@ -1,7 +1,20 @@
 import XCTest
 @testable import NumericAnnex
 
-class RandomXorshiftTests: XCTestCase {
+class RandomXorshiftTests : XCTestCase {
+  func testEntropy() {
+    guard let rng = Random() else {
+      XCTFail()
+      return
+    }
+    XCTAssertTrue(rng.state != (0, 0))
+
+    let mock = {
+      return nil as [UInt64]?
+    }
+    XCTAssertNil(Random(_entropy: mock()))
+  }
+
   func test_1_2() {
     let rng = Random(state: (1, 2))
     // Values were generated using the reference implementation found at:
@@ -427,6 +440,7 @@ class RandomXorshiftTests: XCTestCase {
   }
 
   static var allTests = [
+    ("testEntropy", testEntropy),
     ("test_1_2", test_1_2),
     ("test_42_42", test_42_42),
   ]

Tests/NumericAnnexTests/RationalTests.swift

@@ -1,16 +1,26 @@
 import XCTest
 @testable import NumericAnnex
 
-class RationalTests: XCTestCase {
+class RationalTests : XCTestCase {
   func testRational() {
     let a = 6 / 4 as Rational<Int>
+    XCTAssertEqual(a.magnitude, a)
+    XCTAssertEqual(abs(a), a)
     XCTAssertEqual(a.description, "3/2")
     XCTAssertEqual(a, 3 / 2 as Rational<Int>)
+    XCTAssertTrue(a.isFinite)
+    XCTAssertFalse(a.isProper)
+    XCTAssertTrue(a.reciprocal().isProper)
 
     let b = 5 / 3 as Rational<Int>
+    XCTAssertEqual(b.magnitude, b)
+    XCTAssertEqual(abs(b), b)
     XCTAssertEqual(b.description, "5/3")
-    XCTAssertLessThan(a, b)
+    XCTAssertTrue(b.isFinite)
+    XCTAssertFalse(b.isProper)
+    XCTAssertTrue(b.reciprocal().isProper)
 
+    XCTAssertLessThan(a, b)
     XCTAssertEqual(a + b, 19 / 6 as Ratio)
     XCTAssertEqual(b + a, 19 / 6 as Ratio)
     XCTAssertEqual(a - b, -1 / 6 as Ratio)
@@ -20,13 +30,91 @@ class RationalTests: XCTestCase {
     XCTAssertEqual(a / b, 9 / 10 as Ratio)
     XCTAssertEqual(b / a, 10 / 9 as Ratio)
 
+    var r = a
+    r += b
+    XCTAssertEqual(r, a + b)
+    r = a
+    r -= b
+    XCTAssertEqual(r, a - b)
+    r = b
+    r -= a
+    XCTAssertEqual(r, b - a)
+    r = a
+    r *= b
+    XCTAssertEqual(r, a * b)
+    r = a
+    r /= b
+    XCTAssertEqual(r, a / b)
+    r = b
+    r /= a
+    XCTAssertEqual(r, b / a)
+
+    let c = -6 / 4 as Rational<Int>
+    XCTAssertEqual(c.magnitude, a)
+    XCTAssertEqual(abs(c), a)
+    XCTAssertEqual(c.description, "-3/2")
+    XCTAssertEqual(c, -3 / 2 as Rational<Int>)
+    XCTAssertTrue(c.isFinite)
+    XCTAssertFalse(c.isProper)
+    XCTAssertTrue(c.reciprocal().isProper)
+
+    let d = -5 / 3 as Rational<Int>
+    XCTAssertEqual(d.magnitude, b)
+    XCTAssertEqual(abs(d), b)
+    XCTAssertEqual(d.description, "-5/3")
+    XCTAssertTrue(d.isFinite)
+    XCTAssertFalse(d.isProper)
+    XCTAssertTrue(d.reciprocal().isProper)
+
+    XCTAssertEqual(c + a, 0)
+    XCTAssertEqual(d + b, 0)
+    XCTAssertEqual(c + c, -3)
+    XCTAssertEqual(d + d, -10 / 3)
+
+    r = c
+    r.negate()
+    XCTAssertEqual(r, a)
+    XCTAssertEqual(r, c.magnitude)
+    XCTAssertEqual(r, abs(c))
+    r = d
+    r.negate()
+    XCTAssertEqual(r, b)
+    XCTAssertEqual(r, d.magnitude)
+    XCTAssertEqual(r, abs(d))
+
+    let e = 42 as Rational<Int>
+    XCTAssertEqual(e.description, "42")
+
     XCTAssertEqual((10 / 9 as Ratio).mixed.whole, 1)
     XCTAssertEqual((10 / 9 as Ratio).mixed.fractional, 1 / 9)
 
     // Test special values.
     let pn = Ratio.nan
+    XCTAssertEqual(pn.description, "nan")
+    XCTAssertTrue(pn.isCanonical)
+
     let pi = Ratio.infinity
+    XCTAssertEqual(pi.description, "inf")
+    XCTAssertTrue(pi.isCanonical)
+
+#if !swift(>=4.1)
+    XCTAssertTrue(Ratio(numerator: 2, denominator: 0).isInfinite)
+    XCTAssertFalse(Ratio(numerator: 2, denominator: 0).isCanonical)
+#endif
+
     let ni = -Ratio.infinity
+    XCTAssertEqual(ni.description, "-inf")
+    XCTAssertTrue(ni.isCanonical)
+
+#if !swift(>=4.1)
+    XCTAssertTrue(Ratio(numerator: -2, denominator: 0).isInfinite)
+    XCTAssertFalse(Ratio(numerator: -2, denominator: 0).isCanonical)
+#endif
+
+    let zero = 0 as Ratio
+    XCTAssertEqual(zero.description, "0")
+    XCTAssertTrue(zero.isCanonical)
+    XCTAssertTrue(zero.isZero)
 
     XCTAssertTrue((pn + pn).isNaN)
     XCTAssertTrue((pn - pn).isNaN)
@@ -35,18 +123,29 @@ class RationalTests: XCTestCase {
 
     XCTAssertEqual(pi + pi, .infinity)
     XCTAssertEqual(ni + ni, -.infinity)
+#if !swift(>=4.1)
     XCTAssertEqual(pi + 0, .infinity)
     XCTAssertEqual(ni + 0, -.infinity)
+#endif
     XCTAssertEqual(pi + 42, .infinity)
     XCTAssertEqual(pi - 42, .infinity)
     XCTAssertEqual(ni + 42, -.infinity)
     XCTAssertEqual(ni - 42, -.infinity)
+#if !swift(>=4.1)
+    XCTAssertEqual(0 + pi, .infinity)
+    XCTAssertEqual(0 + ni, -.infinity)
+#endif
+    XCTAssertEqual(42 + pi, .infinity)
+    XCTAssertEqual(42 - pi, -.infinity)
+    XCTAssertEqual(42 + ni, -.infinity)
+    XCTAssertEqual(42 - ni, .infinity)
 
     XCTAssertTrue((pi + ni).isNaN)
     XCTAssertTrue((ni + pi).isNaN)
     XCTAssertTrue((pi - pi).isNaN)
     XCTAssertTrue((-pi + pi).isNaN)
 
+#if !swift(>=4.1)
     XCTAssertTrue((0 / 0 as Ratio).isNaN)
     XCTAssert((0 / 0 as Ratio) != .nan) // NaN compares unequal to everything.
     XCTAssertTrue((42 / 0 as Ratio).isInfinite)
@@ -55,54 +154,106 @@ class RationalTests: XCTestCase {
     XCTAssert(-42 / 0 as Ratio == -.infinity)
     XCTAssertEqual((42 / 0 as Ratio).description, "inf")
     XCTAssertEqual((-42 / 0 as Ratio).description, "-inf")
+#endif
 
     XCTAssertEqual(pi * pi, .infinity)
     XCTAssertEqual(pi * ni, -.infinity)
     XCTAssertEqual(ni * pi, -.infinity)
     XCTAssertEqual(ni * ni, .infinity)
 
+#if !swift(>=4.1)
     XCTAssertTrue((pi * 0).isNaN)
     XCTAssertTrue((ni * 0).isNaN)
     XCTAssertTrue((0 * pi).isNaN)
     XCTAssertTrue((0 * ni).isNaN)
+#endif
     XCTAssertTrue((pn * pi).isNaN)
     XCTAssertTrue((pi * pn).isNaN)
     XCTAssertTrue((pn * ni).isNaN)
     XCTAssertTrue((ni * pn).isNaN)
+
+    XCTAssertTrue(pn.canonical.isNaN)
   }
 
   func testRationalConversion() {
-    var d: Double, r: Ratio
+    var d: Double, r: Ratio, r8: Rational<Int8>, r64: Rational<Int64>
 
     r = Ratio(42)
     XCTAssertEqual(r, 42)
     XCTAssertEqual(Int(r), 42)
+    XCTAssertEqual(Int(exactly: r)!, 42)
+
+    r = Ratio(exactly: 42)!
+    XCTAssertEqual(r, 42)
+    XCTAssertEqual(Int(r), 42)
+    XCTAssertEqual(Int(exactly: r)!, 42)
+
+    r = Ratio(42 as Double)
+    XCTAssertEqual(r, 42)
+    XCTAssertEqual(Int(r), 42)
+    XCTAssertEqual(Int(exactly: r)!, 42)
+
+    r = Ratio(exactly: 42 as Double)!
+    XCTAssertEqual(r, 42)
+    XCTAssertEqual(Int(r), 42)
+    XCTAssertEqual(Int(exactly: r)!, 42)
 
     r = Ratio(UInt8.max)
     XCTAssertEqual(UInt8(r), .max)
+    XCTAssertEqual(UInt8(exactly: r)!, .max)
+
+    r = Ratio(exactly: UInt8.max)!
+    XCTAssertEqual(UInt8(r), .max)
+    XCTAssertEqual(UInt8(exactly: r)!, .max)
 
     r = Ratio(Int16.min)
     XCTAssertEqual(Int16(r), .min)
     XCTAssertNil(UInt16(exactly: r))
 
-    XCTAssertNil(Ratio(exactly: Int.min))
+    r = Ratio(exactly: Int16.min)!
+    XCTAssertEqual(Int16(r), .min)
+    XCTAssertNil(UInt16(exactly: r))
 
-    r = Ratio(Double.pi)
-    XCTAssertEqual(Double(r), .pi)
-    XCTAssertTrue(r.isCanonical)
+    r8 = Rational(exactly: Int8.max)!
+    XCTAssertEqual(Int8(r8), .max)
+    XCTAssertEqual(Int8(exactly: r8)!, .max)
 
-    r = Ratio(Double.e)
-    XCTAssertEqual(Double(r), .e)
-    XCTAssertTrue(r.isCanonical)
+    r8 = Rational(exactly: Double(Int8.max))!
+    XCTAssertEqual(Int8(r8), .max)
+    XCTAssertEqual(Int8(exactly: r8)!, .max)
 
-    r = Ratio(Double.phi)
-    XCTAssertEqual(Double(r), .phi)
-    XCTAssertTrue(r.isCanonical)
+    XCTAssertNil(Rational<Int8>(exactly: Int16.max))
+    XCTAssertNil(Rational<Int8>(exactly: Double(Int8.max) + 1))
+    XCTAssertNil(Rational<Int8>(exactly: Int8.min))
+    XCTAssertNil(Rational<Int8>(exactly: Double(Int8.min)))
+    XCTAssertNil(Rational<Int8>(exactly: 1 / Double(Int8.min)))
+
+    r64 = Rational(Double.pi)
+    XCTAssertEqual(Double(r64), .pi)
+    XCTAssertTrue(r64.isCanonical)
+
+    r64 = Rational(exactly: Double.pi)!
+    XCTAssertEqual(Double(r64), .pi)
+    XCTAssertTrue(r64.isCanonical)
+
+    XCTAssertNil(Rational<Int32>(exactly: Double.pi))
+
+    r64 = Rational(Double.e)
+    XCTAssertEqual(Double(r64), .e)
+    XCTAssertTrue(r64.isCanonical)
+
+    r64 = Rational(Double.phi)
+    XCTAssertEqual(Double(r64), .phi)
+    XCTAssertTrue(r64.isCanonical)
 
     r = Ratio(Float.pi)
     XCTAssertEqual(Float(r), .pi)
     XCTAssertTrue(r.isCanonical)
 
+    r = Ratio(exactly: Float.pi)!
+    XCTAssertEqual(Float(r), .pi)
+    XCTAssertTrue(r.isCanonical)
+
     r = Ratio(Float.e)
     XCTAssertEqual(Float(r), .e)
     XCTAssertTrue(r.isCanonical)
@@ -113,6 +264,12 @@ class RationalTests: XCTestCase {
 
     d = 0
     XCTAssertEqual(Ratio(d), 0)
+    XCTAssertEqual(Ratio(exactly: d)!, 0)
+    XCTAssertEqual(Double(Ratio(d)), 0)
+
+    d = -0.0
+    XCTAssertEqual(Ratio(d), 0)
+    XCTAssertEqual(Ratio(exactly: d)!, 0)
     XCTAssertEqual(Double(Ratio(d)), 0)
 
     d = 0.ulp
@@ -123,15 +280,50 @@ class RationalTests: XCTestCase {
 
     d = .infinity
     XCTAssertEqual(Ratio(d), .infinity)
+    XCTAssertEqual(Ratio(exactly: d)!, .infinity)
     XCTAssertEqual(Double(Ratio(d)), .infinity)
 
     d = -.infinity
     XCTAssertEqual(Ratio(d), -.infinity)
+    XCTAssertEqual(Ratio(exactly: d)!, -.infinity)
     XCTAssertEqual(Double(Ratio(d)), -.infinity)
 
     d = .nan
     XCTAssertTrue(Ratio(d).isNaN)
+    XCTAssertTrue(Ratio(exactly: d)!.isNaN)
     XCTAssertTrue(Double(Ratio(d)).isNaN)
+
+    let a = Rational<Int8>(0.125)
+    XCTAssertEqual(Float(a), 0.125)
+    XCTAssertEqual(Double(a), 0.125)
+
+    let aa = Rational<Int8>(8)
+    XCTAssertEqual(Float(aa), 8)
+    XCTAssertEqual(Double(aa), 8)
+
+    let b = Rational<Int16>(0.0625)
+    XCTAssertEqual(Float(b), 0.0625)
+    XCTAssertEqual(Double(b), 0.0625)
+
+    let bb = Rational<Int16>(16)
+    XCTAssertEqual(Float(bb), 16)
+    XCTAssertEqual(Double(bb), 16)
+
+    let c = Rational<Int32>(0.03125)
+    XCTAssertEqual(Float(c), 0.03125)
+    XCTAssertEqual(Double(c), 0.03125)
+
+    let cc = Rational<Int32>(32)
+    XCTAssertEqual(Float(cc), 32)
+    XCTAssertEqual(Double(cc), 32)
+
+    r64 = Rational<Int64>(0.015625)
+    XCTAssertEqual(Float(r64), 0.015625)
+    XCTAssertEqual(Double(r64), 0.015625)
+
+    r64 = Rational<Int64>(64)
+    XCTAssertEqual(Float(r64), 64)
+    XCTAssertEqual(Double(r64), 64)
   }
 
   func testRationalComparison() {
@@ -139,19 +331,36 @@ class RationalTests: XCTestCase {
     let b = 1 / 4 as Ratio
     XCTAssert(a != b)
     XCTAssert(a == a)
-    XCTAssert(b == b)
-    XCTAssert(a > b)
-    XCTAssert(b < a)
+    XCTAssert(a <= a)
+    XCTAssert(a >= a)
     XCTAssert(!(a < a))
+    XCTAssert(b == b)
+    XCTAssert(b <= b)
+    XCTAssert(b >= b)
+    XCTAssert(a > b)
+    XCTAssert(a >= b)
+    XCTAssert(b < a)
+    XCTAssert(b <= a)
 
     let c = -1 / 2 as Ratio
     let d = -1 / 4 as Ratio
     XCTAssert(c != d)
     XCTAssert(c == c)
-    XCTAssert(d == d)
-    XCTAssert(c < d)
-    XCTAssert(d > c)
+    XCTAssert(c <= c)
+    XCTAssert(c >= c)
     XCTAssert(!(c < c))
+    XCTAssert(d == d)
+    XCTAssert(d <= d)
+    XCTAssert(d >= d)
+    XCTAssert(c < d)
+    XCTAssert(c <= d)
+    XCTAssert(d > c)
+    XCTAssert(d >= c)
+
+    XCTAssert(c < a)
+    XCTAssert(!(a < c))
+    XCTAssert(c <= a)
+    XCTAssert(!(a <= c))
 
     let e = Ratio(numerator: 0, denominator: 1)
     let f = Ratio(numerator: 0, denominator: 2)
@@ -178,13 +387,21 @@ class RationalTests: XCTestCase {
     XCTAssert(!(j < k))
     XCTAssert(j != l)
     XCTAssert(j > l)
+    XCTAssert(j > 42)
+    XCTAssert(j >= 42)
     XCTAssert(!(j < l))
+    XCTAssert(!(j < 42))
+    XCTAssert(!(j <= 42))
     XCTAssert(j != m)
     XCTAssert(j > m)
     XCTAssert(!(j < m))
     XCTAssert(l != j)
     XCTAssert(l < j)
+    XCTAssert(l < 42)
+    XCTAssert(l <= 42)
     XCTAssert(!(l > j))
+    XCTAssert(!(l > 42))
+    XCTAssert(!(l >= 42))
     XCTAssert(l != k)
     XCTAssert(l < k)
     XCTAssert(!(l > k))
@@ -218,6 +435,11 @@ class RationalTests: XCTestCase {
     XCTAssertEqual(a.rounded(.toNearestOrEven), 3)
     XCTAssertEqual(a.rounded(.toNearestOrAwayFromZero), 3)
 
+    XCTAssertEqual(round(a), 3)
+    XCTAssertEqual(ceil(a), 4)
+    XCTAssertEqual(floor(a), 3)
+    XCTAssertEqual(trunc(a), 3)
+
     let b = -a
     XCTAssertEqual(b, -10 / 3)
     XCTAssertEqual(b.rounded(), -3)
@@ -228,6 +450,11 @@ class RationalTests: XCTestCase {
     XCTAssertEqual(b.rounded(.toNearestOrEven), -3)
     XCTAssertEqual(b.rounded(.toNearestOrAwayFromZero), -3)
 
+    XCTAssertEqual(round(b), -3)
+    XCTAssertEqual(ceil(b), -3)
+    XCTAssertEqual(floor(b), -4)
+    XCTAssertEqual(trunc(b), -3)
+
     let c = 5 / 2 as Ratio
     XCTAssertEqual(c.rounded(), 3)
     XCTAssertEqual(c.rounded(.up), 3)
@@ -237,6 +464,11 @@ class RationalTests: XCTestCase {
     XCTAssertEqual(c.rounded(.toNearestOrEven), 2)
     XCTAssertEqual(c.rounded(.toNearestOrAwayFromZero), 3)
 
+    XCTAssertEqual(round(c), 3)
+    XCTAssertEqual(ceil(c), 3)
+    XCTAssertEqual(floor(c), 2)
+    XCTAssertEqual(trunc(c), 2)
+
     let d = -c
     XCTAssertEqual(d, -5 / 2)
     XCTAssertEqual(d.rounded(), -3)
@@ -247,6 +479,11 @@ class RationalTests: XCTestCase {
     XCTAssertEqual(d.rounded(.toNearestOrEven), -2)
     XCTAssertEqual(d.rounded(.toNearestOrAwayFromZero), -3)
 
+    XCTAssertEqual(round(d), -3)
+    XCTAssertEqual(ceil(d), -2)
+    XCTAssertEqual(floor(d), -3)
+    XCTAssertEqual(trunc(d), -2)
+
     let e = 1 / 9 as Ratio
     XCTAssertEqual(e.rounded(), 0)
     XCTAssertEqual(e.rounded(.up), 1)
@@ -255,6 +492,63 @@ class RationalTests: XCTestCase {
     XCTAssertEqual(e.rounded(.awayFromZero), 1)
     XCTAssertEqual(e.rounded(.toNearestOrEven), 0)
     XCTAssertEqual(e.rounded(.toNearestOrAwayFromZero), 0)
+
+    XCTAssertEqual(round(e), 0)
+    XCTAssertEqual(ceil(e), 1)
+    XCTAssertEqual(floor(e), 0)
+    XCTAssertEqual(trunc(e), 0)
+
+    let f = Ratio.infinity
+    XCTAssertEqual(f.rounded(), f)
+
+    let g = -Ratio.infinity
+    XCTAssertEqual(g.rounded(), g)
+
+    let h = Ratio.nan
+    XCTAssertTrue(h.rounded().isNaN)
+  }
+
+  func testRationalStride() {
+    let a = zip(stride(from: 1/3 as Ratio, to: 3, by: 1/3), 1...8)
+    for i in a {
+      XCTAssertEqual(i.0 * 3, Ratio(i.1))
+    }
+
+    let b = zip(stride(from: 1/3 as Ratio, through: 3, by: 1/3), 1...9)
+    for i in b {
+      XCTAssertEqual(i.0 * 3, Ratio(i.1))
+    }
+    XCTAssertEqual(Int(b.map { $0.0 }.last!), 3)
+
+    let c = zip(stride(from: 3 as Ratio, to: 1/3, by: -1/3), (2...9).reversed())
+    for i in c {
+      XCTAssertEqual(i.0 * 3, Ratio(i.1))
+    }
+
+    let d =
+      zip(stride(from: 3 as Ratio, through: 1/3, by: -1/3), (1...9).reversed())
+    for i in d {
+      XCTAssertEqual(i.0 * 3, Ratio(i.1))
+    }
+    XCTAssertEqual(Int(d.map { $0.0 }.last! * 3), 1)
+
+    let e = zip(stride(from: -3 as Ratio, to: -1/3, by: 1/3), -9...(-1))
+    for i in e {
+      XCTAssertEqual(i.0 * 3, Ratio(i.1))
+    }
+
+    let f = zip(stride(from: -3 as Ratio, through: -1/3, by: 1/3), -9...(-1))
+    for i in f {
+      XCTAssertEqual(i.0 * 3, Ratio(i.1))
+    }
+    XCTAssertEqual(Int(f.map { $0.0 }.last! * 3), -1)
+
+    let values: [Ratio] = [6/4, 5/3, 4/6, 3/5, -3/5, -4/6, -5/3, -6/4]
+    for i in values {
+      for j in values {
+        XCTAssertEqual(i.advanced(by: i.distance(to: j)), j)
+      }
+    }
   }
 
   static var allTests = [
@@ -262,5 +556,6 @@ class RationalTests: XCTestCase {
     ("testRationalConversion", testRationalConversion),
     ("testRationalComparison", testRationalComparison),
     ("testRationalRounding", testRationalRounding),
+    ("testRationalStride", testRationalStride),
   ]
 }

Tests/NumericAnnexTests/RealTests.swift

@@ -0,0 +1,778 @@
+import XCTest
+@testable import NumericAnnex
+
+/// A type that wraps `Double` and implements the methods required for
+/// conformance to `Real`.
+struct MockReal {
+  var _value: Double
+
+  init(_ value: Double) {
+    self._value = value
+  }
+}
+
+extension MockReal : ExpressibleByIntegerLiteral {
+  init(integerLiteral value: Double.IntegerLiteralType) {
+    self = MockReal(Double(integerLiteral: value))
+  }
+}
+
+extension MockReal : Hashable {
+  var hashValue: Int {
+    return _value.hashValue
+  }
+}
+
+extension MockReal : Strideable {
+  func distance(to other: MockReal) -> MockReal {
+    return MockReal(_value.distance(to: other._value))
+  }
+
+  func advanced(by n: MockReal) -> MockReal {
+    return MockReal(_value.advanced(by: n._value))
+  }
+}
+
+#if !swift(>=4.1)
+extension MockReal : _Strideable { }
+#endif
+
+extension MockReal : Numeric {
+  var magnitude: MockReal {
+    return MockReal(_value.magnitude)
+  }
+
+  init?<T>(exactly source: T) where T : BinaryInteger {
+    guard let d = Double(exactly: source) else { return nil }
+    self = MockReal(d)
+  }
+
+  static func + (lhs: MockReal, rhs: MockReal) -> MockReal {
+    return MockReal(lhs._value + rhs._value)
+  }
+
+  static func += (lhs: inout MockReal, rhs: MockReal) {
+    lhs._value += rhs._value
+  }
+
+  static func - (lhs: MockReal, rhs: MockReal) -> MockReal {
+    return MockReal(lhs._value - rhs._value)
+  }
+
+  static func -= (lhs: inout MockReal, rhs: MockReal) {
+    lhs._value -= rhs._value
+  }
+
+  static func * (lhs: MockReal, rhs: MockReal) -> MockReal {
+    return MockReal(lhs._value * rhs._value)
+  }
+
+  static func *= (lhs: inout MockReal, rhs: MockReal) {
+    lhs._value *= rhs._value
+  }
+}
+
+extension MockReal : FloatingPoint {
+  init<T : BinaryInteger>(_ value: T) {
+#if false
+    self = MockReal(Double(value))
+#endif
+    fatalError()
+  }
+
+  init(_ value: Int) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: Int8) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: Int16) {
+    self = MockReal(Double(value))
+  }
+  
+  init(_ value: Int32) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: Int64) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: UInt) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: UInt8) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: UInt16) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: UInt32) {
+    self = MockReal(Double(value))
+  }
+
+  init(_ value: UInt64) {
+    self = MockReal(Double(value))
+  }
+
+  init(sign: FloatingPointSign, exponent: Int, significand: MockReal) {
+    self = MockReal(
+      Double(sign: sign, exponent: exponent, significand: significand._value)
+    )
+  }
+
+  init(signOf s: MockReal, magnitudeOf m: MockReal) {
+    self = MockReal(Double(signOf: s._value, magnitudeOf: m._value))
+  }
+
+  static var greatestFiniteMagnitude = MockReal(.greatestFiniteMagnitude)
+
+  static var infinity = MockReal(.infinity)
+
+  static var leastNonzeroMagnitude = MockReal(.leastNonzeroMagnitude)
+
+  static var leastNormalMagnitude = MockReal(.leastNormalMagnitude)
+
+  static var nan = MockReal(.nan)
+
+  static var pi = MockReal(.pi)
+
+  static var radix = Double.radix
+
+  static var signalingNaN = MockReal(.signalingNaN)
+
+  static func / (lhs: MockReal, rhs: MockReal) -> MockReal {
+    return MockReal(lhs._value / rhs._value)
+  }
+
+  static func /= (lhs: inout MockReal, rhs: MockReal) {
+    lhs._value /= rhs._value
+  }
+
+  static prefix func - (x: MockReal) -> MockReal {
+    return MockReal(-x._value)
+  }
+
+  var exponent: Int {
+    return _value.exponent
+  }
+
+  var isCanonical: Bool {
+    return _value.isCanonical
+  }
+
+  var isFinite: Bool {
+    return _value.isFinite
+  }
+
+  var isInfinite: Bool {
+    return _value.isInfinite
+  }
+
+  var isNaN: Bool {
+    return _value.isNaN
+  }
+
+  var isNormal: Bool {
+    return _value.isNormal
+  }
+
+  var isSignalingNaN: Bool {
+    return _value.isSignalingNaN
+  }
+
+  var isSubnormal: Bool {
+    return _value.isSubnormal
+  }
+
+  var isZero: Bool {
+    return _value.isZero
+  }
+
+  var nextUp: MockReal {
+    return MockReal(_value.nextUp)
+  }
+
+  var sign: FloatingPointSign {
+    return _value.sign
+  }
+
+  var significand: MockReal {
+    return MockReal(_value.significand)
+  }
+
+  var ulp: MockReal {
+    return MockReal(_value.ulp)
+  }
+
+  mutating func addProduct(_ lhs: MockReal, _ rhs: MockReal) {
+    _value.addProduct(lhs._value, rhs._value)
+  }
+
+  mutating func formRemainder(dividingBy other: MockReal) {
+    _value.formRemainder(dividingBy: other._value)
+  }
+
+  mutating func formSquareRoot() {
+    _value.formSquareRoot()
+  }
+
+  mutating func formTruncatingRemainder(dividingBy other: MockReal) {
+    _value.formTruncatingRemainder(dividingBy: other._value)
+  }
+
+  func isEqual(to other: MockReal) -> Bool {
+    return _value.isEqual(to: other._value)
+  }
+
+  func isLess(than other: MockReal) -> Bool {
+    return _value.isLess(than: other._value)
+  }
+
+  func isLessThanOrEqualTo(_ other: MockReal) -> Bool {
+    return _value.isLessThanOrEqualTo(other._value)
+  }
+
+  func isTotallyOrdered(belowOrEqualTo other: MockReal) -> Bool {
+    return _value.isTotallyOrdered(belowOrEqualTo: other._value)
+  }
+
+  mutating func negate() {
+    _value.negate()
+  }
+
+  mutating func round(_ rule: FloatingPointRoundingRule) {
+    _value.round(rule)
+  }
+}
+
+extension MockReal : Real {
+  static func ** (lhs: MockReal, rhs: MockReal) -> MockReal {
+    return MockReal(lhs._value ** rhs._value)
+  }
+
+  static func **= (lhs: inout MockReal, rhs: MockReal) {
+    lhs = lhs ** rhs
+  }
+
+  func naturalExponential() -> MockReal {
+    return MockReal(_value.naturalExponential())
+  }
+
+  func naturalLogarithm() -> MockReal {
+    return MockReal(_value.naturalLogarithm())
+  }
+
+  func sine() -> MockReal {
+    return MockReal(_value.sine())
+  }
+
+  func cosine() -> MockReal {
+    return MockReal(_value.cosine())
+  }
+
+  func inverseSine() -> MockReal {
+    return MockReal(_value.inverseSine())
+  }
+
+  func inverseCosine() -> MockReal {
+    return MockReal(_value.inverseCosine())
+  }
+
+  func inverseTangent() -> MockReal {
+    return MockReal(_value.inverseTangent())
+  }
+
+  func hyperbolicSine() -> MockReal {
+    return MockReal(_value.hyperbolicSine())
+  }
+
+  func hyperbolicCosine() -> MockReal {
+    return MockReal(_value.hyperbolicCosine())
+  }
+
+  func inverseHyperbolicSine() -> MockReal {
+    return MockReal(_value.inverseHyperbolicSine())
+  }
+
+  func inverseHyperbolicCosine() -> MockReal {
+    return MockReal(_value.inverseHyperbolicCosine())
+  }
+
+  func inverseHyperbolicTangent() -> MockReal {
+    return MockReal(_value.inverseHyperbolicTangent())
+  }
+
+  func error() -> MockReal {
+    return MockReal(_value.error())
+  }
+
+  func complementaryError() -> MockReal {
+    return MockReal(_value.complementaryError())
+  }
+
+  func gamma() -> MockReal {
+    return MockReal(_value.gamma())
+  }
+
+  func logarithmicGamma() -> MockReal {
+    return MockReal(_value.logarithmicGamma())
+  }
+}
+
+class RealTests: XCTestCase {
+  func testTranscendentalConstants() {
+    XCTAssertEqual(Double.pi, MockReal.pi._value)
+    XCTAssertEqual(Double.e, MockReal.e._value,
+                   accuracy: abs(Double.e * 0.000_000_000_001))
+    XCTAssertEqual(Double.phi, MockReal.phi._value,
+                   accuracy: abs(Double.phi * 0.000_000_000_001))
+  }
+
+  func testExp2() {
+    for i in -42...42 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.exp2(d), MockReal.exp2(m)._value,
+                     accuracy: abs(Double.exp2(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testExp10() {
+    for i in -42...42 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.exp10(d), MockReal.exp10(m)._value,
+                     accuracy: abs(Double.exp10(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testExpm1() {
+    for i in -42...42 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.expm1(d), MockReal.expm1(m)._value,
+                     accuracy: abs(Double.expm1(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testLog2() {
+    for i in 1...85 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.log2(d), MockReal.log2(m)._value,
+                     accuracy: abs(Double.log2(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testLog10() {
+    for i in 1...85 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.log10(d), MockReal.log10(m)._value,
+                     accuracy: abs(Double.log10(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testLog1p() {
+    for i in 1...85 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.log1p(d), MockReal.log1p(m)._value,
+                     accuracy: abs(Double.log1p(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testCbrt() {
+    for i in -42...42 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.cbrt(d), MockReal.cbrt(m)._value,
+                     accuracy: abs(Double.cbrt(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testTan() {
+    for i in -42...42 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.tan(d), MockReal.tan(m)._value,
+                     accuracy: abs(Double.tan(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testTanh() {
+    for i in -42...42 {
+      let d = Double(i)
+      let m = MockReal(d)
+      XCTAssertEqual(Double.tanh(d), MockReal.tanh(m)._value,
+                     accuracy: abs(Double.tanh(d) * 0.000_000_000_001))
+    }
+  }
+
+  func testHypot() {
+    for i in -42...42 {
+      let dx = Double(i)
+      let mx = MockReal(dx)
+      for j in -42...42 {
+        let dy = Double(j)
+        let my = MockReal(dy)
+        XCTAssertEqual(Double.hypot(dx, dy), MockReal.hypot(mx, my)._value,
+                       accuracy: abs(Double.hypot(dx, dy) * 0.000_000_000_001))
+      }
+    }
+    // Test special values.
+    XCTAssertEqual(MockReal.hypot(.infinity, .nan), .infinity)
+    XCTAssertEqual(MockReal.hypot(.infinity, .infinity), .infinity)
+    XCTAssertEqual(MockReal.hypot(.infinity, 42), .infinity)
+    XCTAssertEqual(MockReal.hypot(42, -.infinity), .infinity)
+    XCTAssertEqual(MockReal.hypot(-.infinity, -.infinity), .infinity)
+    XCTAssertEqual(MockReal.hypot(.nan, -.infinity), .infinity)
+    XCTAssertTrue(MockReal.hypot(.nan, .nan).isNaN)
+    XCTAssertTrue(MockReal.hypot(.nan, 42).isNaN)
+    XCTAssertTrue(MockReal.hypot(42, .nan).isNaN)
+  }
+
+  func testAtan2() {
+    for i in -42...42 {
+      let dx = Double(i)
+      let mx = MockReal(dx)
+      for j in -42...42 {
+        let dy = Double(j)
+        let my = MockReal(dy)
+        XCTAssertEqual(Double.atan2(dy, dx), MockReal.atan2(my, mx)._value,
+                       accuracy: abs(Double.atan2(dy, dx) * 0.000_000_000_001))
+      }
+    }
+    // Test special values.
+    XCTAssertTrue(MockReal.atan2(0, 0).isZero)
+    XCTAssertTrue(MockReal.atan2(0, 0).sign == .plus)
+    XCTAssertTrue(MockReal.atan2(MockReal(-0.0), 0).isZero)
+    XCTAssertTrue(MockReal.atan2(MockReal(-0.0), 0).sign == .minus)
+    XCTAssertEqual(MockReal.atan2(0, MockReal(-0.0)), .pi)
+    XCTAssertEqual(MockReal.atan2(MockReal(-0.0), MockReal(-0.0)), -.pi)
+    XCTAssertEqual(MockReal.atan2(0, -42), .pi)
+    XCTAssertEqual(MockReal.atan2(MockReal(-0.0), -42), -.pi)
+    XCTAssertTrue(MockReal.atan2(0, 42).isZero)
+    XCTAssertTrue(MockReal.atan2(0, 42).sign == .plus)
+    XCTAssertTrue(MockReal.atan2(MockReal(-0.0), 42).isZero)
+    XCTAssertTrue(MockReal.atan2(MockReal(-0.0), 42).sign == .minus)
+    XCTAssertEqual(MockReal.atan2(.infinity, 42), .pi / 2)
+    XCTAssertEqual(MockReal.atan2(-.infinity, 42), -.pi / 2)
+    XCTAssertEqual(MockReal.atan2(.infinity, -.infinity), 3 * .pi / 4)
+    XCTAssertEqual(MockReal.atan2(-.infinity, -.infinity), -3 * .pi / 4)
+    XCTAssertEqual(MockReal.atan2(.infinity, .infinity), .pi / 4)
+    XCTAssertEqual(MockReal.atan2(-.infinity, .infinity), -.pi / 4)
+    XCTAssertEqual(MockReal.atan2(-42, 0), -.pi / 2)
+    XCTAssertEqual(MockReal.atan2(-42, MockReal(-0.0)), -.pi / 2)
+    XCTAssertEqual(MockReal.atan2(42, 0), .pi / 2)
+    XCTAssertEqual(MockReal.atan2(42, MockReal(-0.0)), .pi / 2)
+    XCTAssertEqual(MockReal.atan2(42, -.infinity), .pi)
+    XCTAssertEqual(MockReal.atan2(-42, -.infinity), -.pi)
+    XCTAssertTrue(MockReal.atan2(42, .infinity).isZero)
+    XCTAssertTrue(MockReal.atan2(42, .infinity).sign == .plus)
+    XCTAssertTrue(MockReal.atan2(-42, .infinity).isZero)
+    XCTAssertTrue(MockReal.atan2(-42, .infinity).sign == .minus)
+    XCTAssertTrue(MockReal.atan2(.nan, .nan).isNaN)
+    XCTAssertTrue(MockReal.atan2(.nan, 42).isNaN)
+    XCTAssertTrue(MockReal.atan2(42, .nan).isNaN)
+  }
+
+  func testFloat() {
+    XCTAssertEqual(((-2 as Float) ** -3), -0.125)
+    XCTAssertTrue(((0 as Float) ** 42).isZero)
+    XCTAssertTrue(((0 as Float) ** 42).sign == .plus)
+    XCTAssertTrue(((-0.0 as Float) ** 42).isZero)
+    XCTAssertTrue(((-0.0 as Float) ** 42).sign == .plus)
+    XCTAssertTrue(((0 as Float) ** 43).isZero)
+    XCTAssertTrue(((0 as Float) ** 43).sign == .plus)
+    XCTAssertTrue(((-0.0 as Float) ** 43).isZero)
+    XCTAssertTrue(((-0.0 as Float) ** 43).sign == .minus)
+
+    XCTAssertEqual(Float.cbrt(64), 4)
+    XCTAssertEqual(Float.cbrt(-64), -4)
+    XCTAssertEqual(Float.cbrt(27), 3)
+    XCTAssertEqual(Float.cbrt(-27), -3)
+    XCTAssertEqual(Float.cbrt(8), 2)
+    XCTAssertEqual(Float.cbrt(-8), -2)
+    XCTAssertEqual(Float.cbrt(-0.125), -0.5)
+    XCTAssertEqual(Float.cbrt(0.125), 0.5)
+    XCTAssertEqual(Float.cbrt(-0.015625), -0.25)
+    XCTAssertEqual(Float.cbrt(0.015625), 0.25)
+
+    XCTAssertEqual(Float.exp2(0), 1)
+    XCTAssertEqual(Float.exp2(2), 4)
+    XCTAssertEqual(Float.exp2(4), 16)
+
+    XCTAssertEqual(Float.exp10(0), 1)
+    XCTAssertEqual(Float.exp10(2), 100)
+    XCTAssertEqual(Float.exp10(4), 10000)
+
+    XCTAssertEqual(Float.expm1(16), .exp(16) - 1)
+    XCTAssertTrue(Float.expm1(0).isZero)
+    XCTAssertEqual(Float.expm1(-.infinity), -1)
+    XCTAssertEqual(Float.expm1(.infinity), .infinity)
+
+    XCTAssertEqual(Float.log2(16), 4)
+    XCTAssertEqual(Float.log2(4), 2)
+    XCTAssertEqual(Float.log2(1), 0)
+
+    XCTAssertEqual(Float.log10(10000), 4)
+    XCTAssertEqual(Float.log10(100), 2)
+    XCTAssertEqual(Float.log10(1), 0)
+
+    XCTAssertEqual(Float.log1p(16777216), .log(16777215))
+    XCTAssertTrue(Float.log1p(0).isZero)
+    XCTAssertEqual(Float.log1p(-1), -.infinity)
+    XCTAssertEqual(Float.log1p(.infinity), .infinity)
+
+    XCTAssertEqual(Float.sin(.pi / 2), 1)
+    XCTAssertTrue(Float.sin(0).isZero)
+    XCTAssertTrue(Float.sin(0).sign == .plus)
+    XCTAssertTrue(Float.sin(-0.0).isZero)
+    XCTAssertTrue(Float.sin(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.cos(.pi / 2), 0, accuracy: 0.000_001)
+    XCTAssertEqual(Float.cos(0), 1)
+    XCTAssertEqual(Float.cos(-0.0), 1)
+
+    XCTAssertEqual(Float.tan(.pi / 4), 1, accuracy: 0.000_001)
+    XCTAssertEqual(Float.tan(3 * .pi / 4), -1, accuracy: 0.000_001)
+    XCTAssertTrue(Float.tan(0).isZero)
+    XCTAssertTrue(Float.tan(0).sign == .plus)
+    XCTAssertTrue(Float.tan(-0.0).isZero)
+    XCTAssertTrue(Float.tan(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.asin(1), .pi / 2, accuracy: 0.000_001)
+    XCTAssertTrue(Float.asin(0).isZero)
+    XCTAssertTrue(Float.asin(0).sign == .plus)
+    XCTAssertTrue(Float.asin(-0.0).isZero)
+    XCTAssertTrue(Float.asin(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.acos(0), .pi / 2, accuracy: 0.000_001)
+    XCTAssertEqual(Float.acos(1), 0)
+    XCTAssertEqual(Float.acos(-1), .pi, accuracy: 0.000_001)
+
+    XCTAssertEqual(Float.atan(1), .pi / 4, accuracy: 0.000_001)
+    XCTAssertEqual(Float.atan(-1), -.pi / 4, accuracy: 0.000_001)
+    XCTAssertTrue(Float.atan(0).isZero)
+    XCTAssertTrue(Float.atan(0).sign == .plus)
+    XCTAssertTrue(Float.atan(-0.0).isZero)
+    XCTAssertTrue(Float.atan(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.sinh(1), 1.175201194, accuracy: 0.000_001)
+    XCTAssertTrue(Float.sinh(0).isZero)
+    XCTAssertTrue(Float.sinh(0).sign == .plus)
+    XCTAssertTrue(Float.sinh(-0.0).isZero)
+    XCTAssertTrue(Float.sinh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.cosh(1), 1.543080635, accuracy: 0.000_001)
+    XCTAssertEqual(Float.cosh(0), 1)
+    XCTAssertEqual(Float.cosh(-0.0), 1)
+
+    XCTAssertEqual(Float.tanh(1), 0.761594156, accuracy: 0.000_001)
+    XCTAssertTrue(Float.tanh(0).isZero)
+    XCTAssertTrue(Float.tanh(0).sign == .plus)
+    XCTAssertTrue(Float.tanh(-0.0).isZero)
+    XCTAssertTrue(Float.tanh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.asinh(1), 0.881373587, accuracy: 0.000_001)
+    XCTAssertTrue(Float.asinh(0).isZero)
+    XCTAssertTrue(Float.asinh(0).sign == .plus)
+    XCTAssertTrue(Float.asinh(-0.0).isZero)
+    XCTAssertTrue(Float.asinh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.acosh(2), 1.316957897, accuracy: 0.000_001)
+    XCTAssertEqual(Float.acosh(1), 0)
+    XCTAssertEqual(Float.acosh(.infinity), .infinity)
+
+    XCTAssertEqual(Float.atanh(0.5), 0.549306144, accuracy: 0.000_001)
+    XCTAssertTrue(Float.atanh(0).isZero)
+    XCTAssertTrue(Float.atanh(0).sign == .plus)
+    XCTAssertTrue(Float.atanh(-0.0).isZero)
+    XCTAssertTrue(Float.atanh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Float.erf(1), 0.842700793, accuracy: 0.000_001)
+    XCTAssertTrue(Float.erf(0).isZero)
+    XCTAssertTrue(Float.erf(0).sign == .plus)
+    XCTAssertTrue(Float.erf(-0.0).isZero)
+    XCTAssertTrue(Float.erf(-0.0).sign == .minus)
+    XCTAssertEqual(Float.erf(.infinity), 1)
+    XCTAssertEqual(Float.erf(-.infinity), -1)
+
+    XCTAssertEqual(Float.erfc(1), 0.157299207, accuracy: 0.000_001)
+    XCTAssertEqual(Float.erfc(0), 1)
+    XCTAssertEqual(Float.erfc(.infinity), 0)
+    XCTAssertEqual(Float.erfc(-.infinity), 2)
+
+    XCTAssertEqual(Float.tgamma(10), Float((2...9).reduce(1, *)), accuracy: 0.1)
+    XCTAssertEqual(Float.tgamma(1), 1)
+    XCTAssertEqual(Float.tgamma(.infinity), .infinity)
+
+    XCTAssertEqual(Float.lgamma(10), Float.log(Float((2...9).reduce(1, *))),
+                   accuracy: 0.000_001)
+    XCTAssertEqual(Float.lgamma(1), 0)
+    XCTAssertEqual(Float.lgamma(.infinity), .infinity)
+  }
+
+  func testDouble() {
+    XCTAssertEqual(((-2 as Double) ** -3), -0.125)
+    XCTAssertTrue(((0 as Double) ** 42).isZero)
+    XCTAssertTrue(((0 as Double) ** 42).sign == .plus)
+    XCTAssertTrue(((-0.0 as Double) ** 42).isZero)
+    XCTAssertTrue(((-0.0 as Double) ** 42).sign == .plus)
+    XCTAssertTrue(((0 as Double) ** 43).isZero)
+    XCTAssertTrue(((0 as Double) ** 43).sign == .plus)
+    XCTAssertTrue(((-0.0 as Double) ** 43).isZero)
+    XCTAssertTrue(((-0.0 as Double) ** 43).sign == .minus)
+
+    XCTAssertEqual(Double.cbrt(64), 4)
+    XCTAssertEqual(Double.cbrt(-64), -4)
+    XCTAssertEqual(Double.cbrt(27), 3, accuracy: 0.000_001)
+    XCTAssertEqual(Double.cbrt(-27), -3, accuracy: 0.000_001)
+    XCTAssertEqual(Double.cbrt(8), 2)
+    XCTAssertEqual(Double.cbrt(-8), -2)
+    XCTAssertEqual(Double.cbrt(-0.125), -0.5, accuracy: 0.000_001)
+    XCTAssertEqual(Double.cbrt(0.125), 0.5, accuracy: 0.000_001)
+    XCTAssertEqual(Double.cbrt(-0.015625), -0.25, accuracy: 0.000_001)
+    XCTAssertEqual(Double.cbrt(0.015625), 0.25, accuracy: 0.000_001)
+
+    XCTAssertEqual(Double.exp2(0), 1)
+    XCTAssertEqual(Double.exp2(2), 4)
+    XCTAssertEqual(Double.exp2(4), 16)
+
+    XCTAssertEqual(Double.exp10(0), 1)
+    XCTAssertEqual(Double.exp10(2), 100)
+    XCTAssertEqual(Double.exp10(4), 10000, accuracy: 0.000_001)
+
+    XCTAssertEqual(Double.expm1(16), .exp(16) - 1)
+    XCTAssertTrue(Double.expm1(0).isZero)
+    XCTAssertEqual(Double.expm1(-.infinity), -1)
+    XCTAssertEqual(Double.expm1(.infinity), .infinity)
+
+    XCTAssertEqual(Double.log2(16), 4)
+    XCTAssertEqual(Double.log2(4), 2)
+    XCTAssertEqual(Double.log2(1), 0)
+
+    XCTAssertEqual(Double.log10(10000), 4)
+    XCTAssertEqual(Double.log10(100), 2)
+    XCTAssertEqual(Double.log10(1), 0)
+
+    XCTAssertEqual(Double.log1p(16777216), .log(16777215), accuracy: 0.000_001)
+    XCTAssertTrue(Double.log1p(0).isZero)
+    XCTAssertEqual(Double.log1p(-1), -.infinity)
+    XCTAssertEqual(Double.log1p(.infinity), .infinity)
+
+    XCTAssertEqual(Double.sin(.pi / 2), 1)
+    XCTAssertTrue(Double.sin(0).isZero)
+    XCTAssertTrue(Double.sin(0).sign == .plus)
+    XCTAssertTrue(Double.sin(-0.0).isZero)
+    XCTAssertTrue(Double.sin(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.cos(.pi / 2), 0, accuracy: 0.000_000_000_001)
+    XCTAssertEqual(Double.cos(0), 1)
+    XCTAssertEqual(Double.cos(-0.0), 1)
+
+    XCTAssertEqual(Double.tan(.pi / 4), 1, accuracy: 0.000_000_000_001)
+    XCTAssertEqual(Double.tan(3 * .pi / 4), -1, accuracy: 0.000_000_000_001)
+    XCTAssertTrue(Double.tan(0).isZero)
+    XCTAssertTrue(Double.tan(0).sign == .plus)
+    XCTAssertTrue(Double.tan(-0.0).isZero)
+    XCTAssertTrue(Double.tan(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.asin(1), .pi / 2)
+    XCTAssertTrue(Double.asin(0).isZero)
+    XCTAssertTrue(Double.asin(0).sign == .plus)
+    XCTAssertTrue(Double.asin(-0.0).isZero)
+    XCTAssertTrue(Double.asin(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.acos(0), .pi / 2, accuracy: 0.000_000_000_001)
+    XCTAssertEqual(Double.acos(1), 0)
+    XCTAssertEqual(Double.acos(-1), .pi)
+
+    XCTAssertEqual(Double.atan(1), .pi / 4, accuracy: 0.000_000_000_001)
+    XCTAssertEqual(Double.atan(-1), -.pi / 4, accuracy: 0.000_000_000_001)
+    XCTAssertTrue(Double.atan(0).isZero)
+    XCTAssertTrue(Double.atan(0).sign == .plus)
+    XCTAssertTrue(Double.atan(-0.0).isZero)
+    XCTAssertTrue(Double.atan(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.sinh(1), 1.175201193644, accuracy: 0.000_000_000_001)
+    XCTAssertTrue(Double.sinh(0).isZero)
+    XCTAssertTrue(Double.sinh(0).sign == .plus)
+    XCTAssertTrue(Double.sinh(-0.0).isZero)
+    XCTAssertTrue(Double.sinh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.cosh(1), 1.543080634815, accuracy: 0.000_000_000_001)
+    XCTAssertEqual(Double.cosh(0), 1)
+    XCTAssertEqual(Double.cosh(-0.0), 1)
+
+    XCTAssertEqual(Double.tanh(1), 0.761594155956, accuracy: 0.000_000_000_001)
+    XCTAssertTrue(Double.tanh(0).isZero)
+    XCTAssertTrue(Double.tanh(0).sign == .plus)
+    XCTAssertTrue(Double.tanh(-0.0).isZero)
+    XCTAssertTrue(Double.tanh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.asinh(1), 0.881373587020, accuracy: 0.000_000_000_001)
+    XCTAssertTrue(Double.asinh(0).isZero)
+    XCTAssertTrue(Double.asinh(0).sign == .plus)
+    XCTAssertTrue(Double.asinh(-0.0).isZero)
+    XCTAssertTrue(Double.asinh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.acosh(2), 1.316957896925, accuracy: 0.000_000_000_001)
+    XCTAssertEqual(Double.acosh(1), 0)
+    XCTAssertEqual(Double.acosh(.infinity), .infinity)
+
+    XCTAssertEqual(Double.atanh(0.5), 0.549306144334,
+                   accuracy: 0.000_000_000_001)
+    XCTAssertTrue(Double.atanh(0).isZero)
+    XCTAssertTrue(Double.atanh(0).sign == .plus)
+    XCTAssertTrue(Double.atanh(-0.0).isZero)
+    XCTAssertTrue(Double.atanh(-0.0).sign == .minus)
+
+    XCTAssertEqual(Double.erf(1), 0.842700792950, accuracy: 0.000_000_000_001)
+    XCTAssertTrue(Double.erf(0).isZero)
+    XCTAssertTrue(Double.erf(0).sign == .plus)
+    XCTAssertTrue(Double.erf(-0.0).isZero)
+    XCTAssertTrue(Double.erf(-0.0).sign == .minus)
+    XCTAssertEqual(Double.erf(.infinity), 1)
+    XCTAssertEqual(Double.erf(-.infinity), -1)
+
+    XCTAssertEqual(Double.erfc(1), 0.157299207050, accuracy: 0.000_001)
+    XCTAssertEqual(Double.erfc(0), 1)
+    XCTAssertEqual(Double.erfc(.infinity), 0)
+    XCTAssertEqual(Double.erfc(-.infinity), 2)
+
+    XCTAssertEqual(Double.tgamma(10), Double((2...9).reduce(1, *)),
+                   accuracy: 0.000_001)
+    XCTAssertEqual(Double.tgamma(1), 1)
+    XCTAssertEqual(Double.tgamma(.infinity), .infinity)
+
+    XCTAssertEqual(Double.lgamma(10), Double.log(Double((2...9).reduce(1, *))),
+                   accuracy: 0.000_000_000_001)
+    XCTAssertEqual(Double.lgamma(1), 0)
+    XCTAssertEqual(Double.lgamma(.infinity), .infinity)
+  }
+
+  static var allTests = [
+    ("testTranscendentalConstants", testTranscendentalConstants),
+    ("testExp2", testExp2),
+    ("testExp10", testExp10),
+    ("testExpm1", testExpm1),
+    ("testLog2", testLog2),
+    ("testLog10", testLog10),
+    ("testLog1p", testLog1p),
+    ("testCbrt", testCbrt),
+    ("testTan", testTan),
+    ("testTanh", testTanh),
+    ("testHypot", testHypot),
+    ("testAtan2", testAtan2),
+    ("testFloat", testFloat),
+    ("testDouble", testDouble),
+  ]
+}

Tests/NumericAnnexTests/RootExtractionTests.swift

@@ -0,0 +1,56 @@
+import XCTest
+@testable import NumericAnnex
+
+class RootExtractionTests : XCTestCase {
+  func testSqrt() {
+    XCTAssertEqual(Int.sqrt(0), 0)
+    XCTAssertEqual(Int.sqrt(25), 5)
+    XCTAssertEqual(Int.sqrt(27), 5)
+    XCTAssertEqual(Int.sqrt(256), 16)
+    XCTAssertEqual(Int.sqrt(512), 22)
+    XCTAssertEqual(Int.sqrt(1 << 32) * .sqrt(1 << 32), 1 << 32)
+    XCTAssertEqual(Int.sqrt(1 << 48) * .sqrt(1 << 48), 1 << 48)
+    XCTAssertEqual(Int.sqrt(1 << 50) * .sqrt(1 << 50), 1 << 50)
+    XCTAssertEqual(Int.sqrt(1 << 60) * .sqrt(1 << 60), 1 << 60)
+    XCTAssertEqual(Int.sqrt(1 << 62) * .sqrt(1 << 62), 1 << 62)
+
+    XCTAssertLessThanOrEqual(Int.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int8.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int16.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int32.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int64.sqrt(.max) * .sqrt(.max), .max)
+
+    XCTAssertLessThanOrEqual(UInt.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt8.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt16.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt32.sqrt(.max) * .sqrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt64.sqrt(.max) * .sqrt(.max), .max)
+  }
+
+  func testCbrt() {
+    XCTAssertEqual(UInt.cbrt(0), 0)
+    XCTAssertEqual(UInt.cbrt(25), 2)
+    XCTAssertEqual(UInt.cbrt(27), 3)
+    XCTAssertEqual(UInt.cbrt(256), 6)
+    XCTAssertEqual(UInt.cbrt(512), 8)
+
+    XCTAssertLessThanOrEqual(UInt.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt8.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt16.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt32.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(UInt64.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+
+    XCTAssertEqual(Int.cbrt(-27), -3)
+
+    XCTAssertLessThanOrEqual(Int.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int8.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int16.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int32.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+    XCTAssertLessThanOrEqual(Int64.cbrt(.max) * .cbrt(.max) * .cbrt(.max), .max)
+  }
+
+  static var allTests = [
+    ("testSqrt", testSqrt),
+    ("testCbrt", testCbrt),
+  ]
+}