This commit is contained in:
Ingy döt Net 2013-06-05 21:47:54 +00:00
parent 1f1ad49427
commit 6f050a029e
2496 changed files with 37609 additions and 3031 deletions

View file

@ -0,0 +1,71 @@
!-*- mode: compilation; default-directory: "/tmp/" -*-
!Compilation started at Sat May 18 23:25:42
!
!a=./F && make $a && $a < unixdict.txt
!f95 -Wall -ffree-form F.F -o F
! j example, determinant: 7.00000000
! j example, permanent: 5.00000000
! maxima, determinant: -360.000000
! maxima, permanent: 900.000000
!
!Compilation finished at Sat May 18 23:25:43
! NB. example computed by J
! NB. fixed seed random matrix
! _2+3 3?.@$5
! 2 _1 1
!_1 _2 1
!_1 _1 _1
!
! (-/ .*)_2+3 3?.@$5 NB. determinant
!7
! (+/ .*)_2+3 3?.@$5 NB. permanent
!5
!maxima example
!a: matrix([2, 9, 4], [7, 5, 3], [6, 1, 8])$
!determinant(a);
!-360
!
!permanent(a);
!900
! compute permanent or determinant
program f
implicit none
real, dimension(3,3) :: j, m
data j/ 2,-1, 1,-1,-2, 1,-1,-1,-1/
data m/2, 9, 4, 7, 5, 3, 6, 1, 8/
write(6,*) 'j example, determinant: ',det(j,3,-1)
write(6,*) 'j example, permanent: ',det(j,3,1)
write(6,*) 'maxima, determinant: ',det(m,3,-1)
write(6,*) 'maxima, permanent: ',det(m,3,1)
contains
recursive function det(a,n,permanent) result(accumulation)
! setting permanent to 1 computes the permanent.
! setting permanent to -1 computes the determinant.
real, dimension(n,n), intent(in) :: a
integer, intent(in) :: n, permanent
real, dimension(n-1, n-1) :: b
real :: accumulation
integer :: i, sgn
if (n .eq. 1) then
accumulation = a(1,1)
else
accumulation = 0
sgn = 1
do i=1, n
b(:, :(i-1)) = a(2:, :i-1)
b(:, i:) = a(2:, i+1:)
accumulation = accumulation + sgn * a(1, i) * det(b, n-1, permanent)
sgn = sgn * permanent
enddo
endif
end function det
end program f

View file

@ -0,0 +1,53 @@
sub insert( $x, @xs) { [@xs[0..$_-1], $x, @xs[$_..*]] for 0..@xs }
sub order ($sg, @xs) { $sg > 0 ?? @xs !! @xs.reverse }
multi σ_permutations ([]) { [] => 1 }
multi σ_permutations ([$x, *@xs]) {
σ_permutations(@xs).map({ order($_.value, insert($x, $_.key)) }) Z=> (1,-1) xx *
}
sub m_arith ( @a, $op ) {
note "Not a square matrix" and return
if [||] map { @a.elems cmp @a[$_].elems }, ^@a;
[+] map {
my $permutation = .key;
my $term = $op eq 'perm' ?? 1 !! .value;
for $permutation.kv -> $i, $j { $term *= @a[$i][$j] };
$term
}, σ_permutations [^@a];
}
########### Testing ###########
my @tests = (
[
[ 1, 2 ],
[ 3, 4 ]
],
[
[ 1, 2, 3, 4 ],
[ 4, 5, 6, 7 ],
[ 7, 8, 9, 10 ],
[ 10, 11, 12, 13 ]
],
[
[ 0, 1, 2, 3, 4 ],
[ 5, 6, 7, 8, 9 ],
[ 10, 11, 12, 13, 14 ],
[ 15, 16, 17, 18, 19 ],
[ 20, 21, 22, 23, 24 ]
]
);
sub dump (@matrix) {
say $_».fmt: "%3s" for @matrix, '';
}
for @tests -> @matrix {
say 'Matrix:';
@matrix.&dump;
say "Determinant:\t", @matrix.&m_arith: <det>;
say "Permanent: \t", @matrix.&m_arith: <perm>;
say '-' x 25;
}

View file

@ -0,0 +1,39 @@
/* REXX ***************************************************************
* Test the two functions determinant and permanent
* using the matrix specifications shown for other languages
* 21.05.2013 Walter Pachl
**********************************************************************/
Call test ' 1 2',
' 3 4',2
Call test ' 1 2 3 4',
' 4 5 6 7',
' 7 8 9 10',
'10 11 12 13',4
Call test ' 0 1 2 3 4',
' 5 6 7 8 9',
'10 11 12 13 14',
'15 16 17 18 19',
'20 21 22 23 24',5
Exit
test:
/**********************************************************************
* Show the given matrix and compute and show determinant and permanent
**********************************************************************/
Parse Arg as,n
asc=as
Do i=1 To n
ol=''
Do j=1 To n
Parse Var asc a.i.j asc
ol=ol right(a.i.j,3)
End
Say ol
End
Say 'determinant='right(determinant(as),7)
Say ' permanent='right(permanent(as),7)
Say copies('-',50)
Return

View file

