UNPKG

5.77 kBtext/coffeescriptView Raw
1###
2 Symbolic addition
3
4 Terms in a sum are combined if they are identical modulo rational
5 coefficients.
6
7 For example, A + 2A becomes 3A.
8
9 However, the sum A + sqrt(2) A is not modified.
10
11 Combining terms can lead to second-order effects.
12
13 For example, consider the case of
14
15 1/sqrt(2) A + 3/sqrt(2) A + sqrt(2) A
16
17 The first two terms are combined to yield 2 sqrt(2) A.
18
19 This result can now be combined with the third term to yield
20
21 3 sqrt(2) A
22###
23
24
25
26flag = 0
27
28Eval_add = ->
29 h = tos
30 p1 = cdr(p1)
31 while (iscons(p1))
32 push(car(p1))
33 Eval()
34 p2 = pop()
35 push_terms(p2)
36 p1 = cdr(p1)
37 add_terms(tos - h)
38
39# Add n terms, returns one expression on the stack.
40
41stackAddsCount = 0
42add_terms = (n) ->
43 stackAddsCount++
44 i = 0
45
46 h = tos - n
47
48 s = h
49
50 # ensure no infinite loop, use "for"
51
52 if DEBUG then console.log "stack before adding terms #" + stackAddsCount
53 #if stackAddsCount == 137
54 # debugger
55
56 if DEBUG
57 for i in [0...tos]
58 console.log print_list stack[i]
59
60 for i in [0...10]
61
62 if (n < 2)
63 break
64
65 flag = 0
66
67
68 #qsort(s, n, sizeof (U *), cmp_terms)
69 subsetOfStack = stack.slice(h,h+n)
70 subsetOfStack.sort(cmp_terms)
71 stack = stack.slice(0,h).concat(subsetOfStack).concat(stack.slice(h+n))
72
73 if (flag == 0)
74 break
75
76 n = combine_terms(h, n)
77
78 moveTos h + n
79
80 switch (n)
81 when 0
82 if evaluatingAsFloats then push_double(0.0) else push(zero)
83 when 1
84 else
85 list(n)
86 p1 = pop()
87 push_symbol(ADD)
88 push(p1)
89 cons()
90
91 if DEBUG then console.log "stack after adding terms #" + stackAddsCount
92 #if stackAddsCount == 5
93 # debugger
94
95 if DEBUG
96 for i in [0...tos]
97 console.log print_list stack[i]
98
99# Compare terms for order, clobbers p1 and p2.
100
101cmp_terms_count = 0
102cmp_terms = (p1, p2) ->
103 cmp_terms_count++
104 #if cmp_terms_count == 52
105 # debugger
106
107 i = 0
108 # numbers can be combined
109
110 if (isnum(p1) && isnum(p2))
111 flag = 1
112 #if DEBUG then console.log "cmp_terms #" + cmp_terms_count + " returns 0"
113 return 0
114
115 # congruent tensors can be combined
116
117 if (istensor(p1) && istensor(p2))
118 if (p1.tensor.ndim < p2.tensor.ndim)
119 #if DEBUG then console.log "cmp_terms #" + cmp_terms_count + " returns -1"
120 return -1
121 if (p1.tensor.ndim > p2.tensor.ndim)
122 #if DEBUG then console.log "cmp_terms #" + cmp_terms_count + " returns 1"
123 return 1
124 for i in [0...p1.tensor.ndim]
125 if (p1.tensor.dim[i] < p2.tensor.dim[i])
126 #if DEBUG then console.log "cmp_terms #" + cmp_terms_count + " returns -1"
127 return -1
128 if (p1.tensor.dim[i] > p2.tensor.dim[i])
129 #if DEBUG then console.log "cmp_terms #" + cmp_terms_count + " returns 1"
130 return 1
131 flag = 1
132 #if DEBUG then console.log "cmp_terms #" + cmp_terms_count + " returns 0"
133 return 0
134
135 if (car(p1) == symbol(MULTIPLY))
136 p1 = cdr(p1)
137 if (isnum(car(p1)))
138 p1 = cdr(p1)
139 if (cdr(p1) == symbol(NIL))
140 p1 = car(p1)
141
142 if (car(p2) == symbol(MULTIPLY))
143 p2 = cdr(p2)
144 if (isnum(car(p2)))
145 p2 = cdr(p2)
146 if (cdr(p2) == symbol(NIL))
147 p2 = car(p2)
148
149 t = cmp_expr(p1, p2)
150
151 if (t == 0)
152 flag = 1
153
154 #if DEBUG then console.log "cmp_terms #" + cmp_terms_count + " returns " + t
155 return t
156
157###
158 Compare adjacent terms in s[] and combine if possible.
159
160 Returns the number of terms remaining in s[].
161
162 n number of terms in s[] initially
163###
164
165combine_terms = (s, n) ->
166
167 #debugger
168
169 # I had to turn the coffeescript for loop into
170 # a more mundane while loop because the i
171 # variable was changed from within the body,
172 # which is something that is not supposed to
173 # happen in the coffeescript 'vector' form.
174 # Also this means I had to add a 'i++' jus before
175 # the end of the body and before the "continue"s
176 i=0
177 while i < (n-1)
178 check_esc_flag()
179
180 p3 = stack[s+i]
181 p4 = stack[s+i + 1]
182
183 if (istensor(p3) && istensor(p4))
184 push(p3)
185 push(p4)
186 tensor_plus_tensor()
187 p1 = pop()
188 if (p1 != symbol(NIL))
189 stack[s+i] = p1
190 for j in [(i + 1)...(n - 1)]
191 stack[s+j] = stack[s+j + 1]
192 n--
193 i--
194
195 i++
196 continue
197
198 if (istensor(p3) || istensor(p4))
199
200 i++
201 continue
202
203 if (isnum(p3) && isnum(p4))
204 push(p3)
205 push(p4)
206 add_numbers()
207 p1 = pop()
208 if (iszero(p1))
209 for j in [i...(n - 2)]
210 stack[s+j] = stack[s+j + 2]
211 n -= 2
212 else
213 stack[s+i] = p1
214 for j in [(i + 1)...(n - 1)]
215 stack[s+j] = stack[s+j + 1]
216 n--
217 i--
218
219 i++
220 continue
221
222 if (isnum(p3) || isnum(p4))
223
224 i++
225 continue
226
227 if evaluatingAsFloats
228 p1 = one_as_double
229 p2 = one_as_double
230 else
231 p1 = one
232 p2 = one
233
234 t = 0
235
236 if (car(p3) == symbol(MULTIPLY))
237 p3 = cdr(p3)
238 t = 1; # p3 is now denormal
239 if (isnum(car(p3)))
240 p1 = car(p3)
241 p3 = cdr(p3)
242 if (cdr(p3) == symbol(NIL))
243 p3 = car(p3)
244 t = 0
245
246 if (car(p4) == symbol(MULTIPLY))
247 p4 = cdr(p4)
248 if (isnum(car(p4)))
249 p2 = car(p4)
250 p4 = cdr(p4)
251 if (cdr(p4) == symbol(NIL))
252 p4 = car(p4)
253
254 if (!equal(p3, p4))
255
256 i++
257 continue
258
259 push(p1)
260 push(p2)
261 add_numbers()
262
263 p1 = pop()
264
265 if (iszero(p1))
266 for j in [i...(n - 2)]
267 stack[s+j] = stack[s+j + 2]
268 n -= 2
269 i--
270
271 i++
272 continue
273
274 push(p1)
275
276 if (t)
277 push(symbol(MULTIPLY))
278 push(p3)
279 cons()
280 else
281 push(p3)
282
283 multiply()
284
285 stack[s+i] = pop()
286
287 for j in [(i + 1)...(n - 1)]
288 stack[s+j] = stack[s+j + 1]
289
290 n--
291 i--
292
293 # this i++ is to match the while
294 i++
295
296 return n
297
298push_terms = (p) ->
299 if (car(p) == symbol(ADD))
300 p = cdr(p)
301 while (iscons(p))
302 push(car(p))
303 p = cdr(p)
304 else if (!iszero(p))
305 push(p)
306
307# add two expressions
308
309add = ->
310 save()
311 p2 = pop()
312 p1 = pop()
313 h = tos
314 push_terms(p1)
315 push_terms(p2)
316 add_terms(tos - h)
317 restore()
318
319add_all = (k) ->
320 i = 0
321 save()
322 s = tos - k
323 h = tos
324 for i in [0...k]
325 push_terms(stack[s+i])
326 add_terms(tos - h)
327 p1 = pop()
328 moveTos tos - k
329 push(p1)
330 restore()
331
332subtract = ->
333 negate()
334 add()