Repository: sympy/sympy

Issue:
make diophantine output canonical wrt sign
```python
>>> diophantine(y+5-x**2)
{(-t, t**2 - 5)}
```
The solution for `x` might as well be `t` since it appears as a square in the equation.

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:
* solvers
  * DiophantineSolutionSet now makes solution tuples canonical wrt sign ([#25830](https://github.com/sympy/sympy/pull/25830) by [@smichr](https://github.com/smichr))

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
    <!-- 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. -->

    fixes #25775

    #### 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 -->
    * solvers
        * DiophantineSolutionSet now makes solution tuples canonical wrt sign
    <!-- END RELEASE NOTES -->


</details><p>


Patch:
diff --git a/sympy/external/gmpy.py b/sympy/external/gmpy.py
index a4551355ca..a8b2a43bc9 100644
--- a/sympy/external/gmpy.py
+++ b/sympy/external/gmpy.py
@@ -26,13 +26,6 @@
     is_fermat_prp as python_is_fermat_prp,
     is_euler_prp as python_is_euler_prp,
     is_strong_prp as python_is_strong_prp,
-    is_fibonacci_prp as python_is_fibonacci_prp,
-    is_lucas_prp as python_is_lucas_prp,
-    is_selfridge_prp as python_is_selfridge_prp,
-    is_strong_lucas_prp as python_is_strong_lucas_prp,
-    is_strong_selfridge_prp as python_is_strong_selfridge_prp,
-    is_bpsw_prp as python_is_bpsw_prp,
-    is_strong_bpsw_prp as python_is_strong_bpsw_prp,
 )
 
 
@@ -77,13 +70,6 @@
     'is_fermat_prp',
     'is_euler_prp',
     'is_strong_prp',
-    'is_fibonacci_prp',
-    'is_lucas_prp',
-    'is_selfridge_prp',
-    'is_strong_lucas_prp',
-    'is_strong_selfridge_prp',
-    'is_bpsw_prp',
-    'is_strong_bpsw_prp',
 ]
 
 
@@ -194,13 +180,6 @@ def iroot(x, n):
     is_fermat_prp = gmpy.is_fermat_prp
     is_euler_prp = gmpy.is_euler_prp
     is_strong_prp = gmpy.is_strong_prp
-    is_fibonacci_prp = gmpy.is_fibonacci_prp
-    is_lucas_prp = gmpy.is_lucas_prp
-    is_selfridge_prp = gmpy.is_selfridge_prp
-    is_strong_lucas_prp = gmpy.is_strong_lucas_prp
-    is_strong_selfridge_prp = gmpy.is_strong_selfridge_prp
-    is_bpsw_prp = gmpy.is_bpsw_prp
-    is_strong_bpsw_prp = gmpy.is_strong_bpsw_prp
 
 elif GROUND_TYPES == 'flint':
 
@@ -252,13 +231,6 @@ def iroot(x, n):
     is_fermat_prp = python_is_fermat_prp
     is_euler_prp = python_is_euler_prp
     is_strong_prp = python_is_strong_prp
-    is_fibonacci_prp = python_is_fibonacci_prp
-    is_lucas_prp = python_is_lucas_prp
-    is_selfridge_prp = python_is_selfridge_prp
-    is_strong_lucas_prp = python_is_strong_lucas_prp
-    is_strong_selfridge_prp = python_is_strong_selfridge_prp
-    is_bpsw_prp = python_is_bpsw_prp
-    is_strong_bpsw_prp = python_is_strong_bpsw_prp
 
 elif GROUND_TYPES == 'python':
 
@@ -286,13 +258,6 @@ def iroot(x, n):
     is_fermat_prp = python_is_fermat_prp
     is_euler_prp = python_is_euler_prp
     is_strong_prp = python_is_strong_prp
-    is_fibonacci_prp = python_is_fibonacci_prp
-    is_lucas_prp = python_is_lucas_prp
-    is_selfridge_prp = python_is_selfridge_prp
-    is_strong_lucas_prp = python_is_strong_lucas_prp
-    is_strong_selfridge_prp = python_is_strong_selfridge_prp
-    is_bpsw_prp = python_is_bpsw_prp
-    is_strong_bpsw_prp = python_is_strong_bpsw_prp
 
 else:
     assert False
diff --git a/sympy/external/ntheory.py b/sympy/external/ntheory.py
index 296b1331a9..f7bcde8fd9 100644
--- a/sympy/external/ntheory.py
+++ b/sympy/external/ntheory.py
@@ -324,20 +324,6 @@ def is_euler_prp(n, a):
     return pow(a, n >> 1, n) == jacobi(a, n) % n
 
 
-def _is_strong_prp(n, a):
-    s = bit_scan1(n - 1)
-    a = pow(a, n >> s, n)
-    if a == 1 or a == n - 1:
-        return True
-    for _ in range(s - 1):
-        a = pow(a, 2, n)
-        if a == n - 1:
-            return True
-        if a == 1:
-            return False
-    return False
-
-
 def is_strong_prp(n, a):
     if a < 2:
         raise ValueError("is_strong_prp() requires 'a' greater than or equal to 2")
