Repository: sympy/sympy

Issue:
problems with numerical evaluation
I working with a simple statistical mechanics example,
using the Stirling approximation of the factorial in combinatorial functions:

```
from __future__ import division
from sympy.abc import *
from sympy import *
from sympy.physics.units import avogadro_number as N


# Stirling's approximation for factorial
def f(n):
    return sqrt(2*pi*n)*(n/E)**n

# number of k-equipartitions of n particles
def m(n, k):
    return  f(n)/ (f(n/k)**k)

# probability of a k-partition being an equipartition
def p(n,k):
    return m(n, k)/(k**n)

q = (p(n, 2)/p(n, 3)).simplify()

z = log(p(n, k)/p(n, k+1)).expand(force=True)
```

Working tests:

```
>>> print q.n(4)
0.9648*n**0.5
>>> print q.n(4,subs={n: N})
7.487e+11
>>> print z   
k*n*log(n)/(k + 1) - k*n*log(k + 1)/(k + 1) - k*n/(k + 1) + k*log(k)/2 - k*log(k + 1)/2 - n*log(n) + n*log(k + 1) + n + n*log(n)/(k + 1) - n*log(k + 1)/(k + 1) - n/(k + 1) + log(n)/2 - log(k + 1)/2 + log(2)/2 + log(pi)/2
>>> print z.simplify()
k*log(k)/2 - k*log(k + 1)/2 + log(n)/2 - log(k + 1)/2 + log(2)/2 + log(pi)/2
>>> print z.simplify().n(4)
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189
>>> print z.n(4).simplify()
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189
>>> print z.subs(n,N).simplify()
k*log(k)/2 - k*log(k + 1)/2 - log(k + 1)/2 + log(2)/2 + log(pi)/2 + log(602214179000000000000000)/2
>>> print z.subs(n,N).simplify().n(4)
0.5*k*log(k) - 0.5*k*log(k + 1) - 0.5*log(k + 1) + 28.3
```

The problem cases:
- Substitution not working `n` is not replaced by `N` - note that `print q.n(4,subs={n: N})` above works

```
>>> print z.n(4, subs={n: N})
k*n*log(n)/(k + 1.0) - k*n*log(k + 1)/(k + 1.0) - k*n/(k + 1.0) + 0.5*k*log(k) - 0.5*k*log(k + 1) - n*log(n) + n*log(k + 1) + n + n*log(n)/(k + 1.0) - n*log(k + 1)/(k + 1.0) - n/(k + 1.0) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189

>>> print z.simplify().n(4, subs={n: N}) 
0.5*k*log(k) - 0.5*k*log(k + 1) + 0.5*log(n) - 0.5*log(k + 1) + 0.9189
```
- This looses the `log(n)` term completely and performs some other invalid transformations

```
>>> print z.subs(n,N).n(4).simplify()
0.5*k*(k*log(k) - k*log(k + 1) + log(k))/(k + 1)
```


Discussion:

Patch:
diff --git a/doc/src/modules/physics/control/lti.rst b/doc/src/modules/physics/control/lti.rst
index fcd94ef218..5c28064dce 100644
--- a/doc/src/modules/physics/control/lti.rst
+++ b/doc/src/modules/physics/control/lti.rst
@@ -32,5 +32,3 @@ lti
    :members:
 
 .. autofunction:: bilinear
