Ask Your Question
1

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

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

Tony-64 gravatar image

updated 2021-07-05 18:32:48 +0100

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 flag offensive close merge delete

Comments

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
Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2021-07-06 21:31:55 +0100 )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.

Tony-64 gravatar imageTony-64 ( 2021-07-07 15:42:18 +0100 )edit

2 Answers

Sort by ยป oldest newest most voted
0

answered 2021-07-07 22:27:16 +0100

Emmanuel Charpentier gravatar image

updated 2021-07-08 13:27:58 +0100

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.

edit flag offensive delete link more

Comments

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.

Tony-64 gravatar imageTony-64 ( 2021-07-08 17:39:18 +0100 )edit

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

Could you repost it UNformatted ?

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2021-07-08 18:51:33 +0100 )edit

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

Tony-64 gravatar imageTony-64 ( 2021-07-09 16:15:54 +0100 )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.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2021-07-10 00:55:54 +0100 )edit
0

answered 2021-07-09 16:37:56 +0100

Tony-64 gravatar image

updated 2021-07-09 16:43:54 +0100

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())
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

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

Seen: 417 times

Last updated: Jul 09 '21