ported/regression_forecast_family

fn checked(low: f64, high: f64) -> Option<Self> {
    (low.is_finite() && high.is_finite() && low <= high).then_some(Self { low, high })
}

fn subtract(self, other: Self) -> Option<Self> {
    Self::checked((self.low-other.high).next_down(), (self.high-other.low).next_up())
}

fn multiply(self, other: Self) -> Option<Self> {
    let products = [self.low*other.low, self.low*other.high, self.high*other.low, self.high*other.high];
    if products.iter().any(|v| !v.is_finite()) { return None; }
    Self::checked(products.into_iter().fold(f64::INFINITY, f64::min).next_down(),
        products.into_iter().fold(f64::NEG_INFINITY, f64::max).next_up())
}

fn square(self) -> Option<Self> {
    let a = self.low*self.low;
    let b = self.high*self.high;
    let low = if self.low <= 0.0 && self.high >= 0.0 { 0.0 } else { a.min(b).next_down().max(0.0) };
    Self::checked(low, a.max(b).next_up())
}

fn divide_positive(self, other: Self) -> Option<Self> {
    if other.low <= 0.0 { return None; }
    self.multiply(Self::checked((1.0/other.high).next_down(), (1.0/other.low).next_up())?)
}

fn moment(&self, a: usize, sa: f64, b: Option<(usize, f64)>) -> Option<GramInterval> {
    let mut sum = 0.0;
    let mut absolute = 0.0;
    for (i, &raw) in self.columns[a].iter().enumerate() {
        let x = raw/sa;
        let y = b.map_or(1.0, |(index, scale)| self.columns[index][i]/scale);
        let product = x*y;
        // The relative-error bound requires normal nonzero products.
        if x != 0.0 && y != 0.0 && !product.is_normal() { return None; }
        sum += product;
        absolute += product.abs();
    }
    if absolute == 0.0 { return GramInterval::checked(sum, sum); }
    let denominator = (1.0-self.gamma).next_down();
    if denominator <= 0.0 { return None; }
    let absolute_upper = (absolute/denominator).next_up();
    let error = (self.gamma*absolute_upper).next_up();
    GramInterval::checked((sum-error).next_down(), (sum+error).next_up())
}

fn column(&mut self, index: usize) -> Option<GramColumn> {
    if let Some(stats) = self.stats.get(&index) { return Some(*stats); }
    let values = &self.columns[index];
    if values.len() as f64 != self.count || values.iter().any(|v| !v.is_finite()) { return None; }
    let largest = values.iter().map(|v| v.abs()).fold(0.0_f64, f64::max);
    if largest == 0.0 { return None; }
    // Power-of-two normalization is exact when the scaled values stay
    // normal. A representation failure keeps the existing QR path.
    let exponent = largest.to_bits() & 0x7ff0_0000_0000_0000;
    let scale = f64::from_bits(exponent.max(f64::MIN_POSITIVE.to_bits()));
    if values.iter().any(|v| *v != 0.0 && !(v/scale).is_normal()) { return None; }
    let stats = GramColumn { scale, sum: self.moment(index, scale, None)? };
    self.stats.insert(index, stats);
    Some(stats)
}

fn entry(&mut self, a: usize, b: usize) -> Option<GramInterval> {
    let key = (a.min(b), a.max(b));
    if let Some(value) = self.entries.get(&key) { return Some(*value); }
    let x = self.column(a)?;
    let y = self.column(b)?;
    let mut value = self.moment(a, x.scale, Some((b, y.scale)))?;
    if self.centered {
        let count = GramInterval { low: self.count, high: self.count };
        value = value.subtract(x.sum.multiply(y.sum)?.divide_positive(count)?)?;
    }
    self.entries.insert(key, value);
    Some(value)
}

fn factor_column(&mut self, index: usize) -> Option<(Vec<GramInterval>, GramInterval)> {
    let mut row: Vec<GramInterval> = Vec::with_capacity(self.basis.len());
    for k in 0..self.basis.len() {
        let mut value = self.entry(index, self.basis[k])?;
        for r in 0..k {
            value = value.subtract(row[r].multiply(self.lower[k][r])?.multiply(self.pivots[r])?)?;
        }
        row.push(value.divide_positive(self.pivots[k])?);
    }
    let mut pivot = self.entry(index, index)?;
    for (r, coefficient) in row.iter().enumerate() {
        pivot = pivot.subtract(coefficient.square()?.multiply(self.pivots[r])?)?;
    }
    Some((row, pivot))
}

