lib/simd/src/fastmath.zig
daab053ee43316e1809a84551d573ddd1e5bf3d2
1 const std = @import("std");
2 const math = @import("math.zig");
3
4 pub fn fastTan(comptime D: type, value: D.Vector) D.Vector {
5 return math.tan(D, value);
6 }
7
8 pub fn fastAtan(comptime D: type, value: D.Vector) D.Vector {
9 return math.atan(D, value);
10 }
11
12 pub fn fastAtanPositive(comptime D: type, value: D.Vector) D.Vector {
13 return math.atan(D, value);
14 }
15
16 pub fn fastAtan2(comptime D: type, y: D.Vector, x: D.Vector) D.Vector {
17 return math.atan2(D, y, x);
18 }
19
20 pub fn fastTanh(comptime D: type, value: D.Vector) D.Vector {
21 return math.tanh(D, value);
22 }
23
24 pub fn fastLog(comptime D: type, value: D.Vector) D.Vector {
25 return math.log(D, value);
26 }
27
28 pub fn fastExp(comptime D: type, value: D.Vector) D.Vector {
29 return math.exp(D, value);
30 }
31
32 pub fn fastExp2(comptime D: type, value: D.Vector) D.Vector {
33 return math.exp2(D, value);
34 }
35
36 pub fn fastExpMinusOrZero(comptime D: type, value: D.Vector) D.Vector {
37 return flushSubnormals(D, math.exp(D, value));
38 }
39
40 pub fn fastLog2(comptime D: type, value: D.Vector) D.Vector {
41 return math.log2(D, value);
42 }
43
44 pub fn fastLog10(comptime D: type, value: D.Vector) D.Vector {
45 return math.log10(D, value);
46 }
47
48 pub fn fastLog1p(comptime D: type, value: D.Vector) D.Vector {
49 return math.log1p(D, value);
50 }
51
52 pub fn fastPow(comptime D: type, base: D.Vector, exponent: D.Vector) D.Vector {
53 return math.pow(D, base, exponent);
54 }
55
56 pub fn fastExpNormal(comptime D: type, value: D.Vector) D.Vector {
57 return flushSubnormals(D, math.exp(D, value));
58 }
59
60 pub fn fastExp2Normal(comptime D: type, value: D.Vector) D.Vector {
61 return flushSubnormals(D, math.exp2(D, value));
62 }
63
64 pub fn fastLogPositiveNormal(comptime D: type, value: D.Vector) D.Vector {
65 return math.log(D, value);
66 }
67
68 pub fn fastLog2PositiveNormal(comptime D: type, value: D.Vector) D.Vector {
69 return math.log2(D, value);
70 }
71
72 pub fn fastLog10PositiveNormal(comptime D: type, value: D.Vector) D.Vector {
73 return math.log10(D, value);
74 }
75
76 pub fn fastLog1pPositiveNormal(comptime D: type, value: D.Vector) D.Vector {
77 return math.log1p(D, value);
78 }
79
80 pub fn fastPowNormal(comptime D: type, base: D.Vector, exponent: D.Vector) D.Vector {
81 return flushSubnormals(D, math.pow(D, base, exponent));
82 }
83
84 fn flushSubnormals(comptime D: type, value: D.Vector) D.Vector {
85 var result = value;
86 inline for (0..D.lane_count) |index| {
87 if (@abs(result[index]) < std.math.floatMin(D.Lane)) {
88 result[index] = std.math.copysign(@as(D.Lane, 0), result[index]);
89 }
90 }
91 return result;
92 }
93
94 test "Highway fast math stays within the documented relative bounds" {
95 const simd = @import("root.zig");
96 inline for (.{ f32, f64 }) |T| {
97 const D = simd.FixedTag(T, 4);
98 const values: D.Vector = .{ -0.75, -0.125, 0.125, 0.75 };
99 const positive: D.Vector = .{ 0.25, 0.5, 2, 8 };
100 const tangent = fastTan(D, values);
101 const inverse = fastAtan(D, values);
102 const hyperbolic = fastTanh(D, values);
103 const logarithm = fastLog(D, positive);
104 const exponent = fastExp(D, values);
105 inline for (0..D.lane_count) |index| {
106 try std.testing.expectApproxEqRel(@tan(values[index]), tangent[index], @as(T, 0.002));
107 try std.testing.expectApproxEqRel(std.math.atan(values[index]), inverse[index], @as(T, 0.000006));
108 try std.testing.expectApproxEqRel(std.math.tanh(values[index]), hyperbolic[index], @as(T, 0.000006));
109 try std.testing.expectApproxEqRel(@log(positive[index]), logarithm[index], @as(T, 0.000012));
110 try std.testing.expectApproxEqRel(@exp(values[index]), exponent[index], @as(T, 0.000007));
111 }
112 }
113 }
114
115 test "Highway fast logarithm and exponent variants retain identities" {
116 const simd = @import("root.zig");
117 const D = simd.FixedTag(f64, 4);
118 const values: D.Vector = .{ 0.125, 0.5, 2, 8 };
119 const log_two = fastLog2(D, values);
120 const exp_two = fastExp2(D, log_two);
121 inline for (0..D.lane_count) |index| {
122 try std.testing.expectApproxEqRel(values[index], exp_two[index], 0x1p-50);
123 }
124
125 const log_ten = fastLog10(D, @as(D.Vector, .{ 0.1, 1, 10, 100 }));
126 try std.testing.expectApproxEqAbs(@as(f64, -1), log_ten[0], 0x1p-50);
127 try std.testing.expectApproxEqAbs(@as(f64, 2), log_ten[3], 0x1p-50);
128 const log_one = fastLog1p(D, @as(D.Vector, .{ -0.5, -0.25, 0.25, 1 }));
129 try std.testing.expectApproxEqAbs(std.math.log1p(@as(f64, -0.5)), log_one[0], 0x1p-50);
130 }
131
132 test "Highway fast normal modes flush subnormal results" {
133 const simd = @import("root.zig");
134 const D32 = simd.FixedTag(f32, 4);
135 const result32 = fastExpNormal(D32, @as(D32.Vector, .{ -100, -90, -80, 0 }));
136 try std.testing.expectEqual(@as(f32, 0), result32[0]);
137 try std.testing.expectEqual(@as(f32, 0), result32[1]);
138 try std.testing.expect(result32[2] > 0);
139 try std.testing.expectEqual(@as(f32, 1), result32[3]);
140
141 const D64 = simd.FixedTag(f64, 4);
142 const result64 = fastExpMinusOrZero(D64, @as(D64.Vector, .{ -1000, -800, -700, 0 }));
143 try std.testing.expectEqual(@as(f64, 0), result64[0]);
144 try std.testing.expectEqual(@as(f64, 0), result64[1]);
145 try std.testing.expect(result64[2] > 0);
146 try std.testing.expectEqual(@as(f64, 1), result64[3]);
147 }
148
149 test "Highway fast binary and positive-only entry points remain available" {
150 const simd = @import("root.zig");
151 const D = simd.FixedTag(f32, 4);
152 const positive: D.Vector = .{ 0, 0.25, 1, 8 };
153 const signed_values: D.Vector = .{ -2, -0.5, 0.5, 2 };
154 const positive_atan = fastAtanPositive(D, positive);
155 const ordinary_atan = fastAtan(D, positive);
156 try std.testing.expect(@reduce(.And, positive_atan == ordinary_atan));
157
158 const quadrants = fastAtan2(D, signed_values, @as(D.Vector, .{ -1, 1, 1, -1 }));
159 inline for (0..D.lane_count) |index| {
160 try std.testing.expectApproxEqAbs(
161 std.math.atan2(signed_values[index], @as(D.Vector, .{ -1, 1, 1, -1 })[index]),
162 quadrants[index],
163 0.000001,
164 );
165 }
166
167 const powers = fastPow(D, @as(D.Vector, .{ 2, 4, 8, 16 }), @as(D.Vector, .{ 3, 0.5, -1, 0.25 }));
168 const normal_powers = fastPowNormal(D, @as(D.Vector, .{ 2, 4, 8, 16 }), @as(D.Vector, .{ 3, 0.5, -1, 0.25 }));
169 try std.testing.expect(@reduce(.And, powers == normal_powers));
170 _ = fastLogPositiveNormal(D, positive + @as(D.Vector, @splat(1)));
171 _ = fastLog2PositiveNormal(D, positive + @as(D.Vector, @splat(1)));
172 _ = fastLog10PositiveNormal(D, positive + @as(D.Vector, @splat(1)));
173 _ = fastLog1pPositiveNormal(D, positive);
174 _ = fastExp2Normal(D, signed_values);
175 }