290 lines
4.3 KiB
CoffeeScript
290 lines
4.3 KiB
CoffeeScript
### det =====================================================================
|
|
|
|
Tags
|
|
----
|
|
scripting, JS, internal, treenode, general concept
|
|
|
|
Parameters
|
|
----------
|
|
m
|
|
|
|
General description
|
|
-------------------
|
|
Returns the determinant of matrix m.
|
|
Uses Gaussian elimination for numerical matrices.
|
|
|
|
Example:
|
|
|
|
det(((1,2),(3,4)))
|
|
> -2
|
|
|
|
###
|
|
|
|
|
|
DET_check_arg = ->
|
|
if (!istensor(p1))
|
|
return 0
|
|
else if (p1.tensor.ndim != 2)
|
|
return 0
|
|
else if (p1.tensor.dim[0] != p1.tensor.dim[1])
|
|
return 0
|
|
else
|
|
return 1
|
|
|
|
det = ->
|
|
i = 0
|
|
n = 0
|
|
#U **a
|
|
|
|
save()
|
|
|
|
p1 = pop()
|
|
|
|
if (DET_check_arg() == 0)
|
|
push_symbol(DET)
|
|
push(p1)
|
|
list(2)
|
|
restore()
|
|
return
|
|
|
|
n = p1.tensor.nelem
|
|
|
|
a = p1.tensor.elem
|
|
|
|
for i in [0...n]
|
|
if (!isNumericAtom(a[i]))
|
|
break
|
|
|
|
if (i == n)
|
|
yydetg()
|
|
else
|
|
for i in [0...p1.tensor.nelem]
|
|
push(p1.tensor.elem[i])
|
|
determinant(p1.tensor.dim[0])
|
|
|
|
restore()
|
|
|
|
# determinant of n * n matrix elements on the stack
|
|
|
|
determinant = (n) ->
|
|
h = 0
|
|
i = 0
|
|
j = 0
|
|
k = 0
|
|
q = 0
|
|
s = 0
|
|
sign_ = 0
|
|
t = 0
|
|
|
|
a = []
|
|
#int *a, *c, *d
|
|
|
|
h = tos - n * n
|
|
|
|
#a = (int *) malloc(3 * n * sizeof (int))
|
|
|
|
#if (a == NULL)
|
|
# out_of_memory()
|
|
|
|
for i in [0...n]
|
|
a[i] = i
|
|
a[i+n] = 0
|
|
a[i+n+n] = 1
|
|
|
|
sign_ = 1
|
|
|
|
push(zero)
|
|
|
|
while 1
|
|
|
|
if (sign_ == 1)
|
|
push_integer(1)
|
|
else
|
|
push_integer(-1)
|
|
|
|
for i in [0...n]
|
|
k = n * a[i] + i
|
|
push(stack[h + k])
|
|
multiply(); # FIXME -- problem here
|
|
|
|
add()
|
|
|
|
# next permutation (Knuth's algorithm P)
|
|
|
|
j = n - 1
|
|
s = 0
|
|
|
|
breakFromOutherWhile = false
|
|
while 1
|
|
q = a[n+j] + a[n+n+j]
|
|
if (q < 0)
|
|
a[n+n+j] = -a[n+n+j]
|
|
j--
|
|
continue
|
|
if (q == j + 1)
|
|
if (j == 0)
|
|
breakFromOutherWhile = true
|
|
break
|
|
s++
|
|
a[n+n+j] = -a[n+n+j]
|
|
j--
|
|
continue
|
|
break
|
|
|
|
if breakFromOutherWhile
|
|
break
|
|
|
|
t = a[j - a[n+j] + s]
|
|
a[j - a[n+j] + s] = a[j - q + s]
|
|
a[j - q + s] = t
|
|
a[n+j] = q
|
|
|
|
sign_ = -sign_
|
|
|
|
|
|
stack[h] = stack[tos - 1]
|
|
|
|
moveTos h + 1
|
|
|
|
#-----------------------------------------------------------------------------
|
|
#
|
|
# Input: Matrix on stack
|
|
#
|
|
# Output: Determinant on stack
|
|
#
|
|
# Note:
|
|
#
|
|
# Uses Gaussian elimination which is faster for numerical matrices.
|
|
#
|
|
# Gaussian Elimination works by walking down the diagonal and clearing
|
|
# out the columns below it.
|
|
#
|
|
#-----------------------------------------------------------------------------
|
|
|
|
detg = ->
|
|
save()
|
|
|
|
p1 = pop()
|
|
|
|
if (DET_check_arg() == 0)
|
|
push_symbol(DET)
|
|
push(p1)
|
|
list(2)
|
|
restore()
|
|
return
|
|
|
|
yydetg()
|
|
|
|
restore()
|
|
|
|
yydetg = ->
|
|
i = 0
|
|
n = 0
|
|
|
|
n = p1.tensor.dim[0]
|
|
|
|
for i in [0...(n * n)]
|
|
push(p1.tensor.elem[i])
|
|
|
|
lu_decomp(n)
|
|
|
|
moveTos tos - n * n
|
|
|
|
push(p1)
|
|
|
|
#-----------------------------------------------------------------------------
|
|
#
|
|
# Input: n * n matrix elements on stack
|
|
#
|
|
# Output: p1 determinant
|
|
#
|
|
# p2 mangled
|
|
#
|
|
# upper diagonal matrix on stack
|
|
#
|
|
#-----------------------------------------------------------------------------
|
|
|
|
M = (h,n,i, j) ->
|
|
stack[h + n * (i) + (j)]
|
|
|
|
setM = (h,n,i,j,value) ->
|
|
stack[h + n * (i) + (j)] = value
|
|
|
|
lu_decomp = (n) ->
|
|
d = 0
|
|
h = 0
|
|
i = 0
|
|
j = 0
|
|
|
|
h = tos - n * n
|
|
|
|
p1 = one
|
|
|
|
for d in [0...(n - 1)]
|
|
|
|
# diagonal element zero?
|
|
|
|
if (equal(M(h,n,d, d), zero))
|
|
|
|
# find a new row
|
|
|
|
for i in [(d + 1)...n]
|
|
if (!equal(M(h,n,i, d), zero))
|
|
break
|
|
|
|
if (i == n)
|
|
p1 = zero
|
|
break
|
|
|
|
# exchange rows
|
|
|
|
for j in [d...n]
|
|
p2 = M(h,n,d, j)
|
|
setM(h,n,d, j, M(h,n,i, j))
|
|
setM(h,n,i, j, p2)
|
|
|
|
# negate det
|
|
|
|
push(p1)
|
|
negate()
|
|
p1 = pop()
|
|
|
|
# update det
|
|
|
|
push(p1)
|
|
push(M(h,n,d, d))
|
|
multiply()
|
|
p1 = pop()
|
|
|
|
# update lower diagonal matrix
|
|
|
|
for i in [(d + 1)...n]
|
|
|
|
# multiplier
|
|
|
|
push(M(h,n,i, d))
|
|
push(M(h,n,d, d))
|
|
divide()
|
|
negate()
|
|
|
|
p2 = pop()
|
|
|
|
# update one row
|
|
|
|
setM(h,n,i, d, zero); # clear column below pivot d
|
|
|
|
for j in [(d + 1)...n]
|
|
push(M(h,n,d, j))
|
|
push(p2)
|
|
multiply()
|
|
push(M(h,n,i, j))
|
|
add()
|
|
setM(h,n,i, j, pop())
|
|
|
|
# last diagonal element
|
|
|
|
push(p1)
|
|
push(M(h,n,n - 1, n - 1))
|
|
multiply()
|
|
p1 = pop()
|