fn unresolved_inner(&mut self, prefix: &[usize], index: usize) -> Option<bool> {
    if self.basis != prefix {
        self.basis.clear(); self.lower.clear(); self.pivots.clear();
        for &earlier in prefix {
            let (row, pivot) = self.factor_column(earlier)?;
            if pivot.low <= 0.0 { return None; }
            self.basis.push(earlier); self.lower.push(row); self.pivots.push(pivot);
        }
    }
    let (row, pivot) = self.factor_column(index)?;
    if pivot.high < 0.0 { return None; }
    if pivot.low <= 0.0 { return Some(true); }
    self.basis.push(index); self.lower.push(row); self.pivots.push(pivot);
    Some(false)
}

fn unresolved(&mut self, prefix: &[usize], index: usize) -> bool {
    if self.disabled || prefix.is_empty() { return false; }
    match self.unresolved_inner(prefix, index) {
        Some(result) => result,
        None => { self.disabled = true; false }
    }
}

/// An argument as the regression functions read it: the caller has read a reference into an
/// array; an empty or omitted argument is #VALUE!.
fn resolve_arg_eval(arg: &CalcValue) -> Result<CalcValue, WorksheetErrorCode> {
    match arg.core_ref() {
        CoreRef::Missing | CoreRef::Empty => Err(WorksheetErrorCode::Value),
        _ => Ok(arg.clone()),
    }
}

fn optional_arg_value(
    arg: Option<&CalcValue>
) -> Result<Option<CalcValue>, WorksheetErrorCode> {
    match arg {
        None => Ok(None),
        Some(value) if matches!(value.core_ref(), CoreRef::Missing) => Ok(None),
        Some(other) => Ok(Some(resolve_arg_eval(other)?)),
    }
}

fn numeric_from_cell(cell: &CalcValue) -> Result<f64, WorksheetErrorCode> {
    match cell.core_ref() {
        CoreRef::Number(n) => Ok(n),
        CoreRef::Error(code) => Err(code),
        CoreRef::Text(_)
        | CoreRef::Logical(_)
        | CoreRef::Empty
        | CoreRef::Missing
        | CoreRef::Array(_)
        | CoreRef::Reference(_) => Err(WorksheetErrorCode::Value,),
    }
}

fn numeric_scalar_from_eval(value: &CalcValue) -> Result<f64, WorksheetErrorCode> {
    match value.core_ref() {
        CoreRef::Number(n) => Ok(n),
        CoreRef::Logical(flag) => Ok(if flag { 1.0 } else { 0.0 }),
        CoreRef::Error(code) => Err(code),
        CoreRef::Array(array) => {
            let shape = array.shape();
            if shape.rows == 1 && shape.cols == 1 {
                numeric_from_cell(array.get(0, 0).expect("single-cell array"))
            } else {
                Err(WorksheetErrorCode::Value,)
            }
        }
        CoreRef::Text(_) | CoreRef::Reference(_) => Err(WorksheetErrorCode::Value,),
        _ => Err(WorksheetErrorCode::Value,),
    }
}

fn bool_flag_from_eval_or_default(
    value: Option<CalcValue>,
    default: bool,
) -> Result<bool, WorksheetErrorCode> {
    match value {
        None => Ok(default),
        Some(value) if matches!(value.core_ref(), CoreRef::Logical(_)) => match value.core_ref() {
            CoreRef::Logical(flag) => Ok(flag),
            _ => unreachable!("guarded logical"),
        },
        Some(other) => Ok(numeric_scalar_from_eval(&other)? != 0.0),
    }
}

fn collect_numeric_vector_from_eval(
    value: &CalcValue,
) -> Result<NumericVector, WorksheetErrorCode> {
    match value.core_ref() {
        CoreRef::Number(n) => Ok(NumericVector {
            values: vec![n],
            shape: VectorShape::Scalar,
        }),
        CoreRef::Error(code) => Err(code),
        CoreRef::Array(array) => {
            let shape = array.shape();
            if shape.rows > 1 && shape.cols > 1 {
                return Err(WorksheetErrorCode::Ref);
            }
            let mut values = Vec::with_capacity(shape.rows * shape.cols);
            for cell in array.iter_row_major() {
                values.push(numeric_from_cell(cell)?);
            }
            let vector_shape = if shape.rows == 1 && shape.cols == 1 {
                VectorShape::Scalar
            } else if shape.rows == 1 {
                VectorShape::Row(shape.cols)
            } else {
                VectorShape::Column(shape.rows)
            };
            Ok(NumericVector {
                values,
                shape: vector_shape,
            })
        }
        CoreRef::Text(_) | CoreRef::Logical(_) | CoreRef::Reference(_) => Err(
            WorksheetErrorCode::Value,
        ),
        _ => Err(WorksheetErrorCode::Value,),
    }
}