-
-.. autofunction:: backward_diff
diff --git a/pyproject.toml b/pyproject.toml
deleted file mode 100644
index 1f4d34828d..0000000000
--- a/pyproject.toml
+++ /dev/null
@@ -1,51 +0,0 @@
-[tool.ruff]
-# Enable Pyflakes `E` and `F` codes by default.
-select = ["C4", "E", "F"]
-
-# Ignore rules that currently fail on the SymPy codebase
-ignore = [
-    "E401",  # Multiple imports on one line
-    "E402",  # Module level import not at top of file
-    "E501",  # Line too long (<LENGTH> > 88 characters)
-    "E701",  # Multiple statements on one line (colon)
-    "E702",  # Multiple statements on one line (semicolon)
-    "E703",  # Statement ends with an unnecessary semicolon
-    "E711",  # Comparison to `None` should be `cond is not None`
-    "E712",  # Comparison to `<BOOL>` should be `cond is <BOOL>`
-    "E713",  # Test for membership should be `not in`
-    "E714",  # Test for object identity should be `is not`
-    "E721",  # Do not compare types, use `isinstance()`
-    "E722",  # Do not use bare `except`
-    "E731",  # Do not assign a `lambda` expression, use a `def`
-    "E741",  # Ambiguous variable name: `<VARIABLE>`
-    "E743",  # Ambiguous function name: `<FUNCTION>`
-    "F401",  # `<TYPE>` imported but unused
-    "F403",  # `from <MODULE> import *` used; unable to detect undefined names
-    "F405",  # `<TYPE>` may be undefined, or defined from star imports: `<MODULE>`
-    "F523",  # `.format` call has unused arguments at position(s): <INDEX>
-    "F601",  # Dictionary key literal `'<KEY>'` repeated
-    "F811",  # Redefinition of unused `<VARIABLE>` from line <LINE>
-    "F821",  # Undefined name `VARIABLE`
-    "F823",  # Local variable `VARIABLE` referenced before assignment
-    "F841",  # Local variable `VARIABLE` is assigned to but never used
-]
-
-# Exclude paths currently excluded in the flake8 configuration
-exclude = [
-    "sympy/assumptions/*generated.py",
-    "sympy/core/*_generated.py",
-    "sympy/parsing/latex/_antlr/*",
-    "sympy/parsing/autolev/_antlr/*",
-    "sympy/parsing/autolev/test-examples/*",
-    "sympy/integrals/rubi/*",
-]
-
-# Black default, although irrelevant with E501 ignored
-line-length = 88
-
-# Allow unused variables when underscore-prefixed.
-dummy-variable-rgx = "^(_+|(_+[a-zA-Z0-9_]*[a-zA-Z0-9]+?))$"
-
-# Per-file ignores currently specified in the flake8 configuration
-[tool.ruff.per-file-ignores]
-"sympy/interactive/session.py" = ["F821"]
diff --git a/sympy/integrals/laplace.py b/sympy/integrals/laplace.py
index 6515db4673..515afe0bbf 100644
--- a/sympy/integrals/laplace.py
+++ b/sympy/integrals/laplace.py
@@ -1327,15 +1327,7 @@ def _frac(f, s):
     _ILT_rules = [
         (a/s, a, S.true, same, 1),
         (b*(s+a)**(-c), t**(c-1)*exp(-a*t)/gamma(c), c>0, same, 1),
-        (1/(s**2+a**2)**2, (sin(a*t) - a*t*cos(a*t))/(2*a**3),
-         S.true, same, 1),
-        # The next two rules must be there in that order. For the second
-        # one, the condition would be a != 0 or, respectively, to take the
-        # limit a -> 0 after the transform if a == 0. It is much simpler if
-        # the case a == 0 has its own rule.
-        (1/(s**b), t**(b - 1)/gamma(b), S.true, same, 1),
-        (1/(s*(s+a)**b), lowergamma(b, a*t)/(a**b*gamma(b)),
-          S.true, same, 1)
+        (1/(s**2+a**2)**2, (sin(a*t) - a*t*cos(a*t))/(2*a**3), S.true, same, 1)
     ]
     return _ILT_rules, s, t
 
@@ -1424,10 +1416,7 @@ def _inverse_laplace_time_diff(F, s, t, plane):
         debugf('      ma:   %s', (ma1,))
         r, c = _inverse_laplace_transform(ma1[g], s, t, plane)
         r = r.replace(Heaviside(t), 1)
-        if r.has(InverseLaplaceTransform):
-            return diff(r, t, ma1[n]), c
-        else:
-            return Heaviside(t)*diff(r, t, ma1[n]), c
+        return diff(r, t, ma1[n]), c
     return None
 
 
