| ... | @@ -0,0 +1,131 @@ |
| 1 | const std = @import("std"); |
| 2 | const builtin = @import("builtin"); |
| 3 | const is_test = builtin.is_test; |
| 4 | const native_arch = builtin.cpu.arch; |
| 5 | |
| 6 | // AArch64 is the only ABI (at the moment) to support f16 arguments without the |
| 7 | // need for extending them to wider fp types. |
| 8 | pub const F16T = if (native_arch.isAARCH64()) f16 else u16; |
| 9 | |
| 10 | pub fn __extendhfxf2(a: F16T) callconv(.C) f80 { |
| 11 | return extendF80(f16, @bitCast(u16, a)); |
| 12 | } |
| 13 | |
| 14 | pub fn __extendffxf2(a: f32) callconv(.C) f80 { |
| 15 | return extendF80(f32, @bitCast(u32, a)); |
| 16 | } |
| 17 | |
| 18 | pub fn __extenddfxf2(a: f64) callconv(.C) f80 { |
| 19 | return extendF80(f64, @bitCast(u64, a)); |
| 20 | } |
| 21 | |
| 22 | inline fn extendF80(comptime src_t: type, a: std.meta.Int(.unsigned, @typeInfo(src_t).Float.bits)) f80 { |
| 23 | @setRuntimeSafety(builtin.is_test); |
| 24 | |
| 25 | const src_rep_t = std.meta.Int(.unsigned, @typeInfo(src_t).Float.bits); |
| 26 | const src_sig_bits = std.math.floatMantissaBits(src_t); |
| 27 | const dst_int_bit = 0x8000000000000000; |
| 28 | const dst_sig_bits = std.math.floatMantissaBits(f80) - 1; // -1 for the integer bit |
| 29 | |
| 30 | const dst_exp_bias = 16383; |
| 31 | |
| 32 | const src_bits = @bitSizeOf(src_t); |
| 33 | const src_exp_bits = src_bits - src_sig_bits - 1; |
| 34 | const src_inf_exp = (1 << src_exp_bits) - 1; |
| 35 | const src_exp_bias = src_inf_exp >> 1; |
| 36 | |
| 37 | const src_min_normal = 1 << src_sig_bits; |
| 38 | const src_inf = src_inf_exp << src_sig_bits; |
| 39 | const src_sign_mask = 1 << (src_sig_bits + src_exp_bits); |
| 40 | const src_abs_mask = src_sign_mask - 1; |
| 41 | const src_qnan = 1 << (src_sig_bits - 1); |
| 42 | const src_nan_code = src_qnan - 1; |
| 43 | |
| 44 | var dst: std.math.F80Repr align(16) = undefined; |
| 45 | |
| 46 | // Break a into a sign and representation of the absolute value |
| 47 | const a_abs = a & src_abs_mask; |
| 48 | const sign: u16 = if (a & src_sign_mask != 0) 0x8000 else 0; |
| 49 | |
| 50 | if (a_abs -% src_min_normal < src_inf - src_min_normal) { |
| 51 | // a is a normal number. |
| 52 | // Extend to the destination type by shifting the significand and |
| 53 | // exponent into the proper position and rebiasing the exponent. |
| 54 | dst.exp = @intCast(u16, a_abs >> src_sig_bits); |
| 55 | dst.exp += dst_exp_bias - src_exp_bias; |
| 56 | dst.fraction = @as(u64, a_abs) << (dst_sig_bits - src_sig_bits); |
| 57 | dst.fraction |= dst_int_bit; // bit 64 is always set for normal numbers |
| 58 | } else if (a_abs >= src_inf) { |
| 59 | // a is NaN or infinity. |
| 60 | // Conjure the result by beginning with infinity, then setting the qNaN |
| 61 | // bit (if needed) and right-aligning the rest of the trailing NaN |
| 62 | // payload field. |
| 63 | dst.exp = 0x7fff; |
| 64 | dst.fraction = dst_int_bit; |
| 65 | dst.fraction |= @as(u64, a_abs & src_qnan) << (dst_sig_bits - src_sig_bits); |
| 66 | dst.fraction |= @as(u64, a_abs & src_nan_code) << (dst_sig_bits - src_sig_bits); |
| 67 | } else if (a_abs != 0) { |
| 68 | // a is denormal. |
| 69 | // renormalize the significand and clear the leading bit, then insert |
| 70 | // the correct adjusted exponent in the destination type. |
| 71 | const scale: u16 = @clz(src_rep_t, a_abs) - |
| 72 | @clz(src_rep_t, @as(src_rep_t, src_min_normal)); |
| 73 | |
| 74 | dst.fraction = @as(u64, a_abs) << @intCast(u6, dst_sig_bits - src_sig_bits + scale); |
| 75 | dst.fraction |= dst_int_bit; // bit 64 is always set for normal numbers |
| 76 | dst.exp = @truncate(u16, a_abs >> @intCast(u4, src_sig_bits - scale)); |
| 77 | dst.exp ^= 1; |
| 78 | dst.exp |= dst_exp_bias - src_exp_bias - scale + 1; |
| 79 | } else { |
| 80 | // a is zero. |
| 81 | dst.exp = 0; |
| 82 | dst.fraction = 0; |
| 83 | } |
| 84 | |
| 85 | dst.exp |= sign; |
| 86 | return @ptrCast(*const f80, &dst).*; |
| 87 | } |
| 88 | |
| 89 | pub fn __extendxftf2(a: f80) callconv(.C) f128 { |
| 90 | @setRuntimeSafety(builtin.is_test); |
| 91 | |
| 92 | const src_int_bit: u64 = 0x8000000000000000; |
| 93 | const src_sig_mask = ~src_int_bit; |
| 94 | const src_sig_bits = std.math.floatMantissaBits(f80) - 1; // -1 for the integer bit |
| 95 | const dst_sig_bits = std.math.floatMantissaBits(f128); |
| 96 | |
| 97 | const dst_bits = @bitSizeOf(f128); |
| 98 | |
| 99 | const dst_min_normal = @as(u128, 1) << dst_sig_bits; |
| 100 | |
| 101 | // Break a into a sign and representation of the absolute value |
| 102 | var a_rep = @ptrCast(*const std.math.F80Repr, &a).*; |
| 103 | const sign = a_rep.exp & 0x8000; |
| 104 | a_rep.exp &= 0x7FFF; |
| 105 | var abs_result: u128 = undefined; |
| 106 | |
| 107 | if (a_rep.exp == 0 and a_rep.fraction == 0) { |
| 108 | // zero |
| 109 | abs_result = 0; |
| 110 | } else if (a_rep.exp == 0x7FFF) { |
| 111 | // a is nan or infinite |
| 112 | abs_result = @as(u128, a_rep.fraction) << (dst_sig_bits - src_sig_bits); |
| 113 | abs_result |= @as(u128, a_rep.exp) << dst_sig_bits; |
| 114 | } else if (a_rep.fraction & src_int_bit != 0) { |
| 115 | // a is a normal value |
| 116 | abs_result = @as(u128, a_rep.fraction & src_sig_mask) << (dst_sig_bits - src_sig_bits); |
| 117 | abs_result |= @as(u128, a_rep.exp) << dst_sig_bits; |
| 118 | } else { |
| 119 | // a is denormal |
| 120 | // renormalize the significand and clear the leading bit and integer part, |
| 121 | // then insert the correct adjusted exponent in the destination type. |
| 122 | const scale: u32 = @clz(u64, a_rep.fraction); |
| 123 | abs_result = @as(u128, a_rep.fraction) << @intCast(u7, dst_sig_bits - src_sig_bits + scale + 1); |
| 124 | abs_result ^= dst_min_normal; |
| 125 | abs_result |= @as(u128, scale + 1) << dst_sig_bits; |
| 126 | } |
| 127 | |
| 128 | // Apply the signbit to (dst_t)abs(a). |
| 129 | const result: u128 align(@alignOf(f128)) = abs_result | @as(u128, sign) << (dst_bits - 16); |
| 130 | return @bitCast(f128, result); |
| 131 | } |
Veikka Tuominen <git@vexu.eu> 2022-01-20 13:50:30+02:00
lib