# 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. The result for G is clearly incorrect.

418409.917305475 + 1.28757494213663e11*I
1.65784923163565e22
4.77205703148314e29

edit retag close merge delete

Tour copy seems incorrect

:sage: f(1).n()
1.39111866114058e-12 + 6.95559330500736e-6*I
sage: F(1).n()
4.83802782246652e-11
sage: G(1).n()
0.00139260822082924

( 2021-07-06 21:31:55 +0200 )edit

Not sure what 'tour copy seems incorrect' means. You can get the reported output by pasting the above code in SageMathCell.

How have your results been obtained? Was there a mistake copying the code? Comparing with the originals, the real part of f(1) is about 300,77227e15 times smaller, while the imaginary part of f(1) is about 18.51136e15 times smaller. The results for F(1) and G(1) both are 342,67046e30 = (18,51136e15)^2 times smaller. Note that the real part of f(1) is not significant for the calculation of F(1) and (G1). The difference between F(1) and G(1) remains and should not be there.

( 2021-07-07 15:42:18 +0200 )edit

Sort by ยป oldest newest most voted

Not an answer, but a comment needing more room than allowed in a comment :

sage: %cpaste
Pasting code; enter '--' alone on the line to stop or use Ctrl-D.
:# Raw cut and paste :
:
: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))
:
:# deleting blanks and newlines (adding a space before "+" and "-") :
:
:g(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))
:
:--
sage: f(1).n()
418409.917305474 + 1.28757494213663e11*I
sage: g(1).n()
1.39111866114058e-12 + 6.95559330500736e-6*I


Something is amiss in your code entry... or Sage's interpretation.

And the problem you noted is still present :

sage: F = f.real()^2 + f.imag()^2
sage: Ff = (f.real()^2 + f.imag()^2).factor()
sage: G = g.real()^2 + g.imag()^2
sage: Gf = (g.real()^2 + g.imag()^2).factor()
sage: F(1).n()
1.65784923163565e22
sage: Ff(1).n()
4.77205703148314e29
sage: G(1).n()
4.83802782246652e-11
sage: Gf(1).n()
0.00139260822082924


This smells of a bug, or rather two bugs :d

• Code interpretation on entry

• Factoring

with a possible common origin (pynac comes to mind...)

I'll report that on sage support, for confirmation and advice on how to file useful tickets...

EDIT : Reported in sage support

EDIT 2: According to Nils Bruin, the firs problem is a preparser problem : the Sage preparser does not (yet) mimic Python handling of incomplete expression lines.

I'm not convinced by its second explanation (divergence due to numerical ill-conditioning), and I'll explore it further.

more

Thanks for confirming. I didn't run into the preparer issue. f(x) resulted from solving a set of equations, each of which fits on a single line easily. I just formatted it across multiple lines for presentation in my question.

( 2021-07-08 17:39:18 +0200 )edit

I just formatted it across multiple lines for presentation in my question.

Could you repost it UNformatted ?

( 2021-07-08 18:51:33 +0200 )edit

The unformatted version is identical to your g(x) above.

( 2021-07-09 16:15:54 +0200 )edit

Pleas see sage-support : Sage has trouble converting your semi-numerical expression into a pure symbolic one. This might hint the source of possible bugs in SR...

One also notes that Sympy correctly converts your expression to an (almost) purely symbolic one.

( 2021-07-10 00:55:54 +0200 )edit

Not an answer, but here is some additional code that avoids the preparser problem. The numbers used in this code are not identical to the one in the original question, but the cases included should provide more details of the issue.

Case 1 doesn't use factor() at all.

Case 2 uses factor() in the second stage only.

Case 3 uses factor() in the first stage only.

Case 4 uses factor() in both stages.

Cases 5 and 6 are the same as cases 3 and 4, but squash the intermediate results to RR.

Cases 2 and 4, where the second stage has to deal with the complicated constants calculated in stage 1, clearly deviate from the other results. In case 6, there is not issue.

Here is the output.

Case 1: 3.58751258113361
Case 2: 4.69382197605726e85
Case 3: 3.58751258113361
Case 4: 1.03265208389561e8
Case 5: 3.58751258113361
Case 6: 3.58751258113361


