| ... | ... | @@ -0,0 +1,159 @@ |
| 1 | const std = @import("std"); |
| 2 | const builtin = @import("builtin"); |
| 3 | const native_arch = builtin.cpu.arch; |
| 4 | |
| 5 | // AArch64 is the only ABI (at the moment) to support f16 arguments without the |
| 6 | // need for extending them to wider fp types. |
| 7 | pub const F16T = if (native_arch.isAARCH64()) f16 else u16; |
| 8 | |
| 9 | pub fn __truncxfhf2(a: f80) callconv(.C) F16T { |
| 10 | return @bitCast(F16T, trunc(f16, a)); |
| 11 | } |
| 12 | |
| 13 | pub fn __truncxfff2(a: f80) callconv(.C) f32 { |
| 14 | return trunc(f32, a); |
| 15 | } |
| 16 | |
| 17 | pub fn __truncxfdf2(a: f80) callconv(.C) f64 { |
| 18 | return trunc(f64, a); |
| 19 | } |
| 20 | |
| 21 | inline fn trunc(comptime dst_t: type, a: f80) dst_t { |
| 22 | @setRuntimeSafety(builtin.is_test); |
| 23 | |
| 24 | const dst_rep_t = std.meta.Int(.unsigned, @typeInfo(dst_t).Float.bits); |
| 25 | const src_sig_bits = std.math.floatMantissaBits(f80) - 1; // -1 for the integer bit |
| 26 | const dst_sig_bits = std.math.floatMantissaBits(dst_t); |
| 27 | |
| 28 | const src_exp_bias = 16383; |
| 29 | |
| 30 | const round_mask = (1 << (src_sig_bits - dst_sig_bits)) - 1; |
| 31 | const halfway = 1 << (src_sig_bits - dst_sig_bits - 1); |
| 32 | |
| 33 | const dst_bits = @typeInfo(dst_t).Float.bits; |
| 34 | const dst_exp_bits = dst_bits - dst_sig_bits - 1; |
| 35 | const dst_inf_exp = (1 << dst_exp_bits) - 1; |
| 36 | const dst_exp_bias = dst_inf_exp >> 1; |
| 37 | |
| 38 | const underflow = src_exp_bias + 1 - dst_exp_bias; |
| 39 | const overflow = src_exp_bias + dst_inf_exp - dst_exp_bias; |
| 40 | |
| 41 | const dst_qnan = 1 << (dst_sig_bits - 1); |
| 42 | const dst_nan_mask = dst_qnan - 1; |
| 43 | |
| 44 | // Break a into a sign and representation of the absolute value |
| 45 | var a_rep = @ptrCast(*const std.math.F80Repr, &a).*; |
| 46 | const sign = a_rep.exp & 0x8000; |
| 47 | a_rep.exp &= 0x7FFF; |
| 48 | a_rep.fraction &= 0x7FFFFFFFFFFFFFFF; |
| 49 | var abs_result: dst_rep_t = undefined; |
| 50 | |
| 51 | if (a_rep.exp -% underflow < a_rep.exp -% overflow) { |
| 52 | // The exponent of a is within the range of normal numbers in the |
| 53 | // destination format. We can convert by simply right-shifting with |
| 54 | // rounding and adjusting the exponent. |
| 55 | abs_result = @as(dst_rep_t, a_rep.exp) << dst_sig_bits; |
| 56 | abs_result |= @truncate(dst_rep_t, a_rep.fraction >> (src_sig_bits - dst_sig_bits)); |
| 57 | abs_result -%= @as(dst_rep_t, src_exp_bias - dst_exp_bias) << dst_sig_bits; |
| 58 | |
| 59 | const round_bits = a_rep.fraction & round_mask; |
| 60 | if (round_bits > halfway) { |
| 61 | // Round to nearest |
| 62 | abs_result += 1; |
| 63 | } else if (round_bits == halfway) { |
| 64 | // Ties to even |
| 65 | abs_result += abs_result & 1; |
| 66 | } |
| 67 | } else if (a_rep.exp == 0x7FFF and a_rep.fraction != 0) { |
| 68 | // a is NaN. |
| 69 | // Conjure the result by beginning with infinity, setting the qNaN |
| 70 | // bit and inserting the (truncated) trailing NaN field. |
| 71 | abs_result = @intCast(dst_rep_t, dst_inf_exp) << dst_sig_bits; |
| 72 | abs_result |= dst_qnan; |
| 73 | abs_result |= @intCast(dst_rep_t, (a_rep.fraction >> (src_sig_bits - dst_sig_bits)) & dst_nan_mask); |
| 74 | } else if (a_rep.exp >= overflow) { |
| 75 | // a overflows to infinity. |
| 76 | abs_result = @intCast(dst_rep_t, dst_inf_exp) << dst_sig_bits; |
| 77 | } else { |
| 78 | // a underflows on conversion to the destination type or is an exact |
| 79 | // zero. The result may be a denormal or zero. Extract the exponent |
| 80 | // to get the shift amount for the denormalization. |
| 81 | const shift = src_exp_bias - dst_exp_bias - a_rep.exp; |
| 82 | |
| 83 | // Right shift by the denormalization amount with sticky. |
| 84 | if (shift > src_sig_bits) { |
| 85 | abs_result = 0; |
| 86 | } else { |
| 87 | const sticky = @boolToInt(a_rep.fraction << @intCast(u6, shift) != 0); |
| 88 | const denormalized_significand = a_rep.fraction >> @intCast(u6, shift) | sticky; |
| 89 | abs_result = @intCast(dst_rep_t, denormalized_significand >> (src_sig_bits - dst_sig_bits)); |
| 90 | const round_bits = denormalized_significand & round_mask; |
| 91 | if (round_bits > halfway) { |
| 92 | // Round to nearest |
| 93 | abs_result += 1; |
| 94 | } else if (round_bits == halfway) { |
| 95 | // Ties to even |
| 96 | abs_result += abs_result & 1; |
| 97 | } |
| 98 | } |
| 99 | } |
| 100 | |
| 101 | const result align(@alignOf(dst_t)) = abs_result | @as(dst_rep_t, sign) << dst_bits - 16; |
| 102 | return @bitCast(dst_t, result); |
| 103 | } |
| 104 | |
| 105 | pub fn __trunctfxf2(a: f128) callconv(.C) f80 { |
| 106 | const src_sig_bits = std.math.floatMantissaBits(f128); |
| 107 | const dst_sig_bits = std.math.floatMantissaBits(f80) - 1; // -1 for the integer bit |
| 108 | |
| 109 | // Various constants whose values follow from the type parameters. |
| 110 | // Any reasonable optimizer will fold and propagate all of these. |
| 111 | const src_bits = @typeInfo(f128).Float.bits; |
| 112 | const src_exp_bits = src_bits - src_sig_bits - 1; |
| 113 | const src_inf_exp = 0x7FFF; |
| 114 | |
| 115 | const src_inf = src_inf_exp << src_sig_bits; |
| 116 | const src_sign_mask = 1 << (src_sig_bits + src_exp_bits); |
| 117 | const src_abs_mask = src_sign_mask - 1; |
| 118 | const round_mask = (1 << (src_sig_bits - dst_sig_bits)) - 1; |
| 119 | const halfway = 1 << (src_sig_bits - dst_sig_bits - 1); |
| 120 | const src_qnan = 1 << (src_sig_bits - 1); |
| 121 | const src_nan_mask = src_qnan - 1; |
| 122 | |
| 123 | // Break a into a sign and representation of the absolute value |
| 124 | const a_rep = @bitCast(u128, a); |
| 125 | const a_abs = a_rep & src_abs_mask; |
| 126 | const sign: u16 = if (a_rep & src_sign_mask != 0) 0x8000 else 0; |
| 127 | |
| 128 | var res: std.math.F80Repr align(16) = undefined; |
| 129 | |
| 130 | if (a_abs > src_inf) { |
| 131 | // a is NaN. |
| 132 | // Conjure the result by beginning with infinity, setting the qNaN |
| 133 | // bit and inserting the (truncated) trailing NaN field. |
| 134 | res.exp = 0x7fff; |
| 135 | res.fraction = 0x8000000000000000; |
| 136 | res.fraction |= @truncate(u64, (a_abs & src_qnan) << (src_sig_bits - dst_sig_bits)); |
| 137 | res.fraction |= @truncate(u64, (a_abs & src_nan_mask) << (src_sig_bits - dst_sig_bits)); |
| 138 | } else { |
| 139 | // The exponent of a is within the range of normal numbers in the |
| 140 | // destination format. We can convert by simply right-shifting with |
| 141 | // rounding and adjusting the exponent. |
| 142 | res.fraction = @truncate(u64, a_abs >> (src_sig_bits - dst_sig_bits)); |
| 143 | res.exp = @truncate(u16, a_abs >> src_sig_bits); |
| 144 | |
| 145 | const round_bits = a_abs & round_mask; |
| 146 | if (round_bits > halfway) { |
| 147 | // Round to nearest |
| 148 | const exp = @addWithOverflow(u64, res.fraction, 1, &res.fraction); |
| 149 | res.exp += @boolToInt(exp); |
| 150 | } else if (round_bits == halfway) { |
| 151 | // Ties to even |
| 152 | const exp = @addWithOverflow(u64, res.fraction, res.fraction & 1, &res.fraction); |
| 153 | res.exp += @boolToInt(exp); |
| 154 | } |
| 155 | } |
| 156 | |
| 157 | res.exp |= sign; |
| 158 | return @ptrCast(*const f80, &res).*; |
| 159 | } |
Veikka Tuominen <git@vexu.eu> 2022-01-21 12:41:09+02:00
lib