1831 lines
37 KiB
CoffeeScript
1831 lines
37 KiB
CoffeeScript
|
|
|
|
|
|
#define POLY p1
|
|
#define X p2
|
|
#define A p3
|
|
#define B p4
|
|
#define C p5
|
|
#define Y p6
|
|
|
|
show_power_debug = false
|
|
performing_roots = false
|
|
|
|
Eval_roots = ->
|
|
|
|
# this transforms simple "equation" forms into
|
|
# something that can be processed. E.g., say, y = 3x - 2
|
|
# A == B -> A - B
|
|
# A = B -> A - B
|
|
p2 = cadr(p1)
|
|
if (car(p2) == symbol(SETQ) || car(p2) == symbol(TESTEQ))
|
|
push(cadr(p2))
|
|
Eval()
|
|
push(caddr(p2))
|
|
Eval()
|
|
subtract()
|
|
else
|
|
push(p2)
|
|
Eval()
|
|
p2 = pop()
|
|
if (car(p2) == symbol(SETQ) || car(p2) == symbol(TESTEQ))
|
|
push(cadr(p2))
|
|
Eval()
|
|
push(caddr(p2))
|
|
Eval()
|
|
subtract()
|
|
else
|
|
push(p2)
|
|
|
|
# 2nd arg, x
|
|
|
|
push(caddr(p1))
|
|
Eval()
|
|
p2 = pop()
|
|
if (p2 == symbol(NIL))
|
|
guess()
|
|
else
|
|
push(p2)
|
|
|
|
p2 = pop()
|
|
p1 = pop()
|
|
|
|
if (!ispolyexpandedform(p1, p2))
|
|
stop("roots: 1st argument is not a polynomial in the variable " + p2)
|
|
|
|
push(p1)
|
|
push(p2)
|
|
|
|
roots()
|
|
|
|
|
|
hasImaginaryCoeff = (k) ->
|
|
|
|
#polycoeff = tos
|
|
|
|
|
|
imaginaryCoefficients = false
|
|
h = tos
|
|
for i in [k...0] by -1
|
|
#console.log "hasImaginaryCoeff - coeff.:" + stack[tos-i].toString()
|
|
if iscomplexnumber(stack[tos-i])
|
|
imaginaryCoefficients = true
|
|
break
|
|
return imaginaryCoefficients
|
|
|
|
isSimpleRoot = (k) ->
|
|
|
|
#polycoeff = tos
|
|
|
|
|
|
#tos-n Coefficient of x^0
|
|
#tos-1 Coefficient of x^(n-1)
|
|
|
|
if k > 2
|
|
isSimpleRootPolynomial = true
|
|
h = tos
|
|
|
|
if isZeroAtomOrTensor(stack[tos-k])
|
|
isSimpleRootPolynomial = false
|
|
|
|
for i in [(k-1)...1] by -1
|
|
#console.log "hasImaginaryCoeff - coeff.:" + stack[tos-i].toString()
|
|
if !isZeroAtomOrTensor(stack[tos-i])
|
|
isSimpleRootPolynomial = false
|
|
break
|
|
else
|
|
isSimpleRootPolynomial = false
|
|
|
|
return isSimpleRootPolynomial
|
|
|
|
|
|
normalisedCoeff = (variable, polynomial) ->
|
|
k = coeff(variable, polynomial)
|
|
#console.log("->" + tos)
|
|
divideBy = stack[tos-1]
|
|
miniStack = []
|
|
for i in [1..k]
|
|
miniStack.push(pop())
|
|
#console.log(tos)
|
|
|
|
for i in [(k-1)..0]
|
|
push(miniStack[i])
|
|
push(divideBy)
|
|
divide()
|
|
#console.log(tos)
|
|
return k
|
|
|
|
|
|
|
|
# takes the polynomial and the
|
|
# variable on the stack
|
|
|
|
roots = ->
|
|
h = 0
|
|
i = 0
|
|
n = 0
|
|
|
|
if DEBUG then console.log "roots: " + stack[tos-1].toString() + " " + stack[tos-2].toString()
|
|
|
|
save()
|
|
|
|
# the simplification of nested radicals uses
|
|
# "roots", which in turn uses simplification
|
|
# of nested radicals. Usually there is no problem,
|
|
# one level of recursion does the job. Beyond that,
|
|
# we probably got stuck in a strange case of infinite
|
|
# recursion, so bail out and return NIL.
|
|
if recursionLevelNestedRadicalsRemoval > 1
|
|
pop()
|
|
pop()
|
|
push(symbol(NIL))
|
|
restore()
|
|
return
|
|
|
|
performing_roots = true
|
|
h = tos - 2
|
|
|
|
if DEBUG then console.log("roots checking if " + stack[tos-1].toString() + " is a case of simple roots")
|
|
|
|
p2 = pop()
|
|
p1 = pop()
|
|
|
|
push(p1)
|
|
push(p2)
|
|
|
|
k = normalisedCoeff(p2, p1)
|
|
|
|
if isSimpleRoot(k)
|
|
if DEBUG then console.log("yes, " + stack[tos-1].toString() + " is a case of simple roots")
|
|
lastCoeff = stack[tos-k]
|
|
leadingCoeff = stack[tos-1]
|
|
moveTos tos - k
|
|
pop()
|
|
pop()
|
|
getSimpleRoots(k, leadingCoeff, lastCoeff)
|
|
else
|
|
moveTos tos - k
|
|
roots2()
|
|
|
|
n = tos - h
|
|
if (n == 0)
|
|
stop("roots: the polynomial is not factorable, try nroots")
|
|
if (n == 1)
|
|
performing_roots = false
|
|
restore()
|
|
return
|
|
sort_stack(n)
|
|
p1 = alloc_tensor(n)
|
|
p1.tensor.ndim = 1
|
|
p1.tensor.dim[0] = n
|
|
for i in [0...n]
|
|
p1.tensor.elem[i] = stack[h + i]
|
|
moveTos h
|
|
push(p1)
|
|
restore()
|
|
performing_roots = false
|
|
|
|
# ok to generate these roots take a look at their form
|
|
# in the case of even and odd exponents here:
|
|
# http://www.wolframalpha.com/input/?i=roots+x%5E14+%2B+1
|
|
# http://www.wolframalpha.com/input/?i=roots+ax%5E14+%2B+b
|
|
# http://www.wolframalpha.com/input/?i=roots+x%5E15+%2B+1
|
|
# http://www.wolframalpha.com/input/?i=roots+a*x%5E15+%2B+b
|
|
getSimpleRoots = (n, leadingCoeff, lastCoeff) ->
|
|
if DEBUG then console.log("getSimpleRoots")
|
|
save()
|
|
|
|
#tos-n Coefficient of x^0
|
|
#tos-1 Coefficient of x^(n-1)
|
|
|
|
n = n-1
|
|
|
|
push(lastCoeff)
|
|
push_rational(1,n)
|
|
power()
|
|
|
|
push(leadingCoeff)
|
|
push_rational(1,n)
|
|
power()
|
|
divide()
|
|
|
|
commonPart = pop()
|
|
|
|
if n % 2 == 0
|
|
for rootsOfOne in [1..n] by 2
|
|
push(commonPart)
|
|
push_integer(-1)
|
|
push_rational(rootsOfOne,n)
|
|
power()
|
|
multiply()
|
|
aSol = pop()
|
|
push(aSol)
|
|
push(aSol)
|
|
negate()
|
|
else
|
|
for rootsOfOne in [1..n]
|
|
push(commonPart)
|
|
push_integer(-1)
|
|
push_rational(rootsOfOne,n)
|
|
power()
|
|
multiply()
|
|
if rootsOfOne % 2 == 0
|
|
negate()
|
|
|
|
|
|
restore()
|
|
|
|
roots2 = ->
|
|
save()
|
|
|
|
if DEBUG then console.log "roots2: " + stack[tos-1].toString() + " " + stack[tos-2].toString()
|
|
|
|
p2 = pop() # the polynomial variable
|
|
p1 = pop() # the polynomial
|
|
|
|
push(p1)
|
|
push(p2)
|
|
|
|
k = normalisedCoeff(p2, p1)
|
|
|
|
if !hasImaginaryCoeff(k)
|
|
moveTos tos - k
|
|
factorpoly()
|
|
p1 = pop()
|
|
else
|
|
moveTos tos - k
|
|
pop()
|
|
pop()
|
|
|
|
|
|
if (car(p1) == symbol(MULTIPLY))
|
|
p1 = cdr(p1)
|
|
# scan through all the factors
|
|
# and find the roots of each of them
|
|
while (iscons(p1))
|
|
push(car(p1))
|
|
push(p2)
|
|
roots3()
|
|
p1 = cdr(p1)
|
|
else
|
|
push(p1)
|
|
push(p2)
|
|
roots3()
|
|
|
|
restore()
|
|
|
|
roots3 = ->
|
|
save()
|
|
p2 = pop()
|
|
p1 = pop()
|
|
if (car(p1) == symbol(POWER) && ispolyexpandedform(cadr(p1), p2) && isposint(caddr(p1)))
|
|
n = normalisedCoeff(p2, cadr(p1))
|
|
mini_solve(n)
|
|
else if (ispolyexpandedform(p1, p2))
|
|
n = normalisedCoeff(p2, p1)
|
|
mini_solve(n)
|
|
restore()
|
|
|
|
#-----------------------------------------------------------------------------
|
|
#
|
|
# Input: stack[tos - 2] polynomial
|
|
#
|
|
# stack[tos - 1] dependent symbol
|
|
#
|
|
# Output: stack roots on stack
|
|
#
|
|
# (input args are popped first)
|
|
#
|
|
#-----------------------------------------------------------------------------
|
|
|
|
# note that for many quadratic, cubic and quartic polynomials we don't
|
|
# actually end up using the quadratic/cubic/quartic formulas in here,
|
|
# since there is a chance we factored the polynomial and in so
|
|
# doing we found some solutions and lowered the degree.
|
|
mini_solve = (n) ->
|
|
#console.log "mini_solve >>>>>>>>>>>>>>>>>>>>>>>> tos:" + tos
|
|
|
|
save()
|
|
|
|
# AX + B, X = -B/A
|
|
|
|
if (n == 2)
|
|
#console.log "mini_solve >>>>>>>>> 1st degree"
|
|
p3 = pop()
|
|
p4 = pop()
|
|
push(p4)
|
|
push(p3)
|
|
divide()
|
|
negate()
|
|
restore()
|
|
return
|
|
|
|
# AX^2 + BX + C, X = (-B +/- (B^2 - 4AC)^(1/2)) / (2A)
|
|
|
|
if (n == 3)
|
|
#console.log "mini_solve >>>>>>>>> 2nd degree"
|
|
p3 = pop() # A
|
|
p4 = pop() # B
|
|
p5 = pop() # C
|
|
|
|
# B^2
|
|
push(p4)
|
|
push_integer(2)
|
|
power()
|
|
|
|
# 4AC
|
|
push_integer(4)
|
|
push(p3)
|
|
multiply()
|
|
push(p5)
|
|
multiply()
|
|
|
|
# B^2 - 4AC
|
|
subtract()
|
|
|
|
#(B^2 - 4AC)^(1/2)
|
|
push_rational(1, 2)
|
|
power()
|
|
|
|
#p6 is (B^2 - 4AC)^(1/2)
|
|
p6 = pop()
|
|
push(p6);
|
|
|
|
# B
|
|
push(p4)
|
|
subtract() # -B + (B^2 - 4AC)^(1/2)
|
|
|
|
# 1/2A
|
|
push(p3)
|
|
push_integer(2)
|
|
multiply()
|
|
divide()
|
|
#simplify()
|
|
#rationalize()
|
|
# tos - 1 now is 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
|
|
|
|
push(p6);
|
|
# tos - 1 now is (B^2 - 4AC)^(1/2)
|
|
# tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
|
|
|
|
# add B to tos
|
|
push(p4)
|
|
add()
|
|
# tos - 1 now is B + (B^2 - 4AC)^(1/2)
|
|
# tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
|
|
|
|
negate()
|
|
# tos - 1 now is -B -(B^2 - 4AC)^(1/2)
|
|
# tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
|
|
|
|
# 1/2A again
|
|
push(p3)
|
|
divide()
|
|
push_rational(1, 2)
|
|
multiply()
|
|
#simplify()
|
|
#rationalize()
|
|
# tos - 1: 2nd root: (-B - (B^2 - 4AC)^(1/2)) / (2A)
|
|
# tos - 2: 1st root: (-B + (B^2 - 4AC)^(1/2)) / (2A)
|
|
|
|
restore()
|
|
return
|
|
|
|
#if (n == 4)
|
|
if (n == 4 or n == 5)
|
|
|
|
p3 = pop() # A
|
|
p4 = pop() # B
|
|
p5 = pop() # C
|
|
p6 = pop() # D
|
|
|
|
# C - only related calculations
|
|
push(p5)
|
|
push(p5)
|
|
multiply()
|
|
R_c2 = pop()
|
|
|
|
push(R_c2)
|
|
push(p5)
|
|
multiply()
|
|
R_c3 = pop()
|
|
|
|
# B - only related calculations
|
|
push(p4)
|
|
push(p4)
|
|
multiply()
|
|
R_b2 = pop()
|
|
|
|
push(R_b2)
|
|
push(p4)
|
|
multiply()
|
|
R_b3 = pop()
|
|
|
|
push(R_b3)
|
|
push(p6)
|
|
multiply()
|
|
R_b3_d = pop()
|
|
|
|
push(R_b3_d)
|
|
push_integer(-4)
|
|
multiply()
|
|
R_m4_b3_d = pop()
|
|
|
|
|
|
push(R_b3)
|
|
push_integer(2)
|
|
multiply()
|
|
R_2_b3 = pop()
|
|
|
|
# A - only related calculations
|
|
push(p3)
|
|
push(p3)
|
|
multiply()
|
|
R_a2 = pop()
|
|
|
|
push(R_a2)
|
|
push(p3)
|
|
multiply()
|
|
R_a3 = pop()
|
|
|
|
push_integer(3)
|
|
push(p3)
|
|
multiply()
|
|
R_3_a = pop()
|
|
|
|
push(R_a2)
|
|
push(p6)
|
|
multiply()
|
|
R_a2_d = pop()
|
|
|
|
push(R_a2_d)
|
|
push(p6)
|
|
multiply()
|
|
R_a2_d2 = pop()
|
|
|
|
push(R_a2_d)
|
|
push_integer(27)
|
|
multiply()
|
|
R_27_a2_d = pop()
|
|
|
|
push(R_a2_d2)
|
|
push_integer(-27)
|
|
multiply()
|
|
R_m27_a2_d2 = pop()
|
|
|
|
push(R_3_a)
|
|
push_integer(2)
|
|
multiply()
|
|
R_6_a = pop()
|
|
|
|
# mixed calculations
|
|
push(p3)
|
|
push(p5)
|
|
multiply()
|
|
R_a_c = pop()
|
|
|
|
push(R_a_c)
|
|
push(p4)
|
|
multiply()
|
|
R_a_b_c = pop()
|
|
|
|
push(R_a_b_c)
|
|
push(p6)
|
|
multiply()
|
|
R_a_b_c_d = pop()
|
|
|
|
push(R_a_c)
|
|
push_integer(3)
|
|
multiply()
|
|
R_3_a_c = pop()
|
|
|
|
push_integer(-4)
|
|
push(p3)
|
|
push(R_c3)
|
|
multiply()
|
|
multiply()
|
|
R_m4_a_c3 = pop()
|
|
|
|
push(R_a_b_c)
|
|
push_integer(9)
|
|
multiply()
|
|
negate()
|
|
R_m9_a_b_c = pop()
|
|
|
|
push(R_a_b_c_d)
|
|
push_integer(18)
|
|
multiply()
|
|
R_18_a_b_c_d = pop()
|
|
|
|
push(R_b2)
|
|
push(R_3_a_c)
|
|
subtract()
|
|
R_DELTA0 = pop()
|
|
|
|
push(R_b2)
|
|
push(R_c2)
|
|
multiply()
|
|
R_b2_c2 = pop()
|
|
|
|
push(p4)
|
|
negate()
|
|
push(R_3_a)
|
|
divide()
|
|
R_m_b_over_3a = pop()
|
|
|
|
|
|
if (n == 4)
|
|
if DEBUG then console.log ">>>>>>>>>>>>>>>> actually using cubic formula <<<<<<<<<<<<<<< "
|
|
|
|
#console.log ">>>> A:" + p3.toString()
|
|
#console.log ">>>> B:" + p4.toString()
|
|
#console.log ">>>> C:" + p5.toString()
|
|
#console.log ">>>> D:" + p6.toString()
|
|
|
|
|
|
if DEBUG then console.log "cubic: D0: " + R_DELTA0.toString()
|
|
|
|
push(R_DELTA0)
|
|
push_integer(3)
|
|
power()
|
|
push_integer(4)
|
|
multiply()
|
|
R_4_DELTA03 = pop()
|
|
|
|
push(R_DELTA0)
|
|
simplify()
|
|
absValFloat()
|
|
R_DELTA0_toBeCheckedIfZero = pop()
|
|
if DEBUG then console.log "cubic: D0 as float: " + R_DELTA0_toBeCheckedIfZero.toString()
|
|
#if isZeroAtomOrTensor(R_DELTA0_toBeCheckedIfZero)
|
|
# console.log " *********************************** D0 IS ZERO"
|
|
|
|
|
|
# DETERMINANT
|
|
push(R_18_a_b_c_d)
|
|
push(R_m4_b3_d)
|
|
push(R_b2_c2)
|
|
push(R_m4_a_c3)
|
|
push(R_m27_a2_d2)
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
simplify()
|
|
absValFloat()
|
|
R_determinant = pop()
|
|
if DEBUG then console.log "cubic: DETERMINANT: " + R_determinant.toString()
|
|
|
|
# R_DELTA1
|
|
push(R_2_b3)
|
|
push(R_m9_a_b_c)
|
|
push(R_27_a2_d)
|
|
add()
|
|
add()
|
|
R_DELTA1 = pop()
|
|
if DEBUG then console.log "cubic: D1: " + R_DELTA1.toString()
|
|
|
|
# R_Q
|
|
push(R_DELTA1)
|
|
push_integer(2)
|
|
power()
|
|
push(R_4_DELTA03)
|
|
subtract()
|
|
push_rational(1, 2)
|
|
power()
|
|
simplify()
|
|
R_Q = pop()
|
|
|
|
|
|
if isZeroAtomOrTensor(R_determinant)
|
|
if isZeroAtomOrTensor(R_DELTA0_toBeCheckedIfZero)
|
|
if DEBUG then console.log " cubic: DETERMINANT IS ZERO and delta0 is zero"
|
|
push(R_m_b_over_3a) # just same solution three times
|
|
restore()
|
|
return
|
|
else
|
|
if DEBUG then console.log " cubic: DETERMINANT IS ZERO and delta0 is not zero"
|
|
push(p3)
|
|
push(p6)
|
|
push_integer(9)
|
|
multiply()
|
|
multiply()
|
|
push(p4)
|
|
push(p5)
|
|
multiply()
|
|
subtract()
|
|
push(R_DELTA0)
|
|
push_integer(2)
|
|
multiply()
|
|
divide() # first solution
|
|
root_solution = pop()
|
|
push(root_solution) # pushing two of them on the stack
|
|
push(root_solution)
|
|
|
|
# second solution here
|
|
# 4abc
|
|
push(R_a_b_c)
|
|
push_integer(4)
|
|
multiply()
|
|
|
|
# -9a*a*d
|
|
push(p3)
|
|
push(p3)
|
|
push(p6)
|
|
push_integer(9)
|
|
multiply()
|
|
multiply()
|
|
multiply()
|
|
negate()
|
|
|
|
# -9*b^3
|
|
push(R_b3)
|
|
negate()
|
|
|
|
# sum the three terms
|
|
add()
|
|
add()
|
|
|
|
# denominator is a*delta0
|
|
push(p3)
|
|
push(R_DELTA0)
|
|
multiply()
|
|
|
|
# build the fraction
|
|
divide()
|
|
|
|
restore()
|
|
return
|
|
|
|
|
|
|
|
C_CHECKED_AS_NOT_ZERO = false
|
|
flipSignOFQSoCIsNotZero = false
|
|
|
|
# C will go as denominator, we have to check
|
|
# that is not zero
|
|
while !C_CHECKED_AS_NOT_ZERO
|
|
|
|
# R_C
|
|
push(R_Q)
|
|
if flipSignOFQSoCIsNotZero
|
|
negate()
|
|
push(R_DELTA1)
|
|
add()
|
|
push_rational(1, 2)
|
|
multiply()
|
|
push_rational(1, 3)
|
|
power()
|
|
simplify()
|
|
R_C = pop()
|
|
if DEBUG then console.log "cubic: C: " + R_C.toString()
|
|
|
|
push(R_C)
|
|
simplify()
|
|
absValFloat()
|
|
R_C_simplified_toCheckIfZero = pop()
|
|
if DEBUG then console.log "cubic: C as absval and float: " + R_C_simplified_toCheckIfZero.toString()
|
|
|
|
if isZeroAtomOrTensor(R_C_simplified_toCheckIfZero)
|
|
if DEBUG then console.log " cubic: C IS ZERO flipping the sign"
|
|
flipSignOFQSoCIsNotZero = true
|
|
else
|
|
C_CHECKED_AS_NOT_ZERO = true
|
|
|
|
|
|
push(R_C)
|
|
push(R_3_a)
|
|
multiply()
|
|
R_3_a_C = pop()
|
|
|
|
push(R_3_a_C)
|
|
push_integer(2)
|
|
multiply()
|
|
R_6_a_C = pop()
|
|
|
|
|
|
# imaginary parts calculations
|
|
push(imaginaryunit)
|
|
push_integer(3)
|
|
push_rational(1, 2)
|
|
power()
|
|
multiply()
|
|
i_sqrt3 = pop()
|
|
push_integer(1)
|
|
push(i_sqrt3)
|
|
add()
|
|
one_plus_i_sqrt3 = pop()
|
|
push_integer(1)
|
|
push(i_sqrt3)
|
|
subtract()
|
|
one_minus_i_sqrt3 = pop()
|
|
|
|
|
|
push(R_C)
|
|
push(R_3_a)
|
|
divide()
|
|
R_C_over_3a = pop()
|
|
|
|
# first solution
|
|
push(R_m_b_over_3a) # first term
|
|
push(R_C_over_3a)
|
|
negate() # second term
|
|
push(R_DELTA0)
|
|
push(R_3_a_C)
|
|
divide()
|
|
negate() # third term
|
|
# now add the three terms together
|
|
add()
|
|
add()
|
|
simplify()
|
|
|
|
# second solution
|
|
push(R_m_b_over_3a) # first term
|
|
push(R_C_over_3a)
|
|
push(one_plus_i_sqrt3)
|
|
multiply()
|
|
push_integer(2)
|
|
divide() # second term
|
|
push(one_minus_i_sqrt3)
|
|
push(R_DELTA0)
|
|
multiply()
|
|
push(R_6_a_C)
|
|
divide() # third term
|
|
# now add the three terms together
|
|
add()
|
|
add()
|
|
simplify()
|
|
|
|
# third solution
|
|
push(R_m_b_over_3a) # first term
|
|
push(R_C_over_3a)
|
|
push(one_minus_i_sqrt3)
|
|
multiply()
|
|
push_integer(2)
|
|
divide() # second term
|
|
push(one_plus_i_sqrt3)
|
|
push(R_DELTA0)
|
|
multiply()
|
|
push(R_6_a_C)
|
|
divide() # third term
|
|
# now add the three terms together
|
|
add()
|
|
add()
|
|
simplify()
|
|
|
|
restore()
|
|
return
|
|
|
|
# See http://www.sscc.edu/home/jdavidso/Math/Catalog/Polynomials/Fourth/Fourth.html
|
|
# for a description of general shapes and properties of fourth degree polynomials
|
|
if (n == 5)
|
|
if DEBUG then console.log ">>>>>>>>>>>>>>>> actually using quartic formula <<<<<<<<<<<<<<< "
|
|
p7 = pop() # E
|
|
|
|
if isZeroAtomOrTensor(p4) and isZeroAtomOrTensor(p6) and !isZeroAtomOrTensor(p5) and !isZeroAtomOrTensor(p7)
|
|
if DEBUG then console.log("biquadratic case")
|
|
push(p3)
|
|
push(symbol(SECRETX))
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
|
|
push(p5)
|
|
push(symbol(SECRETX))
|
|
multiply()
|
|
|
|
push(p7)
|
|
|
|
add()
|
|
add()
|
|
|
|
push(symbol(SECRETX))
|
|
roots()
|
|
|
|
biquadraticSolutions = pop()
|
|
|
|
for eachSolution in biquadraticSolutions.tensor.elem
|
|
push(eachSolution)
|
|
push_rational(1,2)
|
|
power()
|
|
simplify()
|
|
|
|
push(eachSolution)
|
|
push_rational(1,2)
|
|
power()
|
|
negate()
|
|
simplify()
|
|
|
|
restore()
|
|
return
|
|
|
|
|
|
|
|
|
|
|
|
# D - only related calculations
|
|
push(p6)
|
|
push(p6)
|
|
multiply()
|
|
R_d2 = pop()
|
|
|
|
# E - only related calculations
|
|
push(p7)
|
|
push(p7)
|
|
multiply()
|
|
R_e2 = pop()
|
|
|
|
push(R_e2)
|
|
push(p7)
|
|
multiply()
|
|
R_e3 = pop()
|
|
|
|
# DETERMINANT
|
|
|
|
push_integer(256)
|
|
push(R_a3)
|
|
push(R_e3)
|
|
multiply()
|
|
multiply() # first term 256 a^3 e^3
|
|
|
|
push_integer(-192)
|
|
push(R_a2_d)
|
|
push(R_e2)
|
|
push(p4)
|
|
multiply()
|
|
multiply()
|
|
multiply() # second term -192 a^3 b d e^2
|
|
|
|
push_integer(-128)
|
|
push(R_a2)
|
|
push(R_c2)
|
|
push(R_e2)
|
|
multiply()
|
|
multiply()
|
|
multiply() # third term -128 a^2 c^2 e^2
|
|
|
|
push_integer(144)
|
|
push(R_a2_d2)
|
|
push(p5)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
multiply() # fourth term 144 a^2 c d^2 e
|
|
|
|
push(R_m27_a2_d2)
|
|
push(R_d2)
|
|
multiply() # fifth term -27 a^2 d^4
|
|
|
|
push_integer(144)
|
|
push(R_a_b_c)
|
|
push(p4)
|
|
push(R_e2)
|
|
multiply()
|
|
multiply()
|
|
multiply() # sixth term 144 a b^2 c e^2
|
|
|
|
push_integer(-6)
|
|
push(p3)
|
|
push(R_b2)
|
|
push(R_d2)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
multiply()
|
|
multiply() # seventh term -6 a b^2 d^2 e
|
|
|
|
push_integer(-80)
|
|
push(R_a_b_c_d)
|
|
push(p5)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
multiply() # eigth term -80 a b c^2 d e
|
|
|
|
push_integer(18)
|
|
push(R_a_b_c_d)
|
|
push(R_d2)
|
|
multiply()
|
|
multiply() # ninth term 18 a b c d^3
|
|
|
|
push_integer(16)
|
|
push(R_a_c)
|
|
push(R_c3)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
multiply() # tenth term 16 a c^4 e
|
|
|
|
push_integer(-4)
|
|
push(R_a_c)
|
|
push(R_c2)
|
|
push(R_d2)
|
|
multiply()
|
|
multiply()
|
|
multiply() # eleventh term -4 a c^3 d^2
|
|
|
|
push_integer(-27)
|
|
push(R_b3)
|
|
push(p4)
|
|
push(R_e2)
|
|
multiply()
|
|
multiply()
|
|
multiply() # twelveth term -27 b^4 e^2
|
|
|
|
push_integer(18)
|
|
push(R_b3_d)
|
|
push(p5)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
multiply() # thirteenth term 18 b^3 c d e
|
|
|
|
push(R_m4_b3_d)
|
|
push(R_d2)
|
|
multiply() # fourteenth term -4 b^3 d^3
|
|
|
|
push_integer(-4)
|
|
push(R_b2_c2)
|
|
push(p5)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
multiply() # fifteenth term -4 b^2 c^3 e
|
|
|
|
push(R_b2_c2)
|
|
push(R_d2)
|
|
multiply() # sixteenth term b^2 c^2 d^2
|
|
|
|
# add together the sixteen terms by doing
|
|
# fifteen adds
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
|
|
R_determinant = pop()
|
|
if DEBUG then console.log("R_determinant: " + R_determinant.toString())
|
|
|
|
# DELTA0
|
|
push(R_c2) # term one of DELTA0
|
|
|
|
push_integer(-3)
|
|
push(p4)
|
|
push(p6)
|
|
multiply()
|
|
multiply() # term two of DELTA0
|
|
|
|
push_integer(12)
|
|
push(p3)
|
|
push(p7)
|
|
multiply()
|
|
multiply() # term three of DELTA0
|
|
|
|
# add the three terms together
|
|
add()
|
|
add()
|
|
|
|
R_DELTA0 = pop()
|
|
if DEBUG then console.log("R_DELTA0: " + R_DELTA0.toString())
|
|
|
|
# DELTA1
|
|
push_integer(2)
|
|
push(R_c3)
|
|
multiply()
|
|
|
|
push_integer(-9)
|
|
push(p4)
|
|
push(p5)
|
|
push(p6)
|
|
multiply()
|
|
multiply()
|
|
multiply()
|
|
|
|
push_integer(27)
|
|
push(R_b2)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
|
|
push_integer(27)
|
|
push(p3)
|
|
push(R_d2)
|
|
multiply()
|
|
multiply()
|
|
|
|
push_integer(-72)
|
|
push(R_a_c)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
|
|
# add the five terms together
|
|
add()
|
|
add()
|
|
add()
|
|
add()
|
|
|
|
R_DELTA1 = pop()
|
|
if DEBUG then console.log("R_DELTA1: " + R_DELTA1.toString())
|
|
|
|
# p
|
|
push_integer(8)
|
|
push(R_a_c)
|
|
multiply()
|
|
|
|
push_integer(-3)
|
|
push(R_b2)
|
|
multiply()
|
|
|
|
add()
|
|
|
|
push_integer(8)
|
|
push(R_a2)
|
|
multiply()
|
|
|
|
divide()
|
|
|
|
R_p = pop()
|
|
if DEBUG then console.log("p: " + R_p.toString())
|
|
|
|
# q
|
|
push(R_b3)
|
|
|
|
push_integer(-4)
|
|
push(R_a_b_c)
|
|
multiply()
|
|
|
|
push_integer(8)
|
|
push(R_a2_d)
|
|
multiply()
|
|
|
|
add()
|
|
add()
|
|
|
|
push_integer(8)
|
|
push(R_a3)
|
|
multiply()
|
|
|
|
divide()
|
|
|
|
R_q = pop()
|
|
if DEBUG then console.log("q: " + R_q.toString())
|
|
|
|
if DEBUG then console.log("tos 1 " + tos)
|
|
if !isZeroAtomOrTensor(p4)
|
|
if DEBUG then console.log("tos 2 " + tos)
|
|
|
|
push_integer(8)
|
|
push(p5)
|
|
push(p3)
|
|
multiply()
|
|
multiply()
|
|
|
|
push_integer(-3)
|
|
push(p4)
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
|
|
add()
|
|
|
|
push_integer(8)
|
|
push(p3)
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
|
|
divide()
|
|
|
|
R_p = pop()
|
|
if DEBUG then console.log("p for depressed quartic: " + R_p.toString())
|
|
|
|
push(p4)
|
|
push_integer(3)
|
|
power()
|
|
|
|
push_integer(-4)
|
|
push(p3)
|
|
push(p4)
|
|
push(p5)
|
|
multiply()
|
|
multiply()
|
|
multiply()
|
|
|
|
push_integer(8)
|
|
push(p6)
|
|
push(p3)
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
multiply()
|
|
|
|
add()
|
|
add()
|
|
|
|
|
|
push_integer(8)
|
|
push(p3)
|
|
push_integer(3)
|
|
power()
|
|
multiply()
|
|
|
|
divide()
|
|
|
|
R_q = pop()
|
|
if DEBUG then console.log("q for depressed quartic: " + R_q.toString())
|
|
|
|
|
|
# convert to depressed quartic
|
|
push(p4)
|
|
push_integer(4)
|
|
power()
|
|
push_integer(-3)
|
|
multiply()
|
|
|
|
push_integer(256)
|
|
push(R_a3)
|
|
push(p7)
|
|
multiply()
|
|
multiply()
|
|
|
|
push_integer(-64)
|
|
push(R_a2_d)
|
|
push(p4)
|
|
multiply()
|
|
multiply()
|
|
|
|
push_integer(16)
|
|
push(R_b2)
|
|
push(p3)
|
|
push(p5)
|
|
multiply()
|
|
multiply()
|
|
multiply()
|
|
|
|
add()
|
|
add()
|
|
add()
|
|
|
|
push_integer(256)
|
|
push(p3)
|
|
push_integer(4)
|
|
power()
|
|
multiply()
|
|
divide()
|
|
|
|
R_r = pop()
|
|
if DEBUG then console.log("r for depressed quartic: " + R_r.toString())
|
|
|
|
if DEBUG then console.log("tos 4 " + tos)
|
|
|
|
push(symbol(SECRETX))
|
|
push_integer(4)
|
|
power()
|
|
if DEBUG then console.log("4 * x^4: " + stack[tos-1].toString())
|
|
|
|
push(R_p)
|
|
push(symbol(SECRETX))
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
if DEBUG then console.log("R_p * x^2: " + stack[tos-1].toString())
|
|
|
|
push(R_q)
|
|
push(symbol(SECRETX))
|
|
multiply()
|
|
if DEBUG then console.log("R_q * x: " + stack[tos-1].toString())
|
|
|
|
push(R_r)
|
|
if DEBUG then console.log("R_r: " + stack[tos-1].toString())
|
|
|
|
add()
|
|
add()
|
|
add()
|
|
|
|
simplify()
|
|
if DEBUG then console.log("solving depressed quartic: " + stack[tos-1].toString())
|
|
|
|
push(symbol(SECRETX))
|
|
|
|
roots()
|
|
|
|
depressedSolutions = pop()
|
|
if DEBUG then console.log("depressedSolutions: " + depressedSolutions)
|
|
|
|
for eachSolution in depressedSolutions.tensor.elem
|
|
push(eachSolution)
|
|
push(p4)
|
|
push_integer(4)
|
|
push(p3)
|
|
multiply()
|
|
divide()
|
|
subtract()
|
|
simplify()
|
|
if DEBUG then console.log("solution from depressed: " + stack[tos-1].toString())
|
|
|
|
restore()
|
|
return
|
|
|
|
else
|
|
R_p = p5
|
|
R_q = p6
|
|
R_r = p7
|
|
|
|
###
|
|
# Descartes' solution
|
|
# https://en.wikipedia.org/wiki/Quartic_function#Descartes.27_solution
|
|
# finding the "u" in the depressed equation
|
|
|
|
push_integer(2)
|
|
push(R_p)
|
|
multiply()
|
|
coeff2 = pop()
|
|
|
|
push_integer(-4)
|
|
push(R_p)
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
push(R_r)
|
|
multiply()
|
|
coeff3 = pop()
|
|
|
|
push(R_q)
|
|
push_integer(2)
|
|
power()
|
|
negate()
|
|
coeff4 = pop()
|
|
|
|
# now build the polynomial
|
|
push(symbol(SECRETX))
|
|
push_integer(3)
|
|
power()
|
|
|
|
push(coeff2)
|
|
push(symbol(SECRETX))
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
|
|
push(coeff3)
|
|
push(symbol(SECRETX))
|
|
multiply()
|
|
|
|
push(coeff4)
|
|
|
|
add()
|
|
add()
|
|
add()
|
|
|
|
console.log("Descarte's resolventCubic: " + stack[tos-1].toString())
|
|
push(symbol(SECRETX))
|
|
|
|
roots()
|
|
|
|
resolventCubicSolutions = pop()
|
|
console.log("Descarte's resolventCubic solutions: " + resolventCubicSolutions)
|
|
console.log("tos: " + tos)
|
|
|
|
R_u = null
|
|
#R_u = resolventCubicSolutions.tensor.elem[1]
|
|
for eachSolution in resolventCubicSolutions.tensor.elem
|
|
console.log("examining solution: " + eachSolution)
|
|
push(eachSolution)
|
|
push_integer(2)
|
|
multiply()
|
|
push(R_p)
|
|
add()
|
|
|
|
absValFloat()
|
|
toBeCheckedIFZero = pop()
|
|
console.log("abs value is: " + eachSolution)
|
|
if !isZeroAtomOrTensor(toBeCheckedIFZero)
|
|
R_u = eachSolution
|
|
break
|
|
|
|
console.log("chosen solution: " + R_u)
|
|
|
|
push(R_u)
|
|
negate()
|
|
R_s = pop()
|
|
|
|
push(R_p)
|
|
push(R_u)
|
|
push_integer(2)
|
|
power()
|
|
push(R_q)
|
|
push(R_u)
|
|
divide()
|
|
add()
|
|
add()
|
|
push_integer(2)
|
|
divide()
|
|
R_t = pop()
|
|
|
|
push(R_p)
|
|
push(R_u)
|
|
push_integer(2)
|
|
power()
|
|
push(R_q)
|
|
push(R_u)
|
|
divide()
|
|
subtract()
|
|
add()
|
|
push_integer(2)
|
|
divide()
|
|
R_v = pop()
|
|
|
|
# factoring the quartic into two quadratics:
|
|
|
|
# now build the polynomial
|
|
push(symbol(SECRETX))
|
|
push_integer(2)
|
|
power()
|
|
|
|
push(R_s)
|
|
push(symbol(SECRETX))
|
|
multiply()
|
|
|
|
push(R_t)
|
|
|
|
add()
|
|
add()
|
|
|
|
console.log("factored quartic 1: " + stack[tos-1].toString())
|
|
|
|
push(symbol(SECRETX))
|
|
push_integer(2)
|
|
power()
|
|
|
|
push(R_u)
|
|
push(symbol(SECRETX))
|
|
multiply()
|
|
|
|
push(R_v)
|
|
|
|
add()
|
|
add()
|
|
|
|
console.log("factored quartic 2: " + stack[tos-1].toString())
|
|
pop()
|
|
|
|
restore()
|
|
return
|
|
###
|
|
|
|
# Ferrari's solution
|
|
# https://en.wikipedia.org/wiki/Quartic_function#Ferrari.27s_solution
|
|
# finding the "m" in the depressed equation
|
|
push_rational(5,2)
|
|
push(R_p)
|
|
multiply()
|
|
coeff2 = pop()
|
|
|
|
push_integer(2)
|
|
push(R_p)
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
push(R_r)
|
|
subtract()
|
|
coeff3 = pop()
|
|
|
|
push(R_p)
|
|
push_integer(3)
|
|
power()
|
|
push_integer(2)
|
|
divide()
|
|
|
|
push_rational(-1,2)
|
|
push(R_p)
|
|
push(R_r)
|
|
multiply()
|
|
multiply()
|
|
|
|
push_rational(-1,8)
|
|
push(R_q)
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
|
|
add()
|
|
add()
|
|
|
|
coeff4 = pop()
|
|
|
|
push(symbol(SECRETX))
|
|
push_integer(3)
|
|
power()
|
|
|
|
push(coeff2)
|
|
push(symbol(SECRETX))
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
|
|
push(coeff3)
|
|
push(symbol(SECRETX))
|
|
multiply()
|
|
|
|
push(coeff4)
|
|
|
|
add()
|
|
add()
|
|
add()
|
|
|
|
if DEBUG then console.log("resolventCubic: " + stack[tos-1].toString())
|
|
push(symbol(SECRETX))
|
|
|
|
roots()
|
|
|
|
resolventCubicSolutions = pop()
|
|
if DEBUG then console.log("resolventCubicSolutions: " + resolventCubicSolutions)
|
|
|
|
R_m = null
|
|
#R_m = resolventCubicSolutions.tensor.elem[1]
|
|
for eachSolution in resolventCubicSolutions.tensor.elem
|
|
if DEBUG then console.log("examining solution: " + eachSolution)
|
|
push(eachSolution)
|
|
push_integer(2)
|
|
multiply()
|
|
push(R_p)
|
|
add()
|
|
|
|
absValFloat()
|
|
toBeCheckedIFZero = pop()
|
|
if DEBUG then console.log("abs value is: " + eachSolution)
|
|
if !isZeroAtomOrTensor(toBeCheckedIFZero)
|
|
R_m = eachSolution
|
|
break
|
|
|
|
if DEBUG then console.log("chosen solution: " + R_m)
|
|
push(R_m)
|
|
push_integer(2)
|
|
multiply()
|
|
push(R_p)
|
|
add()
|
|
push_rational(1,2)
|
|
power()
|
|
simplify()
|
|
sqrtPPlus2M = pop()
|
|
|
|
push(R_q)
|
|
push_integer(2)
|
|
multiply()
|
|
push(sqrtPPlus2M)
|
|
divide()
|
|
simplify()
|
|
TwoQOversqrtPPlus2M = pop()
|
|
|
|
push(R_p)
|
|
push_integer(3)
|
|
multiply()
|
|
push(R_m)
|
|
push_integer(2)
|
|
multiply()
|
|
add()
|
|
ThreePPlus2M = pop()
|
|
|
|
# solution1
|
|
push(sqrtPPlus2M)
|
|
push(ThreePPlus2M)
|
|
push(TwoQOversqrtPPlus2M)
|
|
add()
|
|
negate()
|
|
push_rational(1,2)
|
|
power()
|
|
simplify()
|
|
add()
|
|
push_integer(2)
|
|
divide()
|
|
# solution2
|
|
push(sqrtPPlus2M)
|
|
push(ThreePPlus2M)
|
|
push(TwoQOversqrtPPlus2M)
|
|
add()
|
|
negate()
|
|
push_rational(1,2)
|
|
power()
|
|
simplify()
|
|
subtract()
|
|
push_integer(2)
|
|
divide()
|
|
# solution3
|
|
push(sqrtPPlus2M)
|
|
negate()
|
|
push(ThreePPlus2M)
|
|
push(TwoQOversqrtPPlus2M)
|
|
subtract()
|
|
negate()
|
|
push_rational(1,2)
|
|
power()
|
|
simplify()
|
|
add()
|
|
push_integer(2)
|
|
divide()
|
|
# solution4
|
|
push(sqrtPPlus2M)
|
|
negate()
|
|
push(ThreePPlus2M)
|
|
push(TwoQOversqrtPPlus2M)
|
|
subtract()
|
|
negate()
|
|
push_rational(1,2)
|
|
power()
|
|
simplify()
|
|
subtract()
|
|
push_integer(2)
|
|
divide()
|
|
|
|
restore()
|
|
return
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# Q ---------------------------
|
|
|
|
push(R_determinant)
|
|
simplify()
|
|
absValFloat()
|
|
R_determinant_simplified_toCheckIfZero = pop()
|
|
|
|
push(R_DELTA0)
|
|
simplify()
|
|
absValFloat()
|
|
R_DELTA0_simplified_toCheckIfZero = pop()
|
|
|
|
S_CHECKED_AS_NOT_ZERO = false
|
|
choiceOfRadicalInQSoSIsNotZero = 0
|
|
|
|
while !S_CHECKED_AS_NOT_ZERO
|
|
|
|
Q_CHECKED_AS_NOT_ZERO = false
|
|
flipSignOFRadicalSoQIsNotZero = false
|
|
|
|
# Q will go as denominator to calculate S
|
|
# so we have to check
|
|
# that is not zero
|
|
while !Q_CHECKED_AS_NOT_ZERO
|
|
|
|
# D1 under the outer radical
|
|
push(R_DELTA1)
|
|
|
|
# D1^2 under the inner radical
|
|
push(R_DELTA1)
|
|
push_integer(2)
|
|
power()
|
|
|
|
# 4*D0^3 under the inner radical
|
|
push_integer(-4)
|
|
push(R_DELTA0)
|
|
push_integer(3)
|
|
power()
|
|
multiply()
|
|
|
|
# addition under the inner radical
|
|
add()
|
|
|
|
# the second radical
|
|
push_rational(1,2)
|
|
power()
|
|
|
|
if flipSignOFRadicalSoQIsNotZero
|
|
negate()
|
|
|
|
# the addition under the outer radical
|
|
add()
|
|
|
|
# content of outer radical divided by two
|
|
push_integer(2)
|
|
divide()
|
|
|
|
if DEBUG then console.log("content of cubic root: " + stack[tos-1].toString())
|
|
|
|
# outer radical calculation: cubic root
|
|
# now we actually have to find all the roots
|
|
# because we have to pick the one that makes S != 0
|
|
push_rational(1,3)
|
|
power()
|
|
simplify()
|
|
R_principalCubicRoot = pop()
|
|
if DEBUG then console.log("principal cubic root: " + R_principalCubicRoot.toString())
|
|
|
|
if DEBUG then console.log("tos : " + tos)
|
|
|
|
if choiceOfRadicalInQSoSIsNotZero == 0
|
|
if DEBUG then console.log("chosing principal cubic root")
|
|
push(R_principalCubicRoot)
|
|
else if choiceOfRadicalInQSoSIsNotZero == 1
|
|
if DEBUG then console.log("chosing cubic root beyond principal")
|
|
push(R_principalCubicRoot)
|
|
push_rational(-1,2)
|
|
multiply()
|
|
|
|
push_integer(3)
|
|
push_rational(1,2)
|
|
power()
|
|
push(imaginaryunit)
|
|
multiply()
|
|
push_rational(-1,2)
|
|
multiply()
|
|
push(R_principalCubicRoot)
|
|
multiply()
|
|
add()
|
|
else if choiceOfRadicalInQSoSIsNotZero == 1
|
|
if DEBUG then console.log("chosing cubic root beyond beyond principal")
|
|
push(R_principalCubicRoot)
|
|
push_rational(-1,2)
|
|
multiply()
|
|
|
|
push_integer(3)
|
|
push_rational(1,2)
|
|
power()
|
|
push(imaginaryunit)
|
|
multiply()
|
|
push_rational(1,2)
|
|
multiply()
|
|
push(R_principalCubicRoot)
|
|
multiply()
|
|
add()
|
|
|
|
|
|
|
|
|
|
simplify()
|
|
R_Q = pop()
|
|
if DEBUG then console.log "Q " + R_Q.toString()
|
|
if DEBUG then console.log("tos: " + tos)
|
|
|
|
|
|
push(R_Q)
|
|
simplify()
|
|
absValFloat()
|
|
|
|
R_Q_simplified_toCheckIfZero = pop()
|
|
if DEBUG then console.log "Q simplified and abs" + R_Q_simplified_toCheckIfZero.toString()
|
|
if isZeroAtomOrTensor(R_Q_simplified_toCheckIfZero) and (!isZeroAtomOrTensor(R_determinant_simplified_toCheckIfZero) and isZeroAtomOrTensor(R_DELTA0_simplified_toCheckIfZero))
|
|
if DEBUG then console.log " *********************************** Q IS ZERO and it matters, flipping the sign"
|
|
flipSignOFRadicalSoQIsNotZero = true
|
|
else
|
|
Q_CHECKED_AS_NOT_ZERO = true
|
|
|
|
|
|
if DEBUG then console.log("tos: " + tos)
|
|
|
|
# S
|
|
push_rational(-2,3)
|
|
push(R_p)
|
|
multiply()
|
|
|
|
push(R_Q)
|
|
|
|
push(R_DELTA0)
|
|
push(R_Q)
|
|
divide()
|
|
|
|
add()
|
|
#rationalize()
|
|
#console.log("rationalised: " + stack[tos-1].toString())
|
|
#simplify()
|
|
|
|
push(R_3_a)
|
|
divide()
|
|
|
|
add()
|
|
|
|
push_rational(1,2)
|
|
power()
|
|
|
|
push_integer(2)
|
|
divide()
|
|
|
|
show_power_debug = true
|
|
simplify()
|
|
R_S = pop()
|
|
if DEBUG then console.log "S " + R_S.toString()
|
|
|
|
# now check if S is zero
|
|
push(R_S)
|
|
simplify()
|
|
absValFloat()
|
|
|
|
R_S_simplified_toCheckIfZero = pop()
|
|
if DEBUG then console.log "S " + R_S_simplified_toCheckIfZero.toString()
|
|
if isZeroAtomOrTensor(R_S_simplified_toCheckIfZero)
|
|
if DEBUG then console.log " *********************************** S IS ZERO chosing another cubic root"
|
|
choiceOfRadicalInQSoSIsNotZero++
|
|
else
|
|
S_CHECKED_AS_NOT_ZERO = true
|
|
|
|
if DEBUG then console.log("tos: " + tos)
|
|
|
|
|
|
# ----------------------------
|
|
|
|
if DEBUG then console.log("tos: " + tos)
|
|
|
|
push(p4)
|
|
negate()
|
|
push(p3)
|
|
push_integer(4)
|
|
multiply()
|
|
divide()
|
|
R_minus_b_over_4a = pop()
|
|
|
|
push_integer(-4)
|
|
push(R_S)
|
|
push_integer(2)
|
|
power()
|
|
multiply()
|
|
push_integer(2)
|
|
push(R_p)
|
|
multiply()
|
|
subtract()
|
|
R_minus_4S2_minus_2p = pop()
|
|
|
|
push(R_q)
|
|
push(R_S)
|
|
divide()
|
|
R_q_over_S = pop()
|
|
|
|
if DEBUG then console.log("tos before putting together the 4 solutions: " + tos)
|
|
|
|
# first solution
|
|
push(R_minus_b_over_4a) # first term
|
|
push(R_S)
|
|
subtract()
|
|
push(R_minus_4S2_minus_2p)
|
|
push(R_q_over_S)
|
|
add()
|
|
push_rational(1,2)
|
|
power()
|
|
push_integer(2)
|
|
divide()
|
|
add()
|
|
simplify()
|
|
|
|
# second solution
|
|
push(R_minus_b_over_4a) # first term
|
|
push(R_S)
|
|
subtract()
|
|
push(R_minus_4S2_minus_2p)
|
|
push(R_q_over_S)
|
|
add()
|
|
push_rational(1,2)
|
|
power()
|
|
push_integer(2)
|
|
divide()
|
|
subtract()
|
|
simplify()
|
|
|
|
# third solution
|
|
push(R_minus_b_over_4a) # first term
|
|
push(R_S)
|
|
add()
|
|
push(R_minus_4S2_minus_2p)
|
|
push(R_q_over_S)
|
|
subtract()
|
|
push_rational(1,2)
|
|
power()
|
|
push_integer(2)
|
|
divide()
|
|
add()
|
|
simplify()
|
|
|
|
# fourth solution
|
|
push(R_minus_b_over_4a) # first term
|
|
push(R_S)
|
|
add()
|
|
push(R_minus_4S2_minus_2p)
|
|
push(R_q_over_S)
|
|
subtract()
|
|
push_rational(1,2)
|
|
power()
|
|
push_integer(2)
|
|
divide()
|
|
subtract()
|
|
simplify()
|
|
|
|
restore()
|
|
return
|
|
|
|
moveTos tos - n
|
|
|
|
restore()
|
|
|