UNPKG

4.22 kBtext/coffeescriptView Raw
1#-----------------------------------------------------------------------------
2#
3# Input: Matrix on stack (must have two dimensions but
4# it can be non-numerical)
5#
6# Output: Inverse on stack
7#
8# Example:
9#
10# > inv(((1,2),(3,4))
11# ((-2,1),(3/2,-1/2))
12#
13# > inv(((a,b),(c,d))
14# ((d / (a d - b c),-b / (a d - b c)),(-c / (a d - b c),a / (a d - b c)))
15#
16# Note:
17#
18# THIS IS DIFFERENT FROM INVERSE OF AN EXPRESSION (inv)
19# Uses Gaussian elimination for numerical matrices.
20#
21#-----------------------------------------------------------------------------
22
23
24
25INV_check_arg = ->
26 if (!istensor(p1))
27 return 0
28 else if (p1.tensor.ndim != 2)
29 return 0
30 else if (p1.tensor.dim[0] != p1.tensor.dim[1])
31 return 0
32 else
33 return 1
34
35inv = ->
36 i = 0
37 n = 0
38 #U **a
39
40 save()
41
42 p1 = pop()
43
44 # an inv just goes away when
45 # applied to another inv
46 if (isinv(p1))
47 push car(cdr(p1))
48 restore()
49 return
50
51 # inverse goes away in case
52 # of identity matrix
53 if isidentitymatrix(p1)
54 push p1
55 restore()
56 return
57
58 # distribute the inverse of a dot
59 # if in expanding mode
60 # note that the distribution happens
61 # in reverse.
62 # The dot operator is not
63 # commutative, so, it matters.
64 if (expanding && isinnerordot(p1))
65 p1 = cdr(p1)
66 accumulator = []
67 while (iscons(p1))
68 accumulator.push car(p1)
69 p1 = cdr(p1)
70
71 for eachEntry in [accumulator.length-1..0]
72 push(accumulator[eachEntry])
73 inv()
74 if eachEntry != accumulator.length-1
75 inner()
76
77 restore()
78 return
79
80
81 if (INV_check_arg() == 0)
82 push_symbol(INV)
83 push(p1)
84 list(2)
85 restore()
86 return
87
88
89 if isnumerictensor p1
90 yyinvg()
91 else
92 push(p1)
93 adj()
94 push(p1)
95 det()
96 p2 = pop()
97 if (iszero(p2))
98 stop("inverse of singular matrix")
99 push(p2)
100 divide()
101
102 restore()
103
104invg = ->
105 save()
106
107 p1 = pop()
108
109 if (INV_check_arg() == 0)
110 push_symbol(INVG)
111 push(p1)
112 list(2)
113 restore()
114 return
115
116 yyinvg()
117
118 restore()
119
120# inverse using gaussian elimination
121
122yyinvg = ->
123 h = 0
124 i = 0
125 j = 0
126 n = 0
127
128 n = p1.tensor.dim[0]
129
130 h = tos
131
132 for i in [0...n]
133 for j in [0...n]
134 if (i == j)
135 push(one)
136 else
137 push(zero)
138
139 for i in [0...(n * n)]
140 push(p1.tensor.elem[i])
141
142 INV_decomp(n)
143
144 p1 = alloc_tensor(n * n)
145
146 p1.tensor.ndim = 2
147 p1.tensor.dim[0] = n
148 p1.tensor.dim[1] = n
149
150 for i in [0...(n * n)]
151 p1.tensor.elem[i] = stack[h + i]
152
153 moveTos tos - 2 * n * n
154
155 push(p1)
156
157#-----------------------------------------------------------------------------
158#
159# Input: n * n unit matrix on stack
160#
161# n * n operand on stack
162#
163# Output: n * n inverse matrix on stack
164#
165# n * n garbage on stack
166#
167# p2 mangled
168#
169#-----------------------------------------------------------------------------
170
171#define A(i, j) stack[a + n * (i) + (j)]
172#define U(i, j) stack[u + n * (i) + (j)]
173
174INV_decomp = (n) ->
175 a = 0
176 d = 0
177 i = 0
178 j = 0
179 u = 0
180
181 a = tos - n * n
182
183 u = a - n * n
184
185 for d in [0...n]
186
187 # diagonal element zero?
188
189 if (equal( (stack[a + n * (d) + (d)]) , zero))
190
191 # find a new row
192
193 for i in [(d + 1)...n]
194 if (!equal( (stack[a + n * (i) + (d)]) , zero))
195 break
196
197 if (i == n)
198 stop("inverse of singular matrix")
199
200 # exchange rows
201
202 for j in [0...n]
203
204 p2 = stack[a + n * (d) + (j)]
205 stack[a + n * (d) + (j)] = stack[a + n * (i) + (j)]
206 stack[a + n * (i) + (j)] = p2
207
208 p2 = stack[u + n * (d) + (j)]
209 stack[u + n * (d) + (j)] = stack[u + n * (i) + (j)]
210 stack[u + n * (i) + (j)] = p2
211
212 # multiply the pivot row by 1 / pivot
213
214 p2 = stack[a + n * (d) + (d)]
215
216 for j in [0...n]
217
218 if (j > d)
219 push(stack[a + n * (d) + (j)])
220 push(p2)
221 divide()
222 stack[a + n * (d) + (j)] = pop()
223
224 push(stack[u + n * (d) + (j)])
225 push(p2)
226 divide()
227 stack[u + n * (d) + (j)] = pop()
228
229 # clear out the column above and below the pivot
230
231 for i in [0...n]
232
233 if (i == d)
234 continue
235
236 # multiplier
237
238 p2 = stack[a + n * (i) + (d)]
239
240 # add pivot row to i-th row
241
242 for j in [0...n]
243
244 if (j > d)
245 push(stack[a + n * (i) + (j)])
246 push(stack[a + n * (d) + (j)])
247 push(p2)
248 multiply()
249 subtract()
250 stack[a + n * (i) + (j)] = pop()
251
252 push(stack[u + n * (i) + (j)])
253 push(stack[u + n * (d) + (j)])
254 push(p2)
255 multiply()
256 subtract()
257 stack[u + n * (i) + (j)] = pop()