Trac #31137: Fix subs failure involving integer-valued rational exponents

From [https://ask.sagemath.org/question/54986 Ask Sage question 54986]:
`subs` fails with a memory error (or doesn't return) on some
expressions.

The following works fine:
{{{
sage: _ = var('A, L, G, R, f, k, n, q, u, beta, gamma',
domain="positive")
sage: a = I*R^2*f^3*k*q*A*u
sage: b = 2*pi*L*R^2*G*f^4*k^2*q - 2*pi*L*R^2*G*f^4*q -
2*pi*L*R^2*beta^2*G*q
sage: c = (2*I*pi*L*R^2*beta*gamma*q + 2*I*pi*L*R*(beta + q))*G*f^3
sage: d = 2*(pi*(beta^2 + 1)*L*R^2*q + pi*L*R*beta*gamma*q +
pi*L*beta)*G*f^2
sage: e = (-2*I*pi*L*R^2*beta*gamma*q - 2*I*pi*(beta^2*q +
beta)*L*R)*G*f
sage: expr = a / ((b + c + d + e)*n)
sage: R1 = (expr.real()^2 + expr.imag()^2).factor()
sage: R2 = (sqrt(expr.real()^2 + expr.imag()^2)^2).factor()
sage: R3 = ((sqrt(expr.real()^2 + expr.imag()^2).factor())^2).factor()
sage: all((R1 == R2, R1 == R3, R2 == R3))
True
}}}

These expressions have the same string representation:
{{{
sage: [len(str(ex)) for ex in (R1, R2, R3)]
[734, 734, 734]
sage: all((str(R1) == str(R2), str(R1) == str(R3), str(R2) == str(R3)))
}}}

Substituting `f = 2*beta` in `R1` `R2` `R3` works:
{{{
sage: R1s = R1.subs(f = 2*beta)
sage: R2s = R2.subs(f = 2*beta)
sage: R3s = R3.subs(f = 2*beta)
}}}

However, the following **never return**:
{{{
sage: bool(R3s == R1s)  # hangs
}}}
{{{
sage: len(str(R1s))
520
sage: len(str(R2s))
520
sage: len(str(R3s))  # hangs
}}}

Attempting to interrupt the last call can result in a Sage crash;
the original poster reports a "Memory Error" exception.

As a comparison, substituting in `R3` via !SymPy works:
{{{
sage: import sympy
sage: bool(sympy.sympify(R3).subs(f, 2*beta)._sage_() == R1.subs(f =
2*beta))
True
}}}
but seems to fail via Mathematica somehow:
{{{
sage: mma = mathematica
sage: bool(mma.Replace(R3, mma.Rule(f, 2*beta)).sage() == R1.subs(f =
2*beta))
False
}}}

Priority set to critical, because it affects a very basic feature of
Sage.

-----

RESOLUTION: solved by the patch to pynac in #30446. There was a bug in
evaluating an expression (perhaps `sqrt(expr.real())^2`) that has a
rational exponent (such as `2/2`) that simplifies to an integer. This
ticket just adds a doctest to ensure that the problem is not allowed to
return.

URL: https://trac.sagemath.org/31137
Reported by: charpent
Ticket author(s): Dave Morris
Reviewer(s): Dima Pasechnik

Author: Release Manager

src/sage/symbolic/expression.pyx

@@ -1238,6 +1238,13 @@ cdef class Expression(CommutativeRingElement):
 
             sage: latex(1+x^(2/3)+x^(-2/3))
             x^{\frac{2}{3}} + \frac{1}{x^{\frac{2}{3}}} + 1
+
+        Check that pynac understands rational powers (:trac:`30446`)::
+
+            sage: QQ((24*sqrt(3))^(100/50))==1728
+            True
+            sage: float((24*sqrt(3))^(100/51))
+            1493.0092154...
         """
         return self._parent._latex_element_(self)
 
@@ -4162,6 +4169,20 @@ cdef class Expression(CommutativeRingElement):
             4*pi^2
             sage: exp(-3*ln(-9*x)/3)
             -1/9/x
+
+        Check that :trac:`31137` is also fixed::
+
+            sage: _ = var('A, L, G, R, f, k, n, q, u, beta, gamma', domain="positive")
+            sage: a = I*R^2*f^3*k*q*A*u
+            sage: b = 2*pi*L*R^2*G*f^4*k^2*q - 2*pi*L*R^2*G*f^4*q - 2*pi*L*R^2*beta^2*G*q
+            sage: c = (2*I*pi*L*R^2*beta*gamma*q + 2*I*pi*L*R*(beta + q))*G*f^3
+            sage: d = 2*(pi*(beta^2 + 1)*L*R^2*q + pi*L*R*beta*gamma*q + pi*L*beta)*G*f^2
+            sage: e = (-2*I*pi*L*R^2*beta*gamma*q - 2*I*pi*(beta^2*q + beta)*L*R)*G*f
+            sage: expr = a / ((b + c + d + e)*n)
+            sage: R = ((sqrt(expr.real()^2 + expr.imag()^2).factor())^2).factor()
+            sage: Rs = R.subs(f = 2*beta)
+            sage: len(str(Rs))
+            520
         """
         cdef Expression nexp = <Expression>other
         cdef GEx x