Rosetta Code

Cramer's rule: Difference between revisions

Page Discussion

Cramer's rule (view source)

Revision as of 19:34, 6 April 2024

75,846 bytes added ,  2 months ago
Added Easylang
(Added Easylang)
 
(78 intermediate revisions by 41 users not shown)
Line 36:
solve for <big><math>w</math>, <math>x</math>, <math>y</math></big> and <big><math>z</math></big>, using Cramer's rule.
<br><br>
 
=={{header|11l}}==
{{trans|Nim}}
 
<syntaxhighlight lang="11l">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)))</syntaxhighlight>
 
{{out}}
<pre>
2.000
-12.000
-4.000
1.000
</pre>
 
=={{header|Ada}}==
<syntaxhighlight lang="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;</syntaxhighlight>
 
{{out}}
<pre>[ 2.00 -12.00 -4.00 1.00 ]</pre>
 
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
Uses the non-standard DET operator available in Algol 68G.
<langsyntaxhighlight lang="algol68"># returns the solution of a.x = b via Cramer's rule #
# this is for REAL arrays, could define additional operators #
# for INT, COMPL, etc. #
Line 87 ⟶ 217:
OD;
print( ( newline ) )
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 95 ⟶ 225:
=={{header|C}}==
 
<langsyntaxhighlight Clang="c">#include <math.h>
#include <stdio.h>
#include <stdlib.h>
Line 206 ⟶ 336:
deinit_square_matrix(A);
return EXIT_SUCCESS;
}</langsyntaxhighlight>
 
{{out}}
Line 213 ⟶ 343:
-4.000
1.000</pre>
 
=={{header|C sharp|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.
<syntaxhighlight lang="csharp">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();
}
}
 
}</syntaxhighlight>
{{out}}
<pre>
2, -12, -4, 1
</pre>
 
=={{header|C++}}==
{{trans|C#}}
<syntaxhighlight lang="cpp">#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;
}</syntaxhighlight>
{{out}}
<pre>[2, -12, -4, 1]</pre>
 
=={{header|Common Lisp}}==
<syntaxhighlight lang="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))</syntaxhighlight>
{{out}}
<pre>(2 -12 -4 1)</pre>
 
=={{header|D}}==
{{trans|Kotlin}}
<syntaxhighlight lang="d">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]);
}</syntaxhighlight>
{{out}}
<pre>w = 2, x = -12, y = -4, z = 1</pre>
 
=={{header|EasyLang}}==
{{trans|Phix}}
<syntaxhighlight>
proc det . a0[][] res .
res = 1
a[][] = a0[][]
n = len a[][]
for j to n
imax = j
for i = j + 1 to n
if a[i][j] > a[imax][j]
imax = i
.
.
if imax <> j
swap a[imax][] a[j][]
res = -res
.
if abs a[j][j] < 1e-12
print "Singular matrix!"
res = 0 / 0
return
.
for i = j + 1 to n
mult = -a[i][j] / a[j][j]
for k to n
a[i][k] += mult * a[j][k]
.
.
.
for i to n
res *= a[i][i]
.
.
proc cramer_solve . a0[][] deta b[] col res .
a[][] = a0[][]
for i to len a[][]
a[i][col] = b[i]
.
det a[][] d
res = d / deta
.
a[][] = [ [ 2 -1 5 1 ] [ 3 2 2 -6 ] [ 1 3 3 -1 ] [ 5 -2 -3 3 ] ]
b[] = [ -3 -32 -47 49 ]
det a[][] deta
for i to len a[][]
cramer_solve a[][] deta b[] i r
print r
.
</syntaxhighlight>
{{out}}
<pre>
2.00
-12
-4
1.00
</pre>
 
=={{header|EchoLisp}}==
<langsyntaxhighlight lang="scheme">
(lib 'matrix)
(string-delimiter "")
Line 232 ⟶ 807:
(writeln "B = " B)
(writeln "X = " (cramer A B)))
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 243 ⟶ 818:
B = #( -3 -32 -47 49)
X = #( 2 -12 -4 1)
</pre>
 
=={{header|Factor}}==
<syntaxhighlight lang="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</syntaxhighlight>
{{out}}
<pre>
{ 2 -12 -4 1 }
</pre>
 
