1 |
|
2 |
|
3 |
|
4 |
|
5 |
|
6 |
|
7 |
|
8 |
|
9 |
|
10 |
|
11 |
|
12 |
|
13 |
|
14 |
|
15 | polycoeff = 0
|
16 | factpoly_expo = 0
|
17 |
|
18 | factorpoly = ->
|
19 | save()
|
20 |
|
21 | p2 = pop()
|
22 | p1 = pop()
|
23 |
|
24 | if (!Find(p1, p2))
|
25 | push(p1)
|
26 | restore()
|
27 | return
|
28 |
|
29 | if (!ispoly(p1, p2))
|
30 | push(p1)
|
31 | restore()
|
32 | return
|
33 |
|
34 | if (!issymbol(p2))
|
35 | push(p1)
|
36 | restore()
|
37 | return
|
38 |
|
39 | push(p1)
|
40 | push(p2)
|
41 | yyfactorpoly()
|
42 |
|
43 | restore()
|
44 |
|
45 |
|
46 |
|
47 |
|
48 |
|
49 |
|
50 |
|
51 |
|
52 |
|
53 |
|
54 |
|
55 | yyfactorpoly = ->
|
56 | h = 0
|
57 | i = 0
|
58 |
|
59 | save()
|
60 |
|
61 | p2 = pop()
|
62 | p1 = pop()
|
63 |
|
64 | h = tos
|
65 |
|
66 | if (isfloating(p1))
|
67 | stop("floating point numbers in polynomial")
|
68 |
|
69 | polycoeff = tos
|
70 |
|
71 | push(p1)
|
72 | push(p2)
|
73 | factpoly_expo = coeff() - 1
|
74 |
|
75 | rationalize_coefficients(h)
|
76 |
|
77 |
|
78 |
|
79 | whichRootsAreWeFinding = "real"
|
80 | remainingPoly = null
|
81 | while (factpoly_expo > 0)
|
82 |
|
83 | if (iszero(stack[polycoeff+0]))
|
84 | push_integer(1)
|
85 | p4 = pop()
|
86 | push_integer(0)
|
87 | p5 = pop()
|
88 | else
|
89 |
|
90 | if whichRootsAreWeFinding == "real"
|
91 | foundRealRoot = get_factor_from_real_root()
|
92 | else if whichRootsAreWeFinding == "complex"
|
93 | foundComplexRoot = get_factor_from_complex_root(remainingPoly)
|
94 |
|
95 | if whichRootsAreWeFinding == "real"
|
96 | if foundRealRoot == 0
|
97 | whichRootsAreWeFinding = "complex"
|
98 | continue
|
99 | else
|
100 |
|
101 |
|
102 | push(p4)
|
103 | push(p2)
|
104 | multiply()
|
105 | push(p5)
|
106 | add()
|
107 | p8 = pop()
|
108 |
|
109 | if (DEBUG)
|
110 | console.log("success\nFACTOR=" + p8)
|
111 |
|
112 |
|
113 |
|
114 | |
115 |
|
116 |
|
117 |
|
118 |
|
119 |
|
120 |
|
121 |
|
122 |
|
123 |
|
124 |
|
125 |
|
126 |
|
127 |
|
128 |
|
129 | push(p7)
|
130 | push(p8)
|
131 | multiply_noexpand()
|
132 | p7 = pop()
|
133 |
|
134 |
|
135 |
|
136 |
|
137 |
|
138 |
|
139 | yydivpoly()
|
140 |
|
141 | while (factpoly_expo and iszero(stack[polycoeff+factpoly_expo]))
|
142 | factpoly_expo--
|
143 |
|
144 | push(zero)
|
145 | for i in [0..factpoly_expo]
|
146 | push(stack[polycoeff+i])
|
147 | push(p2)
|
148 | push_integer(i)
|
149 | power()
|
150 | multiply()
|
151 | add()
|
152 | remainingPoly = pop()
|
153 |
|
154 |
|
155 | else if whichRootsAreWeFinding == "complex"
|
156 | if foundComplexRoot == 0
|
157 | break
|
158 | else
|
159 |
|
160 |
|
161 | push(p4)
|
162 | push(p2)
|
163 | subtract()
|
164 |
|
165 |
|
166 | push(p4)
|
167 | conjugate()
|
168 | push(p2)
|
169 | subtract()
|
170 |
|
171 |
|
172 | multiply()
|
173 |
|
174 |
|
175 |
|
176 |
|
177 |
|
178 | p8 = pop()
|
179 |
|
180 | if (DEBUG)
|
181 | console.log("success\nFACTOR=" + p8)
|
182 |
|
183 |
|
184 |
|
185 | |
186 |
|
187 |
|
188 |
|
189 |
|
190 |
|
191 |
|
192 |
|
193 |
|
194 |
|
195 |
|
196 |
|
197 |
|
198 |
|
199 |
|
200 |
|
201 | push(p7)
|
202 | previousFactorisation = pop()
|
203 |
|
204 |
|
205 |
|
206 | push(p7)
|
207 | push(p8)
|
208 | multiply_noexpand()
|
209 | p7 = pop()
|
210 |
|
211 |
|
212 |
|
213 |
|
214 |
|
215 |
|
216 |
|
217 | if !remainingPoly?
|
218 | push(zero)
|
219 | for i in [0..factpoly_expo]
|
220 | push(stack[polycoeff+i])
|
221 | push(p2)
|
222 | push_integer(i)
|
223 | power()
|
224 | multiply()
|
225 | add()
|
226 | remainingPoly = pop()
|
227 |
|
228 |
|
229 | dividend = remainingPoly
|
230 |
|
231 |
|
232 |
|
233 | push(dividend)
|
234 |
|
235 |
|
236 | push(p8)
|
237 | push(p2)
|
238 | divpoly()
|
239 | remainingPoly = pop()
|
240 |
|
241 | push(remainingPoly)
|
242 | push(p8)
|
243 | multiply()
|
244 | checkingTheDivision = pop()
|
245 |
|
246 | if !equal(checkingTheDivision, dividend)
|
247 |
|
248 |
|
249 |
|
250 |
|
251 |
|
252 |
|
253 | if DEBUG then console.log("we found a polynomial based on complex root and its conj but it doesn't divide the poly, quitting")
|
254 | if DEBUG then console.log("so just returning previousFactorisation times dividend: " + previousFactorisation + " * " + dividend)
|
255 | push(previousFactorisation)
|
256 | push(dividend)
|
257 |
|
258 | prev_expanding = expanding
|
259 | expanding = 0
|
260 | yycondense()
|
261 | expanding = prev_expanding
|
262 |
|
263 | multiply_noexpand()
|
264 | p7 = pop()
|
265 | stack[h] = p7
|
266 | moveTos h + 1
|
267 | restore()
|
268 | return
|
269 |
|
270 |
|
271 |
|
272 |
|
273 |
|
274 |
|
275 |
|
276 | |
277 |
|
278 |
|
279 |
|
280 |
|
281 |
|
282 |
|
283 |
|
284 |
|
285 |
|
286 |
|
287 | for i in [0..factpoly_expo]
|
288 | pop()
|
289 |
|
290 | push(remainingPoly)
|
291 | push(p2)
|
292 | coeff()
|
293 |
|
294 |
|
295 | factpoly_expo -= 2
|
296 |
|
297 |
|
298 |
|
299 |
|
300 |
|
301 | push(zero)
|
302 | for i in [0..factpoly_expo]
|
303 | push(stack[polycoeff+i])
|
304 | push(p2)
|
305 | push_integer(i)
|
306 | power()
|
307 | multiply()
|
308 | add()
|
309 | p1 = pop()
|
310 |
|
311 | if (DEBUG)
|
312 | console.log("POLY=" + p1)
|
313 |
|
314 | push(p1)
|
315 |
|
316 | prev_expanding = expanding
|
317 | expanding = 0
|
318 | yycondense()
|
319 | expanding = prev_expanding
|
320 |
|
321 | p1 = pop()
|
322 |
|
323 |
|
324 |
|
325 |
|
326 |
|
327 | if (factpoly_expo > 0 && isnegativeterm(stack[polycoeff+factpoly_expo]))
|
328 | push(p1)
|
329 |
|
330 |
|
331 | negate()
|
332 |
|
333 | p1 = pop()
|
334 | push(p7)
|
335 | negate_noexpand()
|
336 | p7 = pop()
|
337 |
|
338 | push(p7)
|
339 | push(p1)
|
340 | multiply_noexpand()
|
341 | p7 = pop()
|
342 |
|
343 | if (DEBUG)
|
344 | console.log("RESULT=" + p7)
|
345 |
|
346 | stack[h] = p7
|
347 |
|
348 | moveTos h + 1
|
349 |
|
350 | restore()
|
351 |
|
352 | rationalize_coefficients = (h) ->
|
353 | i = 0
|
354 |
|
355 |
|
356 |
|
357 | p7 = one
|
358 | for i in [h...tos]
|
359 | push(stack[i])
|
360 | denominator()
|
361 | push(p7)
|
362 | lcm()
|
363 | p7 = pop()
|
364 |
|
365 |
|
366 |
|
367 | for i in [h...tos]
|
368 | push(p7)
|
369 | push(stack[i])
|
370 | multiply()
|
371 | stack[i] = pop()
|
372 |
|
373 |
|
374 |
|
375 | push(p7)
|
376 | reciprocate()
|
377 | p7 = pop()
|
378 | if DEBUG then console.log("rationalize_coefficients result")
|
379 |
|
380 |
|
381 | get_factor_from_real_root = ->
|
382 |
|
383 | i = 0
|
384 | j = 0
|
385 | h = 0
|
386 | a0 = 0
|
387 | an = 0
|
388 | na0 = 0
|
389 | nan = 0
|
390 |
|
391 | if (DEBUG)
|
392 | push(zero)
|
393 | for i in [0..factpoly_expo]
|
394 | push(stack[polycoeff+i])
|
395 | push(p2)
|
396 | push_integer(i)
|
397 | power()
|
398 | multiply()
|
399 | add()
|
400 | p1 = pop()
|
401 | console.log("POLY=" + p1)
|
402 |
|
403 | h = tos
|
404 |
|
405 | an = tos
|
406 | push(stack[polycoeff+factpoly_expo])
|
407 |
|
408 | divisors_onstack()
|
409 |
|
410 | nan = tos - an
|
411 |
|
412 | a0 = tos
|
413 | push(stack[polycoeff+0])
|
414 | divisors_onstack()
|
415 | na0 = tos - a0
|
416 |
|
417 | if (DEBUG)
|
418 | console.log("divisors of base term")
|
419 | for i in [0...na0]
|
420 | console.log(", " + stack[a0 + i])
|
421 | console.log("divisors of leading term")
|
422 | for i in [0...nan]
|
423 | console.log(", " + stack[an + i])
|
424 |
|
425 |
|
426 |
|
427 | for rootsTries_i in [0...nan]
|
428 | for rootsTries_j in [0...na0]
|
429 |
|
430 |
|
431 |
|
432 | p4 = stack[an + rootsTries_i]
|
433 | p5 = stack[a0 + rootsTries_j]
|
434 |
|
435 | push(p5)
|
436 | push(p4)
|
437 | divide()
|
438 | negate()
|
439 | p3 = pop()
|
440 |
|
441 | Evalpoly()
|
442 |
|
443 | if (DEBUG)
|
444 | console.log("try A=" + p4)
|
445 | console.log(", B=" + p5)
|
446 | console.log(", root " + p2)
|
447 | console.log("=-B/A=" + p3)
|
448 | console.log(", POLY(" + p3)
|
449 | console.log(")=" + p6)
|
450 |
|
451 | if (iszero(p6))
|
452 | moveTos h
|
453 | if DEBUG then console.log "get_factor_from_real_root returning 1"
|
454 | return 1
|
455 |
|
456 | push(p5)
|
457 | negate()
|
458 | p5 = pop()
|
459 |
|
460 | push(p3)
|
461 | negate()
|
462 | p3 = pop()
|
463 |
|
464 | Evalpoly()
|
465 |
|
466 | if (DEBUG)
|
467 | console.log("try A=" + p4)
|
468 | console.log(", B=" + p5)
|
469 | console.log(", root " + p2)
|
470 | console.log("=-B/A=" + p3)
|
471 | console.log(", POLY(" + p3)
|
472 | console.log(")=" + p6)
|
473 |
|
474 | if (iszero(p6))
|
475 | moveTos h
|
476 | if DEBUG then console.log "get_factor_from_real_root returning 1"
|
477 | return 1
|
478 |
|
479 | moveTos h
|
480 |
|
481 | if DEBUG then console.log "get_factor_from_real_root returning 0"
|
482 | return 0
|
483 |
|
484 | get_factor_from_complex_root = (remainingPoly) ->
|
485 |
|
486 | i = 0
|
487 | j = 0
|
488 | h = 0
|
489 | a0 = 0
|
490 | an = 0
|
491 | na0 = 0
|
492 | nan = 0
|
493 |
|
494 | if factpoly_expo <= 2
|
495 | if DEBUG then console.log("no more factoring via complex roots to be found in polynomial of degree <= 2")
|
496 | return 0
|
497 |
|
498 | p1 = remainingPoly
|
499 | if DEBUG then console.log("complex root finding for POLY=" + p1)
|
500 |
|
501 | h = tos
|
502 | an = tos
|
503 |
|
504 |
|
505 |
|
506 | push_integer(-1)
|
507 | push_rational(2,3)
|
508 | power()
|
509 | rect()
|
510 | p4 = pop()
|
511 | if DEBUG then console.log("complex root finding: trying with " + p4)
|
512 | push(p4)
|
513 | p3 = pop()
|
514 | push(p3)
|
515 | Evalpoly()
|
516 | if DEBUG then console.log("complex root finding result: " + p6)
|
517 | if (iszero(p6))
|
518 | moveTos h
|
519 | if DEBUG then console.log "get_factor_from_complex_root returning 1"
|
520 | return 1
|
521 |
|
522 |
|
523 |
|
524 |
|
525 | push_integer(1)
|
526 | push_rational(2,3)
|
527 | power()
|
528 | rect()
|
529 | p4 = pop()
|
530 | if DEBUG then console.log("complex root finding: trying with " + p4)
|
531 | push(p4)
|
532 | p3 = pop()
|
533 | push(p3)
|
534 | Evalpoly()
|
535 | if DEBUG then console.log("complex root finding result: " + p6)
|
536 | if (iszero(p6))
|
537 | moveTos h
|
538 | if DEBUG then console.log "get_factor_from_complex_root returning 1"
|
539 | return 1
|
540 |
|
541 |
|
542 |
|
543 |
|
544 | for rootsTries_i in [-10..10]
|
545 | for rootsTries_j in [1..5]
|
546 | push_integer(rootsTries_i)
|
547 | push_integer(rootsTries_j)
|
548 | push(imaginaryunit)
|
549 | multiply()
|
550 | add()
|
551 | rect()
|
552 | p4 = pop()
|
553 |
|
554 |
|
555 | push(p4)
|
556 | p3 = pop()
|
557 |
|
558 |
|
559 | push(p3)
|
560 |
|
561 | Evalpoly()
|
562 |
|
563 |
|
564 | if (iszero(p6))
|
565 | moveTos h
|
566 | if DEBUG then console.log "found complex root: " + p6
|
567 | return 1
|
568 |
|
569 | moveTos h
|
570 |
|
571 | if DEBUG then console.log "get_factor_from_complex_root returning 0"
|
572 | return 0
|
573 |
|
574 |
|
575 |
|
576 |
|
577 |
|
578 |
|
579 |
|
580 |
|
581 |
|
582 |
|
583 |
|
584 |
|
585 |
|
586 |
|
587 |
|
588 |
|
589 |
|
590 | yydivpoly = ->
|
591 | i = 0
|
592 | p6 = zero
|
593 | for i in [factpoly_expo...0]
|
594 | push(stack[polycoeff+i])
|
595 | stack[polycoeff+i] = p6
|
596 | push(p4)
|
597 | divide()
|
598 | p6 = pop()
|
599 | push(stack[polycoeff+i - 1])
|
600 | push(p6)
|
601 | push(p5)
|
602 | multiply()
|
603 | subtract()
|
604 | stack[polycoeff+i - 1] = pop()
|
605 | stack[polycoeff+0] = p6
|
606 | if DEBUG then console.log("yydivpoly Q:")
|
607 |
|
608 |
|
609 | Evalpoly = ->
|
610 | i = 0
|
611 | push(zero)
|
612 | for i in [factpoly_expo..0]
|
613 | push(p3)
|
614 | multiply()
|
615 | push(stack[polycoeff+i])
|
616 | if DEBUG
|
617 | console.log("Evalpoly top of stack:")
|
618 | console.log print_list(stack[tos-i])
|
619 | add()
|
620 | p6 = pop()
|
621 |
|