(x^2).subs({x: sqrt(x)}) returns sqrt(x) instead of x

[egourgoulhon]

As discussed in https://groups.google.com/d/msg/sage-devel/kYZ4N73oiCM/-iDyQLWyBAAJ, we have currently

sage: (x^2).subs({x: sqrt(x)})                                                  
sqrt(x)

A consequence of this bug is

sage: f(x) = x^2                                                                
sage: f(sqrt(x))                                                                
sqrt(x)

Note that the following outputs are correct:

sage: (x^2).subs({x: sqrt(pi)})                                                 
pi
sage: (x^2).subs({x: x^(1/3)})                                                  
x^(2/3)
sage: y = var('y')                                                              
sage: (x^2).subs({x: sqrt(y)})                                                  
y

Moreover, distinct powers of x seem OK:

sage: (x^3).subs({x: sqrt(x)})                                                  
x^(3/2)
sage: (x^(-2)).subs({x: sqrt(x)})                                               
1/x

But a match is problematic:

sage: (x^3).subs({x:x^(1/3)})
x^(1/3)

CC: @nbruin @rwst @kcrisman @EmmanuelCharpentier

Component: symbolics

Keywords: sqrt, substitution, pynac

Author: Dave Morris

Branch: 8d4db46

Reviewer: Matthias Koeppe

Issue created by migration from https://trac.sagemath.org/ticket/30378

[nbruin]

Description changed:

--- 
+++ 
@@ -22,7 +22,7 @@
 sage: (x^2).subs({x: sqrt(y)})                                                  
 y
 ```
-Moreover, powers of `x` distinct from `2` seem OK:
+Moreover, distinct powers of `x` distinct seem OK:
 
 ```
 sage: (x^3).subs({x: sqrt(x)})                                                  
@@ -30,3 +30,10 @@
 sage: (x^(-2)).subs({x: sqrt(x)})                                               
 1/x
 ```
+But a match is problematic:
+
+```
+sage: (x^3).subs({x:x^(1/3)})
+x^(1/3)
+```
+}}}

[nbruin]

Description changed:

--- 
+++ 
@@ -22,7 +22,7 @@
 sage: (x^2).subs({x: sqrt(y)})                                                  
 y
 ```
-Moreover, distinct powers of `x` distinct seem OK:
+Moreover, distinct powers of `x` seem OK:
 
 ```
 sage: (x^3).subs({x: sqrt(x)})                                                  
@@ -36,4 +36,4 @@
 sage: (x^3).subs({x:x^(1/3)})
 x^(1/3)
 ```
-}}}
+

[behackl]

comment:6

I did some debugging and can report the following: this bug happens in https://github.com/pynac/pynac/blob/85c43c230e8050a1f71898889a021fec28318502/ginac/power.cpp#L747; in fact the correct power in the case of the substitution x=sqrt(x) in x^2 is constructed in p, but then, for a reason I am not quite sure of, the substitution gets applied to the correct output again (L748) and then compared to the original expression, x^2, (L749) where it is decided to return the result with the substitution applied once again and thus sqrt(x) is obtained.

Changing the code such that p from L747 is returned causes one doctest to fail, namely:

sage -t src/sage/symbolic/expression.pyx
**********************************************************************
File "src/sage/symbolic/expression.pyx", line 5150, in sage.symbolic.expression.Expression.substitute
Failed example:
    t.subs(w0 == w0^2)
Expected:
    a^8 + b^8 + (x^2 + y^2)^6
Got:
    a^4 + b^4 + (x^2 + y^2)^3
**********************************************************************
1 item had failures:
   1 of  68 in sage.symbolic.expression.Expression.substitute
    [2871 tests, 1 failure, 34.59 s]

(Context: t = a^2 + b^2 + (x+y)^3 and w0 is a wildcard.)

I don't really see a good solution right now, maybe someone else has some thoughts?

[kcrisman]

comment:7

I am not sure how to fix this either - presumably it was precisely for some sort of wildcard situation that this was in there - but I do note the following in expression.pyx:

            sage: t = a^2 + b^2 + (x+y)^3

