def charpoly_GF(n=100, p=16411, system='sage'):
    """
    Given a n x n matrix over GF with random entries, compute the
    charpoly.

    INPUT:

    - ``n`` - matrix dimension (default: 100)
    - ``p`` - prime number (default: ``16411``)
    - ``system`` - either 'magma' or 'sage' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.charpoly_GF(100)
        sage: tm = b.charpoly_GF(100, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(GF(p), n, n)
        t = cputime()
        v = A.charpoly()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := Random(MatrixAlgebra(GF(%s), n));
t := Cputime();
K := CharacteristicPolynomial(A);
s := Cputime(t);
"""%(n,p)
        if verbose: print code
        magma.eval(code)
        return magma.eval('s')
    else:
        raise ValueError('unknown system "%s"'%system)

def rank2_GF(n=500, p=16411, system='sage'):
    """
    Rank over GF(p): Given a (n + 10) x n matrix over GF(p) with
    random entries, compute the rank.

    INPUT:

    - ``n`` - matrix dimension (default: 300)
    - ``p`` - prime number (default: ``16411``)
    - ``system`` - either 'magma' or 'sage' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.rank2_GF(500)
        sage: tm = b.rank2_GF(500, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(GF(p), n+10, n)
        t = cputime()
        v = A.rank()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := Random(MatrixAlgebra(GF(%s), n));
t := Cputime();
K := Rank(A);
s := Cputime(t);
"""%(n,p)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def nullspace_GF(n=300, p=16411, system='sage'):
    """
    Given a n+1 x n  matrix over GF(p) with random
    entries, compute the nullspace.

    INPUT:

    - ``n`` - matrix dimension (default: 300)
    - ``p`` - prime number (default: ``16411``)
    - ``system`` - either 'magma' or 'sage' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.nullspace_GF(300)
        sage: tm = b.nullspace_GF(300, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(GF(p), n, n+1)
        t = cputime()
        v = A.kernel()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := Random(RMatrixSpace(GF(%s), n, n+1));
t := Cputime();
K := Kernel(A);
s := Cputime(t);
"""%(n,p)
        if verbose: print code
        magma.eval(code)
        return magma.eval('s')
    else:
        raise ValueError('unknown system "%s"'%system)

def invert_hilbert_QQ(n=40, system='sage'):
    """
    Runs the benchmark for calculating the inverse of the hilbert
    matrix over rationals of dimension n.

    INPUT:

    - ``n`` - matrix dimension (default: ``300``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.invert_hilbert_QQ(30)
        sage: tm = b.invert_hilbert_QQ(30, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = hilbert_matrix(n)
        t = cputime()
        d = A**(-1)
        return cputime(t)
    elif system == 'magma':
        code = """
