| ... | @@ -1,67 +1,68 @@ | ... | @@ -1,67 +1,68 @@ |
| 1 | const builtin = @import("builtin"); | 1 | const builtin = @import("builtin"); |
| | 2 | const std = @import("std"); |
| | 3 | const math = std.math; |
| 2 | | 4 | |
| 3 | // mulo - multiplication overflow | 5 | // mulo - multiplication overflow |
| 4 | // - muloXi4_generic for unoptimized version | 6 | // * return a*b. |
| | 7 | // * return if a*b overflows => 1 else => 0 |
| | 8 | // - muloXi4_genericSmall as default |
| | 9 | // - muloXi4_genericFast for 2*bitsize <= usize |
| 5 | | 10 | |
| 6 | // return a*b. | 11 | inline fn muloXi4_genericSmall(comptime ST: type, a: ST, b: ST, overflow: *c_int) ST { |
| 7 | // return if a*b overflows => 1 else => 0 | 12 | @setRuntimeSafety(builtin.is_test); |
| 8 | // see https://stackoverflow.com/a/26320664 for possible implementations | 13 | overflow.* = 0; |
| | 14 | const min = math.minInt(ST); |
| | 15 | var res: ST = a *% b; |
| | 16 | // Hacker's Delight section Overflow subsection Multiplication |
| | 17 | // case a=-2^{31}, b=-1 problem, because |
| | 18 | // on some machines a*b = -2^{31} with overflow |
| | 19 | // Then -2^{31}/-1 overflows and any result is possible. |
| | 20 | // => check with a<0 and b=-2^{31} |
| | 21 | if ((a < 0 and b == min) or (a != 0 and @divTrunc(res, a) != b)) |
| | 22 | overflow.* = 1; |
| | 23 | return res; |
| | 24 | } |
| 9 | | 25 | |
| 10 | inline fn muloXi4_generic(comptime ST: type, a: ST, b: ST, overflow: *c_int) ST { | 26 | inline fn muloXi4_genericFast(comptime ST: type, a: ST, b: ST, overflow: *c_int) ST { |
| 11 | @setRuntimeSafety(builtin.is_test); | 27 | @setRuntimeSafety(builtin.is_test); |
| 12 | const BSIZE = @bitSizeOf(ST); | 28 | overflow.* = 0; |
| 13 | comptime var UT = switch (ST) { | 29 | const EST = switch (ST) { |
| 14 | i32 => u32, | 30 | i32 => i64, |
| 15 | i64 => u64, | 31 | i64 => i128, |
| 16 | i128 => u128, | 32 | i128 => i256, |
| 17 | else => unreachable, | 33 | else => unreachable, |
| 18 | }; | 34 | }; |
| 19 | const min = @bitCast(ST, @as(UT, 1 << (BSIZE - 1))); | 35 | const min = math.minInt(ST); |
| 20 | const max = ~min; | 36 | const max = math.maxInt(ST); |
| 21 | overflow.* = 0; | 37 | var res: EST = @as(EST, a) * @as(EST, b); |
| 22 | const result = a *% b; | 38 | //invariant: -2^{bitwidth(EST)} < res < 2^{bitwidth(EST)-1} |
| 23 | | 39 | if (res < min or max < res) |
| 24 | // edge cases | 40 | overflow.* = 1; |
| 25 | if (a == min) { | 41 | return @truncate(ST, res); |
| 26 | if (b != 0 and b != 1) overflow.* = 1; | | |
| 27 | return result; | | |
| 28 | } | | |
| 29 | if (b == min) { | | |
| 30 | if (a != 0 and a != 1) overflow.* = 1; | | |
| 31 | return result; | | |
| 32 | } | | |
| 33 | | | |
| 34 | // take sign of x sx | | |
| 35 | const sa = a >> (BSIZE - 1); | | |
| 36 | const sb = b >> (BSIZE - 1); | | |
| 37 | // take absolute value of a and b via | | |
| 38 | // abs(x) = (x^sx)) - sx | | |
| 39 | const abs_a = (a ^ sa) -% sa; | | |
| 40 | const abs_b = (b ^ sb) -% sb; | | |
| 41 | | | |
| 42 | // unitary magnitude, cannot have overflow | | |
| 43 | if (abs_a < 2 or abs_b < 2) return result; | | |
| 44 | | | |
| 45 | // compare the signs of operands | | |
| 46 | if ((a ^ b) >> (BSIZE - 1) != 0) { | | |
| 47 | if (abs_a > @divTrunc(max, abs_b)) overflow.* = 1; | | |
| 48 | } else { | | |
| 49 | if (abs_a > @divTrunc(min, -abs_b)) overflow.* = 1; | | |
| 50 | } | | |
| 51 | | | |
| 52 | return result; | | |
| 53 | } | 42 | } |
| 54 | | 43 | |
| 55 | pub fn __mulosi4(a: i32, b: i32, overflow: *c_int) callconv(.C) i32 { | 44 | pub fn __mulosi4(a: i32, b: i32, overflow: *c_int) callconv(.C) i32 { |
| 56 | return muloXi4_generic(i32, a, b, overflow); | 45 | if (2 * @bitSizeOf(i32) <= @bitSizeOf(usize)) { |
| | 46 | return muloXi4_genericFast(i32, a, b, overflow); |
| | 47 | } else { |
| | 48 | return muloXi4_genericSmall(i32, a, b, overflow); |
| | 49 | } |
| 57 | } | 50 | } |
| 58 | | 51 | |
| 59 | pub fn __mulodi4(a: i64, b: i64, overflow: *c_int) callconv(.C) i64 { | 52 | pub fn __mulodi4(a: i64, b: i64, overflow: *c_int) callconv(.C) i64 { |
| 60 | return muloXi4_generic(i64, a, b, overflow); | 53 | if (2 * @bitSizeOf(i64) <= @bitSizeOf(usize)) { |
| | 54 | return muloXi4_genericFast(i64, a, b, overflow); |
| | 55 | } else { |
| | 56 | return muloXi4_genericSmall(i64, a, b, overflow); |
| | 57 | } |
| 61 | } | 58 | } |
| 62 | | 59 | |
| 63 | pub fn __muloti4(a: i128, b: i128, overflow: *c_int) callconv(.C) i128 { | 60 | pub fn __muloti4(a: i128, b: i128, overflow: *c_int) callconv(.C) i128 { |
| 64 | return muloXi4_generic(i128, a, b, overflow); | 61 | if (2 * @bitSizeOf(i128) <= @bitSizeOf(usize)) { |
| | 62 | return muloXi4_genericFast(i128, a, b, overflow); |
| | 63 | } else { |
| | 64 | return muloXi4_genericSmall(i128, a, b, overflow); |
| | 65 | } |
| 65 | } | 66 | } |
| 66 | | 67 | |
| 67 | test { | 68 | test { |