Sage seems to incorrectly evaluate fractional powers of complexes

[Emmanuel Charpentier]

Motivation : see [this post], which shows a case where Sage fails to find the roots of a three equations system.

I signalled in this thread that Sympy was able to find these roots. But I stumbled on a difficulty checking these solutions.

Set up the systems :

# Pretext : https://groups.google.com/g/sage-support/c/Nw12vYR0L0U Unks=var('x,y,l') f(x, y) = 10*x^(1/3)*y^(2/3) g(x, y) = 5*x - 6*y h(x, y) = 5*x^2 + 6*y fx = diff(f,x) fy = diff(f, y) gx = diff(g, x) gy = diff(g, y) hx = diff(h, x) hy = diff(h, y) Sys1 = [fx(x, y)==l*gx(x, y),fy(x, y)==l*gy(x, y),g(x, y)==120] Sys2 = [fx(x, y)==l*hx(x, y),fy(x, y)==l*hy(x, y),h(x, y)==120]

Sage’s (default) solution of the first system

DSol1 = solve(Sys1, Unks, solution_dict=True)

As already shown in the original thread, Sage (slowly) fails to solve the second system. Sympy can solve it (anfd the first one too…) :

SSol1 = solve(Sys1, Unks, algorithm="sympy") SSol2 = solve(Sys2, Unks, algorithm="sympy")

For what it’s worth, chech the “competition” :

MSol1 = [{u[1].sage():u[2].sage() for u in s} for s in mathematica(Sys1).Solve(Unks)] MSol2 = [{u[1].sage():u[2].sage() for u in s} for s in mathematica(Sys2).Solve(Unks)]

Both Sympy and Mathematica find one solution to the first system and two for the second. These solutions have different ex^pressins, but that is not problematic. The fly in the ointment is that Sage has trouble checking these solutions. A simple numerical check of these solutions is :

sage: [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys1] for s in DSol1] [[-0.000977, 8.44 - 4.87*I, 0.000]] sage: [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys1] for s in SSol1] [[7.03 - 4.06*I, 0.000488*I, 0.000]] sage: [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys2] for s in SSol2] [[0.000977, 0.000, 0.000], [-0.000977 + 0.00195*I, 0.000122 - 0.000244*I, 0.000]] sage: [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys1] for s in MSol1] [[7.03 - 4.06*I, 0.000, 0.000]] sage: [[(e.lhs()-e.rhs()).subs(s).n(digits=3) for e in Sys2] for s in MSol2] [[-0.00146 - 0.000977*I, 0.000122 + 0.000244*I, 0.000], [0.00195, 0.000244, 0.000]]

None of these solutions seems to check. Uh oh…

However, using sympy to compute the same numerical checks gives different results :

sage: [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e in Sys1] for s in DSol1] [[-7.03 + 4.06*I, 8.43 - 4.87*I, 0]] sage: [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e in Sys1] for s in SSol1] [[0, 0, 0]] sage: [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e in Sys2] for s in SSol2] [[0, 0, 0], [0, 0, 0]] sage: [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e in Sys1] for s in MSol1] [[0, 0, 0]] sage: [[(e.lhs()-e.rhs())._sympy_().subs(sympy.sympify(s)).simplify().n(3) for e in Sys2] for s in MSol2] [[0, 0, 0], [0, 0, 0]]

Sage’s solution to the first system does not check ; Both Sympy’s and Mathematica’s solutions of both systems check.

This tells us that

• Sage’s solution of the first system is wrong, and that

• Sage’s computaion of the numerical checks of knows solutins is also wrong.

This suggests a bug in Sage’s computations ; this problem may be bound to the use of fractional powers.

Any hint on the “right” way to report this issue eficiently welcome.

​

[Emmanuel Charpentier]

The problem seems Sage-specific : the same systems solve correctly (up to numerical noise) in “pure” Maxima :