fn collect_numeric_vector_arg(
    arg: &CalcValue
) -> Result<NumericVector, WorksheetErrorCode> {
    let eval = resolve_arg_eval(arg)?;
    collect_numeric_vector_from_eval(&eval)
}

fn cell(&self, row: usize, col: usize) -> f64 {
    self.values[row * self.cols + col]
}

fn transpose(&self) -> Self {
    let mut values = vec![0.0; self.rows * self.cols];
    for row in 0..self.rows {
        for col in 0..self.cols {
            values[col * self.rows + row] = self.cell(row, col);
        }
    }
    Self {
        rows: self.cols,
        cols: self.rows,
        values,
    }
}

fn row_values(&self, row: usize) -> &[f64] {
    let start = row * self.cols;
    &self.values[start..(start + self.cols)]
}

fn collect_numeric_matrix_from_eval(
    value: &CalcValue,
) -> Result<NumericMatrix, WorksheetErrorCode> {
    match value.core_ref() {
        CoreRef::Number(n) => Ok(NumericMatrix {
            rows: 1,
            cols: 1,
            values: vec![n],
        }),
        CoreRef::Error(code) => Err(code),
        CoreRef::Array(array) => {
            let shape = array.shape();
            let mut values = Vec::with_capacity(shape.rows * shape.cols);
            for cell in array.iter_row_major() {
                values.push(numeric_from_cell(cell)?);
            }
            Ok(NumericMatrix {
                rows: shape.rows,
                cols: shape.cols,
                values,
            })
        }
        CoreRef::Text(_) | CoreRef::Logical(_) | CoreRef::Reference(_) => Err(
            WorksheetErrorCode::Value,
        ),
        _ => Err(WorksheetErrorCode::Value,),
    }
}

fn collect_numeric_matrix_arg(
    arg: &CalcValue
) -> Result<NumericMatrix, WorksheetErrorCode> {
    let eval = resolve_arg_eval(arg)?;
    collect_numeric_matrix_from_eval(&eval)
}

fn sequence_predictors(len: usize) -> NumericMatrix {
    NumericMatrix {
        rows: len,
        cols: 1,
        values: (1..=len).map(|n| n as f64).collect(),
    }
}

fn default_output_shape_for_known_y(shape: VectorShape) -> OutputShape {
    match shape {
        VectorShape::Scalar => OutputShape::Scalar,
        VectorShape::Row(cols) => OutputShape::Row(cols),
        VectorShape::Column(rows) => OutputShape::Column(rows),
    }
}

fn normalize_known_x_matrix(
    raw: NumericMatrix,
    observation_count: usize,
) -> Result<NumericMatrix, WorksheetErrorCode> {
    if raw.rows == observation_count {
        Ok(raw)
    } else if raw.rows == 1 && raw.cols == observation_count {
        Ok(raw.transpose())
    } else if observation_count == 1 {
        // One known y: every known x is a predictor of that one observation (Excel 16.0.20430,
        // 2026-10-25, dc-6au: `=LINEST(1,{1;2;3})` and `=LINEST(1,{1,2,3})` spill 0,0,0,1, three
        // zero slopes and the intercept; `=TREND(1,{1;2;3})` and `=GROWTH(1,{1;2;3})` are 1,
        // `=LOGEST(1,{1;2;3})` 1,1,1,1; `=TREND(5,{1;2;3},4)` is #REF!).
        Ok(NumericMatrix {
            rows: 1,
            cols: raw.rows * raw.cols,
            values: raw.values,
        })
    } else {
        Err(WorksheetErrorCode::Ref)
    }
}