h := HilbertMatrix(%s);
tinit := Cputime();
d := h^(-1);
s := Cputime(tinit);
delete h;
"""%n
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))

def EllipticCurve_from_cubic(F, P):
    r"""
    Construct an elliptic curve from a ternary cubic with a rational point.
    
    INPUT:

    - ``F`` -- a homogeneous cubic in three variables with rational
      coefficients (either as a polynomial ring element or as a
      string) defining a smooth plane cubic curve.

    - ``P`` -- a 3-tuple `(x,y,z)` defining a projective point on the
      curve `F=0`.

    OUTPUT:

    (elliptic curve) An elliptic curve (in minimal Weierstrass form)
    isomorphic to the curve `F=0`.

    .. note::

       USES MAGMA - This function will not work on computers that
       do not have magma installed.

    TO DO: implement this without using MAGMA.

    For a more general version, see the function
    ``EllipticCurve_from_plane_curve()``.
    
    EXAMPLES:

    First we find that the Fermat cubic is isomorphic to the curve
    with Cremona label 27a1::
    
        sage: E = EllipticCurve_from_cubic('x^3 + y^3 + z^3', [1,-1,0])  # optional - magma
        sage: E         # optional - magma
        Elliptic Curve defined by y^2 + y = x^3 - 7 over Rational Field
        sage: E.cremona_label()     # optional - magma
        '27a1'
    
    Next we find the minimal model and conductor of the Jacobian of the
    Selmer curve.
    
    ::
    
        sage: E = EllipticCurve_from_cubic('u^3 + v^3 + 60*w^3', [1,-1,0])   # optional - magma
        sage: E                # optional - magma
        Elliptic Curve defined by y^2  = x^3 - 24300 over Rational Field
        sage: E.conductor()    # optional - magma
        24300
    """
    from sage.interfaces.all import magma

    cmd = "P<%s,%s,%s> := ProjectivePlane(RationalField());" % SR(F).variables()
    magma.eval(cmd)
    cmd = "aInvariants(MinimalModel(EllipticCurve(Curve(Scheme(P, %s)),P!%s)));" % (F, P)
    s = magma.eval(cmd)
    return EllipticCurve(rings.RationalField(), eval(s))

def vecmat_ZZ(n=300, min=-9, max=9, system='sage', times=200):
    """
    Vector matrix multiplication over ZZ.

    Given an n x n  matrix A over ZZ with random entries
    between min and max, inclusive, and v the first row of A,
    compute the product v * A.

    INPUT:

    - ``n`` - matrix dimension (default: ``300``)
    - ``min`` - minimal value for entries of matrix (default: ``-9``)
    - ``max`` - maximal value for entries of matrix (default: ``9``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')
    - ``times`` - number of runs (default: ``200``)

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.vecmat_ZZ(300)  # long time
        sage: tm = b.vecmat_ZZ(300, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(ZZ, n, n, x=min, y=max+1)
        v = A.row(0)
        t = cputime()
        for z in range(times):
            w = v * A
        return cputime(t)/times
    elif system == 'magma':
        code = """
n := %s;
A := MatrixAlgebra(IntegerRing(), n)![Random(%s,%s) : i in [1..n^2]];
v := A[1];
t := Cputime();
for z in [1..%s] do
    K := v * A;
end for;
s := Cputime(t);
"""%(n,min,max,times)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))/times
    else:
        raise ValueError('unknown system "%s"'%system)

def matrix_add_ZZ(n=200, min=-9, max=9, system='sage', times=50):
    """
    Matrix addition over ZZ
    Given an n x n matrix A and B over ZZ with random entries between
    ``min`` and ``max``, inclusive, compute A + B ``times`` times.

    INPUT:

    - ``n`` - matrix dimension (default: ``200``)
    - ``min`` - minimal value for entries of matrix (default: ``-9``)
    - ``max`` - maximal value for entries of matrix (default: ``9``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')
    - ``times`` - number of experiments (default: ``50``)

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.matrix_add_ZZ(200)
        sage: tm = b.matrix_add_ZZ(200, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(ZZ, n, n, x=min, y=max+1)
        B = random_matrix(ZZ, n, n, x=min, y=max+1)
        t = cputime()
        for z in range(times):
            v = A + B
        return cputime(t)/times
    elif system == 'magma':
        code = """
n := %s;
min := %s;
max := %s;
A := MatrixAlgebra(IntegerRing(), n)![Random(min,max) : i in [1..n^2]];
B := MatrixAlgebra(IntegerRing(), n)![Random(min,max) : i in [1..n^2]];
t := Cputime();
for z in [1..%s] do
    K := A + B;