;;; Loading #P"/usr/lib/x86_64-linux-gnu/ecl-21.2.1/sb-bsd-sockets.fas" ;;; Loading #P"/usr/lib/x86_64-linux-gnu/ecl-21.2.1/sockets.fas" Maxima 5.46.0 https://maxima.sourceforge.io using Lisp ECL 21.2.1 Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) display2d:false; (%o1) false (%i2) Unks:[x,y,l]; (%o2) [x,y,l] (%i3) f:10*x^(1/3)*y^(2/3); (%o3) 10*x^(1/3)*y^(2/3) (%i4) g:5*x - 6*y; (%o4) 5*x-6*y (%i5) h:5*x^2 + 6*y; (%o5) 6*y+5*x^2 (%i6) fx:diff(f,x); (%o6) (10*y^(2/3))/(3*x^(2/3)) (%i7) fy:diff(f, y); (%o7) (20*x^(1/3))/(3*y^(1/3)) (%i8) gx:diff(g, x); (%o8) 5 (%i9) gy:diff(g, y); (%o9) -6 (%i10) hx:diff(h, x); (%o10) 10*x (%i11) hy:diff(h, y); (%o11) 6 (%i12) Sys1:[fx=l*gx,fy=l*gy,g=120]; (%o12) [(10*y^(2/3))/(3*x^(2/3)) = 5*l,(20*x^(1/3))/(3*y^(1/3)) = -6*l, 5*x-6*y = 120] (%i13) Sys2:[fx=l*hx,fy=l*hy,h=120]; (%o13) [(10*y^(2/3))/(3*x^(2/3)) = 10*l*x,(20*x^(1/3))/(3*y^(1/3)) = 6*l, 6*y+5*x^2 = 120] (%i14) Sol1:solve(Sys1, Unks); (%o14) [[x = 8,y = -40/3,l = (2*5^(2/3))/3^(5/3)]] (%i15) Sol2:solve(Sys2, Unks); (%o15) [[x = (2*sqrt(6))/sqrt(5),y = 16,l = 18750^(1/6)/9], [x = -(2*sqrt(6))/sqrt(5),y = 16,l = -18750^(1/6)/9]] (%i16) map(lambda([w], map(lambda([v], subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))), Sys1))), Sol1); (%o16) [[0,0,0]] (%i17) map(lambda([w], map(lambda([v], subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))), Sys2))), Sol2); (%o17) [[(5^(1/3)*(5*2^(11/3)-(2^(8/3)*5^(1/6)*6^(5/6)*18750^(1/6))/3)) /(3*2^(2/3)*6^(1/3)), -(2^(7/3)*18750^(1/6)-2^(7/3)*5^(5/6)*6^(1/6))/(3*2^(4/3)),0], [(5^(1/3)*(5*2^(11/3)-(2^(8/3)*5^(1/6)*6^(5/6)*18750^(1/6))/3)) /(3*2^(2/3)*6^(1/3)), -(2^(7/3)*5^(5/6)*6^(1/6)-2^(7/3)*18750^(1/6))/(3*2^(4/3)),0]] (%i18) %, numer; (%o18) [[1.404069277432862e-15,0.0,0],[1.404069277432862e-15,0.0,0]]

HTH,

​

[Dima Pasechnik]

You can get the same errors from pure Maxima if you set domain to "complex",

no?

[Emmanuel Charpentier]

Indeed :

(%i29) domain:complex; (%o29) complex (%i30) Sol1C:solve(Sys1, Unks); (%o30) [[x = 8,y = -40/3,l = (2*25^(1/3))/(3*9^(1/3))]] (%i31) Sol2C:solve(Sys2, Unks); (%o31) [] (%i32) map(lambda([w], map(lambda([v], subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))), Sys1))), Sol1C); (%o32) [[((8*(-1)^(2/3)*5^(5/3))/3^(2/3)-(40*25^(1/3))/9^(1/3))/12, ((8*(-1)^(1/3)*3^(2/3)*5^(1/3)*25^(1/3))/9^(1/3)+40) /(2*(-1)^(1/3)*3^(2/3)*5^(1/3)),0]] (%i33) %, numer; (%o33) [[0.08333333333333333*(56.22884435344994 *(0.8660254037844387*%i-0.4999999999999998) -56.22884435344994), 0.1405721108836249*(39.99999999999999 *(0.8660254037844386*%i+0.5000000000000001) +40)*(0.5000000000000001-0.8660254037844386*%i), 0]] (%i34) map(lambda([w], map(lambda([v], subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))), Sys2))), Sol2C); (%o34) [] (%i35) %, numer; (%o35) []

