IRR

18 irr_fn "IRR" (1, 2) [R "values", V "guess"] |s, a, sv| {
    let mut values = resolved(s, a, 0, sv)?;
    // Fn7 scalar readings: held cells keep the collector's held-value rule,
    // even when the range is one cell. Direct text and logicals are ignored
    // too. Admit the guess before these scalar collections; a direct error
    // and larger arrays retain their existing validation order.
    let scalar = values.callable_value().is_none()
        && (matches!(a[0].core_ref(), CoreRef::Reference(_))
            && !matches!(values.core_ref(), CoreRef::Array(_))
            || matches!(values.core_ref(), CoreRef::Text(_) | CoreRef::Logical(_) | CoreRef::Empty));
    let mut g = None;
    if scalar {
        g = guess(s, a, 1, sv)?;
        values = CalcValue::array(CalcArray::from_scalar(values));
    }
    let cashflows = match cash::collect_irr_values_from_eval(&values) {
        Ok(cashflows) => cashflows,
        Err(error) => return cash_err(error),
    };
    if !scalar {
        g = guess(s, a, 1, sv)?;
    }
    out(s, cash::irr_kernel(&cashflows, g))
}

/// Excel `IRR` collects the values array/range with these rules (verified vs live
/// Excel 16.0 b20026, identical for array constants and ranges):
/// numbers are included; text, logical, and blank/empty cells are **ignored**;
/// and **any** error cell maps to `#VALUE!` (the specific code is not propagated).
/// A whole-argument scalar error still propagates as usual. This differs from the
/// strict positional collector used by XNPV/XIRR, whose values must stay aligned
/// with the parallel dates array.
pub(crate) fn collect_irr_values_from_eval(value: &CalcValue) -> Result<Vec<f64>, CashflowRateEvalError> {
    match value.core_ref() {
        CoreRef::Number(n) => Ok(vec![n]),
        CoreRef::Error(code) => Err(CashflowRateEvalError::Domain(code)),
        CoreRef::Array(array) => {
            let shape = array.shape();
            if shape.rows > 1 && shape.cols > 1 {
                return Err(CashflowRateEvalError::Domain(WorksheetErrorCode::Ref));
            }
            let mut out = Vec::with_capacity(shape.rows * shape.cols);
            for cell in array.iter_row_major() {
                match cell.core_ref() {
                    CoreRef::Number(n) => out.push(n),
                    CoreRef::Error(_) => {
                        return Err(CashflowRateEvalError::Domain(WorksheetErrorCode::Value));
                    }
                    CoreRef::Text(_)
                    | CoreRef::Logical(_)
                    | CoreRef::Empty
                    | CoreRef::Missing => {} // ignored, like Excel
                    CoreRef::Array(_) | CoreRef::Reference(_) => {
                        return Err(CashflowRateEvalError::Domain(WorksheetErrorCode::Value));
                    }
                }
            }
            Ok(out)
        }
        CoreRef::Text(_) | CoreRef::Logical(_) | CoreRef::Reference(_) => {
            Err(CashflowRateEvalError::Domain(WorksheetErrorCode::Value))
        }
        _ => Err(CashflowRateEvalError::Domain(WorksheetErrorCode::Value)),
    }
}

fn periodic_npv_with_t0(rate: f64, cashflows: &[f64]) -> Result<f64, WorksheetErrorCode> {
    validate_cashflows(cashflows)?;
    if !rate.is_finite() || rate <= MIN_VALID_RATE {
        return Err(WorksheetErrorCode::Num);
    }
    let base = 1.0 + rate;
    let mut total = 0.0;
    for (idx, cashflow) in cashflows.iter().enumerate() {
        total += *cashflow / crt_pow(base, idx as f64);
    }
    if total.is_finite() {
        Ok(total)
    } else {
        Err(WorksheetErrorCode::Num)
    }
}

fn periodic_npv_derivative(rate: f64, cashflows: &[f64]) -> Result<f64, WorksheetErrorCode> {
    validate_cashflows(cashflows)?;
    if !rate.is_finite() || rate <= MIN_VALID_RATE {
        return Err(WorksheetErrorCode::Num);
    }
    let base = 1.0 + rate;
    let mut total = 0.0;
    for (idx, cashflow) in cashflows.iter().enumerate().skip(1) {
        let period = idx as f64;
        total += -period * *cashflow / crt_pow(base, period + 1.0);
    }
    if total.is_finite() {
        Ok(total)
    } else {
        Err(WorksheetErrorCode::Num)
    }
}

fn ordered_rate(rate: f64) -> u64 {
    let bits = rate.to_bits();
    if rate.is_sign_negative() { !bits } else { bits | (1 << 63) }
}

fn rate_from_ordered(key: u64) -> f64 {
    f64::from_bits(if key & (1 << 63) == 0 { !key } else { key & !(1 << 63) })
}