@@ -1539,10 +1528,10 @@ def _inverse_laplace_transform(fn, s_, t_, plane, simplify=True, dorational=True
             pass
         elif (r := _inverse_laplace_apply_simple_rules(f, s_, t_)) is not None:
             pass
-        elif (r := _inverse_laplace_expand(f, s_, t_, plane)) is not None:
-            pass
         elif (r := _inverse_laplace_apply_prog_rules(f, s_, t_, plane)) is not None:
             pass
+        elif (r := _inverse_laplace_expand(f, s_, t_, plane)) is not None:
+            pass
         elif any(undef.has(s_) for undef in f.atoms(AppliedUndef)):
             # If there are undefined functions f(t) then integration is
             # unlikely to do anything useful so we skip it and given an
diff --git a/sympy/physics/control/__init__.py b/sympy/physics/control/__init__.py
index 26d8a87de2..9bdddcb2e3 100644
--- a/sympy/physics/control/__init__.py
+++ b/sympy/physics/control/__init__.py
@@ -1,5 +1,5 @@
 from .lti import (TransferFunction, Series, MIMOSeries, Parallel, MIMOParallel,
-    Feedback, MIMOFeedback, TransferFunctionMatrix, bilinear, backward_diff)
+    Feedback, MIMOFeedback, TransferFunctionMatrix, bilinear)
 from .control_plots import (pole_zero_numerical_data, pole_zero_plot, step_response_numerical_data,
     step_response_plot, impulse_response_numerical_data, impulse_response_plot, ramp_response_numerical_data,
     ramp_response_plot, bode_magnitude_numerical_data, bode_phase_numerical_data, bode_magnitude_plot,
@@ -7,7 +7,7 @@
 
 __all__ = ['TransferFunction', 'Series', 'MIMOSeries', 'Parallel',
     'MIMOParallel', 'Feedback', 'MIMOFeedback', 'TransferFunctionMatrix','bilinear',
-    'backward_diff', 'pole_zero_numerical_data',
+    'pole_zero_numerical_data',
     'pole_zero_plot', 'step_response_numerical_data', 'step_response_plot',
     'impulse_response_numerical_data', 'impulse_response_plot',
     'ramp_response_numerical_data', 'ramp_response_plot',
diff --git a/sympy/physics/control/lti.py b/sympy/physics/control/lti.py
index 71b13803e8..3b87d28c7a 100644
--- a/sympy/physics/control/lti.py
+++ b/sympy/physics/control/lti.py
@@ -22,7 +22,7 @@
 from mpmath.libmp.libmpf import prec_to_dps
 
 __all__ = ['TransferFunction', 'Series', 'MIMOSeries', 'Parallel', 'MIMOParallel',
-    'Feedback', 'MIMOFeedback', 'TransferFunctionMatrix', 'bilinear', 'backward_diff']
+    'Feedback', 'MIMOFeedback', 'TransferFunctionMatrix', 'bilinear']
 
 
 def _roots(poly, var):
@@ -35,8 +35,8 @@ def _roots(poly, var):
 
 def bilinear(tf, sample_per):
         """
-        Returns falling coefficients of H(z) from numerator and denominator.
-        Where H(z) is the corresponding discretized transfer function,
+        Returns falling coeffs of H(z) from numerator and denominator.
+        H(z) is the corresponding discretized transfer function,
         discretized with the bilinear transform method.
         H(z) is obtained from the continuous transfer function H(s)
         by substituting s(z) = 2/T * (z-1)/(z+1) into H(s), where T is the
@@ -44,6 +44,7 @@ def bilinear(tf, sample_per):
         Coefficients are falling, i.e. H(z) = (az+b)/(cz+d) is returned
         as [a, b], [c, d].
 
+
         Examples
         ========
 
@@ -57,10 +58,9 @@ def bilinear(tf, sample_per):
         [2*L + R*T, -2*L + R*T]
         """
 
-
+        z = Symbol('z') # discrete variable z
         T = sample_per  # and sample period T
         s = tf.var
-        z =  s # dummy discrete variable z
 
         np = tf.num.as_poly(s).all_coeffs()
         dp = tf.den.as_poly(s).all_coeffs()
@@ -76,49 +76,6 @@ def bilinear(tf, sample_per):
         return num_coefs, den_coefs
 
 
-def backward_diff(tf, sample_per):
-        """
-        Returns falling coefficients of H(z) from numerator and denominator.
-        Where H(z) is the corresponding discretized transfer function,
-        discretized with the backward difference transform method.
-        H(z) is obtained from the continuous transfer function H(s)
-        by substituting s(z) =  (z-1)/(T*z) into H(s), where T is the
-        sample period.
-        Coefficients are falling, i.e. H(z) = (az+b)/(cz+d) is returned
-        as [a, b], [c, d].
-
-        Examples
-        ========
-
-        >>> from sympy.physics.control.lti import TransferFunction, backward_diff
-        >>> from sympy.abc import s, L, R, T
-        >>> tf = TransferFunction(1, s*L + R, s)
-        >>> numZ, denZ = backward_diff(tf, T)
-        >>> numZ
-        [T, 0]
-        >>> denZ
-        [L + R*T, -L]
-        """
-
-        T = sample_per  # and sample period T
-        s = tf.var
-        z =  s         # dummy discrete variable z
-
-        np = tf.num.as_poly(s).all_coeffs()
-        dp = tf.den.as_poly(s).all_coeffs()
-
-        # The next line results from multiplying H(z) with z^N/z^N
-
-        N = max(len(np), len(dp)) - 1
-        num = Add(*[ T**(N-i)*c*(z-1)**i*(z)**(N-i) for c, i in zip(np[::-1], range(len(np))) ])
-        den = Add(*[ T**(N-i)*c*(z-1)**i*(z)**(N-i) for c, i in zip(dp[::-1], range(len(dp))) ])
-
-        num_coefs = num.as_poly(z).all_coeffs()
-        den_coefs = den.as_poly(z).all_coeffs()
-
-        return num_coefs, den_coefs
-
-
 class LinearTimeInvariant(Basic, EvalfMixin):
     """A common class for all the Linear Time-Invariant Dynamical Systems."""
 
diff --git a/sympy/polys/numberfields/basis.py b/sympy/polys/numberfields/basis.py
index 7c9cb41925..a4ff2df0b2 100644
--- a/sympy/polys/numberfields/basis.py
+++ b/sympy/polys/numberfields/basis.py
@@ -4,6 +4,7 @@
 from sympy.polys.domains.algebraicfield import AlgebraicField
 from sympy.polys.domains.integerring import ZZ
 from sympy.polys.domains.rationalfield import QQ
+from sympy.polys.polyerrors import CoercionFailed
 from sympy.utilities.decorator import public
 from .modules import ModuleEndomorphism, ModuleHomomorphism, PowerBasis
 from .utilities import extract_fundamental_discriminant
@@ -99,8 +100,8 @@ def round_two(T, radicals=None):
     Explanation
     ===========
 
-    Carry out Zassenhaus's "Round 2" algorithm on an irreducible polynomial
-    *T* over :ref:`ZZ` or :ref:`QQ`. This computes an integral basis and the
+    Carry out Zassenhaus's "Round 2" algorithm on a monic irreducible
+    polynomial *T* over :ref:`ZZ`. This computes an integral basis and the
     discriminant for the field $K = \mathbb{Q}[x]/(T(x))$.
 
     Alternatively, you may pass an :py:class:`~.AlgebraicField` instance, in
@@ -152,9 +153,9 @@ def round_two(T, radicals=None):
     ==========
 
     T : :py:class:`~.Poly`, :py:class:`~.AlgebraicField`
-        Either (1) the irreducible polynomial over :ref:`ZZ` or :ref:`QQ`
-        defining the number field, or (2) an :py:class:`~.AlgebraicField`
-        representing the number field itself.
+        Either (1) the irreducible monic polynomial over :ref:`ZZ` defining the
+        number field, or (2) an :py:class:`~.AlgebraicField` representing the
+        number field itself.
 
     radicals : dict, optional
         This is a way for any $p$-radicals (if computed) to be returned by
@@ -190,11 +191,16 @@ def round_two(T, radicals=None):
     K = None
     if isinstance(T, AlgebraicField):
         K, T = T, T.ext.minpoly_of_element()
+    if T.domain == QQ:
+        try:
+            T = Poly(T, domain=ZZ)
+        except CoercionFailed:
+            pass  # Let the error be raised by the next clause.
     if (   not T.is_univariate
         or not T.is_irreducible
-        or T.domain not in [ZZ, QQ]):
-        raise ValueError('Round 2 requires an irreducible univariate polynomial over ZZ or QQ.')
-    T, _ = T.make_monic_over_integers_by_scaling_roots()
+        or not T.is_monic
+        or not T.domain == ZZ):
+        raise ValueError('Round 2 requires a monic irreducible univariate polynomial over ZZ.')
     n = T.degree()
     D = T.discriminant()
     D_modulus = ZZ.from_sympy(abs(D))
