C#で行列式を計算するコードを書いてみた
ネットで見つけた順列を生成するC#のコードを参考にして、行列式を計算するコードを書いてみた。
２行２列で要素の値が０，１，２のどれかで行列式の値が１のものを表示するプログラム。２行２列で要素の値が０，１，２のどれかである行列は重複ありの順列を使って生成させてみた。
>|cs|
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Collections;
static class Extensions
{
public static IEnumerable<IEnumerable<T>> Perm<T>(this IEnumerable<T> items, int? k = null)
{
if (k == null)
k = items.Count();
if (k == 0)
{
yield return Enumerable.Empty<T>();
}
else
{
foreach (var x in items)
{
var xs = items.ToList();
xs.Remove(x);
foreach (var c in Perm(xs,k - 1)){
var newelem = (new T {x}).Concat(c).ToArray();
yield return newelem;
}
}
}
}
public static IEnumerable<IEnumerable<T>> RepeatedPerm<T>(this IEnumerable<T> items, int? k = null)
{
if (k == null)
k = items.Count();
if (k == 0)
{
yield return Enumerable.Empty<T>();
}
else
{
foreach (var x in items)
{
var xs = items.ToList();
foreach (var c in RepeatedPerm(xs,k - 1)){
var newelem = (new T {x}).Concat(c).ToArray();
yield return newelem;
}
}
}
}
}
public class Matrix
{
int dim;
public int[,] elements;
public Matrix(int order)
{
dim = order;
elements = new int[dim,dim];
}
public Matrix(int order,int[,] elm)
{
dim = order;
elements = elm;
}
public static Matrix operator +(Matrix a, Matrix b)
{
int[,] elm;
int n, i, j;
n = a.dim;
elm = new int[n,n];
for(i=0;i<n;i++){
for(j=0;j<n;j++){
elm[i,j] = a.elements[i,j] + b.elements[i,j];
}
}
return new Matrix(n,elm);
}
public static Matrix operator -(Matrix a, Matrix b)
{
int[,] elm;
int n, i, j;
n = a.dim;
elm = new int[n,n];
for(i=0;i<n;i++){
for(j=0;j<n;j++){
elm[i,j] = a.elements[i,j] - b.elements[i,j];
}
}
return new Matrix(n,elm);
}
public static Matrix operator *(Matrix a, Matrix b)
{
int[,] elm;
int n, i, j, k, sum;
n = a.dim;
elm = new int[n,n];
for(i=0;i<n;i++){
for(j=0;j<n;j++){
sum = 0;
for(k=0;k<n;k++){
sum += a.elements[i,k] * b.elements[k,j];
}
elm[i,j] = sum;
}
}
return new Matrix(n,elm);
}
public void WriteLine()
{
int i, j;
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
Console.Write(" {0,6:d0}",elements[i, j]);
}
Console.WriteLine();
}
}
//public static int signature(int src, int ar)
public static int signature(int src, IEnumerable<int> plist)
{
int cnt = 0, swp, i, j;
int ar = plist.ToArray();
int br = new int[ar.Length];
for(i=0;i<ar.Length;i++){
br[i] = ar[i];
}
for(i=0;i<src.Length;i++){
if ( src[i] != br[i] ){
for(j=i;j<src.Length;j++){
if ( br[j] == src[i] ){
swp = br[i];
br[i] = src[i];
br[j] = swp;
cnt++;
}
}
}
}
return cnt%2==0?1:(-1);
}
public int Determinant()
{
int src;
int sgn, sum, pro, i, j;
src = new int[dim];
for(i=0;i<dim;i++) src[i]=i;
sum = 0;
var perm = src.Perm();
foreach( var plist in perm )
{
pro = 1;
j = 0;
foreach( var p in plist )
{
pro *= elements[j,p];
j++;
}
sgn = signature(src,plist);
pro *= sgn;
sum += pro;
}
//Console.WriteLine("----- << det >> ----");
//Console.WriteLine("sum = {0}",sum);
return sum;
}
public Matrix Small(int n, int m)
{
int[,] small;
int i,j, c_row, c_col;
Console.WriteLine("--- small ---");
small = new int[dim-1,dim-1];
c_row = 0;
for ( i=0; i< dim; i++ ){
if ( i != n ){
c_col = 0;
for ( j=0; j< dim; j++ ){
if ( j != m ){
small[c_row,c_col] = elements[i,j];
//Console.Write(" {0}",small[c_row,c_col]);
c_col++;
}
}
//Console.WriteLine();
c_row++;
}
}
//Console.WriteLine("--------");
for ( i=0; i< dim -1; i++ ){
for ( j=0; j< dim -1; j++ ){
//Console.Write(" {0}",small[i,j]);
}
//Console.WriteLine();
}
Matrix a = new Matrix(dim-1,small);
return a;
}
public int cofactor(int n, int m)
{
int[,] small;
int sgn,det,i,j, c_row, c_col;
if ( dim < 2 ){
return 0;
}
//Console.WriteLine("--- cofactor ---");
small = new int[dim-1,dim-1];
c_row = 0;
for ( i=0; i< dim; i++ ){
if ( i != n ){
c_col = 0;
for ( j=0; j< dim; j++ ){
if ( j != m ){
small[c_row,c_col] = elements[i,j];
//Console.Write(" {0}",small[c_row,c_col]);
c_col++;
}
}
//Console.WriteLine();
c_row++;
}
}
//Console.WriteLine("--------");
for ( i=0; i< dim -1; i++ ){
for ( j=0; j< dim -1; j++ ){
//Console.Write(" {0}",small[i,j]);
}
//Console.WriteLine();
}
Matrix a = new Matrix(dim-1,small);
det = a.Determinant();
sgn = (n+m)%2==0?1:(-1);
//Console.WriteLine("< {0} : {1} > det = {2} sgn={3}",n,m,det,sgn);
return sgn*det;
}
public Matrix Cofactorize()
{
int[,] inv = new int[dim,dim];
int i,j;
for(i = 0; i < dim; i++) {
for(j = 0; j < dim; j++)
{
inv[i,j] = cofactor(j,i);
}
}
return new Matrix(dim,inv);
}
public int Rank()
{
Matrix s;
int max =0, rk, det, i, j;
if ( dim == 1 ){
if ( elements[0,0] == 0 ){
return 0;
}
else {
return 1;
}
}
det = Determinant();
Console.WriteLine("det = {0}",det);
if ( det != 0 ){
return dim;
}
else{
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
s = Small(i,j);
Console.WriteLine("< {0} : {1} > ",i,j);
s.WriteLine();
//det = s.Determinant();
//Console.WriteLine("det = {0}",det);
rk = s.Rank();
if ( rk > max ){
max = rk;
}
}
}
}
Console.WriteLine("max = {0}",max);
return max;
}
public void Unit()
{
for(int i = 0; i < dim; i++) {
elements[i,i] = 1;
}
}
public void ZeroClear()
{
int i,j;
for(i = 0; i < dim; i++) {
for(j = 0; j < dim; j++)
{
elements[i,j] = 0;
}
}
}
public void Jordan()
{
int i,j;
for(i = 0; i < dim; i++) {
for(j = 0; j < dim; j++)
{
if ( j == (i+1) ){
elements[i,j] = 1;
}
else {
elements[i,j] = 0;
}
}
}
}
}
class Program
{
static void Main(string args)
{
int order = 2;
int i, j;
int cnt = 0;
Matrix A;
int[,] matrix = new int[order, order];
var source = new int { 0, 1, 2, };
//var source = new int { -1, 0, 1 };
//var source = new int { 1, 2, };
var duplex = source.RepeatedPerm( order * order );
foreach( var dlist in duplex )
{
i = 0;
j = 0;
foreach( var d in dlist )
{
if ( j > 0 && j % order == 0 )
{
i++;
j=0;
}
matrix[i,j] = d;
j++;
}
A = new Matrix(order,matrix);
if ( A.Determinant() == 1 ){
Console.WriteLine("[{0}]",++cnt);
A.WriteLine();
Console.WriteLine();
}
}
Console.WriteLine("--- det 1 matrix :{0} / case num = {1} ---",cnt,duplex.Count());
}
}
||<
◎実行例
２行２列で要素の値が０，１，２のどれかで行列式の値が１のものを表示する。
PS >.\matrix.exe
[1]
1 0
0 1
[2]
1 0
1 1
[3]
1 0
2 1
[4]
1 1
0 1
[5]
1 1
1 2
[6]
1 2
0 1
[7]
2 1
1 1
--- det 1 matrix :7 / case num = 81 ---