And here is the code.

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

k = sqrt(3/10)
K = 45*sqrt(448*pi + 2025)

eq1 = x*pi*18607995147 == 4*I*k*(93523937500*q/pi - (11107995147*I + 236271000)*pi*y)
eq2 = y*18607995147 == -10*I*k*((3702665049*I + 157514000) + 3937850000*r)*x

Xc = solve([eq1, eq2], x, y)[0][0].rhs()

X = (Xc.real()^2 + Xc.imag()^2) # WITHOUT FACTOR()

A = 4*pi*sqrt(X*r^2*pi^2)/38
B = 24*pi^4*(1 + 25*r)*X/45125

R = solve(A^2/B == 135/224/pi, r)[0].rhs()
S = solve(B(r=R) == 84*pi^3/1805, q)[0].rhs()

eq1 = 38*S/pi + 4*I*k*pi*p*x == (12*pi + 600*I*pi*p - 600*I*pi/p)*y/125
eq2 = p*y == -I*k*(224*pi + 5600*I*pi*p + K - 5600*I*pi/p + 2025)*x/pi/1680

Xc = solve([eq1, eq2], x, y)[0][0].rhs()

X = (Xc.real()^2 + Xc.imag()^2) # WITHOUT FACTOR()

print('Case 1:', X(1).n())

X = (Xc.real()^2 + Xc.imag()^2).factor()

print('Case 2:', X(1).n())

eq1 = x*pi*18607995147 == 4*I*k*(93523937500*q/pi - (11107995147*I + 236271000)*pi*y)
eq2 = y*18607995147 == -10*I*k*((3702665049*I + 157514000) + 3937850000*r)*x

Xc = solve([eq1, eq2], x, y)[0][0].rhs()

X = (Xc.real()^2 + Xc.imag()^2).factor() # WITH FACTOR()

A = 4*pi*sqrt(X*r^2*pi^2)/38
B = 24*pi^4*(1 + 25*r)*X/45125

R = solve(A^2/B == 135/224/pi, r)[0].rhs()
S = solve(B(r=R) == 84*pi^3/1805, q)[0].rhs()

eq1 = 38*S/pi + 4*I*k*pi*p*x == (12*pi + 600*I*pi*p - 600*I*pi/p)*y/125
eq2 = p*y == -I*k*(224*pi + 5600*I*pi*p + K - 5600*I*pi/p + 2025)*x/pi/1680

Xc = solve([eq1, eq2], x, y)[0][0].rhs()

X = (Xc.real()^2 + Xc.imag()^2) # WITHOUT FACTOR()

print('Case 3:', X(1).n())

X = (Xc.real()^2 + Xc.imag()^2).factor()

print('Case 4:', X(1).n())

eq1 = x*pi*18607995147 == 4*I*k*(93523937500*q/pi - (11107995147*I + 236271000)*pi*y)
eq2 = y*18607995147 == -10*I*k*((3702665049*I + 157514000) + 3937850000*r)*x

Xc = solve([eq1, eq2], x, y)[0][0].rhs()

X = (Xc.real()^2 + Xc.imag()^2).factor() # WITH FACTOR()

A = 4*pi*sqrt(X*r^2*pi^2)/38
B = 24*pi^4*(1 + 25*r)*X/45125

R = (solve(A^2/B == 135/224/pi, r)[0].rhs()).n()
S = (solve(B(r=R) == 84*pi^3/1805, q)[0].rhs()).n()

eq1 = 38*S/pi + 4*I*k*pi*p*x == (12*pi + 600*I*pi*p - 600*I*pi/p)*y/125
eq2 = p*y == -I*k*(224*pi + 5600*I*pi*p + K - 5600*I*pi/p + 2025)*x/pi/1680

Xc = solve([eq1, eq2], x, y)[0][0].rhs()

X = (Xc.real()^2 + Xc.imag()^2) # WITHOUT FACTOR()

print('Case 5:', X(1).n())

X = (Xc.real()^2 + Xc.imag()^2).factor()

print('Case 6:', X(1).n())

more