Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

asked 2021-07-03 22:26:48 +0100

Tony-64 gravatar image

Bug using factor() and complicated constants?

Here is the output of the code below.

Case 1 - original code
10.4887875965689
3.32958130339336
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246651e-11
1.45140834701674e-31

Case 3 - use N(Q), N(R), and do not use factor()
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246651e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the results in case 1 are widely different from those in cases 2-4. Am I doing something wrong or is this indeed a bug in factor()? I am using Sage 9.3.

var('x, y', domain='positive')
var('p, q, r', domain='positive')



P = 1

Q = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R = 9/280*sqrt(448*pi + 2025) + 81/56



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

Q = N(Q)
R = N(R)

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

Q = N(Q)
R = N(R)

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - use N(Q), N(R), and do not use factor()")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))

Bug using factor() and complicated constants?

Here is the output of the code below.

Case 1 - original code
10.4887875965689
3.32958130339336
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246651e-11
1.45140834701674e-31

Case 3 - use N(Q), N(R), and do not use factor()
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246651e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the results in case 1 are widely different from those in cases 2-4. Am I doing something wrong or is this indeed a bug in factor()? I am using Sage 9.3.

var('x, y', domain='positive')
var('p, q, r', domain='positive')



P = 1

Q = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R = 9/280*sqrt(448*pi + 2025) + 81/56



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

Q = N(Q)
R = N(R)

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

Q = N(Q)
R = N(R)

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - use N(Q), N(R), and do not use factor()")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))

Bug using factor() and complicated constants?constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

It looks like factor() doesn't play well with complicated constants like Q in the code below. Here is the output of the code below.output.

Case 1 - original code
10.4887875965689
3.32958130339336

Case 1 - original code; use Q, R
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
10.4887875965689
3.32958130339336
1.39111866114058e-12 1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 4 - use N(Q), R
1.39111866114059e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246651e-11
4.83802782246652e-11
1.45140834701674e-31

Case 3 - use N(Q), N(R), 5 - use Q, N(R)
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
1.45140834701674e-31

Case 6 - use q, r and substitute Q, R after solving
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 7 - use q, r; substitute Q, R after solving, do not use factor()
10.4887875965689
3.32958130339336
1.39111866114058e-12 + 6.95559330500736e-6*I
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246651e-11
4.83802782246652e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the results in output says differently. Going with the simplest case, I would say case 1 are widely different from those in cases 2-4. 2 yields the correct answer. Cases 3, 4 and 7 yield the same answer. Case 4 simplifies Q by taking its a numeric value. Cases 3 and 7 don't use factor(). Case 6 shows that solving the equations is not the issue. Interestingly, case 5 shows an incorrect result for X2 only.

What is going on here? Am I doing something wrong or is this indeed a bug in factor()? I am using Sage 9.3.

var('x, y', domain='positive')
var('p, q, r', domain='positive')



P = 1

Q = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R = 9/280*sqrt(448*pi + 2025) + 81/56


print(N(Q))
print(N(R))



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code")
print(N(Q))
print(N(R))
code; use Q, R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

Q = N(Q)
R = N(R)

eqn_1 = Q N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



 ### Case 3 4 ###

Q = N(Q)
R = N(R)

eqn_1 = Q N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 4 - use N(Q), R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 5 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 5 - use Q, N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 6 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 6 - use q, r and substitute Q, R after solving")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 7 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - use N(Q), N(R), and 7 - use q, r; substitute Q, R after solving, do not use factor()")
print(N(Q))
print(N(R))

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))

Bug using factor() and complicated constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

It looks like factor() doesn't play well with complicated constants like Q in the code below. Here is the output.

10.4887875965689
3.32958130339336

Case 1 - original code; use Q, R
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 4 - use N(Q), R
1.39111866114059e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 5 - use Q, N(R)
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
1.45140834701674e-31