...
            sage: t.subs({w0^2: w0^3})
            a^3 + b^3 + (x + y)^3

        Substitute with one or more relational expressions::

            sage: t.subs(w0^2 == w0^3)
            a^3 + b^3 + (x + y)^3

so far, so good

            sage: t.subs(w0 == w0^2)
            a^8 + b^8 + (x^2 + y^2)^6

Wait, what? Shouldn't these be fourth powers? That x and y get replaced by their squares is fine, as is the sixth power, but I'm not at all sure about the eighth powers. (Also note that the wildcard only found literal squares in the {w0^2: w0^3} case, where the cube was not replaced by 1.5 squares - which is presumably good.)

Note that the corresponding Ginac code always returns the subs_one_level business in this situation where either the exponent or the base are 'not trivially equal', whatever that means in this context. Unfortunately I can't find out from the commit history where this changed in Pynac - it seems to go back to Burcin's original rewrite back when we changed to !Pynac/Ginac in the first place for basic symbolics a decade ago.

[behackl]

comment:8

Replying to @kcrisman:

so far, so good

            sage: t.subs(w0 == w0^2)
            a^8 + b^8 + (x^2 + y^2)^6

Wait, what? Shouldn't these be fourth powers? That x and y get replaced by their squares is fine, as is the sixth power, but I'm not at all sure about the eighth powers. (Also note that the wildcard only found literal squares in the {w0^2: w0^3} case, where the cube was not replaced by 1.5 squares - which is presumably good.)

I can only make sense of the eighth powers in a setting where the substitution acts on a power like a^2 in such a way that both the base and the exponent get squared. Whether or not this is the intended behavior for a general wildcard substitution like this seems debateable; one could argue that it is actually inconsistent as the exponent of the parenthesis is not treated in the same way.

Unfortunately, make ptestlong is currently broken for my system, but as far as I can tell from manual testing, there are no other obvious doctest errors other than this wildcard substitution one.

Note that the corresponding Ginac code always returns the subs_one_level business in this situation where either the exponent or the base are 'not trivially equal', whatever that means in this context. Unfortunately I can't find out from the commit history where this changed in Pynac - it seems to go back to Burcin's original rewrite back when we changed to !Pynac/Ginac in the first place for basic symbolics a decade ago.

Interesting! I only saw that something concerning subs was mentioned in #24262, and that the particular segment around L747 has been introduced in the upgrade on that ticket.

[DaveWitteMorris]

comment:9

In the w0 == w0^2 example, I suspect that w0 matches a, so a gets squared, but w0 also matches the entire term a^2, so it gets squared again: ((a<sup>2)</sup>2)^2 = a^8. See, for example:

sage: t.subs(w0 == e^w0)                                                                                                  
e^((e^x + e^y)^3) + e^(e^(2*a)) + e^(e^(2*b))
sage: (cos(x)).subs(w0 == e^w0)                                                                                           
e^(cos(e^x))

[mkoeppe]

comment:10

In the case of substitution where all keys are variables rather than general expressions, shouldn't we be able to write robust code that does not rely on pattern matching?

[kcrisman]

comment:11

Presumably, but just keep in mind the "not reinventing the wheel" guideline. w0 in this example is an intentional wildcard. I'm still not 100% sure what is going on with the simpler case that originated the ticket - again, note the default Ginac behavior.

[mkoeppe]

comment:12

Several things to note here:

• pynac has flags that are supposed to disable pattern matching (subs_options::no_pattern, https://github.com/pynac/pynac/blob/85c43c230e8050a1f71898889a021fec28318502/ginac/flags.h#L61)
• but power::subs does not look at the pattern matching flag
• also Sage's Expression.substitute never passes this or any other flags to Pynac's subs method
• I think any form of wildcarding (and iterated substitution) should be disabled in CallableSymbolicExpressionRing_class._call_element_ (I think it's quite shocking that Sage's callable expressions are built upon general subs)