end for;
s := Cputime(t);
"""%(n,min,max,times)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))/times
    else:
        raise ValueError('unknown system "%s"'%system)

def MatrixVector_QQ(n=1000,h=100,system='sage',times=1):
    """
    Compute product of square ``n`` matrix by random vector with num and
    denom bounded by ``h`` the given number of ``times``.

    INPUT:

    - ``n`` - matrix dimension (default: ``300``)
    - ``h`` - numerator and denominator bound (default: ``bnd``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')
    - ``times`` - number of experiments (default: ``1``)

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.MatrixVector_QQ(500)
        sage: tm = b.MatrixVector_QQ(500, system='magma')  # optional - magma
    """
    if system=='sage':
        V=QQ**n
        v=V.random_element(h)
        M=random_matrix(QQ,n)
        t=cputime()
        for i in range(times):
            w=M*v
        return cputime(t)
    elif system == 'magma':
        code = """
            n:=%s;
            h:=%s;
            times:=%s;
            v:=VectorSpace(RationalField(),n)![Random(h)/(Random(h)+1) : i in [1..n]];
            M:=MatrixAlgebra(RationalField(),n)![Random(h)/(Random(h)+1) : i in [1..n^2]];
            t := Cputime();
            for z in [1..times] do
                W:=v*M;
            end for;
            s := Cputime(t);
        """%(n,h,times)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def matrix_multiply_QQ(n=100, bnd=2, system='sage', times=1):
    """
    Given an n x n matrix A over QQ with random entries
    whose numerators and denominators are bounded by bnd,
    compute A * (A+1).

    INPUT:

    - ``n`` - matrix dimension (default: ``300``)
    - ``bnd`` - numerator and denominator bound (default: ``bnd``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')
    - ``times`` - number of experiments (default: ``1``)

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.matrix_multiply_QQ(100)
        sage: tm = b.matrix_multiply_QQ(100, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(QQ, n, n, num_bound=bnd, den_bound=bnd)
        B = A + 1
        t = cputime()
        for z in range(times):
            v = A * B
        return cputime(t)/times
    elif system == 'magma':
        A = magma(random_matrix(QQ, n, n, num_bound=bnd, den_bound=bnd))
        code = """
n := %s;
A := %s;
B := A + 1;
t := Cputime();
for z in [1..%s] do
    K := A * B;
end for;
s := Cputime(t);
"""%(n, A.name(), times)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))/times
    else:
        raise ValueError('unknown system "%s"'%system)

def matrix_multiply_GF(n=100, p=16411, system='sage', times=3):
    """
    Given an n x n matrix A over GF(p) with random entries, compute
    A * (A+1).

    INPUT:

    - ``n`` - matrix dimension (default: 100)
    - ``p`` - prime number (default: ``16411``)
    - ``system`` - either 'magma' or 'sage' (default: 'sage')
    - ``times`` - number of experiments (default: ``3``)

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.matrix_multiply_GF(100, p=19)
        sage: tm = b.matrix_multiply_GF(100, p=19, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(GF(p), n)
        B = A + 1
        t = cputime()
        for n in range(times):
            v = A * B
        return cputime(t) / times
    elif system == 'magma':
        code = """
n := %s;
A := Random(MatrixAlgebra(GF(%s), n));
B := A + 1;
t := Cputime();
for z in [1..%s] do
    K := A * B;
end for;
s := Cputime(t);
"""%(n,p,times)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))/times
    else:
        raise ValueError('unknown system "%s"'%system)

def matrix_add_GF(n=1000, p=16411, system='sage',times=100):
    """
    Given two n x n matrix over GF(p) with random entries, add them.

    INPUT:

    - ``n`` - matrix dimension (default: 300)
    - ``p`` - prime number (default: ``16411``)
    - ``system`` - either 'magma' or 'sage' (default: 'sage')
    - ``times`` - number of experiments (default: ``100``)

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.matrix_add_GF(500, p=19)
        sage: tm = b.matrix_add_GF(500, p=19, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(GF(p), n, n)
        B = random_matrix(GF(p), n, n)
        t = cputime()
        for n in range(times):
            v = A + B
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := Random(MatrixAlgebra(GF(%s), n));
B := Random(MatrixAlgebra(GF(%s), n));
t := Cputime();
for z in [1..%s] do
    K := A + B;
end for;
s := Cputime(t);
"""%(n,p,p,times)
        if verbose: print(code)
        magma.eval(code)
        return magma.eval('s')
    else:
        raise ValueError('unknown system "%s"'%system)

def nullspace_RDF(n=300, min=0, max=10, system='sage'):
    """
    Nullspace over RDF:
    Given a n+1 x n  matrix over RDF with random entries
    between min and max, compute the nullspace.

    INPUT:

    - ``n`` - matrix dimension (default: ``300``)
    - ``min`` - minimal value for entries of matrix (default: ``0``)
    - ``max`` - maximal value for entries of matrix (default: `10``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.nullspace_RDF(100)  # long time
        sage: tm = b.nullspace_RDF(100, system='magma')  # optional - magma
    """
    if system == 'sage':
        from sage.rings.real_double import RDF
        A = random_matrix(ZZ, n+1, n, x=min, y=max+1).change_ring(RDF)
        t = cputime()
        v = A.kernel()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := RMatrixSpace(RealField(16), n+1,n)![Random(%s,%s) : i in [1..n*(n+1)]];
t := Cputime();
K := Kernel(A);
s := Cputime(t);
"""%(n,min,max)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def det_QQ(n=300, num_bound=10, den_bound=10, system='sage'):
    """
    Dense rational determinant over QQ.
    Given an n x n matrix A over QQ with random entries
    with numerator bound and denominator bound, compute det(A).

    INPUT:

    - ``n`` - matrix dimension (default: ``200``)
    - ``num_bound`` - numerator bound, inclusive (default: ``10``)
    - ``den_bound`` - denominator bound, inclusive (default: ``10``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.det_QQ(200)
        sage: ts = b.det_QQ(10, num_bound=100000, den_bound=10000)
        sage: tm = b.det_QQ(200, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(QQ, n, n, num_bound=num_bound, den_bound=den_bound)
        t = cputime()
        d = A.determinant()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := MatrixAlgebra(RationalField(), n)![Random(%s,%s)/Random(1,%s) : i in [1..n^2]];
t := Cputime();
d := Determinant(A);
s := Cputime(t);
"""%(n,-num_bound, num_bound, den_bound)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def rank2_ZZ(n=400, min=0, max=2**64, system='sage'):
    """
    Rank 2 over ZZ:
    Given a (n + 10) x n matrix over ZZ with random entries
    between min and max, compute the rank.

    INPUT:

    - ``n`` - matrix dimension (default: ``400``)
    - ``min`` - minimal value for entries of matrix (default: ``0``)
    - ``max`` - maximal value for entries of matrix (default: ``2**64``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.rank2_ZZ(300)
        sage: tm = b.rank2_ZZ(300, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(ZZ, n+10, n, x=min, y=max+1)
        t = cputime()
        v = A.rank()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := RMatrixSpace(IntegerRing(), n+10, n)![Random(%s,%s) : i in [1..n*(n+10)]];
t := Cputime();
K := Rank(A);
s := Cputime(t);
"""%(n,min,max)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def det_ZZ(n=200, min=1, max=100, system='sage'):
    """
    Dense integer determinant over ZZ.
    Given an n x n matrix A over ZZ with random entries
    between min and max, inclusive, compute det(A).

    INPUT:

    - ``n`` - matrix dimension (default: ``200``)
    - ``min`` - minimal value for entries of matrix (default: ``1``)
    - ``max`` - maximal value for entries of matrix (default: ``100``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.det_ZZ(200)
        sage: tm = b.det_ZZ(200, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(ZZ, n, n, x=min, y=max+1)
        t = cputime()
        d = A.determinant()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := MatrixAlgebra(IntegerRing(), n)![Random(%s,%s) : i in [1..n^2]];
t := Cputime();
d := Determinant(A);
s := Cputime(t);
"""%(n,min,max)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def smithform_ZZ(n=128, min=0, max=9, system='sage'):
    """
    Smith Form over ZZ:
    Given a n x n matrix over ZZ with random entries
    between min and max, compute the Smith normal form.

    INPUT:

    - ``n`` - matrix dimension (default: ``128``)
    - ``min`` - minimal value for entries of matrix (default: ``0``)
    - ``max`` - maximal value for entries of matrix (default: ``9``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.smithform_ZZ(100)
        sage: tm = b.smithform_ZZ(100, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(ZZ, n, n, x=min, y=max+1)
        t = cputime()
        v = A.elementary_divisors()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := MatrixAlgebra(IntegerRing(), n)![Random(%s,%s) : i in [1..n^2]];
t := Cputime();
K := ElementaryDivisors(A);
s := Cputime(t);
"""%(n,min,max)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def inverse_QQ(n=100, min=0, max=9, system='sage'):
    """
    Given a n x n matrix over QQ with random integer entries
    between min and max, compute the reduced row echelon form.

    INPUT:

    - ``n`` - matrix dimension (default: ``300``)
    - ``min`` - minimal value for entries of matrix (default: ``-9``)
    - ``max`` - maximal value for entries of matrix (default: ``9``)
    - ``system`` - either 'sage' or 'magma' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.inverse_QQ(100)
        sage: tm = b.inverse_QQ(100, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(ZZ, n, n, x=min, y=max+1).change_ring(QQ)
        t = cputime()
        v = ~A
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := MatrixAlgebra(RationalField(), n)![Random(%s,%s) : i in [1..n*n]];
t := Cputime();
K := A^(-1);
s := Cputime(t);
"""%(n,min,max)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def det_GF(n=400, p=16411 , system='sage'):
    """
    Dense determinant over GF(p).
    Given an n x n matrix A over GF with random entries compute
    det(A).

    INPUT:

    - ``n`` - matrix dimension (default: 300)
    - ``p`` - prime number (default: ``16411``)
    - ``system`` - either 'magma' or 'sage' (default: 'sage')

    EXAMPLES::

        sage: import sage.matrix.benchmark as b
        sage: ts = b.det_GF(1000)
        sage: tm = b.det_GF(1000, system='magma')  # optional - magma
    """
    if system == 'sage':
        A = random_matrix(GF(p), n, n)
        t = cputime()
        d = A.determinant()
        return cputime(t)
    elif system == 'magma':
        code = """
n := %s;
A := Random(MatrixAlgebra(GF(%s), n));
t := Cputime();
d := Determinant(A);
s := Cputime(t);
"""%(n,p)
        if verbose: print code
        magma.eval(code)
        return float(magma.eval('s'))
    else:
        raise ValueError('unknown system "%s"'%system)

def hilbert_class_polynomial(D, algorithm=None):
    r"""
    Returns the Hilbert class polynomial for discriminant `D`.

    INPUT:

    - ``D`` (int) -- a negative integer congruent to 0 or 1 modulo 4.

    - ``algorithm`` (string, default None).

    OUTPUT:

    (integer polynomial) The Hilbert class polynomial for the
    discriminant `D`.

    ALGORITHM:

    - If ``algorithm`` = "arb" (default): Use Arb's implementation which uses complex interval arithmetic.

    - If ``algorithm`` = "sage": Use complex approximations to the roots.

    - If ``algorithm`` = "magma": Call the appropriate Magma function (if available).

    AUTHORS:

    - Sage implementation originally by Eduardo Ocampo Alvarez and
      AndreyTimofeev

    - Sage implementation corrected by John Cremona (using corrected precision bounds from Andreas Enge)

    - Magma implementation by David Kohel

    EXAMPLES::

        sage: hilbert_class_polynomial(-4)
        x - 1728
        sage: hilbert_class_polynomial(-7)
        x + 3375
        sage: hilbert_class_polynomial(-23)
        x^3 + 3491750*x^2 - 5151296875*x + 12771880859375
        sage: hilbert_class_polynomial(-37*4)
        x^2 - 39660183801072000*x - 7898242515936467904000000
        sage: hilbert_class_polynomial(-37*4, algorithm="magma") # optional - magma
        x^2 - 39660183801072000*x - 7898242515936467904000000
        sage: hilbert_class_polynomial(-163)
        x + 262537412640768000
        sage: hilbert_class_polynomial(-163, algorithm="sage")
        x + 262537412640768000
        sage: hilbert_class_polynomial(-163, algorithm="magma") # optional - magma
        x + 262537412640768000

    TESTS::

        sage: all([hilbert_class_polynomial(d, algorithm="arb") == \
        ....:      hilbert_class_polynomial(d, algorithm="sage") \
        ....:        for d in range(-1,-100,-1) if d%4 in [0,1]])
        True

    """
    if algorithm is None:
        algorithm = "arb"

    D = Integer(D)
    if D >= 0:
        raise ValueError("D (=%s) must be negative"%D)
    if not (D%4 in [0,1]):
         raise ValueError("D (=%s) must be a discriminant"%D)

    if algorithm == "arb":
        import sage.libs.arb.arith
        return sage.libs.arb.arith.hilbert_class_polynomial(D)

    if algorithm == "magma":
        magma.eval("R<x> := PolynomialRing(IntegerRing())")
        f = str(magma.eval("HilbertClassPolynomial(%s)"%D))
        return IntegerRing()['x'](f)

    if algorithm != "sage":
        raise ValueError("%s is not a valid algorithm"%algorithm)

    from sage.quadratic_forms.binary_qf import BinaryQF_reduced_representatives
    from sage.rings.all import RR, ZZ, ComplexField
    from sage.functions.all import elliptic_j

    # get all primitive reduced quadratic forms, (necessary to exclude
    # imprimitive forms when D is not a fundamental discriminant):

    rqf = BinaryQF_reduced_representatives(D, primitive_only=True)

    # compute needed precision
    #
    # NB: [http://arxiv.org/abs/0802.0979v1], quoting Enge (2006), is
    # incorrect.  Enge writes (2009-04-20 email to John Cremona) "The
    # source is my paper on class polynomials
    # [http://hal.inria.fr/inria-00001040] It was pointed out to me by
    # the referee after ANTS that the constant given there was
    # wrong. The final version contains a corrected constant on p.7
    # which is consistent with your example. It says:

    # "The logarithm of the absolute value of the coefficient in front
#.........这里部分代码省略.........

def EllipticCurve_from_plane_curve(C, P):
    r"""
    Construct an elliptic curve from a smooth plane cubic with a rational point.
    
    INPUT:

    - ``C`` -- a plane curve of genus one.

    - ``P`` -- a 3-tuple `(x,y,z)` defining a projective point on the
      curve ``C``.

    OUTPUT:

    (elliptic curve) An elliptic curve (in minimal Weierstrass form)
    isomorphic to ``C``.


    .. note::

       USES MAGMA - This function will not work on computers that
       do not have magma installed.

    TO DO: implement this without using MAGMA.

    EXAMPLES:

    First we check that the Fermat cubic is isomorphic to the curve 
    with Cremona label '27a1'::

        sage: x,y,z=PolynomialRing(QQ,3,'xyz').gens() # optional - magma 
        sage: C=Curve(x^3+y^3+z^3) # optional - magma
        sage: P=C(1,-1,0) # optional - magma
        sage: E=EllipticCurve_from_plane_curve(C,P) # optional - magma
        sage: E # optional - magma
        Elliptic Curve defined by y^2 + y = x^3 - 7 over Rational Field
        sage: E.label() # optional - magma
        '27a1'

    Now we try a quartic example::
    
        sage: u,v,w=PolynomialRing(QQ,3,'uvw').gens() # optional - magma 
        sage: C=Curve(u^4+u^2*v^2-w^4) # optional - magma
        sage: P=C(1,0,1) # optional - magma
        sage: E=EllipticCurve_from_plane_curve(C,P) # optional - magma
        sage: E # optional - magma
        Elliptic Curve defined by y^2  = x^3 + 4*x over Rational Field
        sage: E.label() # optional - magma
        '32a1'

        """
    from sage.interfaces.all import magma
    if C.genus()!=1:
        raise TypeError, "The curve C must have genus 1"
    elif P.parent()!=C.point_set(C.base_ring()):
        raise TypeError, "The point P must be on the curve C"
    dp=C.defining_polynomial()
    x,y,z = dp.parent().variable_names()
    cmd = "PR<%s,%s,%s>:=ProjectivePlane(RationalField());"%(x,y,z)
    magma.eval(cmd)
    cmd = 'CC:=Curve(PR, %s);'%(dp)
    magma.eval(cmd)
    cmd='aInvariants(MinimalModel(EllipticCurve(CC,CC!%s)));'%([P[0],P[1],P[2]])
    s=magma.eval(cmd)
    return EllipticCurve(rings.RationalField(), eval(s))