Case 6 - use q, r and substitute Q, R after solving
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 7 - use q, r; substitute Q, R after solving, do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the output says differently. Going with the simplest case, I would say case 2 yields the correct answer. Cases 3, 4 and 7 yield the same answer. Case 4 simplifies Q by taking its a numeric value. Cases 3 and 7 don't use factor(). Case 6 shows that solving the equations is not the issue. Interestingly, case 5 shows an incorrect result for X2 only.

What is going on here? Am I doing something wrong or is this indeed a bug in factor()? I am using Sage 9.3.

var('x, y', domain='positive')
var('p, q, r', domain='positive')



P = 1

Q = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R = 9/280*sqrt(448*pi + 2025) + 81/56


print(N(Q))
print(N(R))



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code; use Q, R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))





### Case 4 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 4 - use N(Q), R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 5 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 5 - use Q, N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 6 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 6 - use q, r and substitute Q, R after solving")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 7 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 7 - use q, r; substitute Q, R after solving, do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))

Bug using factor() and complicated constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

It looks like factor() doesn't play well with complicated constants like Q in the code below. Here is the output.

10.4887875965689
3.32958130339336

Case 1 - original code; use Q, R
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 4 - use N(Q), R
1.39111866114059e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 5 - use Q, N(R)
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
1.45140834701674e-31

Case 6 - use q, r and substitute Q, R after solving
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 7 - use q, r; substitute Q, R after solving, do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the output says differently. Going with the simplest case, I would say case 2 yields the correct answer. Cases 3, 4 and 7 yield the same answer. Case 4 simplifies Q by taking its numeric value. Cases 3 and 7 don't use factor(). Case 6 shows that solving the equations is not the issue. The results of X and Y basically do that as well. Interestingly, case 5 shows an incorrect result for X2 only.

What is going on here? Am I doing something wrong or is this indeed a bug in factor()? I am using Sage 9.3.

var('x, y', domain='positive')
var('p, q, r', domain='positive')



P = 1

Q = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R = 9/280*sqrt(448*pi + 2025) + 81/56


print(N(Q))
print(N(R))



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code; use Q, R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))





### Case 4 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 4 - use N(Q), R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 5 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 5 - use Q, N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 6 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 6 - use q, r and substitute Q, R after solving")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 7 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 7 - use q, r; substitute Q, R after solving, do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))

Bug using factor() and complicated constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

It looks like factor() doesn't play well with complicated constants like Q in the code below. Here is the output.

10.4887875965689
3.32958130339336

Case 1 - original code; use Q, R
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 4 - use N(Q), R
1.39111866114059e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 5 - use Q, N(R)
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
1.45140834701674e-31

Case 6 - use q, r and substitute Q, R after solving
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 7 - use q, r; substitute Q, R after solving, do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the output says differently. Going with the simplest case, I would say case 2 yields the correct answer. Cases 3, 4 and 7 yield the same answer. Case 4 simplifies Q by taking its numeric value. Cases 3 and 7 don't use factor(). Case 6 shows that solving the equations is not the issue. The results of X and Y basically do that so as well. Interestingly, case 5 shows an incorrect result for X2 only.

What is going on here? Am I doing something wrong or is this indeed a bug in factor()? I am using Sage 9.3.

var('x, y', domain='positive')
var('p, q, r', domain='positive')



P = 1

Q = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R = 9/280*sqrt(448*pi + 2025) + 81/56


print(N(Q))
print(N(R))



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code; use Q, R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))





### Case 4 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 4 - use N(Q), R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 5 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 5 - use Q, N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 6 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 6 - use q, r and substitute Q, R after solving")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 7 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 7 - use q, r; substitute Q, R after solving, do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))

Bug using factor() and complicated constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

It looks like factor() doesn't play well with complicated constants like Q in the code below. Here is the output.

10.4887875965689
3.32958130339336

Case 1 - original code; use Q, R
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 4 - use N(Q), R
1.39111866114059e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 5 - use Q, N(R)
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
1.45140834701674e-31

