Yet another (SR) solver(s) issue

EmmanuelCharpentier commented 1 year ago

Is there an existing issue for this?

[X] I have searched the existing issues for a bug report that matches the one I want to file, without success.

Did you read the documentation and troubleshoot guide?

[X] I have read the documentation and troubleshoot guide

Environment

- **OS**: Linux (Debian testing, `Linux p-202-021 6.1.0-7-amd64 #1 SMP PREEMPT_DYNAMIC Debian 6.1.20-2 (2023-04-08) x86_64 GNU/Linux`)
- **Sage Version**: `SageMath version 10.0.beta9, Release Date: 2023-04-13 Using Python 3.11.2` compiled from source

Steps To Reproduce

See this ask.sagemath.org question.

sage: x, m, n = var('x, m, n')
sage: hh = (exp(x)*(1-x))^n/x^m
sage: Ex = hh.diff(x)
sage: S1 = Ex.solve(x)

Expected Behavior

One can expect the solutions set of the equation.

Actual Behavior

sage: S1
[x^2 == (m*x^2 - m)/n]

Implicit solution. Let's try a recursive call :

sage: S2 = flatten([s.solve(x) for s in S1]) ; S2
[x == -sqrt(m/(m - n)), x == sqrt(m/(m - n))]

Explicit solution, indeed... But :

sage: [Ex.subs(s).is_zero() for s in S2]
[False, False]

doesn't check.

However :

sage: S3 = Ex.factor().solve(x) ; S3
[x == 1, x == 1/2*(m - sqrt(m^2 - 4*m*n))/n, x == 1/2*(m + sqrt(m^2 - 4*m*n))/n, x == 0]
sage: [Ex.subs(s).is_zero() for s in S3]
[True, True, True, True]

Much better...

Additional Information

One notes that S1 doesn't check with S3 (checked) solutions :

sage: [[bool(u.subs(s).simplify_full()) for u in S1] for s in S3]
[[False], [False], [False], [False]]

FWIW, "the competition" doesn't make the same error, bit misses two solutions :

sage: S4 = [{v[1].sage():v[2].sage() for v in s} for s in mathematica.Solve(Ex==0, x)] ; S4
[{x: -1/2*(sqrt(m - 4*n)*sqrt(m) - m)/n},
 {x: 1/2*(sqrt(m - 4*n)*sqrt(m) + m)/n}]
sage: [Ex.subs(s).is_zero() for s in S4]
[True, True]

Sympy isn't much better :

sage: S5 = Ex.solve(x, algorithm="sympy") ; S5
[Complement(ConditionSet(x, Eq(x**m, 0), Complexes), ConditionSet(x, Eq(x**(2*m), 0), Complexes)),
 Complement(ConditionSet(x, Eq((-x*exp(x) + exp(x))**n, 0), Complexes), ConditionSet(x, Eq(x**(2*m), 0), Complexes)),
 Complement(Intersection(Complement({m/(2*n) - sqrt(m*(m - 4*n))/(2*n), m/(2*n) + sqrt(m*(m - 4*n))/(2*n)}, ConditionSet(x, Eq(x**(2*m), 0), Complexes)), {m/(2*n) - sqrt(m*(m - 4*n))/(2*n), m/(2*n) + sqrt(m*(m - 4*n))/(2*n)}), {0, 1})]
sage: [type(u) for u in S5]
[<class 'sympy.sets.sets.Complement'>,
 <class 'sympy.sets.sets.Complement'>,
 <class 'sympy.sets.sets.Complement'>]

and factoring doesn't help it :

sage: S6 = Ex.factor().solve(x, algorithm="sympy") ; S6
[Complement(ConditionSet(x, Eq(x**(-m - 1), 0), Complexes), {1}),
 Complement(ConditionSet(x, Eq((-(x - 1)*exp(x))**n, 0), Complexes), {1}),
 Complement({m/(2*n) - sqrt(m*(m - 4*n))/(2*n), m/(2*n) + sqrt(m*(m - 4*n))/(2*n)}, {1})]

Giac gets the same two roots :