fn prepare_predictor_input(
    raw: Option<NumericMatrix>,
    known_y: &NumericVector,
    treat_single_predictor_cells_as_observations: bool,
) -> Result<PreparedPredictorInput, WorksheetErrorCode> {
    let observation_count = known_y.values.len();
    match raw {
        None => Ok(PreparedPredictorInput {
            predictors: sequence_predictors(observation_count),
            output_shape: default_output_shape_for_known_y(known_y.shape),
        }),
        Some(matrix) => {
            let normalized = normalize_known_x_matrix(matrix.clone(), observation_count)?;
            if treat_single_predictor_cells_as_observations && normalized.cols == 1 {
                let output_shape = if matrix.rows == 1 && matrix.cols == 1 {
                    OutputShape::Scalar
                } else if matrix.rows == 1 {
                    OutputShape::Row(matrix.cols)
                } else if matrix.cols == 1 {
                    OutputShape::Column(matrix.rows)
                } else {
                    OutputShape::Matrix {
                        rows: matrix.rows,
                        cols: matrix.cols,
                    }
                };
                return Ok(PreparedPredictorInput {
                    predictors: NumericMatrix {
                        rows: matrix.rows * matrix.cols,
                        cols: 1,
                        values: matrix.values,
                    },
                    output_shape,
                });
            }

            if normalized.cols == 1 {
                let output_shape = if matrix.rows == 1 && matrix.cols == 1 {
                    OutputShape::Scalar
                } else if matrix.rows == 1 {
                    OutputShape::Row(matrix.cols)
                } else {
                    OutputShape::Column(matrix.rows)
                };
                return Ok(PreparedPredictorInput {
                    predictors: normalized,
                    output_shape,
                });
            }

            if matrix.cols != normalized.cols && observation_count != 1 {
                return Err(WorksheetErrorCode::Ref);
            }
            let output_shape = if matrix.rows == 1 || observation_count == 1 {
                OutputShape::Scalar
            } else {
                OutputShape::Column(matrix.rows)
            };
            Ok(PreparedPredictorInput {
                predictors: normalized,
                output_shape,
            })
        }
    }
}

fn predictor_input_from_optional(
    raw: Option<NumericMatrix>,
    known_y: &NumericVector,
) -> Result<PreparedPredictorInput, WorksheetErrorCode> {
    prepare_predictor_input(raw, known_y, false)
}

fn prediction_input_from_optional(
    raw: Option<NumericMatrix>,
    known_y: &NumericVector,
    known_x: &PreparedPredictorInput,
) -> Result<PreparedPredictorInput, WorksheetErrorCode> {
    match raw {
        None => Ok(PreparedPredictorInput {
            predictors: known_x.predictors.clone(),
            output_shape: default_output_shape_for_known_y(known_y.shape),
        }),
        Some(matrix) => {
            let predictor_count = known_x.predictors.cols;
            if predictor_count == 1 {
                let output_shape = if matrix.rows == 1 && matrix.cols == 1 {
                    OutputShape::Scalar
                } else if matrix.rows == 1 {
                    OutputShape::Row(matrix.cols)
                } else if matrix.cols == 1 {
                    OutputShape::Column(matrix.rows)
                } else {
                    OutputShape::Matrix {
                        rows: matrix.rows,
                        cols: matrix.cols,
                    }
                };
                return Ok(PreparedPredictorInput {
                    predictors: NumericMatrix {
                        rows: matrix.rows * matrix.cols,
                        cols: 1,
                        values: matrix.values,
                    },
                    output_shape,
                });
            }

            if matrix.cols != predictor_count {
                return Err(WorksheetErrorCode::Ref);
            }
            let output_shape = if matrix.rows == 1 {
                OutputShape::Scalar
            } else {
                OutputShape::Column(matrix.rows)
            };
            Ok(PreparedPredictorInput {
                predictors: matrix,
                output_shape,
            })
        }
    }
}

