Repository: sympy/sympy

Issue:
solveset issue with floats/trig
from https://groups.google.com/g/sympy/c/8sZ_pbh3QPQ
```
solveset(Eq(sin(x),1/2),x,S.Reals)
∅

solveset(Eq(sin(x),1/2),x,S.Complexes)
{−i(i(2nπ+0.523598775598299)−1.11022302462516⋅10−16)∣∣n∈Z}∪{−i(i(2nπ−0.523598775598299+π)−1.11022302462516⋅10−16)∣∣n∈Z}
```
Oscar writes:
```
The correct answer is given if you use exact rational numbers (S(1)/2
or Rational(1, 2)) rather than the float 1/2:

In [6]: solveset(Eq(sin(x), S(1)/2), x, Reals)
Out[6]:
⎧ 5⋅π │ ⎫ ⎧ π │ ⎫
⎨2⋅n⋅π + ─── │ n ∊ ℤ⎬ ∪ ⎨2⋅n⋅π + ─ │ n ∊ ℤ⎬
⎩ 6 │ ⎭ ⎩ 6 │ ⎭

The problem is that solveset solves this in an overly complicated way
using complex numbers even if the domain is reals. Then a small
rounding error makes it look like the solutions are not real.
```

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).



* No release notes entry will be added for this pull request.

<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
    <!-- 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. -->

    closes #25768

    #### Brief description of what is fixed or changed

    Now detect when a single trig function is present that can be inverted in the real domain.

    #### 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 -->
    NO ENTRY
    <!-- END RELEASE NOTES -->


</details><p>

smichr: @MaqOwais , do you have an idea of what to do in the area that you worked on before in #21212? Not all values in the dictionary are ImageSets so the XFAIL test fails because it doesn't handle those properly (formerly) and had non-Expr in the `rep` dictionary.
github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
| Change   | Before [8059df73] <sympy-1.12^0>   | After [368eb0f4]    |   Ratio | Benchmark (Parameter)                                                |
|----------|------------------------------------|---------------------|---------|----------------------------------------------------------------------|
| -        | 111±3ms                            | 72.4±1ms            |    0.65 | integrate.TimeIntegrationRisch02.time_doit_risch(10)                 |
| +        | 25.9±0.1μs                         | 42.2±0.8μs          |    1.63 | integrate.TimeIntegrationRisch03.time_doit(1)                        |
| -        | 8.47±0.04ms                        | 4.65±0.06ms         |    0.55 | logic.LogicSuite.time_load_file                                      |
| -        | 2.55±0.02ms                        | 792±3μs             |    0.31 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(3, 'sparse')          |
| -        | 12.3±0.1ms                         | 2.38±0.01ms         |    0.19 | polys.TimePREM_LinearDenseQuadraticGCD.time_op(5, 'sparse')          |
| -        | 459±4μs                            | 99.6±2μs            |    0.22 | polys.TimePREM_QuadraticNonMonicGCD.time_op(1, 'sparse')             |
| -        | 5.81±0.1ms                         | 440±1μs             |    0.08 | polys.TimePREM_QuadraticNonMonicGCD.time_op(3, 'sparse')             |
| -        | 12.7±0.2ms                         | 1.33±0ms            |    0.11 | polys.TimePREM_QuadraticNonMonicGCD.time_op(5, 'sparse')             |
| -        | 7.86±0.07ms                        | 4.83±0.07ms         |    0.61 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(2, 'sparse') |
| -        | 33.4±0.2ms                         | 14.6±0.04ms         |    0.44 | polys.TimeSUBRESULTANTS_LinearDenseQuadraticGCD.time_op(3, 'sparse') |
| -        | 8.74±0.09ms                        | 1.42±0.01ms         |    0.16 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(1, 'sparse')    |
| -        | 19.9±0.3ms                         | 11.3±0.04ms         |    0.57 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(2, 'sparse')    |
| -        | 250±5ms                            | 85.7±0.09ms         |    0.34 | polys.TimeSUBRESULTANTS_QuadraticNonMonicGCD.time_op(3, 'sparse')    |
| -        | 7.74±0.1ms                         | 644±2μs             |    0.08 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(3, 'sparse')     |
| -        | 33.5±0.5ms                         | 1.04±0ms            |    0.03 | polys.TimeSUBRESULTANTS_SparseGCDHighDegree.time_op(5, 'sparse')     |
| -        | 772±5μs                            | 239±1μs             |    0.31 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(1, 'sparse') |
| -        | 7.91±0.06ms                        | 240±0.6μs           |    0.03 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(3, 'sparse') |
| -        | 19.5±0.1ms                         | 246±1μs             |    0.01 | polys.TimeSUBRESULTANTS_SparseNonMonicQuadratic.time_op(5, 'sparse') |
| -        | 195±5μs                            | 106±0.3μs           |    0.54 | solve.TimeMatrixOperations.time_rref(3, 0)                           |
| -        | 370±2μs                            | 130±1μs             |    0.35 | solve.TimeMatrixOperations.time_rref(4, 0)                           |
| -        | 38.2±0.07ms                        | 16.7±0.3ms          |    0.44 | solve.TimeSolveLinSys189x49.time_solve_lin_sys                       |

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

Patch:
diff --git a/.mailmap b/.mailmap
index 1f673809c7..2475033668 100644
--- a/.mailmap
+++ b/.mailmap
@@ -466,7 +466,6 @@ Colin B. Macdonald <cbm@m.fsf.org> <macdonald@maths.ox.ac.uk>
 Colin Marquardt <github@marquardt-home.de>
 Colleen Lee <colleenclee@gmail.com> <clee@coursera.org>
 Comer Duncan <comer.duncan@gmail.com>
-Congxu Yang <u7189828@anu.edu.au>
 Constantin Mateescu <costica1234@me.com>
 Costor <pcs2009@web.de>
 Craig A. Stoudt <craig.stoudt@gmail.com>