/// Find the end of the residual plateau by doubling the distance, then bisecting its last
/// interval. Near zero, `1 + rate` can be constant over millions of rates: walking them one
/// by one is unbounded in practice. There are only 64 bits in the ordered representation,
/// so each end needs at most 64 distance probes and 64 bisections. A plateau with no finite
/// end returns #NUM!.
fn residual_plateau_end<F>(
    root: u64,
    residual: u64,
    lower: bool,
    f: &mut F,
) -> Result<u64, WorksheetErrorCode>
where
    F: FnMut(f64) -> Result<f64, WorksheetErrorCode>,
{
    let mut inside = root;
    let mut on_plateau = |rate| match f(rate) {
        Ok(value) => Ok(value.abs().to_bits() == residual),
        // A probe outside the NPV domain bounds the plateau just as a different residual
        // does. Bisect back into the valid interval instead of rejecting the root.
        Err(WorksheetErrorCode::Num) => Ok(false),
        Err(code) => Err(code),
    };
    for bit in 0..64 {
        let distance = 1u64 << bit;
        let outside = if lower { root.saturating_sub(distance) } else { root.saturating_add(distance) };
        let candidate = rate_from_ordered(outside);
        if candidate.is_finite() && on_plateau(candidate)? {
            inside = outside;
            continue;
        }
        let mut outside = outside;
        for _ in 0..64 {
            if inside.abs_diff(outside) <= 1 {
                // No finite end exists when even the last finite rate is on the plateau.
                if !rate_from_ordered(outside).is_finite() {
                    return Err(WorksheetErrorCode::Num);
                }
                return Ok(inside);
            }
            let midpoint = inside.min(outside) + inside.abs_diff(outside) / 2;
            let candidate = rate_from_ordered(midpoint);
            if candidate.is_finite() && on_plateau(candidate)? {
                inside = midpoint;
            } else {
                outside = midpoint;
            }
        }
        return Err(WorksheetErrorCode::Num);
    }
    Err(WorksheetErrorCode::Num)
}

fn midpoint_of_local_min_abs_residual_plateau<F>(
    root: f64,
    mut f: F,
) -> Result<f64, WorksheetErrorCode>
where
    F: FnMut(f64) -> Result<f64, WorksheetErrorCode>,
{
    let mut best_rate = root;
    let mut best_abs_bits = f(root)?.abs().to_bits();

    let mut lower_probe = root;
    let mut upper_probe = root;
    for _ in 0..IRR_PUBLICATION_ULP_SCAN_RADIUS {
        lower_probe = lower_probe.next_down();
        if lower_probe.is_finite() {
            let abs_bits = f(lower_probe)?.abs().to_bits();
            if abs_bits < best_abs_bits {
                best_abs_bits = abs_bits;
                best_rate = lower_probe;
            }
        }

        upper_probe = upper_probe.next_up();
        if upper_probe.is_finite() {
            let abs_bits = f(upper_probe)?.abs().to_bits();
            if abs_bits < best_abs_bits {
                best_abs_bits = abs_bits;
                best_rate = upper_probe;
            }
        }
    }

    let root = ordered_rate(best_rate);
    let lower = residual_plateau_end(root, best_abs_bits, true, &mut f)?;
    let upper = residual_plateau_end(root, best_abs_bits, false, &mut f)?;
    // Both zero representations are in the same plateau. Publish the exact zero root rather
    // than letting the different spacings on either side choose a subnormal rate.
    if rate_from_ordered(lower) <= 0.0 && rate_from_ordered(upper) >= 0.0 {
        return Ok(0.0);
    }
    Ok(rate_from_ordered(lower + (upper - lower) / 2))
}

/// Two-period cashflows have a single pole in q = 1 + rate. The recorded
/// rows favor a q-space secant fixpoint, including the final subtraction.
/// Equal residuals use the analytic Newton step, including an initial guess of
/// zero. A bounded solve keeps nonconvergence as the worksheet numeric error.
fn irr_two_cashflows(cashflows: &[f64], guess: f64) -> Result<f64, WorksheetErrorCode> {
    let mut current = 1.0 + guess;
    let mut previous = 1.0;
    let mut previous_value = periodic_npv_with_t0(0.0, cashflows)?;
    for _ in 0..20 {
        let rate = current - 1.0;
        let value = periodic_npv_with_t0(rate, cashflows)?;
        let step = if value == previous_value {
            let derivative = periodic_npv_derivative(rate, cashflows)?;
            if derivative == 0.0 {
                return Err(WorksheetErrorCode::Num);
            }
            value / derivative
        } else {
            value * (current - previous) / (value - previous_value)
        };
        let next = current - step;
        if next.to_bits() == current.to_bits() {
            return Ok(next - 1.0);
        }
        previous = current;
        previous_value = value;
        current = next;
    }
    Err(WorksheetErrorCode::Num)
}

pub(crate) fn irr_kernel(cashflows: &[f64], guess: Option<f64>) -> Result<f64, WorksheetErrorCode> {
    validate_cashflows(cashflows)?;
    if cashflows.len() == 2 {
        return irr_two_cashflows(cashflows, guess.unwrap_or(0.1));
    }
    let root = bounded_rate_solve(
        guess.unwrap_or(0.1),
        |rate| periodic_npv_with_t0(rate, cashflows),
        |rate| periodic_npv_derivative(rate, cashflows),
    )?;
    midpoint_of_local_min_abs_residual_plateau(root, |rate| periodic_npv_with_t0(rate, cashflows))
}