catalogue/rounding

/// An eager named call that takes from `min` to as many arguments as it names parameters.
const fn call(
    id: FunctionId,
    spelling: &'static str,
    min: usize,
    parameters: &'static [Parameter],
    kernel: EagerKernel,
) -> CatalogueEntry {
    CatalogueEntry {
        id,
        spelling,
        notation: Notation::Call,
        arity: CallableArityShape::range(min, parameters.len()),
        parameters,
        demand: ArgumentDemand::Eager,
        kernel: Kernel::Eager(kernel),
    }
}

/// Runs `compute` on the call's arguments, read as numbers in order, and publishes its number.
fn numeric(
    id: FunctionId,
    arguments: &[CalcValue],
    services: &EvalServices<'_>,
    compute: fn(&Site<'_>, &[f64]) -> Step<f64>,
) -> Result<CalcValue, FunctionFault> {
    let site = Site::new(id, services);
    let flush_input = matches!(id, TRUNC | ROUND | ROUNDUP | ROUNDDOWN | MOD);
    Site::finish(numbers(&site, arguments, flush_input).and_then(|read| {
        let result = compute(&site, &read)?;
        // Only `stored` makes a NaN here: Excel's own bit pattern for a result past DBL_MAX,
        // which the worksheet holds as a number.
        if result.is_nan() && matches!(id, TRUNC | ROUND | ROUNDDOWN) {
            Ok(CalcValue::excel_overflow_number(result).expect("the stored rounding result"))
        } else if result.is_subnormal() && matches!(id, TRUNC | ROUNDUP | ROUNDDOWN) {
            // Only the measured initial-reading return below can reach this publication.
            super::kit::out_with_subnormal(&site, Ok(result))
        } else {
            site.number(result)
        }
    }))
}

/// The arguments as numbers: an array argument is refused first, then each argument is read in
/// order, so the first one's error is the result.
fn numbers(site: &Site<'_>, arguments: &[CalcValue], flush_input: bool) -> Step<Vec<f64>> {
    site.scalar_operands(arguments)?;
    arguments
        .iter()
        .enumerate()
        .map(|(index, argument)| {
            // W111 typed MROUND/QUOTIENT readings: logicals are #VALUE!, an explicit
            // omission #N/A. A referenced blank remains Empty and coerces to zero.
            if matches!(site.id(), MROUND | QUOTIENT) {
                match argument.core_ref() {
                    CoreRef::Logical(_) => return Err(Exit::Error(WorksheetErrorCode::Value)),
                    CoreRef::Missing => return Err(Exit::Error(WorksheetErrorCode::NA)),
                    _ => {}
                }
            }
            if flush_input && argument.as_number().is_some_and(f64::is_subnormal) {
                Ok(0.0)
            } else {
                site.number_argument(index, argument)
            }
        })
        .collect()
}

/// ABS: the number without its sign.
fn abs(
    arguments: &[CalcValue],
    services: &mut EvalServices<'_>,
) -> Result<CalcValue, FunctionFault> {
    if arguments[0].is_excel_overflow() {
        let magnitude = arguments[0].as_number().expect("qualified overflow number").abs();
        return Ok(CalcValue::excel_overflow_number(magnitude).expect("sign removal retains overflow encoding"));
    }
    numeric(ABS, arguments, services, |_, n| Ok(n[0].abs()))
}

/// The fifteen significant digits the rounding functions round: the x87 reading where it is
/// known ([`ExactDecimal::fifteen_x87`]), and otherwise the number read to fifteen significant
/// digits, a half away from zero ([`ExactDecimal::fifteen`]).
fn reading(number: f64, ties_up: bool) -> ExactDecimal {
    ExactDecimal::fifteen_x87(number, ties_up).unwrap_or_else(|| ExactDecimal::fifteen(number).0)
}

/// A result below the normal range is 0, as Excel gives it (dc-bd5:
/// ROUNDDOWN(2.5E-308,308) is 0 in Excel 16.0.20430).
fn flush(result: f64) -> f64 {
    if result.is_subnormal() { 0.0 } else { result }
}

/// The reading itself as the result: the number unchanged, with no rounding, when the reading
/// is its exact value, and otherwise the reading converted at an arithmetic site.
fn whole_reading(site: &Site<'_>, number: f64, reading: &ExactDecimal) -> Step<f64> {
    if reading.is_value_of(number) {
        Ok(number)
    } else {
        Ok(flush(stored(site, reading)?))
    }
}

/// An exact decimal result as Excel stores it in binary64: the nearest binary64, except that a
/// result of 2^1024 or more keeps its 53-bit significand and carries its exponent into the
/// exponent field without a check, so the field reads 2047 and the bits are a NaN. So
/// ROUND(1.7976931348623157E+308,-308), 2E+308, is 0x7ff1ccf385ebc8a0 and the reading of
/// DBL_MAX, 1.79769313486232E+308, is 0x7ff000000000000a, as Excel 16.0.20430 and 16.0.20527 give
/// them (627 of 627 such bank rows). Excel shows such a value as 2E+308, ISNUMBER is TRUE, and
/// arithmetic on it is #NUM!. A decimal of 2^1025 or more, which no reading reaches, is #NUM!.
fn stored(site: &Site<'_>, decimal: &ExactDecimal) -> Step<f64> {
    let nearest = site.nearest(decimal)?;
    if nearest.is_finite() {
        return Ok(nearest);
    }
    // Half the decimal, exactly: the same significand one binade lower.
    let mut coefficient = decimal.coefficient.clone();
    coefficient.add(&decimal.coefficient.times(&decimal::Natural::from_u64(4)));
    let half = site.nearest(&ExactDecimal {
        negative: decimal.negative,
        coefficient,
        exponent: decimal.exponent - 1,
    })?;
    if half.is_finite() {
        Ok(f64::from_bits(half.to_bits() + (1 << 52)))
    } else {
        Err(Exit::Error(WorksheetErrorCode::Num))
    }
}

/// How ROUNDUP, ROUNDDOWN and TRUNC count num_digits: the 16 low bits of the integer part of its
/// magnitude, with its sign. A magnitude of 2^32 or more, or one whose 16 low bits reach 32768,
/// counts as 100 places either way and is out of range (false in the second place). So 32768 to
/// 65535 places count as 100, 65537 as 1, and -2147483748 as -100, as Excel 16.0.20430 gives
/// them (the OxFunc-Evidence bank's count-refinement and held-out runs).
fn directed_places(places: f64) -> (i32, bool) {
    let magnitude = places.abs();
    let low = if magnitude >= 4_294_967_296.0 {
        None
    } else {
        // Below 2^32, so the conversion is exact; the mask keeps 16 bits.
        let low = (decimal::truncate(magnitude) as u64) & 0xFFFF;
        (low < 0x8000).then(|| i32::try_from(low).expect("15 bits"))
    };
    let (count, in_range) = match low {
        Some(low) => (low, true),
        None => (100, false),
    };
    (if places < 0.0 { -count } else { count }, in_range)
}

/// ROUNDUP's, ROUNDDOWN's and TRUNC's kernel: the reading truncated at 10^-places, and for
/// ROUNDUP, when a digit was dropped, that truncation plus one unit of the place, the two
/// converted to binary64 and added at an arithmetic site (so ROUNDUP(-1.125,2) is
/// -1.1300000000000001 and ROUNDUP(3.535,2) 3.5399999999999996, as Excel gives them). With
/// num_digits out of range ROUNDUP adds no unit for a positive count and is #NUM! for a negative
/// one whenever a digit is dropped.
fn directed(site: &Site<'_>, number: f64, places: f64, direction: Direction) -> Step<f64> {
    if number == 0.0 {
        return Ok(0.0);
    }
    let (places, in_range) = directed_places(places);
    let position = (-places).clamp(-PLACES, PLACES);
    let reading = reading(number, true);
    if let Some(initial) = directed_initial_underflow(site, number, &reading)? {
        return Ok(initial);
    }
    let Some(down) = reading.round_at(position, Direction::TowardZero) else {
        let whole = whole_reading(site, number, &reading)?;
        // ROUNDUP gives #NUM! where ROUNDDOWN and TRUNC keep the bits past DBL_MAX (Excel
        // 16.0.20430: 102 of 102 bank rows, ROUNDUP(1.7976931348623153E+308,0) among them).
        return if whole.is_nan() && direction == Direction::AwayFromZero {
            Err(Exit::Error(WorksheetErrorCode::Num))
        } else {
            Ok(whole)
        };
    };
    let lower = stored(site, &down)?;
    let dropped = direction == Direction::AwayFromZero
        && reading
            .round_at(position, Direction::AwayFromZero)
            .is_some_and(|up| up.coefficient != down.coefficient);
    if !dropped {
        return Ok(flush(lower));
    }
    if !in_range {
        return if places > 0 {
            Ok(flush(lower))
        } else {
            Err(Exit::Error(WorksheetErrorCode::Num))
        };
    }
    if lower.is_nan() {
        // Past DBL_MAX, adding the unit is an overflow: #NUM!, as Excel gives it.
        return Err(Exit::Error(WorksheetErrorCode::Num));
    }
    let unit = site.nearest(&ExactDecimal::unit(position))?;
    // A unit below the normal range is not added, whichever the sign, and a positive number whose
    // reading lies less than 2^-1026 above the truncation, both as binary64, stays at it (Excel
    // 16.0.20527: ROUNDUP(1.0001086461844975E-307,307) is 1E-307, while
    // ROUNDUP(-1.0001086461844975E-307,307) is -2E-307). The negative rule has no such distance.
    if unit.is_subnormal() || (number > 0.0 && site.nearest(&reading)? - lower < SMALLEST_DISTANCE)
    {
        return Ok(flush(lower));
    }
    let magnitude = site.add(lower.abs(), unit)?;
    Ok(flush(if number < 0.0 { -magnitude } else { magnitude }))
}

/// A normal input whose initial15-digit reading falls below MIN_NORMAL returns before
/// digit rounding in the directed functions. Excel20430 and20527 retain its sign and the
/// normalized significand with exponent field zero, without shifting for a subnormal.
/// Doubling the exact decimal gives a normal value; decrementing its encoded exponent
/// reproduces that observed packing. This is an output model, not an internal-code claim.
/// The other directed paths still flush final subnormal results; ROUND keeps its own policy.
fn directed_initial_underflow(
    site: &Site<'_>,
    number: f64,
    reading: &ExactDecimal,
) -> Step<Option<f64>> {
    // Only the lowest normal binade can have a15-digit reading below MIN_NORMAL.
    // Keep arithmetic-free exact calls on their existing path, including under Exact mode.
    if number.abs() >= 2.0 * f64::MIN_POSITIVE || reading.is_value_of(number) {
        return Ok(None);
    }
    if !site.nearest(reading)?.is_subnormal() {
        return Ok(None);
    }
    let mut twice = reading.clone();
    twice.coefficient.add(&reading.coefficient);
    let normal = site.nearest(&twice)?;
    debug_assert!(normal.is_normal());
    Ok(Some(f64::from_bits(normal.to_bits() - (1 << 52))))
}

fn ceiling_for(
    id: FunctionId,
    arguments: &[CalcValue],
    services: &mut EvalServices<'_>,
) -> Result<CalcValue, FunctionFault> {
    numeric(id, arguments, services, |site, n| {
        let (number, significance) = (n[0], n[1]);
        if significance == 0.0 || number == 0.0 {
            Ok(0.0)
        } else if number > 0.0 && significance < 0.0 {
            Err(Exit::Error(WorksheetErrorCode::Num))
        } else {
            to_multiple(site, number, significance, Toward::Up)
        }
    })
}

/// The multiple of `significance` that `toward` chooses for `number`, both non-zero, in binary64:
/// significance * ceil(number / significance) for CEILING and significance *
/// floor(number / significance) for FLOOR, whatever the signs. The division and the product are
/// arithmetic sites, each rounding once; the integer above or below a binary64 is exact. The
/// quotient is first rounded to 15 significant digits, so one that rounds to just below or above
/// an integer counts as that integer (0.3 / 0.1 is 2.9999999999999996, counted as 3), and a
/// quotient beyond binary64's range gives an infinite product, which the number rule publishes
/// as #NUM!.
fn to_multiple(site: &Site<'_>, number: f64, significance: f64, toward: Toward) -> Step<f64> {
    let quotient = site.divide(number, significance)?;
    // Counted from the quotient rounded to 15 significant digits, as Excel counts (dc-j9k; see
    // ceiling_floor_family::excel_quotient).
    let quotient = if quotient.is_finite() && quotient != 0.0 {
        crate::ported::ceiling_floor_family::excel_quotient(number, significance)
    } else {
        quotient
    };
    let count = match toward {
        Toward::Up => decimal::ceil(quotient),
        Toward::Down => decimal::floor(quotient),
    };
    site.multiply(significance, count)
}

/// QUOTIENT: the integer part of the binary64 quotient numerator / denominator, toward zero. A
/// zero denominator is #DIV/0!.
fn quotient(
    arguments: &[CalcValue],
    services: &mut EvalServices<'_>,
) -> Result<CalcValue, FunctionFault> {
    numeric(QUOTIENT, arguments, services, |site, n| {
        let (numerator, denominator) = (n[0], n[1]);
        if denominator == 0.0 {
            return Err(Exit::Error(WorksheetErrorCode::Div0));
        }
        if numerator == 0.0 {
            return Ok(0.0);
        }
        // The binary64 quotient, toward zero (dc-bd5: QUOTIENT(0.3,0.1) is 2 in Excel
        // 16.0.20430, the quotient being 2.9999999999999996).
        let quotient = site.divide(numerator, denominator)?;
        Ok(decimal::truncate(quotient))
    })
}

/// ODD: rounded away from zero to an odd integer; an odd integer is unchanged, and 0 gives 1.
fn odd(
    arguments: &[CalcValue],
    services: &mut EvalServices<'_>,
) -> Result<CalcValue, FunctionFault> {
    numeric(ODD, arguments, services, |site, n| {
        to_parity(site, n[0], true)
    })
}

/// The integer of the wanted parity at or beyond |number|, with the number's sign. The integer
/// above |number| is exact; stepping on to the next integer is an addition at an arithmetic site,
/// which from 2^53 on, where every binary64 is even, rounds the odd integer to a binary64.
fn to_parity(site: &Site<'_>, number: f64, odd: bool) -> Step<f64> {
    if number == 0.0 {
        return Ok(if odd { 1.0 } else { 0.0 });
    }
    let magnitude = decimal::ceil(number.abs());
    // Halving an integer of at least 1 is exact.
    let half = magnitude * 0.5;
    let is_even = decimal::truncate(half) == half;
    let rounded = if is_even == odd {
        site.add(magnitude, 1.0)?
    } else {
        magnitude
    };
    Ok(if number < 0.0 { -rounded } else { rounded })
}