1const builtin = @import("builtin");
2const std = @import("std");
3const math = std.math;
4const Log2Int = std.math.Log2Int;
5const assert = std.debug.assert;
6const expect = std.testing.expect;
7
8/// Returns x * 2^n.
9pub fn ldexp(x: anytype, n: i32) @TypeOf(x) {
10 const T = @TypeOf(x);
11 const TBits = @Int(.unsigned, @typeInfo(T).float.bits);
12
13 const exponent_bits = math.floatExponentBits(T);
14 const mantissa_bits = math.floatMantissaBits(T);
15 const fractional_bits = math.floatFractionalBits(T);
16
17 const max_biased_exponent = 2 * math.floatExponentMax(T);
18 const mantissa_mask = @as(TBits, (1 << mantissa_bits) - 1);
19
20 const repr = @as(TBits, @bitCast(x));
21 const sign_bit = repr & (1 << (exponent_bits + mantissa_bits));
22
23 if (math.isNan(x) or !math.isFinite(x))
24 return x;
25
26 var exponent: i32 = @as(i32, @intCast((repr << 1) >> (mantissa_bits + 1)));
27 if (exponent == 0)
28 exponent += (@as(i32, exponent_bits) + @intFromBool(T == f80)) - @clz(repr << 1);
29
30 if (n >= 0) {
31 if (n > max_biased_exponent - exponent) {
32 // Overflow. Return +/- inf
33 return @as(T, @bitCast(@as(TBits, @bitCast(math.inf(T))) | sign_bit));
34 } else if (exponent + n <= 0) {
35 // Result is subnormal
36 return @as(T, @bitCast((repr << @as(Log2Int(TBits), @intCast(n))) | sign_bit));
37 } else if (exponent <= 0) {
38 // Result is normal, but needs shifting
39 var result = @as(TBits, @intCast(n + exponent)) << mantissa_bits;
40 result |= (repr << @as(Log2Int(TBits), @intCast(1 - exponent))) & mantissa_mask;
41 return @as(T, @bitCast(result | sign_bit));
42 }
43
44 // Result needs no shifting
45 return @as(T, @bitCast(repr + (@as(TBits, @intCast(n)) << mantissa_bits)));
46 } else {
47 if (n <= -exponent) {
48 if (n < -(mantissa_bits + exponent))
49 return @as(T, @bitCast(sign_bit)); // Severe underflow. Return +/- 0
50
51 // Result underflowed, we need to shift and round
52 const shift = @as(Log2Int(TBits), @intCast(@min(-n, -(exponent + n) + 1)));
53 const exact_tie: bool = @ctz(repr) == shift - 1;
54 var result = repr & mantissa_mask;
55
56 if (T != f80) // Include integer bit
57 result |= @as(TBits, @intFromBool(exponent > 0)) << fractional_bits;
58 result = @as(TBits, @intCast((result >> (shift - 1))));
59
60 // Round result, including round-to-even for exact ties
61 result = ((result + 1) >> 1) & ~@as(TBits, @intFromBool(exact_tie));
62 return @as(T, @bitCast(result | sign_bit));
63 }
64
65 // Result is exact, and needs no shifting
66 return @as(T, @bitCast(repr - (@as(TBits, @intCast(-n)) << mantissa_bits)));
67 }
68}
69
70test ldexp {
71 if (builtin.zig_backend == .stage2_c and builtin.cpu.arch.isAarch64() and builtin.os.tag == .netbsd) return error.SkipZigTest; // https://codeberg.org/ziglang/zig/issues/36765
72 if (builtin.zig_backend == .stage2_llvm and builtin.cpu.arch.isAarch64() and builtin.os.tag == .netbsd) return error.SkipZigTest; // https://codeberg.org/ziglang/zig/issues/36765
73
74 // subnormals
75 try expect(ldexp(@as(f16, 0x1.1FFp14), -14 - 9 - 15) == math.floatTrueMin(f16));
76 try expect(ldexp(@as(f32, 0x1.3FFFFFp-1), -126 - 22) == math.floatTrueMin(f32));
77 try expect(ldexp(@as(f64, 0x1.7FFFFFFFFFFFFp-1), -1022 - 51) == math.floatTrueMin(f64));
78 try expect(ldexp(@as(f80, 0x1.7FFFFFFFFFFFFFFEp-1), -16382 - 62) == math.floatTrueMin(f80));
79 try expect(ldexp(@as(f128, 0x1.7FFFFFFFFFFFFFFFFFFFFFFFFFFFp-1), -16382 - 111) == math.floatTrueMin(f128));
80
81 try expect(ldexp(math.floatMax(f32), -128 - 149) > 0.0);
82 try expect(ldexp(math.floatMax(f32), -128 - 149 - 1) == 0.0);
83
84 @setEvalBranchQuota(10_000);
85
86 inline for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
87 const fractional_bits = math.floatFractionalBits(T);
88
89 const min_exponent = math.floatExponentMin(T);
90 const max_exponent = math.floatExponentMax(T);
91 const exponent_bias = max_exponent;
92
93 // basic usage
94 try expect(ldexp(@as(T, 1.5), 4) == 24.0);
95
96 // normals -> subnormals
97 try expect(math.isNormal(ldexp(@as(T, 1.0), min_exponent)));
98 try expect(!math.isNormal(ldexp(@as(T, 1.0), min_exponent - 1)));
99
100 // normals -> zero
101 try expect(ldexp(@as(T, 1.0), min_exponent - fractional_bits) > 0.0);
102 try expect(ldexp(@as(T, 1.0), min_exponent - fractional_bits - 1) == 0.0);
103
104 // subnormals -> zero
105 try expect(ldexp(math.floatTrueMin(T), 0) > 0.0);
106 try expect(ldexp(math.floatTrueMin(T), -1) == 0.0);
107
108 // Multiplications might flush the denormals to zero, esp. at
109 // runtime, so we manually construct the constants here instead.
110 const Z = @Int(.unsigned, @bitSizeOf(T));
111 const EightTimesTrueMin = @as(T, @bitCast(@as(Z, 8)));
112 const TwoTimesTrueMin = @as(T, @bitCast(@as(Z, 2)));
113
114 // subnormals -> subnormals
115 try expect(ldexp(math.floatTrueMin(T), 3) == EightTimesTrueMin);
116 try expect(ldexp(EightTimesTrueMin, -2) == TwoTimesTrueMin);
117 try expect(ldexp(EightTimesTrueMin, -3) == math.floatTrueMin(T));
118
119 // subnormals -> normals (+)
120 try expect(ldexp(math.floatTrueMin(T), fractional_bits) == math.floatMin(T));
121 try expect(ldexp(math.floatTrueMin(T), fractional_bits - 1) == math.floatMin(T) * 0.5);
122
123 // subnormals -> normals (-)
124 try expect(ldexp(-math.floatTrueMin(T), fractional_bits) == -math.floatMin(T));
125 try expect(ldexp(-math.floatTrueMin(T), fractional_bits - 1) == -math.floatMin(T) * 0.5);
126
127 // subnormals -> float limits (+inf)
128 try expect(math.isFinite(ldexp(math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits - 1)));
129 try expect(ldexp(math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits) == math.inf(T));
130
131 // subnormals -> float limits (-inf)
132 try expect(math.isFinite(ldexp(-math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits - 1)));
133 try expect(ldexp(-math.floatTrueMin(T), max_exponent + exponent_bias + fractional_bits) == -math.inf(T));
134
135 // infinity -> infinity
136 try expect(ldexp(math.inf(T), math.maxInt(i32)) == math.inf(T));
137 try expect(ldexp(math.inf(T), math.minInt(i32)) == math.inf(T));
138 try expect(ldexp(math.inf(T), max_exponent) == math.inf(T));
139 try expect(ldexp(math.inf(T), min_exponent) == math.inf(T));
140 try expect(ldexp(-math.inf(T), math.maxInt(i32)) == -math.inf(T));
141 try expect(ldexp(-math.inf(T), math.minInt(i32)) == -math.inf(T));
142
143 // extremely large n
144 try expect(ldexp(math.floatMax(T), math.maxInt(i32)) == math.inf(T));
145 try expect(ldexp(math.floatMax(T), -math.maxInt(i32)) == 0.0);
146 try expect(ldexp(math.floatMax(T), math.minInt(i32)) == 0.0);
147 try expect(ldexp(math.floatTrueMin(T), math.maxInt(i32)) == math.inf(T));
148 try expect(ldexp(math.floatTrueMin(T), -math.maxInt(i32)) == 0.0);
149 try expect(ldexp(math.floatTrueMin(T), math.minInt(i32)) == 0.0);
150 }
151}