ported/power_fn

fn exact_integer_magnitude(power: f64) -> Option<u32> {
    if !power.is_finite() {
        return None;
    }
    let magnitude = power.abs();
    // W111 width controls distinguish the unsigned32 sentinel itself from
    // its predecessor. Out-of-range integers use the noninteger body below.
    if magnitude.fract() != 0.0 || magnitude >= u32::MAX as f64 {
        None
    } else {
        Some(magnitude as u32)
    }
}

/// Excel's real-root acceptance for a negative base with a fractional exponent.
/// Excel accepts `POWER(x<0, y)` only when `1/y`, rounded to **15 significant
/// decimal digits**, is an odd integer within signed `int32`; the result is then
/// `-POWER(|x|, y)`, evaluated at the ORIGINAL `y`. (Round-1 live-Excel finding:
/// the acceptance band around `n=101` is exactly 100x the band around `n=3` — a
/// decimal signature no binary tolerance reproduces. BUG-FUNC-042 / W108.)
fn accept_negative_base_odd_root(power: f64) -> bool {
    // 15 significant decimal digits == 14 after the point in scientific form,
    // matching Excel's `%.14e` fuzz on 1/y.
    let n: f64 = match format!("{:.14e}", 1.0 / power).parse() {
        Ok(v) => v,
        Err(_) => return false,
    };
    n.fract() == 0.0 && (-2_147_483_648.0..=2_147_483_647.0).contains(&n) && (n as i64) % 2 != 0
}

/// LSB-first (right-to-left) square-and-multiply `base^n` for `n >= 0`, with
/// staged products followed by binary64 stores (W111 numerical observations,
/// distinct from `powf`). Returns the POSITIVE power; the caller reciprocates.
fn binexp_positive(base: f64, mut n: u64) -> f64 {
    let mut acc = 1.0;
    let mut b = base;
    while n > 0 {
        if n & 1 == 1 {
            acc = crate::ported::excel_numeric::excel_x87_mul(acc, b);
        }
        n >>= 1;
        if n > 0 {
            b = crate::ported::excel_numeric::excel_x87_mul(b, b);
        }
    }
    acc
}

/// `POWER(x, y)` candidate graph (the `^` operator shares this kernel).
/// Extends the black-box BUG-FUNC-042 / W108 graph with W111 integer-power
/// observations: positive base ten has a decimal-scale branch; other integer
/// exponents use staged binary exponentiation. Fractional exponents use
/// → `exp(y·ln x)` via the x87 backend with the `y<0` reciprocal staging (positive
/// power stored to `binary64`, then one x87 double-rounded reciprocal) that
/// resolved the earlier 1-ULP residual — `exp(y·ln x)` was the right function
/// evaluated at the wrong point for `y < 0`.
pub fn power_kernel(number: f64, power: f64) -> Result<f64, WorksheetErrorCode> {
    // 1. zero exponent
    if power == 0.0 {
        return if number == 0.0 {
            Err(WorksheetErrorCode::Num) // POWER(0, 0)
        } else {
            Ok(1.0)
        };
    }
    // 2. zero base
    if number == 0.0 {
        return if power < 0.0 {
            Err(WorksheetErrorCode::Div0) // POWER(0, negative)
        } else {
            Ok(0.0)
        };
    }

    // 3. integer exponent → binary-exponentiation publication. This is the
    // declared precision quirk on `POWER_META`, also used by `^`. Consumers
    // retain their separately assessed call-site graphs. Dispatch is EXACT —
    // `POWER(2, 3+ulp)` takes the fractional
    // path. The positive power is formed by `binexp`, reciprocated for `y<0`.
    // POWER_META declares PrecisionRoundingProfile::IntegerExponentPublication.
    let integer_publication = true;
    if let Some(magnitude) = exact_integer_magnitude(power).filter(|_| integer_publication) {
        // This branch is specific to positive ten; negative ten retains the
        // ordinary repeated-product graph and differs at many large powers.
        let mut p = if number == 10.0 {
            // The observed decimal scale wraps through signed32 inside the
            // admitted unsigned32 magnitude. E.g. 2^32-2 selects scale -2.
            let scale = magnitude as i32;
            if scale > 308 {
                f64::INFINITY
            } else if scale < -308 {
                0.0
            } else {
                format!("1e{scale}").parse::<f64>().unwrap_or(f64::INFINITY)
            }
        } else {
            binexp_positive(number, u64::from(magnitude))
        };
        if power < 0.0 {
            if p.abs() < f64::MIN_POSITIVE {
                return Err(WorksheetErrorCode::Div0); // reciprocal sees tiny magnitude as zero
            }
            p = excel_x87_recip(p); // 1/inf -> 0, published below
        }
        if p.is_nan() || p.is_infinite() {
            return Err(WorksheetErrorCode::Num);
        }
        // Integer path flushes subnormal / -0 / 0 to +0 (Excel: POWER(2,-1023)=0).
        return Ok(if p.abs() < f64::MIN_POSITIVE { 0.0 } else { p });
    }

    // 4. fractional exponent.
    let mut sign = 1.0;
    let mut base = number;
    if number < 0.0 {
        if !accept_negative_base_odd_root(power) {
            return Err(WorksheetErrorCode::Num); // no real root
        }
        sign = -1.0;
        base = -number; // continue with |x|
    }

    // p = exp_x87( RN53(RN64( |y| * ln_x87(base) )) ) — the positive power.
    let p = excel_pow_positive(base, power.abs());

    if power > 0.0 {
        if p.is_infinite() {
            return Err(WorksheetErrorCode::Num); // overflow
        }
        if p == 0.0 {
            return Ok(0.0);
        }
        return Ok(sign * p); // subnormals ARE published on the y>0 path
    }

    // y < 0: reciprocal staging (the resolved 1-ULP puzzle).
    if p < f64::MIN_POSITIVE {
        return Err(WorksheetErrorCode::Div0); // subnormal/zero positive power (DAZ-style)
    }
    if p.is_infinite() {
        return Ok(0.0); // 1 / inf
    }
    let r = excel_x87_recip(p); // RN53(RN64(1/p))
    if r < f64::MIN_POSITIVE {
        return Ok(0.0); // subnormal quotient flushes to +0
    }
    if r.is_infinite() {
        return Err(WorksheetErrorCode::Num);
    }
    Ok(sign * r)
}