106 lines
3.8 KiB
Rexx
106 lines
3.8 KiB
Rexx
/* 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
|