Line 255 ⟶ 856:
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 <code>WRITE (6,*) A</code>, 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 <code>READ(INF,*) A</code>
 
Similarly, if instead a DATA statement were used to initialise the array for the example problem, and it looked something like <langsyntaxhighlight Fortranlang="fortran"> DATA A/2, -1, 5, 1
1 3, 2, 2, -6
2 1, 3, 3, -1
3 5, -2, -3, 3/</langsyntaxhighlight> (ignoring integer/floating-point issues) thus corresponding to the layout of the example problem, there would need to be a statement <code>A = TRANSPOSE(A)</code> to obtain the required order.
 
I have never seen an explanation of why this choice was made for Fortran.<langsyntaxhighlight Fortranlang="fortran"> MODULE BAD IDEA !Employ Cramer's rule for solving A.x = b...
INTEGER MSG !Might as well be here.
CONTAINS !The details.
Line 367 ⟶ 968:
Closedown.
100 WRITE (6,*) "That was interesting." !Quite.
END !Open files are closed, allocated memory is released. </langsyntaxhighlight>
 
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:
Line 394 ⟶ 995:
</pre>
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.
 
 
=={{header|FreeBASIC}}==
<syntaxhighlight lang="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
</syntaxhighlight>
{{out}}
<pre>
2
-12
-4
1
</pre>
 
=={{header|Go}}==
'''Library gonum:'''
<syntaxhighlight lang="go">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)
}</syntaxhighlight>
{{out}}
<pre>
[2 -12.000000000000007 -4.000000000000001 1.0000000000000009]
</pre>
'''Library go.matrix:'''
<syntaxhighlight lang="go">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)
}</syntaxhighlight>
{{out}}
<pre>
[2.0000000000000004 -11.999999999999998 -4 0.9999999999999999]
</pre>
 
=={{header|Groovy}}==
{{trans|Java}}
<syntaxhighlight lang="groovy">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)
}
}
}</syntaxhighlight>
{{out}}
<pre>Solution = [2.0, -12.0, -4.0, 1.0]</pre>
 
=={{header|Haskell}}==
===Version 1===
<langsyntaxhighlight lang="haskell">import Data.Matrix
 
solveCramer :: (Ord a, Fractional a) => Matrix a -> Matrix a -> Maybe [a]
Line 413 ⟶ 1,273:
,[1, 3, 3,-1]
,[5,-2,-3, 3]]
y = fromLists [[-3], [-32], [-47], [49]]</langsyntaxhighlight>
{{out}}
 
Line 419 ⟶ 1,279:
Just [2.0000000000000004,-11.999999999999998,-4.0,0.9999999999999999]
===Version 2===
We use Rational numbers for having more precision.
<lang Haskell>
a % b is the rational a / b.
s_permutations :: [a] -> [([a], Int)]
<syntaxhighlight lang="haskell">s_permutations :: [a] -> [([a], Int)]
s_permutations = flip zip (cycle [1, -1]) . (foldl aux [[]])
where aux items x = do
Line 427 ⟶ 1,288:
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..]
Line 436 ⟶ 1,303:
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)
 
determinantelemCramerPos::Int Num-> Int -> [[a]] =-> [[a]] -> Int -> Int -> a
determinantelemCramerPos l k ks ms i j = s_determinantif j /= l then elemPos ms $i s_permutationsj [0..pred.lengthelse $elemPos ms]ks i k
 
elemCramerPossolveCramer::Int -> [a[Rational]] -> [[aRational]] -> Int -> Int -> a[[Rational]]
solveCramer ms ks = xs
elemCramerPos l ks ms i j = if j /= l then elemPos ms i j else ks !! i
 
solveCramer::Fractional a => [[a]] -> [a] -> [a]
solveCramer ms ks = map (\l -> (/d) $ s_determinant (elemCramerPos l ks) ms ps) $ ls
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
task = solveCramer a y
wherelet a = [[2,-1, 5, 1]
,[3, 2, 2,-6]
,[1, 3, 3,-1]
,[5,-2,-3, 3]]
let b y = [[-3], [-32], [-47], [49]]
task a b
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
 
