(D, N, t=symbols("t", integer=True))
Solves the equation `x^2 - Dy^2 = N`.

    Mainly concerned in the case `D > 0, D` is not a perfect square, which is
    the same as generalized Pell equation. To solve the generalized Pell
    equation this function Uses LMM algorithm. Refer [1]_ for more details on
    the algorithm.
    Returns one solution for each class of the solutions. Other solutions of
    the class can be constructed according to the values of ``D`` and ``N``.
    Returns a list containing the solution tuples `(x, y)`.

    Usage
    =====

    ``diop_DN(D, N, t)``: D and N are integers as in `x^2 - Dy^2 = N` and
    ``t`` is the parameter to be used in the solutions.

    Details
    =======

    ``D`` and ``N`` correspond to D and N in the equation.
    ``t`` is the parameter to be used in the solutions.

    Examples
    ========

    >>> from sympy.solvers.diophantine import diop_DN
    >>> diop_DN(13, -4) # Solves equation x**2 - 13*y**2 = -4
    [(3, 1), (393, 109), (36, 10)]

    The output can be interpreted as follows: There are three fundamental
    solutions to the equation `x^2 - 13y^2 = -4` given by (3, 1), (393, 109)
    and (36, 10). Each tuple is in the form (x, y), i. e solution (3, 1) means
    that `x = 3` and `y = 1`.

    >>> diop_DN(986, 1) # Solves equation x**2 - 986*y**2 = 1
    [(49299, 1570)]

    See Also
    ========

    find_DN(), diop_bf_DN()

    References
    ==========

    .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
        Robertson, July 31, 2004, Pages 16 - 17. [online], Available:
        http://www.jpr2718.org/pell.pdf