Case 6 - use q, r and substitute Q, R after solving
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 7 - use q, r; substitute Q, R after solving, do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the output says differently. Going with the simplest case, I would say case 2 yields the correct answer. Cases 3, 4 and 7 yield the same answer. Case 4 simplifies Q by taking its numeric value. Cases 3 and 7 don't use factor(). Case 6 shows that solving the equations is not the issue. The results of for X and Y basically do so as well. Interestingly, case 5 shows an incorrect result for X2 only.

What is going on here? Am I doing something wrong or is this indeed a bug in factor()? I am using Sage 9.3.

var('x, y', domain='positive')
var('p, q, r', domain='positive')



P = 1

Q = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R = 9/280*sqrt(448*pi + 2025) + 81/56


print(N(Q))
print(N(R))



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code; use Q, R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))





### Case 4 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 4 - use N(Q), R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 5 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 5 - use Q, N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 6 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 6 - use q, r and substitute Q, R after solving")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 7 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 7 - use q, r; substitute Q, R after solving, do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))

Bug using factor() and complicated constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

Edit 2: simplified to a single example.

It looks like factor() doesn't play well with complicated constants expressions like Q in the code below. Here is the output.

10.4887875965689
3.32958130339336

Case 1 - original code; use Q, R
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 2 - use N(Q), N(R)
1.39111866114058e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 3 - do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 4 - use N(Q), R
1.39111866114059e-12 + 6.95559330500736e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

Case 5 - use Q, N(R)
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
1.45140834701674e-31

Case 6 - use q, r and substitute Q, R after solving
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
0.00139260822082924
4.17782466328442e-24

Case 7 - use q, r; substitute Q, R after solving, do not use factor()
1.39111866114059e-12 + 6.95559330500737e-6*I
-1.23801031100746e-21 - 3.80973535432766e-16*I
4.83802782246652e-11
1.45140834701674e-31

I expect the results for X2 and Y2 (last two lines of each case) to be (nearly) the same in all cases. However, the output says differently. Going with the simplest case, I would say case 2 yields the correct answer. Cases 3, 4 and 7 yield the same answer. Case 4 simplifies Q by taking its numeric value. Cases 3 and 7 don't use factor(). Case 6 shows that solving the equations is not the issue. The results for X and Y basically do so as well. Interestingly, case 5 shows an incorrect result for X2 only.

What is going on here? Am I doing something wrong or is below.

Is this indeed a bug in factor()? known bug? I am using Sage 9.3.

var('x, y', var('x', domain='positive')
var('p, q, r', domain='positive')



P 
f(x) = 1

Q (3/174465461165747500*pi*(-1750000*I*pi*x^3 - 31250000*(224*pi + 45*sqrt(448*pi + 2025) + 2025)*x^2
                                + 17500000000000000*I*pi*x)
        *sqrt(92821652156334811582567480952850314403/10*pi^2/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
              + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 7713517620898636162808584411766250000*pi
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 659225266976959904108326638192187500*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 29665137013963195684874698718648437500
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025))
        /(63*pi^2*x^4 - (504000*I*pi^2 + 67500*I*pi*(sqrt(448*pi + 2025) + 45))*x^3
          - 3000000*(560224*pi^2 + 45*pi*(sqrt(448*pi + 2025) + 45))*x^2 + 8400000000000000000000*pi^2
          - (-5040000000000000*I*pi^2 - 675000000000000*I*pi*(sqrt(448*pi + 2025) + 45))*x))

F = 3/174465461165747500*pi*sqrt(92821652156334811582567480952850314403/10*pi^2
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 7713517620898636162808584411766250000*pi
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 659225266976959904108326638192187500
                                 *sqrt(448*pi + 2025)/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
                                 + 29665137013963195684874698718648437500
                                 /(224*pi + 45*sqrt(448*pi + 2025) + 2025))

R f.real()^2 + f.imag()^2
G = 9/280*sqrt(448*pi + 2025) + 81/56


print(N(Q))
print(N(R))