@ -0,0 +1,60 @@
/* REXX ***************************************************************
* determinant.rex
* compute the determinant of the given square matrix
* Input: as: the representation of the matrix as vector (n**2 elements)
* 21.05.2013 Walter Pachl
**********************************************************************/
Parse Arg as
n=sqrt(words(as))
Do i=1 To n
Do j=1 To n
Parse Var as a.i.j as
End
End
Select
When n=2 Then det=a.1.1*a.2.2-a.1.2*a.2.1
When n=3 Then det= a.1.1*a.2.2*a.3.3,
+a.1.2*a.2.3*a.3.1,
+a.1.3*a.2.1*a.3.2,
-a.1.3*a.2.2*a.3.1,
-a.1.2*a.2.1*a.3.3,
-a.1.1*a.2.3*a.3.2
Otherwise Do
det=0
Do k=1 To n
det=det+((-1)**(k+1))*a.1.k*determinant(subm(k))
End
End
End
Return det
subm: Procedure Expose a. n
/**********************************************************************
* compute the submatrix resulting when row 1 and column k are removed
* Input: a.*.*, k
* Output: bs the representation of the submatrix as vector
**********************************************************************/
Parse Arg k
bs=''
do i=2 To n
Do j=1 To n
If j=k Then Iterate
bs=bs a.i.j
End
End
Return bs
sqrt: Procedure
/**********************************************************************
* compute and return the (integer) square root of the given argument
* terminate the program if the argument is not a square
**********************************************************************/
Parse Arg nn
Do n=1 By 1 while n*n<nn
End
If n*n=nn Then
Return n
Else Do
Say 'invalid number of elements:' nn 'is not a square.'
Exit
End

View file

@ -0,0 +1,105 @@
/* REXX ***************************************************************
* permanent.rex
* compute the permanent of a matrix
* I found an algorithm here:
* http://www.codeproject.com/Articles/21282/Compute-Permanent-of-a-Matrix-with-Ryser-s-Algorit
* see there for the original author.
* translated it to REXX (hopefully correctly) to REXX
* and believe that I can "publish" it here, on rosettacode
* when I look at the copyright rules shown there:
* http://www.codeproject.com/info/cpol10.aspx
* 20.05.2013 Walter Pachl
**********************************************************************/
Call init arg(1) /* initialize the matrix (n and a.* */
sum=0
rowsumprod=0
rowsum=0
chi.=0
c=2**n
Do k=1 To c-1 /* loop all 2^n submatrices of A */
rowsumprod = 1
chis=dec2binarr(k,n) /* characteristic vector */
Do ci=0 By 1 While chis<>''
Parse Var chis chi.ci chis
End
Do m=0 To n-1 /* loop columns of submatrix #k */
rowsum = 0
Do p=0 To n-1 /* loop rows and compute rowsum */
mnp=m*n+p
rowsum=rowsum+chi.p*A.mnp
End
rowsumprod=rowsumprod*rowsum /* update product of rowsums */
/* (optional -- use for sparse matrices) */
/* if (rowsumprod == 0) break; */
End
sum=sum+((-1)**(n-chi.n))*rowsumprod
End
Return sum
/**********************************************************************
* Notes
* 1.The submatrices are chosen by use of a characteristic vector chi
* (only the columns are considered, where chi[p] == 1).
* To retrieve the t from Ryser's formula, we need to save the number
* n-t, as is done in chi[n]. Then we get t = n - chi[n].
* 2.The matrix parameter A is expected to be a one-dimensional integer
* array -- should the matrix be encoded row-wise or column-wise?
* -- It doesn't matter. The permanent is invariant under
* row-switching and column-switching, and it is Screenshot
* - per_inv.gif .
* 3.Further enhancements: If any rowsum equals zero,
* the entire rowsumprod becomes zero, and thus the m-loop can be broken.
* Since if-statements are relatively expensive compared to integer
* operations, this might save time only for sparse matrices
* (where most entries are zeros).
* 4.If anyone finds a polynomial algorithm for permanents,
* he will get rich and famous (at least in the computer science world).
**********************************************************************/
/**********************************************************************
* At first, we need to transform a decimal to a binary array
* with an additional element
* (the last one) saving the number of ones in the array:
**********************************************************************/
dec2binarr: Procedure
Parse Arg n,dim
ol='n='n 'dim='dim
res.=0
pos=dim-1
Do While n>0
res.pos=n//2
res.dim=res.dim+res.pos
n=n%2
pos=pos-1
End
res_s=''
Do i=0 To dim
res_s=res_s res.i
End
Return res_s
init: Procedure Expose a. n
/**********************************************************************
* a.* (starting with index 0) contains all array elements
* n is the dimension of the square matrix
**********************************************************************/
Parse Arg as
n=sqrt(words(as))
a.=0
Do ai=0 By 1 While as<>''
Parse Var as a.ai as
End
Return
sqrt: Procedure
/**********************************************************************
* compute and return the (integer) square root of the given argument
* terminate the program if the argument is not a square
**********************************************************************/
Parse Arg nn
Do n=1 By 1 while n*n<nn
End
If n*n=nn Then
Return n
Else Do
Say 'invalid number of elements:' nn 'is not a square.'
Exit
End

View file

@ -0,0 +1,9 @@
#lang racket
(require math)
(define determinant matrix-determinant)
(define (permanent M)
(define n (matrix-num-rows M))
(for/sum ([σ (in-permutations (range n))])
(for/product ([i n] [σi σ])
(matrix-ref M i σi))))