diff --git a/sympy/polys/numberfields/galoisgroups.py b/sympy/polys/numberfields/galoisgroups.py
index fb2d778b77..7741fcc60f 100644
--- a/sympy/polys/numberfields/galoisgroups.py
+++ b/sympy/polys/numberfields/galoisgroups.py
@@ -563,7 +563,7 @@ def galois_group(f, *gens, by_name=False, max_tries=30, randomize=False, **args)
     ==========
 
     f : Expr
-        Irreducible polynomial over :ref:`ZZ` or :ref:`QQ`, whose Galois group
+        Irreducible, monic polynomial over :ref:`ZZ`, whose Galois group
         is to be determined.
     gens : optional list of symbols
         For converting *f* to Poly, and will be passed on to the
diff --git a/sympy/polys/rings.py b/sympy/polys/rings.py
index e54b6dbe12..93f3617564 100644
--- a/sympy/polys/rings.py
+++ b/sympy/polys/rings.py
@@ -2441,13 +2441,6 @@ def symmetrize(self):
 
         sympy.polys.polyfuncs.symmetrize
 
-        References
-        ==========
-
-        .. [1] Lauer, E. Algorithms for symmetrical polynomials, Proc. 1976
-            ACM Symp. on Symbolic and Algebraic Computing, NY 242-247.
-            https://dl.acm.org/doi/pdf/10.1145/800205.806342
-
         """
         f = self.copy()
         ring = f.ring
diff --git a/sympy/printing/aesaracode.py b/sympy/printing/aesaracode.py
index 87117e06fa..1806a07213 100644
--- a/sympy/printing/aesaracode.py
+++ b/sympy/printing/aesaracode.py
@@ -17,13 +17,6 @@
     from aesara.tensor.elemwise import Elemwise
     from aesara.tensor.elemwise import DimShuffle
 
-    # `true_divide` replaced `true_div` in Aesara 2.8.11 (released 2023) to
-    # match NumPy
-    # XXX: Remove this when not needed to support older versions.
-    true_divide = getattr(aet, 'true_divide', None)
-    if true_divide is None:
-        true_divide = aet.true_div
-
     mapping = {
             sympy.Add: aet.add,
             sympy.Mul: aet.mul,
@@ -245,8 +238,8 @@ def _print_Piecewise(self, expr, **kwargs):
         return aet.switch(p_cond, p_e, p_remaining)
 
     def _print_Rational(self, expr, **kwargs):
-        return true_divide(self._print(expr.p, **kwargs),
-                           self._print(expr.q, **kwargs))
+        return aet.true_div(self._print(expr.p, **kwargs),
+                            self._print(expr.q, **kwargs))
 
     def _print_Integer(self, expr, **kwargs):
         return expr.p


Test patch:
diff --git a/.github/workflows/runtests.yml b/.github/workflows/runtests.yml
index 981bf09f29..8e93a72f4a 100644
--- a/.github/workflows/runtests.yml
+++ b/.github/workflows/runtests.yml
@@ -22,7 +22,7 @@ jobs:
         with:
           python-version: '3.10'
       - run: python -m pip install --upgrade pip
-      - run: pip install mpmath flake8 flake8-comprehensions ruff
+      - run: pip install mpmath flake8 flake8-comprehensions
 
       - name: Basic code quality tests
         run: bin/test quality
@@ -30,9 +30,6 @@ jobs:
       - name: Run flake8 on the sympy package
         run: flake8 sympy
 
-      - name: Run Ruff on the sympy package
-        run: ruff check .
-
       - name: Detect invalid escapes like '\e'
         run: python -We:invalid -We::SyntaxWarning -m compileall -f -q sympy/
 
diff --git a/sympy/core/tests/test_arit.py b/sympy/core/tests/test_arit.py
index cbce59da8b..abf0f39637 100644
--- a/sympy/core/tests/test_arit.py
+++ b/sympy/core/tests/test_arit.py
@@ -1724,11 +1724,6 @@ def test_Pow_as_coeff_mul_doesnt_expand():
     assert exp(x + y).as_coeff_mul() == (1, (exp(x + y),))
     assert exp(x + exp(x + y)) != exp(x + exp(x)*exp(y))
 
-def test_issue_24751():
-    expr = Add(-2, -3, evaluate=False)
-    expr1 = Add(-1, expr, evaluate=False)
-    assert int(expr1) == int((-3 - 2) - 1)
-
 
 def test_issue_3514_18626():
     assert sqrt(S.Half) * sqrt(6) == 2 * sqrt(3)/2
diff --git a/sympy/integrals/tests/test_laplace.py b/sympy/integrals/tests/test_laplace.py
index 9f5a20d603..85749d03f3 100644
--- a/sympy/integrals/tests/test_laplace.py
+++ b/sympy/integrals/tests/test_laplace.py
@@ -402,11 +402,6 @@ def simp_hyp(expr):
     assert ILT(1/s, s, t) == Heaviside(t)
     assert ILT(a/(a + s), s, t) == a*exp(-a*t)*Heaviside(t)
     assert ILT(s/(a + s), s, t) == -a*exp(-a*t)*Heaviside(t) + DiracDelta(t)
-    assert ILT(s/(a + s)**3, s, t) == t*(-a*t + 4)*exp(-a*t)*Heaviside(t)/2
-    assert ILT(1/(s*(a + s)**3), s, t) == (
-        -a**2*t**2 - 4*a*t + 4*exp(a*t) - 4)*exp(-a*t)*Heaviside(t)/(2*a**3)
-    assert ILT(1/(s*(a + s)**n), s, t) == (
-        Heaviside(t)*lowergamma(n, a*t)/(a**n*gamma(n)))
     assert ILT((s-a)**(-b), s, t) == t**(b - 1)*exp(a*t)*Heaviside(t)/gamma(b)
     assert ILT((a + s)**(-2), s, t) == t*exp(-a*t)*Heaviside(t)
     assert ILT((a + s)**(-5), s, t) == t**4*exp(-a*t)*Heaviside(t)/24
@@ -470,9 +465,10 @@ def simp_hyp(expr):
     # Test time_diff rule
     assert ILT(s**42*f(s), s, t) ==\
         Derivative(InverseLaplaceTransform(f(s), s, t, None), (t, 42))
-    assert ILT((b*s**2 + d)/(a**2 + s**2)**2, s, t) == (
-        a**3*b*t*cos(a*t) + 5*a**2*b*sin(a*t) - a*d*t*cos(a*t) +
-        d*sin(a*t))*Heaviside(t)/(2*a**3)
+    assert ILT((b*s**2 + d)/(a**2 + s**2)**2, s, t) ==\
+        (a**3*b*t*cos(a*t) + 4*a**2*b*sin(a*t)*Heaviside(t) +\
+         a**2*b*sin(a*t) - a*d*t*cos(a*t)*Heaviside(t) +\
+             d*sin(a*t)*Heaviside(t))/(2*a**3)
     assert ILT(cos(s), s, t) == InverseLaplaceTransform(cos(s), s, t, None)
     # Rules for testing different DiracDelta cases
     assert ILT(2, s, t) == 2*DiracDelta(t)
@@ -483,12 +479,6 @@ def simp_hyp(expr):
     assert ILT(exp(2*s), s, t) == InverseLaplaceTransform(exp(2*s), s, t, None)
     r = Symbol('r', real=True)
     assert ILT(exp(r*s), s, t) == InverseLaplaceTransform(exp(r*s), s, t, None)
-    # Rules for testing whether Heaviside(t) is treated properly in diff rule
-    assert ILT(s**2/(a**2 + s**2), s, t) == (
-        -a*sin(a*t)*Heaviside(t) + DiracDelta(t))
-    assert ILT(s**2*(f(s) + 1/(a**2 + s**2)), s, t) == (
-        -a*sin(a*t)*Heaviside(t) + DiracDelta(t) +
-        Derivative(InverseLaplaceTransform(f(s), s, t, None), (t, 2)))
     # Rules from the previous test_inverse_laplace_transform_delta_cond():
     assert ILT(exp(r*s), s, t, noconds=False) ==\
         (InverseLaplaceTransform(exp(r*s), s, t, None), True)
diff --git a/sympy/physics/control/tests/test_lti.py b/sympy/physics/control/tests/test_lti.py
index 5d0f4b67e2..eeea0d4638 100644
--- a/sympy/physics/control/tests/test_lti.py
+++ b/sympy/physics/control/tests/test_lti.py
@@ -14,8 +14,7 @@
 from sympy.core.containers import Tuple
 from sympy.matrices import ImmutableMatrix, Matrix
 from sympy.physics.control import (TransferFunction, Series, Parallel,
-    Feedback, TransferFunctionMatrix, MIMOSeries, MIMOParallel, MIMOFeedback,
-    bilinear, backward_diff)
+    Feedback, TransferFunctionMatrix, MIMOSeries, MIMOParallel, MIMOFeedback, bilinear)
 from sympy.testing.pytest import raises
 
 a, x, b, s, g, d, p, k, a0, a1, a2, b0, b1, b2, tau, zeta, wn, T = symbols('a, x, b, s, g, d, p, k,\
@@ -1232,14 +1231,3 @@ def test_TransferFunction_bilinear():
     tf_test_manual = TransferFunction(s*T+T, s*(T*b+2*a)+T*b-2*a, s)
 
     assert S.Zero == (tf_test_bilinear-tf_test_manual).simplify().num
-
-def test_TransferFunction_backward_diff():
-    # simple transfer function, e.g. ohms law
-    tf = TransferFunction(1, a*s+b, s)
-    numZ, denZ = backward_diff(tf, T)
-    # discretized transfer function with coefs from tf.bilinear()
-    tf_test_bilinear = TransferFunction(s*numZ[0]+numZ[1], s*denZ[0]+denZ[1], s)
-    # corresponding tf with manually calculated coefs
-    tf_test_manual = TransferFunction(s*T, s*(T*b+a)-a, s)
-
-    assert S.Zero == (tf_test_bilinear-tf_test_manual).simplify().num
diff --git a/sympy/polys/numberfields/tests/test_basis.py b/sympy/polys/numberfields/tests/test_basis.py
index c0ed017936..1647608d97 100644
--- a/sympy/polys/numberfields/tests/test_basis.py
+++ b/sympy/polys/numberfields/tests/test_basis.py
@@ -10,9 +10,10 @@
 
 
 def test_round_two():
-    # Poly must be irreducible, and over ZZ or QQ:
+    # Poly must be monic, irreducible, and over ZZ:
+    raises(ValueError, lambda: round_two(Poly(3 * x ** 2 + 1)))
     raises(ValueError, lambda: round_two(Poly(x ** 2 - 1)))
-    raises(ValueError, lambda: round_two(Poly(x ** 2 + sqrt(2))))
+    raises(ValueError, lambda: round_two(Poly(x ** 2 + QQ(1, 2))))
 
     # Test on many fields:
     cases = (
@@ -58,10 +59,6 @@ def test_round_two():
         (x**3 + 9 * x**2 + 6 * x - 8, DM([(1, 0, 0), (0, (1, 2), (1, 2)), (0, 0, 1)], QQ).transpose(), 3969),
         # F = 2^2 * 3^2 * 7
         (x**3 + 15 * x**2 - 9 * x + 13, DM([((1, 6), (1, 3), (1, 6)), (0, 1, 0), (0, 0, 1)], QQ).transpose(), -5292),
-        # Polynomial need not be monic
-        (5*x**3 + 5*x**2 - 10 * x + 40, DM([[1, 0, 0], [0, 1, 0], [0, (1, 2), (1, 2)]], QQ).transpose(), -503),
-        # Polynomial can have non-integer rational coeffs
-        (QQ(5, 3)*x**3 + QQ(5, 3)*x**2 - QQ(10, 3)*x + QQ(40, 3), DM([[1, 0, 0], [0, 1, 0], [0, (1, 2), (1, 2)]], QQ).transpose(), -503),
     )
     for f, B_exp, d_exp in cases:
         K = QQ.alg_field_from_poly(f)
diff --git a/sympy/printing/tests/test_aesaracode.py b/sympy/printing/tests/test_aesaracode.py
index 21484626dc..d820b9fc97 100644
--- a/sympy/printing/tests/test_aesaracode.py
+++ b/sympy/printing/tests/test_aesaracode.py
@@ -30,8 +30,6 @@
     from aesara.tensor.elemwise import Elemwise, DimShuffle
     from aesara.tensor.math import Dot
 
-    from sympy.printing.aesaracode import true_divide
-
     xt, yt, zt = [aet.scalar(name, 'floatX') for name in 'xyz']
     Xt, Yt, Zt = [aet.tensor('floatX', (False, False), name=n) for n in 'XYZ']
 else:
@@ -266,8 +264,8 @@ def test_MatAdd():
 
 
 def test_Rationals():
-    assert theq(aesara_code_(sy.Integer(2) / 3), true_divide(2, 3))
-    assert theq(aesara_code_(S.Half), true_divide(1, 2))
+    assert theq(aesara_code_(sy.Integer(2) / 3), aet.true_div(2, 3))
+    assert theq(aesara_code_(S.Half), aet.true_div(1, 2))
 
 def test_Integers():
     assert aesara_code_(sy.Integer(3)) == 3
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index 4ed4fc47f6..a26e8e33a2 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -1037,6 +1037,30 @@ def test_issue_23543():
     assert (x*(y + z/2)).simplify() == x*(2*y + z)/2
 
 
+def test_issue_11004():
+
+    def f(n):
+        return sqrt(2*pi*n) * (n/E)**n
+
+    def m(n, k):
+        return  f(n) / (f(n/k)**k)
+
+    def p(n,k):
+        return m(n, k) / (k**n)
+
+    N, k = symbols('N k')
+    half = Float('0.5', 4)
+    z = log(p(n, k) / p(n, k + 1)).expand(force=True)
+    r = simplify(z.subs(n, N).n(4))
+    assert r == (
+        half*k*log(k)
+        - half*k*log(k + 1)
+        + half*log(N)
+        - half*log(k + 1)
+        + Float(0.9189224, 4)
+    )
+
+
 def test_issue_19161():
     polynomial = Poly('x**2').simplify()
     assert (polynomial-x**2).simplify() == 0