/// The least-squares fit of LINEST, LOGEST, TREND and GROWTH as Excel computes it (dc-j9k,
/// identified on captured coefficients and predictions of one, two and three predictors, with and
/// without a constant; the catalogue oracle's sheet dcj9k_regression holds them):
///
/// 1. With a constant, each predictor column and `y` are centred on their means (each mean the
///    forward sum divided by the count), and a column of ones leads the matrix; without one, the
///    raw columns are the matrix.
/// 2. Householder QR. Before each predictor column's step, the predictor column with the largest
///    sum of squares over the remaining rows comes next (the first such column on a tie); the
///    ones column always goes first. At step `j`, with `x0` the diagonal entry, `rest` the sum
///    of squares below it and `s` the sum of squares from it down (each a forward sum from 0):
///    `alpha = -sign(x0)*sqrt(s)`, `v = (x0 - alpha, the entries below)`, `vv = rest + v0*v0`, and
///    each later column and `y` become `a - (2*(v.a)/vv)*v`.
/// 3. Back substitution: `b_j = (y_j - d)/alpha_j`, with `d` the forward sum of `r_jk*b_k` over
///    the later columns.
/// 4. The intercept is `my - b_1*m_1 - b_2*m_2 - ...`, subtracted in predictor order.
///
/// 5. Exactly proportional predictors are removed after the larger column has
///    been selected. Zero and surplus predictors retain their structural handling.
///    A further predictor is numerically eliminated when outward Gram-entry
///    bounds leave its LDL pivot unresolved. Coefficients still use this QR.
fn least_squares(
    columns: &[Vec<f64>],
    y: &[f64],
    use_const: bool,
) -> Result<RegressionFit, WorksheetErrorCode> {
    let n = y.len();
    let mean = |values: &[f64]| values.iter().fold(0.0, |sum, v| sum + v) / n as f64;
    let means: Vec<f64> = columns.iter().map(|column| mean(column)).collect();
    let y_mean = if use_const { mean(y) } else { 0.0 };
    // The matrix's columns, each with the predictor it holds (None for the ones column).
    let mut matrix: Vec<(Option<usize>, Vec<f64>)> = Vec::with_capacity(columns.len() + 1);
    if use_const {
        matrix.push((None, vec![1.0; n]));
    }
    for (index, column) in columns.iter().enumerate() {
        let values = if use_const {
            column.iter().map(|v| v - means[index]).collect()
        } else {
            column.clone()
        };
        matrix.push((Some(index), values));
    }
    let mut target: Vec<f64> = if use_const {
        y.iter().map(|v| v - y_mean).collect()
    } else {
        y.to_vec()
    };
    let centred_y = target.clone();
    let p = matrix.len();
    if n == 1 && use_const {
        // One observation with a constant: every centred predictor is 0, so each coefficient is
        // 0 and the intercept is the one y (Excel 16.0.20430, 2026-10-25, dc-6au:
        // `=LINEST(1,{1;2;3})` 0,0,0,1, `=LINEST(5,2)` 0,5, `=TREND(1,{1;2;3})` 1).
        return Ok(RegressionFit {
            model: LinearRegressionModel {
                coefficients: vec![0.0; columns.len()],
                intercept: y_mean,
            },
            order: Vec::new(),
            r: Vec::new(),
            qty: Vec::new(),
            means,
            y: Vec::new(),
            observations: n,
        });
    }
    // R25 full-rank orthogonal probes require PC64 square-product arithmetic.
    // Use the shared stored product; the exact accumulator/store layout stays
    // unresolved, so keep the existing addition.
    let squares_from = |values: &[f64], from: usize| {
        values[from..].iter().fold(0.0, |sum, v| sum + excel_numeric::excel_x87_mul(*v, *v))
    };
    // Each column's sum of squares before the steps, the scale its remainder is judged against.
    let initial: Vec<f64> = matrix.iter().map(|(_, column)| squares_from(column, 0)).collect();
    let mut scale: Vec<(Option<usize>, f64)> =
        matrix.iter().map(|(predictor, _)| *predictor).zip(initial).collect();
    let finite_norms = scale.iter().all(|(predictor, norm)| predictor.is_none() || norm.is_finite());
    let mut gram_rank = GramRank::new(columns, use_const, n, finite_norms);
    let mut diagonal = vec![0.0; p];
    // Pivot by the remaining sum of squares; drop an exactly eliminated column
    // only when it was nonzero before the reflections. Initially zero columns
    // retain their zero coefficient and count in the degrees of freedom.
    let mut j = 0;
    while j < matrix.len() {
        if j >= n {
            matrix.truncate(j);
            scale.truncate(j);
            break;
        }
        if matrix[j].0.is_some() {
            let mut best = j;
            let mut best_squares = -1.0;
            for (k, (_, column)) in matrix.iter().enumerate().skip(j) {
                let squares = squares_from(column, j);
                if squares > best_squares {
                    best_squares = squares;
                    best = k;
                }
            }
            matrix.swap(j, best);
            scale.swap(j, best);
        }
        let x0 = matrix[j].1[j];
        let rest = squares_from(&matrix[j].1, j + 1);
        let squares = squares_from(&matrix[j].1, j);
        // Reflections can leave a nonzero floating-point remainder even when
        // two original columns are exactly proportional. Test that dependence
        // in the original binary inputs, without a numerical rank tolerance.
        let proportional = scale[j].1 > 0.0 && matrix[j].0.is_some_and(|index| {
            matrix[..j].iter().zip(&scale[..j]).any(|((earlier, _), (_, original_squares))| {
                *original_squares > 0.0
                    && earlier.is_some_and(|other| exactly_proportional(&columns[index], &columns[other]))
            })
        });
        if proportional {
            matrix.remove(j);
            scale.remove(j);
            continue;
        }
        let unresolved = scale[j].1 > 0.0 && matrix[j].0.is_some_and(|index| {
            let prefix: Vec<usize> = matrix[..j].iter().zip(&scale[..j])
                .filter_map(|((predictor, _), (_, norm))| if *norm > 0.0 { *predictor } else { None })
                .collect();
            gram_rank.unresolved(&prefix, index)
        });
        if scale[j].1 > 0.0 && matrix[j].0.is_some() && (squares == 0.0 || unresolved) {
            matrix.remove(j);
            scale.remove(j);
            continue;
        }
        if squares == 0.0 {
            diagonal[j] = 0.0;
            j += 1;
            continue;
        }
        let norm = squares.sqrt();
        let alpha = -norm.copysign(x0);
        let mut v = matrix[j].1[j..].to_vec();
        // Householder's I - 2vv'/(v'v), with every reflector operation
        // stored as RN53(RN64(operation)). Use the shared native/software
        // arithmetic, including each dot-product term and vector update.
        v[0] = excel_numeric::excel_x87_sub(x0, alpha);
        let vv = excel_numeric::excel_x87_add(rest, excel_numeric::excel_x87_mul(v[0], v[0]));
        let reflect = |a: &mut [f64]| {
            let t = v
                .iter()
                .zip(a[j..].iter())
                .fold(0.0, |sum, (vi, ai)| {
                    excel_numeric::excel_x87_add(sum, excel_numeric::excel_x87_mul(*vi, *ai))
                });
            let f = excel_numeric::excel_x87_div(excel_numeric::excel_x87_mul(2.0, t), vv);
            for (ai, vi) in a[j..].iter_mut().zip(v.iter()) {
                *ai = excel_numeric::excel_x87_sub(*ai, excel_numeric::excel_x87_mul(f, *vi));
            }
        };
        for (_, column) in matrix.iter_mut().skip(j + 1) {
            reflect(column);
        }
        reflect(&mut target);
        diagonal[j] = alpha;
        j += 1;
    }
    let rank = matrix.len();
    let mut solution = vec![0.0; rank];
    for j in (0..rank).rev() {
        let d = ((j + 1)..rank).fold(0.0, |sum, k| sum + matrix[k].1[j] * solution[k]);
        solution[j] = if diagonal[j] == 0.0 { 0.0 } else { (target[j] - d) / diagonal[j] };
    }
    let mut coefficients = vec![0.0; columns.len()];
    for (j, (predictor, _)) in matrix.iter().take(rank).enumerate() {
        if let Some(index) = predictor {
            coefficients[*index] = solution[j];
        }
    }
    let intercept = if use_const {
        coefficients
            .iter()
            .zip(means.iter())
            .fold(y_mean, |a, (b, m)| a - b * m)
    } else {
        0.0
    };
    let r = (0..rank)
        .map(|j| {
            (0..rank)
                .map(|k| if k < j { 0.0 } else if k == j { diagonal[j] } else { matrix[k].1[j] })
                .collect()
        })
        .collect();
    Ok(RegressionFit {
        model: LinearRegressionModel {
            coefficients,
            intercept,
        },
        order: matrix.iter().map(|(predictor, _)| *predictor).collect(),
        r,
        qty: target,
        means,
        y: centred_y,
        observations: n,
    })
}