Solves the equation `x^2 - Dy^2 = N`.
Solves the equation x^2 - Dy^2 = N .
def diop_DN(D, N, t=symbols("t", integer=True)):
    """
    Solves the equation `x^2 - Dy^2 = N`.

    Mainly concerned in the case `D > 0, D` is not a perfect square, which is
    the same as generalized Pell equation. To solve the generalized Pell
    equation this function Uses LMM algorithm. Refer [1]_ for more details on
    the algorithm.
    Returns one solution for each class of the solutions. Other solutions of
    the class can be constructed according to the values of ``D`` and ``N``.
    Returns a list containing the solution tuples `(x, y)`.

    Usage
    =====

    ``diop_DN(D, N, t)``: D and N are integers as in `x^2 - Dy^2 = N` and
    ``t`` is the parameter to be used in the solutions.

    Details
    =======

    ``D`` and ``N`` correspond to D and N in the equation.
    ``t`` is the parameter to be used in the solutions.

    Examples
    ========

    >>> from sympy.solvers.diophantine import diop_DN
    >>> diop_DN(13, -4) # Solves equation x**2 - 13*y**2 = -4
    [(3, 1), (393, 109), (36, 10)]

    The output can be interpreted as follows: There are three fundamental
    solutions to the equation `x^2 - 13y^2 = -4` given by (3, 1), (393, 109)
    and (36, 10). Each tuple is in the form (x, y), i. e solution (3, 1) means
    that `x = 3` and `y = 1`.

    >>> diop_DN(986, 1) # Solves equation x**2 - 986*y**2 = 1
    [(49299, 1570)]

    See Also
    ========

    find_DN(), diop_bf_DN()

    References
    ==========

    .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
        Robertson, July 31, 2004, Pages 16 - 17. [online], Available:
        http://www.jpr2718.org/pell.pdf
    """
    if D < 0:
        if N == 0:
            return [(S.Zero, S.Zero)]
        elif N < 0:
            return []
        elif N > 0:
            d = divisors(square_factor(N))
            sol = []

            for divisor in d:
                sols = cornacchia(1, -D, N // divisor**2)
                if sols:
                    for x, y in sols:
                        sol.append((divisor*x, divisor*y))

            return sol

    elif D == 0:
        if N < 0 or not isinstance(sqrt(N), Integer):
            return []
        if N == 0:
            return [(S.Zero, t)]
        if isinstance(sqrt(N), Integer):
            return [(sqrt(N), t)]

    else: # D > 0
        if isinstance(sqrt(D), Integer):
            r = sqrt(D)

            if N == 0:
                return [(r*t, t)]
            else:
                sol = []

                for y in xrange(floor(sign(N)*(N - 1)/(2*r)) + 1):
                    if isinstance(sqrt(D*y**2 + N), Integer):
                        sol.append((sqrt(D*y**2 + N), y))

                return sol
        else:
            if N == 0:
                return [(S.Zero, S.Zero)]

            elif abs(N) == 1:

                pqa = PQa(0, 1, D)
                a_0 = floor(sqrt(D))
                l = 0
                G = []
                B = []

                for i in pqa:

                    a = i[2]
                    G.append(i[5])
                    B.append(i[4])

                    if l != 0 and a == 2*a_0:
                        break
                    l = l + 1

                if l % 2 == 1:

                    if N == -1:
                        x = G[l-1]
                        y = B[l-1]
                    else:
                        count = l
                        while count < 2*l - 1:
                            i = next(pqa)
                            G.append(i[5])
                            B.append(i[4])
                            count = count + 1

                        x = G[count]
                        y = B[count]
                else:
                    if N == 1:
                        x = G[l-1]
                        y = B[l-1]
                    else:
                        return []

                return [(x, y)]

            else:

                fs = []
                sol = []
                div = divisors(N)

                for d in div:
                    if divisible(N, d**2):
                        fs.append(d)

                for f in fs:
                    m = N // f**2
                    zs = sqrt_mod(D, abs(m), True)

                    zs = [i for i in zs if i <= abs(m) // 2 ]
                    if abs(m) != 2:
                        zs = zs + [-i for i in zs]
                        if S.Zero in zs:
                            zs.remove(S.Zero) # Remove duplicate zero

                    for z in zs:

                        pqa = PQa(z, abs(m), D)
                        l = 0
                        G = []
                        B = []

                        for i in pqa:

                            a = i[2]
                            G.append(i[5])
                            B.append(i[4])

                            if l != 0 and abs(i[1]) == 1:
                                r = G[l-1]
                                s = B[l-1]

                                if r**2 - D*s**2 == m:
                                    sol.append((f*r, f*s))

                                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

                            l = l + 1
                            if l == length(z, abs(m), D):
                                break

                return sol
def diop_DN ( D , N , t = symbols ( "t" , integer = True ) ) : if D < 0 : if N == 0 : return [ ( S . Zero , S . Zero ) ] elif N < 0 : return [ ] elif N > 0 : d = divisors ( square_factor ( N ) ) sol = [ ] for divisor in d : sols = cornacchia ( 1 , - D , N // divisor ** 2 ) if sols : for x , y in sols : sol . append ( ( divisor * x , divisor * y ) ) return sol elif D == 0 : if N < 0 or not isinstance ( sqrt ( N ) , Integer ) : return [ ] if N == 0 : return [ ( S . Zero , t ) ] if isinstance ( sqrt ( N ) , Integer ) : return [ ( sqrt ( N ) , t ) ] else : # D > 0 if isinstance ( sqrt ( D ) , Integer ) : r = sqrt ( D ) if N == 0 : return [ ( r * t , t ) ] else : sol = [ ] for y in xrange ( floor ( sign ( N ) * ( N - 1 ) / ( 2 * r ) ) + 1 ) : if isinstance ( sqrt ( D * y ** 2 + N ) , Integer ) : sol . append ( ( sqrt ( D * y ** 2 + N ) , y ) ) return sol else : if N == 0 : return [ ( S . Zero , S . Zero ) ] elif abs ( N ) == 1 : pqa = PQa ( 0 , 1 , D ) a_0 = floor ( sqrt ( D ) ) l = 0 G = [ ] B = [ ] for i in pqa : a = i [ 2 ] G . append ( i [ 5 ] ) B . append ( i [ 4 ] ) if l != 0 and a == 2 * a_0 : break l = l + 1 if l % 2 == 1 : if N == - 1 : x = G [ l - 1 ] y = B [ l - 1 ] else : count = l while count < 2 * l - 1 : i = next ( pqa ) G . append ( i [ 5 ] ) B . append ( i [ 4 ] ) count = count + 1 x = G [ count ] y = B [ count ] else : if N == 1 : x = G [ l - 1 ] y = B [ l - 1 ] else : return [ ] return [ ( x , y ) ] else : fs = [ ] sol = [ ] div = divisors ( N ) for d in div : if divisible ( N , d ** 2 ) : fs . append ( d ) for f in fs : m = N // f ** 2 zs = sqrt_mod ( D , abs ( m ) , True ) zs = [ i for i in zs if i <= abs ( m ) // 2 ] if abs ( m ) != 2 : zs = zs + [ - i for i in zs ] if S . Zero in zs : zs . remove ( S . Zero ) # Remove duplicate zero for z in zs : pqa = PQa ( z , abs ( m ) , D ) l = 0 G = [ ] B = [ ] for i in pqa : a = i [ 2 ] G . append ( i [ 5 ] ) B . append ( i [ 4 ] ) if l != 0 and abs ( i [ 1 ] ) == 1 : r = G [ l - 1 ] s = B [ l - 1 ] if r ** 2 - D * s ** 2 == m : sol . append ( ( f * r , f * s ) ) 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 l = l + 1 if l == length ( z , abs ( m ) , D ) : break return sol