UNPKG

8.48 kBtext/coffeescriptView Raw
1# Partial fraction expansion
2#
3# Example
4#
5# expand(1/(x^3+x^2),x)
6#
7# 1 1 1
8# ---- - --- + -------
9# 2 x x + 1
10# x
11
12
13
14Eval_expand = ->
15
16 # 1st arg
17
18 push(cadr(p1))
19 Eval()
20
21 # 2nd arg
22
23 push(caddr(p1))
24 Eval()
25
26 p2 = pop()
27 if (p2 == symbol(NIL))
28 guess()
29 else
30 push(p2)
31
32 expand()
33
34#define A p2
35#define B p3
36#define C p4
37#define F p5
38#define P p6
39#define Q p7
40#define T p8
41#define X p9
42
43expand = ->
44
45 save()
46
47 p9 = pop()
48 p5 = pop()
49
50 if (istensor(p5))
51 expand_tensor()
52 restore()
53 return
54
55 # if sum of terms then sum over the expansion of each term
56
57 if (car(p5) == symbol(ADD))
58 push_integer(0)
59 p1 = cdr(p5)
60 while (iscons(p1))
61 push(car(p1))
62 push(p9)
63 expand()
64 add()
65 p1 = cdr(p1)
66 restore()
67 return
68
69 # B = numerator
70
71 push(p5)
72 numerator()
73 p3 = pop()
74
75 # A = denominator
76
77 push(p5)
78 denominator()
79 p2 = pop()
80
81 remove_negative_exponents()
82
83 # Q = quotient
84
85 push(p3)
86 push(p2)
87 push(p9)
88
89 # if the denominator is one then always bail out
90 # also bail out if the denominator is not one but
91 # it's not anything recognizable as a polynomial.
92 if isone(p3) || isone(p2)
93 if !ispoly(p2,p9) || isone(p2)
94 pop()
95 pop()
96 pop()
97 push(p5)
98 # p5 is the original input, leave unchanged
99 restore()
100 return
101
102 divpoly()
103 p7 = pop()
104
105 # remainder B = B - A * Q
106
107 push(p3)
108 push(p2)
109 push(p7)
110 multiply()
111 subtract()
112 p3 = pop()
113
114 # if the remainder is zero then we're done
115
116 if (iszero(p3))
117 push(p7)
118 restore()
119 return
120
121 # A = factor(A)
122
123 #console.log("expand - to be factored: " + p2)
124 push(p2)
125 push(p9)
126 factorpoly()
127 p2 = pop()
128 #console.log("expand - factored to: " + p2)
129
130 expand_get_C()
131 expand_get_B()
132 expand_get_A()
133
134 if (istensor(p4))
135 push(p4)
136 prev_expanding = expanding
137 expanding = 1
138 inv()
139 expanding = prev_expanding
140 push(p3)
141 inner()
142 push(p2)
143 inner()
144 else
145 push(p3)
146 push(p4)
147 prev_expanding = expanding
148 expanding = 1
149 divide()
150 expanding = prev_expanding
151 push(p2)
152 multiply()
153
154 push(p7)
155 add()
156
157 restore()
158
159expand_tensor = ->
160 i = 0
161 push(p5)
162 copy_tensor()
163 p5 = pop()
164 for i in [0...p5.tensor.nelem]
165 push(p5.tensor.elem[i])
166 push(p9)
167 expand()
168 p5.tensor.elem[i] = pop()
169 push(p5)
170
171remove_negative_exponents = ->
172 h = 0
173 i = 0
174 j = 0
175 k = 0
176 n = 0
177
178 h = tos
179 factors(p2)
180 factors(p3)
181 n = tos - h
182
183 # find the smallest exponent
184
185 j = 0
186 for i in [0...n]
187 p1 = stack[h + i]
188 if (car(p1) != symbol(POWER))
189 continue
190 if (cadr(p1) != p9)
191 continue
192 push(caddr(p1))
193 k = pop_integer()
194 if (isNaN(k))
195 continue
196 if (k < j)
197 j = k
198
199 moveTos h
200
201 if (j == 0)
202 return
203
204 # A = A / X^j
205
206 push(p2)
207 push(p9)
208 push_integer(-j)
209 power()
210 multiply()
211 p2 = pop()
212
213 # B = B / X^j
214
215 push(p3)
216 push(p9)
217 push_integer(-j)
218 power()
219 multiply()
220 p3 = pop()
221
222# Returns the expansion coefficient matrix C.
223#
224# Example:
225#
226# B 1
227# --- = -----------
228# A 2
229# x (x + 1)
230#
231# We have
232#
233# B Y1 Y2 Y3
234# --- = ---- + ---- + -------
235# A 2 x x + 1
236# x
237#
238# Our task is to solve for the unknowns Y1, Y2, and Y3.
239#
240# Multiplying both sides by A yields
241#
242# AY1 AY2 AY3
243# B = ----- + ----- + -------
244# 2 x x + 1
245# x
246#
247# Let
248#
249# A A A
250# W1 = ---- W2 = --- W3 = -------
251# 2 x x + 1
252# x
253#
254# Then the coefficient matrix C is
255#
256# coeff(W1,x,0) coeff(W2,x,0) coeff(W3,x,0)
257#
258# C = coeff(W1,x,1) coeff(W2,x,1) coeff(W3,x,1)
259#
260# coeff(W1,x,2) coeff(W2,x,2) coeff(W3,x,2)
261#
262# It follows that
263#
264# coeff(B,x,0) Y1
265#
266# coeff(B,x,1) = C Y2
267#
268# coeff(B,x,2) = Y3
269#
270# Hence
271#
272# Y1 coeff(B,x,0)
273# -1
274# Y2 = C coeff(B,x,1)
275#
276# Y3 coeff(B,x,2)
277
278expand_get_C = ->
279 h = 0
280 i = 0
281 j = 0
282 n = 0
283 #U **a
284 h = tos
285 if (car(p2) == symbol(MULTIPLY))
286 p1 = cdr(p2)
287 while (iscons(p1))
288 p5 = car(p1)
289 expand_get_CF()
290 p1 = cdr(p1)
291 else
292 p5 = p2
293 expand_get_CF()
294 n = tos - h
295 if (n == 1)
296 p4 = pop()
297 return
298 p4 = alloc_tensor(n * n)
299 p4.tensor.ndim = 2
300 p4.tensor.dim[0] = n
301 p4.tensor.dim[1] = n
302 a = h
303 for i in [0...n]
304 for j in [0...n]
305 push(stack[a+j])
306 push(p9)
307 push_integer(i)
308 power()
309 prev_expanding = expanding
310 expanding = 1
311 divide()
312 expanding = prev_expanding
313 push(p9)
314 filter()
315 p4.tensor.elem[n * i + j] = pop()
316 moveTos tos - n
317
318# The following table shows the push order for simple roots, repeated roots,
319# and inrreducible factors.
320#
321# Factor F Push 1st Push 2nd Push 3rd Push 4th
322#
323#
324# A
325# x ---
326# x
327#
328#
329# 2 A A
330# x ---- ---
331# 2 x
332# x
333#
334#
335# A
336# x + 1 -------
337# x + 1
338#
339#
340# 2 A A
341# (x + 1) ---------- -------
342# 2 x + 1
343# (x + 1)
344#
345#
346# 2 A Ax
347# x + x + 1 ------------ ------------
348# 2 2
349# x + x + 1 x + x + 1
350#
351#
352# 2 2 A Ax A Ax
353# (x + x + 1) --------------- --------------- ------------ ------------
354# 2 2 2 2 2 2
355# (x + x + 1) (x + x + 1) x + x + 1 x + x + 1
356#
357#
358# For T = A/F and F = P^N we have
359#
360#
361# Factor F Push 1st Push 2nd Push 3rd Push 4th
362#
363# x T
364#
365# 2
366# x T TP
367#
368#
369# x + 1 T
370#
371# 2
372# (x + 1) T TP
373#
374# 2
375# x + x + 1 T TX
376#
377# 2 2
378# (x + x + 1) T TX TP TPX
379#
380#
381# Hence we want to push in the order
382#
383# T * (P ^ i) * (X ^ j)
384#
385# for all i, j such that
386#
387# i = 0, 1, ..., N - 1
388#
389# j = 0, 1, ..., deg(P) - 1
390#
391# where index j runs first.
392
393expand_get_CF = ->
394 d = 0
395 i = 0
396 j = 0
397 n = 0
398
399 if (!Find(p5, p9))
400 return
401 prev_expanding = expanding
402 expanding = 1
403 trivial_divide()
404 expanding = prev_expanding
405 if (car(p5) == symbol(POWER))
406 push(caddr(p5))
407 n = pop_integer()
408 p6 = cadr(p5)
409 else
410 n = 1
411 p6 = p5
412
413 push(p6)
414 push(p9)
415 degree()
416 d = pop_integer()
417 for i in [0...n]
418 for j in [0...d]
419 push(p8)
420 push(p6)
421 push_integer(i)
422 power()
423 prev_expanding = expanding
424 expanding = 1
425 multiply()
426 expanding = prev_expanding
427 push(p9)
428 push_integer(j)
429 power()
430 prev_expanding = expanding
431 expanding = 1
432 multiply()
433 expanding = prev_expanding
434
435# Returns T = A/F where F is a factor of A.
436
437trivial_divide = ->
438 h = 0
439 if (car(p2) == symbol(MULTIPLY))
440 h = tos
441 p0 = cdr(p2)
442 while (iscons(p0))
443 if (!equal(car(p0), p5))
444 push(car(p0))
445 Eval(); # force expansion of (x+1)^2, f.e.
446 p0 = cdr(p0)
447 multiply_all(tos - h)
448 else
449 push_integer(1)
450 p8 = pop()
451
452# Returns the expansion coefficient vector B.
453
454expand_get_B = ->
455 i = 0
456 n = 0
457 if (!istensor(p4))
458 return
459 n = p4.tensor.dim[0]
460 p8 = alloc_tensor(n)
461 p8.tensor.ndim = 1
462 p8.tensor.dim[0] = n
463 for i in [0...n]
464 push(p3)
465 push(p9)
466 push_integer(i)
467 power()
468 prev_expanding = expanding
469 expanding = 1
470 divide()
471 expanding = prev_expanding
472 push(p9)
473 filter()
474 p8.tensor.elem[i] = pop()
475 p3 = p8
476
477# Returns the expansion fractions in A.
478
479expand_get_A = ->
480 h = 0
481 i = 0
482 n = 0
483 if (!istensor(p4))
484 push(p2)
485 reciprocate()
486 p2 = pop()
487 return
488 h = tos
489 if (car(p2) == symbol(MULTIPLY))
490 p8 = cdr(p2)
491 while (iscons(p8))
492 p5 = car(p8)
493 expand_get_AF()
494 p8 = cdr(p8)
495 else
496 p5 = p2
497 expand_get_AF()
498 n = tos - h
499 p8 = alloc_tensor(n)
500 p8.tensor.ndim = 1
501 p8.tensor.dim[0] = n
502 for i in [0...n]
503 p8.tensor.elem[i] = stack[h + i]
504 moveTos h
505 p2 = p8
506
507expand_get_AF = ->
508 d = 0
509 i = 0
510 j = 0
511 n = 1
512 if (!Find(p5, p9))
513 return
514 if (car(p5) == symbol(POWER))
515 push(caddr(p5))
516 n = pop_integer()
517 p5 = cadr(p5)
518 push(p5)
519 push(p9)
520 degree()
521 d = pop_integer()
522 for i in [n...0]
523 for j in [0...d]
524 push(p5)
525 push_integer(i)
526 power()
527 reciprocate()
528 push(p9)
529 push_integer(j)
530 power()
531 multiply()
532
533