/// Equality of cross products of finite binary64 values, without rounding or
/// overflow. A significand product needs at most 106 bits, so u128 is sufficient.
fn exact_product(a: f64, b: f64) -> (bool, u128, i32) {
    let parts = |x: f64| {
        let bits = x.to_bits();
        let exponent = ((bits >> 52) & 0x7ff) as i32;
        let fraction = bits & ((1u64 << 52) - 1);
        if exponent == 0 {
            (fraction, -1074)
        } else {
            (fraction | (1u64 << 52), exponent - 1075)
        }
    };
    let (am, ae) = parts(a);
    let (bm, be) = parts(b);
    let product = u128::from(am) * u128::from(bm);
    if product == 0 {
        return (false, 0, 0);
    }
    let zeros = product.trailing_zeros();
    (a.is_sign_negative() != b.is_sign_negative(), product >> zeros, ae + be + zeros as i32)
}

fn exactly_proportional(a: &[f64], b: &[f64]) -> bool {
    if a.len() != b.len() || !a.iter().chain(b).all(|x| x.is_finite()) {
        return false;
    }
    // Initially zero columns have their own recorded degrees-of-freedom policy.
    let Some((&ax, &bx)) = a.iter().zip(b).find(|(x, y)| **x != 0.0 && **y != 0.0) else {
        return false;
    };
    a.iter().zip(b).all(|(&x, &y)| exact_product(x, bx) == exact_product(y, ax))
}

