Raise to the power of an _integer_

[drinkcat]

Slightly simpler than #45 and #74.

I'd also need pow function for uutils (https://github.com/uutils/coreutils/issues/7708), but only for integer exponents, and I want a fixed limited precision.

The use case is hexadecimal floats, where to get the value for a number like 0x1p100000, we need to compute 0x1 * 2**100000). My first naive solution would do the operation on a BigUint, but the number of digits totally explodes for large exponents.

Ended up using exponentiation by squaring basic algorithm from Wikipedia: https://en.wikipedia.org/wiki/Exponentiation_by_squaring . Came up with this basic function for now:

/// Compute bd**exp using exponentiation by squaring algorithm, while maintaining the
/// precision specified in ctx (the number of digits would otherwise explode).
// TODO: We do lose a little bit of precision, and the last digits may not be correct.
// TODO: Upstream this to bigdecimal-rs.
fn pow_with_context(bd: BigDecimal, exp: u32, ctx: &bigdecimal::Context) -> BigDecimal {
    if exp == 0 {
        return 1.into();
    }

    fn trim_precision(bd: BigDecimal, ctx: &bigdecimal::Context) -> BigDecimal {
        if bd.digits() > ctx.precision().get() {
            bd.with_precision_round(ctx.precision(), ctx.rounding_mode())
        } else {
            bd
        }
    }

    let bd = trim_precision(bd, ctx);
    let ret = if exp % 2 == 0 {
        pow_with_context(bd.square(), exp / 2, ctx)
    } else {
        &bd * pow_with_context(bd.square(), (exp - 1) / 2, ctx)
    };
    trim_precision(ret, ctx)
}

I'll try to refine this and submit a PR here. Still a bit unclear how to get the exact required precision (e.g. 100 digits):

// Input
0x1p3000000000
// Wolfram Alpha says
9.8162042336235053508313854078782835648991393286913072670026492205522618203568834202759669215027003865...e+903089986
// We get:               9.816204233623505350831385407878283564899139328691307267002649220552261820356883420275966921514831318e+903089986

I could avoid losing precision by rounding intermediate values to 2*ctx.precision(), but that feels like a big hammer.