הפתרון השלם של Ax=b

function rref(input, relativeTolerance=1e-12, scaleColumnCount) {
  const matrix = input.map((row) => row.map(Number));
  const columns = matrix[0].length;
  const referenceColumns = scaleColumnCount ?? columns;
  const rowScales = matrix.map((row) =>
    Math.max(0, ...row.slice(0, referenceColumns).map(Math.abs)) ||
    Math.max(0, ...row.map(Math.abs)),
  );
  const pivots = [];
  let pivotRow = 0;

  for (let column = 0;
       column < columns && pivotRow < matrix.length;
       column += 1) {
    let best = -1;
    let bestRatio = 0;
    for (let row = pivotRow; row < matrix.length; row += 1) {
      const ratio = rowScales[row] === 0
        ? 0
        : Math.abs(matrix[row][column])/rowScales[row];
      if (ratio > bestRatio) {
        bestRatio = ratio;
        best = row;
      }
    }
    if (best === -1 || bestRatio <= relativeTolerance) continue;
    [matrix[pivotRow], matrix[best]] = [matrix[best], matrix[pivotRow]];
    [rowScales[pivotRow], rowScales[best]] = [rowScales[best], rowScales[pivotRow]];

    const scale = matrix[pivotRow][column];
    matrix[pivotRow] = matrix[pivotRow].map((value) => value/scale);
    for (let row = 0; row < matrix.length; row += 1) {
      if (row === pivotRow) continue;
      const factor = matrix[row][column];
      matrix[row] = matrix[row].map(
        (value, j) => value - factor*matrix[pivotRow][j],
      );
    }
    pivots.push(column);
    pivotRow += 1;
  }

  return {
    matrix: matrix.map((row) =>
      row.map((value) => Math.abs(value) <= relativeTolerance ? 0 : value),
    ),
    pivots,
  };
}

function completeSolution(A, b, relativeTolerance=1e-12) {
  const variables = A[0].length;
  const augmented = A.map((row, i) => [...row, b[i]]);
  const { matrix, pivots } = rref(
    augmented,
    relativeTolerance,
    variables,
  );
  const coefficientPivots = pivots.filter((column) => column < variables);
  const inconsistent = pivots.includes(variables);

  if (inconsistent) {
    return {
      consistent: false,
      rankA: coefficientPivots.length,
      rankAugmented: pivots.length,
      rref: matrix,
    };
  }

  const free = Array.from({ length: variables }, (_, j) => j)
    .filter((j) => !coefficientPivots.includes(j));
  const particular = Array(variables).fill(0);
  coefficientPivots.forEach((column, row) => {
    particular[column] = matrix[row][variables];
  });
  const nullspaceBasis = free.map((freeColumn) => {
    const vector = Array(variables).fill(0);
    vector[freeColumn] = 1;
    coefficientPivots.forEach((column, row) => {
      vector[column] = -matrix[row][freeColumn];
    });
    return vector;
  });

  return {
    consistent: true,
    rankA: coefficientPivots.length,
    rankAugmented: coefficientPivots.length,
    particular,
    nullspaceBasis,
    rref: matrix,
  };
}

const A = [[2, 4, 6, 4], [2, 5, 7, 6], [2, 3, 5, 2]];
console.log(completeSolution(A, [4, 3, 5]));
console.log(completeSolution(A, [4, 3, 4]));

def rref(source, relative_tolerance=1e-12, scale_column_count=None):
    matrix = [[float(value) for value in row] for row in source]
    columns = len(matrix[0])
    reference_columns = columns if scale_column_count is None else scale_column_count
    row_scales = [
        max((abs(value) for value in row[:reference_columns]), default=0)
        or max((abs(value) for value in row), default=0)
        for row in matrix
    ]
    pivots = []
    pivot_row = 0

    for column in range(columns):
        if pivot_row == len(matrix):
            break
        best = max(
            range(pivot_row, len(matrix)),
            key=lambda row: (
                abs(matrix[row][column])/row_scales[row]
                if row_scales[row] else 0.0
            ),
        )
        best_ratio = (
            abs(matrix[best][column])/row_scales[best]
            if row_scales[best] else 0.0
        )
        if best_ratio <= relative_tolerance:
            continue
        matrix[pivot_row], matrix[best] = matrix[best], matrix[pivot_row]
        row_scales[pivot_row], row_scales[best] = (
            row_scales[best], row_scales[pivot_row]
        )
        scale = matrix[pivot_row][column]
        matrix[pivot_row] = [value/scale for value in matrix[pivot_row]]
        for row in range(len(matrix)):
            if row == pivot_row:
                continue
            factor = matrix[row][column]
            matrix[row] = [
                value - factor*pivot_value
                for value, pivot_value
                in zip(matrix[row], matrix[pivot_row])
            ]
        pivots.append(column)
        pivot_row += 1

    matrix = [
        [0.0 if abs(value) <= relative_tolerance else value for value in row]
        for row in matrix
    ]
    return matrix, pivots

def complete_solution(A, b, relative_tolerance=1e-12):
    variables = len(A[0])
    augmented = [list(row) + [b[i]] for i, row in enumerate(A)]
    matrix, pivots = rref(
        augmented,
        relative_tolerance,
        variables,
    )
    coefficient_pivots = [
        column for column in pivots if column < variables
    ]
    inconsistent = variables in pivots

    if inconsistent:
        return {
            "consistent": False,
            "rank_A": len(coefficient_pivots),
            "rank_augmented": len(pivots),
            "rref": matrix,
        }

    free = [
        column for column in range(variables)
        if column not in coefficient_pivots
    ]
    particular = [0.0]*variables
    for row, column in enumerate(coefficient_pivots):
        particular[column] = matrix[row][variables]

    nullspace_basis = []
    for free_column in free:
        vector = [0.0]*variables
        vector[free_column] = 1.0
        for row, column in enumerate(coefficient_pivots):
            vector[column] = -matrix[row][free_column]
        nullspace_basis.append(vector)

    return {
        "consistent": True,
        "rank_A": len(coefficient_pivots),
        "rank_augmented": len(coefficient_pivots),
        "particular": particular,
        "nullspace_basis": nullspace_basis,
        "rref": matrix,
    }

A = [[2, 4, 6, 4], [2, 5, 7, 6], [2, 3, 5, 2]]
print(complete_solution(A, [4, 3, 5]))
print(complete_solution(A, [4, 3, 4]))