"""
Functions related to the congruent number problem for elliptic curves.

AUTHOR: William Stein, 2005-10-17

EXAMPLES:
    sage: attach "cong.sage"
    sage: cong_number_sets(101)
    ([], [])
    sage: is_conj_congruent_number(101)
    True
    sage: is_congruent_number(101)
    True
    sage: print conj_congruent_number_list(30)
    [5, 6, 7, 13, 14, 15, 20, 21, 22, 23, 24, 28, 29, 30]
    sage: cong_number_curve(101)
    Elliptic Curve defined by y^2  = x^3 - 10201*x over Rational Field
    sage: congruent_curve_gens(101)
    Entering qsieve::search: y^2 = 1x^3 + 0x^2 + -10201x^1 + 0
     ...
    sage: P = _[0]
    sage: congruent_triangle(P)
    (44538033219/1326635050,
     267980280100/44538033219,
     2015242462949760001961/59085715926389725950)
    sage: congruent_triangle(2*P)
    (-3808424574270625821973109905486908673845521/238144087377294983538196567899844505175900,
     -48105105650213586674715706715768590045531800/3808424574270625821973109905486908673845521,
     18482628628473862825682187406224713731933103369206814890394108346849598187226371761441/906953794584941364348309838871418957694907021279988390325859242498283526441532143900)
"""

def cong_number_sets(n):
    """
    Given a positive integer n, returns the two sets appearing in the
    conjectural criterion for when a number is congruent.
    """
    n = ZZ(n)
    n = ZZ(prod([p for p, e in n.factor() if e%2 == 1]))
    
    if n % 2 == 0:   # even case
        E = []       # with c even
        O = []       # with c odd
        a_bound = (sqrt(n/8) + 1).integer_part()
        c_bound = (sqrt(n)/4 + 1).integer_part()
        half_n = n//2
        for c in range(-c_bound, c_bound+1):
            c_square = c^2
            for a in range(-a_bound, a_bound+1):
                a_square = a^2
                z = half_n - 4*a_square - 8*c_square
                try:
                    b = z.square_root()
                    if b.denominator() == 1:
                        b = b.numerator()
                        assert 4*a^2 + b^2 + 8*c^2 == n/2
                        if c % 2 == 0:
                            E.append((a,b,c))
                        else:
                            O.append((a,b,c))
                except ValueError:
                    pass
    else:
        E = []       # with c even
        O = []       # with c odd
        a_bound = (sqrt(n/2)+1).integer_part()
        c_bound = (sqrt(n/8) + 1).integer_part()
        for c in range(-c_bound, c_bound+1):
            c_square = c^2
            for a in range(-a_bound, a_bound+1):
                a_square = a^2
                z = n - 2*a_square - 8*c_square
                try:
                    b = z.square_root()
                    if b.denominator() == 1:
                        b = b.numerator()
                        assert 2*a^2 + b^2 + 8*c^2 == n
                        if c % 2 == 0:
                            E.append((a,b,c))
                        else:
                            O.append((a,b,c))
                except ValueError:
                    pass
    return E, O

def is_conj_congruent_number(n):
    """
    Returns True if n is conjecturally a congruent number, according
    to the Birch and Swinnerton-Dyer conjecture.
    """
    E, O = cong_number_sets(n)
    return len(E) == len(O)

def is_congruent_number(n):
    """
    Returns True if n is provably a congruent number.  This computes
    the rank of the corresponding elliptic curve.
    """
    return cong_number_curve(n).rank() > 0

def conj_congruent_number_list(bound):
    """
    Returns the conjectural congruent numbers up to a given bound.
    """
    return [n for n in range(1,bound+1) if is_conj_congruent_number(n)]

def cong_number_curve(n):
    """
    Returns the elliptic curve corresponding to the integer n.  This
    curve has positive rank if and only if n is a congruent number.
    """
    return EllipticCurve([-n^2,0])

def congruent_curve_gens(n, verbose=True):
    """
    Returns generators for the Mordell-Weil group of the elliptic
    curve corresponding to n.
    """
    E = cong_number_curve(n)
    return E.gens(verbose=verbose)

def congruent_triangle(P):
    """
    Returns the triangle corresponding to the point P on a congruent
    number elliptic curve.  The point P must have nonzero y coordinate.
    """
    m = P.curve().a_invariants()[3]
    n = (-m).square_root()
    (x, y) = (P[0],P[1])
    if y == 0:
        raise ArithmeticError, "point does not define a rational right triangle"
    T = ((n^2-x^2)/y, -2*x*n/y, (n^2+x^2)/y)
    assert (T[0]^2 + T[1]^2 == T[2]^2), "bug in program."
    return T

def point_from_triangle(a,b,c):
    """
    Return the point on $y^2=x^3-n^2x$ corresponding to the
    right triangle with side lengths $a,b,c$.
    """
    n = a*b/2
    E = EllipticCurve([-n^2, 0])
    return E([-n*b/(a+c), 2*n^2/(a+c)])