public fun max(a: u128, b: u128): u128 {
    if (a >= b) a else b
}

public fun min(a: u128, b: u128): u128 {
    if (a < b) a else b
}

public fun average(a: u128, b: u128): u128 {
    if (a < b) {
        a + (b - a) / 2
    } else {
        b + (a - b) / 2
    }
}

public inline fun gcd(a: u128, b: u128): u128 {
    let (large, small) = if (a > b) (a, b) else (b, a);
    while (small != 0) {
        let tmp = small;
        small = large % small;
        large = tmp;
    };
    large
}

public inline fun lcm(a: u128, b: u128): u128 {
    if (a == 0 || b == 0) {
        0
    } else {
        a / gcd(a, b) * b
    }
}

public inline fun mul_div(a: u128, b: u128, c: u128): u128 {
    // Inline functions cannot take constants, as then every module using it needs the constant
    assert!(c != 0, std::error::invalid_argument(4));
    (((a as u256) * (b as u256) / (c as u256)) as u128)
}

public fun clamp(x: u128, lower: u128, upper: u128): u128 {
    min(upper, max(lower, x))
}

public fun pow(n: u128, e: u128): u128 {
    if (e == 0) {
        1
    } else {
        let p = 1;
        while (e > 1) {
            if (e % 2 == 1) {
                p *= n;
            };
            e /= 2;
            n *= n;
        };
        p * n
    }
}

public fun floor_log2(x: u128): u8 {
    let res = 0;
    assert!(x != 0, std::error::invalid_argument(EINVALID_ARG_FLOOR_LOG2));
    // Effectively the position of the most significant set bit
    let n = 64;
    while ({
        spec {
            invariant n <= 64;
            invariant (res as u64) + 2 * (n as u64) <= 128;
            invariant n == 0 ==> (res as u64) <= 127;
        };
        n > 0
    }) {
        if (x >= (1 << n)) {
            x >>= n;
            res += n;
        };
        n >>= 1;
    };
    res
}

public fun log2(x: u128): FixedPoint32 {
    let integer_part = floor_log2(x);
    // Normalize x to [1, 2) in fixed point 32.
    if (x >= 1 << 32) {
        x >>= (integer_part - 32);
    } else {
        x <<= (32 - integer_part);
    };
    let frac = 0;
    let delta = 1 << 31;
    while (delta != 0) {
        // log x = 1/2 log x^2
        // x in [1, 2)
        x = (x * x) >> 32;
        // x is now in [1, 4)
        // if x in [2, 4) then log x = 1 + log (x / 2)
        if (x >= (2 << 32)) { frac += delta; x >>= 1; };
        delta >>= 1;
    };
    fixed_point32::create_from_raw_value (((integer_part as u64) << 32) + frac)
}

public fun log2_64(x: u128): FixedPoint64 {
    let integer_part = floor_log2(x);
    // Normalize x to [1, 2) in fixed point 63. To ensure x is smaller then 1<<64
    if (x >= 1 << 63) {
        x >>= (integer_part - 63);
    } else {
        x <<= (63 - integer_part);
    };
    let frac = 0;
    let delta = 1 << 63;
    while (delta != 0) {
        // log x = 1/2 log x^2
        // x in [1, 2)
        x = (x * x) >> 63;
        // x is now in [1, 4)
        // if x in [2, 4) then log x = 1 + log (x / 2)
        if (x >= (2 << 63)) { frac += delta; x >>= 1; };
        delta >>= 1;
    };
    fixed_point64::create_from_raw_value (((integer_part as u128) << 64) + frac)
}

public fun sqrt(x: u128): u128 {
    if (x == 0) return 0;
    // Note the plus 1 in the expression. Let n = floor_lg2(x) we have x in [2^n, 2^{n+1}) and thus the answer in
    // the half-open interval [2^(n/2), 2^{(n+1)/2}). For even n we can write this as [2^(n/2), sqrt(2) 2^{n/2})
    // for odd n [2^((n+1)/2)/sqrt(2), 2^((n+1)/2). For even n the left end point is integer for odd the right
    // end point is integer. If we choose as our first approximation the integer end point we have as maximum
    // relative error either (sqrt(2) - 1) or (1 - 1/sqrt(2)) both are smaller then 1/2.
    let res = 1 << ((floor_log2(x) + 1) >> 1);
    // We use standard newton-rhapson iteration to improve the initial approximation.
    // The error term evolves as delta_i+1 = delta_i^2 / 2 (quadratic convergence).
    // It turns out that after 5 iterations the delta is smaller than 2^-64 and thus below the treshold.
    res = (res + x / res) >> 1;
    res = (res + x / res) >> 1;
    res = (res + x / res) >> 1;
    res = (res + x / res) >> 1;
    res = (res + x / res) >> 1;
    min(res, x / res)
}

public inline fun ceil_div(x: u128, y: u128): u128 {
    // ceil_div(x, y) = floor((x + y - 1) / y) = floor((x - 1) / y) + 1
    // (x + y - 1) could spuriously overflow. so we use the later version
    if (x == 0) {
        // Inline functions cannot take constants, as then every module using it needs the constant
        assert!(y != 0, std::error::invalid_argument(4));
        0
    }
    else (x - 1) / y + 1
}