1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
//===--- Real.swift -------------------------------------------*- swift -*-===//
//
// This source file is part of the Swift Numerics open source project
//
// Copyright (c) 2019 Apple Inc. and the Swift Numerics project authors
// Licensed under Apache License v2.0 with Runtime Library Exception
//
// See https://swift.org/LICENSE.txt for license information
//
//===----------------------------------------------------------------------===//
/// A type that models the real numbers.
///
/// Types conforming to this protocol provide the arithmetic and utility
/// operations defined by the `FloatingPoint` protocol, and provide all of the
/// math functions defined by the `ElementaryFunctions` and `RealFunctions`
/// protocols. This protocol does not add any additional conformances itself,
/// but is very useful as a protocol against which to write generic code. For
/// example, we can naturally write a generic implementation of a sigmoid
/// function:
/// ```
/// func sigmoid<T: Real>(_ x: T) -> T {
/// return 1/(1 + .exp(-x))
/// }
/// ```
/// See Also:
/// -
/// - `ElementaryFunctions`
/// - `RealFunctions`
/// - `AlgebraicField`
public protocol Real: FloatingPoint, RealFunctions, AlgebraicField {
}
// While `Real` does not provide any additional customization points,
// it does allow us to default the implementation of a few operations,
// and also provides `signGamma`.
extension Real {
// Most math libraries do not provide exp10, so we need a default
// implementation.
@_transparent
public static func exp10(_ x: Self) -> Self {
return pow(10, x)
}
/// cos(x) - 1, computed in such a way as to maintain accuracy for small x.
///
/// See also:
/// -
/// - `ElementaryFunctions.expMinusOne()`
@_transparent
public static func cosMinusOne(_ x: Self) -> Self {
let sinxOver2 = sin(x/2)
return -2*sinxOver2*sinxOver2
}
#if !os(Windows)
public static func signGamma(_ x: Self) -> FloatingPointSign {
// Gamma is strictly positive for x >= 0.
if x >= 0 { return .plus }
// For negative x, we arbitrarily choose to assign a sign of .plus to the
// poles.
let trunc = x.rounded(.towardZero)
if x == trunc { return .plus }
// Otherwise, signGamma is .minus if the integral part of x is even.
return trunc.isEven ? .minus : .plus
}
// Determines if this value is even, assuming that it is an integer.
@inline(__always)
private var isEven: Bool {
if Self.radix == 2 {
// For binary types, we can just check if x/2 is an integer. This works
// because x/2 is always computed exactly.
let half = self/2
return half == half.rounded(.towardZero)
} else {
// For decimal types, it's not quite that simple, because x/2 is not
// necessarily computed exactly. As an example, suppose that we had a
// decimal type with a one digit significand, and self = 7. Then self/2
// would round to 4, and we would (wrongly) conclude that it was an
// integer, and hence that self was even.
//
// Instead, for decimal types, we check if 2*trunc(self/2) == self,
// using an FMA; this is always correct; this approach works for any
// radix, but the previous method is more efficient for radix == 2.
let half = self/2
return self.addingProduct(-2, half.rounded(.towardZero)) == 0
}
}
#endif
@_transparent
public static func _mulAdd(_ a: Self, _ b: Self, _ c: Self) -> Self {
a*b + c
}
@_transparent
public static func sqrt(_ x: Self) -> Self {
return x.squareRoot()
}
/// The (approximate) reciprocal (multiplicative inverse) of this number,
/// if it is representable.
///
/// If `x` if finite and nonzero, and `1/x` overflows or underflows,
/// then `x.reciprocal` is `nil`. Otherwise, `a.reciprocal` is `1/x`.
///
/// If `x.reciprocal` is non-nil, you may be able to replace division by `x`
/// with multiplication by this value. It is not advantageous to do this
/// for an isolated division unless it is a compile-time constant visible
/// to the compiler, but if you are dividing many values by a single
/// denominator, this will often be a significant performance win.
///
/// A typical use case looks something like this:
/// ```
/// func divide<T: Real>(data: [T], by divisor: T) -> [T] {
/// // If divisor is well-scaled, multiply by reciprocal.
/// if let recip = divisor.reciprocal {
/// return data.map { $0 * recip }
/// }
/// // Fallback on using division.
/// return data.map { $0 / divisor }
/// }
/// ```
///
/// Error Bounds:
///
/// Multiplying by the reciprocal instead of dividing will slightly
/// perturb results. For example `5.0 / 3` is 1.6666666666666667, but
/// `5.0 * 3.reciprocal!` is 1.6666666666666665.
///
/// The error of a normal division is bounded by half an ulp of the
/// result; we can derive a quick error bound for multiplication by
/// the real reciprocal (when it exists) as follows (I will use circle
/// operators to denote real-number arithmetic, and normal operators
/// for floating-point arithmetic):
///
/// ```
/// a * b.reciprocal! = a * (1/b)
/// = a * (1 ⊘ b)(1 + δ₁)
/// = (a ⊘ b)(1 + δ₁)(1 + δ₂)
/// = (a ⊘ b)(1 + δ₁ + δ₂ + δ₁δ₂)
/// ```
///
/// where `0 < δᵢ <= ulpOfOne/2`. This gives a roughly 1-ulp error,
/// about twice the error bound we get using division. For most
/// purposes this is an acceptable error, but if you need to match
/// results obtained using division, you should not use this.
@inlinable
public var reciprocal: Self? {
let recip = 1/self
if recip.isNormal || isZero || !isFinite {
return recip
}
return nil
}
}
|
