128 bits (a * b) / c en Rust avec protection contre le débordement fantôme

Question

Soit a, b et c des nombres de 128 bits. Un débordement fantôme se produit lorsque la valeur intermédiaire a * b déborde de 128 bits. Cependant, le produit final après la division numerator / c tient dans 128 bits.

Pour les nombres de 64 bits et moins, on peut passer à 128 bits ou utiliser la fonction caisse de multivision . Quelle est la meilleure façon de procéder pour des nombres de 128 bits ? Le problème est le suivant :

• Rust n'a pas de types natifs de plus de 128 bits.
• muldiv crate n'a pas le support 128 bit

Référence- muldiv en solidité par Remco Bloemen

Answer 1

A moins qu'il n'existe une méthode plus simple/plus facile (dont je n'ai pas connaissance), ce n'est PAS aussi trivial qu'il n'y paraît. Cela a pris plus de temps que prévu pour le faire, mais j'ai trouvé le code suivant. Je ne l'ai pas bien testé, je ne peux pas garantir qu'il soit exact, et il n'est presque certainement pas optimal. Cependant, il a semblé fonctionner avec quelques cas de test. Utilisez-le à vos risques et périls :

Terrain de jeux

// 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
}

Alternativement, si vous savez que les nombres se divisent de façon égale, vous pouvez simplement faire (a/gcd(a, c)) * (b/gcd(b, c)) * gcd(a, b, c) .

Answer 2

La multiplication de deux nombres de 128 bits donne un résultat de 256 bits. Il existe peut-être une instruction d'assemblage qui vous permet d'accéder à ces bits.

Vous pouvez également découper chaque entier de 128 bits en une moitié supérieure de 64 bits et une moitié inférieure de 64 bits. La multiplication de ces moitiés donne un maximum de 128 bits de sortie et peut donc être effectuée dans un entier de 128 bits.

(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

Par conséquent, les 128 bits supérieurs sont a_up*b_up*2^128 + (a_up*b_lo*2^64 + a_lo*b_up*2^64)_up et les 128 bits inférieurs sont (a_up*b_lo*2^64 + a_lo*b_up*2^64 + a_lo*b_lo)_lo .

Il suffit ensuite de diviser ce résultat de 256 bits par c...

Answer 3

J'essayais juste de produire l'exemple ci-dessous lorsque le commentaire sur l'utilisation de l'option num-bigint est apparu... (rien de plus que le commentaire à expliquer)

// [dependencies]
// num-bigint = "0.4"

use num_bigint::BigUint;

fn muldiv(
    a: u128,
    b: u128,
    c: u128,
) -> u128 {
    let mut result = 0_u128;
    for (i, d) in (BigUint::from(a) * b / c).iter_u64_digits().enumerate() {
        match i {
            0 | 1 => result |= (d as u128) << (64 * i),
            _ => panic!("u128 overflow"),
        }
    }
    result
}

fn main() {
    let a = 123 * ((1_u128 << 100) + 1);
    let b = 456 * ((1_u128 << 100) + 1);
    let c = 789 * ((1_u128 << 100) + 1);
    let result = muldiv(a, b, c);
    println!("{:?} × {:?} / {:?} --> {:?}", a, b, c, result);
}