===Version 3===
<syntaxhighlight lang="haskell">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 = print taskdo
let a = [[2,-1, 5, 1]
</lang>
,[3, 2, 2,-6]
,[1, 3, 3,-1]
,[5,-2,-3, 3]]
let b = [[-3], [-32], [-47], [49]]
task a b
</syntaxhighlight>
{{out}}
<pre>
a =
[2.0,-12.0,-4.0,1.0]
[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
</pre>
 
Line 467 ⟶ 1,548:
Implementation:
 
<langsyntaxhighlight Jlang="j">cramer=:4 :0
A=. x [ b=. y
det=. -/ .*
A %~&det (i.#A) b"_`[`]}&.|:"0 2 A
)</langsyntaxhighlight>
 
Task data:
 
<langsyntaxhighlight Jlang="j">A=: _&".;._2]t=: 0 :0
2 -1 5 1
3 2 2 -6
Line 482 ⟶ 1,563:
)
 
b=: _3 _32 _47 49</langsyntaxhighlight>
 
Task example:
 
<langsyntaxhighlight Jlang="j"> A cramer b
2 _12 _4 1</langsyntaxhighlight>
 
=={{header|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.
 
<syntaxhighlight lang="java">
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);
}
}
 
}
</syntaxhighlight>
{{out}}
<pre>
Solution = [2.0, -12.0, -4.0, 1.0]
</pre>
 
=={{header|JavaScript}}==
 
<syntaxhighlight lang="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();});
}
 
</syntaxhighlight>
{{out}}
<pre>[ 2, -12, -4, 1 ]</pre>
 
