Cramer's rule

From Rosetta Code
Task
Cramer's rule
You are encouraged to solve this task according to the task description, using any language you may know.
In linear algebra, Cramer's rule is an explicit formula for the solution of a system of linear equations with as many equations as unknowns, valid whenever the system has a unique solution. It expresses the solution in terms of the determinants of the (square) coefficient matrix and of matrices obtained from it by replacing one column by the vector of right hand sides of the equations.


Given

which in matrix format is

Then the values of and can be found as follows:


Task

Given the following system of equations:

solve for , , and , using Cramer's rule.

11l

Translation of: Nim
F det(mm)
   V m = copy(mm)
   V result = 1.0

   L(j) 0 .< m.len
      V imax = j
      L(i) j + 1 .< m.len
         I m[i][j] > m[imax][j]
            imax = i

      I imax != j
         swap(&m[imax], &m[j])
         result = -result

      I abs(m[j][j]) < 1e-12
         R Float.infinity

      L(i) j + 1 .< m.len
         V mult = -m[i][j] / m[j][j]
         L(k) 0 .< m.len
            m[i][k] += mult * m[j][k]

   L(i) 0 .< m.len
      result *= m[i][i]
   R result

F cramerSolve(aa, detA, b, col)
   V a = copy(aa)
   L(i) 0 .< a.len
      a[i][col] = b[i]
   R det(a) / detA

V A = [[2.0, -1.0,  5.0,  1.0],
       [3.0,  2.0,  2.0, -6.0],
       [1.0,  3.0,  3.0, -1.0],
       [5.0, -2.0, -3.0,  3.0]]

V B = [-3.0, -32.0, -47.0, 49.0]

V detA = det(A)