@@ -350,288 +336,14 @@ def is_strong_prp(n, a):
     a %= n
     if gcd(n, a) != 1:
         raise ValueError("is_strong_prp() requires gcd(n,a) == 1")
-    return _is_strong_prp(n, a)
-
-
-def _lucas_sequence(n, P, Q, k):
-    r"""Return the modular Lucas sequence (U_k, V_k, Q_k).
-
-    Explanation
-    ===========
-
-    Given a Lucas sequence defined by P, Q, returns the kth values for
-    U and V, along with Q^k, all modulo n. This is intended for use with
-    possibly very large values of n and k, where the combinatorial functions
-    would be completely unusable.
-
-    .. math ::
-        U_k = \begin{cases}
-             0 & \text{if } k = 0\\
-             1 & \text{if } k = 1\\
-             PU_{k-1} - QU_{k-2} & \text{if } k > 1
-        \end{cases}\\
-        V_k = \begin{cases}
-             2 & \text{if } k = 0\\
-             P & \text{if } k = 1\\
-             PV_{k-1} - QV_{k-2} & \text{if } k > 1
-        \end{cases}
-
-    The modular Lucas sequences are used in numerous places in number theory,
-    especially in the Lucas compositeness tests and the various n + 1 proofs.
-
-    Parameters
-    ==========
-
-    n : int
-        n is an odd number greater than or equal to 3
-    P : int
-    Q : int
-        D determined by D = P**2 - 4*Q is non-zero
-    k : int
-        k is a nonnegative integer
-
-    Returns
-    =======
-
-    U, V, Qk : (int, int, int)
-        `(U_k \bmod{n}, V_k \bmod{n}, Q^k \bmod{n})`
-
-    Examples
-    ========
-
-    >>> from sympy.external.ntheory import _lucas_sequence
-    >>> N = 10**2000 + 4561
-    >>> sol = U, V, Qk = _lucas_sequence(N, 3, 1, N//2); sol
-    (0, 2, 1)
-
-    References
-    ==========
-
-    .. [1] https://en.wikipedia.org/wiki/Lucas_sequence
-
-    """
-    if k == 0:
-        return (0, 2, 1)
-    D = P**2 - 4*Q
-    U = 1
-    V = P
-    Qk = Q % n
-    if Q == 1:
-        # Optimization for extra strong tests.
-        for b in bin(k)[3:]:
-            U = (U*V) % n
-            V = (V*V - 2) % n
-            if b == "1":
-                U, V = U*P + V, V*P + U*D
-                if U & 1:
-                    U += n
-                if V & 1:
-                    V += n
-                U, V = U >> 1, V >> 1
-    elif P == 1 and Q == -1:
-        # Small optimization for 50% of Selfridge parameters.
-        for b in bin(k)[3:]:
-            U = (U*V) % n
-            if Qk == 1:
-                V = (V*V - 2) % n
-            else:
-                V = (V*V + 2) % n
-                Qk = 1
-            if b == "1":
-                # new_U = (U + V) // 2
-                # new_V = (5*U + V) // 2 = 2*U + new_U
-                U, V  = U + V, U << 1
-                if U & 1:
-                    U += n
-                U >>= 1
-                V += U
-                Qk = -1
-        Qk %= n
-    elif P == 1:
-        for b in bin(k)[3:]:
-            U = (U*V) % n
-            V = (V*V - 2*Qk) % n
-            Qk *= Qk
-            if b == "1":
-                # new_U = (U + V) // 2
-                # new_V = new_U - 2*Q*U
-                U, V  = U + V, (Q*U) << 1
-                if U & 1:
-                    U += n
-                U >>= 1
-                V = U - V
-                Qk *= Q
-            Qk %= n
-    else:
-        # The general case with any P and Q.
-        for b in bin(k)[3:]:
-            U = (U*V) % n
-            V = (V*V - 2*Qk) % n
-            Qk *= Qk
-            if b == "1":
-                U, V = U*P + V, V*P + U*D
-                if U & 1:
-                    U += n
-                if V & 1:
-                    V += n
-                U, V = U >> 1, V >> 1
-                Qk *= Q
-            Qk %= n
-    return (U % n, V % n, Qk)
-
-
-def is_fibonacci_prp(n, p, q):
-    d = p**2 - 4*q
-    if d == 0 or p <= 0 or q not in [1, -1]:
-        raise ValueError("invalid values for p,q in is_fibonacci_prp()")
-    if n < 1:
-        raise ValueError("is_fibonacci_prp() requires 'n' be greater than 0")
-    if n == 1:
-        return False
-    if n % 2 == 0:
-        return n == 2
-    return _lucas_sequence(n, p, q, n)[1] == p % n
-
-
-def is_lucas_prp(n, p, q):
-    d = p**2 - 4*q
-    if d == 0:
-        raise ValueError("invalid values for p,q in is_lucas_prp()")
-    if n < 1:
-        raise ValueError("is_lucas_prp() requires 'n' be greater than 0")
-    if n == 1:
-        return False
-    if n % 2 == 0:
-        return n == 2
-    if gcd(n, q*d) not in [1, n]:
-        raise ValueError("is_lucas_prp() requires gcd(n,2*q*D) == 1")
-    return _lucas_sequence(n, p, q, n - jacobi(d, n))[0] == 0
-
-
-def _is_selfridge_prp(n):
-    """Lucas compositeness test with the Selfridge parameters for n.
-
-    Explanation
-    ===========
-
-    The Lucas compositeness test checks whether n is a prime number.
-    The test can be run with arbitrary parameters ``P`` and ``Q``, which also change the performance of the test.
-    So, which parameters are most effective for running the Lucas compositeness test?
-    As an algorithm for determining ``P`` and ``Q``, Selfridge proposed method A [1]_ page 1401
-    (Since two methods were proposed, referred to simply as A and B in the paper,
-    we will refer to one of them as "method A").
-
-    method A fixes ``P = 1``. Then, ``D`` defined by ``D = P**2 - 4Q`` is varied from 5, -7, 9, -11, 13, and so on,
-    with the first ``D`` being ``jacobi(D, n) == -1``. Once ``D`` is determined,
-    ``Q`` is determined to be ``(P**2 - D)//4``.
-
-    References
-    ==========
-
-    .. [1] Robert Baillie, Samuel S. Wagstaff, Lucas Pseudoprimes,
-           Math. Comp. Vol 35, Number 152 (1980), pp. 1391-1417,
-           https://doi.org/10.1090%2FS0025-5718-1980-0583518-6
-           http://mpqs.free.fr/LucasPseudoprimes.pdf
-
-    """
-    for D in range(5, 1_000_000, 2):
-        if D & 2: # if D % 4 == 3
-            D = -D
-        j = jacobi(D, n)
-        if j == -1:
-            return _lucas_sequence(n, 1, (1-D) // 4, n + 1)[0] == 0
-        if j == 0 and D % n:
-            return False
-        # When j == -1 is hard to find, suspect a square number
-        if D == 13 and is_square(n):
-            return False
-    raise ValueError("appropriate value for D cannot be found in is_selfridge_prp()")
-
-
-def is_selfridge_prp(n):
-    if n < 1:
-        raise ValueError("is_selfridge_prp() requires 'n' be greater than 0")
-    if n == 1:
-        return False
-    if n % 2 == 0:
-        return n == 2
-    return _is_selfridge_prp(n)
-
-
-def is_strong_lucas_prp(n, p, q):
-    D = p**2 - 4*q
-    if D == 0:
-        raise ValueError("invalid values for p,q in is_strong_lucas_prp()")
-    if n < 1:
-        raise ValueError("is_selfridge_prp() requires 'n' be greater than 0")
-    if n == 1:
-        return False
-    if n % 2 == 0:
-        return n == 2
-    if gcd(n, q*D) not in [1, n]:
-        raise ValueError("is_strong_lucas_prp() requires gcd(n,2*q*D) == 1")
-    j = jacobi(D, n)
-    s = bit_scan1(n - j)
-    U, V, Qk = _lucas_sequence(n, p, q, (n - j) >> s)
-    if U == 0 or V == 0:
+    s = bit_scan1(n - 1)
+    a = pow(a, n >> s, n)
+    if a == 1 or a == n - 1:
         return True
     for _ in range(s - 1):
-        V = (V*V - 2*Qk) % n
-        if V == 0:
+        a = pow(a, 2, n)
+        if a == n - 1:
             return True
-        Qk = pow(Qk, 2, n)
-    return False
-
-
-def _is_strong_selfridge_prp(n):
-    for D in range(5, 1_000_000, 2):
-        if D & 2: # if D % 4 == 3
-            D = -D
-        j = jacobi(D, n)
-        if j == -1:
-            s = bit_scan1(n + 1)
-            U, V, Qk = _lucas_sequence(n, 1, (1-D) // 4, (n + 1) >> s)
-            if U == 0 or V == 0:
-                return True
-            for _ in range(s - 1):
-                V = (V*V - 2*Qk) % n
-                if V == 0:
-                    return True
-                Qk = pow(Qk, 2, n)
-            return False
-        if j == 0 and D % n:
-            return False
-        # When j == -1 is hard to find, suspect a square number
-        if D == 13 and is_square(n):
+        if a == 1:
             return False
-    raise ValueError("appropriate value for D cannot be found in is_strong_selfridge_prp()")
-
-
-def is_strong_selfridge_prp(n):
-    if n < 1:
-        raise ValueError("is_strong_selfridge_prp() requires 'n' be greater than 0")
-    if n == 1:
-        return False
-    if n % 2 == 0:
-        return n == 2
-    return _is_strong_selfridge_prp(n)
-
-
-def is_bpsw_prp(n):
-    if n < 1:
-        raise ValueError("is_bpsw_prp() requires 'n' be greater than 0")
-    if n == 1:
-        return False
-    if n % 2 == 0:
-        return n == 2
-    return _is_strong_prp(n, 2) and _is_selfridge_prp(n)
-
-
-def is_strong_bpsw_prp(n):
-    if n < 1:
-        raise ValueError("is_strong_bpsw_prp() requires 'n' be greater than 0")
-    if n == 1:
-        return False
-    if n % 2 == 0:
-        return n == 2
-    return _is_strong_prp(n, 2) and _is_strong_selfridge_prp(n)
+    return False
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index e913717093..bda9a8d7bf 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -111,6 +111,12 @@ def __init__(self, symbols_seq, parameters):
     def add(self, solution):
         if len(solution) != len(self.symbols):
             raise ValueError("Solution should have a length of %s, not %s" % (len(self.symbols), len(solution)))
+        # make solution canonical wrt sign (i.e. no -x unless x is also present as an arg)
+        args = set(solution)
+        for i in range(len(solution)):
+            x = solution[i]
+            if not type(x) is int and (-x).is_Symbol and -x not in args:
+                solution = [_.subs(-x, x) for _ in solution]
         super().add(Tuple(*solution))
 
     def update(self, *solutions):
@@ -1322,7 +1328,7 @@ def diophantine(eq, param=symbols("t", integer=True), syms=None,
     >>> diophantine(x*(2*x + 3*y - z))
     {(0, n1, n2), (t_0, t_1, 2*t_0 + 3*t_1)}
     >>> diophantine(x**2 + 3*x*y + 4*x)
-    {(0, n1), (3*t_0 - 4, -t_0)}
+    {(0, n1), (-3*t_0 - 4, t_0)}
 
     See Also
     ========
diff --git a/sympy/stats/joint_rv_types.py b/sympy/stats/joint_rv_types.py
index b2bd86f918..c1a24f934e 100644
--- a/sympy/stats/joint_rv_types.py
+++ b/sympy/stats/joint_rv_types.py
@@ -633,7 +633,8 @@ def MultivariateEwens(syms, n, theta):
     ==========
 
     .. [1] https://en.wikipedia.org/wiki/Ewens%27s_sampling_formula
-    .. [2] https://www.jstor.org/stable/24780825
+    .. [2] https://www.researchgate.net/publication/280311472_The_Ubiquitous_Ewens_Sampling_Formula
+
     """
     return multivariate_rv(MultivariateEwensDistribution, syms, n, theta)
 


Test patch:
diff --git a/sympy/external/tests/test_ntheory.py b/sympy/external/tests/test_ntheory.py
index 00824481ad..437b4d81be 100644
--- a/sympy/external/tests/test_ntheory.py
+++ b/sympy/external/tests/test_ntheory.py
@@ -1,10 +1,7 @@
 from itertools import permutations
 
 from sympy.external.ntheory import (bit_scan1, remove, bit_scan0, is_fermat_prp,
-                                    is_euler_prp, is_strong_prp, gcdext, _lucas_sequence,
-                                    is_fibonacci_prp, is_lucas_prp, is_selfridge_prp,
-                                    is_strong_lucas_prp, is_strong_selfridge_prp,
-                                    is_bpsw_prp, is_strong_bpsw_prp)
+                                    is_euler_prp, is_strong_prp, gcdext)
 from sympy.testing.pytest import raises
 
 
@@ -144,164 +141,3 @@ def test_is_strong_prp():
                    16531, 18721, 19345, 23521, 31621, 44287, 47197]
     for n in pseudorpime:
         assert is_strong_prp(n, 3)
-
-
-def test_lucas_sequence():
-    def lucas_u(P, Q, length):
-        array = [0] * length
-        array[1] = 1
-        for k in range(2, length):
-            array[k] = P * array[k - 1] - Q * array[k - 2]
-        return array
-
-    def lucas_v(P, Q, length):
-        array = [0] * length
-        array[0] = 2
-        array[1] = P
-        for k in range(2, length):
-            array[k] = P * array[k - 1] - Q * array[k - 2]
-        return array
-
-    length = 20
-    for P in range(-10, 10):
-        for Q in range(-10, 10):
-            D = P**2 - 4*Q
-            if D == 0:
-                continue
-            us = lucas_u(P, Q, length)
-            vs = lucas_v(P, Q, length)
-            for n in range(3, 100, 2):
-                for k in range(length):
-                    U, V, Qk = _lucas_sequence(n, P, Q, k)
-                    assert U == us[k] % n
-                    assert V == vs[k] % n
-                    assert pow(Q, k, n) == Qk
-
-
-def test_is_fibonacci_prp():
-    # invalid input
-    raises(ValueError, lambda: is_fibonacci_prp(3, 2, 1))
-    raises(ValueError, lambda: is_fibonacci_prp(3, -5, 1))
-    raises(ValueError, lambda: is_fibonacci_prp(3, 5, 2))
-    raises(ValueError, lambda: is_fibonacci_prp(0, 5, -1))
-
-    # n = 1
-    assert not is_fibonacci_prp(1, 3, 1)
-
-    # n is prime
-    assert is_fibonacci_prp(2, 5, 1)
-    assert is_fibonacci_prp(3, 6, -1)
-    assert is_fibonacci_prp(11, 7, 1)
-    assert is_fibonacci_prp(2**31-1, 8, -1)
-
-    # A005845
-    pseudorpime = [705, 2465, 2737, 3745, 4181, 5777, 6721,
-                   10877, 13201, 15251, 24465, 29281, 34561]
-    for n in pseudorpime:
-        assert is_fibonacci_prp(n, 1, -1)
-
-
-def test_is_lucas_prp():
-    # invalid input
-    raises(ValueError, lambda: is_lucas_prp(3, 2, 1))
-    raises(ValueError, lambda: is_lucas_prp(0, 5, -1))
-    raises(ValueError, lambda: is_lucas_prp(15, 3, 1))
-
-    # n = 1
-    assert not is_lucas_prp(1, 3, 1)
-
-    # n is prime
-    assert is_lucas_prp(2, 5, 2)
-    assert is_lucas_prp(3, 6, -1)
-    assert is_lucas_prp(11, 7, 5)
-    assert is_lucas_prp(2**31-1, 8, -3)
-
-    # A081264
-    pseudorpime = [323, 377, 1891, 3827, 4181, 5777, 6601, 6721,
-                   8149, 10877, 11663, 13201, 13981, 15251, 17119]
-    for n in pseudorpime:
-        assert is_lucas_prp(n, 1, -1)
-
-
-def test_is_selfridge_prp():
-    # invalid input
-    raises(ValueError, lambda: is_selfridge_prp(0))
-
-    # n = 1
-    assert not is_selfridge_prp(1)
-
-    # n is prime
-    assert is_selfridge_prp(2)
-    assert is_selfridge_prp(3)
-    assert is_selfridge_prp(11)
-    assert is_selfridge_prp(2**31-1)
-
-    # A217120
-    pseudorpime = [323, 377, 1159, 1829, 3827, 5459, 5777, 9071,
-                   9179, 10877, 11419, 11663, 13919, 14839, 16109]
-    for n in pseudorpime:
-        assert is_selfridge_prp(n)
-
-
-def test_is_strong_lucas_prp():
-    # invalid input
-    raises(ValueError, lambda: is_strong_lucas_prp(3, 2, 1))
-    raises(ValueError, lambda: is_strong_lucas_prp(0, 5, -1))
-    raises(ValueError, lambda: is_strong_lucas_prp(15, 3, 1))
-
-    # n = 1
-    assert not is_strong_lucas_prp(1, 3, 1)
-
-    # n is prime
-    assert is_strong_lucas_prp(2, 5, 2)
-    assert is_strong_lucas_prp(3, 6, -1)
-    assert is_strong_lucas_prp(11, 7, 5)
-    assert is_strong_lucas_prp(2**31-1, 8, -3)
-
-
-def test_is_strong_selfridge_prp():
-    # invalid input
-    raises(ValueError, lambda: is_strong_selfridge_prp(0))
-
-    # n = 1
-    assert not is_strong_selfridge_prp(1)
-
-    # n is prime
-    assert is_strong_selfridge_prp(2)
-    assert is_strong_selfridge_prp(3)
-    assert is_strong_selfridge_prp(11)
-    assert is_strong_selfridge_prp(2**31-1)
-
-    # A217255
-    pseudorpime = [5459, 5777, 10877, 16109, 18971, 22499, 24569,
-                   25199, 40309, 58519, 75077, 97439, 100127, 113573]
-    for n in pseudorpime:
-        assert is_strong_selfridge_prp(n)
-
-
-def test_is_bpsw_prp():
-    # invalid input
-    raises(ValueError, lambda: is_bpsw_prp(0))
-
-    # n = 1
-    assert not is_bpsw_prp(1)
-
-    # n is prime
-    assert is_bpsw_prp(2)
-    assert is_bpsw_prp(3)
-    assert is_bpsw_prp(11)
-    assert is_bpsw_prp(2**31-1)
-
-
-def test_is_strong_bpsw_prp():
-    # invalid input
-    raises(ValueError, lambda: is_strong_bpsw_prp(0))
-
-    # n = 1
-    assert not is_strong_bpsw_prp(1)
-
-    # n is prime
-    assert is_strong_bpsw_prp(2)
-    assert is_strong_bpsw_prp(3)
-    assert is_strong_bpsw_prp(11)
-    assert is_strong_bpsw_prp(2**31-1)
diff --git a/sympy/ntheory/primetest.py b/sympy/ntheory/primetest.py
index 99b3f60501..c853ec6d6d 100644
--- a/sympy/ntheory/primetest.py
+++ b/sympy/ntheory/primetest.py
@@ -8,10 +8,7 @@
 from sympy.core.sympify import sympify
 from sympy.external.gmpy import (gmpy as _gmpy, gcd, jacobi,
                                  is_square as gmpy_is_square,
-                                 bit_scan1, is_fermat_prp, is_euler_prp,
-                                 is_selfridge_prp, is_strong_selfridge_prp,
-                                 is_strong_bpsw_prp)
-from sympy.external.ntheory import _lucas_sequence
+                                 bit_scan1, is_fermat_prp, is_euler_prp)
 from sympy.utilities.misc import as_int
 
 
@@ -257,6 +254,190 @@ def mr(n, bases):
     return True
 
 
+def _lucas_sequence(n, P, Q, k):
+    r"""Return the modular Lucas sequence (U_k, V_k, Q_k).
+
+    Explanation
+    ===========
+
+    Given a Lucas sequence defined by P, Q, returns the kth values for
+    U and V, along with Q^k, all modulo n. This is intended for use with
+    possibly very large values of n and k, where the combinatorial functions
+    would be completely unusable.
+
+    .. math ::
+        U_k = \begin{cases}
+             0 & \text{if } k = 0\\
+             1 & \text{if } k = 1\\
+             PU_{k-1} - QU_{k-2} & \text{if } k > 1
+        \end{cases}\\
+        V_k = \begin{cases}
+             2 & \text{if } k = 0\\
+             P & \text{if } k = 1\\
+             PV_{k-1} - QV_{k-2} & \text{if } k > 1
+        \end{cases}
+
+    The modular Lucas sequences are used in numerous places in number theory,
+    especially in the Lucas compositeness tests and the various n + 1 proofs.
+
+    Parameters
+    ==========
+
+    n : int
+        n is an odd number greater than or equal to 3
+    P : int
+    Q : int
+        D determined by D = P**2 - 4*Q is non-zero
+    k : int
+        k is a nonnegative integer
+
+    Returns
+    =======
+
+    U, V, Qk : (int, int, int)
+        `(U_k \bmod{n}, V_k \bmod{n}, Q^k \bmod{n})`
+
+    Examples
+    ========
+
+    >>> from sympy.ntheory.primetest import _lucas_sequence
+    >>> N = 10**2000 + 4561
+    >>> sol = U, V, Qk = _lucas_sequence(N, 3, 1, N//2); sol
+    (0, 2, 1)
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Lucas_sequence
+
+    """
+    if k == 0:
+        return (0, 2, 1)
+    D = P**2 - 4*Q
+    U = 1
+    V = P
+    Qk = Q % n
+    if Q == 1:
+        # Optimization for extra strong tests.
+        for b in bin(k)[3:]:
+            U = (U*V) % n
+            V = (V*V - 2) % n
+            if b == "1":
+                U, V = U*P + V, V*P + U*D
+                if U & 1:
+                    U += n
+                if V & 1:
+                    V += n
+                U, V = U >> 1, V >> 1
+    elif P == 1 and Q == -1:
+        # Small optimization for 50% of Selfridge parameters.
+        for b in bin(k)[3:]:
+            U = (U*V) % n
+            if Qk == 1:
+                V = (V*V - 2) % n
+            else:
+                V = (V*V + 2) % n
+                Qk = 1
+            if b == "1":
+                # new_U = (U + V) // 2
+                # new_V = (5*U + V) // 2 = 2*U + new_U
+                U, V  = U + V, U << 1
+                if U & 1:
+                    U += n
+                U >>= 1
+                V += U
+                Qk = -1
+        Qk %= n
+    elif P == 1:
+        for b in bin(k)[3:]:
+            U = (U*V) % n
+            V = (V*V - 2*Qk) % n
+            Qk *= Qk
+            if b == "1":
+                # new_U = (U + V) // 2
+                # new_V = new_U - 2*Q*U
+                U, V  = U + V, (Q*U) << 1
+                if U & 1:
+                    U += n
+                U >>= 1
+                V = U - V
+                Qk *= Q
+            Qk %= n
+    else:
+        # The general case with any P and Q.
+        for b in bin(k)[3:]:
+            U = (U*V) % n
+            V = (V*V - 2*Qk) % n
+            Qk *= Qk
+            if b == "1":
+                U, V = U*P + V, V*P + U*D
+                if U & 1:
+                    U += n
+                if V & 1:
+                    V += n
+                U, V = U >> 1, V >> 1
+                Qk *= Q
+            Qk %= n
+    return (U % n, V % n, Qk)
+
+
+def _lucas_selfridge_params(n):
+    """Calculates the Selfridge parameters (D, P, Q) for n.
+
+    Explanation
+    ===========
+
+    The Lucas compositeness test checks whether n is a prime number.
+    The test can be run with arbitrary parameters ``P`` and ``Q``, which also change the performance of the test.
+    So, which parameters are most effective for running the Lucas compositeness test?
+    As an algorithm for determining ``P`` and ``Q``, Selfridge proposed method A [1]_ page 1401
+    (Since two methods were proposed, referred to simply as A and B in the paper,
+    we will refer to one of them as "method A").
+
+    method A fixes ``P = 1``. Then, ``D`` defined by ``D = P**2 - 4Q`` is varied from 5, -7, 9, -11, 13, and so on,
+    with the first ``D`` being ``jacobi(D, n) == -1``. Once ``D`` is determined,
+    ``Q`` is determined to be ``(P**2 - D)//4``. Here is an implementation of the method A.
+
+    Parameters
+    ==========
+
+    n : int
+        positive odd integer
+
+    Returns
+    =======
+
+    D, P, Q: Selfridge parameters for Lucas compositeness test.
+             ``(0, 0, 0)`` if we find a nontrivial divisor of ``n``.
+
+    Examples
+    ========
+
+    >>> from sympy.ntheory.primetest import _lucas_selfridge_params
+    >>> _lucas_selfridge_params(101)
+    (-7, 1, 2)
+    >>> _lucas_selfridge_params(15)
+    (0, 0, 0)
+
+    References
+    ==========
+
+    .. [1] Robert Baillie, Samuel S. Wagstaff, Lucas Pseudoprimes,
+           Math. Comp. Vol 35, Number 152 (1980), pp. 1391-1417,
+           https://doi.org/10.1090%2FS0025-5718-1980-0583518-6
+           http://mpqs.free.fr/LucasPseudoprimes.pdf
+
+    """
+    for D in count(5, 2):
+        if D & 2: # if D % 4 == 3
+            D = -D
+        j = jacobi(D, n)
+        if j == -1:
+            return (D, 1, (1 - D)//4)
+        elif j == 0 and D % n:
+            return (0, 0, 0)
+
+
 def _lucas_extrastrong_params(n):
     """Calculates the "extra strong" parameters (D, P, Q) for n.
 
@@ -332,9 +513,17 @@ def is_lucas_prp(n):
     9179
     """
     n = as_int(n)
-    if n < 2:
+    if n == 2:
+        return True
+    if n < 2 or (n % 2) == 0:
+        return False
+    if gmpy_is_square(n):
+        return False
+
+    D, P, Q = _lucas_selfridge_params(n)
+    if D == 0:
         return False
-    return is_selfridge_prp(n)
+    return _lucas_sequence(n, P, Q, n + 1)[0] == 0
 
 
 def is_strong_lucas_prp(n):
@@ -370,9 +559,31 @@ def is_strong_lucas_prp(n):
     18971
     """
     n = as_int(n)
-    if n < 2:
+    if n == 2:
+        return True
+    if n < 2 or (n % 2) == 0:
         return False
-    return is_strong_selfridge_prp(n)
+    if gmpy_is_square(n):
+        return False
+
+    D, P, Q = _lucas_selfridge_params(n)
+    if D == 0:
+        return False
+
+    # remove powers of 2 from n+1 (= k * 2**s)
+    s = bit_scan1(n + 1)
+    k = (n + 1) >> s
+
+    U, V, Qk = _lucas_sequence(n, P, Q, k)
+
+    if U == 0 or V == 0:
+        return True
+    for _ in range(1, s):
+        V = (V*V - 2*Qk) % n
+        if V == 0:
+            return True
+        Qk = pow(Qk, 2, n)
+    return False
 
 
 def is_extra_strong_lucas_prp(n):
@@ -551,7 +762,7 @@ def isprime(n):
     # This should be a bit faster than our step 2, and for large values will
     # be a lot faster than our step 3 (C+GMP vs. Python).
     if _gmpy is not None:
-        return is_strong_bpsw_prp(n)
+        return _gmpy.is_strong_prp(n, 2) and _gmpy.is_strong_selfridge_prp(n)
 
 
     # Step 2: deterministic Miller-Rabin testing for numbers < 2^64.  See:
@@ -615,7 +826,7 @@ def isprime(n):
     #      3.2s   extra strong BPSW
 
     # Classic BPSW from page 1401 of the paper.  See alternate ideas below.
-    return is_strong_bpsw_prp(n)
+    return mr(n, [2]) and is_strong_lucas_prp(n)
 
     # Using extra strong test, which is somewhat faster
     #return mr(n, [2]) and is_extra_strong_lucas_prp(n)
diff --git a/sympy/ntheory/tests/test_primetest.py b/sympy/ntheory/tests/test_primetest.py
index 46061fe426..0f23f9a85f 100644
--- a/sympy/ntheory/tests/test_primetest.py
+++ b/sympy/ntheory/tests/test_primetest.py
@@ -1,7 +1,8 @@
 from math import gcd
 
 from sympy.ntheory.generate import Sieve, sieve
-from sympy.ntheory.primetest import (mr, _lucas_extrastrong_params, is_lucas_prp, is_square,
+from sympy.ntheory.primetest import (mr, _lucas_sequence, _lucas_selfridge_params, _lucas_extrastrong_params,
+                                     is_lucas_prp, is_square,
                                      is_strong_lucas_prp, is_extra_strong_lucas_prp, isprime, is_euler_pseudoprime,
                                      is_gaussian_prime, is_fermat_pseudoprime, is_euler_jacobi_pseudoprime)
 
@@ -53,6 +54,48 @@ def test_is_euler_jacobi_pseudoprime():
     assert is_euler_jacobi_pseudoprime(15, 1)
 
 
+def test_lucas_sequence():
+    def lucas_u(P, Q, length):
+        array = [0] * length
+        array[1] = 1
+        for k in range(2, length):
+            array[k] = P * array[k - 1] - Q * array[k - 2]
+        return array
+
+    def lucas_v(P, Q, length):
+        array = [0] * length
+        array[0] = 2
+        array[1] = P
+        for k in range(2, length):
+            array[k] = P * array[k - 1] - Q * array[k - 2]
+        return array
+
+    length = 20
+    for P in range(-10, 10):
+        for Q in range(-10, 10):
+            D = P**2 - 4*Q
+            if D == 0:
+                continue
+            us = lucas_u(P, Q, length)
+            vs = lucas_v(P, Q, length)
+            for n in range(3, 100, 2):
+                for k in range(length):
+                    U, V, Qk = _lucas_sequence(n, P, Q, k)
+                    assert U == us[k] % n
+                    assert V == vs[k] % n
+                    assert pow(Q, k, n) == Qk
+
+
+def test_lucas_selfridge_params():
+    assert _lucas_selfridge_params(3) == (5, 1, -1)
+    assert _lucas_selfridge_params(5) == (-7, 1, 2)
+    assert _lucas_selfridge_params(7) == (5, 1, -1)
+    assert _lucas_selfridge_params(9) == (0, 0, 0)
+    assert _lucas_selfridge_params(11) == (13, 1, -3)
+    assert _lucas_selfridge_params(19) == (-7, 1, 2)
+    assert _lucas_selfridge_params(29) == (-11, 1, 3)
+
+
 def test_lucas_extrastrong_params():
     assert _lucas_extrastrong_params(3) == (5, 3, 1)
     assert _lucas_extrastrong_params(5) == (12, 4, 1)
diff --git a/sympy/solvers/diophantine/tests/test_diophantine.py b/sympy/solvers/diophantine/tests/test_diophantine.py
index 067b68f93b..3b3fd3522b 100644
--- a/sympy/solvers/diophantine/tests/test_diophantine.py
+++ b/sympy/solvers/diophantine/tests/test_diophantine.py
@@ -96,7 +96,7 @@ def test_linear():
     assert diop_solve(2*x - 3*y - 5) == (3*t_0 - 5, 2*t_0 - 5)
     assert diop_solve(-2*x - 3*y - 5) == (3*t_0 + 5, -2*t_0 - 5)
     assert diop_solve(7*x + 5*y) == (5*t_0, -7*t_0)
-    assert diop_solve(2*x + 4*y) == (2*t_0, -t_0)
+    assert diop_solve(2*x + 4*y) == (-2*t_0, t_0)
     assert diop_solve(4*x + 6*y - 4) == (3*t_0 - 2, -2*t_0 + 2)
     assert diop_solve(4*x + 6*y - 3) == (None, None)
     assert diop_solve(0*x + 3*y - 4*z + 5) == (4*t_0 + 5, 3*t_0 + 5)
@@ -569,10 +569,10 @@ def test_diophantine():
            {(-3, -2), (-3, 2), (-2, -3), (-2, 3), (2, -3), (2, 3), (3, -2), (3, 2)}
 
     # issue 18122
-    assert check_solutions(x**2-y)
-    assert check_solutions(y**2-x)
-    assert diophantine((x**2-y), t) == {(t, t**2)}
-    assert diophantine((y**2-x), t) == {(t**2, -t)}
+    assert check_solutions(x**2 - y)
+    assert check_solutions(y**2 - x)
+    assert diophantine((x**2 - y), t) == {(t, t**2)}
+    assert diophantine((y**2 - x), t) == {(t**2, t)}
 
 
 def test_general_pythagorean():
@@ -769,8 +769,8 @@ def test_diopcoverage():
     eq = (2*x + y + 1)**2
     assert diop_solve(eq) == {(t_0, -2*t_0 - 1)}
     eq = 2*x**2 + 6*x*y + 12*x + 4*y**2 + 18*y + 18
-    assert diop_solve(eq) == {(t, -t - 3), (2*t - 3, -t)}
-    assert diop_quadratic(x + y**2 - 3) == {(-t**2 + 3, -t)}
+    assert diop_solve(eq) == {(t, -t - 3), (-2*t - 3, t)}
+    assert diop_quadratic(x + y**2 - 3) == {(-t**2 + 3, t)}
 
     assert diop_linear(x + y - 3) == (t_0, 3 - t_0)
 
@@ -1033,5 +1033,5 @@ def test_quadratic_parameter_passing():
     eq = -33*x*y + 3*y**2
     solution = BinaryQuadratic(eq).solve(parameters=[t, u])
     # test that parameters are passed all the way to the final solution
-    assert solution == {(t, 11*t), (-t, 22*t)}
+    assert solution == {(t, 11*t), (t, -22*t)}
     assert solution(0, 0) == {(0, 0)}