=={{header|jq}}==
'''Adapted from [[#Wren|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.
<syntaxhighlight lang="jq"># 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))"</syntaxhighlight>
{{out}}
<pre>
Solution is [2,-12,-4,1]
</pre>
 
=={{header|Julia}}==
{{works with|Julia|0.6}}
<lang Julia>function cramer_solve(A, v)
<syntaxhighlight lang="julia">function cramersolve(A::Matrix, b::Vector)
return [ begin
return collect(begin B = copy(A); B[:, i] = b; det(B) end for i Bin =eachindex(b)) copy./ det(A)
det((B[:, i] = v; B))
end for i = 1 : length(v) ] ./ det(A)
end
 
A = [2 -1 5 1
2 -13 52 1;2 -6
3 1 2 3 2 3 -6;1
1 5 -2 -3 3 -1;]
5 -2 -3 3;
]
v = [-3, -32, -47, 49]
 
b = [-3, -32, -47, 49]
map(println, round(cramer_solve(A, v), 14))</lang>
 
@show cramersolve(A, b)</syntaxhighlight>
 
{{out}}
<pre>cramersolve(A, b) = [2.0, -12.0, -4.0, 1.0]</pre>
 
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:
 
<lang Julia>A \ v</lang>
<syntaxhighlight lang="julia">@show A \ b</syntaxhighlight>
 
=={{header|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:
<syntaxhighlight lang="scala">// 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")
}</syntaxhighlight>
 
{{out}}
<pre>
w = 2.0, x = -12.0, y = -4.0, z = 1.0
</pre>
 
=={{header|Lua}}==
<pre>2.0
<syntaxhighlight lang="lua">local matrix = require "matrix" -- https://github.com/davidm/lua-matrix
-12.0
 
-4.0
local function cramer(mat, vec)
1.0</pre>
-- 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, ", "))
</syntaxhighlight>
 
{{out}}
<pre>Result: 2, -12, -4, 1</pre>
 
=={{header|Maple}}==
<syntaxhighlight lang="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));</syntaxhighlight>
 
{{out}}
 
<pre>Vector(4, [2,-12,-4,1])</pre>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">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}]</syntaxhighlight>
 
{{out}}
<pre>{2,-12,-4,1}</pre>
 
=={{header|Maxima}}==
<syntaxhighlight lang="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]
</syntaxhighlight>
 
=={{header|Nim}}==
{{trans|Phix}}
 
<syntaxhighlight lang="nim">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}"</syntaxhighlight>
 
{{out}}
<pre> 2.000
-12.000
-4.000
1.000</pre>
 
=={{header|PARI/GP}}==
 
<langsyntaxhighlight lang="parigp">
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)
</syntaxhighlight>
</lang>
 
Output:
Line 530 ⟶ 2,164:
 
=={{header|Perl}}==
<langsyntaxhighlight lang="perl">use Math::Matrix;
 
sub cramers_rule {
Line 559 ⟶ 2,193:
print "x = $x\n";
print "y = $y\n";
print "z = $z\n";</langsyntaxhighlight>
{{out}}
<pre>
Line 568 ⟶ 2,202:
</pre>
 
=={{header|Perl 6Phix}}==
{{trans|C}}
<lang perl6>sub det(@matrix) {
<del>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.</del><br>
my @a = @matrix.map: { [|$_] };
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.
my $sign = +1;
<!--<syntaxhighlight lang="phix">(phixonline)-->
my $pivot = 1;
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"0.8.4"</span><span style="color: #0000FF;">)</span>
for ^@a -> $k {
<span style="color: #008080;">with</span> <span style="color: #000000;">javascript_semantics</span>
my @r = ($k+1 .. @a.end);
my $previous-pivot = $pivot;
<span style="color: #008080;">constant</span> <span style="color: #000000;">inf</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1e300</span><span style="color: #0000FF;">*</span><span style="color: #000000;">1e300</span><span style="color: #0000FF;">,</span>
if 0 == ($pivot = @a[$k][$k]) {
<span style="color: #000000;">nan</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-(</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">/</span><span style="color: #000000;">inf</span><span style="color: #0000FF;">)</span>
(my $s = @r.first: { @a[$_][$k] != 0 }) // return 0;
(@a[$s],@a[$k]) = (@a[$k], @a[$s]);
<span style="color: #008080;">function</span> <span style="color: #000000;">det</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
my $pivot = @a[$k][$k];
<span style="color: #004080;">atom</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
$sign = -$sign;
}
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
for @r X @r -> ($i, $j) {
<span style="color: #004080;">integer</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
((@a[$i][$j] *= $pivot) -= @a[$i][$k]*@a[$k][$j]) /= $previous-pivot;
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
}
<span style="color: #004080;">integer</span> <span style="color: #000000;">i_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">j</span>
}
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">j</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
$sign * $pivot
<span style="color: #008080;">if</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">></span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i_max</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
}
<span style="color: #000000;">i_max</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">i</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
sub cramers_rule(@A, @terms) {
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
gather for ^@A -> $i {
<span style="color: #008080;">if</span> <span style="color: #000000;">i_max</span> <span style="color: #0000FF;">!=</span> <span style="color: #000000;">j</span> <span style="color: #008080;">then</span>
my @Ai = @A.map: { [|$_] };
<span style="color: #004080;">sequence</span> <span style="color: #000000;">aim</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i_max</span><span style="color: #0000FF;">]</span>
for ^@terms -> $j {
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i_max</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
@Ai[$j][$i] = @terms[$j];
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">aim</span>
}
<span style="color: #000000;">res</span> <span style="color: #0000FF;">*=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span>
take det(@Ai);
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
} »/» det(@A);
}
<span style="color: #008080;">if</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">])</span> <span style="color: #0000FF;"><</span> <span style="color: #000000;">1e-12</span> <span style="color: #008080;">then</span>
 
<span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Singular matrix!"</span><span style="color: #0000FF;">)</span>
my @matrix = (
<span style="color: #008080;">return</span> <span style="color: #000000;">nan</span>
[2, -1, 5, 1],
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
[3, 2, 2, -6],
[1, 3, 3, -1],
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">j</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
[5, -2, -3, 3],
<span style="color: #004080;">atom</span> <span style="color: #000000;">mult</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
);
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
 
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">mult</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span>
my @free_terms = (-3, -32, -47, 49);
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
my ($w, $x, $y, $z) = |cramers_rule(@matrix, @free_terms);
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
 
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
say "w = $w";
say "x = $x";
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">l</span> <span style="color: #008080;">do</span>
say "y = $y";
<span style="color: #000000;">res</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
say "z = $z";</lang>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
{{out}}
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">cramer_solve</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">det_a</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">v</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">det</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">det_a</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">5</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">6</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">5</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">3</span><span style="color: #0000FF;">}},</span>
<span style="color: #000000;">b</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">32</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">47</span><span style="color: #0000FF;">,</span><span style="color: #000000;">49</span><span style="color: #0000FF;">}</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">det_a</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">det</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span> <span style="color: #008000;">"%7.3f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cramer_solve</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">det_a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<!--</syntaxhighlight>-->
{{Out}}
<pre>
w = 2.000
x = -12.000
-4.000
y = -4
z = 1.000
</pre>
 
=={{header|Prolog}}==
{{works with|GNU Prolog}}
<langsyntaxhighlight lang="prolog">removeElement([_|Tail], 0, Tail).
removeElement([Head|Tail], J, [Head|X]) :-
J_2 is J - 1,
Line 686 ⟶ 2,342:
],
B = [-3, -32, -47, 49],
cramer(A, B, X).</langsyntaxhighlight>
 
{{out}}
Line 696 ⟶ 2,352:
 
=={{header|Python}}==
<langsyntaxhighlight lang="python">
def det(m,n):
# a simple implementation using numpy
if n==1: return m[0][0]
from numpy import linalg
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)
</syntaxhighlight>
 
=={{header|Racket}}==
A=[[2,-1,5,1],[3,2,2,-6],[1,3,3,-1],[5,-2,-3,3]]
<syntaxhighlight lang="racket">
B=[-3,-32,-47,49]
#lang racket
C=[[2,-1,5,1],[3,2,2,-6],[1,3,3,-1],[5,-2,-3,3]]
(require math/matrix)
X=[]
 
for i in range(0,len(B)):
(define sys
for j in range(0,len(B)):
(matrix [[2 -1 5 C[j][i]=B[j1]
if i>0: [3 2 2 -6]
C[j][i-1]=A[j][i 3 3 -1]
[5 -2 -3 3]]))
X.append(round(linalg.det(C)/linalg.det(A),1))
 
(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)
</syntaxhighlight>
 
{{out}}
print('w=%s'%X[0],'x=%s'%X[1],'y=%s'%X[2],'z=%s'%X[3])
 
</lang>
<pre>'((w . 2) (x . -12) (y . -4) (z . 1))</pre>
 
=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku" line>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;</syntaxhighlight>
{{out}}
<pre>
w = 2
w=2.0 x=-12.0 y=-4.0 z=1.0
x = -12
y = -4
z = 1
</pre>
 
=={{header|REXX}}==
=== version 1 ===
<lang rexx>/* Use Cramer's rule to compute solutions of given linear equations */
<syntaxhighlight lang="rexx">/* REXX Use Cramer's rule to compute solutions of given linear equations */
Numeric Digits 20
names='w x y z'
Line 833 ⟶ 2,570:
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</langsyntaxhighlight>
{{out}}
<pre> 2 -1 5 1 -3
Line 846 ⟶ 2,605:
z = 1</pre>
 
=== version 2 ===
The REXX version is based on the REXX version 1 program with the following improvements:
* &nbsp; aligns all output
* &nbsp; shows the &nbsp; values &nbsp; of the linear equations
* &nbsp; uses a PARSE for finding some matrix elements
* &nbsp; allows larger matrices to be used
* &nbsp; finds the widest decimal numbers for better formatting instead of assuming six digits
* &nbsp; use true variable names, doesn't assume that there are only four variables
* &nbsp; uses exact comparisons where appropriate
* &nbsp; added a check to see if the matrix had all its elements specified, added an error message
* &nbsp; uses a faster form of &nbsp; '''DO''' &nbsp; loops
* &nbsp; elided dead code and superfluous statements
* &nbsp; elided the need for a high precision '''sqrt''' function
* &nbsp; eschewed the use of a variable name with a function with the same name &nbsp; (bad practice)
* &nbsp; eschewed the deprecated use of: &nbsp; &nbsp; ''call func(args)'' &nbsp; &nbsp; syntax
* &nbsp; automatically used the minimum width when showing the matrix elements and equation values
<syntaxhighlight lang="rexx">/*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</syntaxhighlight>
{{out|output|text=&nbsp; when using the internal default inputs:}}
<pre>
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
</pre>
 
=={{header|RPL}}==
{{works with|Halcyon Calc|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
{{out}}
<pre>
1: [ 2 -12 -4 1 ]
</pre>
 
===Implicit Cramer's rule===
RPL use Cramer's rule for its built-in equation system resolution feature, performed by the <code>/</code> instruction.
[ -3 -32 -47 49 ]
[[ 2 -1 5 1 ][ 3 2 2 -6 ][ 1 3 3 -1 ][ 5 -2 -3 3 ]]
/
{{out}}
<pre>
1: [ 2 -12 -4 1 ]
</pre>
=={{header|Ruby}}==
<langsyntaxhighlight lang="ruby">require 'matrix'
def cramers_rule(a, terms)
raise ArgumentError, " Matrix not square" unless a.square?
det = a.det
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
Line 868 ⟶ 2,736:
 
vector = [-3, -32, -47, 49]
puts cramers_rule(matrix, vector)</langsyntaxhighlight>
{{out}}
<pre>
Line 877 ⟶ 2,745:
</pre>
 
=={{header|SidefRust}}==
<syntaxhighlight lang="rust">use std::ops::{Index, IndexMut};
<lang ruby>func det(a) {
a = a.map{.map{_}}
var sign = +1
var pivot = 1
 
fn main() {
for k in ^a {
let var rm = matrix(k+1 .. a.end)
var previous_pivot = pivotvec![
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)]
if ((pivot = a[k][k]) == 0) {
struct Matrix {
a.swap(r.first_by { a[_][k] != 0 } \\ return 0, k)
elts: Vec<f64>,
pivot = a[k][k]
dim: sign.neg!usize,
}
}
 
impl Matrix {
for i,j in (r ~X r) {
// Compute determinant using cofactor method
a[i][j] *= pivot ->
// Using Gaussian elimination would have been more efficient, but it also solves the linear
-= a[i][k]*a[k][j] ->
// system, so…
/= previous_pivot
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)
}
sign * pivot
}
 
fn matrix(elts: Vec<f64>, dim: usize) -> Matrix {
func cramers_rule(A, terms) {
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)]
}
}</syntaxhighlight>
 
Which outputs:
<pre>[2.0, -12.0, -4.0, 1.0]</pre>
 
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">func cramers_rule(A, terms) {
gather {
for i in ^A {
Line 909 ⟶ 2,846:
Ai[j][i] = terms[j]
}
take(det(Ai).det)
}
} »/» det(A).det
}
 
Line 927 ⟶ 2,864:
say "x = #{x}"
say "y = #{y}"
say "z = #{z}"</langsyntaxhighlight>
{{out}}
<pre>
Line 934 ⟶ 2,871:
y = -4
z = 1
</pre>
 
=={{header|Tcl}}==
<syntaxhighlight lang="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
}
</syntaxhighlight>
{{out}}
<pre>
2.0000
-12.0000
-4.0000
1.0000
</pre>
 
=={{header|VBA}}==
<syntaxhighlight lang="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
 
</syntaxhighlight>
{{out}}
<pre>
1: 2
2: -12
3: -4
4: 1
</pre>
 
=={{header|Visual Basic .NET}}==
{{trans|C#}}
<syntaxhighlight lang="vbnet">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</syntaxhighlight>
{{out}}
<pre>2, -12, -4, 1</pre>
 
=={{header|Wren}}==
{{libheader|Wren-matrix}}
<syntaxhighlight lang="wren">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)")</syntaxhighlight>
 
{{out}}
<pre>
Solution is [2, -12, -4, 1]
</pre>
 
=={{header|XPL0}}==
<syntaxhighlight lang="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);
]</syntaxhighlight>
 
{{out}}
<pre>
2.00000 -12.00000 -4.00000 1.00000
</pre>
 
=={{header|zkl}}==
Using the GNU Scientific Library, we define the values:
<syntaxhighlight lang="zkl">var [const] GSL=Import("zklGSL"); // libGSL (GNU Scientific Library)
<lang zkl>A:=GSL.Matrix(4,4).set(2,-1, 5, 1,
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);</langsyntaxhighlight>
First, just let GSL solve:
<langsyntaxhighlight lang="zkl">A.AxEQb(b).format().println();</langsyntaxhighlight>
Actually implement Cramer's rule:
{{Trans|Julia}}
<langsyntaxhighlight lang="zkl">fcn cramersRule(A,b){
b.len().pump(GSL.Vector(b.len()).sink(),'wrap(i){ // put calculations into new Vector
A.copy().setColumn(i,b).det();
}).close()/A.det();
}
cramersRule(A,b).format().println();</langsyntaxhighlight>
{{out}}
<pre>
Chkas
2,064

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/263502...363338"