Let a, b, and c be 128 bit numbers. A phantom overflow happens when the intermediary value a * b overflows 128 bits. However the final product after division numerator / c fits inside 128 bits.
For 64 bit and smaller numbers, one can cast to 128 bits or use the muldiv crate. What is the best way to do this for 128 bit numbers? Issue is that:
Reference- muldiv in solidity by Remco Bloemen
Unless there exists simpler/easier method (of which I am unaware), this is NOT as trivial as it looks. It took longer than expected to make, but I came up with the following code. I have not tested it well, I cannot guarantee it is accurate, and it almost certainly isn't optimal. However, it seemed to work with a few test cases. Use at your own risk ♂️:
// Inspired by: https://medium.com/wicketh/mathemagic-512-bit-division-in-solidity-afa55870a65
/// Returns the least significant 64 bits of a
fn lo128(a: u128) -> u128 {
a & ((1<<64)-1)
}
/// Returns the most significant 64 bits of a
fn hi128(a: u128) -> u128 {
a >> 64
}
/// Returns 2^128 - a (two's complement)
fn neg128(a: u128) -> u128 {
(!a).wrapping_add(1)
}
/// Returns 2^128 / a
fn div128(a: u128) -> u128 {
(neg128(a)/a).wrapping_add(1)
}
/// Returns 2^128 % a
fn mod128(a: u128) -> u128 {
neg128(a) % a
}
/// Returns a+b (wrapping)
fn add256(a0: u128, a1: u128, b0: u128, b1: u128) -> (u128, u128) {
let (r0, overflow) = a0.overflowing_add(b0);
let r1 = a1.wrapping_add(b1).wrapping_add(overflow as u128);
(r0, r1)
}
/// Returns a*b (in 256 bits)
fn mul128(a: u128, b: u128) -> (u128, u128) {
// Split a and b into hi and lo 64-bit parts
// a*b = (a1<<64 + a0)*(b1<<64 + b0)
// = (a1*b1)<<128 + (a1*b0 + b1*a0)<<64 + a0*b0
let (a0, a1) = (lo128(a), hi128(a));
let (b0, b1) = (lo128(b), hi128(b));
let (x, y) = (a1*b0, b1*a0);
let (r0, r1) = (a0*b0, a1*b1);
let (r0, r1) = add256(r0, r1, lo128(x)<<64, hi128(x));
let (r0, r1) = add256(r0, r1, lo128(y)<<64, hi128(y));
(r0, r1)
}
/// Returns a/b
fn div256_128(mut a0: u128, mut a1: u128, b: u128) -> (u128, u128) {
if b == 1 {
return (a0, a1);
}
// Calculate a/b
// = (a0<<128 + a1)/b
// let (q, r) = (div128(b), mod128(b));
// = (a1*(q*b + r)) + a0)/b
// = (a1*q*b + a1*r + a0)/b
// = (a1*r + a0)/b + a1*q
let (q, r) = (div128(b), mod128(b));
// x = current result
// a = next number
// t = temp
let (mut x0, mut x1) = (0, 0);
while a1 != 0 {
// x += a1*q
let (t0, t1) = mul128(a1, q);
let (new_x0, new_x1) = add256(x0, x1, t0, t1);
x0 = new_x0;
x1 = new_x1;
// a = a1*r + a0
let (t0, t1) = mul128(a1, r);
let (new_a0, new_a1) = add256(t0, t1, a0, 0);
a0 = new_a0;
a1 = new_a1;
}
add256(x0, x1, a0/b, 0)
}
/// Returns a*b/c (wrapping to 128 bits)
fn muldiv128(a: u128, b: u128, c: u128) -> u128 {
let (t0, t1) = mul128(a, b);
div256_128(t0, t1, c).0
}
fn main() {
let a = 54994939174662544931549714329375118023;
let b = 100391518367866121278432408876359705009;
let c = 183619791236921872854293928964252165930;
let result = muldiv128(a, b, c);
let expected = 30067703536211509756706998560804057558;
println!("{}", result == expected); // true
}
Alternatively, if you know that the numbers divide evenly, you can simply do (a/gcd(a, c)) * (b/gcd(b, c)) * gcd(a, b, c).
Multiplying two 128 bit numbers result in 256 bits of output. There may be an assembly instruction that lets you get at those bits.
Alternatively you can chop each 128 bit integer into an upper 64-bit half and a lower 64-bit half. Multiplying those halves results in a maximum of 128 bits of output and thus can be done in an 128-bit integer.
(a_up*2^64 + a_lo) * (b_up*2^64 + b_lo) = a_up*b_up*2^128 + a_up*b_lo*2^64 + a_lo*b_up*2^64 + a_lo*b_lo
Therefore the upper 128 bits are a_up*b_up*2^128 + (a_up*b_lo*2^64 + a_lo*b_up*2^64)_up and the lower 128 bits are (a_up*b_lo*2^64 + a_lo*b_up*2^64 + a_lo*b_lo)_lo.
Then you only need to divide this 256 bit result by c...
