UNPKG

3.84 kBtext/coffeescriptView Raw
1###
2 Simplify factorials
3
4The following script
5
6 F(n,k) = k binomial(n,k)
7 (F(n,k) + F(n,k-1)) / F(n+1,k)
8
9generates
10
11 k! n! n! (1 - k + n)! k! n!
12 -------------------- + -------------------- - ----------------------
13 (-1 + k)! (1 + n)! (1 + n)! (-k + n)! k (-1 + k)! (1 + n)!
14
15Simplify each term to get
16
17 k 1 - k + n 1
18 ------- + ----------- - -------
19 1 + n 1 + n 1 + n
20
21Then simplify the sum to get
22
23 n
24 -------
25 1 + n
26
27###
28
29
30
31# simplify factorials term-by-term
32
33Eval_simfac = ->
34 push(cadr(p1))
35 Eval()
36 simfac()
37
38#if 1
39
40simfac = ->
41 h = 0
42 save()
43 p1 = pop()
44 if (car(p1) == symbol(ADD))
45 h = tos
46 p1 = cdr(p1)
47 while (p1 != symbol(NIL))
48 push(car(p1))
49 simfac_term()
50 p1 = cdr(p1)
51 add_all(tos - h)
52 else
53 push(p1)
54 simfac_term()
55 restore()
56
57
58#else
59###
60void
61simfac(void)
62{
63 int h
64 save()
65 p1 = pop()
66 if (car(p1) == symbol(ADD)) {
67 h = tos
68 p1 = cdr(p1)
69 while (p1 != symbol(NIL)) {
70 push(car(p1))
71 simfac_term()
72 p1 = cdr(p1)
73 }
74 addk(tos - h)
75 p1 = pop()
76 if (find(p1, symbol(FACTORIAL))) {
77 push(p1)
78 if (car(p1) == symbol(ADD)) {
79 Condense()
80 simfac_term()
81 }
82 }
83 } else {
84 push(p1)
85 simfac_term()
86 }
87 restore()
88}
89
90#endif
91###
92
93simfac_term = ->
94 h = 0
95
96 save()
97
98 p1 = pop()
99
100 # if not a product of factors then done
101
102 if (car(p1) != symbol(MULTIPLY))
103 push(p1)
104 restore()
105 return
106
107 # push all factors
108
109 h = tos
110 p1 = cdr(p1)
111 while (p1 != symbol(NIL))
112 push(car(p1))
113 p1 = cdr(p1)
114
115 # keep trying until no more to do
116
117 while (yysimfac(h))
118 doNothing = 1
119
120 multiply_all_noexpand(tos - h)
121
122 restore()
123
124# try all pairs of factors
125
126yysimfac = (h) ->
127 i = 0
128 j = 0
129
130 for i in [h...tos]
131 p1 = stack[i]
132 for j in [h...tos]
133 if (i == j)
134 continue
135 p2 = stack[j]
136
137 # n! / n -> (n - 1)!
138
139 if (car(p1) == symbol(FACTORIAL)\
140 && car(p2) == symbol(POWER)\
141 && isminusone(caddr(p2))\
142 && equal(cadr(p1), cadr(p2)))
143 push(cadr(p1))
144 push(one)
145 subtract()
146 factorial()
147 stack[i] = pop()
148 stack[j] = one
149 return 1
150
151 # n / n! -> 1 / (n - 1)!
152
153 if (car(p2) == symbol(POWER)\
154 && isminusone(caddr(p2))\
155 && caadr(p2) == symbol(FACTORIAL)\
156 && equal(p1, cadadr(p2)))
157 push(p1)
158 push_integer(-1)
159 add()
160 factorial()
161 reciprocate()
162 stack[i] = pop()
163 stack[j] = one
164 return 1
165
166 # (n + 1) n! -> (n + 1)!
167
168 if (car(p2) == symbol(FACTORIAL))
169 push(p1)
170 push(cadr(p2))
171 subtract()
172 p3 = pop()
173 if (isplusone(p3))
174 push(p1)
175 factorial()
176 stack[i] = pop()
177 stack[j] = one
178 return 1
179
180 # 1 / ((n + 1) n!) -> 1 / (n + 1)!
181
182 if (car(p1) == symbol(POWER)\
183 && isminusone(caddr(p1))\
184 && car(p2) == symbol(POWER)\
185 && isminusone(caddr(p2))\
186 && caadr(p2) == symbol(FACTORIAL))
187 push(cadr(p1))
188 push(cadr(cadr(p2)))
189 subtract()
190 p3 = pop()
191 if (isplusone(p3))
192 push(cadr(p1))
193 factorial()
194 reciprocate()
195 stack[i] = pop()
196 stack[j] = one
197 return 1
198
199 # (n + 1)! / n! -> n + 1
200
201 # n! / (n + 1)! -> 1 / (n + 1)
202
203 if (car(p1) == symbol(FACTORIAL)\
204 && car(p2) == symbol(POWER)\
205 && isminusone(caddr(p2))\
206 && caadr(p2) == symbol(FACTORIAL))
207 push(cadr(p1))
208 push(cadr(cadr(p2)))
209 subtract()
210 p3 = pop()
211 if (isplusone(p3))
212 stack[i] = cadr(p1)
213 stack[j] = one
214 return 1
215 if (isminusone(p3))
216 push(cadr(cadr(p2)))
217 reciprocate()
218 stack[i] = pop()
219 stack[j] = one
220 return 1
221 if (equaln(p3, 2))
222 stack[i] = cadr(p1)
223 push(cadr(p1))
224 push_integer(-1)
225 add()
226 stack[j] = pop()
227 return 1
228 if (equaln(p3, -2))
229 push(cadr(cadr(p2)))
230 reciprocate()
231 stack[i] = pop()
232 push(cadr(cadr(p2)))
233 push_integer(-1)
234 add()
235 reciprocate()
236 stack[j] = pop()
237 return 1
238 return 0