@@ -1051,7 +1050,6 @@ Oscar Benjamin <oscar.j.benjamin@gmail.com> <oscar@kar-wench.(none)>
 Oscar Gerardo Lazo Arjona <oscar.lazoarjona@physics.ox.ac.uk> oscarlazoarjona <oscar.lazoarjona@physics.ox.ac.uk>
 Oscar Gustafsson <oscar.gustafsson@gmail.com>
 P. Sai Prasanth <psai.prasanth.min16@itbhu.ac.in>
-Pablo <48098178+PabloRuizCuevas@users.noreply.github.com>
 Pablo Galindo Salgado <pablogsal@gmail.com> Pablo Galindo <pablogsal@gmail.com>
 Pablo Galindo Salgado <pablogsal@gmail.com> Pablo Galindo Salgado <Pablogsal@gmail.com>
 Pablo Galindo Salgado <pablogsal@gmail.com> pablogsal <Pablogsal@gmail.com>
diff --git a/doc/src/explanation/active-deprecations.md b/doc/src/explanation/active-deprecations.md
index eea06df81b..794b91b52c 100644
--- a/doc/src/explanation/active-deprecations.md
+++ b/doc/src/explanation/active-deprecations.md
@@ -265,30 +265,8 @@ Now they should be imported from ``sympy.physics.mechanics``:
 >>> from sympy.physics.mechanics.loads import gravity
 ```
 
-(modularinteger-compare)=
-### Ordered comparisons like ``a < b`` with modular integers
 
-SymPy's ``GF`` domains represent modular integers. Previously it was possible
-to compare these with ordered comparisons like ``a < b``:
-```py
->>> from sympy import GF
->>> F5 = GF(5)
->>> F5(2) < F5(3) # doctest: +SKIP
-True
-```
-This will now fail with ``TypeError`` when the ground types are set to
-``flint``. When the ground types are not ``flint`` these comparisons are now
-deprecated: they will still work but will give a deprecation warning when used.
-
-Ordered comparisons of modular integer or finite fields do not make sense
-because these are not ordered fields:
-```
->>> e = F5(4)
->>> e + 1 > e # doctest: +SKIP
-False
-```
-
-(modularinteger-to-int)=
+modularinteger-to-int=
 ### The ``ModularInteger.to_int()`` method
 
 SymPy's ``GF`` domains are for modular integers e.g. ``GF(n)`` is for the
diff --git a/sympy/polys/densetools.py b/sympy/polys/densetools.py
index b56a4a773b..ff3ea5f2ea 100644
--- a/sympy/polys/densetools.py
+++ b/sympy/polys/densetools.py
@@ -453,10 +453,6 @@ def dup_trunc(f, p, K):
                 g.append(c - p)
             else:
                 g.append(c)
-    elif K.is_FiniteField:
-        # XXX: python-flint's nmod does not support %
-        pi = int(p)
-        g = [ K(int(c) % pi) for c in f ]
     else:
         g = [ c % p for c in f ]
 
diff --git a/sympy/polys/domains/finitefield.py b/sympy/polys/domains/finitefield.py
index 9569e7f92c..27801f00b4 100644
--- a/sympy/polys/domains/finitefield.py
+++ b/sympy/polys/domains/finitefield.py
@@ -1,6 +1,5 @@
 """Implementation of :class:`FiniteField` class. """
 
-from sympy.external.gmpy import GROUND_TYPES
 
 from sympy.core.numbers import int_valued
 from sympy.polys.domains.field import Field
@@ -12,43 +11,6 @@
 from sympy.utilities import public
 from sympy.polys.domains.groundtypes import SymPyInteger
 
-
-if GROUND_TYPES == 'flint':
-    import flint
-    # Don't use python-flint < 0.5.0 because nmod was missing some features in
-    # previous versions of python-flint and fmpz_mod was not yet added.
-    _major, _minor, *_ = flint.__version__.split('.')
-    if (int(_major), int(_minor)) < (0, 5):
-        flint = None
-else:
-    flint = None
-
-
-def _modular_int_factory(mod, dom, symmetric, self):
-
-    # Use flint if available
-    if flint is not None:
-        try:
-            mod = dom.convert(mod)
-        except CoercionFailed:
-            raise ValueError('modulus must be an integer, got %s' % mod)
-
-        # flint's nmod is only for moduli up to 2^64-1 (on a 64-bit machine)
-        try:
-            flint.nmod(0, mod)
-        except OverflowError:
-            # Use fmpz_mod
-            ctx = flint.fmpz_mod_ctx(mod)
-        else:
-            # Use nmod
-            ctx = lambda x: flint.nmod(x, mod)
-
-        return ctx
-
-    # Use the Python implementation
-    return ModularIntegerFactory(mod, dom, symmetric, self)
-
-
 @public
 class FiniteField(Field, SimpleDomain):
     r"""Finite field of prime order :ref:`GF(p)`
@@ -158,17 +120,12 @@ def __init__(self, mod, symmetric=True):
         if mod <= 0:
             raise ValueError('modulus must be a positive integer, got %s' % mod)
 
-        self.dtype = _modular_int_factory(mod, dom, symmetric, self)
+        self.dtype = ModularIntegerFactory(mod, dom, symmetric, self)
         self.zero = self.dtype(0)
         self.one = self.dtype(1)
         self.dom = dom
         self.mod = mod
         self.sym = symmetric
-        self._tp = type(self.zero)
-
-    @property
-    def tp(self):
-        return self._tp
 
     def __str__(self):
         return 'GF(%s)' % self.mod
@@ -204,34 +161,18 @@ def from_sympy(self, a):
 
     def to_int(self, a):
         """Convert ``val`` to a Python ``int`` object. """
-        aval = int(a)
+        aval = a.val
         if self.sym and aval > self.mod // 2:
             aval -= self.mod
         return aval
 
-    def is_positive(self, a):
-        """Returns True if ``a`` is positive. """
-        return bool(a)
-
-    def is_nonnegative(self, a):
-        """Returns True if ``a`` is non-negative. """
-        return True
-
-    def is_negative(self, a):
-        """Returns True if ``a`` is negative. """
-        return False
-
-    def is_nonpositive(self, a):
-        """Returns True if ``a`` is non-positive. """
-        return not a
-
     def from_FF(K1, a, K0=None):
         """Convert ``ModularInteger(int)`` to ``dtype``. """
-        return K1.dtype(K1.dom.from_ZZ(int(a), K0.dom))
+        return K1.dtype(K1.dom.from_ZZ(a.val, K0.dom))
 
     def from_FF_python(K1, a, K0=None):
         """Convert ``ModularInteger(int)`` to ``dtype``. """
-        return K1.dtype(K1.dom.from_ZZ_python(int(a), K0.dom))
+        return K1.dtype(K1.dom.from_ZZ_python(a.val, K0.dom))
 
     def from_ZZ(K1, a, K0=None):
         """Convert Python's ``int`` to ``dtype``. """
@@ -274,7 +215,7 @@ def from_RealField(K1, a, K0):
     def is_square(self, a):
         """Returns True if ``a`` is a quadratic residue modulo p. """
         # a is not a square <=> x**2-a is irreducible
-        poly = [int(x) for x in [self.one, self.zero, -a]]
+        poly = [x.val for x in [self.one, self.zero, -a]]
         return not gf_irred_p_rabin(poly, self.mod, self.dom)
 
     def exsqrt(self, a):
@@ -288,7 +229,7 @@ def exsqrt(self, a):
         if self.mod == 2 or a == 0:
             return a
         # Otherwise, use square-free factorization routine to factorize x**2-a
-        poly = [int(x) for x in [self.one, self.zero, -a]]
+        poly = [x.val for x in [self.one, self.zero, -a]]
         for factor in gf_zassenhaus(poly, self.mod, self.dom):
             if len(factor) == 2 and factor[1] <= self.mod // 2:
                 return self.dtype(factor[1])
diff --git a/sympy/polys/domains/modularinteger.py b/sympy/polys/domains/modularinteger.py
index 98ff7b80f4..42bb09fb9a 100644
--- a/sympy/polys/domains/modularinteger.py
+++ b/sympy/polys/domains/modularinteger.py
@@ -150,29 +150,11 @@ def __pow__(self, exp):
     def _compare(self, other, op):
         val = self._get_val(other)
 
-        if val is None:
-            return NotImplemented
-
-        return op(self.val, val % self.mod)
-
-    def _compare_deprecated(self, other, op):
-        val = self._get_val(other)
-
-        if val is None:
+        if val is not None:
+            return op(self.val, val % self.mod)
+        else:
             return NotImplemented
 
-        sympy_deprecation_warning(
-            """Ordered comparisons with modular integers are deprecated.
-
-            Use e.g. int(a) < int(b) instead of a < b.
-            """,
-            deprecated_since_version="1.13",
-            active_deprecations_target="modularinteger-compare",
-            stacklevel=4,
-        )
-
-        return op(self.val, val % self.mod)
-
     def __eq__(self, other):
         return self._compare(other, operator.eq)
 
@@ -180,16 +162,16 @@ def __ne__(self, other):
         return self._compare(other, operator.ne)
 
     def __lt__(self, other):
-        return self._compare_deprecated(other, operator.lt)
+        return self._compare(other, operator.lt)
 
     def __le__(self, other):
-        return self._compare_deprecated(other, operator.le)
+        return self._compare(other, operator.le)
 
     def __gt__(self, other):
-        return self._compare_deprecated(other, operator.gt)
+        return self._compare(other, operator.gt)
 
     def __ge__(self, other):
-        return self._compare_deprecated(other, operator.ge)
+        return self._compare(other, operator.ge)
 
     def __bool__(self):
         return bool(self.val)
diff --git a/sympy/polys/polyutils.py b/sympy/polys/polyutils.py
index b5c49fcc9d..d754d9c413 100644
--- a/sympy/polys/polyutils.py
+++ b/sympy/polys/polyutils.py
@@ -1,8 +1,5 @@
 """Useful utilities for higher level polynomial classes. """
 
-from __future__ import annotations
-
-from sympy.external.gmpy import GROUND_TYPES
 
 from sympy.core import (S, Add, Mul, Pow, Eq, Expr,
     expand_mul, expand_multinomial)
@@ -11,8 +8,8 @@
 from sympy.polys.polyerrors import PolynomialError, GeneratorsError
 from sympy.polys.polyoptions import build_options
 
-import re
 
+import re
 
 _gens_order = {
     'a': 301, 'b': 302, 'c': 303, 'd': 304,
@@ -157,35 +154,20 @@ def _analyze_gens(gens):
 
 def _sort_factors(factors, **args):
     """Sort low-level factors in increasing 'complexity' order. """
-
-    # XXX: GF(p) does not support comparisons so we need a key function to sort
-    # the factors if python-flint is being used. A better solution might be to
-    # add a sort key method to each domain.
-    def order_key(factor):
-        if isinstance(factor, _GF_types):
-            return int(factor)
-        elif isinstance(factor, list):
-            return [order_key(f) for f in factor]
-        else:
-            return factor
-
     def order_if_multiple_key(factor):
         (f, n) = factor
-        return (len(f), n, order_key(f))
+        return (len(f), n, f)
 
     def order_no_multiple_key(f):
-        return (len(f), order_key(f))
+        return (len(f), f)
 
     if args.get('multiple', True):
         return sorted(factors, key=order_if_multiple_key)
     else:
         return sorted(factors, key=order_no_multiple_key)
 
-
 illegal_types = [type(obj) for obj in _illegal]
 finf = [float(i) for i in _illegal[1:3]]
-
-
 def _not_a_coeff(expr):
     """Do not treat NaN and infinities as valid polynomial coefficients. """
     if type(expr) in illegal_types or expr in finf:
@@ -568,15 +550,3 @@ def _zeroth_power(self):
     def _first_power(self):
         """Return a copy of self."""
         raise NotImplementedError
-
-
-_GF_types: tuple[type, ...]
-
-
-if GROUND_TYPES == 'flint':
-    import flint
-    _GF_types = (flint.nmod, flint.fmpz_mod)
-else:
-    from sympy.polys.domains.modularinteger import ModularInteger
-    flint = None
-    _GF_types = (ModularInteger,)
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index f5650e1116..985450284f 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -58,6 +58,7 @@
                                        is_sequence, iterable)
 from sympy.calculus.util import periodicity, continuous_domain, function_range
 
+
 from types import GeneratorType
 
 
@@ -313,12 +314,14 @@ def _invert_real(f, g_ys, symbol):
             def inv(trig):
                 if isinstance(trig, (sin, csc)):
                     F = asin if isinstance(trig, sin) else acsc
-                    return (lambda a: n*pi + S.NegativeOne**n*F(a),)
+                    return (
+                        lambda a: 2*n*pi + F(a),
+                        lambda a: 2*n*pi + pi - F(a))
                 if isinstance(trig, (cos, sec)):
                     F = acos if isinstance(trig, cos) else asec
                     return (
                         lambda a: 2*n*pi + F(a),
-                        lambda a: 2*n*pi - F(a),)
+                        lambda a: 2*n*pi - F(a))
                 if isinstance(trig, (tan, cot)):
                     return (lambda a: n*pi + trig.inverse()(a),)
 
@@ -373,7 +376,7 @@ def _invert_complex(f, g_ys, symbol):
 
     if isinstance(f, exp) or (f.is_Pow and f.base == S.Exp1):
         if isinstance(g_ys, ImageSet):
-            # can solve upto `(d*exp(exp(...(exp(a*x + b))...) + c)` format.
+            # can solve up to `(d*exp(exp(...(exp(a*x + b))...) + c)` format.
             # Further can be improved to `(d*exp(exp(...(exp(a*x**n + b*x**(n-1) + ... + f))...) + c)`.
             g_ys_expr = g_ys.lamda.expr
             g_ys_vars = g_ys.lamda.variables
@@ -1020,6 +1023,7 @@ def _solveset(f, symbol, domain, _check=False):
         return domain
 
     orig_f = f
+    invert_trig = False  # True if we will use inversion to solve
     if f.is_Mul:
         coeff, f = f.as_independent(symbol, as_Add=False)
         if coeff in {S.ComplexInfinity, S.NegativeInfinity, S.Infinity}:
@@ -1030,6 +1034,10 @@ def _solveset(f, symbol, domain, _check=False):
         if m not in {S.ComplexInfinity, S.Zero, S.Infinity,
                               S.NegativeInfinity}:
             f = a/m + h  # XXX condition `m != 0` should be added to soln
+        if isinstance(h, TrigonometricFunction) and (a and a.is_number
+                and a.is_real and domain.is_subset(S.Reals)):
+            # solve this by inversion
+            invert_trig = True
 
     # assign the solvers to use
     solver = lambda f, x, domain=domain: _solveset(f, x, domain)
@@ -1050,8 +1058,8 @@ def _solveset(f, symbol, domain, _check=False):
         # wrong solutions we are using this technique only if both f and g are
         # finite for a finite input.
         result = Union(*[solver(m, symbol) for m in f.args])
-    elif _is_function_class_equation(TrigonometricFunction, f, symbol) or \
-            _is_function_class_equation(HyperbolicFunction, f, symbol):
+    elif not invert_trig and (_is_function_class_equation(TrigonometricFunction, f, symbol) or \
+            _is_function_class_equation(HyperbolicFunction, f, symbol)):
         result = _solve_trig(f, symbol, domain)
     elif isinstance(f, arg):
         a = f.args[0]
@@ -2605,11 +2613,6 @@ def linear_eq_to_matrix(equations, *symbols):
             are to be found.
             '''))
 
-    # Check if 'symbols' is a set and raise an error if it is
-    if isinstance(symbols[0], set):
-        raise TypeError(
-            "Unordered 'set' type is not supported as input for symbols.")
-
     if hasattr(symbols[0], '__iter__'):
         symbols = symbols[0]
 
@@ -3065,7 +3068,6 @@ def _unsolved_syms(eq, sort=False):
             unsolved = list(unsolved)
             unsolved.sort(key=default_sort_key)
         return unsolved
-    # end of _unsolved_syms()
 
     # sort such that equation with the fewest potential symbols is first.
     # means eq with less number of variable first in the list.
@@ -3103,13 +3105,12 @@ def add_intersection_complement(result, intersection_dict, complement_dict):
             if res_copy is not None:
                 final_result.append(res_copy)
         return final_result
-    # end of def add_intersection_complement()
 
     def _extract_main_soln(sym, sol, soln_imageset):
         """Separate the Complements, Intersections, ImageSet lambda expr and
-        its base_set. This function returns the unmasks sol from different classes
+        its base_set. This function returns the unmasked sol from different classes
         of sets and also returns the appended ImageSet elements in a
-        soln_imageset (dict: where key as unmasked element and value as ImageSet).
+        soln_imageset dict: `{unmasked element: ImageSet}`.
         """
         # if there is union, then need to check
         # Complement, Intersection, Imageset.
@@ -3159,9 +3160,7 @@ def _extract_main_soln(sym, sol, soln_imageset):
         if not isinstance(sol, FiniteSet):
             sol = FiniteSet(sol)
         return sol, soln_imageset
-    # end of def _extract_main_soln()
 
-    # helper function for _append_new_soln
     def _check_exclude(rnew, imgset_yes):
         rnew_ = rnew
         if imgset_yes:
@@ -3184,9 +3183,7 @@ def _check_exclude(rnew, imgset_yes):
         except TypeError:
             satisfy_exclude = None
         return satisfy_exclude
-    # end of def _check_exclude()
 
-    # helper function for _append_new_soln
     def _restore_imgset(rnew, original_imageset, newresult):
         restore_sym = set(rnew.keys()) & \
             set(original_imageset.keys())
@@ -3195,7 +3192,6 @@ def _restore_imgset(rnew, original_imageset, newresult):
             rnew[key_sym] = img
         if rnew not in newresult:
             newresult.append(rnew)
-    # end of def _restore_imgset()
 
     def _append_eq(eq, result, res, delete_soln, n=None):
         u = Dummy('u')
@@ -3259,7 +3255,6 @@ def _append_new_soln(rnew, sym, sol, imgset_yes, soln_imageset,
             rnew = {}
         _restore_imgset(rnew, original_imageset, newresult)
         return newresult, delete_soln
-    # end of def _append_new_soln()
 
     def _new_order_result(result, eq):
         # separate first, second priority. `res` that makes `eq` value equals
@@ -3287,24 +3282,24 @@ def _solve_using_known_values(result, solver):
         total_solvest_call = 0
         total_conditionst = 0
 
-        # sort such that equation with the fewest potential symbols is first.
-        # means eq with less variable first
+        # sort equations so the one with the fewest potential
+        # symbols appears first
         for index, eq in enumerate(eqs_in_better_order):
             newresult = []
             original_imageset = {}
-            # if imageset expr is used to solve other symbol
+            # if imageset, expr is used to solve for other symbol
             imgset_yes = False
-            result = _new_order_result(result, eq)
             for res in result:
                 got_symbol = set()  # symbols solved in one iteration
                 # find the imageset and use its expr.
-                for key_res, value_res in res.items():
-                    if isinstance(value_res, ImageSet):
-                        res[key_res] = value_res.lamda.expr
-                        original_imageset[key_res] = value_res
-                        dummy_n = value_res.lamda.expr.atoms(Dummy).pop()
-                        (base,) = value_res.base_sets
+                for k, v in res.items():
+                    if isinstance(v, ImageSet):
+                        res[k] = v.lamda.expr
+                        original_imageset[k] = v
+                        dummy_n = v.lamda.expr.atoms(Dummy).pop()
+                        (base,) = v.base_sets
                         imgset_yes = (dummy_n, base)
+                    assert not isinstance(v, FiniteSet)  # if so, internal error
                 # update eq with everything that is known so far
                 eq2 = eq.subs(res).expand()
                 unsolved_syms = _unsolved_syms(eq2, sort=True)
@@ -3316,7 +3311,7 @@ def _solve_using_known_values(result, solver):
                         if delete_res:
                             # `delete_res` is true, means substituting `res` in
                             # eq2 doesn't return `zero` or deleting the `res`
-                            # (a soln) since it staisfies expr of `exclude`
+                            # (a soln) since it satisfies expr of `exclude`
                             # list.
                             result.remove(res)
                     continue  # skip as it's independent of desired symbols
@@ -3374,7 +3369,7 @@ def _solve_using_known_values(result, solver):
                         # before this loop
                         sol, soln_imageset = _extract_main_soln(
                             sym, sol, soln_imageset)
-                        sol = set(sol).pop()
+                        sol = set(sol).pop()  # XXX what if there are more solutions?
                         free = sol.free_symbols
                         if got_symbol and any(
                             ss in free for ss in got_symbol
@@ -3404,7 +3399,7 @@ def _solve_using_known_values(result, solver):
                             rnew, sym, sol, imgset_yes, soln_imageset,
                             original_imageset, newresult)
                         if delete_res:
-                            # deleting the `res` (a soln) since it staisfies
+                            # deleting the `res` (a soln) since it satisfies
                             # eq of `exclude` list
                             result.remove(res)
                     # solution got for sym
@@ -3414,7 +3409,6 @@ def _solve_using_known_values(result, solver):
             if newresult:
                 result = newresult
         return result, total_solvest_call, total_conditionst
-    # end def _solve_using_know_values()
 
     new_result_real, solve_call1, cnd_call1 = _solve_using_known_values(
         old_result, solveset_real)
@@ -3478,7 +3472,6 @@ def _solve_using_known_values(result, solver):
         temp = [r[symb] for symb in all_symbols]
         result += FiniteSet(tuple(temp))
     return result
-# end of def substitution()
 
 
 def _solveset_work(system, symbols):
@@ -3504,7 +3497,6 @@ def _handle_positive_dimensional(polys, symbols, denominators):
         new_system, symbols, result, [],
         denominators)
     return result
-# end of def _handle_positive_dimensional()
 
 
 def _handle_zero_dimensional(polys, symbols, system):
@@ -3519,7 +3511,6 @@ def _handle_zero_dimensional(polys, symbols, system):
         if all(checksol(eq, dict_sym_value) for eq in system):
             result_update += FiniteSet(res)
     return result_update
-# end of def _handle_zero_dimensional()
 
 
 def _separate_poly_nonpoly(system, symbols):
@@ -3556,7 +3547,6 @@ def _separate_poly_nonpoly(system, symbols):
         else:
             nonpolys.append(eq)
     return polys, polys_expr, nonpolys, denominators, unrad_changed
-# end of def _separate_poly_nonpoly()
 
 
 def _handle_poly(polys, symbols):
diff --git a/sympy/stats/frv_types.py b/sympy/stats/frv_types.py
index f0a392a02e..70795a76aa 100644
--- a/sympy/stats/frv_types.py
+++ b/sympy/stats/frv_types.py
@@ -561,11 +561,11 @@ class HypergeometricDistribution(SingleFiniteDistribution):
     @staticmethod
     def check(n, N, m):
         _value_check((N.is_integer, N.is_nonnegative),
-                     "'N' must be nonnegative integer. N = %s." % str(N))
+                     "'N' must be nonnegative integer. N = %s." % str(n))
         _value_check((n.is_integer, n.is_nonnegative),
                      "'n' must be nonnegative integer. n = %s." % str(n))
         _value_check((m.is_integer, m.is_nonnegative),
-                     "'m' must be nonnegative integer. m = %s." % str(m))
+                     "'m' must be nonnegative integer. m = %s." % str(n))
 
     @property
     def is_symbolic(self):


Test patch:
diff --git a/sympy/ntheory/tests/test_factor_.py b/sympy/ntheory/tests/test_factor_.py
index 9ff6fbb290..959dc6fd96 100644
--- a/sympy/ntheory/tests/test_factor_.py
+++ b/sympy/ntheory/tests/test_factor_.py
@@ -6,7 +6,7 @@
 from sympy.core.symbol import Symbol
 from sympy.functions.combinatorial.factorials import factorial as fac
 from sympy.core.numbers import Integer, Rational
-from sympy.external.gmpy import gcd
+from sympy.external.gmpy import gcd, GROUND_TYPES, SYMPY_INTS
 
 from sympy.ntheory import (totient,
     factorint, primefactors, divisors, nextprime,
@@ -301,6 +301,18 @@ def test_factorint():
     assert str(factorint(Dict(n))) == sans
 
 
+def test_factorint_flint_int_issue():
+    """
+    Because of a bug in flint, we are temporarily wrapping int in `factorint`.
+    When this fails, the int wrapper is no longer needed in factorint.
+
+    https://github.com/flintlib/python-flint/issues/92
+    https://github.com/sympy/sympy/pull/25749
+    """
+    if GROUND_TYPES == 'flint':
+        assert pow(2, 5, SYMPY_INTS[1](1000)) == 1
+
+
 def test_divisors_and_divisor_count():
     assert divisors(-1) == [1]
     assert divisors(0) == []
diff --git a/sympy/polys/domains/tests/test_domains.py b/sympy/polys/domains/tests/test_domains.py
index a569f44f65..f220eca753 100644
--- a/sympy/polys/domains/tests/test_domains.py
+++ b/sympy/polys/domains/tests/test_domains.py
@@ -1,7 +1,5 @@
 """Tests for classes defining properties of ground domains, e.g. ZZ, QQ, ZZ[x] ... """
 
-from sympy.external.gmpy import GROUND_TYPES
-
 from sympy.core.numbers import (AlgebraicNumber, E, Float, I, Integer,
     Rational, oo, pi, _illegal)
 from sympy.core.singleton import S
@@ -11,6 +9,8 @@
 from sympy.polys.polytools import Poly
 from sympy.abc import x, y, z
 
+from sympy.external.gmpy import GROUND_TYPES
+
 from sympy.polys.domains import (ZZ, QQ, RR, CC, FF, GF, EX, EXRAW, ZZ_gmpy,
     ZZ_python, QQ_gmpy, QQ_python)
 from sympy.polys.domains.algebraicfield import AlgebraicField
@@ -32,7 +32,7 @@
     NotInvertible,
     DomainError)
 
-from sympy.testing.pytest import raises, warns_deprecated_sympy
+from sympy.testing.pytest import raises
 
 from itertools import product
 
@@ -811,7 +811,6 @@ def test_PolynomialRing_from_FractionField():
     assert R.to_domain().from_FractionField(g, F.to_domain()) == X**2/4 + Y**2/4
     assert R.to_domain().from_FractionField(h, F.to_domain()) == X**2 + Y**2
 
-
 def test_FractionField_from_PolynomialRing():
     R, x,y = ring("x,y", QQ)
     F, X,Y = field("x,y", ZZ)
@@ -822,12 +821,10 @@ def test_FractionField_from_PolynomialRing():
     assert F.to_domain().from_PolynomialRing(f, R.to_domain()) == 3*X**2 + 5*Y**2
     assert F.to_domain().from_PolynomialRing(g, R.to_domain()) == (5*X**2 + 3*Y**2)/15
 
-
 def test_FF_of_type():
-    # XXX: of_type is not very useful here because in the case of ground types
-    # = flint all elements are of type nmod.
     assert FF(3).of_type(FF(3)(1)) is True
     assert FF(5).of_type(FF(5)(3)) is True
+    assert FF(5).of_type(FF(7)(3)) is False
 
 
 def test___eq__():
@@ -858,81 +855,92 @@ def test_ModularInteger():
     F3 = FF(3)
 
     a = F3(0)
-    assert F3.of_type(a) and a == 0
+    assert isinstance(a, F3.dtype) and a == 0
     a = F3(1)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)
-    assert F3.of_type(a) and a == 2
+    assert isinstance(a, F3.dtype) and a == 2
     a = F3(3)
-    assert F3.of_type(a) and a == 0
+    assert isinstance(a, F3.dtype) and a == 0
     a = F3(4)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
 
     a = F3(F3(0))
-    assert F3.of_type(a) and a == 0
+    assert isinstance(a, F3.dtype) and a == 0
     a = F3(F3(1))
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(F3(2))
-    assert F3.of_type(a) and a == 2
+    assert isinstance(a, F3.dtype) and a == 2
     a = F3(F3(3))
-    assert F3.of_type(a) and a == 0
+    assert isinstance(a, F3.dtype) and a == 0
     a = F3(F3(4))
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
 
     a = -F3(1)
-    assert F3.of_type(a) and a == 2
+    assert isinstance(a, F3.dtype) and a == 2
     a = -F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
 
     a = 2 + F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2) + 2
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2) + F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2) + F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
 
     a = 3 - F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(3) - 2
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(3) - F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(3) - F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
 
     a = 2*F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)*2
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)*F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)*F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
 
     a = 2/F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)/2
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)/F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)/F3(2)
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
+
+    a = 1 % F3(2)
+    assert isinstance(a, F3.dtype) and a == 1
+    a = F3(1) % 2
+    assert isinstance(a, F3.dtype) and a == 1
+    a = F3(1) % F3(2)
+    assert isinstance(a, F3.dtype) and a == 1
+    a = F3(1) % F3(2)
+    assert isinstance(a, F3.dtype) and a == 1
 
     a = F3(2)**0
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
     a = F3(2)**1
-    assert F3.of_type(a) and a == 2
+    assert isinstance(a, F3.dtype) and a == 2
     a = F3(2)**2
-    assert F3.of_type(a) and a == 1
+    assert isinstance(a, F3.dtype) and a == 1
 
     F7 = FF(7)
 
     a = F7(3)**100000000000
-    assert F7.of_type(a) and a == 4
+    assert isinstance(a, F7.dtype) and a == 4
     a = F7(3)**-100000000000
-    assert F7.of_type(a) and a == 2
+    assert isinstance(a, F7.dtype) and a == 2
+    a = F7(3)**S(2)
+    assert isinstance(a, F7.dtype) and a == 2
 
     assert bool(F3(3)) is False
     assert bool(F3(4)) is True
@@ -940,39 +948,60 @@ def test_ModularInteger():
     F5 = FF(5)
 
     a = F5(1)**(-1)
-    assert F5.of_type(a) and a == 1
+    assert isinstance(a, F5.dtype) and a == 1
     a = F5(2)**(-1)
-    assert F5.of_type(a) and a == 3
+    assert isinstance(a, F5.dtype) and a == 3
     a = F5(3)**(-1)
-    assert F5.of_type(a) and a == 2
+    assert isinstance(a, F5.dtype) and a == 2
     a = F5(4)**(-1)
-    assert F5.of_type(a) and a == 4
-
-    if GROUND_TYPES != 'flint':
-        # XXX: This gives a core dump with python-flint...
-        raises(NotInvertible, lambda: F5(0)**(-1))
-        raises(NotInvertible, lambda: F5(5)**(-1))
+    assert isinstance(a, F5.dtype) and a == 4
+
+    assert (F5(1) < F5(2)) is True
+    assert (F5(1) <= F5(2)) is True
+    assert (F5(1) > F5(2)) is False
+    assert (F5(1) >= F5(2)) is False
+
+    assert (F5(3) < F5(2)) is False
+    assert (F5(3) <= F5(2)) is False
+    assert (F5(3) > F5(2)) is True
+    assert (F5(3) >= F5(2)) is True
+
+    assert (F5(1) < F5(7)) is True
+    assert (F5(1) <= F5(7)) is True
+    assert (F5(1) > F5(7)) is False
+    assert (F5(1) >= F5(7)) is False
+
+    assert (F5(3) < F5(7)) is False
+    assert (F5(3) <= F5(7)) is False
+    assert (F5(3) > F5(7)) is True
+    assert (F5(3) >= F5(7)) is True
+
+    assert (F5(1) < 2) is True
+    assert (F5(1) <= 2) is True
+    assert (F5(1) > 2) is False
+    assert (F5(1) >= 2) is False
+
+    assert (F5(3) < 2) is False
+    assert (F5(3) <= 2) is False
+    assert (F5(3) > 2) is True
+    assert (F5(3) >= 2) is True
+
+    assert (F5(1) < 7) is True
+    assert (F5(1) <= 7) is True
+    assert (F5(1) > 7) is False
+    assert (F5(1) >= 7) is False
+
+    assert (F5(3) < 7) is False
+    assert (F5(3) <= 7) is False
+    assert (F5(3) > 7) is True
+    assert (F5(3) >= 7) is True
+
+    raises(NotInvertible, lambda: F5(0)**(-1))
+    raises(NotInvertible, lambda: F5(5)**(-1))
 
     raises(ValueError, lambda: FF(0))
     raises(ValueError, lambda: FF(2.1))
 
-    for n1 in range(5):
-        for n2 in range(5):
-            if GROUND_TYPES != 'flint':
-                with warns_deprecated_sympy():
-                    assert (F5(n1) < F5(n2)) is (n1 < n2)
-                with warns_deprecated_sympy():
-                    assert (F5(n1) <= F5(n2)) is (n1 <= n2)
-                with warns_deprecated_sympy():
-                    assert (F5(n1) > F5(n2)) is (n1 > n2)
-                with warns_deprecated_sympy():
-                    assert (F5(n1) >= F5(n2)) is (n1 >= n2)
-            else:
-                raises(TypeError, lambda: F5(n1) < F5(n2))
-                raises(TypeError, lambda: F5(n1) <= F5(n2))
-                raises(TypeError, lambda: F5(n1) > F5(n2))
-                raises(TypeError, lambda: F5(n1) >= F5(n2))
-
 def test_QQ_int():
     assert int(QQ(2**2000, 3**1250)) == 455431
     assert int(QQ(2**100, 3)) == 422550200076076467165567735125
diff --git a/sympy/sets/tests/test_fancysets.py b/sympy/sets/tests/test_fancysets.py
index 3582d5aabb..08f097b934 100644
--- a/sympy/sets/tests/test_fancysets.py
+++ b/sympy/sets/tests/test_fancysets.py
@@ -168,12 +168,17 @@ def test_halfcircle():
     assert (0, -1) not in halfcircle
     assert (0, 0) in halfcircle
     assert halfcircle._contains((r, 0)) is None
-    # This one doesn't work:
-    #assert (r, 2*pi) not in halfcircle
-
     assert not halfcircle.is_iterable
 
 
+@XFAIL
+def test_halfcircle_fail():
+    r, th = symbols('r, theta', real=True)
+    L = Lambda(((r, th),), (r*cos(th), r*sin(th)))
+    halfcircle = ImageSet(L, Interval(0, 1)*Interval(0, pi))
+    assert (r, 2*pi) not in halfcircle
+
+
 def test_ImageSet_iterator_not_injective():
     L = Lambda(x, x - x % 2)  # produces 0, 2, 2, 4, 4, 6, 6, ...
     evens = ImageSet(L, S.Naturals)
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index ed9a337be9..83778eb197 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -129,16 +129,24 @@ def ireal(x, s=S.Reals):
     assert invert_real(Abs(x**31 + x + 1), y, x) == (lhs, base_values)
 
     assert dumeq(invert_real(sin(x), y, x),
-        (x, imageset(Lambda(n, n*pi + (-1)**n*asin(y)), S.Integers)))
+        (x, Union(
+            ImageSet(Lambda(n, 2*n*pi + asin(y)), S.Integers),
+            ImageSet(Lambda(n, pi*2*n + pi - asin(y)), S.Integers))))
 
     assert dumeq(invert_real(sin(exp(x)), y, x),
-        (x, imageset(Lambda(n, log((-1)**n*asin(y) + n*pi)), S.Integers)))
+        (x, Union(
+            ImageSet(Lambda(n, log(2*n*pi + asin(y))), S.Integers),
+            ImageSet(Lambda(n, log(pi*2*n + pi - asin(y))), S.Integers))))
 
     assert dumeq(invert_real(csc(x), y, x),
-        (x, imageset(Lambda(n, n*pi + (-1)**n*acsc(y)), S.Integers)))
+        (x, Union(
+            ImageSet(Lambda(n, 2*n*pi + acsc(y)), S.Integers),
+            ImageSet(Lambda(n, pi*2*n + pi - acsc(y)), S.Integers))))
 
     assert dumeq(invert_real(csc(exp(x)), y, x),
-        (x, imageset(Lambda(n, log((-1)**n*acsc(y) + n*pi)), S.Integers)))
+        (x, Union(
+            ImageSet(Lambda(n, log(2*n*pi + acsc(y))), S.Integers),
+            ImageSet(Lambda(n, log(pi*2*n + pi - acsc(y))), S.Integers))))
 
     assert dumeq(invert_real(cos(x), y, x),
         (x, Union(imageset(Lambda(n, 2*n*pi + acos(y)), S.Integers), \
@@ -349,6 +357,19 @@ def test_solve_invert():
     assert solveset_real(2**x - 10, x) == FiniteSet(1 + log(5)/log(2))
 
 
+def test_issue_25768():
+    assert dumeq(solveset_real(sin(x) - S.Half, x), Union(
+        ImageSet(Lambda(n, pi*2*n + pi/6), S.Integers),
+        ImageSet(Lambda(n, pi*2*n + pi*5/6), S.Integers)))
+    n1 = solveset_real(sin(x) - 0.5, x).n(5)
+    n2 = solveset_real(sin(x) - S.Half, x).n(5)
+    # help pass despite fp differences
+    eq = [i.replace(
+        lambda x:x.is_Float,
+        lambda x:Rational(x).limit_denominator(1000)) for i in (n1, n2)]
+    assert dumeq(*eq),(n1,n2)
+
+
 def test_errorinverses():
     assert solveset_real(erf(x) - S.Half, x) == \
         FiniteSet(erfinv(S.Half))
@@ -657,10 +678,11 @@ def test_solve_abs():
     for si in sol.subs(reps):
         assert not eqab.subs(x, si)
     assert dumeq(solveset(Eq(sin(Abs(x)), 1), x, domain=S.Reals), Union(
-        Intersection(Interval(0, oo),
-            ImageSet(Lambda(n, (-1)**n*pi/2 + n*pi), S.Integers)),
-        Intersection(Interval(-oo, 0),
-            ImageSet(Lambda(n, n*pi - (-1)**(-n)*pi/2), S.Integers))))
+        Intersection(Interval(0, oo), Union(
+        Intersection(ImageSet(Lambda(n, 2*n*pi + 3*pi/2), S.Integers),
+            Interval(-oo, 0)),
+        Intersection(ImageSet(Lambda(n, 2*n*pi + pi/2), S.Integers),
+            Interval(0, oo))))))
 
 
 def test_issue_9824():
@@ -849,9 +871,7 @@ def test_solve_trig():
 
     assert dumeq(solveset(sin(y + a) - sin(y), a, domain=S.Reals),
         Union(ImageSet(Lambda(n, 2*n*pi), S.Integers),
-        Intersection(ImageSet(Lambda(n, -I*(I*(
-        2*n*pi + arg(-exp(-2*I*y))) +
-        2*im(y))), S.Integers), S.Reals)))
+        Intersection(ImageSet(Lambda(n, -I*(I*(2*n*pi + arg(-exp(-2*I*y))) + 2*im(y))), S.Integers), S.Reals)))
 
     assert dumeq(solveset_real(sin(2*x)*cos(x) + cos(2*x)*sin(x)-1, x),
         ImageSet(Lambda(n, n*pi*Rational(2, 3) + pi/6), S.Integers))
@@ -1272,12 +1292,12 @@ def test__solveset_multi():
             [Interval(0, pi), Interval(-1, 1)]) == FiniteSet((0, 1), (pi, -1))
     assert _solveset_multi([r*cos(theta)-1, r*sin(theta)], [r, theta],
             [Interval(0, 1), Interval(0, pi)]) == FiniteSet((1, 0))
-    #assert _solveset_multi([r*cos(theta)-r, r*sin(theta)], [r, theta],
-    #        [Interval(0, 1), Interval(0, pi)]) == ?
-    assert dumeq(_solveset_multi([r*cos(theta)-r, r*sin(theta)], [r, theta],
-            [Interval(0, 1), Interval(0, pi)]), Union(
-            ImageSet(Lambda(((r,),), (r, 0)), ImageSet(Lambda(r, (r,)), Interval(0, 1))),
-            ImageSet(Lambda(((theta,),), (0, theta)), ImageSet(Lambda(theta, (theta,)), Interval(0, pi)))))
+    assert _solveset_multi([r*cos(theta)-r, r*sin(theta)], [r, theta],
+           [Interval(0, 1), Interval(0, pi)]) == Union(
+           ImageSet(Lambda(((r,),), (r, 0)),
+           ImageSet(Lambda(r, (r,)), Interval(0, 1))),
+           ImageSet(Lambda(((theta,),), (0, theta)),
+           ImageSet(Lambda(theta, (theta,)), Interval(0, pi))))
 
 
 def test_conditionset():
@@ -1408,11 +1428,6 @@ def test_linear_eq_to_matrix():
     assert linear_eq_to_matrix(Eq(x + 2, 1), x) == (
         Matrix([[1]]), Matrix([[-1]]))
 
-    # issue 25423
-    raises(TypeError, lambda: linear_eq_to_matrix([], {x, y}))
-    raises(TypeError, lambda: linear_eq_to_matrix([x + y], {x, y}))
-    raises(ValueError, lambda: linear_eq_to_matrix({x + y}, (x, y)))
-
 
 def test_issue_16577():
     assert linear_eq_to_matrix(Eq(a*(2*x + 3*y) + 4*y, 5), x, y) == (
@@ -1589,6 +1604,7 @@ def test_solve_decomposition():
     assert solve_decomposition(f7, x, S.Reals) == S.EmptySet
     assert solve_decomposition(x, x, Interval(1, 2)) == S.EmptySet
 
+
 # nonlinsolve testcases
 def test_nonlinsolve_basic():
     assert nonlinsolve([],[]) == S.EmptySet
@@ -2249,7 +2265,6 @@ def test_issue_17933():
     eq2 = x*cos(45) - y*sin(q)
     eq3 = 9*x*sin(45)/10 + y*cos(q)
     eq4 = 9*x*cos(45)/10 + y*sin(z) - z
-
     assert nonlinsolve([eq1, eq2, eq3, eq4], x, y, z, q) ==\
         FiniteSet((0, 0, 0, q))