| ... | ... | @@ -225,6 +225,177 @@ fn addXf3(comptime T: type, a: T, b: T) T { |
| 225 | 225 | return @bitCast(T, result); |
| 226 | 226 | } |
| 227 | 227 | |
| 228 | fn normalize_f80(exp: *i32, significand: *u80) void { |
| 229 | const shift = @clz(u64, @truncate(u64, significand.*)); |
| 230 | significand.* = (significand.* << shift); |
| 231 | exp.* += -@as(i8, shift); |
| 232 | } |
| 233 | |
| 234 | pub fn __addxf3(a: f80, b: f80) callconv(.C) f80 { |
| 235 | var a_rep align(16) = @ptrCast(*const std.math.F80Repr, &a).*; |
| 236 | var b_rep align(16) = @ptrCast(*const std.math.F80Repr, &b).*; |
| 237 | var a_exp: i32 = a_rep.exp & 0x7FFF; |
| 238 | var b_exp: i32 = b_rep.exp & 0x7FFF; |
| 239 | |
| 240 | const significand_bits = std.math.floatMantissaBits(f80); |
| 241 | const int_bit = 0x8000000000000000; |
| 242 | const significand_mask = 0x7FFFFFFFFFFFFFFF; |
| 243 | const qnan_bit = 0xC000000000000000; |
| 244 | const max_exp = 0x7FFF; |
| 245 | const sign_bit = 0x8000; |
| 246 | |
| 247 | // Detect if a or b is infinity, or NaN. |
| 248 | if (a_exp == max_exp) { |
| 249 | if (a_rep.fraction ^ int_bit == 0) { |
| 250 | if (b_exp == max_exp and (b_rep.fraction ^ int_bit == 0)) { |
| 251 | // +/-infinity + -/+infinity = qNaN |
| 252 | return std.math.qnan_f80; |
| 253 | } |
| 254 | // +/-infinity + anything = +/- infinity |
| 255 | return a; |
| 256 | } else { |
| 257 | std.debug.assert(a_rep.fraction & significand_mask != 0); |
| 258 | // NaN + anything = qNaN |
| 259 | a_rep.fraction |= qnan_bit; |
| 260 | return @ptrCast(*const f80, &a_rep).*; |
| 261 | } |
| 262 | } |
| 263 | if (b_exp == max_exp) { |
| 264 | if (b_rep.fraction ^ int_bit == 0) { |
| 265 | // anything + +/-infinity = +/-infinity |
| 266 | return b; |
| 267 | } else { |
| 268 | std.debug.assert(b_rep.fraction & significand_mask != 0); |
| 269 | // anything + NaN = qNaN |
| 270 | b_rep.fraction |= qnan_bit; |
| 271 | return @ptrCast(*const f80, &b_rep).*; |
| 272 | } |
| 273 | } |
| 274 | |
| 275 | const a_zero = (a_rep.fraction | @bitCast(u32, a_exp)) == 0; |
| 276 | const b_zero = (b_rep.fraction | @bitCast(u32, b_exp)) == 0; |
| 277 | if (a_zero) { |
| 278 | // zero + anything = anything |
| 279 | if (b_zero) { |
| 280 | // but we need to get the sign right for zero + zero |
| 281 | a_rep.exp &= b_rep.exp; |
| 282 | return @ptrCast(*const f80, &a_rep).*; |
| 283 | } else { |
| 284 | return b; |
| 285 | } |
| 286 | } else if (b_zero) { |
| 287 | // anything + zero = anything |
| 288 | return a; |
| 289 | } |
| 290 | |
| 291 | var a_int: u80 = a_rep.fraction | (@as(u80, a_rep.exp & max_exp) << significand_bits); |
| 292 | var b_int: u80 = b_rep.fraction | (@as(u80, b_rep.exp & max_exp) << significand_bits); |
| 293 | |
| 294 | // Swap a and b if necessary so that a has the larger absolute value. |
| 295 | if (b_int > a_int) { |
| 296 | const temp = a_rep; |
| 297 | a_rep = b_rep; |
| 298 | b_rep = temp; |
| 299 | } |
| 300 | |
| 301 | // Extract the exponent and significand from the (possibly swapped) a and b. |
| 302 | a_exp = a_rep.exp & max_exp; |
| 303 | b_exp = b_rep.exp & max_exp; |
| 304 | a_int = a_rep.fraction; |
| 305 | b_int = b_rep.fraction; |
| 306 | |
| 307 | // Normalize any denormals, and adjust the exponent accordingly. |
| 308 | normalize_f80(&a_exp, &a_int); |
| 309 | normalize_f80(&b_exp, &b_int); |
| 310 | |
| 311 | // The sign of the result is the sign of the larger operand, a. If they |
| 312 | // have opposite signs, we are performing a subtraction; otherwise addition. |
| 313 | const result_sign = a_rep.exp & sign_bit; |
| 314 | const subtraction = (a_rep.exp ^ b_rep.exp) & sign_bit != 0; |
| 315 | |
| 316 | // Shift the significands to give us round, guard and sticky, and or in the |
| 317 | // implicit significand bit. (If we fell through from the denormal path it |
| 318 | // was already set by normalize( ), but setting it twice won't hurt |
| 319 | // anything.) |
| 320 | a_int = a_int << 3; |
| 321 | b_int = b_int << 3; |
| 322 | |
| 323 | // Shift the significand of b by the difference in exponents, with a sticky |
| 324 | // bottom bit to get rounding correct. |
| 325 | const @"align" = @intCast(u80, a_exp - b_exp); |
| 326 | if (@"align" != 0) { |
| 327 | if (@"align" < 80) { |
| 328 | const sticky = if (b_int << @intCast(u7, 80 - @"align") != 0) @as(u80, 1) else 0; |
| 329 | b_int = (b_int >> @truncate(u7, @"align")) | sticky; |
| 330 | } else { |
| 331 | b_int = 1; // sticky; b is known to be non-zero. |
| 332 | } |
| 333 | } |
| 334 | if (subtraction) { |
| 335 | a_int -= b_int; |
| 336 | // If a == -b, return +zero. |
| 337 | if (a_int == 0) return 0.0; |
| 338 | |
| 339 | // If partial cancellation occurred, we need to left-shift the result |
| 340 | // and adjust the exponent: |
| 341 | if (a_int < int_bit << 3) { |
| 342 | const shift = @intCast(i32, @clz(u80, a_int)) - @intCast(i32, @clz(u80, int_bit << 3)); |
| 343 | a_int <<= @intCast(u7, shift); |
| 344 | a_exp -= shift; |
| 345 | } |
| 346 | } else { // addition |
| 347 | a_int += b_int; |
| 348 | |
| 349 | // If the addition carried up, we need to right-shift the result and |
| 350 | // adjust the exponent: |
| 351 | if (a_int & (int_bit << 4) != 0) { |
| 352 | const sticky = a_int & 1; |
| 353 | a_int = a_int >> 1 | sticky; |
| 354 | a_exp += 1; |
| 355 | } |
| 356 | } |
| 357 | |
| 358 | // If we have overflowed the type, return +/- infinity: |
| 359 | if (a_exp >= max_exp) { |
| 360 | a_rep.exp = max_exp | result_sign; |
| 361 | a_rep.fraction = int_bit; // integer bit is set for +/-inf |
| 362 | return @ptrCast(*const f80, &a_rep).*; |
| 363 | } |
| 364 | |
| 365 | if (a_exp <= 0) { |
| 366 | // Result is denormal before rounding; the exponent is zero and we |
| 367 | // need to shift the significand. |
| 368 | const shift = @intCast(u80, 1 - a_exp); |
| 369 | const sticky = if (a_int << @intCast(u7, 80 - shift) != 0) @as(u1, 1) else 0; |
| 370 | a_int = a_int >> @intCast(u7, shift | sticky); |
| 371 | a_exp = 0; |
| 372 | } |
| 373 | |
| 374 | // Low three bits are round, guard, and sticky. |
| 375 | const round_guard_sticky = @truncate(u3, a_int); |
| 376 | |
| 377 | // Shift the significand into place. |
| 378 | a_int = @truncate(u64, a_int >> 3); |
| 379 | |
| 380 | // // Insert the exponent and sign. |
| 381 | a_int |= (@intCast(u80, a_exp) | result_sign) << significand_bits; |
| 382 | |
| 383 | // Final rounding. The result may overflow to infinity, but that is the |
| 384 | // correct result in that case. |
| 385 | if (round_guard_sticky > 0x4) a_int += 1; |
| 386 | if (round_guard_sticky == 0x4) a_int += a_int & 1; |
| 387 | |
| 388 | a_rep.fraction = @truncate(u64, a_int); |
| 389 | a_rep.exp = @truncate(u16, a_int >> significand_bits); |
| 390 | return @ptrCast(*const f80, &a_rep).*; |
| 391 | } |
| 392 | |
| 393 | pub fn __subxf3(a: f80, b: f80) callconv(.C) f80 { |
| 394 | var b_rep align(16) = @ptrCast(*const std.math.F80Repr, &b).*; |
| 395 | b_rep.exp ^= 0x8000; |
| 396 | return __addxf3(a, @ptrCast(*const f80, &b_rep).*); |
| 397 | } |
| 398 | |
| 228 | 399 | test { |
| 229 | 400 | _ = @import("addXf3_test.zig"); |
| 230 | 401 | } |
Veikka Tuominen <git@vexu.eu> 2022-01-21 21:49:02+02:00
lib