מרחבים ותתי־מרחבים

function rref(input, tolerance = 1e-10, scaleColumnCount) {
  if (!Array.isArray(input) || input.length === 0 || !Array.isArray(input[0]) || input[0].length === 0) {
    throw new TypeError("matrix must be a nonempty rectangular array");
  }
  const columnCount = input[0].length;
  if (!input.every((row) => Array.isArray(row) && row.length === columnCount && row.every(Number.isFinite))) {
    throw new TypeError("matrix entries must be finite and rows must have equal length");
  }
  if (!Number.isFinite(tolerance) || tolerance <= 0) {
    throw new RangeError("tolerance must be positive and finite");
  }
  const referenceColumns = scaleColumnCount ?? columnCount;
  if (!Number.isInteger(referenceColumns) || referenceColumns < 1 || referenceColumns > columnCount) {
    throw new RangeError("scaleColumnCount must select existing columns");
  }
  const rowScales = input.map((row) =>
    row.slice(0, referenceColumns)
      .reduce((largest, value) => Math.max(largest, Math.abs(value)), 0) ||
      row.reduce((largest, value) => Math.max(largest, Math.abs(value)), 0) || 1,
  );
  const R = input.map((row, index) => row.map((value) => value / rowScales[index]));
  const rows = R.length, columns = R[0].length;
  const pivotColumns = [];
  let pivotRow = 0;

  for (let column = 0; column < columns && pivotRow < rows; column += 1) {
    let bestRow = pivotRow;
    for (let row = pivotRow + 1; row < rows; row += 1) {
      if (Math.abs(R[row][column]) > Math.abs(R[bestRow][column])) bestRow = row;
    }
    if (Math.abs(R[bestRow][column]) <= tolerance) continue;
    [R[pivotRow], R[bestRow]] = [R[bestRow], R[pivotRow]];
    const pivot = R[pivotRow][column];
    R[pivotRow] = R[pivotRow].map((value) => value / pivot);
    for (let row = 0; row < rows; row += 1) {
      if (row === pivotRow) continue;
      const factor = R[row][column];
      R[row] = R[row].map((value, index) => value - factor * R[pivotRow][index]);
    }
    pivotColumns.push(column);
    pivotRow += 1;
  }

  for (const row of R) {
    for (let column = 0; column < columns; column += 1) {
      if (Math.abs(row[column]) <= tolerance) {
        row[column] = 0;
      } else {
        const nearestInteger = Math.round(row[column]);
        if (Math.abs(row[column] - nearestInteger) <= tolerance * Math.max(1, Math.abs(row[column]))) {
          row[column] = nearestInteger;
        }
      }
    }
  }
  return { matrix: R, pivotColumns, numericalRank: pivotColumns.length, tolerance };
}

function analyzeMatrix(A, tolerance = 1e-10) {
  const reduced = rref(A, tolerance);
  const columnCount = A[0].length;
  const freeColumns = Array.from({ length: columnCount }, (_, index) => index)
    .filter((index) => !reduced.pivotColumns.includes(index));
  const nullspaceBasis = freeColumns.map((freeColumn) => {
    const vector = Array(columnCount).fill(0);
    vector[freeColumn] = 1;
    reduced.pivotColumns.forEach((pivotColumn, row) => {
      vector[pivotColumn] = -reduced.matrix[row][freeColumn];
    });
    return vector;
  });
  const columnBasis = reduced.pivotColumns.map((column) => A.map((row) => row[column]));
  return { ...reduced, freeColumns, columnBasis, nullspaceBasis };
}

function particularSolution(A, b, tolerance = 1e-10) {
  if (!Array.isArray(A) || !Array.isArray(b) || b.length !== A.length || !b.every(Number.isFinite)) {
    throw new TypeError("b must contain one finite value per row of A");
  }
  const columnCount = A[0].length;
  const reduced = rref(
    A.map((row, index) => [...row, b[index]]),
    tolerance,
    columnCount,
  );
  if (reduced.pivotColumns.includes(columnCount)) return null;
  const solution = Array(columnCount).fill(0);
  reduced.pivotColumns.forEach((column, row) => {
    if (column < columnCount) solution[column] = reduced.matrix[row][columnCount];
  });
  return solution;
}

