Repository: sympy/sympy

Issue:
Incorrect Domain of Inverse Trigonometric Functions
Domain of the arc sin function is [-1, 1]. 
Although, the following code returns unhelpful answers. 
Is there a workaround?
Does this happen for other kinds of functions as well?

Code:
```
import sympy as sp
from sympy.calculus.util import continuous_domain, function_range

x = sp.symbols("x")
func = sp.asin(x)
dom = continuous_domain(func, x, sp.Reals)
print(dom)
```

Output
```
Reals
```





Discussion:

Patch:
diff --git a/sympy/calculus/singularities.py b/sympy/calculus/singularities.py
index 3a2ea1a17c..5adafc59ef 100644
--- a/sympy/calculus/singularities.py
+++ b/sympy/calculus/singularities.py
@@ -21,6 +21,8 @@
 from sympy.core.sympify import sympify
 from sympy.functions.elementary.exponential import log
 from sympy.functions.elementary.trigonometric import sec, csc, cot, tan, cos
+from sympy.functions.elementary.hyperbolic import (
+    sech, csch, coth, tanh, cosh, asech, acsch, atanh, acoth)
 from sympy.utilities.misc import filldedent
 
 
@@ -90,13 +92,19 @@ def singularities(expression, symbol, domain=None):
         domain = S.Reals if symbol.is_real else S.Complexes
     try:
         sings = S.EmptySet
-        for i in expression.rewrite([sec, csc, cot, tan], cos).atoms(Pow):
+        e = expression.rewrite([sec, csc, cot, tan], cos)
+        e = e.rewrite([sech, csch, coth, tanh], cosh)
+        for i in e.atoms(Pow):
             if i.exp.is_infinite:
                 raise NotImplementedError
             if i.exp.is_negative:
+                # XXX: exponent of varying sign not handled
                 sings += solveset(i.base, symbol, domain)
-        for i in expression.atoms(log):
+        for i in expression.atoms(log, asech, acsch):
             sings += solveset(i.args[0], symbol, domain)
