1 |
|
2 |
|
3 |
|
4 |
|
5 |
|
6 |
|
7 |
|
8 |
|
9 |
|
10 |
|
11 |
|
12 |
|
13 |
|
14 |
|
15 |
|
16 |
|
17 |
|
18 |
|
19 |
|
20 |
|
21 |
|
22 |
|
23 |
|
24 |
|
25 | INV_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 |
|
35 | inv = ->
|
36 | i = 0
|
37 | n = 0
|
38 |
|
39 |
|
40 | save()
|
41 |
|
42 | p1 = pop()
|
43 |
|
44 |
|
45 |
|
46 | if (isinv(p1))
|
47 | push car(cdr(p1))
|
48 | restore()
|
49 | return
|
50 |
|
51 |
|
52 |
|
53 | if isidentitymatrix(p1)
|
54 | push p1
|
55 | restore()
|
56 | return
|
57 |
|
58 |
|
59 |
|
60 |
|
61 |
|
62 |
|
63 |
|
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 |
|
104 | invg = ->
|
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 |
|
121 |
|
122 | yyinvg = ->
|
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 |
|
160 |
|
161 |
|
162 |
|
163 |
|
164 |
|
165 |
|
166 |
|
167 |
|
168 |
|
169 |
|
170 |
|
171 |
|
172 |
|
173 |
|
174 | INV_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 |
|
188 |
|
189 | if (equal( (stack[a + n * (d) + (d)]) , zero))
|
190 |
|
191 |
|
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 |
|
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 |
|
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 |
|
230 |
|
231 | for i in [0...n]
|
232 |
|
233 | if (i == d)
|
234 | continue
|
235 |
|
236 |
|
237 |
|
238 | p2 = stack[a + n * (i) + (d)]
|
239 |
|
240 |
|
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()
|