const A = [[2, 4, 6, 4], [2, 5, 7, 6], [2, 3, 5, 2]];
const b = [4, 3, 5];
console.log(analyzeMatrix(A));
console.log({ particular: particularSolution(A, b) });

from math import isfinite

def rref(values, tolerance=1e-10, scale_column_count=None):
    if (not isinstance(values, list) or not values or
            not isinstance(values[0], list) or not values[0]):
        raise TypeError("matrix must be a nonempty rectangular list")
    column_count = len(values[0])
    if (any(not isinstance(row, list) or len(row) != column_count or
            any(not isinstance(value, (int, float)) or not isfinite(value) for value in row)
            for row in values)):
        raise TypeError("matrix entries must be finite and rows must have equal length")
    if not isinstance(tolerance, (int, float)) or not isfinite(tolerance) or tolerance <= 0:
        raise ValueError("tolerance must be positive and finite")
    reference_columns = column_count if scale_column_count is None else scale_column_count
    if (not isinstance(reference_columns, int) or
            reference_columns < 1 or reference_columns > column_count):
        raise ValueError("scale_column_count must select existing columns")
    row_scales = [
        max((abs(value) for value in row[:reference_columns]), default=0)
        or max((abs(value) for value in row), default=0)
        or 1
        for row in values
    ]
    R = [[value / row_scales[index] for value in row]
         for index, row in enumerate(values)]
    rows, columns = len(R), len(R[0])
    pivot_columns = []
    pivot_row = 0

    for column in range(columns):
        if pivot_row == rows:
            break
        best_row = max(range(pivot_row, rows), key=lambda row: abs(R[row][column]))
        if abs(R[best_row][column]) <= tolerance:
            continue
        R[pivot_row], R[best_row] = R[best_row], R[pivot_row]
        pivot = R[pivot_row][column]
        R[pivot_row] = [value / pivot for value in R[pivot_row]]
        for row in range(rows):
            if row == pivot_row:
                continue
            factor = R[row][column]
            R[row] = [value - factor * pivot_value
                      for value, pivot_value in zip(R[row], R[pivot_row])]
        pivot_columns.append(column)
        pivot_row += 1

    def clean(value):
        if abs(value) <= tolerance:
            return 0
        nearest_integer = round(value)
        if abs(value - nearest_integer) <= tolerance * max(1, abs(value)):
            return nearest_integer
        return value

    R = [[clean(value) for value in row] for row in R]
    return {"matrix": R, "pivot_columns": pivot_columns,
            "numerical_rank": len(pivot_columns), "tolerance": tolerance}

def analyze_matrix(A, tolerance=1e-10):
    reduced = rref(A, tolerance)
    column_count = len(A[0])
    free_columns = [column for column in range(column_count)
                    if column not in reduced["pivot_columns"]]
    nullspace_basis = []
    for free_column in free_columns:
        vector = [0] * column_count
        vector[free_column] = 1
        for row, pivot_column in enumerate(reduced["pivot_columns"]):
            vector[pivot_column] = -reduced["matrix"][row][free_column]
        nullspace_basis.append(vector)
    column_basis = [[row[column] for row in A]
                    for column in reduced["pivot_columns"]]
    return {**reduced, "free_columns": free_columns,
            "column_basis": column_basis, "nullspace_basis": nullspace_basis}

def particular_solution(A, b, tolerance=1e-10):
    if (not isinstance(A, list) or not isinstance(b, list) or len(b) != len(A) or
            any(not isinstance(value, (int, float)) or not isfinite(value) for value in b)):
        raise TypeError("b must contain one finite value per row of A")
    column_count = len(A[0])
    reduced = rref(
        [row + [b[index]] for index, row in enumerate(A)],
        tolerance,
        column_count,
    )
    if column_count in reduced["pivot_columns"]:
        return None
    solution = [0] * column_count
    for row, column in enumerate(reduced["pivot_columns"]):
        if column < column_count:
            solution[column] = reduced["matrix"][row][column_count]
    return solution

A = [[2, 4, 6, 4], [2, 5, 7, 6], [2, 3, 5, 2]]
b = [4, 3, 5]
print(analyze_matrix(A))
print({"particular": particular_solution(A, b)})