### Case 1 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 1 - original code; use Q, R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 2 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 2 - use N(Q), N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 3 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 3 - do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))





### Case 4 ###

eqn_1 = N(Q) + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*R - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 4 - use N(Q), R")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 5 ###

eqn_1 = Q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*N(R) - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x]
Y = S[y]

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 5 - use Q, N(R)")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))




### Case 6 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = (X.real()^2 + X.imag()^2).factor()
Y2 = (Y.real()^2 + Y.imag()^2).factor()

print("\nCase 6 - use q, r and substitute Q, R after solving")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))



### Case 7 ###

eqn_1 = q + 1/25000*I*sqrt(3/10)*pi*p*y  == 1/62500*(6000*pi + 3*I*pi*p - 30000000000*I*pi/p)*x
eqn_2 = I*sqrt(3/10)*pi*p*x == (4000*pi + I*pi*p + 25000*r - 10000000000*I*pi/p)*y

S = solve([eqn_1, eqn_2], x, y, solution_dict=True)[0]

X = S[x].subs(q=Q, r=R)
Y = S[y].subs(q=Q, r=R)

X2 = X.real()^2 + X.imag()^2
Y2 = Y.real()^2 + Y.imag()^2

print("\nCase 7 - use q, r; substitute Q, R after solving, do not use factor()")

print(N(X(p=P)))
print(N(Y(p=P)))

print(N(X2(p=P)))
print(N(Y2(p=P)))
(f.real()^2 + f.imag()^2).factor()

print(N(f(1)))
print(N(F(1)))
print(N(G(1)))

Here is the output.

418409.917305475 + 1.28757494213663e11*I
1.65784923163565e22
4.77205703148314e29

Bug using factor() and complicated constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

Edit 2: simplified to a single example.

It looks like factor() doesn't play well with complicated expressions like in the code below.

Is this a known bug? I am using Sage 9.3.

var('x', domain='positive')

f(x) = (3/174465461165747500*pi*(-1750000*I*pi*x^3 - 31250000*(224*pi + 45*sqrt(448*pi + 2025) + 2025)*x^2
                                + 17500000000000000*I*pi*x)
        *sqrt(92821652156334811582567480952850314403/10*pi^2/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
              + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 7713517620898636162808584411766250000*pi
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 659225266976959904108326638192187500*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 29665137013963195684874698718648437500
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025))
        /(63*pi^2*x^4 - (504000*I*pi^2 + 67500*I*pi*(sqrt(448*pi + 2025) + 45))*x^3
          - 3000000*(560224*pi^2 + 45*pi*(sqrt(448*pi + 2025) + 45))*x^2 + 8400000000000000000000*pi^2
          - (-5040000000000000*I*pi^2 - 675000000000000*I*pi*(sqrt(448*pi + 2025) + 45))*x))

F = f.real()^2 + f.imag()^2
G = (f.real()^2 + f.imag()^2).factor()

print(N(f(1)))
print(N(F(1)))
print(N(G(1)))

Here is the output.

418409.917305475 + 1.28757494213663e11*I
1.65784923163565e22
4.77205703148314e29

Bug when using factor() and complicated constants? (UPDATED)

Edit: solved some issues with the code, added a few more cases and improved the description.

Edit 2: simplified to a single example.

It looks like factor() doesn't play well with complicated expressions like in the code below.

Is this a known bug? I am using Sage 9.3.

var('x', domain='positive')

f(x) = (3/174465461165747500*pi*(-1750000*I*pi*x^3 - 31250000*(224*pi + 45*sqrt(448*pi + 2025) + 2025)*x^2
                                + 17500000000000000*I*pi*x)
        *sqrt(92821652156334811582567480952850314403/10*pi^2/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
              + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 7713517620898636162808584411766250000*pi
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 659225266976959904108326638192187500*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 29665137013963195684874698718648437500
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025))
        /(63*pi^2*x^4 - (504000*I*pi^2 + 67500*I*pi*(sqrt(448*pi + 2025) + 45))*x^3
          - 3000000*(560224*pi^2 + 45*pi*(sqrt(448*pi + 2025) + 45))*x^2 + 8400000000000000000000*pi^2
          - (-5040000000000000*I*pi^2 - 675000000000000*I*pi*(sqrt(448*pi + 2025) + 45))*x))