But this does not explain whiy Sage is uneble o check (numericalmlmy or otherwise) the solutions given by Sympy or Mathematica, which check in Sympy (I didn’t yet try to check them in Mathematica, the limitations of the current Mathematica interface make this bothersome…).

IMHO, this would deserve a critiocal ticket (possibly a blocker one), but I do not know how to report this efficiently. Suggestions welcome…

HTH,

​

[Oscar Benjamin]

On Thu, 4 Jan 2024 at 12:15, Emmanuel Charpentier
<emanuel.c...@gmail.com> wrote:
>
> But this does not explain whiy Sage is uneble o check (numericalmlmy or otherwise) the solutions given by Sympy or Mathematica, which check in Sympy (I didn’t yet try to check them in Mathematica, the limitations of the current Mathematica interface make this bothersome…).

The problem is with Sage's subs method:

sage: e1 = (-40/3)**(2/3)
sage: e2 = (y**(2/3)).subs({y:-40/3})
sage: e1
4*(-5/3)^(2/3)
sage: e2
4*(5/3)^(2/3)
sage: e1.n()
-2.81144221767250 + 4.86956076355288*I
sage: e2.n()
5.62288443534499

When the number to be substituted has an exact cube factor subs loses
the sign under the radical:

sage: (y^(2/3)).subs({y:-12})
(-12)^(2/3)
sage: (y^(2/3)).subs({y:-24})
4*3^(2/3)

This is most likely due to confusing y^(n/m) as (y^n)^(1/m) when it
should be (y^(1/m))^n. These are not equivalent if y^(1/m) means the
principal mth root.

--
Oscar

[Dima Pasechnik]

On 4 January 2024 12:15:07 WET, Emmanuel Charpentier <emanuel.c...@gmail.com> wrote:
>
>
>Indeed :
>(%i29) domain:complex; (%o29) complex (%i30) Sol1C:solve(Sys1, Unks);
>(%o30) [[x = 8,y = -40/3,l = (2*25^(1/3))/(3*9^(1/3))]] (%i31)
>Sol2C:solve(Sys2, Unks); (%o31) [] (%i32) map(lambda([w], map(lambda([v],
>subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))), Sys1))), Sol1C);
>(%o32) [[((8*(-1)^(2/3)*5^(5/3))/3^(2/3)-(40*25^(1/3))/9^(1/3))/12,
>((8*(-1)^(1/3)*3^(2/3)*5^(1/3)*25^(1/3))/9^(1/3)+40)
>/(2*(-1)^(1/3)*3^(2/3)*5^(1/3)),0]] (%i33) %, numer; (%o33)
>[[0.08333333333333333*(56.22884435344994
>*(0.8660254037844387*%i-0.4999999999999998) -56.22884435344994),
>0.1405721108836249*(39.99999999999999
>*(0.8660254037844386*%i+0.5000000000000001)
>+40)*(0.5000000000000001-0.8660254037844386*%i), 0]] (%i34) map(lambda([w],
>map(lambda([v], subst(w, v)), map(lambda([u], ratsimp(lhs(u)-rhs(u))),
>Sys2))), Sol2C); (%o34) [] (%i35) %, numer; (%o35) []
>
>But this does not explain whiy Sage is uneble o check (numericalmlmy or
>otherwise) the solutions given by Sympy or Mathematica, which check in
>Sympy (I didn’t yet try to check them in Mathematica, the limitations of
>the current Mathematica interface make this bothersome…).
>
>IMHO, this would deserve a critiocal ticket (possibly a blocker one), but I
>do not know how to report this efficiently. Suggestions welcome…

open an issue with this data there

>>>> <https://groups.google.com/g/sage-support/c/Nw12vYR0L0UUnks=var('x,y,l>')

>>>> -

>>>>
>>>> Sage’s solution of the first system is wrong, and that

>>>> -
>>>>
>>>> Sage’s computaion of the numerical checks of knows solutins is *also*
>>>> wrong.
>>>>
>>>> This suggests a bug in Sage’s computations ; this problem *may* be