L(i) 0 .< A.len
   print(‘#3.3’.format(cramerSolve(A, detA, B, i)))
Output:
  2.000
-12.000
 -4.000
  1.000

Ada

with Ada.Text_IO;
with Ada.Numerics.Generic_Real_Arrays;

procedure Cramers_Rules is

   type Real is new Float;
   --  This is the type we want to use in the matrix and vector

   package Real_Arrays is
      new Ada.Numerics.Generic_Real_Arrays (Real);

   use Real_Arrays;

   function Solve_Cramer (M : in Real_Matrix;
                          V : in Real_Vector)
                         return Real_Vector
   is
      Denominator : Real;
      Nom_Matrix  : Real_Matrix (M'Range (1),
                                 M'Range (2));
      Numerator   : Real;
      Result      : Real_Vector (M'Range (1));
   begin
      if
        M'Length (2) /= V'Length or
        M'Length (1) /= M'Length (2)
      then
         raise Constraint_Error with "Dimensions does not match";
      end if;

      Denominator := Determinant (M);

      for Col in V'Range loop
         Nom_Matrix := M;

         --  Substitute column
         for Row in V'Range loop
            Nom_Matrix (Row, Col) := V (Row);
         end loop;

         Numerator    := Determinant (Nom_Matrix);
         Result (Col) := Numerator / Denominator;
      end loop;

      return Result;
   end Solve_Cramer;

   procedure Put (V : Real_Vector) is
      use Ada.Text_IO;
      package Real_IO is
         new Ada.Text_IO.Float_IO (Real);
   begin
      Put ("[");
      for E of V loop
         Real_IO.Put (E, Exp => 0, Aft => 2);
         Put (" ");
      end loop;
      Put ("]");
      New_Line;
   end Put;

   M : constant Real_Matrix := ((2.0, -1.0,  5.0,  1.0),
                                (3.0,  2.0,  2.0, -6.0),
                                (1.0,  3.0,  3.0, -1.0),
                                (5.0, -2.0, -3.0,  3.0));
   V : constant Real_Vector := (-3.0, -32.0, -47.0, 49.0);
   R : constant Real_Vector := Solve_Cramer (M, V);
begin
   Put (R);
end Cramers_Rules;
Output:
[ 2.00 -12.00 -4.00  1.00 ]

ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Uses the non-standard DET operator available in Algol 68G.

# returns the solution of a.x = b via Cramer's rule                    #
#         this is for REAL arrays, could define additional operators   #
#         for INT, COMPL, etc.                                         #
PRIO CRAMER = 1;
OP   CRAMER = ( [,]REAL a, []REAL b )[]REAL:
     IF 1 UPB a /= 2 UPB a
     OR 1 LWB a /= 2 LWB a
     OR 1 UPB a /=   UPB b
     THEN
        # the array sizes and bounds do not match                       #
        print( ( "Invaid parameters to CRAMER", newline ) );
        stop
     ELIF REAL deta = DET a;
          det a = 0
     THEN
        # a is singular                                                 #
        print( ( "Singular matrix for CRAMER", newline ) );
        stop
     ELSE
        # the arrays have matching bounds                               #
        [ LWB b : UPB b ]REAL result;
        FOR col FROM LWB b TO UPB b DO
            # form a matrix from a with the col'th column replaced by b #
            [ 1 LWB a : 1 UPB a, 2 LWB a : 2 UPB a ]REAL m := a;
            m[ : , col ] := b[ : AT 1 ];
            # col'th result elemet as per Cramer's rule                 #
            result[ col ] := DET m / det a
        OD;
        result
     FI; # CRAMER #

# test CRAMER using the matrix and column vector specified in the task  #
[,]REAL a = ( (  2, -1,  5,  1 )
            , (  3,  2,  2, -6 )
            , (  1,  3,  3, -1 )
            , (  5, -2, -3,  3 )
            );
[]REAL  b = (  -3
            , -32
            , -47
            ,  49
            );
[]REAL  solution = a CRAMER b;
FOR c FROM LWB solution TO UPB solution DO
    print( ( " ", fixed( solution[ c ], -8, 4 ) ) )
OD;
print( ( newline ) )
Output:
   2.0000 -12.0000  -4.0000   1.0000

C

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

typedef struct {
    int n;
    double **elems;
} SquareMatrix;

SquareMatrix init_square_matrix(int n, double elems[n][n]) {
    SquareMatrix A = {
        .n = n,
        .elems = malloc(n * sizeof(double *))
    };
    for(int i = 0; i < n; ++i) {
        A.elems[i] = malloc(n * sizeof(double));
        for(int j = 0; j < n; ++j)
            A.elems[i][j] = elems[i][j];
    }

    return A;
}

SquareMatrix copy_square_matrix(SquareMatrix src) {
    SquareMatrix dest;
    dest.n = src.n;
    dest.elems = malloc(dest.n * sizeof(double *));
    for(int i = 0; i < dest.n; ++i) {
        dest.elems[i] = malloc(dest.n * sizeof(double));
        for(int j = 0; j < dest.n; ++j)
            dest.elems[i][j] = src.elems[i][j];
    }

    return dest;
}

double det(SquareMatrix A) {
    double det = 1;

    for(int j = 0; j < A.n; ++j) {
        int i_max = j;
        for(int i = j; i < A.n; ++i)
            if(A.elems[i][j] > A.elems[i_max][j])
                i_max = i;

        if(i_max != j) {
            for(int k = 0; k < A.n; ++k) {
                double tmp = A.elems[i_max][k];
                A.elems[i_max][k] = A.elems[j][k];
                A.elems[j][k]     = tmp;
            }

            det *= -1;
        }

        if(abs(A.elems[j][j]) < 1e-12) {
            puts("Singular matrix!");
            return NAN;
        }

        for(int i = j + 1; i < A.n; ++i) {
            double mult = -A.elems[i][j] / A.elems[j][j];
            for(int k = 0; k < A.n; ++k)
                A.elems[i][k] += mult * A.elems[j][k];
        }
    }

    for(int i = 0; i < A.n; ++i)
        det *= A.elems[i][i];

    return det;
}

void deinit_square_matrix(SquareMatrix A) {
    for(int i = 0; i < A.n; ++i)
        free(A.elems[i]);
    free(A.elems);
}

double cramer_solve(SquareMatrix A, double det_A, double *b, int var) {
    SquareMatrix tmp = copy_square_matrix(A);
    for(int i = 0; i < tmp.n; ++i)
        tmp.elems[i][var] = b[i];

    double det_tmp = det(tmp);
    deinit_square_matrix(tmp);

    return det_tmp / det_A;
}

int main(int argc, char **argv) {
#define N 4
    double elems[N][N] = {
        { 2, -1,  5,  1},
        { 3,  2,  2, -6},
        { 1,  3,  3, -1},
        { 5, -2, -3,  3}
    };
    SquareMatrix A = init_square_matrix(N, elems);

    SquareMatrix tmp = copy_square_matrix(A);
    int det_A = det(tmp);
    deinit_square_matrix(tmp);

    double b[] = {-3, -32, -47, 49};

    for(int i = 0; i < N; ++i)
        printf("%7.3lf\n", cramer_solve(A, det_A, b, i));

    deinit_square_matrix(A);
    return EXIT_SUCCESS;
}
Output:
  2.000
-12.000
 -4.000
  1.000

C#

Instead of copying a bunch of arrays, I made a class with an indexer that simply accesses the correct items in the original matrix.

using System;
using System.Collections.Generic;
using static System.Linq.Enumerable;

public static class CramersRule
{
    public static void Main() {
        var equations = new [] {
            new [] { 2, -1,  5,  1,  -3 },
            new [] { 3,  2,  2, -6, -32 },
            new [] { 1,  3,  3, -1, -47 },
            new [] { 5, -2, -3,  3,  49 }
        };
        var solution = SolveCramer(equations);
        Console.WriteLine(solution.DelimitWith(", "));
    }

    public static int[] SolveCramer(int[][] equations) {
        int size = equations.Length;
        if (equations.Any(eq => eq.Length != size + 1)) throw new ArgumentException($"Each equation must have {size+1} terms.");
        int[,] matrix = new int[size, size];
        int[] column = new int[size];
        for (int r = 0; r < size; r++) {
            column[r] = equations[r][size];
            for (int c = 0; c < size; c++) {
                matrix[r, c] = equations[r][c];
            }
        }
        return Solve(new SubMatrix(matrix, column));
    }

    private static int[] Solve(SubMatrix matrix) {
        int det = matrix.Det();
        if (det == 0) throw new ArgumentException("The determinant is zero.");

        int[] answer = new int[matrix.Size];
        for (int i = 0; i < matrix.Size; i++) {
            matrix.ColumnIndex = i;
            answer[i] = matrix.Det() / det;
        }
        return answer;
    }

    //Extension method from library.
    static string DelimitWith<T>(this IEnumerable<T> source, string separator = " ") =>
        string.Join(separator ?? " ", source ?? Empty<T>());

    private class SubMatrix
    {
        private int[,] source;
        private SubMatrix prev;
        private int[] replaceColumn;

        public SubMatrix(int[,] source, int[] replaceColumn) {
            this.source = source;
            this.replaceColumn = replaceColumn;
            this.prev = null;
            this.ColumnIndex = -1;
            Size = replaceColumn.Length;
        }

        private SubMatrix(SubMatrix prev, int deletedColumnIndex = -1) {
            this.source = null;
            this.prev = prev;
            this.ColumnIndex = deletedColumnIndex;
            Size = prev.Size - 1;
        }

        public int ColumnIndex { get; set; }
        public int Size { get; }

        public int this[int row, int column] {
            get {
                if (source != null) return column == ColumnIndex ? replaceColumn[row] : source[row, column];
                return prev[row + 1, column < ColumnIndex ? column : column + 1];
            }
        }

        public int Det() {
            if (Size == 1) return this[0, 0];
            if (Size == 2) return this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0];
            SubMatrix m = new SubMatrix(this);
            int det = 0;
            int sign = 1;
            for (int c = 0; c < Size; c++) {
                m.ColumnIndex = c;
                int d = m.Det();
                det += this[0, c] * d * sign;
                sign = -sign;
            }
            return det;
        }

        public void Print() {
            for (int r = 0; r < Size; r++) {
                Console.WriteLine(Range(0, Size).Select(c => this[r, c]).DelimitWith(", "));
            }
            Console.WriteLine();
        }
    }

}
Output:
2, -12, -4, 1

C++

Translation of: C#
#include <algorithm>
#include <iostream>
#include <vector>

class SubMatrix {
    const std::vector<std::vector<double>> *source;
    std::vector<double> replaceColumn;
    const SubMatrix *prev;
    size_t sz;
    int colIndex = -1;

public:
    SubMatrix(const std::vector<std::vector<double>> &src, const std::vector<double> &rc) : source(&src), replaceColumn(rc), prev(nullptr), colIndex(-1) {
        sz = replaceColumn.size();
    }

    SubMatrix(const SubMatrix &p) : source(nullptr), prev(&p), colIndex(-1) {
        sz = p.size() - 1;
    }

    SubMatrix(const SubMatrix &p, int deletedColumnIndex) : source(nullptr), prev(&p), colIndex(deletedColumnIndex) {
        sz = p.size() - 1;
    }

    int columnIndex() const {
        return colIndex;
    }
    void columnIndex(int index) {
        colIndex = index;
    }

    size_t size() const {
        return sz;
    }

    double index(int row, int col) const {
        if (source != nullptr) {
            if (col == colIndex) {
                return replaceColumn[row];
            } else {
                return (*source)[row][col];
            }
        } else {
            if (col < colIndex) {
                return prev->index(row + 1, col);
            } else {
                return prev->index(row + 1, col + 1);
            }
        }
    }

    double det() const {
        if (sz == 1) {
            return index(0, 0);
        }
        if (sz == 2) {
            return index(0, 0) * index(1, 1) - index(0, 1) * index(1, 0);
        }
        SubMatrix m(*this);
        double det = 0.0;
        int sign = 1;
        for (size_t c = 0; c < sz; ++c) {
            m.columnIndex(c);
            double d = m.det();
            det += index(0, c) * d * sign;
            sign = -sign;
        }
        return det;
    }
};

std::vector<double> solve(SubMatrix &matrix) {
    double det = matrix.det();
    if (det == 0.0) {
        throw std::runtime_error("The determinant is zero.");
    }

    std::vector<double> answer(matrix.size());
    for (int i = 0; i < matrix.size(); ++i) {
        matrix.columnIndex(i);
        answer[i] = matrix.det() / det;
    }
    return answer;
}

std::vector<double> solveCramer(const std::vector<std::vector<double>> &equations) {
    int size = equations.size();
    if (std::any_of(
        equations.cbegin(), equations.cend(),
        [size](const std::vector<double> &a) { return a.size() != size + 1; }
    )) {
        throw std::runtime_error("Each equation must have the expected size.");
    }

    std::vector<std::vector<double>> matrix(size);
    std::vector<double> column(size);
    for (int r = 0; r < size; ++r) {
        column[r] = equations[r][size];
        matrix[r].resize(size);
        for (int c = 0; c < size; ++c) {
            matrix[r][c] = equations[r][c];
        }
    }

    SubMatrix sm(matrix, column);
    return solve(sm);
}

template<typename T>
std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
    auto it = v.cbegin();
    auto end = v.cend();

    os << '[';
    if (it != end) {
        os << *it++;
    }
    while (it != end) {
        os << ", " << *it++;
    }

    return os << ']';
}

int main() {
    std::vector<std::vector<double>> equations = {
        { 2, -1,  5,  1,  -3},
        { 3,  2,  2, -6, -32},
        { 1,  3,  3, -1, -47},
        { 5, -2, -3,  3,  49},
    };

    auto solution = solveCramer(equations);
    std::cout << solution << '\n';

    return 0;
}
Output:
[2, -12, -4, 1]

Common Lisp

(defun minor (m col)
  (loop with dim = (1- (array-dimension m 0))
        with result = (make-array (list dim dim))
        for i from 1 to dim
        for r = (1- i)
        do (loop with c = 0
                 for j to dim
                 when (/= j col)
                   do (setf (aref result r c) (aref m i j))
                      (incf c))
        finally (return result)))

(defun det (m)
  (assert (= (array-rank m) 2))
  (assert (= (array-dimension m 0) (array-dimension m 1)))
  (let ((dim (array-dimension m 0)))
    (if (= dim 1)
        (aref m 0 0)
        (loop for col below dim
              for sign = 1 then (- sign)
              sum (* sign (aref m 0 col) (det (minor m col)))))))

(defun replace-column (m col values)
  (let* ((dim (array-dimension m 0))
         (result (make-array (list dim dim))))
    (dotimes (r dim result)
      (dotimes (c dim)
        (setf (aref result r c)
              (if (= c col) (aref values r) (aref m r c)))))))

(defun solve (m v)
  (loop with dim = (array-dimension m 0)
        with det = (det m)
        for col below dim
        collect (/ (det (replace-column m col v)) det)))

(solve #2A((2 -1  5  1)
           (3  2  2 -6)
           (1  3  3 -1)
           (5 -2 -3  3))
       #(-3 -32 -47 49))
Output:
(2 -12 -4 1)

D

Translation of: Kotlin
import std.array : array, uninitializedArray;
import std.range : iota;
import std.stdio : writeln;
import std.typecons : tuple;

alias vector = double[4];
alias matrix = vector[4];

auto johnsonTrotter(int n) {
    auto p = iota(n).array;
    auto q = iota(n).array;
    auto d = uninitializedArray!(int[])(n);
    d[] = -1;
    auto sign = 1;
    int[][] perms;
    int[] signs;

    void permute(int k) {
        if (k >= n) {
            perms ~= p.dup;
            signs ~= sign;
            sign *= -1;
            return;
        }
        permute(k + 1);
        foreach (i; 0..k) {
            auto z = p[q[k] + d[k]];
            p[q[k]] = z;
            p[q[k] + d[k]] = k;
            q[z] = q[k];
            q[k] += d[k];
            permute(k + 1);
        }
        d[k] *= -1;
    }

    permute(0);
    return tuple!("sigmas", "signs")(perms, signs);
}

auto determinant(matrix m) {
    auto jt = johnsonTrotter(m.length);
    auto sum = 0.0;
    foreach (i,sigma; jt.sigmas) {
        auto prod = 1.0;
        foreach (j,s; sigma) {
            prod *= m[j][s];
        }
        sum += jt.signs[i] * prod;
    }
    return sum;
}

auto cramer(matrix m, vector d) {
    auto divisor = determinant(m);
    auto numerators = uninitializedArray!(matrix[])(m.length);
    foreach (i; 0..m.length) {
        foreach (j; 0..m.length) {
            foreach (k; 0..m.length) {
                numerators[i][j][k] = m[j][k];
            }
        }
    }
    vector v;
    foreach (i; 0..m.length) {
        foreach (j; 0..m.length) {
            numerators[i][j][i] = d[j];
        }
    }
    foreach (i; 0..m.length) {
        v[i] = determinant(numerators[i]) / divisor;
    }
    return v;
}

void main() {
    matrix m = [
        [2.0, -1.0,  5.0,  1.0],
        [3.0,  2.0,  2.0, -6.0],
        [1.0,  3.0,  3.0, -1.0],
        [5.0, -2.0, -3.0,  3.0]
    ];
    vector d = [-3.0, -32.0, -47.0, 49.0];
    auto wxyz = cramer(m, d);
    writeln("w = ", wxyz[0], ", x = ", wxyz[1], ", y = ", wxyz[2], ", z = ", wxyz[3]);
}
Output:
w = 2, x = -12, y = -4, z = 1

EchoLisp

(lib 'matrix)
(string-delimiter "")
(define (cramer A B (X)) ;; --> vector X
	(define  (determinant A))
	(for/vector [(i (matrix-col-num A))]
		(set! X (matrix-set-col! (array-copy A) i B))
		(// (determinant X) )))
		
(define (task)
	(define A (list->array 
  	'( 2 -1 5 1 3 2 2 -6 1 3 3 -1 5 -2 -3 3) 4 4))
	(define B #(-3 -32 -47 49))
	(writeln "Solving A * X = B")
	(array-print A)
	(writeln "B = " B)
	(writeln "X = " (cramer A B)))
Output:
(task)
Solving A * X = B    
  2   -1   5    1  
  3   2    2    -6 
  1   3    3    -1 
  5   -2   -3   3  
B =      #( -3 -32 -47 49)    
X =      #( 2 -12 -4 1)    

Factor

USING: kernel math math.matrices.laplace prettyprint sequences ;
IN: rosetta-code.cramers-rule

: replace-col ( elt n seq -- seq' ) flip [ set-nth ] keep flip ;

: solve ( m v -- seq )
    dup length <iota> [
        rot [ replace-col ] keep [ determinant ] bi@ /
    ] 2with map ;

: cramers-rule-demo ( -- )
    {
        { 2 -1  5  1 }
        { 3  2  2 -6 }
        { 1  3  3 -1 }
        { 5 -2 -3  3 }
    }
    { -3 -32 -47 49 } solve . ;

MAIN: cramers-rule-demo
Output:
{ 2 -12 -4 1 }

Fortran

In Numerical Methods That Work (Usually), in the section What not to compute, F. S. Acton remarks "...perhaps we should be glad he didn't resort to Cramer's rule (still taught as the practical method in some high schools) and solve his equations as the ratios of determinants - a process that requires labor proportional to N! if done in the schoolboy manner. The contrast with N3 can be startling!" And further on, "Having hinted darkly at my computational fundamentalism, it is probably time to commit to a public heresy by denouncing recursive calculations. I have never seen a numerical problem arising from the physical world that was best calculated by a recursive subroutine..."

Since this problem requires use of Cramer's rule, one might as well be hung for a sheep instead of a lamb, so the traditions of Old Fortran and heavy computation will be ignored and the fearsome RECURSIVE specification employed so that the determinants will be calculated recursively, all the way down to N = 1 even though the N = 2 case is easy. This requires F90 and later. Similarly, the MODULE protocol will be employed, even though there is no significant context to share. The alternative method for calculating a determinant involves generating permutations, a tiresome process.

Array passing via the modern arrangements of F90 is a source of novel difficulty to set against the slight convenience of not having to pass an additional parameter, N. Explicitly, at least. There are "secret" additional parameters when an array is being passed in the modern way, which are referred to by the new SIZE function. Anyway, for an order N square matrix, the array must be declared A(N,N), and specifically not something like A(100,100) with usage only of elements up to N = 7, say, because the locations in storage of elements in use would be quite different from those used by an array declared A(7,7). This means that the array must be re-declared for each different size usage, a tiresome and error-inviting task. One-dimensional arrays do not have this problem, but they do have to be "long enough" so B and X might as well be included. This also means that the auxiliary matrices needed within the routines have to be made the right size, and fortunately they can be declared in a way that requests this without the blather of ALLOCATE, this being a protocol introduced by Algol in the 1960s. Unfortunately, there is no scheme such as in pl/i to declare AUX "like" A, so some grotesquery results. And in the case of function DET which needs an array of order N - 1, when its recursion bottoms out with N = 1 it will have declared MINOR(0,0), a rather odd situation that fortunately evokes no complaint, and a test run in which its "value" was written out by WRITE (6,*) MINOR produced a blank line: no complaint there either, presumably because zero elements were being sent forth and so there was no improper access of ... nothing.

With matrices, there is a problem all the way from the start in 1958. Everyone agrees that a matrix should be indexed as A(row,column) and that when written out, rows should run down the page and columns across. This is unexceptional and the F90 function MATMUL follows this usage. However, Fortran has always stored its array elements in what is called "column major" order, which is to say that element A(1,1) is followed by element A(2,1) in storage, not A(1,2). Thus, if an array is written (or read) by something like WRITE (6,*) A, consecutive elements, written along a line, will be A(1,1), A(2,1), A(3,1), ... So, subroutine SHOWMATRIX is employed to write the matrix out in the desired form, and to read the values into the array, an explicit loop is used to place them where expected rather than just READ(INF,*) A

Similarly, if instead a DATA statement were used to initialise the array for the example problem, and it looked something like
      DATA A/2, -1,  5,  1
     1       3,  2,  2, -6
     2       1,  3,  3, -1
     3       5, -2, -3,  3/
(ignoring integer/floating-point issues) thus corresponding to the layout of the example problem, there would need to be a statement A = TRANSPOSE(A) to obtain the required order. I have never seen an explanation of why this choice was made for Fortran.
      MODULE BAD IDEA	!Employ Cramer's rule for solving A.x = b...
       INTEGER MSG	!Might as well be here.
       CONTAINS		!The details.
        SUBROUTINE SHOWMATRIX(A)	!With nice vertical bars.
         DOUBLE PRECISION A(:,:)	!The matrix.
         INTEGER R,N			!Assistants.
          N = SIZE(A, DIM = 1)		!Instead of passing an explicit parameter.
          DO R = 1,N			!Work down the rows.
            WRITE (MSG,1) A(R,:)		!Fling forth a row at a go.
    1       FORMAT ("|",<N>F12.3,"|")		!Bounded by bars.
          END DO			!On to the next row.
        END SUBROUTINE SHOWMATRIX	!Furrytran's default order is the transpose.

        RECURSIVE DOUBLE PRECISION FUNCTION DET(A)	!Determine the determinant.
         DOUBLE PRECISION A(:,:)	!The square matrix, order N.
         DOUBLE PRECISION MINOR(SIZE(A,DIM=1) - 1,SIZE(A,DIM=1) - 1)	!Order N - 1.
         DOUBLE PRECISION D		!A waystation.
         INTEGER C,N			!Steppers.
          N = SIZE(A, DIM = 1)		!Suplied via secret parameters.
          IF (N .LE. 0) THEN		!This really ought not happen.
            STOP "DET: null array!"		!But I'm endlessly suspicious.
          ELSE IF (N .NE. SIZE(A, DIM = 2)) THEN	!And I'd rather have a decent message
            STOP "DET: not a square array!"	!In place of a crashed run.
          ELSE IF (N .EQ. 1) THEN	!Alright, now get on with it.
            DET = A(1,1)		!This is really easy.
          ELSE				!But otherwise,
            D = 0			!Here we go.
            DO C = 1,N			!Step along the columns of the first row.
              CALL FILLMINOR(C)			!Produce the auxiliary array for each column.
              IF (MOD(C,2) .EQ. 0) THEN		!Odd or even case?
                D = D - A(1,C)*DET(MINOR)	!Even: subtract.
               ELSE				!Otherwise,
                D = D + A(1,C)*DET(MINOR)	!Odd: add.
              END IF				!So much for that term.
            END DO			!On to the next.
            DET = D		!Declare the result.
          END IF	!So much for that.
         CONTAINS	!An assistant.
           SUBROUTINE FILLMINOR(CC)	!Corresponding to A(1,CC).
           INTEGER CC	!The column being omitted.
           INTEGER R	!A stepper.
           DO R = 2,N		!Ignoring the first row,
             MINOR(R - 1,1:CC - 1) = A(R,1:CC - 1)	!Copy columns 1 to CC - 1. There may be none.
             MINOR(R - 1,CC:) = A(R,CC + 1:)		!And from CC + 1 to N. There may be none.
           END DO		!On to the next row.
          END SUBROUTINE FILLMINOR	!Divide and conquer.
        END FUNCTION DET	!Rather than mess with permutations.

        SUBROUTINE CRAMER(A,X,B)	!Solve A.x = b, where A is a matrix...
Careful! The array must be A(N,N), and not say A(100,100) of which only up to N = 6 are in use.
         DOUBLE PRECISION A(:,:)	!A square matrix. I hope.
         DOUBLE PRECISION X(:),B(:)	!Column matrices look rather like 1-D arrays.
         DOUBLE PRECISION AUX(SIZE(A,DIM=1),SIZE(A,DIM=1))	!Can't say "LIKE A", as in pl/i, alas.
         DOUBLE PRECISION D	!To be calculated once.
         INTEGER N		!The order of the square matrix. I hope.
         INTEGER C		!A stepper.
          N = SIZE(A, DIM = 1)	!Alright, what's the order of battle?
          D = DET(A)		!Once only.
          IF (D.EQ.0) STOP "Cramer: zero determinant!"	!Surely, this won't happen...
          AUX = A		!Prepare the assistant.
          DO C = 1,N		!Step across the columns.
            IF (C.GT.1) AUX(1:N,C - 1) = A(1:N,C - 1)	!Repair previous damage.
            AUX(1:N,C) = B(1:N)		!Place the current damage.
            X(C) = DET(AUX)/D		!The result!
          END DO		!On to the next column.
        END SUBROUTINE CRAMER	!This looks really easy!
      END MODULE BAD IDEA	!But actually, it is a bad idea for N > 2.

      PROGRAM TEST	!Try it and see.
      USE BAD IDEA	!Just so.
      DOUBLE PRECISION, ALLOCATABLE ::A(:,:), B(:), X(:)	!Custom work areas.
      INTEGER N,R	!Assistants..
      INTEGER INF	!An I/O unit.

      MSG = 6		!Output.
      INF = 10		!For reading test data.
      OPEN (INF,FILE="Test.dat",STATUS="OLD",ACTION="READ")	!As in this file..

Chew into the next problem.
   10 IF (ALLOCATED(A)) DEALLOCATE(A)	!First,
      IF (ALLOCATED(B)) DEALLOCATE(B)	!Get rid of
      IF (ALLOCATED(X)) DEALLOCATE(X)	!The hired help.
      READ (INF,*,END = 100) N		!Since there is a new order.
      IF (N.LE.0) GO TO 100		!Perhaps a final order.
      WRITE (MSG,11) N			!Othewise, announce prior to acting.
   11 FORMAT ("Order ",I0," matrix A, as follows...")	!In case something goes wrong.
      ALLOCATE(A(N,N))			!For instance,
      ALLOCATE(B(N))			!Out of memory.
      ALLOCATE(X(N))			!But otherwise, a tailored fit.
      DO R = 1,N			!Now read in the values for the matrix.
        READ(INF,*,END=667,ERR=665) A(R,:),B(R)	!One row of A at a go, followed by B's value.
      END DO				!In free format.
      CALL SHOWMATRIX(A)		!Show what we have managed to obtain.
      WRITE (MSG,12) "In the scheme A.x = b, b = ",B	!In case something goes wrong.
   12 FORMAT (A,<N>F12.6)		!How many would be too many?
      CALL CRAMER(A,X,B)		!The deed!
      WRITE (MSG,12) "    Via Cramer's rule, x = ",X	!The result!
      GO TO 10				!And try for another test problem.

Complaints.
  665 WRITE (MSG,666) "Format error",R	!I know where I came from.
  666 FORMAT (A," while reading row ",I0,"!")	!So I can refer to R.
      GO TO 100		!So much for that.
  667 WRITE (MSG,666) "End-of-file", R		!Some hint as to where.

Closedown.
  100 WRITE (6,*) "That was interesting."	!Quite.
      END	!Open files are closed, allocated memory is released.

Oddly, the Compaq Visual Fortran F90/95 compiler is confused by the usage "BAD IDEA" instead of "BADIDEA" - spaces are not normally relevant in Fortran source. Anyway, file Test.dat has been filled with the numbers of the example, as follows:

4           /The order, for A.x = b.
2  -1   5   1,  -3    /First row of A, b
3   2   2  -6, -32    /Second row...
1   3   3  -1, -47     third row.
5  -2  -3   3,  49    /Last row.

Fortran's free-form allows a comma, a tab, and spaces between numbers, and regards the / as starting a comment, but, because each row is read separately, once the required five (N + 1) values have been read, no further scanning of the line takes place and the next READ statement will start with a new line of input. So the / isn't needed for the third row, as shown. Omitted values lead to confusion as the input process would read additional lines to fill the required count and everything gets out of step. Echoing input very soon after it is obtained is helpful in making sense of such mistakes.

For more practical use it would probably be better to constrain the freedom somewhat, perhaps requiring that all the N + 1 values for a row appear on one input record. In such a case, the record could first be read into a text variable (from which the data would be read) so that if a problem arises the text could be printed as a part of the error message. But, this requires guessing a suitably large length for the text variable so as to accommodate the longest possible input line.

Output:

Order 4 matrix A, as follows...
|       2.000      -1.000       5.000       1.000|
|       3.000       2.000       2.000      -6.000|
|       1.000       3.000       3.000      -1.000|
|       5.000      -2.000      -3.000       3.000|
In the scheme A.x = b, b =    -3.000000  -32.000000  -47.000000   49.000000
    Via Cramer's rule, x =     2.000000  -12.000000   -4.000000    1.000000
 That was interesting.

And at this point I suddenly noticed that the habits of Old Fortran are not so easily suppressed: all calculations are done with double precision. Curiously enough, for the specific example data, the same results are obtained if all variables are integer.


FreeBASIC

Function determinant(matrix() As Double) As Double
      Dim As long n=Ubound(matrix,1),sign=1
      Dim As Double det=1,s=1
      Dim As Double b(1 To n,1 To n)
      For c As long=1 To n
            For d As long=1 To n
                  b(c,d)=matrix(c,d)
            Next d
      Next c
      #macro pivot(num)
      For p1 As long  = num To n - 1
            For p2 As long  = p1 + 1 To n  
                  If Abs(b(p1,num))<Abs(b(p2,num)) Then
                        sign=-sign
                        For g As long=1 To n
                              Swap b(p1,g),b(p2,g)
                        Next g
                  End If
            Next p2
      Next p1
      #endmacro
      For k As long=1 To n-1
            pivot(k)
            For row As long =k To n-1
                  If b(row+1,k)=0 Then Exit For
                  Var f=b(k,k)/b(row+1,k)
                  s=s*f
                  For g As long=1 To n
                        b((row+1),g)=(b((row+1),g)*f-b(k,g))/f
                  Next g
            Next row
      Next k
      For z As long=1 To n
            det=det*b(z,z)
      Next z
      Return sign*det
End Function

'CRAMER COLUMN SWAPS
Sub swapcolumn(m() As Double,c() As Double,_new() As Double,column As long)
      Redim _new(1 To Ubound(m,1),1 To Ubound(m,1))
      For x As long=1 To Ubound(m,1)
            For y As long=1 To Ubound(m,1)
                  _new(x,y)=m(x,y)
            Next y
      Next x
      For z As long=1 To Ubound(m,1)
            _new(z,column)=c(z)
      Next z
End Sub

Sub solve(mat() As Double,rhs() As Double,_out() As Double)
      redim _out(Lbound(mat,1) To Ubound(mat,1))
      Redim As Double result(Lbound(mat,1) To Ubound(mat,1),Lbound(mat,1) To Ubound(mat,1))
      Dim As Double maindeterminant=determinant(mat())
      If Abs(maindeterminant) < 1e-12 Then Print "singular":Exit Sub 
      For column As Long=1 To Ubound(mat,1)
            swapcolumn(mat(),rhs(),result(),column)
            _out(column)= determinant(result())/maindeterminant
      Next
End Sub



Dim As Double MainMat(1 To ...,1 To ...)={{2,-1,5,1}, _
                                          {3,2,2,-6}, _
                                          {1,3,3,-1}, _
                                          {5,-2,-3,3}}

Dim As Double rhs(1 To ...)={-3,-32,-47,49}

Redim ans() As Double
solve(MainMat(),rhs(),ans())

For n As Long=1 To Ubound(ans)
      Print Csng(ans(n))
Next
Sleep
Output:
 2
-12
-4
 1

Go

Library gonum:

package main

import (
    "fmt"

    "gonum.org/v1/gonum/mat"
)

var m = mat.NewDense(4, 4, []float64{
    2, -1, 5, 1,
    3, 2, 2, -6,
    1, 3, 3, -1,
    5, -2, -3, 3,
})

var v = []float64{-3, -32, -47, 49}

func main() {
    x := make([]float64, len(v))
    b := make([]float64, len(v))
    d := mat.Det(m)
    for c := range v {
        mat.Col(b, c, m)
        m.SetCol(c, v)
        x[c] = mat.Det(m) / d
        m.SetCol(c, b)
    }
    fmt.Println(x)
}
Output:
[2 -12.000000000000007 -4.000000000000001 1.0000000000000009]

Library go.matrix:

package main

import (
    "fmt"

    "github.com/skelterjohn/go.matrix"
)

var m = matrix.MakeDenseMatrixStacked([][]float64{
    {2, -1, 5, 1},
    {3, 2, 2, -6},
    {1, 3, 3, -1},
    {5, -2, -3, 3},
})

var v = []float64{-3, -32, -47, 49}

func main() {
    x := make([]float64, len(v))
    b := make([]float64, len(v))
    d := m.Det()
    for c := range v {
        m.BufferCol(c, b)
        m.FillCol(c, v)
        x[c] = m.Det() / d
        m.FillCol(c, b)
    }
    fmt.Println(x)
}
Output:
[2.0000000000000004 -11.999999999999998 -4 0.9999999999999999]

Groovy

Translation of: Java
class CramersRule {
    static void main(String[] args) {
        Matrix mat = new Matrix(Arrays.asList(2d, -1d, 5d, 1d),
                Arrays.asList(3d, 2d, 2d, -6d),
                Arrays.asList(1d, 3d, 3d, -1d),
                Arrays.asList(5d, -2d, -3d, 3d))
        List<Double> b = Arrays.asList(-3d, -32d, -47d, 49d)
        println("Solution = " + cramersRule(mat, b))
    }

    private static List<Double> cramersRule(Matrix matrix, List<Double> b) {
        double denominator = matrix.determinant()
        List<Double> result = new ArrayList<>()
        for (int i = 0; i < b.size(); i++) {
            result.add(matrix.replaceColumn(b, i).determinant() / denominator)
        }
        return result
    }

    private static class Matrix {
        private List<List<Double>> matrix

        @Override
        String toString() {
            return matrix.toString()
        }

        @SafeVarargs
        Matrix(List<Double>... lists) {
            matrix = new ArrayList<>()
            for (List<Double> list : lists) {
                matrix.add(list)
            }
        }

        Matrix(List<List<Double>> mat) {
            matrix = mat
        }

        double determinant() {
            if (matrix.size() == 1) {
                return get(0, 0)
            }
            if (matrix.size() == 2) {
                return get(0, 0) * get(1, 1) - get(0, 1) * get(1, 0)
            }
            double sum = 0
            double sign = 1
            for (int i = 0; i < matrix.size(); i++) {
                sum += sign * get(0, i) * coFactor(0, i).determinant()
                sign *= -1
            }
            return sum
        }

        private Matrix coFactor(int row, int col) {
            List<List<Double>> mat = new ArrayList<>()
            for (int i = 0; i < matrix.size(); i++) {
                if (i == row) {
                    continue
                }
                List<Double> list = new ArrayList<>()
                for (int j = 0; j < matrix.size(); j++) {
                    if (j == col) {
                        continue
                    }
                    list.add(get(i, j))
                }
                mat.add(list)
            }
            return new Matrix(mat)
        }

        private Matrix replaceColumn(List<Double> b, int column) {
            List<List<Double>> mat = new ArrayList<>()
            for (int row = 0; row < matrix.size(); row++) {
                List<Double> list = new ArrayList<>()
                for (int col = 0; col < matrix.size(); col++) {
                    double value = get(row, col)
                    if (col == column) {
                        value = b.get(row)
                    }
                    list.add(value)
                }
                mat.add(list)
            }
            return new Matrix(mat)
        }

        private double get(int row, int col) {
            return matrix.get(row).get(col)
        }
    }
}
Output:
Solution = [2.0, -12.0, -4.0, 1.0]

Haskell

Version 1

import Data.Matrix

solveCramer :: (Ord a, Fractional a) => Matrix a -> Matrix a -> Maybe [a]
solveCramer a y
  | da == 0 = Nothing
  | otherwise = Just $ map (\i -> d i / da) [1..n]
  where da = detLU a
        d i = detLU $ submatrix 1 n 1 n $ switchCols i (n+1) ay
        ay = a <|> y
        n = ncols a

task = solveCramer a y
  where a = fromLists [[2,-1, 5, 1]
                      ,[3, 2, 2,-6]
                      ,[1, 3, 3,-1]
                      ,[5,-2,-3, 3]]
        y = fromLists [[-3], [-32], [-47], [49]]
Output:
λ> task
Just [2.0000000000000004,-11.999999999999998,-4.0,0.9999999999999999]

Version 2

We use Rational numbers for having more precision. a % b is the rational a / b.

s_permutations :: [a] -> [([a], Int)]
s_permutations = flip zip (cycle [1, -1]) . (foldl aux [[]])
  where aux items x = do
          (f,item) <- zip (cycle [reverse,id]) items
          f (insertEv x item)
        insertEv x [] = [[x]]
        insertEv x l@(y:ys) = (x:l) :  map (y:) (insertEv x ys)
 
mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs). zipWith (\vs u -> map (u*) vs) vss) uss

matI::(Num a) => Int -> [[a]]
matI n = [ [fromIntegral.fromEnum $ i == j | i <- [1..n]] | j <- [1..n]]
 
elemPos::[[a]] -> Int -> Int -> a
elemPos ms i j = (ms !! i) !! j
 
prod:: Num a => ([[a]] -> Int -> Int -> a) -> [[a]] -> [Int] -> a
prod f ms = product.zipWith (f ms) [0..]
 
s_determinant:: Num a => ([[a]] -> Int -> Int -> a) -> [[a]] -> [([Int],Int)] -> a
s_determinant f ms = sum.map (\(is,s) -> fromIntegral s * prod f ms is)
 
elemCramerPos::Int -> Int -> [[a]] -> [[a]] -> Int -> Int -> a
elemCramerPos l k ks ms i j = if j /= l then elemPos ms i j else elemPos ks i k
 
solveCramer:: [[Rational]] -> [[Rational]] -> [[Rational]]
solveCramer ms ks = xs
  where
  xs | d /= 0    = go (reverse [0..pred.length.head $ ks])
     | otherwise = []
  go (u:us) = foldl glue (col u) us
  glue us u = zipWith (\ys (y:_) -> y:ys) us (col u)
  col k = map (\l -> [(/d) $ s_determinant (elemCramerPos l k ks) ms ps]) $ ls
  ls = [0..pred.length $ ms]
  ps = s_permutations ls
  d = s_determinant elemPos ms ps
 
task::[[Rational]] -> [[Rational]] -> IO()
task a b = do
  let x         = solveCramer a b
  let u         = map (map fromRational) x
  let y         = mult a x
  let identity  = matI (length x)
  let a1        = solveCramer a identity
  let h         = mult a a1
  let z         = mult a1 b
  putStrLn "a ="
  mapM_ print a
  putStrLn "b ="
  mapM_ print b
  putStrLn "solve: a * x = b => x = solveCramer a b ="
  mapM_ print x
  putStrLn "u = fromRationaltoDouble x ="
  mapM_ print u
  putStrLn "verification: y = a * x = mult a x ="
  mapM_ print y
  putStrLn $ "test: y == b = "
  print $ y == b
  putStrLn "identity matrix: identity ="
  mapM_ print identity
  putStrLn "find: a1 = inv(a) => solve: a * a1 = identity => a1 = solveCramer a identity ="
  mapM_ print a1
  putStrLn "verification: h = a * a1 = mult a a1 ="
  mapM_ print h
  putStrLn $ "test: h == identity = "
  print $ h == identity
  putStrLn "z = a1 * b = mult a1 b ="
  mapM_ print z
  putStrLn "test: z == x ="
  print $ z == x
  
main = do
  let a  = [[2,-1, 5, 1]
           ,[3, 2, 2,-6]
           ,[1, 3, 3,-1]
           ,[5,-2,-3, 3]]
  let b   =  [[-3], [-32], [-47], [49]]
  task a b
Output:
a =
[2 % 1,(-1) % 1,5 % 1,1 % 1]
[3 % 1,2 % 1,2 % 1,(-6) % 1]
[1 % 1,3 % 1,3 % 1,(-1) % 1]
[5 % 1,(-2) % 1,(-3) % 1,3 % 1]
b =
[(-3) % 1]
[(-32) % 1]
[(-47) % 1]
[49 % 1]
solve: a * x = b => x = solveCramer a b =
[2 % 1]
[(-12) % 1]
[(-4) % 1]
[1 % 1]
u = fromRationaltoDouble x =
[2.0]
[-12.0]
[-4.0]
[1.0]
verification: y = a * x = mult a x =
[(-3) % 1]
[(-32) % 1]
[(-47) % 1]
[49 % 1]
test: y == b = 
True
identity matrix: identity =
[1 % 1,0 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
find: a1 = inv(a) => solve: a * a1 = identity => a1 = solveCramer a identity =
[4 % 171,11 % 171,10 % 171,8 % 57]
[(-55) % 342,(-23) % 342,119 % 342,2 % 57]
[107 % 684,(-5) % 684,11 % 684,(-7) % 114]
[7 % 684,(-109) % 684,103 % 684,7 % 114]
verification: h = a * a1 = mult a a1 =
[1 % 1,0 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
test: h == identity = 
True
z = a1 * b = mult a1 b =
[2 % 1]
[(-12) % 1]
[(-4) % 1]
[1 % 1]
test: z == x =
True

Version 3

import Data.List

determinant::(Fractional a, Ord a) => [[a]] -> a
determinant ls = if null ls then 0 else pivot 1 (zip ls [(0::Int)..])
  where
  good rs ts = (abs.head.fst $ ts) <= (abs.head.fst $ rs)
  go us (vs,i) = if v == 0 then (ws,i) else (zipWith (\x y -> y - x*v) us ws,i)
    where (v,ws) = (head $ vs,tail vs)
  change i (ys:zs) = map (\xs -> if (==i).snd $ xs then ys else xs) zs
  pivot d [] = d
  pivot d zs@((_,j):ys) = if 0 == u then 0 else pivot e ws
    where
    e = if i == j then u*d else -u*d
    ((u:us),i) = foldl1 (\rs ts ->  if good rs ts then rs else ts) zs
    ws = map (go (map (/u) us)) $ if i == j then ys else change i zs

solveCramer::(Fractional a, Ord a) => [[a]] -> [[a]] -> [[a]]
solveCramer as bs = if 0 == d then [] else ans bs
  where 
  d = determinant as
  ans = transpose.map go.transpose
    where
    ms = zip [0..] (transpose as)
    go us =  [ (/d) $ determinant [if i /= j then vs else us | (j,vs) <- ms] | (i,_) <- ms]
    
matI::(Num a) => Int -> [[a]]
matI n = [ [fromIntegral.fromEnum $ i == j | i <- [1..n]] | j <- [1..n]]

mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs). zipWith (\vs u -> map (u*) vs) vss) uss

task::[[Rational]] -> [[Rational]] -> IO()
task a b = do
  let x         = solveCramer a b
  let u         = map (map fromRational) x
  let y         = mult a x
  let identity  = matI (length x)
  let a1        = solveCramer a identity
  let h         = mult a a1
  let z         = mult a1 b
  putStrLn "a ="
  mapM_ print a
  putStrLn "b ="
  mapM_ print b
  putStrLn "solve: a * x = b => x = solveCramer a b ="
  mapM_ print x
  putStrLn "u = fromRationaltoDouble x ="
  mapM_ print u
  putStrLn "verification: y = a * x = mult a x ="
  mapM_ print y
  putStrLn $ "test: y == b = "
  print $ y == b
  putStrLn "identity matrix: identity ="
  mapM_ print identity
  putStrLn "find: a1 = inv(a) => solve: a * a1 = identity => a1 = solveCramer a identity ="
  mapM_ print a1
  putStrLn "verification: h = a * a1 = mult a a1 ="
  mapM_ print h
  putStrLn $ "test: h == identity = "
  print $ h == identity
  putStrLn "z = a1 * b = mult a1 b ="
  mapM_ print z
  putStrLn "test: z == x ="
  print $ z == x
 
main = do
  let a  = [[2,-1, 5, 1]
           ,[3, 2, 2,-6]
           ,[1, 3, 3,-1]
           ,[5,-2,-3, 3]]
  let b  = [[-3], [-32], [-47], [49]]
  task a b
Output:
a =
[2 % 1,(-1) % 1,5 % 1,1 % 1]
[3 % 1,2 % 1,2 % 1,(-6) % 1]
[1 % 1,3 % 1,3 % 1,(-1) % 1]
[5 % 1,(-2) % 1,(-3) % 1,3 % 1]
b =
[(-3) % 1]
[(-32) % 1]
[(-47) % 1]
[49 % 1]
solve: a * x = b => x = solveCramer a b =
[2 % 1]
[(-12) % 1]
[(-4) % 1]
[1 % 1]
u = fromRationaltoDouble x =
[2.0]
[-12.0]
[-4.0]
[1.0]
verification: y = a * x = mult a x =
[(-3) % 1]
[(-32) % 1]
[(-47) % 1]
[49 % 1]
test: y == b = 
True
identity matrix: identity =
[1 % 1,0 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
find: a1 = inv(a) => solve: a * a1 = identity => a1 = solveCramer a identity =
[4 % 171,11 % 171,10 % 171,8 % 57]
[(-55) % 342,(-23) % 342,119 % 342,2 % 57]
[107 % 684,(-5) % 684,11 % 684,(-7) % 114]
[7 % 684,(-109) % 684,103 % 684,7 % 114]
verification: h = a * a1 = mult a a1 =
[1 % 1,0 % 1,0 % 1,0 % 1]
[0 % 1,1 % 1,0 % 1,0 % 1]
[0 % 1,0 % 1,1 % 1,0 % 1]
[0 % 1,0 % 1,0 % 1,1 % 1]
test: h == identity = 
True
z = a1 * b = mult a1 b =
[2 % 1]
[(-12) % 1]
[(-4) % 1]
[1 % 1]
test: z == x =
True

J

Implementation:

cramer=:4 :0
  A=. x [ b=. y
  det=. -/ .*
  A %~&det (i.#A) b"_`[`]}&.|:"0 2 A
)

Task data:

A=: _&".;._2]t=: 0 :0
  2 -1  5  1
  3  2  2 -6
  1  3  3 -1
  5 -2 -3  3
)

b=: _3 _32 _47 49

Task example:

   A cramer b
2 _12 _4 1

Java

Supports double data type. A more robust solution would support arbitrary precision integers, arbitrary precision decimals, arbitrary precision rationals, or even arbitrary precision algebraic numbers.

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class CramersRule {

    public static void main(String[] args) {
        Matrix mat = new Matrix(Arrays.asList(2d, -1d, 5d, 1d), 
                                Arrays.asList(3d, 2d, 2d, -6d), 
                                Arrays.asList(1d, 3d, 3d, -1d),
                                Arrays.asList(5d, -2d, -3d, 3d));
        List<Double> b = Arrays.asList(-3d, -32d, -47d, 49d);
        System.out.println("Solution = " + cramersRule(mat, b));
    }
    
    private static List<Double> cramersRule(Matrix matrix, List<Double> b) {
        double denominator = matrix.determinant();
        List<Double> result = new ArrayList<>();
        for ( int i = 0 ; i < b.size() ; i++ ) {
            result.add(matrix.replaceColumn(b, i).determinant() / denominator);
        }
        return result;
    }
        
    private static class Matrix {
        
        private List<List<Double>> matrix;
        
        @Override
        public String toString() {
            return matrix.toString();
        }
        
        @SafeVarargs
        public Matrix(List<Double> ... lists) {
            matrix = new ArrayList<>();
            for ( List<Double> list : lists) {
                matrix.add(list);
            }
        }
        
        public Matrix(List<List<Double>> mat) {
            matrix = mat;
        }
        
        public double determinant() {
            if ( matrix.size() == 1 ) {
                return get(0, 0);
            }
            if ( matrix.size() == 2 ) {
                return get(0, 0) * get(1, 1) - get(0, 1) * get(1, 0);
            }
            double sum = 0;
            double sign = 1;
            for ( int i = 0 ; i < matrix.size() ; i++ ) {
                sum += sign * get(0, i) * coFactor(0, i).determinant();
                sign *= -1;
            }
            return sum;
        }
        
        private Matrix coFactor(int row, int col) {
            List<List<Double>> mat = new ArrayList<>();
            for ( int i = 0 ; i < matrix.size() ; i++ ) {
                if ( i == row ) {
                    continue;
                }
                List<Double> list = new ArrayList<>();
                for ( int j = 0 ; j < matrix.size() ; j++ ) {
                    if ( j == col ) {
                        continue;
                    }
                    list.add(get(i, j));
                }
                mat.add(list);
            }
            return new Matrix(mat);
        }

        private Matrix replaceColumn(List<Double> b, int column) {
            List<List<Double>> mat = new ArrayList<>();
            for ( int row = 0 ; row < matrix.size() ; row++ ) {
                List<Double> list = new ArrayList<>();
                for ( int col = 0 ; col < matrix.size() ; col++ ) {
                    double value = get(row, col);
                    if ( col == column ) {
                        value = b.get(row);
                    }
                    list.add(value);
                }
                mat.add(list);
            }
            return new Matrix(mat);
        }

        private double get(int row, int col) {
            return matrix.get(row).get(col);
        }
        
    }

}
Output:
Solution = [2.0, -12.0, -4.0, 1.0]

JavaScript

var matrix = [
	[2, -1,  5,  1],
	[3,  2,  2, -6],
	[1,  3,  3, -1],
	[5, -2, -3,  3]
];
var freeTerms = [-3, -32, -47, 49];

var result = cramersRule(matrix,freeTerms);
console.log(result);

/**
 * Compute Cramer's Rule
 * @param  {array} matrix    x,y,z, etc. terms
 * @param  {array} freeTerms
 * @return {array}           solution for x,y,z, etc.
 */
function cramersRule(matrix,freeTerms) {
	var det = detr(matrix),
		returnArray = [],
		i,
		tmpMatrix;

	for(i=0; i < matrix[0].length; i++) {
		var tmpMatrix = insertInTerms(matrix, freeTerms,i)
		returnArray.push(detr(tmpMatrix)/det)
	}
	return returnArray;
}

/**
 * Inserts single dimensional array into
 * @param  {array} matrix multidimensional array to have ins inserted into
 * @param  {array} ins single dimensional array to be inserted vertically into matrix
 * @param  {array} at  zero based offset for ins to be inserted into matrix
 * @return {array}     New multidimensional array with ins replacing the at column in matrix
 */
function insertInTerms(matrix, ins, at) {
	var tmpMatrix = clone(matrix),
		i;
	for(i=0; i < matrix.length; i++) {
		tmpMatrix[i][at] = ins[i];
	}
	return tmpMatrix;
}
/**
 * Compute the determinate of a matrix.  No protection, assumes square matrix
 * function borrowed, and adapted from MIT Licensed numericjs library (www.numericjs.com)
 * @param  {array} m Input Matrix (multidimensional array)
 * @return {number}   result rounded to 2 decimal
 */
function detr(m) {
	var ret = 1,
		k,
		A=clone(m),
		n=m[0].length,
		alpha;

	for(var j =0; j < n-1; j++) {
		k=j;
		for(i=j+1;i<n;i++) { if(Math.abs(A[i][j]) > Math.abs(A[k][j])) { k = i; } }
		if(k !== j) {
		    temp = A[k]; A[k] = A[j]; A[j] = temp;
		    ret *= -1;
		}
		Aj = A[j];
		for(i=j+1;i<n;i++) {
			Ai = A[i];
            alpha = Ai[j]/Aj[j];
            for(k=j+1;k<n-1;k+=2) {
                k1 = k+1;
                Ai[k] -= Aj[k]*alpha;
                Ai[k1] -= Aj[k1]*alpha;
            }
            if(k!==n) { Ai[k] -= Aj[k]*alpha; }
        }
        if(Aj[j] === 0) { return 0; }
        ret *= Aj[j];
	    }
    return Math.round(ret*A[j][j]*100)/100;
}

/**
 * Clone two dimensional Array using ECMAScript 5 map function and EcmaScript 3 slice
 * @param  {array} m Input matrix (multidimensional array) to clone
 * @return {array}   New matrix copy
 */
function clone(m) {
	return m.map(function(a){return a.slice();});
}
Output:
[ 2, -12, -4, 1 ]

jq

Adapted from Wren (*)

Works with: jq

Works with gojq, the Go implementation of jq

(*) Note that cramer(a;d) as defined here assumes that d is a 1-d vector in accordance with the task description and common practice.

# The minor of the input matrix after removing the specified row and column.
# Assumptions: the input is square and the indices are hunky dory.
def minor(rowNum; colNum):
  . as $in
  | (length - 1) as $len
  | reduce range(0; $len) as $i (null;
      reduce range(0; $len) as $j (.;
        if $i < rowNum and $j < colNum
        then .[$i][$j] = $in[$i][$j]
        elif $i >= rowNum and $j < colNum
        then .[$i][$j] = $in[$i+1][$j]
        elif $i < rowNum and $j >= colNum
        then .[$i][$j] = $in[$i][$j+1]
        else .[$i][$j] = $in[$i+1][$j+1]
        end) );

# The determinant using Laplace expansion.
# Assumption: the matrix is square
def det:
 . as $a
 | length as $nr
 | if $nr == 1 then .[0][0]
   elif $nr == 2 then .[1][1] * .[0][0] - .[0][1] * .[1][0]
   else reduce range(0; $nr) as $i (
     { sign: 1, sum: 0 };
     ($a|minor(0; $i)) as $m
     | .sum += .sign * $a[0][$i] * ($m|det)
     | .sign *= -1 )
     | .sum
   end ;

# Solve A X = D using Cramer's method
# a is assumed to be a JSON array representing the 2-d square matrix A
# d is assumed to be a JSON array representing the 1-d vector D
def cramer(a; d):
  (a | length) as $n
  | (a | det) as $ad
  | if $ad == 0 then "matrix determinant is 0" | error
    else reduce range(0; $n) as $c (null;
      (reduce range(0; $n) as $r (a; .[$r][$c] = d[$r])) as $aa
      | .[$c] = ($aa|det) / $ad )
    end ;
 
def a: [
    [2, -1,  5,  1],
    [3,  2,  2, -6],
    [1,  3,  3, -1],
    [5, -2, -3,  3]
];
 
def d:
  [ -3, -32, -47, 49 ] ;

"Solution is \(cramer(a; d))"
Output:
Solution is [2,-12,-4,1]

Julia

Works with: Julia version 0.6
function cramersolve(A::Matrix, b::Vector)
    return collect(begin B = copy(A); B[:, i] = b; det(B) end for i in eachindex(b)) ./ det(A)
end

A = [2 -1  5  1
     3  2  2 -6
     1  3  3 -1
     5 -2 -3  3]

b = [-3, -32, -47, 49]

@show cramersolve(A, b)
Output:
cramersolve(A, b) = [2.0, -12.0, -4.0, 1.0]

Note that it is entirely impractical to use Cramer's rule in this situation. It would be much better to use the built-in operator for solving linear systems. Assuming that the coefficient matrix and constant vector are defined as before, the solution vector is given by:

@show A \ b

Kotlin

As in the case of the Matrix arithmetic task, I've used the Johnson-Trotter permutations generator to assist with the calculation of the determinants for the various matrices:

// version 1.1.3

typealias Vector = DoubleArray
typealias Matrix = Array<Vector>

fun johnsonTrotter(n: Int): Pair<List<IntArray>, List<Int>> {
    val p = IntArray(n) { it }  // permutation
    val q = IntArray(n) { it }  // inverse permutation
    val d = IntArray(n) { -1 }  // direction = 1 or -1
    var sign = 1
    val perms = mutableListOf<IntArray>()
    val signs = mutableListOf<Int>()

    fun permute(k: Int) {
        if (k >= n) {
            perms.add(p.copyOf())
            signs.add(sign)
            sign *= -1
            return
        }
        permute(k + 1)
        for (i in 0 until k) {
            val z = p[q[k] + d[k]]
            p[q[k]] = z
            p[q[k] + d[k]] = k
            q[z] = q[k]
            q[k] += d[k]
            permute(k + 1)
        }
        d[k] *= -1
    }

    permute(0)
    return perms to signs
}

fun determinant(m: Matrix): Double {
    val (sigmas, signs) = johnsonTrotter(m.size)
    var sum = 0.0
    for ((i, sigma) in sigmas.withIndex()) {
        var prod = 1.0
        for ((j, s) in sigma.withIndex()) prod *= m[j][s]
        sum += signs[i] * prod
    }
    return sum
}

fun cramer(m: Matrix, d: Vector): Vector {
    val divisor = determinant(m)
    val numerators = Array(m.size) { Matrix(m.size) { m[it].copyOf() } }
    val v = Vector(m.size)
    for (i in 0 until m.size) {
        for (j in 0 until m.size) numerators[i][j][i] = d[j]
    }
    for (i in 0 until m.size) v[i] = determinant(numerators[i]) / divisor
    return v
}

fun main(args: Array<String>) {
    val m = arrayOf(
        doubleArrayOf(2.0, -1.0,  5.0,  1.0),
        doubleArrayOf(3.0,  2.0,  2.0, -6.0),
        doubleArrayOf(1.0,  3.0,  3.0, -1.0),
        doubleArrayOf(5.0, -2.0, -3.0,  3.0)
    )
    val d = doubleArrayOf(-3.0, -32.0, -47.0, 49.0)
    val (w, x, y, z) = cramer(m, d)
    println("w = $w, x = $x, y = $y, z = $z")
}
Output:
w = 2.0, x = -12.0, y = -4.0, z = 1.0

Lua

local matrix = require "matrix" -- https://github.com/davidm/lua-matrix

local function cramer(mat, vec)
  -- Check if matrix is quadratic
  assert(#mat == #mat[1], "Matrix is not square!")
  -- Check if vector has the same size of the matrix
  assert(#mat == #vec, "Vector has not the same size of the matrix!")
	
  local size = #mat
  local main_det = matrix.det(mat)
  
  local aux_mats = {}
  local dets = {}
  local result = {}
  for i = 1, size do
    -- Construct the auxiliary matrixes
    aux_mats[i] = matrix.copy(mat)
    for j = 1, size do
      aux_mats[i][j][i] = vec[j]
    end
    
    -- Calculate the auxiliary determinants
    dets[i] = matrix.det(aux_mats[i])
    
    -- Calculate results
    result[i] = dets[i]/main_det
  end
  
  return result
end

-----------------------------------------------

local A = {{ 2, -1,  5,  1},
           { 3,  2,  2, -6},
           { 1,  3,  3, -1},
           { 5, -2, -3,  3}}
local b = {-3, -32, -47, 49}

local result = cramer(A, b)
print("Result: " .. table.concat(result, ", "))
Output:
Result: 2, -12, -4, 1

Maple

with(LinearAlgebra):
cramer:=proc(A,B)
  local n,d,X,V,i;
  n:=upperbound(A,2);
  d:=Determinant(A);
  X:=Vector(n,0);
  for i from 1 to n do
    V:=A(1..-1,i);
    A(1..-1,i):=B;
    X[i]:=Determinant(A)/d;
    A(1..-1,i):=V;
  od;
  X;
end:

A:=Matrix([[2,-1,5,1],[3,2,2,-6],[1,3,3,-1],[5,-2,-3,3]]):
B:=Vector([-3,-32,-47,49]):
printf("%a",cramer(A,B));
Output:
Vector(4, [2,-12,-4,1])

Mathematica/Wolfram Language

crule[m_, b_] := Module[{d = Det[m], a},
  Table[a = m; a[[All, k]] = b; Det[a]/d, {k, Length[m]}]]

crule[{
  {2, -1, 5, 1},
  {3, 2, 2, -6},
  {1, 3, 3, -1},
  {5, -2, -3, 3}
 } , {-3, -32, -47, 49}]
Output:
{2,-12,-4,1}

Maxima

(%i1) eqns: [ 2*w-x+5*y+z=-3, 3*w+2*x+2*y-6*z=-32, w+3*x+3*y-z=-47, 5*w-2*x-3*y+3*z=49];
(%o1) [z + 5 y - x + 2 w = - 3, (- 6 z) + 2 y + 2 x + 3 w = - 32, 
                      (- z) + 3 y + 3 x + w = - 47, 3 z - 3 y - 2 x + 5 w = 49]
(%i2) A: augcoefmatrix (eqns, [w,x,y,z]);
                          [ 2  - 1   5    1    3   ]
                          [                        ]
                          [ 3   2    2   - 6   32  ]
(%o2)                     [                        ]
                          [ 1   3    3   - 1   47  ]
                          [                        ]
                          [ 5  - 2  - 3   3   - 49 ]
(%i3) C: coefmatrix(eqns, [w,x,y,z]);
                             [ 2  - 1   5    1  ]
                             [                  ]
                             [ 3   2    2   - 6 ]
(%o3)                        [                  ]
                             [ 1   3    3   - 1 ]
                             [                  ]
                             [ 5  - 2  - 3   3  ]
(%i4) c[n]:= (-1)^(n+1) * determinant (submatrix (A,n))/determinant (C);
                            n + 1
                       (- 1)      determinant(submatrix(A, n))
(%o4)            c  := ---------------------------------------
                  n                determinant(C)
(%i5) makelist (c[n],n,1,4);
(%o5)                          [2, - 12, - 4, 1]
(%i6) linsolve(eqns, [w,x,y,z]);
(%o6)                  [w = 2, x = - 12, y = - 4, z = 1]

Nim

Translation of: Phix
type
  SquareMatrix[N: static Positive] = array[N, array[N, float]]
  Vector[N: static Positive] = array[N, float]


####################################################################################################
# Templates.

template `[]`(m: SquareMatrix; i, j: Natural): float =
  ## Allow to get value of an element using m[i, j] syntax.
  m[i][j]

template `[]=`(m: var SquareMatrix; i, j: Natural; val: float) =
  ## Allow to set value of an element using m[i, j] syntax.
  m[i][j] = val

#---------------------------------------------------------------------------------------------------

func det(m: SquareMatrix): float =
  ## Return the determinant of matrix "m".

  var m = m
  result = 1

  for j in 0..m.high:
    var imax = j
    for i in (j + 1)..m.high:
      if m[i, j] > m[imax, j]:
        imax = i

    if imax != j:
      swap m[iMax], m[j]
      result = -result

    if abs(m[j, j]) < 1e-12:
      return NaN

    for i in (j + 1)..m.high:
      let mult = -m[i, j] / m[j, j]
      for k in 0..m.high:
        m[i, k] += mult * m[j, k]

  for i in 0..m.high:
    result *= m[i, i]

#---------------------------------------------------------------------------------------------------

func cramerSolve(a: SquareMatrix; detA: float; b: Vector; col: Natural): float =
  ## Apply Cramer rule on matrix "a", using vector "b" to replace column "col".

  when a.N != b.N:
    {.error: "incompatible matrix and vector sizes".}

  else:
    var a = a
    for i in 0..a.high:
      a[i, col] = b[i]
    result = det(a) / detA

#———————————————————————————————————————————————————————————————————————————————————————————————————

import strformat

const

  A: SquareMatrix[4] = [[2.0, -1.0,  5.0,  1.0],
                        [3.0,  2.0,  2.0, -6.0],
                        [1.0,  3.0,  3.0, -1.0],
                        [5.0, -2.0, -3.0,  3.0]]

  B: Vector[4] = [-3.0, -32.0, -47.0, 49.0]

let detA = det(A)
if detA == NaN:
  echo "Singular matrix!"
  quit(QuitFailure)

for i in 0..A.high:
  echo &"{cramerSolve(A, detA, B, i):7.3f}"
Output:
  2.000
-12.000
 -4.000
  1.000

PARI/GP

M = [2,-1,5,1;3,2,2,-6;1,3,3,-1;5,-2,-3,3];
V = Col([-3,-32,-47,49]);

matadjoint(M) * V / matdet(M)

Output:

[2, -12, -4, 1]~

Perl

use Math::Matrix;

sub cramers_rule {
    my ($A, $terms) = @_;
    my @solutions;
    my $det = $A->determinant;
    foreach my $i (0 .. $#{$A}) {
        my $Ai = $A->clone;
        foreach my $j (0 .. $#{$terms}) {
            $Ai->[$j][$i] = $terms->[$j];
        }
        push @solutions, $Ai->determinant / $det;
    }
    @solutions;
}

my $matrix = Math::Matrix->new(
    [2, -1,  5,  1],
    [3,  2,  2, -6],
    [1,  3,  3, -1],
    [5, -2, -3,  3],
);

my $free_terms = [-3, -32, -47, 49];
my ($w, $x, $y, $z) = cramers_rule($matrix, $free_terms);

print "w = $w\n";
print "x = $x\n";
print "y = $y\n";
print "z = $z\n";
Output:
w = 2
x = -12
y = -4
z = 1

Phix

Translation of: C

The copy-on-write semantics of Phix really shine here; because there is no explicit return/re-assign, you can treat parameters as a private workspace, confident in the knowledge that the updated version will be quietly discarded; all the copying and freeing of the C version is automatic/unnecessary here.
UPDATE: For the next release, "with js" (or "with javascript_semantics") diables said copy-on-write semantics, so this now needs a couple of deep_copy() calls.

requires("0.8.4")
with javascript_semantics

constant inf = 1e300*1e300,
         nan = -(inf/inf)
 
function det(sequence a)
    atom res = 1
 
    a = deep_copy(a)
    integer l = length(a)
    for j=1 to l do
        integer i_max = j
        for i=j+1 to l do
            if a[i][j] > a[i_max][j] then
                i_max = i
            end if
        end for
        if i_max != j then
            sequence aim = a[i_max]
            a[i_max] = a[j]
            a[j] = aim
            res *= -1
        end if
 
        if abs(a[j][j]) < 1e-12 then
            puts(1,"Singular matrix!")
            return nan
        end if
 
        for i=j+1 to l do
            atom mult = -a[i][j] / a[j][j]
            for k=1 to l do
                a[i][k] += mult * a[j][k]
            end for
        end for
    end for
 
    for i=1 to l do
        res *= a[i][i]
    end for
    return res
end function
 
function cramer_solve(sequence a, atom det_a, sequence b, integer v)
    a = deep_copy(a)
    for i=1 to length(a) do
      a[i][v] = b[i]
    end for
    return det(a)/det_a
end function
 
sequence a = {{2,-1, 5, 1},
              {3, 2, 2,-6},
              {1, 3, 3,-1},
              {5,-2,-3, 3}},
         b = {-3,-32,-47,49}
integer det_a = det(a)
for i=1 to length(a) do
    printf(1, "%7.3f\n", cramer_solve(a, det_a, b, i))
end for
Output:
  2.000
-12.000
 -4.000
  1.000

Prolog

Works with: GNU Prolog
removeElement([_|Tail], 0, Tail).
removeElement([Head|Tail], J, [Head|X]) :-
    J_2 is J - 1,
    removeElement(Tail, J_2, X).

removeColumn([], _, []).
removeColumn([Matrix_head|Matrix_tail], J, [X|Y]) :-
    removeElement(Matrix_head, J, X),
    removeColumn(Matrix_tail, J, Y).

removeRow([_|Matrix_tail], 0, Matrix_tail).
removeRow([Matrix_head|Matrix_tail], I, [Matrix_head|X]) :-
    I_2 is I - 1,
    removeRow(Matrix_tail, I_2, X).

cofactor(Matrix, I, J, X) :-
    removeRow(Matrix, I, Matrix_2),
    removeColumn(Matrix_2, J, Matrix_3),
    det(Matrix_3, Y),
    X is (-1) ** (I + J) * Y.

det_summand(_, _, [], 0).
det_summand(Matrix, J, B, X) :-
    B = [B_head|B_tail],
    cofactor(Matrix, 0, J, Z),
    J_2 is J + 1,
    det_summand(Matrix, J_2, B_tail, Y),
    X is B_head * Z + Y.

det([[X]], X).
det(Matrix, X) :-
    Matrix = [Matrix_head|_],
    det_summand(Matrix, 0, Matrix_head, X).

replaceElement([_|Tail], 0, New, [New|Tail]).
replaceElement([Head|Tail], J, New, [Head|Y]) :-
    J_2 is J - 1,
    replaceElement(Tail, J_2, New, Y).

replaceColumn([], _, _, []).
replaceColumn([Matrix_head|Matrix_tail], J, [Column_head|Column_tail], [X|Y]) :-
    replaceElement(Matrix_head, J, Column_head, X),
    replaceColumn(Matrix_tail, J, Column_tail, Y).

cramerElements(_, B, L, []) :- length(B, L).
cramerElements(A, B, J, [X_J|Others]) :-
    replaceColumn(A, J, B, A_J),
    det(A_J, Det_A_J),
    det(A, Det_A),
    X_J is Det_A_J / Det_A,
    J_2 is J + 1,
    cramerElements(A, B, J_2, Others).

cramer(A, B, X) :- cramerElements(A, B, 0, X).

results(X) :-
    A = [
            [2, -1,  5,  1],
            [3,  2,  2, -6],
            [1,  3,  3, -1],
            [5, -2, -3,  3]
        ],
    B = [-3, -32, -47, 49],
    cramer(A, B, X).
Output:
| ?- results(X).
X = [2.0,-12.0,-4.0,1.0] ? 
yes

Python

def det(m,n):
 if n==1: return m[0][0]
 z=0
 for r in range(n):
  k=m[:]
  del k[r]
  z+=m[r][0]*(-1)**r*det([p[1:]for p in k],n-1)
 return z
w=len(t)
d=det(h,w)
if d==0:r=[]
else:r=[det([r[0:i]+[s]+r[i+1:]for r,s in zip(h,t)],w)/d for i in range(w)]
print(r)

Racket

#lang racket
(require math/matrix)

(define sys
  (matrix [[2 -1 5 1]
           [3 2 2 -6]
           [1 3 3 -1]
           [5 -2 -3 3]]))

(define soln
  (col-matrix [-3 -32 -47 49]))

(define (matrix-set-column M new-col idx)
  (matrix-augment (list-set (matrix-cols M) idx new-col)))

(define (cramers-rule M soln)
  (let ([denom (matrix-determinant M)]
        [nvars (matrix-num-cols M)])
    (letrec ([roots (λ (position)
                      (if (>= position nvars)
                          '()
                          (cons (/ (matrix-determinant
                                    (matrix-set-column M soln position))
                                   denom)
                                (roots (add1 position)))))])
      (map cons '(w x y z) (roots 0)))))

(cramers-rule sys soln)
Output:
'((w . 2) (x . -12) (y . -4) (z . 1))

Raku

(formerly Perl 6)

sub det(@matrix) {
    my @a     = @matrix.map: { [|$_] };
    my $sign  = 1;
    my $pivot = 1;
    for ^@a -> \k {
      my @r = (k+1 ..^ @a);
      my $previous-pivot = $pivot;
      if 0 == ($pivot = @a[k;k]) {
        (my \s = @r.first: { @a[$_;k] != 0 }) // return 0;
        (@a[s], @a[k]) = (@a[k], @a[s]);
        my $pivot = @a[k;k];
        $sign = -$sign;
      }
      for @r X @r -> (\i,\j) {
        ((@a[i;j] ×= $pivot) -= @a[i;k@a[k;j]) /= $previous-pivot;
      }
    }
    $sign × $pivot
}

sub cramers_rule(@A, @terms) {
    gather for ^@A -> \i {
        my @Ai = @A.map: { [|$_] };
        for ^@terms -> \j {
            @Ai[j;i] = @terms[j];
        }
        take det(@Ai);
    } »/» det(@A);
}

my @matrix = (
    [2, -1,  5,  1],
    [3,  2,  2, -6],
    [1,  3,  3, -1],
    [5, -2, -3,  3],
);

my @free_terms = <-3 -32 -47 49>;
my ($w, $x, $y, $z) = cramers_rule(@matrix, @free_terms);
("w = $w", "x = $x", "y = $y", "z = $z").join("\n").say;
Output:
w = 2
x = -12
y = -4
z = 1

REXX

version 1

/* REXX Use Cramer's rule to compute solutions of given linear equations  */
Numeric Digits 20
names='w x y z'
M='  2  -1   5   1',
  '  3   2   2  -6',
  '  1   3   3  -1',
  '  5  -2  -3   3'
v=' -3',
  '-32',
  '-47',
  ' 49'
Call mk_mat(m)                      /* M -> a.i.j                    */
Do j=1 To dim                       /* Show the input                */
  ol=''
  Do i=1 To dim
    ol=ol format(a.i.j,6)
    End
  ol=ol format(word(v,j),6)
  Say ol
  End
Say copies('-',35)

d=det(m)                            /* denominator determinant       */

Do k=1 To dim                       /* construct nominator matrix    */
  Do j=1 To dim
    Do i=1 To dim
      If i=k Then
        b.i.j=word(v,j)
      Else
        b.i.j=a.i.j
      End
    End
  Call show_b
  d.k=det(mk_str())                 /* numerator determinant         */
  Say word(names,k) '=' d.k/d       /* compute value of variable k   */
  End
Exit

mk_mat: Procedure Expose a. dim     /* Turn list into matrix a.i.j */
  Parse Arg list
  dim=sqrt(words(list))
  k=0
  Do j=1 To dim
    Do i=1 To dim
      k=k+1
      a.i.j=word(list,k)
      End
    End
  Return

mk_str: Procedure Expose b. dim     /* Turn matrix b.i.j into list   */
  str=''
Do j=1 To dim
  Do i=1 To dim
    str=str b.i.j
    End
  End
Return str

show_b: Procedure Expose b. dim     /* show numerator matrix         */
  do j=1 To dim
    ol=''
    Do i=1 To dim
      ol=ol format(b.i.j,6)
      end
    Call dbg ol
    end
  Return

det: Procedure                      /* compute determinant           */
Parse Arg list
n=words(list)
call dbg 'det:' list
do dim=1 To 10
  If dim**2=n Then Leave
  End
call dbg 'dim='dim
If dim=2 Then Do
  det=word(list,1)*word(list,4)-word(list,2)*word(list,3)
  call dbg 'det=>'det
  Return det
  End
k=0
Do j=1 To dim
  Do i=1 To dim
    k=k+1
    a.i.j=word(list,k)
    End
  End
Do j=1 To dim
  ol=j
  Do i=1 To dim
    ol=ol format(a.i.j,6)
    End
  call dbg ol
  End
det=0
Do i=1 To dim
  ol=''
  Do j=2 To dim
    Do ii=1 To dim
      If ii<>i Then
        ol=ol a.ii.j
      End
    End
  call dbg 'i='i 'ol='ol
  If i//2 Then
    det=det+a.i.1*det(ol)
  Else
    det=det-a.i.1*det(ol)
  End
Call dbg 'det=>>>'det
Return det
sqrt: Procedure
/* REXX ***************************************************************
* EXEC to calculate the square root of a = 2 with high precision
**********************************************************************/
  Parse Arg x,prec
  If prec<9 Then prec=9
  prec1=2*prec
  eps=10**(-prec1)
  k = 1
  Numeric Digits 3
  r0= x
  r = 1
  Do i=1 By 1 Until r=r0 | (abs(r*r-x)<eps)
    r0 = r
    r  = (r + x/r) / 2
    k  = min(prec1,2*k)
    Numeric Digits (k + 5)
    End
  Numeric Digits prec
  r=r+0
  Return r


dbg: Return
Output:
      2     -1      5      1     -3
      3      2      2     -6    -32
      1      3      3     -1    -47
      5     -2     -3      3     49
-----------------------------------
w = 2
x = -12
y = -4
z = 1

version 2

The REXX version is based on the REXX version 1 program with the following improvements:

  •   aligns all output
  •   shows the   values   of the linear equations
  •   uses a PARSE for finding some matrix elements
  •   allows larger matrices to be used
  •   finds the widest decimal numbers for better formatting instead of assuming six digits
  •   use true variable names, doesn't assume that there are only four variables
  •   uses exact comparisons where appropriate
  •   added a check to see if the matrix had all its elements specified, added an error message
  •   uses a faster form of   DO   loops
  •   elided dead code and superfluous statements
  •   elided the need for a high precision sqrt function
  •   eschewed the use of a variable name with a function with the same name   (bad practice)
  •   eschewed the deprecated use of:     call func(args)     syntax
  •   automatically used the minimum width when showing the matrix elements and equation values
/*REXX program uses Cramer's rule to find and display solution of given linear equations*/
values=     '-3 -32 -47 49'                      /*values of each matrix row of numbers.*/
variables= substr('abcdefghijklmnopqrstuvwxyz', 27 - words(values) )   /*variable names.*/
call makeM  ' 2  -1   5  1        3   2   2  -6        1   3   3  -1        5  -2  -3   3'
              do   y=1  for sz;  $=              /*display the matrix (linear equations)*/
                do x=1  for sz;  $= $ right(psign(@.x.y), w)'*'substr(variables, x, 1)
                end   /*y*/                      /* [↑]   right─justify matrix elements.*/
              pad= left('', length($) - 2);    say $   ' = '   right( word(values, y), wv)
              end     /*x*/                      /* [↑]   obtain value of the equation. */
say; say
              do     k=1  for sz                 /*construct the nominator matrix.      */
                do   j=1  for sz
                  do i=1  for sz;  if i==k  then !.i.j= word(values, j)
                                            else !.i.j= @.i.j
                  end   /*i*/
                end     /*j*/
              say pad substr(variables,k,1) ' = ' right(det(makeL())/det(mat), digits()+2)
              end       /*k*/
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
makeL: $=; do x=1  for sz; do y=1  for sz; $= $ !.x.y; end; end; return $ /*matrix─►list*/
mSize: arg _; do sz=0 for 1e3; if sz*sz==_ then return; end; say 'error,bad matrix';exit 9
psign: parse arg num;  if left(num, 1)\=='-'  &  x>1  then return "+"num;   return num
/*──────────────────────────────────────────────────────────────────────────────────────*/
det: procedure;  parse arg a b c d 1 nums;        call mSize words(nums);    _= 0
     if sz==2  then return a*d  -  b*c
              do   j=1  for sz
                do i=1  for sz;    _= _ + 1;      @.i.j= word(nums, _)
                end   /*i*/
              end
     aa= 0
              do     i=1  for sz;  odd= i//2;     $=
                do   j=2  for sz-1
                  do k=1  for sz;  if k\==i  then $= $  @.k.j
                  end   /*k*/
                end     /*j*/
              aa= aa   -   (-1 ** odd)  *  @.i.1  *  det($)
              end;      /*i*/;                                               return aa
/*──────────────────────────────────────────────────────────────────────────────────────*/
makeM: procedure expose @. values mat sz w wv;  parse arg mat;    call mSize words(mat)
       #= 0;                     wv= 0;                           w= 0
              do   j=1  for sz;  wv= max(wv, length( word( values, j) ) )
                do k=1  for sz;  #= #+1;  @.k.j= word(mat, #);    w= max(w, length(@.k.j))
                end   /*k*/
              end;    /*j*/;     w= w + 1;                                   return
output   when using the internal default inputs:
   2*w  -1*x  +5*y  +1*z  =   -3
   3*w  +2*x  +2*y  -6*z  =  -32
   1*w  +3*x  +3*y  -1*z  =  -47
   5*w  -2*x  -3*y  +3*z  =   49


                       w  =            2
                       x  =          -12
                       y  =           -4
                       z  =            1

RPL

Works with: Halcyon Calc version 4.2.7

Explicit Cramer's rule

IF OVER DET
   THEN 
     LAST ROT DUP SIZE 1 GET → vect det mtx dim 
     ≪ 1 dim FOR c 
          mtx 
          1 dim FOR r 
             r c 2 →LIST vect r GET PUT 
          NEXT 
          DET det / 
       NEXT 
       dim →ARRY 
     ≫ 
  END 
≫
‘CRAMR’ STO
[[ 2 -1 5 1 ][ 3 2 2 -6 ][ 1 3 3 -1 ][ 5 -2 -3 3 ]]
[ -3 -32 -47 49 ] CRAMR
Output:
1: [ 2 -12 -4 1 ]

Implicit Cramer's rule

RPL use Cramer's rule for its built-in equation system resolution feature, performed by the / instruction.

[ -3 -32 -47 49 ]
[[ 2 -1 5 1 ][ 3 2 2 -6 ][ 1 3 3 -1 ][ 5 -2 -3 3 ]]
/
Output:
1: [ 2 -12 -4 1 ]

Ruby

require 'matrix'
 
def cramers_rule(a, terms)
  raise ArgumentError, " Matrix not square"  unless a.square?
  cols = a.to_a.transpose
  cols.each_index.map do |i|
    c = cols.dup
    c[i] = terms
    Matrix.columns(c).det / a.det
  end
end

matrix = Matrix[
    [2, -1,  5,  1],
    [3,  2,  2, -6],
    [1,  3,  3, -1],
    [5, -2, -3,  3],
]

vector = [-3, -32, -47, 49]
puts cramers_rule(matrix, vector)
Output:
2
-12
-4
1

Rust

use std::ops::{Index, IndexMut};

fn main() {
    let m = matrix(
        vec![
            2., -1., 5., 1., 3., 2., 2., -6., 1., 3., 3., -1., 5., -2., -3., 3.,
        ],
        4,
    );
    let mm = m.solve(&vec![-3., -32., -47., 49.]);
    println!("{:?}", mm);
}

#[derive(Clone)]
struct Matrix {
    elts: Vec<f64>,
    dim: usize,
}

impl Matrix {
    // Compute determinant using cofactor method
    // Using Gaussian elimination would have been more efficient, but it also solves the linear
    // system, so…
    fn det(&self) -> f64 {
        match self.dim {
            0 => 0.,
            1 => self[0][0],
            2 => self[0][0] * self[1][1] - self[0][1] * self[1][0],
            d => {
                let mut acc = 0.;
                let mut signature = 1.;
                for k in 0..d {
                    acc += signature * self[0][k] * self.comatrix(0, k).det();
                    signature *= -1.
                }
                acc
            }
        }
    }

    // Solve linear systems using Cramer's method
    fn solve(&self, target: &Vec<f64>) -> Vec<f64> {
        let mut solution: Vec<f64> = vec![0.; self.dim];
        let denominator = self.det();
        for j in 0..self.dim {
            let mut mm = self.clone();
            for i in 0..self.dim {
                mm[i][j] = target[i]
            }
            solution[j] = mm.det() / denominator
        }
        solution
    }

    // Compute the cofactor matrix for determinant computations
    fn comatrix(&self, k: usize, l: usize) -> Matrix {
        let mut v: Vec<f64> = vec![];
        for i in 0..self.dim {
            for j in 0..self.dim {
                if i != k && j != l {
                    v.push(self[i][j])
                }
            }
        }
        matrix(v, self.dim - 1)
    }
}

fn matrix(elts: Vec<f64>, dim: usize) -> Matrix {
    assert_eq!(elts.len(), dim * dim);
    Matrix { elts, dim }
}

impl Index<usize> for Matrix {
    type Output = [f64];

    fn index(&self, i: usize) -> &Self::Output {
        let m = self.dim;
        &self.elts[m * i..m * (i + 1)]
    }
}

impl IndexMut<usize> for Matrix {
    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
        let m = self.dim;
        &mut self.elts[m * i..m * (i + 1)]
    }
}

Which outputs:

[2.0, -12.0, -4.0, 1.0]

Sidef

func cramers_rule(A, terms) {
    gather {
        for i in ^A {
            var Ai = A.map{.map{_}}
            for j in ^terms {
                Ai[j][i] = terms[j]
            }
            take(Ai.det)
        }
    } »/» A.det
}

var matrix = [
    [2, -1,  5,  1],
    [3,  2,  2, -6],
    [1,  3,  3, -1],
    [5, -2, -3,  3],
]

var free_terms = [-3, -32, -47, 49]
var (w, x, y, z) = cramers_rule(matrix, free_terms)...

say "w = #{w}"
say "x = #{x}"
say "y = #{y}"
say "z = #{z}"
Output:
w = 2
x = -12
y = -4
z = 1

Tcl

package require math::linearalgebra
namespace path ::math::linearalgebra

# Setting matrix to variable A and size to n
set A [list { 2 -1  5  1} { 3  2  2 -6} { 1  3  3 -1} { 5 -2 -3  3}]
set n [llength $A]
# Setting right side of equation
set right {-3 -32 -47 49}

# Calculating determinant of A
set detA [det $A]

# Apply Cramer's rule
for {set i 0} {$i < $n} {incr i} {
  set tmp $A                    ;# copy A to tmp
  setcol tmp $i $right          ;# replace column i with right side vector
  set detTmp [det $tmp]         ;# calculate determinant of tmp
  set v [expr $detTmp / $detA]  ;# divide two determinants
  puts [format "%0.4f" $v]      ;# format and display result
}
Output:
2.0000
-12.0000
-4.0000
1.0000

VBA

Sub CramersRule()
    OrigM = [{2,  -1, 5, 1; 3,2,2,-6;1,3,3,-1;5,-2,-3,3}]
    OrigD = [{-3;-32;-47;49}]
    
    MatrixSize = UBound(OrigM)
    DetOrigM = WorksheetFunction.MDeterm(OrigM)
    
    For i = 1 To MatrixSize
        ChangeM = OrigM
        
        For j = 1 To MatrixSize
            ChangeM(j, i) = OrigD(j, 1)
        Next j
        
        DetChangeM = WorksheetFunction.MDeterm(ChangeM)
        Debug.Print i & ": " & DetChangeM / DetOrigM
    Next i
End Sub
Output:
1: 2
2: -12
3: -4
4: 1

Visual Basic .NET

Translation of: C#
Imports System.Runtime.CompilerServices
Imports System.Linq.Enumerable

Module Module1
    <Extension()>
    Function DelimitWith(Of T)(source As IEnumerable(Of T), Optional seperator As String = " ") As String
        Return String.Join(seperator, source)
    End Function

    Private Class SubMatrix
        Private ReadOnly source As Integer(,)
        Private ReadOnly prev As SubMatrix
        Private ReadOnly replaceColumn As Integer()

        Public Sub New(source As Integer(,), replaceColumn As Integer())
            Me.source = source
            Me.replaceColumn = replaceColumn
            prev = Nothing
            ColumnIndex = -1
            Size = replaceColumn.Length
        End Sub

        Public Sub New(prev As SubMatrix, Optional deletedColumnIndex As Integer = -1)
            source = Nothing
            replaceColumn = Nothing
            Me.prev = prev
            ColumnIndex = deletedColumnIndex
            Size = prev.Size - 1
        End Sub

        Public Property ColumnIndex As Integer
        Public ReadOnly Property Size As Integer

        Default Public ReadOnly Property Index(row As Integer, column As Integer) As Integer
            Get
                If Not IsNothing(source) Then
                    Return If(column = ColumnIndex, replaceColumn(row), source(row, column))
                Else
                    Return prev(row + 1, If(column < ColumnIndex, column, column + 1))
                End If
            End Get
        End Property

        Public Function Det() As Integer
            If Size = 1 Then Return Me(0, 0)
            If Size = 2 Then Return Me(0, 0) * Me(1, 1) - Me(0, 1) * Me(1, 0)
            Dim m As New SubMatrix(Me)
            Dim detVal = 0
            Dim sign = 1
            For c = 0 To Size - 1
                m.ColumnIndex = c
                Dim d = m.Det()
                detVal += Me(0, c) * d * sign
                sign = -sign
            Next
            Return detVal
        End Function

        Public Sub Print()
            For r = 0 To Size - 1
                Dim rl = r
                Console.WriteLine(Range(0, Size).Select(Function(c) Me(rl, c)).DelimitWith(", "))
            Next
            Console.WriteLine()
        End Sub
    End Class

    Private Function Solve(matrix As SubMatrix) As Integer()
        Dim det = matrix.Det()
        If det = 0 Then Throw New ArgumentException("The determinant is zero.")

        Dim answer(matrix.Size - 1) As Integer
        For i = 0 To matrix.Size - 1
            matrix.ColumnIndex = i
            answer(i) = matrix.Det() / det
        Next
        Return answer
    End Function

    Public Function SolveCramer(equations As Integer()()) As Integer()
        Dim size = equations.Length
        If equations.Any(Function(eq) eq.Length <> size + 1) Then Throw New ArgumentException($"Each equation must have {size + 1} terms.")
        Dim matrix(size - 1, size - 1) As Integer
        Dim column(size - 1) As Integer
        For r = 0 To size - 1
            column(r) = equations(r)(size)
            For c = 0 To size - 1
                matrix(r, c) = equations(r)(c)
            Next
        Next
        Return Solve(New SubMatrix(matrix, column))
    End Function

    Sub Main()
        Dim equations = {
            ({2, -1, 5, 1, -3}),
            ({3, 2, 2, -6, -32}),
            ({1, 3, 3, -1, -47}),
            ({5, -2, -3, 3, 49})
        }
        Dim solution = SolveCramer(equations)
        Console.WriteLine(solution.DelimitWith(", "))
    End Sub

End Module
Output:
2, -12, -4, 1

Wren

Library: Wren-matrix
import "./matrix" for Matrix

var cramer = Fn.new { |a, d|
    var n = a.numRows
    var x = List.filled(n, 0)
    var ad = a.det
    for (c in 0...n) {
        var aa = a.copy()
        for (r in 0...n) aa[r, c] = d[r, 0]
        x[c] = aa.det/ad
    }
    return x
}

var a = Matrix.new([
    [2, -1,  5,  1],
    [3,  2,  2, -6],
    [1,  3,  3, -1],
    [5, -2, -3,  3]
])

var d = Matrix.new([
    [- 3],
    [-32],
    [-47],
    [ 49]
])

var x = cramer.call(a, d)
System.print("Solution is %(x)")
Output:
Solution is [2, -12, -4, 1]

XPL0

func Det(A, N);         \Return value of determinate A, order N
int  A, N;
int  B, Sum, I, K, L, Term;
[if N = 1 then return A(0, 0);
B:= Reserve((N-1)*4\IntSize\);
Sum:= 0;
for I:= 0 to N-1 do
    [L:= 0;
    for K:= 0 to N-1 do
      if K # I then
        [B(L):= @A(K, 1);  L:= L+1];
    Term:= A(I, 0) * Det(B, N-1);
    if I & 1 then Term:= -Term;
    Sum:= Sum + Term;
    ];
return Sum;
];

real D;
[D:= float(Det([[2,-1,5,1], [3,2,2,-6], [1,3,3,-1], [5,-2,-3,3]], 4));
RlOut(0, float(Det([[-3,-1,5,1], [-32,2,2,-6], [-47,3,3,-1], [49,-2,-3,3]], 4)) / D);
RlOut(0, float(Det([[2,-3,5,1], [3,-32,2,-6], [1,-47,3,-1], [5,49,-3,3]], 4)) / D);
RlOut(0, float(Det([[2,-1,-3,1], [3,2,-32,-6], [1,3,-47,-1], [5,-2,49,3]], 4)) / D);
RlOut(0, float(Det([[2,-1,5,-3], [3,2,2,-32], [1,3,3,-47], [5,-2,-3,49]], 4)) / D);
]
Output:
    2.00000  -12.00000   -4.00000    1.00000

zkl

Using the GNU Scientific Library, we define the values:

var [const] GSL=Import("zklGSL");	// libGSL (GNU Scientific Library)
A:=GSL.Matrix(4,4).set(2,-1, 5, 1,
		       3, 2, 2,-6,
		       1, 3, 3,-1,
		       5,-2,-3, 3);
b:=GSL.Vector(4).set(-3,-32,-47,49);

First, just let GSL solve:

A.AxEQb(b).format().println();

Actually implement Cramer's rule:

Translation of: Julia
fcn cramersRule(A,b){
   b.len().pump(GSL.Vector(b.len()),'wrap(i){ // put calculations into new Vector
      A.copy().setColumn(i,b).det();
   }).close()/A.det();
}
cramersRule(A,b).format().println();
Output:
2.00,-12.00,-4.00,1.00