Repository: sympy/sympy

Issue:
Asymptotic series expansion for Cosine Integral function
I was just testing out the asymptotic series & limits at infinity for some error functions and I found this 
```
>>> Si(x).series(x, oo)
pi/2 - (-6/x**3 + 1/x + O(x**(-7), (x, oo)))*sin(x)/x - (24/x**4 - 2/x**2 + 1 + O(x**(-7), (x, oo)))*cos(x)/x
>>>
>>> Si(x).series(x, -oo)
-24*cos(x)/x**5 + 6*sin(x)/x**4 + 2*cos(x)/x**3 - sin(x)/x**2 - cos(x)/x - pi/2 + O(x**(-6), (x, -oo))
>>>
>>> Ci(x).series(x, oo)
-(-6/x**3 + 1/x + O(x**(-7), (x, oo)))*cos(x)/x + (24/x**4 - 2/x**2 + 1 + O(x**(-7), (x, oo)))*sin(x)/x
>>>
>>> Ci(x).series(x, -oo)
Traceback (most recent call last):
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 2994, in series
    s = self.subs(x, cdir/x).series(x, n=n, dir='+', cdir=1)
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3015, in series
    rv = self.subs(x, xpos).series(xpos, x0, n, dir, logx=logx, cdir=cdir)
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3023, in series
    s1 = self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
  File "C:\Users\91989\sympy\sympy\sympy\functions\special\error_functions.py", line 1800, in _eval_nseries
    return super()._eval_nseries(x, n, logx)
  File "C:\Users\91989\sympy\sympy\sympy\core\function.py", line 690, in _eval_nseries
    return self._eval_aseries(n, args0, x, logx)
  File "C:\Users\91989\sympy\sympy\sympy\functions\special\error_functions.py", line 2082, in _eval_aseries
    return super(Ci, self)._eval_aseries(n, args0, x, logx)
  File "C:\Users\91989\sympy\sympy\sympy\core\function.py", line 653, in _eval_aseries
    raise PoleError(filldedent('''
sympy.core.function.PoleError:
Asymptotic expansion of Ci around [-oo] is not implemented.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 2999, in series
    s = self.subs(x, cdir*x).aseries(x, n=n)
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3198, in aseries
    return self.subs(x, xpos).aseries(xpos, n, bound, hir).subs(xpos, x)
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3216, in aseries
    s = self.subs(x, exp(x)).aseries(x, n, bound, hir).subs(x, log(x))
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3238, in aseries
    s = func.series(k, 0, n)
  File "C:\Users\91989\sympy\sympy\sympy\core\expr.py", line 3023, in series
    s1 = self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
  File "C:\Users\91989\sympy\sympy\sympy\functions\special\error_functions.py", line 1800, in _eval_nseries
    return super()._eval_nseries(x, n, logx)
  File "C:\Users\91989\sympy\sympy\sympy\core\function.py", line 690, in _eval_nseries
    return self._eval_aseries(n, args0, x, logx)
  File "C:\Users\91989\sympy\sympy\sympy\functions\special\error_functions.py", line 2082, in _eval_aseries
    return super(Ci, self)._eval_aseries(n, args0, x, logx)
  File "C:\Users\91989\sympy\sympy\sympy\core\function.py", line 653, in _eval_aseries
    raise PoleError(filldedent('''
sympy.core.function.PoleError:
Asymptotic expansion of Ci around [-oo] is not implemented.
```

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* functions
  * Fixes for order and missing terms in trigonometric integrals in error functions ([#26290](https://github.com/sympy/sympy/pull/26290) by [@arnabnandikgp](https://github.com/arnabnandikgp))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.13.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    Fixes #26207 
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->


    #### Brief description of what is fixed or changed


    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below between the BEGIN and END
    statements. The basic format is a bulleted list with the name of the subpackage
    and the release note for this PR. For example:

    * solvers
      * Added a new solver for logarithmic equations.

    * functions
      * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

    * physics.units
      * Corrected a semantical error in the conversion between volt and statvolt which
        reported the volt as being larger than the statvolt.

    or if no release note(s) should be included use:

    NO ENTRY

    See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
    information on how to write release notes. The bot will check your release
    notes automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->

    * functions
      * Fixes for order and missing terms in trigonometric integrals in error functions
 
    <!-- END RELEASE NOTES -->


</details><p>


Patch:
diff --git a/doc/src/modules/physics/vector/advanced.rst b/doc/src/modules/physics/vector/advanced.rst
index 316487f64a..4f0f922d62 100644
--- a/doc/src/modules/physics/vector/advanced.rst
+++ b/doc/src/modules/physics/vector/advanced.rst
@@ -198,34 +198,3 @@ so dynamic symbols created before or after will print the same way.
 Also note that ``Vector``'s ``.dt`` method uses the ``._t`` attribute of
 ``dynamicsymbols``, along with a number of other important functions and
 methods. Don't mix and match symbols representing time.
-
-Solving Vector Equations
-========================
-
-To solve equations involving vectors, you cannot directly use the solve
-functions on a vector. Instead, you must convert the vector to a set of scalar
-equations.
-
-Suppose that we have two frames ``N`` and ``A``, where ``A`` is rotated 30
-degrees about the z-axis with respect to ``N``. ::
-
-  >>> from sympy import pi, symbols, solve
-  >>> from sympy.physics.vector import ReferenceFrame
-  >>> N = ReferenceFrame("N")
-  >>> A = ReferenceFrame("A")
-  >>> A.orient_axis(N, pi / 6, N.z)
-
-Suppose that we have two vectors ``v1`` and ``v2``, which represent the same
-vector using different symbols. ::
-
-  >>> v1x, v1y, v1z = symbols("v1x v1y v1z")
-  >>> v2x, v2y, v2z = symbols("v2x v2y v2z")
-  >>> v1 = v1x * N.x + v1y * N.y + v1z * N.z
-  >>> v2 = v2x * A.x + v2y * A.y + v2z * A.z
-
-Our goal is to find the relationship between the symbols used in ``v2`` and the
-symbols used in ``v1``. We can achieve this by converting the vector to a matrix
-and then solving the matrix using :meth:`sympy.solvers.solvers.solve`. ::
-
-  >>> solve((v1 - v2).to_matrix(N), [v2x, v2y, v2z])
-  {v2x: sqrt(3)*v1x/2 + v1y/2, v2y: -v1x/2 + sqrt(3)*v1y/2, v2z: v1z}
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index 808866dea3..292ca7089b 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -1795,7 +1795,6 @@ def _eval_rewrite_as_uppergamma(self, z, **kwargs):
 
     def _eval_nseries(self, x, n, logx, cdir=0):
         # NOTE this is fairly inefficient
-        n += 1
         if self.args[0].subs(x, 0) != 0:
             return super()._eval_nseries(x, n, logx)
         baseseries = self._trigfunc(x)._eval_nseries(x, n, logx)
@@ -1928,11 +1927,11 @@ def _eval_aseries(self, n, args0, x, logx):
         # Expansion at oo
         if point is S.Infinity:
             z = self.args[0]
-            p = [S.NegativeOne**k * factorial(2*k) / z**(2*k)
-                    for k in range(int((n - 1)/2))] + [Order(1/z**n, x)]
-            q = [S.NegativeOne**k * factorial(2*k + 1) / z**(2*k + 1)
-                    for k in range(int(n/2) - 1)] + [Order(1/z**n, x)]
-            return pi/2 - (cos(z)/z)*Add(*p) - (sin(z)/z)*Add(*q)
+            p = [S.NegativeOne**k * factorial(2*k) / z**(2*k + 1)
+                    for k in range(n//2 + 1)] + [Order(1/z**n, x)]
+            q = [S.NegativeOne**k * factorial(2*k + 1) / z**(2*(k + 1))
+                    for k in range(n//2)] + [Order(1/z**n, x)]
+            return pi/2 - cos(z)*Add(*p) - sin(z)*Add(*q)
 
         # All other points are not handled
         return super(Si, self)._eval_aseries(n, args0, x, logx)
@@ -2071,11 +2070,11 @@ def _eval_aseries(self, n, args0, x, logx):
 
         if point in (S.Infinity, S.NegativeInfinity):
             z = self.args[0]
-            p = [S.NegativeOne**k * factorial(2*k) / z**(2*k)
-                    for k in range(int((n - 1)/2))] + [Order(1/z**n, x)]
-            q = [S.NegativeOne**k * factorial(2*k + 1) / z**(2*k + 1)
-                    for k in range(int(n/2) - 1)] + [Order(1/z**n, x)]
-            result = (sin(z)/z)*Add(*p) - (cos(z)/z)*Add(*q)
+            p = [S.NegativeOne**k * factorial(2*k) / z**(2*k + 1)
+                    for k in range(n//2 + 1)] + [Order(1/z**n, x)]
+            q = [S.NegativeOne**k * factorial(2*k + 1) / z**(2*(k + 1))
+                    for k in range(n//2)] + [Order(1/z**n, x)]
+            result = sin(z)*(Add(*p)) - cos(z)*(Add(*q))
 
             if point is S.NegativeInfinity:
                 result += I*pi
diff --git a/sympy/physics/optics/utils.py b/sympy/physics/optics/utils.py
index 881f1cdfc0..4f1f612b21 100644
--- a/sympy/physics/optics/utils.py
+++ b/sympy/physics/optics/utils.py
@@ -671,8 +671,8 @@ def hyperfocal_distance(f, N, c):
 def transverse_magnification(si, so):
     """
 
-    Calculates the transverse magnification upon reflection in a mirror,
-    which is the ratio of the image size to the object size.
+    Calculates the transverse magnification, which is the ratio of the
+    image size to the object size.
 
     Parameters
     ==========
diff --git a/sympy/polys/matrices/sdm.py b/sympy/polys/matrices/sdm.py
index 3fc3643cde..a465fe66b7 100644
--- a/sympy/polys/matrices/sdm.py
+++ b/sympy/polys/matrices/sdm.py
@@ -1805,13 +1805,6 @@ def sdm_rref_den(A, K):
         Aij = Ai[j]
         return ({0: Ai.copy()}, Aij, [j])
 
-    # For inexact domains like RR[x] we use quo and discard the remainder.
-    # Maybe it would be better for K.exquo to do this automatically.
-    if K.is_Exact:
-        exquo = K.exquo
-    else:
-        exquo = K.quo
-
     # Make sure we have the rows in order to make this deterministic from the
     # outset.
     _, rows_in_order = zip(*sorted(A.items()))
@@ -1908,7 +1901,7 @@ def sdm_rref_den(A, K):
                 for l, Akl in Ak.items():
                     Akl = Akl * Aij
                     if divisor is not None:
-                        Akl = exquo(Akl, divisor)
+                        Akl = K.exquo(Akl, divisor)
                     Ak[l] = Akl
                 continue
 
@@ -1919,19 +1912,19 @@ def sdm_rref_den(A, K):
             for l in Ai_nz - Ak_nz:
                 Ak[l] = - Akj * Ai[l]
                 if divisor is not None:
-                    Ak[l] = exquo(Ak[l], divisor)
+                    Ak[l] = K.exquo(Ak[l], divisor)
 
             # This loop also not needed in sdm_irref.
             for l in Ak_nz - Ai_nz:
                 Ak[l] = Aij * Ak[l]
                 if divisor is not None:
-                    Ak[l] = exquo(Ak[l], divisor)
+                    Ak[l] = K.exquo(Ak[l], divisor)
 
             for l in Ai_nz & Ak_nz:
                 Akl = Aij * Ak[l] - Akj * Ai[l]
                 if Akl:
                     if divisor is not None:
-                        Akl = exquo(Akl, divisor)
+                        Akl = K.exquo(Akl, divisor)
                     Ak[l] = Akl
                 else:
                     Ak.pop(l)
@@ -1955,7 +1948,7 @@ def sdm_rref_den(A, K):
                 denom *= Aij
 
         if divisor is not None:
-            denom = exquo(denom, divisor)
+            denom = K.exquo(denom, divisor)
 
         # Update the divisor.
         divisor = denom
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index 215c5e8f3b..cba48dabc7 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -1,5 +1,3 @@
-from __future__ import annotations
-
 from sympy.core.add import Add
 from sympy.core.assumptions import check_assumptions
 from sympy.core.containers import Tuple
@@ -3265,45 +3263,15 @@ def descent(A, B):
     return _remove_gcd(x_0*x_1 + A*z_0*z_1, z_0*x_1 + x_0*z_1, t_1*t_2*y_1)
 
 
-def gaussian_reduce(w:int, a:int, b:int) -> tuple[int, int]:
+def gaussian_reduce(w, a, b):
     r"""
     Returns a reduced solution `(x, z)` to the congruence
-    `X^2 - aZ^2 \equiv 0 \pmod{b}` so that `x^2 + |a|z^2` is minimal.
-    Here ``w`` is a solution of the congruence `x^2 \equiv a \pmod{b}`.
-
-    Explanation
-    ===========
-
-    We define a special dot product of the vectors `u = (u_{1}, u_{2})` and
-    `v = (v_{1}, v_{2})` which is defined in order to reduce solution of
-    the congruence equation `X^2 - aZ^2 \equiv 0 \pmod{b}`.
-
-    .. math ::
-        u \cdot v := (wu_1 + bu_2)(wv_1 + bv_2) + |a|u_1v_1
-
-    Here, if `u_3` is defined as `wu_1 + bu_2`, `u \cdot v` can be computed in `u_3v_3 + |a|u_1v_1`.
-    Therefore, it is kept in pairs of ``(u_3, u_1)``.
+    `X^2 - aZ^2 \equiv 0 \ (mod \ b)` so that `x^2 + |a|z^2` is minimal.
 
-    Parameters
-    ==========
-
-    w : int
-        ``w`` s.t. `w^2 \equiv a \pmod{b}`
-    a : int
-    b : int
+    Details
+    =======
 
-    Examples
-    ========
-
-    >>> from sympy.solvers.diophantine.diophantine import gaussian_reduce
-    >>> from sympy.ntheory.residue_ntheory import sqrt_mod
-    >>> a, b = 19, 101
-    >>> gaussian_reduce(sqrt_mod(a, b), a, b) # 1**2 - 19*(-4)**2 = -303
-    (1, -4)
-    >>> a, b = 11, 14
-    >>> x, z = gaussian_reduce(sqrt_mod(a, b), a, b)
-    >>> (x**2 - a*z**2) % b == 0
-    True
+    Here ``w`` is a solution of the congruence `x^2 \equiv a \ (mod \ b)`
 
     References
     ==========
@@ -3314,25 +3282,48 @@ def gaussian_reduce(w:int, a:int, b:int) -> tuple[int, int]:
            Mathematics of Computation, 72(243), 1417-1441.
            https://doi.org/10.1090/S0025-5718-02-01480-1
     """
-    a = abs(a)
-    def _dot(u, v):
-        return u[0]*v[0] + a*u[1]*v[1]
+    u = (0, 1)
+    v = (1, 0)
 
-    u = (b, 0)
-    v = (w, 1) if b*w >= 0 else (-w, -1)
-    # i.e., _dot(u, v) >= 0
+    if dot(u, v, w, a, b) < 0:
+        v = (-v[0], -v[1])
 
-    if b**2 < w**2 + a:
+    if norm(u, w, a, b) < norm(v, w, a, b):
         u, v = v, u
-    # i.e., norm(u) >= norm(v), where norm(u) := sqrt(_dot(u, u))
-
-    while _dot(u, u) > (dv := _dot(v, v)):
-        k = _dot(u, v) // dv
-        u, v = v, (u[0] - k*v[0], u[1] - k*v[1])
-    c = (v[0] - u[0], v[1] - u[1])
-    if _dot(c, c) <= _dot(u, u) <= 2*_dot(u, v):
-        return c
-    return u
+
+    while norm(u, w, a, b) > norm(v, w, a, b):
+        k = dot(u, v, w, a, b) // dot(v, v, w, a, b)
+        u, v = v, (u[0]- k*v[0], u[1]- k*v[1])
+
+    u, v = v, u
+
+    if dot(u, v, w, a, b) < dot(v, v, w, a, b)/2 or norm((u[0]-v[0], u[1]-v[1]), w, a, b) > norm(v, w, a, b):
+        c = v
+    else:
+        c = (u[0] - v[0], u[1] - v[1])
+
+    return c[0]*w + b*c[1], c[0]
+
+
+def dot(u, v, w, a, b):
+    r"""
+    Returns a special dot product of the vectors `u = (u_{1}, u_{2})` and
+    `v = (v_{1}, v_{2})` which is defined in order to reduce solution of
+    the congruence equation `X^2 - aZ^2 \equiv 0 \ (mod \ b)`.
+    """
+    u_1, u_2 = u
+    v_1, v_2 = v
+    return (w*u_1 + b*u_2)*(w*v_1 + b*v_2) + abs(a)*u_1*v_1
+
+
+def norm(u, w, a, b):
+    r"""
+    Returns the norm of the vector `u = (u_{1}, u_{2})` under the dot product
+    defined by `u \cdot v = (wu_{1} + bu_{2})(w*v_{1} + bv_{2}) + |a|*u_{1}*v_{1}`
+    where `u = (u_{1}, u_{2})` and `v = (v_{1}, v_{2})`.
+    """
+    u_1, u_2 = u
+    return sqrt(dot((u_1, u_2), (u_1, u_2), w, a, b))
 
 
 def holzer(x, y, z, a, b, c):


Test patch:
diff --git a/.github/workflows/runtests.yml b/.github/workflows/runtests.yml
index 1552e84a60..fe46ad79c5 100644
--- a/.github/workflows/runtests.yml
+++ b/.github/workflows/runtests.yml
@@ -29,12 +29,12 @@ jobs:
       - name: Basic code quality tests
         run: bin/test quality
 
-      - name: Run Ruff on the sympy package
-        run: ruff check .
-
       - 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/functions/special/tests/test_error_functions.py b/sympy/functions/special/tests/test_error_functions.py
index 2ad877cc9a..8fb6645761 100644
--- a/sympy/functions/special/tests/test_error_functions.py
+++ b/sympy/functions/special/tests/test_error_functions.py
@@ -621,16 +621,15 @@ def test_si():
         Si(1/x)._eval_as_leading_term(x, cdir=-1) == Si(1/x)
 
     assert Si(x).nseries(x, n=8) == \
-        x - x**3/18 + x**5/600 - x**7/35280 + O(x**9)
+        x - x**3/18 + x**5/600 - x**7/35280 + O(x**8)
     assert Shi(x).nseries(x, n=8) == \
-        x + x**3/18 + x**5/600 + x**7/35280 + O(x**9)
-    assert Si(sin(x)).nseries(x, n=5) == x - 2*x**3/9 + 17*x**5/450 + O(x**6)
+        x + x**3/18 + x**5/600 + x**7/35280 + O(x**8)
+    assert Si(sin(x)).nseries(x, n=5) == x - 2*x**3/9 + O(x**5)
     assert Si(x).nseries(x, 1, n=3) == \
         Si(1) + (x - 1)*sin(1) + (x - 1)**2*(-sin(1)/2 + cos(1)/2) + O((x - 1)**3, (x, 1))
 
-    assert Si(x).series(x, oo) == pi/2 - (- 6/x**3 + 1/x \
-        + O(x**(-7), (x, oo)))*sin(x)/x - (24/x**4 - 2/x**2 + 1 \
-        + O(x**(-7), (x, oo)))*cos(x)/x
+    assert Si(x).series(x, oo) == -sin(x)*(-6/x**4 + x**(-2) + O(x**(-6), (x, oo))) - \
+        cos(x)*(24/x**5 - 2/x**3 + 1/x + O(x**(-6), (x, oo))) + pi/2
 
     t = Symbol('t', Dummy=True)
     assert Si(x).rewrite(sinc).dummy_eq(Integral(sinc(t), (t, 0, x)))
@@ -675,16 +674,15 @@ def test_ci():
     assert tn_arg(Chi)
 
     assert Ci(x).nseries(x, n=4) == \
-        EulerGamma + log(x) - x**2/4 + x**4/96 + O(x**5)
+        EulerGamma + log(x) - x**2/4 + O(x**4)
     assert Chi(x).nseries(x, n=4) == \
-        EulerGamma + log(x) + x**2/4 + x**4/96 + O(x**5)
+        EulerGamma + log(x) + x**2/4 + O(x**4)
 
-    assert Ci(x).series(x, oo) == -cos(x)*(-6/x**3 + 1/x \
-        + O(x**(-7), (x, oo)))/x + (24/x**4 - 2/x**2 + 1 \
-        + O(x**(-7), (x, oo)))*sin(x)/x
+    assert Ci(x).series(x, oo) == -cos(x)*(-6/x**4 + x**(-2) + O(x**(-6), (x, oo))) + \
+        sin(x)*(24/x**5 - 2/x**3 + 1/x + O(x**(-6), (x, oo)))
 
-    assert Ci(x).series(x, -oo) == -cos(x)*(-6/x**3 + 1/x + O(-1/x**7, (x, -oo)))/x + \
-            (24/x**4 - 2/x**2 + 1 + O(-1/x**7, (x, -oo)))*sin(x)/x + I*pi
+    assert Ci(x).series(x, -oo) == -cos(x)*(-6/x**4 + x**(-2) + O(x**(-6), (x, -oo))) + \
+        sin(x)*(24/x**5 - 2/x**3 + 1/x + O(x**(-6), (x, -oo))) + I*pi
 
     assert limit(log(x) - Ci(2*x), x, 0) == -log(2) - EulerGamma
     assert Ci(x).rewrite(uppergamma) == -expint(1, x*exp_polar(-I*pi/2))/2 -\
diff --git a/sympy/polys/matrices/tests/test_inverse.py b/sympy/polys/matrices/tests/test_inverse.py
index 47c8279932..9b3415469e 100644
--- a/sympy/polys/matrices/tests/test_inverse.py
+++ b/sympy/polys/matrices/tests/test_inverse.py
@@ -1,4 +1,4 @@
-from sympy import ZZ, Matrix
+from sympy import ZZ
 from sympy.polys.matrices import DM, DomainMatrix
 from sympy.polys.matrices.dense import ddm_iinv
 from sympy.polys.matrices.exceptions import DMNonInvertibleMatrixError
@@ -6,9 +6,6 @@
 
 import pytest
 from sympy.testing.pytest import raises
-from sympy.core.numbers import all_close
-
-from sympy.abc import x
 
 
 # Examples are given as adjugate matrix and determinant adj_det should match
@@ -154,40 +151,3 @@ def test_Matrix_adjugate(name, A, A_inv, den):
 @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
 def test_dm_adj_det(name, A, A_inv, den):
     assert A.adj_det() == (A_inv, den)
-
-
-def test_inverse_inexact():
-
-    M = Matrix([[x-0.3, -0.06, -0.22],
-                [-0.46, x-0.48, -0.41],
-                [-0.14, -0.39, x-0.64]])
-
-    Mn = Matrix([[1.0*x**2 - 1.12*x + 0.1473, 0.06*x + 0.0474, 0.22*x - 0.081],
-                 [0.46*x - 0.237, 1.0*x**2 - 0.94*x + 0.1612, 0.41*x - 0.0218],
-                 [0.14*x + 0.1122, 0.39*x - 0.1086, 1.0*x**2 - 0.78*x + 0.1164]])
-
-    d = 1.0*x**3 - 1.42*x**2 + 0.4249*x - 0.0546540000000002
-
-    Mi = Mn / d
-
-    M_dm = M.to_DM()
-    M_dmd = M_dm.to_dense()
-    M_dm_num, M_dm_den = M_dm.inv_den()
-    M_dmd_num, M_dmd_den = M_dmd.inv_den()
-
-    # XXX: We don't check M_dm().to_field().inv() which currently uses division
-    # and produces a more complicate result from gcd cancellation failing.
-    # DomainMatrix.inv() over RR(x) should be changed to clear denominators and
-    # use DomainMatrix.inv_den().
-
-    Minvs = [
-        M.inv(),
-        (M_dm_num.to_field() / M_dm_den).to_Matrix(),
-        (M_dmd_num.to_field() / M_dmd_den).to_Matrix(),
-        M_dm_num.to_Matrix() / M_dm_den.as_expr(),
-        M_dmd_num.to_Matrix() / M_dmd_den.as_expr(),
-    ]
-
-    for Minv in Minvs:
-        for Mi1, Mi2 in zip(Minv.flat(), Mi.flat()):
-            assert all_close(Mi2, Mi1)