F = f.real()^2 + f.imag()^2
G = (f.real()^2 + f.imag()^2).factor()

print(N(f(1)))
print(N(F(1)))
print(N(G(1)))

Here is the output.

418409.917305475 + 1.28757494213663e11*I
1.65784923163565e22
4.77205703148314e29

Bug when Wrong result after using factor() (UPDATED)factor(): known Sage bug?

Edit: solved some issues with the code, added a few more cases and improved the description.

Edit 2: simplified to a single example.

It looks like factor() doesn't play well with complicated expressions like in the code below.

Is this a known bug? I am using Sage 9.3.

var('x', domain='positive')

f(x) = (3/174465461165747500*pi*(-1750000*I*pi*x^3 - 31250000*(224*pi + 45*sqrt(448*pi + 2025) + 2025)*x^2
                                + 17500000000000000*I*pi*x)
        *sqrt(92821652156334811582567480952850314403/10*pi^2/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
              + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 7713517620898636162808584411766250000*pi
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 659225266976959904108326638192187500*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 29665137013963195684874698718648437500
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025))
        /(63*pi^2*x^4 - (504000*I*pi^2 + 67500*I*pi*(sqrt(448*pi + 2025) + 45))*x^3
          - 3000000*(560224*pi^2 + 45*pi*(sqrt(448*pi + 2025) + 45))*x^2 + 8400000000000000000000*pi^2
          - (-5040000000000000*I*pi^2 - 675000000000000*I*pi*(sqrt(448*pi + 2025) + 45))*x))

F = f.real()^2 + f.imag()^2
G = (f.real()^2 + f.imag()^2).factor()

print(N(f(1)))
print(N(F(1)))
print(N(G(1)))

Here is the output.

418409.917305475 + 1.28757494213663e11*I
1.65784923163565e22
4.77205703148314e29

Wrong result after using factor(): known Sage bug?

It looks like factor() doesn't play well with complicated expressions like in the code below.

Is this a known bug? I am using Sage 9.3.

var('x', domain='positive')

f(x) = (3/174465461165747500*pi*(-1750000*I*pi*x^3 - 31250000*(224*pi + 45*sqrt(448*pi + 2025) + 2025)*x^2
                                + 17500000000000000*I*pi*x)
        *sqrt(92821652156334811582567480952850314403/10*pi^2/(224*pi + 45*sqrt(448*pi + 2025) + 2025)
              + 98489794142024498175862287197250000*pi*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 7713517620898636162808584411766250000*pi
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 659225266976959904108326638192187500*sqrt(448*pi + 2025)
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025) + 29665137013963195684874698718648437500
              /(224*pi + 45*sqrt(448*pi + 2025) + 2025))
        /(63*pi^2*x^4 - (504000*I*pi^2 + 67500*I*pi*(sqrt(448*pi + 2025) + 45))*x^3
          - 3000000*(560224*pi^2 + 45*pi*(sqrt(448*pi + 2025) + 45))*x^2 + 8400000000000000000000*pi^2
          - (-5040000000000000*I*pi^2 - 675000000000000*I*pi*(sqrt(448*pi + 2025) + 45))*x))

F = f.real()^2 + f.imag()^2
G = (f.real()^2 + f.imag()^2).factor()

print(N(f(1)))
print(N(F(1)))
print(N(G(1)))

Here is the output.output. The result for G is clearly incorrect.

418409.917305475 + 1.28757494213663e11*I
1.65784923163565e22
4.77205703148314e29
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
about | faq | help | privacy policy | terms of service
Powered by Askbot version 0.7.59