fn linear_fit_from_data(
    known_y: &NumericVector,
    known_x: &PreparedPredictorInput,
    use_const: bool,
) -> Result<RegressionFit, WorksheetErrorCode> {
    let observation_count = known_y.values.len();
    if known_x.predictors.rows != observation_count {
        return Err(WorksheetErrorCode::Ref);
    }
    let columns: Vec<Vec<f64>> = (0..known_x.predictors.cols)
        .map(|col| {
            (0..observation_count)
                .map(|row| known_x.predictors.cell(row, col))
                .collect()
        })
        .collect();
    least_squares(&columns, &known_y.values, use_const)
}

fn exponential_fit_from_data(
    known_y: &NumericVector,
    known_x: &PreparedPredictorInput,
    use_const: bool,
) -> Result<RegressionFit, WorksheetErrorCode> {
    let mut log_y = Vec::with_capacity(known_y.values.len());
    for value in &known_y.values {
        if *value <= 0.0 || !value.is_finite() {
            return Err(WorksheetErrorCode::Num);
        }
        log_y.push(excel_numeric::excel_log(*value));
    }
    let logged = NumericVector {
        values: log_y,
        shape: known_y.shape,
    };
    linear_fit_from_data(&logged, known_x, use_const)
}

fn eval_value_from_output_shape(
    values: &[f64],
    shape: OutputShape,
) -> Result<CalcValue, WorksheetErrorCode> {
    match shape {
        OutputShape::Scalar => Ok(CalcValue::number(values[0])),
        OutputShape::Row(cols) => CalcArray::from_rows(vec![
            values
                .iter()
                .copied()
                .take(cols)
                .map(CalcValue::number)
                .collect::<Vec<_>>(),
        ])
        .map(CalcValue::array)
        .ok_or(WorksheetErrorCode::Value),
        OutputShape::Column(rows) => CalcArray::from_rows(
            values
                .iter()
                .copied()
                .take(rows)
                .map(|value| vec![CalcValue::number(value)])
                .collect::<Vec<_>>(),
        )
        .map(CalcValue::array)
        .ok_or(WorksheetErrorCode::Value),
        OutputShape::Matrix { rows, cols } => {
            let row_major = values
                .chunks(cols)
                .take(rows)
                .map(|chunk| chunk.iter().copied().map(CalcValue::number).collect())
                .collect::<Vec<Vec<CalcValue>>>();
            CalcArray::from_rows(row_major)
                .map(CalcValue::array)
                .ok_or(WorksheetErrorCode::Value)
        }
    }
}

fn row_array(values: &[f64]) -> CalcValue {
    CalcValue::array(
        CalcArray::from_rows(vec![
            values.iter().copied().map(CalcValue::number).collect(),
        ])
        .expect("non-empty coefficient row"),
    )
}

/// The statistics of a fit to one observation with a constant: the first row `first` (the
/// coefficients, then the constant), every standard error 0, R squared 1 and its standard error
/// 0, F #NUM! with 0 degrees of freedom, and both sums of squares 0, the cells beyond the second
/// column of rows 3 to 5 #N/A (Excel 16.0.20430, 2026-10-25, dc-6au: `=LINEST(5,{1;2;3},TRUE,TRUE)`
/// spills 0,0,0,5 / 0,0,0,0 / 1,0,#N/A,#N/A / #NUM!,0,#N/A,#N/A / 0,0,#N/A,#N/A, `=LINEST(5,2,TRUE,TRUE)`
/// 0,5 / 0,0 / 1,0 / #NUM!,0 / 0,0, and LOGEST's the same below its first row).
fn one_observation_statistics(first: Vec<f64>) -> CalcValue {
    let cols = first.len();
    let pad = |head: Vec<CalcValue>| {
        let mut row = head;
        while row.len() < cols {
            row.push(CalcValue::error(WorksheetErrorCode::NA));
        }
        row
    };
    let n = CalcValue::number;
    CalcValue::array(
        CalcArray::from_rows(vec![
            first.into_iter().map(n).collect(),
            vec![n(0.0); cols],
            pad(vec![n(1.0), n(0.0)]),
            pad(vec![CalcValue::error(WorksheetErrorCode::Num), n(0.0)]),
            pad(vec![n(0.0), n(0.0)]),
        ])
        .expect("five rows of one length"),
    )
}

