UNPKG

32.3 kBtext/coffeescriptView Raw
1
2
3
4#define POLY p1
5#define X p2
6#define A p3
7#define B p4
8#define C p5
9#define Y p6
10
11show_power_debug = false
12performing_roots = false
13
14Eval_roots = ->
15 # A == B -> A - B
16
17 p2 = cadr(p1)
18
19 if (car(p2) == symbol(SETQ) || car(p2) == symbol(TESTEQ))
20 push(cadr(p2))
21 Eval()
22 push(caddr(p2))
23 Eval()
24 subtract()
25 else
26 push(p2)
27 Eval()
28 p2 = pop()
29 if (car(p2) == symbol(SETQ) || car(p2) == symbol(TESTEQ))
30 push(cadr(p2))
31 Eval()
32 push(caddr(p2))
33 Eval()
34 subtract()
35 else
36 push(p2)
37
38 # 2nd arg, x
39
40 push(caddr(p1))
41 Eval()
42 p2 = pop()
43 if (p2 == symbol(NIL))
44 guess()
45 else
46 push(p2)
47
48 p2 = pop()
49 p1 = pop()
50
51 if (!ispoly(p1, p2))
52 stop("roots: 1st argument is not a polynomial")
53
54 push(p1)
55 push(p2)
56
57 roots()
58
59
60hasImaginaryCoeff = (k) ->
61
62 #polycoeff = tos
63
64
65 imaginaryCoefficients = false
66 h = tos
67 for i in [k...0] by -1
68 #console.log "hasImaginaryCoeff - coeff.:" + stack[tos-i].toString()
69 if iscomplexnumber(stack[tos-i])
70 imaginaryCoefficients = true
71 break
72 return imaginaryCoefficients
73
74isSimpleRoot = (k) ->
75
76 #polycoeff = tos
77
78
79 #tos-n Coefficient of x^0
80 #tos-1 Coefficient of x^(n-1)
81
82 if k > 2
83 isSimpleRootPolynomial = true
84 h = tos
85
86 if iszero(stack[tos-k])
87 isSimpleRootPolynomial = false
88
89 for i in [(k-1)...1] by -1
90 #console.log "hasImaginaryCoeff - coeff.:" + stack[tos-i].toString()
91 if !iszero(stack[tos-i])
92 isSimpleRootPolynomial = false
93 break
94 else
95 isSimpleRootPolynomial = false
96
97 return isSimpleRootPolynomial
98
99
100normalisedCoeff = ->
101 k = coeff()
102 #console.log("->" + tos)
103 divideBy = stack[tos-1]
104 miniStack = []
105 for i in [1..k]
106 miniStack.push(pop())
107 #console.log(tos)
108
109 for i in [(k-1)..0]
110 push(miniStack[i])
111 push(divideBy)
112 divide()
113 #console.log(tos)
114 return k
115
116
117
118# takes the polynomial and the
119# variable on the stack
120
121roots = ->
122 h = 0
123 i = 0
124 n = 0
125
126 save()
127
128 # the simplification of nested radicals uses
129 # "roots", which in turn uses simplification
130 # of nested radicals. Usually there is no problem,
131 # one level of recursion does the job. Beyond that,
132 # we probably got stuck in a strange case of infinite
133 # recursion, so bail out and return NIL.
134 if recursionLevelNestedRadicalsRemoval > 1
135 pop()
136 pop()
137 push(symbol(NIL))
138 restore()
139 return
140
141 performing_roots = true
142 h = tos - 2
143
144 if DEBUG then console.log("checking if " + stack[tos-1].toString() + " is a case of simple roots")
145
146 p2 = pop()
147 p1 = pop()
148
149 push(p1)
150 push(p2)
151
152 push(p1)
153 push(p2)
154 k = normalisedCoeff()
155
156 if isSimpleRoot(k)
157 if DEBUG then console.log("yes, " + stack[tos-1].toString() + " is a case of simple roots")
158 lastCoeff = stack[tos-k]
159 leadingCoeff = stack[tos-1]
160 moveTos tos - k
161 pop()
162 pop()
163 getSimpleRoots(k, leadingCoeff, lastCoeff)
164 else
165 moveTos tos - k
166 roots2()
167
168 n = tos - h
169 if (n == 0)
170 stop("roots: the polynomial is not factorable, try nroots")
171 if (n == 1)
172 performing_roots = false
173 restore()
174 return
175 sort_stack(n)
176 p1 = alloc_tensor(n)
177 p1.tensor.ndim = 1
178 p1.tensor.dim[0] = n
179 for i in [0...n]
180 p1.tensor.elem[i] = stack[h + i]
181 moveTos h
182 push(p1)
183 restore()
184 performing_roots = false
185
186# ok to generate these roots take a look at their form
187# in the case of even and odd exponents here:
188# http://www.wolframalpha.com/input/?i=roots+x%5E14+%2B+1
189# http://www.wolframalpha.com/input/?i=roots+ax%5E14+%2B+b
190# http://www.wolframalpha.com/input/?i=roots+x%5E15+%2B+1
191# http://www.wolframalpha.com/input/?i=roots+a*x%5E15+%2B+b
192getSimpleRoots = (n, leadingCoeff, lastCoeff) ->
193 if DEBUG then console.log("getSimpleRoots")
194 save()
195
196 #tos-n Coefficient of x^0
197 #tos-1 Coefficient of x^(n-1)
198
199 n = n-1
200
201 push(lastCoeff)
202 push_rational(1,n)
203 power()
204
205 push(leadingCoeff)
206 push_rational(1,n)
207 power()
208 divide()
209
210 commonPart = pop()
211
212 if n % 2 == 0
213 for rootsOfOne in [1..n] by 2
214 push(commonPart)
215 push_integer(-1)
216 push_rational(rootsOfOne,n)
217 power()
218 multiply()
219 aSol = pop()
220 push(aSol)
221 push(aSol)
222 negate()
223 else
224 for rootsOfOne in [1..n]
225 push(commonPart)
226 push_integer(-1)
227 push_rational(rootsOfOne,n)
228 power()
229 multiply()
230 if rootsOfOne % 2 == 0
231 negate()
232
233
234 restore()
235
236roots2 = ->
237 save()
238
239 p2 = pop() # the polynomial variable
240 p1 = pop() # the polynomial
241
242 push(p1)
243 push(p2)
244
245 push(p1)
246 push(p2)
247 k = normalisedCoeff()
248
249 if !hasImaginaryCoeff(k)
250 moveTos tos - k
251 factorpoly()
252 p1 = pop()
253 else
254 moveTos tos - k
255 pop()
256 pop()
257
258
259 if (car(p1) == symbol(MULTIPLY))
260 p1 = cdr(p1)
261 # scan through all the factors
262 # and find the roots of each of them
263 while (iscons(p1))
264 push(car(p1))
265 push(p2)
266 roots3()
267 p1 = cdr(p1)
268 else
269 push(p1)
270 push(p2)
271 roots3()
272
273 restore()
274
275roots3 = ->
276 save()
277 p2 = pop()
278 p1 = pop()
279 if (car(p1) == symbol(POWER) && ispoly(cadr(p1), p2) && isposint(caddr(p1)))
280 push(cadr(p1))
281 push(p2)
282 n = normalisedCoeff()
283 mini_solve(n)
284 else if (ispoly(p1, p2))
285 push(p1)
286 push(p2)
287 n = normalisedCoeff()
288 mini_solve(n)
289 restore()
290
291#-----------------------------------------------------------------------------
292#
293# Input: stack[tos - 2] polynomial
294#
295# stack[tos - 1] dependent symbol
296#
297# Output: stack roots on stack
298#
299# (input args are popped first)
300#
301#-----------------------------------------------------------------------------
302
303# note that for many quadratic, cubic and quartic polynomials we don't
304# actually end up using the quadratic/cubic/quartic formulas in here,
305# since there is a chance we factored the polynomial and in so
306# doing we found some solutions and lowered the degree.
307mini_solve = (n) ->
308 #console.log "mini_solve >>>>>>>>>>>>>>>>>>>>>>>> tos:" + tos
309
310 save()
311
312 # AX + B, X = -B/A
313
314 if (n == 2)
315 #console.log "mini_solve >>>>>>>>> 1st degree"
316 p3 = pop()
317 p4 = pop()
318 push(p4)
319 push(p3)
320 divide()
321 negate()
322 restore()
323 return
324
325 # AX^2 + BX + C, X = (-B +/- (B^2 - 4AC)^(1/2)) / (2A)
326
327 if (n == 3)
328 #console.log "mini_solve >>>>>>>>> 2nd degree"
329 p3 = pop() # A
330 p4 = pop() # B
331 p5 = pop() # C
332
333 # B^2
334 push(p4)
335 push_integer(2)
336 power()
337
338 # 4AC
339 push_integer(4)
340 push(p3)
341 multiply()
342 push(p5)
343 multiply()
344
345 # B^2 - 4AC
346 subtract()
347
348 #(B^2 - 4AC)^(1/2)
349 push_rational(1, 2)
350 power()
351
352 #p6 is (B^2 - 4AC)^(1/2)
353 p6 = pop()
354 push(p6);
355
356 # B
357 push(p4)
358 subtract() # -B + (B^2 - 4AC)^(1/2)
359
360 # 1/2A
361 push(p3)
362 push_integer(2)
363 multiply()
364 divide()
365 #simplify()
366 #rationalize()
367 # tos - 1 now is 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
368
369 push(p6);
370 # tos - 1 now is (B^2 - 4AC)^(1/2)
371 # tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
372
373 # add B to tos
374 push(p4)
375 add()
376 # tos - 1 now is B + (B^2 - 4AC)^(1/2)
377 # tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
378
379 negate()
380 # tos - 1 now is -B -(B^2 - 4AC)^(1/2)
381 # tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
382
383 # 1/2A again
384 push(p3)
385 divide()
386 push_rational(1, 2)
387 multiply()
388 #simplify()
389 #rationalize()
390 # tos - 1: 2nd root: (-B - (B^2 - 4AC)^(1/2)) / (2A)
391 # tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
392
393 restore()
394 return
395
396 #if (n == 4)
397 if (n == 4 or n == 5)
398
399 p3 = pop() # A
400 p4 = pop() # B
401 p5 = pop() # C
402 p6 = pop() # D
403
404 # C - only related calculations
405 push(p5)
406 push(p5)
407 multiply()
408 R_c2 = pop()
409
410 push(R_c2)
411 push(p5)
412 multiply()
413 R_c3 = pop()
414
415 # B - only related calculations
416 push(p4)
417 push(p4)
418 multiply()
419 R_b2 = pop()
420
421 push(R_b2)
422 push(p4)
423 multiply()
424 R_b3 = pop()
425
426 push(R_b3)
427 push(p6)
428 multiply()
429 R_b3_d = pop()
430
431 push(R_b3_d)
432 push_integer(-4)
433 multiply()
434 R_m4_b3_d = pop()
435
436
437 push(R_b3)
438 push_integer(2)
439 multiply()
440 R_2_b3 = pop()
441
442 # A - only related calculations
443 push(p3)
444 push(p3)
445 multiply()
446 R_a2 = pop()
447
448 push(R_a2)
449 push(p3)
450 multiply()
451 R_a3 = pop()
452
453 push_integer(3)
454 push(p3)
455 multiply()
456 R_3_a = pop()
457
458 push(R_a2)
459 push(p6)
460 multiply()
461 R_a2_d = pop()
462
463 push(R_a2_d)
464 push(p6)
465 multiply()
466 R_a2_d2 = pop()
467
468 push(R_a2_d)
469 push_integer(27)
470 multiply()
471 R_27_a2_d = pop()
472
473 push(R_a2_d2)
474 push_integer(-27)
475 multiply()
476 R_m27_a2_d2 = pop()
477
478 push(R_3_a)
479 push_integer(2)
480 multiply()
481 R_6_a = pop()
482
483 # mixed calculations
484 push(p3)
485 push(p5)
486 multiply()
487 R_a_c = pop()
488
489 push(R_a_c)
490 push(p4)
491 multiply()
492 R_a_b_c = pop()
493
494 push(R_a_b_c)
495 push(p6)
496 multiply()
497 R_a_b_c_d = pop()
498
499 push(R_a_c)
500 push_integer(3)
501 multiply()
502 R_3_a_c = pop()
503
504 push_integer(-4)
505 push(p3)
506 push(R_c3)
507 multiply()
508 multiply()
509 R_m4_a_c3 = pop()
510
511 push(R_a_b_c)
512 push_integer(9)
513 multiply()
514 negate()
515 R_m9_a_b_c = pop()
516
517 push(R_a_b_c_d)
518 push_integer(18)
519 multiply()
520 R_18_a_b_c_d = pop()
521
522 push(R_b2)
523 push(R_3_a_c)
524 subtract()
525 R_DELTA0 = pop()
526
527 push(R_b2)
528 push(R_c2)
529 multiply()
530 R_b2_c2 = pop()
531
532 push(p4)
533 negate()
534 push(R_3_a)
535 divide()
536 R_m_b_over_3a = pop()
537
538
539 if (n == 4)
540 if DEBUG then console.log ">>>>>>>>>>>>>>>> actually using cubic formula <<<<<<<<<<<<<<< "
541
542 #console.log ">>>> A:" + p3.toString()
543 #console.log ">>>> B:" + p4.toString()
544 #console.log ">>>> C:" + p5.toString()
545 #console.log ">>>> D:" + p6.toString()
546
547
548 if DEBUG then console.log "cubic: D0: " + R_DELTA0.toString()
549
550 push(R_DELTA0)
551 push_integer(3)
552 power()
553 push_integer(4)
554 multiply()
555 R_4_DELTA03 = pop()
556
557 push(R_DELTA0)
558 simplify()
559 absValFloat()
560 R_DELTA0_toBeCheckedIfZero = pop()
561 if DEBUG then console.log "cubic: D0 as float: " + R_DELTA0_toBeCheckedIfZero.toString()
562 #if iszero(R_DELTA0_toBeCheckedIfZero)
563 # console.log " *********************************** D0 IS ZERO"
564
565
566 # DETERMINANT
567 push(R_18_a_b_c_d)
568 push(R_m4_b3_d)
569 push(R_b2_c2)
570 push(R_m4_a_c3)
571 push(R_m27_a2_d2)
572 add()
573 add()
574 add()
575 add()
576 simplify()
577 absValFloat()
578 R_determinant = pop()
579 if DEBUG then console.log "cubic: DETERMINANT: " + R_determinant.toString()
580
581 # R_DELTA1
582 push(R_2_b3)
583 push(R_m9_a_b_c)
584 push(R_27_a2_d)
585 add()
586 add()
587 R_DELTA1 = pop()
588 if DEBUG then console.log "cubic: D1: " + R_DELTA1.toString()
589
590 # R_Q
591 push(R_DELTA1)
592 push_integer(2)
593 power()
594 push(R_4_DELTA03)
595 subtract()
596 push_rational(1, 2)
597 power()
598 simplify()
599 R_Q = pop()
600
601
602 if iszero(R_determinant)
603 if iszero(R_DELTA0_toBeCheckedIfZero)
604 if DEBUG then console.log " cubic: DETERMINANT IS ZERO and delta0 is zero"
605 push(R_m_b_over_3a) # just same solution three times
606 restore()
607 return
608 else
609 if DEBUG then console.log " cubic: DETERMINANT IS ZERO and delta0 is not zero"
610 push(p3)
611 push(p6)
612 push_integer(9)
613 multiply()
614 multiply()
615 push(p4)
616 push(p5)
617 multiply()
618 subtract()
619 push(R_DELTA0)
620 push_integer(2)
621 multiply()
622 divide() # first solution
623 root_solution = pop()
624 push(root_solution) # pushing two of them on the stack
625 push(root_solution)
626
627 # second solution here
628 # 4abc
629 push(R_a_b_c)
630 push_integer(4)
631 multiply()
632
633 # -9a*a*d
634 push(p3)
635 push(p3)
636 push(p6)
637 push_integer(9)
638 multiply()
639 multiply()
640 multiply()
641 negate()
642
643 # -9*b^3
644 push(R_b3)
645 negate()
646
647 # sum the three terms
648 add()
649 add()
650
651 # denominator is a*delta0
652 push(p3)
653 push(R_DELTA0)
654 multiply()
655
656 # build the fraction
657 divide()
658
659 restore()
660 return
661
662
663
664 C_CHECKED_AS_NOT_ZERO = false
665 flipSignOFQSoCIsNotZero = false
666
667 # C will go as denominator, we have to check
668 # that is not zero
669 while !C_CHECKED_AS_NOT_ZERO
670
671 # R_C
672 push(R_Q)
673 if flipSignOFQSoCIsNotZero
674 negate()
675 push(R_DELTA1)
676 add()
677 push_rational(1, 2)
678 multiply()
679 push_rational(1, 3)
680 power()
681 simplify()
682 R_C = pop()
683 if DEBUG then console.log "cubic: C: " + R_C.toString()
684
685 push(R_C)
686 simplify()
687 absValFloat()
688 R_C_simplified_toCheckIfZero = pop()
689 if DEBUG then console.log "cubic: C as absval and float: " + R_C_simplified_toCheckIfZero.toString()
690
691 if iszero(R_C_simplified_toCheckIfZero)
692 if DEBUG then console.log " cubic: C IS ZERO flipping the sign"
693 flipSignOFQSoCIsNotZero = true
694 else
695 C_CHECKED_AS_NOT_ZERO = true
696
697
698 push(R_C)
699 push(R_3_a)
700 multiply()
701 R_3_a_C = pop()
702
703 push(R_3_a_C)
704 push_integer(2)
705 multiply()
706 R_6_a_C = pop()
707
708
709 # imaginary parts calculations
710 push(imaginaryunit)
711 push_integer(3)
712 push_rational(1, 2)
713 power()
714 multiply()
715 i_sqrt3 = pop()
716 push_integer(1)
717 push(i_sqrt3)
718 add()
719 one_plus_i_sqrt3 = pop()
720 push_integer(1)
721 push(i_sqrt3)
722 subtract()
723 one_minus_i_sqrt3 = pop()
724
725
726 push(R_C)
727 push(R_3_a)
728 divide()
729 R_C_over_3a = pop()
730
731 # first solution
732 push(R_m_b_over_3a) # first term
733 push(R_C_over_3a)
734 negate() # second term
735 push(R_DELTA0)
736 push(R_3_a_C)
737 divide()
738 negate() # third term
739 # now add the three terms together
740 add()
741 add()
742 simplify()
743
744 # second solution
745 push(R_m_b_over_3a) # first term
746 push(R_C_over_3a)
747 push(one_plus_i_sqrt3)
748 multiply()
749 push_integer(2)
750 divide() # second term
751 push(one_minus_i_sqrt3)
752 push(R_DELTA0)
753 multiply()
754 push(R_6_a_C)
755 divide() # third term
756 # now add the three terms together
757 add()
758 add()
759 simplify()
760
761 # third solution
762 push(R_m_b_over_3a) # first term
763 push(R_C_over_3a)
764 push(one_minus_i_sqrt3)
765 multiply()
766 push_integer(2)
767 divide() # second term
768 push(one_plus_i_sqrt3)
769 push(R_DELTA0)
770 multiply()
771 push(R_6_a_C)
772 divide() # third term
773 # now add the three terms together
774 add()
775 add()
776 simplify()
777
778 restore()
779 return
780
781 # See http://www.sscc.edu/home/jdavidso/Math/Catalog/Polynomials/Fourth/Fourth.html
782 # for a description of general shapes and properties of fourth degree polynomials
783 if (n == 5)
784 if DEBUG then console.log ">>>>>>>>>>>>>>>> actually using quartic formula <<<<<<<<<<<<<<< "
785 p7 = pop() # E
786
787 if iszero(p4) and iszero(p6) and !iszero(p5) and !iszero(p7)
788 if DEBUG then console.log("biquadratic case")
789 push(p3)
790 push(symbol(SECRETX))
791 push_integer(2)
792 power()
793 multiply()
794
795 push(p5)
796 push(symbol(SECRETX))
797 multiply()
798
799 push(p7)
800
801 add()
802 add()
803
804 push(symbol(SECRETX))
805 roots()
806
807 biquadraticSolutions = pop()
808
809 for eachSolution in biquadraticSolutions.tensor.elem
810 push(eachSolution)
811 push_rational(1,2)
812 power()
813 simplify()
814
815 push(eachSolution)
816 push_rational(1,2)
817 power()
818 negate()
819 simplify()
820
821 restore()
822 return
823
824
825
826
827
828 # D - only related calculations
829 push(p6)
830 push(p6)
831 multiply()
832 R_d2 = pop()
833
834 # E - only related calculations
835 push(p7)
836 push(p7)
837 multiply()
838 R_e2 = pop()
839
840 push(R_e2)
841 push(p7)
842 multiply()
843 R_e3 = pop()
844
845 # DETERMINANT
846
847 push_integer(256)
848 push(R_a3)
849 push(R_e3)
850 multiply()
851 multiply() # first term 256 a^3 e^3
852
853 push_integer(-192)
854 push(R_a2_d)
855 push(R_e2)
856 push(p4)
857 multiply()
858 multiply()
859 multiply() # second term -192 a^3 b d e^2
860
861 push_integer(-128)
862 push(R_a2)
863 push(R_c2)
864 push(R_e2)
865 multiply()
866 multiply()
867 multiply() # third term -128 a^2 c^2 e^2
868
869 push_integer(144)
870 push(R_a2_d2)
871 push(p5)
872 push(p7)
873 multiply()
874 multiply()
875 multiply() # fourth term 144 a^2 c d^2 e
876
877 push(R_m27_a2_d2)
878 push(R_d2)
879 multiply() # fifth term -27 a^2 d^4
880
881 push_integer(144)
882 push(R_a_b_c)
883 push(p4)
884 push(R_e2)
885 multiply()
886 multiply()
887 multiply() # sixth term 144 a b^2 c e^2
888
889 push_integer(-6)
890 push(p3)
891 push(R_b2)
892 push(R_d2)
893 push(p7)
894 multiply()
895 multiply()
896 multiply()
897 multiply() # seventh term -6 a b^2 d^2 e
898
899 push_integer(-80)
900 push(R_a_b_c_d)
901 push(p5)
902 push(p7)
903 multiply()
904 multiply()
905 multiply() # eigth term -80 a b c^2 d e
906
907 push_integer(18)
908 push(R_a_b_c_d)
909 push(R_d2)
910 multiply()
911 multiply() # ninth term 18 a b c d^3
912
913 push_integer(16)
914 push(R_a_c)
915 push(R_c3)
916 push(p7)
917 multiply()
918 multiply()
919 multiply() # tenth term 16 a c^4 e
920
921 push_integer(-4)
922 push(R_a_c)
923 push(R_c2)
924 push(R_d2)
925 multiply()
926 multiply()
927 multiply() # eleventh term -4 a c^3 d^2
928
929 push_integer(-27)
930 push(R_b3)
931 push(p4)
932 push(R_e2)
933 multiply()
934 multiply()
935 multiply() # twelveth term -27 b^4 e^2
936
937 push_integer(18)
938 push(R_b3_d)
939 push(p5)
940 push(p7)
941 multiply()
942 multiply()
943 multiply() # thirteenth term 18 b^3 c d e
944
945 push(R_m4_b3_d)
946 push(R_d2)
947 multiply() # fourteenth term -4 b^3 d^3
948
949 push_integer(-4)
950 push(R_b2_c2)
951 push(p5)
952 push(p7)
953 multiply()
954 multiply()
955 multiply() # fifteenth term -4 b^2 c^3 e
956
957 push(R_b2_c2)
958 push(R_d2)
959 multiply() # sixteenth term b^2 c^2 d^2
960
961 # add together the sixteen terms by doing
962 # fifteen adds
963 add()
964 add()
965 add()
966 add()
967 add()
968 add()
969 add()
970 add()
971 add()
972 add()
973 add()
974 add()
975 add()
976 add()
977 add()
978
979 R_determinant = pop()
980 if DEBUG then console.log("R_determinant: " + R_determinant.toString())
981
982 # DELTA0
983 push(R_c2) # term one of DELTA0
984
985 push_integer(-3)
986 push(p4)
987 push(p6)
988 multiply()
989 multiply() # term two of DELTA0
990
991 push_integer(12)
992 push(p3)
993 push(p7)
994 multiply()
995 multiply() # term three of DELTA0
996
997 # add the three terms together
998 add()
999 add()
1000
1001 R_DELTA0 = pop()
1002 if DEBUG then console.log("R_DELTA0: " + R_DELTA0.toString())
1003
1004 # DELTA1
1005 push_integer(2)
1006 push(R_c3)
1007 multiply()
1008
1009 push_integer(-9)
1010 push(p4)
1011 push(p5)
1012 push(p6)
1013 multiply()
1014 multiply()
1015 multiply()
1016
1017 push_integer(27)
1018 push(R_b2)
1019 push(p7)
1020 multiply()
1021 multiply()
1022
1023 push_integer(27)
1024 push(p3)
1025 push(R_d2)
1026 multiply()
1027 multiply()
1028
1029 push_integer(-72)
1030 push(R_a_c)
1031 push(p7)
1032 multiply()
1033 multiply()
1034
1035 # add the five terms together
1036 add()
1037 add()
1038 add()
1039 add()
1040
1041 R_DELTA1 = pop()
1042 if DEBUG then console.log("R_DELTA1: " + R_DELTA1.toString())
1043
1044 # p
1045 push_integer(8)
1046 push(R_a_c)
1047 multiply()
1048
1049 push_integer(-3)
1050 push(R_b2)
1051 multiply()
1052
1053 add()
1054
1055 push_integer(8)
1056 push(R_a2)
1057 multiply()
1058
1059 divide()
1060
1061 R_p = pop()
1062 if DEBUG then console.log("p: " + R_p.toString())
1063
1064 # q
1065 push(R_b3)
1066
1067 push_integer(-4)
1068 push(R_a_b_c)
1069 multiply()
1070
1071 push_integer(8)
1072 push(R_a2_d)
1073 multiply()
1074
1075 add()
1076 add()
1077
1078 push_integer(8)
1079 push(R_a3)
1080 multiply()
1081
1082 divide()
1083
1084 R_q = pop()
1085 if DEBUG then console.log("q: " + R_q.toString())
1086
1087 if DEBUG then console.log("tos 1 " + tos)
1088 if !iszero(p4)
1089 if DEBUG then console.log("tos 2 " + tos)
1090
1091 push_integer(8)
1092 push(p5)
1093 push(p3)
1094 multiply()
1095 multiply()
1096
1097 push_integer(-3)
1098 push(p4)
1099 push_integer(2)
1100 power()
1101 multiply()
1102
1103 add()
1104
1105 push_integer(8)
1106 push(p3)
1107 push_integer(2)
1108 power()
1109 multiply()
1110
1111 divide()
1112
1113 R_p = pop()
1114 if DEBUG then console.log("p for depressed quartic: " + R_p.toString())
1115
1116 push(p4)
1117 push_integer(3)
1118 power()
1119
1120 push_integer(-4)
1121 push(p3)
1122 push(p4)
1123 push(p5)
1124 multiply()
1125 multiply()
1126 multiply()
1127
1128 push_integer(8)
1129 push(p6)
1130 push(p3)
1131 push_integer(2)
1132 power()
1133 multiply()
1134 multiply()
1135
1136 add()
1137 add()
1138
1139
1140 push_integer(8)
1141 push(p3)
1142 push_integer(3)
1143 power()
1144 multiply()
1145
1146 divide()
1147
1148 R_q = pop()
1149 if DEBUG then console.log("q for depressed quartic: " + R_q.toString())
1150
1151
1152 # convert to depressed quartic
1153 push(p4)
1154 push_integer(4)
1155 power()
1156 push_integer(-3)
1157 multiply()
1158
1159 push_integer(256)
1160 push(R_a3)
1161 push(p7)
1162 multiply()
1163 multiply()
1164
1165 push_integer(-64)
1166 push(R_a2_d)
1167 push(p4)
1168 multiply()
1169 multiply()
1170
1171 push_integer(16)
1172 push(R_b2)
1173 push(p3)
1174 push(p5)
1175 multiply()
1176 multiply()
1177 multiply()
1178
1179 add()
1180 add()
1181 add()
1182
1183 push_integer(256)
1184 push(p3)
1185 push_integer(4)
1186 power()
1187 multiply()
1188 divide()
1189
1190 R_r = pop()
1191 if DEBUG then console.log("r for depressed quartic: " + R_r.toString())
1192
1193 if DEBUG then console.log("tos 4 " + tos)
1194
1195 push(symbol(SECRETX))
1196 push_integer(4)
1197 power()
1198 if DEBUG then console.log("4 * x^4: " + stack[tos-1].toString())
1199
1200 push(R_p)
1201 push(symbol(SECRETX))
1202 push_integer(2)
1203 power()
1204 multiply()
1205 if DEBUG then console.log("R_p * x^2: " + stack[tos-1].toString())
1206
1207 push(R_q)
1208 push(symbol(SECRETX))
1209 multiply()
1210 if DEBUG then console.log("R_q * x: " + stack[tos-1].toString())
1211
1212 push(R_r)
1213 if DEBUG then console.log("R_r: " + stack[tos-1].toString())
1214
1215 add()
1216 add()
1217 add()
1218
1219 simplify()
1220 if DEBUG then console.log("solving depressed quartic: " + stack[tos-1].toString())
1221
1222 push(symbol(SECRETX))
1223
1224 roots()
1225
1226 depressedSolutions = pop()
1227 if DEBUG then console.log("depressedSolutions: " + depressedSolutions)
1228
1229 for eachSolution in depressedSolutions.tensor.elem
1230 push(eachSolution)
1231 push(p4)
1232 push_integer(4)
1233 push(p3)
1234 multiply()
1235 divide()
1236 subtract()
1237 simplify()
1238 if DEBUG then console.log("solution from depressed: " + stack[tos-1].toString())
1239
1240 restore()
1241 return
1242
1243 else
1244 R_p = p5
1245 R_q = p6
1246 R_r = p7
1247
1248 ###
1249 # Descartes' solution
1250 # https://en.wikipedia.org/wiki/Quartic_function#Descartes.27_solution
1251 # finding the "u" in the depressed equation
1252
1253 push_integer(2)
1254 push(R_p)
1255 multiply()
1256 coeff2 = pop()
1257
1258 push_integer(-4)
1259 push(R_p)
1260 push_integer(2)
1261 power()
1262 multiply()
1263 push(R_r)
1264 multiply()
1265 coeff3 = pop()
1266
1267 push(R_q)
1268 push_integer(2)
1269 power()
1270 negate()
1271 coeff4 = pop()
1272
1273 # now build the polynomial
1274 push(symbol(SECRETX))
1275 push_integer(3)
1276 power()
1277
1278 push(coeff2)
1279 push(symbol(SECRETX))
1280 push_integer(2)
1281 power()
1282 multiply()
1283
1284 push(coeff3)
1285 push(symbol(SECRETX))
1286 multiply()
1287
1288 push(coeff4)
1289
1290 add()
1291 add()
1292 add()
1293
1294 console.log("Descarte's resolventCubic: " + stack[tos-1].toString())
1295 push(symbol(SECRETX))
1296
1297 roots()
1298
1299 resolventCubicSolutions = pop()
1300 console.log("Descarte's resolventCubic solutions: " + resolventCubicSolutions)
1301 console.log("tos: " + tos)
1302
1303 R_u = null
1304 #R_u = resolventCubicSolutions.tensor.elem[1]
1305 for eachSolution in resolventCubicSolutions.tensor.elem
1306 console.log("examining solution: " + eachSolution)
1307 push(eachSolution)
1308 push_integer(2)
1309 multiply()
1310 push(R_p)
1311 add()
1312
1313 absValFloat()
1314 toBeCheckedIFZero = pop()
1315 console.log("abs value is: " + eachSolution)
1316 if !iszero(toBeCheckedIFZero)
1317 R_u = eachSolution
1318 break
1319
1320 console.log("chosen solution: " + R_u)
1321
1322 push(R_u)
1323 negate()
1324 R_s = pop()
1325
1326 push(R_p)
1327 push(R_u)
1328 push_integer(2)
1329 power()
1330 push(R_q)
1331 push(R_u)
1332 divide()
1333 add()
1334 add()
1335 push_integer(2)
1336 divide()
1337 R_t = pop()
1338
1339 push(R_p)
1340 push(R_u)
1341 push_integer(2)
1342 power()
1343 push(R_q)
1344 push(R_u)
1345 divide()
1346 subtract()
1347 add()
1348 push_integer(2)
1349 divide()
1350 R_v = pop()
1351
1352 # factoring the quartic into two quadratics:
1353
1354 # now build the polynomial
1355 push(symbol(SECRETX))
1356 push_integer(2)
1357 power()
1358
1359 push(R_s)
1360 push(symbol(SECRETX))
1361 multiply()
1362
1363 push(R_t)
1364
1365 add()
1366 add()
1367
1368 console.log("factored quartic 1: " + stack[tos-1].toString())
1369
1370 push(symbol(SECRETX))
1371 push_integer(2)
1372 power()
1373
1374 push(R_u)
1375 push(symbol(SECRETX))
1376 multiply()
1377
1378 push(R_v)
1379
1380 add()
1381 add()
1382
1383 console.log("factored quartic 2: " + stack[tos-1].toString())
1384 pop()
1385
1386 restore()
1387 return
1388 ###
1389
1390 # Ferrari's solution
1391 # https://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution
1392 # finding the "m" in the depressed equation
1393 push_rational(5,2)
1394 push(R_p)
1395 multiply()
1396 coeff2 = pop()
1397
1398 push_integer(2)
1399 push(R_p)
1400 push_integer(2)
1401 power()
1402 multiply()
1403 push(R_r)
1404 subtract()
1405 coeff3 = pop()
1406
1407 push(R_p)
1408 push_integer(3)
1409 power()
1410 push_integer(2)
1411 divide()
1412
1413 push_rational(-1,2)
1414 push(R_p)
1415 push(R_r)
1416 multiply()
1417 multiply()
1418
1419 push_rational(-1,8)
1420 push(R_q)
1421 push_integer(2)
1422 power()
1423 multiply()
1424
1425 add()
1426 add()
1427
1428 coeff4 = pop()
1429
1430 push(symbol(SECRETX))
1431 push_integer(3)
1432 power()
1433
1434 push(coeff2)
1435 push(symbol(SECRETX))
1436 push_integer(2)
1437 power()
1438 multiply()
1439
1440 push(coeff3)
1441 push(symbol(SECRETX))
1442 multiply()
1443
1444 push(coeff4)
1445
1446 add()
1447 add()
1448 add()
1449
1450 if DEBUG then console.log("resolventCubic: " + stack[tos-1].toString())
1451 push(symbol(SECRETX))
1452
1453 roots()
1454
1455 resolventCubicSolutions = pop()
1456 if DEBUG then console.log("resolventCubicSolutions: " + resolventCubicSolutions)
1457
1458 R_m = null
1459 #R_m = resolventCubicSolutions.tensor.elem[1]
1460 for eachSolution in resolventCubicSolutions.tensor.elem
1461 if DEBUG then console.log("examining solution: " + eachSolution)
1462 push(eachSolution)
1463 push_integer(2)
1464 multiply()
1465 push(R_p)
1466 add()
1467
1468 absValFloat()
1469 toBeCheckedIFZero = pop()
1470 if DEBUG then console.log("abs value is: " + eachSolution)
1471 if !iszero(toBeCheckedIFZero)
1472 R_m = eachSolution
1473 break
1474
1475 if DEBUG then console.log("chosen solution: " + R_m)
1476 push(R_m)
1477 push_integer(2)
1478 multiply()
1479 push(R_p)
1480 add()
1481 push_rational(1,2)
1482 power()
1483 simplify()
1484 sqrtPPlus2M = pop()
1485
1486 push(R_q)
1487 push_integer(2)
1488 multiply()
1489 push(sqrtPPlus2M)
1490 divide()
1491 simplify()
1492 TwoQOversqrtPPlus2M = pop()
1493
1494 push(R_p)
1495 push_integer(3)
1496 multiply()
1497 push(R_m)
1498 push_integer(2)
1499 multiply()
1500 add()
1501 ThreePPlus2M = pop()
1502
1503 # solution1
1504 push(sqrtPPlus2M)
1505 push(ThreePPlus2M)
1506 push(TwoQOversqrtPPlus2M)
1507 add()
1508 negate()
1509 push_rational(1,2)
1510 power()
1511 simplify()
1512 add()
1513 push_integer(2)
1514 divide()
1515 # solution2
1516 push(sqrtPPlus2M)
1517 push(ThreePPlus2M)
1518 push(TwoQOversqrtPPlus2M)
1519 add()
1520 negate()
1521 push_rational(1,2)
1522 power()
1523 simplify()
1524 subtract()
1525 push_integer(2)
1526 divide()
1527 # solution3
1528 push(sqrtPPlus2M)
1529 negate()
1530 push(ThreePPlus2M)
1531 push(TwoQOversqrtPPlus2M)
1532 subtract()
1533 negate()
1534 push_rational(1,2)
1535 power()
1536 simplify()
1537 add()
1538 push_integer(2)
1539 divide()
1540 # solution4
1541 push(sqrtPPlus2M)
1542 negate()
1543 push(ThreePPlus2M)
1544 push(TwoQOversqrtPPlus2M)
1545 subtract()
1546 negate()
1547 push_rational(1,2)
1548 power()
1549 simplify()
1550 subtract()
1551 push_integer(2)
1552 divide()
1553
1554 restore()
1555 return
1556
1557
1558
1559
1560
1561
1562 # Q ---------------------------
1563
1564 push(R_determinant)
1565 simplify()
1566 absValFloat()
1567 R_determinant_simplified_toCheckIfZero = pop()
1568
1569 push(R_DELTA0)
1570 simplify()
1571 absValFloat()
1572 R_DELTA0_simplified_toCheckIfZero = pop()
1573
1574 S_CHECKED_AS_NOT_ZERO = false
1575 choiceOfRadicalInQSoSIsNotZero = 0
1576
1577 while !S_CHECKED_AS_NOT_ZERO
1578
1579 Q_CHECKED_AS_NOT_ZERO = false
1580 flipSignOFRadicalSoQIsNotZero = false
1581
1582 # Q will go as denominator to calculate S
1583 # so we have to check
1584 # that is not zero
1585 while !Q_CHECKED_AS_NOT_ZERO
1586
1587 # D1 under the outer radical
1588 push(R_DELTA1)
1589
1590 # D1^2 under the inner radical
1591 push(R_DELTA1)
1592 push_integer(2)
1593 power()
1594
1595 # 4*D0^3 under the inner radical
1596 push_integer(-4)
1597 push(R_DELTA0)
1598 push_integer(3)
1599 power()
1600 multiply()
1601
1602 # addition under the inner radical
1603 add()
1604
1605 # the second radical
1606 push_rational(1,2)
1607 power()
1608
1609 if flipSignOFRadicalSoQIsNotZero
1610 negate()
1611
1612 # the addition under the outer radical
1613 add()
1614
1615 # content of outer radical divided by two
1616 push_integer(2)
1617 divide()
1618
1619 if DEBUG then console.log("content of cubic root: " + stack[tos-1].toString())
1620
1621 # outer radical calculation: cubic root
1622 # now we actually have to find all the roots
1623 # because we have to pick the one that makes S != 0
1624 push_rational(1,3)
1625 power()
1626 simplify()
1627 R_principalCubicRoot = pop()
1628 if DEBUG then console.log("principal cubic root: " + R_principalCubicRoot.toString())
1629
1630 if DEBUG then console.log("tos : " + tos)
1631
1632 if choiceOfRadicalInQSoSIsNotZero == 0
1633 if DEBUG then console.log("chosing principal cubic root")
1634 push(R_principalCubicRoot)
1635 else if choiceOfRadicalInQSoSIsNotZero == 1
1636 if DEBUG then console.log("chosing cubic root beyond principal")
1637 push(R_principalCubicRoot)
1638 push_rational(-1,2)
1639 multiply()
1640
1641 push_integer(3)
1642 push_rational(1,2)
1643 power()
1644 push(imaginaryunit)
1645 multiply()
1646 push_rational(-1,2)
1647 multiply()
1648 push(R_principalCubicRoot)
1649 multiply()
1650 add()
1651 else if choiceOfRadicalInQSoSIsNotZero == 1
1652 if DEBUG then console.log("chosing cubic root beyond beyond principal")
1653 push(R_principalCubicRoot)
1654 push_rational(-1,2)
1655 multiply()
1656
1657 push_integer(3)
1658 push_rational(1,2)
1659 power()
1660 push(imaginaryunit)
1661 multiply()
1662 push_rational(1,2)
1663 multiply()
1664 push(R_principalCubicRoot)
1665 multiply()
1666 add()
1667
1668
1669
1670
1671 simplify()
1672 R_Q = pop()
1673 if DEBUG then console.log "Q " + R_Q.toString()
1674 if DEBUG then console.log("tos: " + tos)
1675
1676
1677 push(R_Q)
1678 simplify()
1679 absValFloat()
1680
1681 R_Q_simplified_toCheckIfZero = pop()
1682 if DEBUG then console.log "Q simplified and abs" + R_Q_simplified_toCheckIfZero.toString()
1683 if iszero(R_Q_simplified_toCheckIfZero) and (!iszero(R_determinant_simplified_toCheckIfZero) and iszero(R_DELTA0_simplified_toCheckIfZero))
1684 if DEBUG then console.log " *********************************** Q IS ZERO and it matters, flipping the sign"
1685 flipSignOFRadicalSoQIsNotZero = true
1686 else
1687 Q_CHECKED_AS_NOT_ZERO = true
1688
1689
1690 if DEBUG then console.log("tos: " + tos)
1691
1692 # S
1693 push_rational(-2,3)
1694 push(R_p)
1695 multiply()
1696
1697 push(R_Q)
1698
1699 push(R_DELTA0)
1700 push(R_Q)
1701 divide()
1702
1703 add()
1704 #rationalize()
1705 #console.log("rationalised: " + stack[tos-1].toString())
1706 #simplify()
1707
1708 push(R_3_a)
1709 divide()
1710
1711 add()
1712
1713 push_rational(1,2)
1714 power()
1715
1716 push_integer(2)
1717 divide()
1718
1719 show_power_debug = true
1720 simplify()
1721 R_S = pop()
1722 if DEBUG then console.log "S " + R_S.toString()
1723
1724 # now check if S is zero
1725 push(R_S)
1726 simplify()
1727 absValFloat()
1728
1729 R_S_simplified_toCheckIfZero = pop()
1730 if DEBUG then console.log "S " + R_S_simplified_toCheckIfZero.toString()
1731 if iszero(R_S_simplified_toCheckIfZero)
1732 if DEBUG then console.log " *********************************** S IS ZERO chosing another cubic root"
1733 choiceOfRadicalInQSoSIsNotZero++
1734 else
1735 S_CHECKED_AS_NOT_ZERO = true
1736
1737 if DEBUG then console.log("tos: " + tos)
1738
1739
1740 # ----------------------------
1741
1742 if DEBUG then console.log("tos: " + tos)
1743
1744 push(p4)
1745 negate()
1746 push(p3)
1747 push_integer(4)
1748 multiply()
1749 divide()
1750 R_minus_b_over_4a = pop()
1751
1752 push_integer(-4)
1753 push(R_S)
1754 push_integer(2)
1755 power()
1756 multiply()
1757 push_integer(2)
1758 push(R_p)
1759 multiply()
1760 subtract()
1761 R_minus_4S2_minus_2p = pop()
1762
1763 push(R_q)
1764 push(R_S)
1765 divide()
1766 R_q_over_S = pop()
1767
1768 if DEBUG then console.log("tos before putting together the 4 solutions: " + tos)
1769
1770 # first solution
1771 push(R_minus_b_over_4a) # first term
1772 push(R_S)
1773 subtract()
1774 push(R_minus_4S2_minus_2p)
1775 push(R_q_over_S)
1776 add()
1777 push_rational(1,2)
1778 power()
1779 push_integer(2)
1780 divide()
1781 add()
1782 simplify()
1783
1784 # second solution
1785 push(R_minus_b_over_4a) # first term
1786 push(R_S)
1787 subtract()
1788 push(R_minus_4S2_minus_2p)
1789 push(R_q_over_S)
1790 add()
1791 push_rational(1,2)
1792 power()
1793 push_integer(2)
1794 divide()
1795 subtract()
1796 simplify()
1797
1798 # third solution
1799 push(R_minus_b_over_4a) # first term
1800 push(R_S)
1801 add()
1802 push(R_minus_4S2_minus_2p)
1803 push(R_q_over_S)
1804 subtract()
1805 push_rational(1,2)
1806 power()
1807 push_integer(2)
1808 divide()
1809 add()
1810 simplify()
1811
1812 # fourth solution
1813 push(R_minus_b_over_4a) # first term
1814 push(R_S)
1815 add()
1816 push(R_minus_4S2_minus_2p)
1817 push(R_q_over_S)
1818 subtract()
1819 push_rational(1,2)
1820 power()
1821 push_integer(2)
1822 divide()
1823 subtract()
1824 simplify()
1825
1826 restore()
1827 return
1828
1829 moveTos tos - n
1830
1831 restore()
1832