sage: S7 = Ex.solve(x, algorithm="giac") ; S7
[1/2*(m - sqrt(m^2 - 4*m*n))/n, 1/2*(m + sqrt(m^2 - 4*m*n))/n]
sage: [Ex.subs(x==s).is_zero() for s in S7]
[True, True]

Trying algorithm="fricas" gives us the same problematic answer as the default algorithm :

sage: S8 = Ex.solve(x, algorithm="fricas") ; S8
[x^2 == (m*x^2 - m)/n]

Interestingly, this problematic answer is also returned by :

sage: Ex.solve(x, to_poly_solve=True)
[x == -sqrt(m/(m - n)), x == sqrt(m/(m - n))]

HTH,

EmmanuelCharpentier commented 1 year ago

The problem seems to reside in Sage :

sage: reset()
sage: x, m, n = var("x, m, n")
sage: hh = (e^x*(1-x))^n/x^m
sage: maxima_calculus("hh: (%e^x*(1-x))^n/x^m;")
((1-x)*%e^x)^n/x^m
sage: Ex = diff(hh, x);
sage: maxima_calculus("Ex: diff(hh, x);")
(n*%e^-x*((1-x)*%e^x)^n*((1-x)*%e^x-%e^x))/((1-x)*x^m)-m*x^((-m)-1)*((1-x)*%e^x)^n
sage: S1 = solve(Ex, x)
sage: maxima_calculus("S1: solve(Ex, x);")
[x = 1,x = -(sqrt(m^2-4*m*n)-m)/(2*n),x = (sqrt(m^2-4*m*n)+m)/(2*n)]
sage: bool(hh==maxima_calculus("hh").sage())
True
sage: bool(Ex==maxima_calculus("Ex").sage())
True

So far, so good. But :

sage: bool(S1==maxima_calculus("S1").sage())
False

Uh, oh... Let's see :

sage: S1
[x^2 == (m*x^2 - m)/n]
sage: maxima_calculus("S1").sage()
[x == 1,
 x == 1/2*(m - sqrt(m^2 - 4*m*n))/n,
 x == 1/2*(m + sqrt(m^2 - 4*m*n))/n]
sage: [Ex.subs(s).is_zero() for s in S1]
[False]
sage: [Ex.subs(s).is_zero() for s in maxima_calculus("S1").sage()]
[True, True, True]

Here, Sage is wrong and Maxima, used directly through the maxima_calculus interface is (almost) right (it misses the x=0 solution found in the initial report).

This is quite serious, since Sage silently reports wrong results. I raise the severity to critical for this reason.

EmmanuelCharpentier commented 1 year ago

The problem appears to be independent of the interface (library or pexpect) and the way the interface is used :

sage: Ex.maxima_methods().solve(x)
[x^2 == (m*x^2 - m)/n]
sage: maxima_calculus.solve(Ex, x)
[_SAGE_VAR_x^2 = (_SAGE_VAR_m*_SAGE_VAR_x^2-_SAGE_VAR_m)/_SAGE_VAR_n]
sage: [u.sage() for u in  maxima_calculus.solve(Ex, x)]
[x^2 == (m*x^2 - m)/n]
sage: maxima.solve(Ex, x).sage()
[x^2 == (m*x^2 - m)/n]

Attaching the interfaces ~~keyword~~ component...

nbruin commented 1 year ago

No mystery here. If maxima is directly fed the form of Ex as computed it produces the same error. Just a normal error in maxima that should probably be reported. I'm not sure if garbage output from solve should be interpreted as an "implicit solution". Shouldn't solve just solve for x? If it doesn't then I think you have your error there already.

Maxima code that shows the problem:

E : (-m*x^((-m)-1) *(-(x-1)*%e^x)^n) -(n*(-(x-1)*%e^x)^(n-1) *((x-1)*%e^x+%e^x)) /x^m;
solve(E,x);

This is just maxima(Ex) with the variables backsubstituted, to confirm it's not some weird conversion bug due to the variable name decorations. Note that hh.diff(x) does NOT call maxima, so the derivative may not be in the identical format that maxima itself would generate.

sagemath / sage