/// The statistics block of LINEST and LOGEST below the first row `first` (the coefficients in
/// reverse predictor order, then the constant), from the fit's QR factors: the standard errors
/// of the coefficients and of the constant (#N/A without one), R squared and the standard error
/// of y, F and the degrees of freedom, and the regression and residual sums of squares. LOGEST's
/// statistics are those of the fit to ln y.
fn statistics_block(fit: &RegressionFit, use_const: bool, first: Vec<f64>) -> CalcValue {
    if fit.order.is_empty() {
        return one_observation_statistics(first);
    }
    let n = fit.observations;
    let p = fit.order.len();
    let predictors = p - usize::from(use_const);
    let squares = |values: &[f64]| values.iter().fold(0.0, |sum, v| sum + v * v);
    let ss_resid = squares(&fit.qty[p..]);
    let ss_total = squares(&fit.y);
    let ss_reg = ss_total - ss_resid;
    let df = (n - p) as f64;
    // With no degrees of freedom left Excel gives R squared 1, every standard error 0 and F
    // #NUM! (an exact fit of as many columns as observations).
    let exact = n == p;
    let r2 = if exact { 1.0 } else { ss_reg / ss_total };
    let se_y = if exact { 0.0 } else { (ss_resid / df).sqrt() };
    let f = if exact { f64::NAN } else { (ss_reg / predictors as f64) / (ss_resid / df) };
    // The inverse of the triangular factor over the columns with a nonzero diagonal (a column
    // kept as zero comes last and has standard error 0), column by column.
    let nz = (0..p).take_while(|j| fit.r[*j][*j] != 0.0).count();
    let mut inverse = vec![vec![0.0; nz]; nz];
    for j in 0..nz {
        inverse[j][j] = 1.0 / fit.r[j][j];
        for i in (0..j).rev() {
            let sum = ((i + 1)..=j).fold(0.0, |sum, k| sum + fit.r[i][k] * inverse[k][j]);
            inverse[i][j] = -sum / fit.r[i][i];
        }
    }
    // (R'R)^-1 = R^-1 R^-T.
    let covariance =
        |a: usize, b: usize| (0..nz).fold(0.0, |sum, k| sum + inverse[a][k] * inverse[b][k]);
    let mut se = vec![0.0; fit.model.coefficients.len()];
    for (position, predictor) in fit.order.iter().enumerate().take(nz) {
        if let Some(index) = predictor {
            if !exact {
                se[*index] = (ss_resid * covariance(position, position) / df).sqrt();
            }
        }
    }
    let n_a = || CalcValue::error(WorksheetErrorCode::NA);
    let number = |v: f64| {
        if v.is_finite() { CalcValue::number(v) } else { CalcValue::error(WorksheetErrorCode::Num) }
    };
    let cols = first.len();
    let mut second: Vec<CalcValue> = se.iter().rev().map(|v| number(*v)).collect();
    if use_const {
        // The constant is my - sum(b_j m_j): its variance is w'(R'R)^-1 w with w = 1 at the
        // ones column and -m_j at predictor j's column.
        let w: Vec<f64> = fit
            .order
            .iter()
            .take(nz)
            .map(|predictor| match predictor {
                None => 1.0,
                Some(index) => -fit.means[*index],
            })
            .collect();
        let mut variance = 0.0;
        for a in 0..nz {
            for b in 0..nz {
                variance += w[a] * w[b] * covariance(a, b);
            }
        }
        let se_b = if exact { 0.0 } else { (ss_resid / df * variance).sqrt() };
        second.push(number(se_b));
    } else {
        second.push(n_a());
    }
    let pad = |head: Vec<CalcValue>| {
        let mut row = head;
        while row.len() < cols {
            row.push(n_a());
        }
        row
    };
    CalcValue::array(
        CalcArray::from_rows(vec![
            first.into_iter().map(CalcValue::number).collect(),
            second,
            pad(vec![number(r2), number(se_y)]),
            pad(vec![number(f), number(df)]),
            pad(vec![number(ss_reg), number(ss_resid)]),
        ])
        .expect("five rows of one length"),
    )
}