+        for i in expression.atoms(atanh, acoth):
+            sings += solveset(i.args[0] - 1, symbol, domain)
+            sings += solveset(i.args[0] + 1, symbol, domain)
         return sings
     except NotImplementedError:
         raise NotImplementedError(filldedent('''
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index e140396d8c..e12eb1cc19 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -1,7 +1,7 @@
 from .accumulationbounds import AccumBounds, AccumulationBounds # noqa: F401
 from .singularities import singularities
 from sympy.core import Pow, S
-from sympy.core.function import diff, expand_mul
+from sympy.core.function import diff, expand_mul, Function
 from sympy.core.kind import NumberKind
 from sympy.core.mod import Mod
 from sympy.core.numbers import equal_valued
@@ -10,23 +10,30 @@
 from sympy.core.sympify import _sympify
 from sympy.functions.elementary.complexes import Abs, im, re
 from sympy.functions.elementary.exponential import exp, log
+from sympy.functions.elementary.integers import frac
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.functions.elementary.trigonometric import (
-    TrigonometricFunction, sin, cos, csc, sec)
+    TrigonometricFunction, sin, cos, tan, cot, csc, sec,
+    asin, acos, acot, atan, asec, acsc)
+from sympy.functions.elementary.hyperbolic import (sinh, cosh, tanh, coth,
+    sech, csch, asinh, acosh, atanh, acoth, asech, acsch)
 from sympy.polys.polytools import degree, lcm_list
 from sympy.sets.sets import (Interval, Intersection, FiniteSet, Union,
                              Complement)
 from sympy.sets.fancysets import ImageSet
+from sympy.sets.conditionset import ConditionSet
 from sympy.utilities import filldedent
 from sympy.utilities.iterables import iterable
 
 
 def continuous_domain(f, symbol, domain):
     """
-    Returns the intervals in the given domain for which the function
-    is continuous.
-    This method is limited by the ability to determine the various
+    Returns the domain on which the function expression f is continuous.
+
+    This function is limited by the ability to determine the various
     singularities and discontinuities of the given function.
+    The result is either given as a union of intervals or constructed using
+    other set operations.
 
     Parameters
     ==========
@@ -69,24 +76,85 @@ def continuous_domain(f, symbol, domain):
     """
     from sympy.solvers.inequalities import solve_univariate_inequality
 
-    if domain.is_subset(S.Reals):
-        constrained_interval = domain
-        for atom in f.atoms(Pow):
-            den = atom.exp.as_numer_denom()[1]
-            if den.is_even and den.is_nonzero:
-                constraint = solve_univariate_inequality(atom.base >= 0,
-                                                         symbol).as_set()
-                constrained_interval = Intersection(constraint,
-                                                    constrained_interval)
-
-        for atom in f.atoms(log):
-            constraint = solve_univariate_inequality(atom.args[0] > 0,
-                                                     symbol).as_set()
-            constrained_interval = Intersection(constraint,
-                                                constrained_interval)
-
+    if not domain.is_subset(S.Reals):
+        raise NotImplementedError(filldedent('''
+            Domain must be a subset of S.Reals.
+            '''))
+    implemented = [Pow, exp, log, Abs, frac,
+                   sin, cos, tan, cot, sec, csc,
+                   asin, acos, atan, acot, asec, acsc,
+                   sinh, cosh, tanh, coth, sech, csch,
+                   asinh, acosh, atanh, acoth, asech, acsch]
+    used = [fct.func for fct in f.atoms(Function) if fct.has(symbol)]
+    if any(func not in implemented for func in used):
+        raise NotImplementedError(filldedent('''
+            Unable to determine the domain of the given function.
+            '''))
+
+    x = Symbol('x')
+    constraints = {
+        log: (x > 0,),
+        asin: (x >= -1, x <= 1),
+        acos: (x >= -1, x <= 1),
+        acosh: (x >= 1,),
+        atanh: (x > -1, x < 1),
+        asech: (x > 0, x <= 1)
+    }
+    constraints_union = {
+        asec: (x <= -1, x >= 1),
+        acsc: (x <= -1, x >= 1),
+        acoth: (x < -1, x > 1)
+    }
+
+    cont_domain = domain
+    for atom in f.atoms(Pow):
+        den = atom.exp.as_numer_denom()[1]
+        if atom.exp.is_rational and den.is_odd:
+            pass    # 0**negative handled by singularities()
+        else:
+            constraint = solve_univariate_inequality(atom.base >= 0,
+                                                        symbol).as_set()
+            cont_domain = Intersection(constraint, cont_domain)
+
+    for atom in f.atoms(Function):
+        if atom.func in constraints:
+            for c in constraints[atom.func]:
+                constraint_relational = c.subs(x, atom.args[0])
+                constraint_set = solve_univariate_inequality(
+                    constraint_relational, symbol).as_set()
+                cont_domain = Intersection(constraint_set, cont_domain)
+        elif atom.func in constraints_union:
+            constraint_set = S.EmptySet
+            for c in constraints_union[atom.func]:
+                constraint_relational = c.subs(x, atom.args[0])
+                constraint_set += solve_univariate_inequality(
+                    constraint_relational, symbol).as_set()
+            cont_domain = Intersection(constraint_set, cont_domain)
+        # XXX: the discontinuities below could be factored out in
+        # a new "discontinuities()".
+        elif atom.func == acot:
+            from sympy.solvers.solveset import solveset_real
+            # Sympy's acot() has a step discontinuity at 0. Since it's
+            # neither an essential singularity nor a pole, singularities()
+            # will not report it. But it's still relevant for determining
+            # the continuity of the function f.
+            cont_domain -= solveset_real(atom.args[0], symbol)
+            # Note that the above may introduce spurious discontinuities, e.g.
+            # for abs(acot(x)) at 0.
+        elif atom.func == frac:
+            from sympy.solvers.solveset import solveset_real
+            r = function_range(atom.args[0], symbol, domain)
+            r = Intersection(r, S.Integers)
+            if r.is_finite_set:
+                discont = S.EmptySet
+                for n in r:
+                    discont += solveset_real(atom.args[0]-n, symbol)
+            else:
+                discont = ConditionSet(
+                    symbol, S.Integers.contains(atom.args[0]), cont_domain)
+            cont_domain -= discont
 
-    return constrained_interval - singularities(f, symbol, domain)
+    return cont_domain - singularities(f, symbol, domain)
 
 
 def function_range(f, symbol, domain):
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index a25e18435f..bda9a8d7bf 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -11,7 +11,6 @@
 from sympy.core.sorting import default_sort_key, ordered
 from sympy.core.symbol import Symbol, symbols
 from sympy.core.sympify import _sympify
-from sympy.external.gmpy import jacobi, remove
 from sympy.functions.elementary.complexes import sign
 from sympy.functions.elementary.integers import floor
 from sympy.functions.elementary.miscellaneous import sqrt
@@ -302,7 +301,9 @@ def solve(self, parameters=None, limit=None):
             q, r = divmod(c, coeff[var[0]])
             if not r:
                 result.add((q,))
-            return result
+                return result
+            else:
+                return result
 
         '''
         base_solution_linear() can solve diophantine equations of the form:
@@ -847,11 +848,12 @@ def unpack_sol(sol):
                 result.add(_remove_gcd(s[0], -coeff[x*z], s[1]))
                 return result
 
-            var[0], var[1] = _var[1], _var[0]
-            y_0, x_0, z_0 = unpack_sol(_diop_ternary_quadratic(var, coeff))
-            if x_0 is not None:
-                result.add((x_0, y_0, z_0))
-            return result
+            else:
+                var[0], var[1] = _var[1], _var[0]
+                y_0, x_0, z_0 = unpack_sol(_diop_ternary_quadratic(var, coeff))
+                if x_0 is not None:
+                    result.add((x_0, y_0, z_0))
+                return result
 
         if coeff[x**2] == 0:
             # If the coefficient of x is zero change the variables
@@ -938,8 +940,11 @@ def matches(self):
             return False
         if not self.homogeneous_order:
             return True
-        # there may be Pow keys like x**2 or Mul keys like x*y
-        return any(k.is_Mul for k in self.coeff) and not self.homogeneous
+        else:
+            # there may be Pow keys like x**2 or Mul keys like x*y
+            if any(k.is_Mul for k in self.coeff): # cross terms
+                return not self.homogeneous
+        return False
 
 
 class HomogeneousGeneralQuadratic(DiophantineEquationType):
@@ -958,8 +963,11 @@ def matches(self):
             return False
         if not self.homogeneous_order:
             return False
-        # there may be Pow keys like x**2 or Mul keys like x*y
-        return any(k.is_Mul for k in self.coeff) and self.homogeneous
+        else:
+            # there may be Pow keys like x**2 or Mul keys like x*y
+            if any(k.is_Mul for k in self.coeff): # cross terms
+                return self.homogeneous
+        return False
 
 
 class GeneralSumOfSquares(DiophantineEquationType):
@@ -1354,7 +1362,9 @@ def diophantine(eq, param=symbols("t", integer=True), syms=None,
             dsol = diophantine(d)
             good = diophantine(n) - dsol
             return {s for s in good if _mexpand(d.subs(zip(var, s)))}
-        eq = factor_terms(n)
+        else:
+            eq = n
+        eq = factor_terms(eq)
         assert not eq.is_number
         eq = eq.as_independent(*var, as_Add=False)[1]
         p = Poly(eq)
@@ -1804,22 +1814,27 @@ def base_solution_linear(c, a, b, t=None):
     a, b, c = _remove_gcd(a, b, c)
 
     if c == 0:
-        if t is None:
+        if t is not None:
+            if b < 0:
+                t = -t
+            return (b*t, -a*t)
+        else:
             return (0, 0)
-        if b < 0:
-            t = -t
-        return (b*t, -a*t)
-
-    x0, y0, d = igcdex(abs(a), abs(b))
-    x0 *= sign(a)
-    y0 *= sign(b)
-    if c % d:
-        return (None, None)
-    if t is None:
-        return (c*x0, c*y0)
-    if b < 0:
-        t = -t
-    return (c*x0 + b*t, c*y0 - a*t)
+    else:
+        x0, y0, d = igcdex(abs(a), abs(b))
+
+        x0 *= sign(a)
+        y0 *= sign(b)
+
+        if divisible(c, d):
+            if t is not None:
+                if b < 0:
+                    t = -t
+                return (c*x0 + b*t, c*y0 - a*t)
+            else:
+                return (c*x0, c*y0)
+        else:
+            return (None, None)
 
 
 def diop_univariate(eq):
@@ -1990,20 +2005,20 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
     if D < 0:
         if N == 0:
             return [(0, 0)]
-        if N < 0:
+        elif N < 0:
             return []
-        # N > 0:
-        sol = []
-        for d in divisors(square_factor(N)):
-            sols = cornacchia(1, -D, N // d**2)
-            if sols:
-                for x, y in sols:
-                    sol.append((d*x, d*y))
-                    if D == -1:
-                        sol.append((d*y, d*x))
-        return sol
-
-    if D == 0:
+        elif N > 0:
+            sol = []
+            for d in divisors(square_factor(N)):
+                sols = cornacchia(1, -D, N // d**2)
+                if sols:
+                    for x, y in sols:
+                        sol.append((d*x, d*y))
+                        if D == -1:
+                            sol.append((d*y, d*x))
+            return sol
+
+    elif D == 0:
         if N < 0:
             return []
         if N == 0:
@@ -2011,118 +2026,125 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
         sN, _exact = integer_nthroot(N, 2)
         if _exact:
             return [(sN, t)]
-        return []
+        else:
+            return []
 
-    # D > 0
-    sD, _exact = integer_nthroot(D, 2)
-    if _exact:
-        if N == 0:
-            return [(sD*t, t)]
-
-        sol = []
-        for y in range(floor(sign(N)*(N - 1)/(2*sD)) + 1):
-            try:
-                sq, _exact = integer_nthroot(D*y**2 + N, 2)
-            except ValueError:
-                _exact = False
-            if _exact:
-                sol.append((sq, y))
-        return sol
+    else:  # D > 0
+        sD, _exact = integer_nthroot(D, 2)
+        if _exact:
+            if N == 0:
+                return [(sD*t, t)]
+            else:
+                sol = []
 
-    if 1 < N**2 < D:
-        # It is much faster to call `_special_diop_DN`.
-        return _special_diop_DN(D, N)
+                for y in range(floor(sign(N)*(N - 1)/(2*sD)) + 1):
+                    try:
+                        sq, _exact = integer_nthroot(D*y**2 + N, 2)
+                    except ValueError:
+                        _exact = False
+                    if _exact:
+                        sol.append((sq, y))
 
-    if N == 0:
-        return [(0, 0)]
+                return sol
 
-    if abs(N) == 1:
-        pqa = PQa(0, 1, D)
-        j = 0
-        G = []
-        B = []
+        elif 1 < N**2 < D:
+            # It is much faster to call `_special_diop_DN`.
+            return _special_diop_DN(D, N)
 
-        for i in pqa:
+        else:
+            if N == 0:
+                return [(0, 0)]
 
-            a = i[2]
-            G.append(i[5])
-            B.append(i[4])
+            elif abs(N) == 1:
 
-            if j != 0 and a == 2*sD:
-                break
-            j = j + 1
+                pqa = PQa(0, 1, D)
+                j = 0
+                G = []
+                B = []
 
-        if _odd(j):
+                for i in pqa:
 
-            if N == -1:
-                x = G[j - 1]
-                y = B[j - 1]
-            else:
-                count = j
-                while count < 2*j - 1:
-                    i = next(pqa)
+                    a = i[2]
                     G.append(i[5])
                     B.append(i[4])
-                    count += 1
 
-                x = G[count]
-                y = B[count]
-        else:
-            if N == 1:
-                x = G[j - 1]
-                y = B[j - 1]
-            else:
-                return []
+                    if j != 0 and a == 2*sD:
+                        break
+                    j = j + 1
 
-        return [(x, y)]
+                if _odd(j):
 
-    fs = []
-    sol = []
-    div = divisors(N)
+                    if N == -1:
+                        x = G[j - 1]
+                        y = B[j - 1]
+                    else:
+                        count = j
+                        while count < 2*j - 1:
+                            i = next(pqa)
+                            G.append(i[5])
+                            B.append(i[4])
+                            count += 1
+
+                        x = G[count]
+                        y = B[count]
+                else:
+                    if N == 1:
+                        x = G[j - 1]
+                        y = B[j - 1]
+                    else:
+                        return []
 
-    for d in div:
-        if divisible(N, d**2):
-            fs.append(d)
+                return [(x, y)]
 
-    for f in fs:
-        m = N // f**2
+            else:
 
-        zs = sqrt_mod(D, abs(m), all_roots=True)
-        zs = [i for i in zs if i <= abs(m) // 2 ]
+                fs = []
+                sol = []
+                div = divisors(N)
 
-        if abs(m) != 2:
-            zs = zs + [-i for i in zs if i]  # omit dupl 0
+                for d in div:
+                    if divisible(N, d**2):
+                        fs.append(d)
 
-        for z in zs:
+                for f in fs:
+                    m = N // f**2
 
-            pqa = PQa(z, abs(m), D)
-            j = 0
-            G = []
-            B = []
+                    zs = sqrt_mod(D, abs(m), all_roots=True)
+                    zs = [i for i in zs if i <= abs(m) // 2 ]
 
-            for i in pqa:
+                    if abs(m) != 2:
+                        zs = zs + [-i for i in zs if i]  # omit dupl 0
 
-                G.append(i[5])
-                B.append(i[4])
+                    for z in zs:
 
-                if j != 0 and abs(i[1]) == 1:
-                    r = G[j-1]
-                    s = B[j-1]
+                        pqa = PQa(z, abs(m), D)
+                        j = 0
+                        G = []
+                        B = []
 
-                    if r**2 - D*s**2 == m:
-                        sol.append((f*r, f*s))
+                        for i in pqa:
 
-                    elif diop_DN(D, -1) != []:
-                        a = diop_DN(D, -1)
-                        sol.append((f*(r*a[0][0] + a[0][1]*s*D), f*(r*a[0][1] + s*a[0][0])))
+                            G.append(i[5])
+                            B.append(i[4])
 
-                    break
+                            if j != 0 and abs(i[1]) == 1:
+                                r = G[j-1]
+                                s = B[j-1]
 
-                j = j + 1
-                if j == length(z, abs(m), D):
-                    break
+                                if r**2 - D*s**2 == m:
+                                    sol.append((f*r, f*s))
 
-    return sol
+                                elif diop_DN(D, -1) != []:
+                                    a = diop_DN(D, -1)
+                                    sol.append((f*(r*a[0][0] + a[0][1]*s*D), f*(r*a[0][1] + s*a[0][0])))
+
+                                break
+
+                            j = j + 1
+                            if j == length(z, abs(m), D):
+                                break
+
+                return sol
 
 
 def _special_diop_DN(D, N):
@@ -2401,28 +2423,32 @@ def diop_bf_DN(D, N, t=symbols("t", integer=True)):
     a = diop_DN(D, 1)
     u = a[0][0]
 
-    if N == 0:
-        if D < 0:
-            return [(0, 0)]
-        if D == 0:
-            return [(0, t)]
-        sD, _exact = integer_nthroot(D, 2)
-        if _exact:
-            return [(sD*t, t), (-sD*t, t)]
-        return [(0, 0)]
-
     if abs(N) == 1:
         return diop_DN(D, N)
 
-    if N > 1:
+    elif N > 1:
         L1 = 0
         L2 = integer_nthroot(int(N*(u - 1)/(2*D)), 2)[0] + 1
-    else: # N < -1
+
+    elif N < -1:
         L1, _exact = integer_nthroot(-int(N/D), 2)
         if not _exact:
             L1 += 1
         L2 = integer_nthroot(-int(N*(u + 1)/(2*D)), 2)[0] + 1
 
+    else:  # N = 0
+        if D < 0:
+            return [(0, 0)]
+        elif D == 0:
+            return [(0, t)]
+        else:
+            sD, _exact = integer_nthroot(D, 2)
+            if _exact:
+                return [(sD*t, t), (-sD*t, t)]
+            else:
+                return [(0, 0)]
+
+
     for y in range(L1, L2):
         try:
             x, _exact = integer_nthroot(N + D*y**2, 2)
@@ -2627,27 +2653,31 @@ def _transformation_to_DN(var, coeff):
         A_0, B_0 = _transformation_to_DN([X, Y], coeff)
         return Matrix(2, 2, [S.One/B, -S(C)/B, 0, 1])*A_0, Matrix(2, 2, [S.One/B, -S(C)/B, 0, 1])*B_0
 
-    if d:
-        B, C = _rational_pq(2*a, d)
-        A, T = _rational_pq(a, B**2)
+    else:
+        if d:
+            B, C = _rational_pq(2*a, d)
+            A, T = _rational_pq(a, B**2)
 
-        # eq_2 = A*X**2 + c*T*Y**2 + e*T*Y + f*T - A*C**2
-        coeff = {X**2: A, X*Y: 0, Y**2: c*T, X: 0, Y: e*T, 1: f*T - A*C**2}
-        A_0, B_0 = _transformation_to_DN([X, Y], coeff)
-        return Matrix(2, 2, [S.One/B, 0, 0, 1])*A_0, Matrix(2, 2, [S.One/B, 0, 0, 1])*B_0 + Matrix([-S(C)/B, 0])
+            # eq_2 = A*X**2 + c*T*Y**2 + e*T*Y + f*T - A*C**2
+            coeff = {X**2: A, X*Y: 0, Y**2: c*T, X: 0, Y: e*T, 1: f*T - A*C**2}
+            A_0, B_0 = _transformation_to_DN([X, Y], coeff)
+            return Matrix(2, 2, [S.One/B, 0, 0, 1])*A_0, Matrix(2, 2, [S.One/B, 0, 0, 1])*B_0 + Matrix([-S(C)/B, 0])
 
-    if e:
-        B, C = _rational_pq(2*c, e)
-        A, T = _rational_pq(c, B**2)
+        else:
+            if e:
+                B, C = _rational_pq(2*c, e)
+                A, T = _rational_pq(c, B**2)
 
-        # eq_3 = a*T*X**2 + A*Y**2 + f*T - A*C**2
-        coeff = {X**2: a*T, X*Y: 0, Y**2: A, X: 0, Y: 0, 1: f*T - A*C**2}
-        A_0, B_0 = _transformation_to_DN([X, Y], coeff)
-        return Matrix(2, 2, [1, 0, 0, S.One/B])*A_0, Matrix(2, 2, [1, 0, 0, S.One/B])*B_0 + Matrix([0, -S(C)/B])
+                # eq_3 = a*T*X**2 + A*Y**2 + f*T - A*C**2
+                coeff = {X**2: a*T, X*Y: 0, Y**2: A, X: 0, Y: 0, 1: f*T - A*C**2}
+                A_0, B_0 = _transformation_to_DN([X, Y], coeff)
+                return Matrix(2, 2, [1, 0, 0, S.One/B])*A_0, Matrix(2, 2, [1, 0, 0, S.One/B])*B_0 + Matrix([0, -S(C)/B])
 
-    # TODO: pre-simplification: Not necessary but may simplify
-    # the equation.
-    return Matrix(2, 2, [S.One/a, 0, 0, 1]), Matrix([0, 0])
+            else:
+                # TODO: pre-simplification: Not necessary but may simplify
+                # the equation.
+
+                return Matrix(2, 2, [S.One/a, 0, 0, 1]), Matrix([0, 0])
 
 
 def find_DN(eq):
@@ -2843,11 +2873,12 @@ def _transformation_to_normal(var, coeff):
             T.col_swap(0, 2)
             return T
 
-        _var[0], _var[1] = var[1], var[0]
-        T = _transformation_to_normal(_var, coeff)
-        T.row_swap(0, 1)
-        T.col_swap(0, 1)
-        return T
+        else:
+            _var[0], _var[1] = var[1], var[0]
+            T = _transformation_to_normal(_var, coeff)
+            T.row_swap(0, 1)
+            T.col_swap(0, 1)
+            return T
 
     # Apply the transformation x --> X - (B*Y + C*Z)/(2*A)
     if coeff[x*y] != 0 or coeff[x*z] != 0:
@@ -2877,21 +2908,24 @@ def _transformation_to_normal(var, coeff):
                 # Apply transformation y -> Y + Z ans z -> Y - Z
                 return Matrix(3, 3, [1, 0, 0, 0, 1, 1, 0, 1, -1])
 
-            # Ax**2 + E*y*z + F*z**2  = 0
-            _var[0], _var[2] = var[2], var[0]
+            else:
+                # Ax**2 + E*y*z + F*z**2  = 0
+                _var[0], _var[2] = var[2], var[0]
+                T = _transformation_to_normal(_var, coeff)
+                T.row_swap(0, 2)
+                T.col_swap(0, 2)
+                return T
+
+        else:
+            # A*x**2 + D*y**2 + E*y*z + F*z**2 = 0, F may be zero
+            _var[0], _var[1] = var[1], var[0]
             T = _transformation_to_normal(_var, coeff)
-            T.row_swap(0, 2)
-            T.col_swap(0, 2)
+            T.row_swap(0, 1)
+            T.col_swap(0, 1)
             return T
 
-        # A*x**2 + D*y**2 + E*y*z + F*z**2 = 0, F may be zero
-        _var[0], _var[1] = var[1], var[0]
-        T = _transformation_to_normal(_var, coeff)
-        T.row_swap(0, 1)
-        T.col_swap(0, 1)
-        return T
-
-    return Matrix.eye(3)
+    else:
+        return Matrix.eye(3)
 
 
 def parametrize_ternary_quadratic(eq):
@@ -3568,24 +3602,6 @@ def prime_as_sum_of_two_squares(p):
     Represent a prime `p` as a unique sum of two squares; this can
     only be done if the prime is congruent to 1 mod 4.
 
-    Parameters
-    ==========
-
-    p : Integer
-        A prime that is congruent to 1 mod 4
-
-    Returns
-    =======
-
-    (int, int) | None : Pair of positive integers ``(x, y)`` satisfying ``x**2 + y**2 = p``.
-                        None if ``p`` is not congruent to 1 mod 4.
-
-    Raises
-    ======
-
-    ValueError
-        If ``p`` is not prime number
-
     Examples
     ========
 
@@ -3598,38 +3614,29 @@ def prime_as_sum_of_two_squares(p):
     =========
 
     .. [1] Representing a number as a sum of four squares, [online],
-           Available: https://schorn.ch/lagrange.html
+        Available: https://schorn.ch/lagrange.html
 
     See Also
     ========
-
-    sum_of_squares
-
+    sum_of_squares()
     """
-    p = as_int(p)
-    if p % 4 != 1:
+    if not p % 4 == 1:
         return
-    if not isprime(p):
-        raise ValueError("p should be a prime number")
 
     if p % 8 == 5:
-        # Legendre symbol (2/p) == -1 if p % 8 in [3, 5]
         b = 2
-    elif p % 12 == 5:
-        # Legendre symbol (3/p) == -1 if p % 12 in [5, 7]
-        b = 3
-    elif p % 5 in [2, 3]:
-        # Legendre symbol (5/p) == -1 if p % 5 in [2, 3]
-        b = 5
     else:
-        b = 7
-        while jacobi(b, p) == 1:
+        b = 3
+
+        while pow(b, (p - 1) // 2, p) == 1:
             b = nextprime(b)
 
-    b = pow(b, p >> 2, p)
+    b = pow(b, (p - 1) // 4, p)
     a = p
+
     while b**2 > p:
         a, b = b, a % b
+
     return (int(a % b), int(b))  # convert from long
 
 
@@ -3841,7 +3848,8 @@ def power_representation(n, p, k, zeros=False):
         return
 
     if p == 1:
-        yield from partition(n, k, zeros=zeros)
+        for t in partition(n, k, zeros=zeros):
+            yield t
         return
 
     if p == 2:
@@ -3978,12 +3986,17 @@ def _can_do_sum_of_squares(n, k):
             if n % 4 == 1:
                 return 1  # signal that it was prime
             return False
-        # n is a composite number
-        # we can proceed iff no prime factor in the form 4*k + 3
-        # has an odd multiplicity
-        return all(p % 4 !=3 or m % 2 == 0 for p, m in factorint(n).items())
+        else:
+            f = factorint(n)
+            for p, m in f.items():
+                # we can proceed iff no prime factor in the form 4*k + 3
+                # has an odd multiplicity
+                if (p % 4 == 3) and m % 2:
+                    return False
+            return True
     if k == 3:
-        return remove(n, 4)[0] % 8 != 7
+        if (n//4**multiplicity(4, n)) % 8 == 7:
+            return False
     # every number can be written as a sum of 4 squares; for k > 4 partitions
     # can be 0
     return True


Test patch:
diff --git a/sympy/calculus/tests/test_singularities.py b/sympy/calculus/tests/test_singularities.py
index 26e1cb875b..74428dc269 100644
--- a/sympy/calculus/tests/test_singularities.py
+++ b/sympy/calculus/tests/test_singularities.py
@@ -1,7 +1,11 @@
-from sympy.core.numbers import (I, Rational, oo)
+from sympy.core.numbers import (I, Rational, pi, oo)
 from sympy.core.singleton import S
-from sympy.core.symbol import Symbol
+from sympy.core.symbol import Symbol, Dummy
+from sympy.core.function import Lambda
 from sympy.functions.elementary.exponential import (exp, log)
+from sympy.functions.elementary.trigonometric import sec, csc
+from sympy.functions.elementary.hyperbolic import (coth, sech,
+                                                   atanh, asech, acoth, acsch)
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.calculus.singularities import (
     singularities,
@@ -11,7 +15,7 @@
     is_strictly_decreasing,
     is_monotonic
 )
-from sympy.sets import Interval, FiniteSet
+from sympy.sets import Interval, FiniteSet, Union, ImageSet
 from sympy.testing.pytest import raises
 from sympy.abc import x, y
 
@@ -25,6 +29,17 @@ def test_singularities():
         FiniteSet(-1, (1 - sqrt(3) * I) / 2, (1 + sqrt(3) * I) / 2)
     assert singularities(1/(y**2 + 2*I*y + 1), y) == \
         FiniteSet(-I + sqrt(2)*I, -I - sqrt(2)*I)
+    _n = Dummy('n')
+    assert singularities(sech(x), x).dummy_eq(Union(
+        ImageSet(Lambda(_n, I*(2*_n*pi + pi/2)), S.Integers),
+        ImageSet(Lambda(_n, I*(2*_n*pi - pi/2)), S.Integers)))
+    assert singularities(coth(x), x).dummy_eq(Union(
+        ImageSet(Lambda(_n, I*(2*_n*pi + pi)), S.Integers),
+        ImageSet(Lambda(_n, 2*_n*I*pi), S.Integers)))
+    assert singularities(atanh(x), x) == FiniteSet(-1, 1)
+    assert singularities(acoth(x), x) == FiniteSet(-1, 1)
+    assert singularities(asech(x), x) == FiniteSet(0)
+    assert singularities(acsch(x), x) == FiniteSet(0)
 
     x = Symbol('x', real=True)
     assert singularities(1/(x**2 + 1), x) == S.EmptySet
@@ -32,6 +47,10 @@ def test_singularities():
     assert singularities(exp(1/x), x, Interval(1, 2)) == S.EmptySet
     assert singularities(log((x - 2)**2), x, Interval(1, 3)) == FiniteSet(2)
     raises(NotImplementedError, lambda: singularities(x**-oo, x))
+    assert singularities(sec(x), x, Interval(0, 3*pi)) == FiniteSet(
+        pi/2, 3*pi/2, 5*pi/2)
+    assert singularities(csc(x), x, Interval(0, 3*pi)) == FiniteSet(
+        0, pi, 2*pi, 3*pi)
 
 
 def test_is_increasing():
diff --git a/sympy/calculus/tests/test_util.py b/sympy/calculus/tests/test_util.py
index 6367fb98b1..02c2db89aa 100644
--- a/sympy/calculus/tests/test_util.py
+++ b/sympy/calculus/tests/test_util.py
@@ -1,11 +1,18 @@
+from sympy.core.function import Lambda
 from sympy.core.numbers import (E, I, Rational, oo, pi)
+from sympy.core.relational import Eq
 from sympy.core.singleton import S
-from sympy.core.symbol import (Symbol, symbols)
+from sympy.core.symbol import (Dummy, Symbol, symbols)
 from sympy.functions.elementary.complexes import (Abs, re)
 from sympy.functions.elementary.exponential import (exp, log)
+from sympy.functions.elementary.integers import frac
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.piecewise import Piecewise
-from sympy.functions.elementary.trigonometric import (cos, cot, csc, sec, sin, tan)
+from sympy.functions.elementary.trigonometric import (
+    cos, cot, csc, sec, sin, tan, asin, acos, atan, acot, asec, acsc)
+from sympy.functions.elementary.hyperbolic import (sinh, cosh, tanh, coth,
+    sech, csch, asinh, acosh, atanh, acoth, asech, acsch)
+from sympy.functions.special.gamma_functions import gamma
 from sympy.functions.special.error_functions import expint
 from sympy.matrices.expressions.matexpr import MatrixSymbol
 from sympy.simplify.simplify import simplify
@@ -13,7 +20,9 @@
                                  periodicity, lcim, is_convex,
                                  stationary_points, minimum, maximum)
 from sympy.sets.sets import (Interval, FiniteSet, Complement, Union)
-from sympy.testing.pytest import raises, _both_exp_pow
+from sympy.sets.fancysets import ImageSet
+from sympy.sets.conditionset import ConditionSet
+from sympy.testing.pytest import XFAIL, raises, _both_exp_pow
 from sympy.abc import x
 
 a = Symbol('a', real=True)
@@ -64,6 +73,8 @@ def test_continuous_domain():
     assert continuous_domain(tan(x), x, Interval(0, 2*pi)) == \
         Union(Interval(0, pi/2, False, True), Interval(pi/2, pi*Rational(3, 2), True, True),
               Interval(pi*Rational(3, 2), 2*pi, True, False))
+    assert continuous_domain(cot(x), x, Interval(0, 2*pi)) == Union(
+        Interval.open(0, pi), Interval.open(pi, 2*pi))
     assert continuous_domain((x - 1)/((x - 1)**2), x, S.Reals) == \
         Union(Interval(-oo, 1, True, True), Interval(1, oo, True, True))
     assert continuous_domain(log(x) + log(4*x - 1), x, S.Reals) == \
@@ -73,11 +84,66 @@ def test_continuous_domain():
         Union(Interval.open(-oo, 0), Interval.open(0, oo))
     assert continuous_domain(1/(x**2 - 4) + 2, x, S.Reals) == \
         Union(Interval.open(-oo, -2), Interval.open(-2, 2), Interval.open(2, oo))
+    assert continuous_domain((x+1)**pi, x, S.Reals) == Interval(-1, oo)
+    assert continuous_domain((x+1)**(pi/2), x, S.Reals) == Interval(-1, oo)
+    assert continuous_domain(x**x, x, S.Reals) == Interval(0, oo)
+    assert continuous_domain((x+1)**log(x**2), x, S.Reals) == Union(
+        Interval.Ropen(-1, 0), Interval.open(0, oo))
     domain = continuous_domain(log(tan(x)**2 + 1), x, S.Reals)
     assert not domain.contains(3*pi/2)
     assert domain.contains(5)
     d = Symbol('d', even=True, zero=False)
     assert continuous_domain(x**(1/d), x, S.Reals) == Interval(0, oo)
+    n = Dummy('n')
+    assert continuous_domain(1/sin(x), x, S.Reals).dummy_eq(Complement(
+        S.Reals, Union(ImageSet(Lambda(n, 2*n*pi + pi), S.Integers),
+                       ImageSet(Lambda(n, 2*n*pi), S.Integers))))
+    assert continuous_domain(sin(x) + cos(x), x, S.Reals) == S.Reals
+    assert continuous_domain(asin(x), x, S.Reals) == Interval(-1, 1) # issue #21786
+    assert continuous_domain(1/acos(log(x)), x, S.Reals) == Interval.Ropen(exp(-1), E)
+    assert continuous_domain(sinh(x)+cosh(x), x, S.Reals) == S.Reals
+    assert continuous_domain(tanh(x)+sech(x), x, S.Reals) == S.Reals
+    assert continuous_domain(atan(x)+asinh(x), x, S.Reals) == S.Reals
+    assert continuous_domain(acosh(x), x, S.Reals) == Interval(1, oo)
+    assert continuous_domain(atanh(x), x, S.Reals) == Interval.open(-1, 1)
+    assert continuous_domain(atanh(x)+acosh(x), x, S.Reals) == S.EmptySet
+    assert continuous_domain(asech(x), x, S.Reals) == Interval.Lopen(0, 1)
+    assert continuous_domain(acoth(x), x, S.Reals) == Union(
+        Interval.open(-oo, -1), Interval.open(1, oo))
+    assert continuous_domain(asec(x), x, S.Reals) == Union(
+        Interval(-oo, -1), Interval(1, oo))
+    assert continuous_domain(acsc(x), x, S.Reals) == Union(
+        Interval(-oo, -1), Interval(1, oo))
+    for f in (coth, acsch, csch):
+        assert continuous_domain(f(x), x, S.Reals) == Union(
+            Interval.open(-oo, 0), Interval.open(0, oo))
+    assert continuous_domain(acot(x), x, S.Reals).contains(0) == False
+    assert continuous_domain(1/(exp(x) - x), x, S.Reals) == Complement(
+        S.Reals, ConditionSet(x, Eq(-x + exp(x), 0), S.Reals))
+    assert continuous_domain(frac(x**2), x, Interval(-2,-1)) == Union(
+        Interval.open(-2, -sqrt(3)), Interval.open(-sqrt(2), -1),
+        Interval.open(-sqrt(3), -sqrt(2)))
+    assert continuous_domain(frac(x), x, S.Reals) == Complement(
+        S.Reals, S.Integers)
+    raises(NotImplementedError, lambda : continuous_domain(
+        1/(x**2+1), x, S.Complexes))
+    raises(NotImplementedError, lambda : continuous_domain(
+        gamma(x), x, Interval(-5,0)))
+    assert continuous_domain(x + gamma(pi), x, S.Reals) == S.Reals
+
+
+@XFAIL
+def test_continuous_domain_acot():
+    acot_cont = Piecewise((pi+acot(x), x<0), (acot(x), True))
+    assert continuous_domain(acot_cont, x, S.Reals) == S.Reals
+
+@XFAIL
+def test_continuous_domain_gamma():
+    assert continuous_domain(gamma(x), x, S.Reals).contains(-1) == False
+
+@XFAIL
+def test_continuous_domain_neg_power():
+    assert continuous_domain((x-2)**(1-x), x, S.Reals) == Interval.open(2, oo)
 
 
 def test_not_empty_in():
@@ -329,3 +395,8 @@ def test_issue_16469():
 @_both_exp_pow
 def test_issue_18747():
     assert periodicity(exp(pi*I*(x/4+S.Half/2)), x) == 8
+
+
+def test_issue_25942():
+    x = Symbol("x")
+    assert (acos(x) > pi/3).as_set() == Interval.Ropen(-1, S(1)/2)
diff --git a/sympy/plotting/tests/test_series.py b/sympy/plotting/tests/test_series.py
index d362737281..e23aa71915 100644
--- a/sympy/plotting/tests/test_series.py
+++ b/sympy/plotting/tests/test_series.py
@@ -90,7 +90,8 @@ def test_detect_poles():
         assert not np.any(np.isnan(yy1))
         assert np.any(np.isnan(yy2)) and np.any(np.isnan(yy2))
         assert not np.allclose(yy1, yy2, equal_nan=True)
-        assert len(s3.poles_locations) == 0
+        # The poles below are actually step discontinuities.
+        assert len(s3.poles_locations) == 21
 
     s1 = LineOver1DRangeSeries(tan(u * x), (x, -pi, pi), params={u: 1},
         adaptive=False, n=1000, detect_poles=False)