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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
//! Ported from:
//!
//! https://github.com/llvm/llvm-project/commit/d674d96bc56c0f377879d01c9d8dfdaaa7859cdb/compiler-rt/lib/builtins/divsf3.c
const std = @import("std");
const builtin = @import("builtin");
const arch = builtin.cpu.arch;
const common = @import("common.zig");
const normalize = common.normalize;
pub const panic = common.panic;
comptime {
if (common.want_aeabi) {
@export(&__aeabi_fdiv, .{ .name = "__aeabi_fdiv", .linkage = common.linkage, .visibility = common.visibility });
} else {
@export(&__divsf3, .{ .name = "__divsf3", .linkage = common.linkage, .visibility = common.visibility });
}
}
pub fn __divsf3(a: f32, b: f32) callconv(.c) f32 {
return div(a, b);
}
fn __aeabi_fdiv(a: f32, b: f32) callconv(.{ .arm_aapcs = .{} }) f32 {
return div(a, b);
}
inline fn div(a: f32, b: f32) f32 {
const Z = std.meta.Int(.unsigned, 32);
const significandBits = std.math.floatMantissaBits(f32);
const exponentBits = std.math.floatExponentBits(f32);
const signBit = (@as(Z, 1) << (significandBits + exponentBits));
const maxExponent = ((1 << exponentBits) - 1);
const exponentBias = (maxExponent >> 1);
const implicitBit = (@as(Z, 1) << significandBits);
const quietBit = implicitBit >> 1;
const significandMask = implicitBit - 1;
const absMask = signBit - 1;
const exponentMask = absMask ^ significandMask;
const qnanRep = exponentMask | quietBit;
const infRep: Z = @bitCast(std.math.inf(f32));
const aExponent: u32 = @truncate((@as(Z, @bitCast(a)) >> significandBits) & maxExponent);
const bExponent: u32 = @truncate((@as(Z, @bitCast(b)) >> significandBits) & maxExponent);
const quotientSign: Z = (@as(Z, @bitCast(a)) ^ @as(Z, @bitCast(b))) & signBit;
var aSignificand: Z = @as(Z, @bitCast(a)) & significandMask;
var bSignificand: Z = @as(Z, @bitCast(b)) & significandMask;
var scale: i32 = 0;
// Detect if a or b is zero, denormal, infinity, or NaN.
if (aExponent -% 1 >= maxExponent - 1 or bExponent -% 1 >= maxExponent - 1) {
const aAbs: Z = @as(Z, @bitCast(a)) & absMask;
const bAbs: Z = @as(Z, @bitCast(b)) & absMask;
// NaN / anything = qNaN
if (aAbs > infRep) return @bitCast(@as(Z, @bitCast(a)) | quietBit);
// anything / NaN = qNaN
if (bAbs > infRep) return @bitCast(@as(Z, @bitCast(b)) | quietBit);
if (aAbs == infRep) {
// infinity / infinity = NaN
if (bAbs == infRep) {
return @bitCast(qnanRep);
}
// infinity / anything else = +/- infinity
else {
return @bitCast(aAbs | quotientSign);
}
}
// anything else / infinity = +/- 0
if (bAbs == infRep) return @bitCast(quotientSign);
if (aAbs == 0) {
// zero / zero = NaN
if (bAbs == 0) {
return @bitCast(qnanRep);
}
// zero / anything else = +/- zero
else {
return @bitCast(quotientSign);
}
}
// anything else / zero = +/- infinity
if (bAbs == 0) return @bitCast(infRep | quotientSign);
// one or both of a or b is denormal, the other (if applicable) is a
// normal number. Renormalize one or both of a and b, and set scale to
// include the necessary exponent adjustment.
if (aAbs < implicitBit) scale +%= normalize(f32, &aSignificand);
if (bAbs < implicitBit) scale -%= normalize(f32, &bSignificand);
}
// Or in the implicit significand bit. (If we fell through from the
// denormal path it was already set by normalize( ), but setting it twice
// won't hurt anything.)
aSignificand |= implicitBit;
bSignificand |= implicitBit;
var quotientExponent: i32 = @as(i32, @bitCast(aExponent -% bExponent)) +% scale;
// Align the significand of b as a Q31 fixed-point number in the range
// [1, 2.0) and get a Q32 approximate reciprocal using a small minimax
// polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
// is accurate to about 3.5 binary digits.
const q31b = bSignificand << 8;
var reciprocal = @as(u32, 0x7504f333) -% q31b;
// Now refine the reciprocal estimate using a Newton-Raphson iteration:
//
// x1 = x0 * (2 - x0 * b)
//
// This doubles the number of correct binary digits in the approximation
// with each iteration, so after three iterations, we have about 28 binary
// digits of accuracy.
var correction: u32 = undefined;
correction = @truncate(~(@as(u64, reciprocal) *% q31b >> 32) +% 1);
reciprocal = @truncate(@as(u64, reciprocal) *% correction >> 31);
correction = @truncate(~(@as(u64, reciprocal) *% q31b >> 32) +% 1);
reciprocal = @truncate(@as(u64, reciprocal) *% correction >> 31);
correction = @truncate(~(@as(u64, reciprocal) *% q31b >> 32) +% 1);
reciprocal = @truncate(@as(u64, reciprocal) *% correction >> 31);
// Exhaustive testing shows that the error in reciprocal after three steps
// is in the interval [-0x1.f58108p-31, 0x1.d0e48cp-29], in line with our
// expectations. We bump the reciprocal by a tiny value to force the error
// to be strictly positive (in the range [0x1.4fdfp-37,0x1.287246p-29], to
// be specific). This also causes 1/1 to give a sensible approximation
// instead of zero (due to overflow).
reciprocal -%= 2;
// The numerical reciprocal is accurate to within 2^-28, lies in the
// interval [0x1.000000eep-1, 0x1.fffffffcp-1], and is strictly smaller
// than the true reciprocal of b. Multiplying a by this reciprocal thus
// gives a numerical q = a/b in Q24 with the following properties:
//
// 1. q < a/b
// 2. q is in the interval [0x1.000000eep-1, 0x1.fffffffcp0)
// 3. the error in q is at most 2^-24 + 2^-27 -- the 2^24 term comes
// from the fact that we truncate the product, and the 2^27 term
// is the error in the reciprocal of b scaled by the maximum
// possible value of a. As a consequence of this error bound,
// either q or nextafter(q) is the correctly rounded
var quotient: Z = @truncate(@as(u64, reciprocal) *% (aSignificand << 1) >> 32);
// Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
// In either case, we are going to compute a residual of the form
//
// r = a - q*b
//
// We know from the construction of q that r satisfies:
//
// 0 <= r < ulp(q)*b
//
// if r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
// already have the correct result. The exact halfway case cannot occur.
// We also take this time to right shift quotient if it falls in the [1,2)
// range and adjust the exponent accordingly.
var residual: Z = undefined;
if (quotient < (implicitBit << 1)) {
residual = (aSignificand << 24) -% quotient *% bSignificand;
quotientExponent -%= 1;
} else {
quotient >>= 1;
residual = (aSignificand << 23) -% quotient *% bSignificand;
}
const writtenExponent = quotientExponent +% exponentBias;
if (writtenExponent >= maxExponent) {
// If we have overflowed the exponent, return infinity.
return @bitCast(infRep | quotientSign);
} else if (writtenExponent < 1) {
if (writtenExponent == 0) {
// Check whether the rounded result is normal.
const round = @intFromBool((residual << 1) > bSignificand);
// Clear the implicit bit.
var absResult = quotient & significandMask;
// Round.
absResult += round;
if ((absResult & ~significandMask) > 0) {
// The rounded result is normal; return it.
return @bitCast(absResult | quotientSign);
}
}
// Flush denormals to zero. In the future, it would be nice to add
// code to round them correctly.
return @bitCast(quotientSign);
} else {
const round = @intFromBool((residual << 1) > bSignificand);
// Clear the implicit bit
var absResult = quotient & significandMask;
// Insert the exponent
absResult |= @as(Z, @bitCast(writtenExponent)) << significandBits;
// Round
absResult +%= round;
// Insert the sign and return
return @bitCast(absResult | quotientSign);
}
}
test {
_ = @import("divsf3_test.zig");
}
|
