One Hat Cyber Team

View File Name : RationalPoints2.m2

-------------------------------------------------------------------------------
-- licensed under GPL v2 or any later version
-------------------------------------------------------------------------------
newPackage(
    "RationalPoints2",
    Version => "0.5",
    Date => "Mar 18, 2021",
    Authors => { {Name => "Jieao Song", Email => "jieao.song@imj-prg.fr"} },
    Headline => "find the rational points on a variety",
    Keywords => {"Commutative Algebra"},
    DebuggingMode => false,
    PackageImports => {"Elimination", "Varieties"}
)
-------------------------------------------------------------------------------
-- This program takes in an ideal generated by a list of polynomials, and finds
-- all of the common zeros of the polynomials. Over a finite field the ideal
-- can be arbitrary; over a number field either the ideal should be
-- 0-dimensional, or a height bound should be specified so that the field
-- elements can be enumerated.
-------------------------------------------------------------------------------
-- The program naively tests points coordinate by coordinate. There are however
-- several tricks used to improve the performance.
--
-- * At each coordinate we use elimination to compute the possible values
-- instead of enumerating over the entire field. The elimination is done by
-- computing a Groebner basis for the Eliminate monomial order.
--
-- * If we get a prime ideal of degree 1, we solve a linear system of equations
-- to quickly get all the points contained in the corresponding linear
-- subspace.
--
-- * If the ideal is homogeneous, we will enumerate the rational points in the
-- projective space first, then reconstruct the points in the affine space;
-- this can also be done recursively.
-------------------------------------------------------------------------------
-- Some utility functions are provided for handling field extensions.
-------------------------------------------------------------------------------
export{
    "integers", "hermiteNormalForm", "ProjectivePoint", "globalHeight",
    "baseChange", "extField", "splittingField", "charpoly", "minpoly", "zeros",
    "rationalPoints", "Amount", "Split", "Bound", "KeepAll", "setPariSize"
}
-------------------------------------------------------------------------------
-- ELS: a list of elements in the finite field k
-- When over a number field/a very large finite field, we set ELS = {}
local ELS;
-- the options
local AMOUNT;
local PROJECTIVE;
local VERBOSE;
local SPLIT;
local BOUND;
local KEEPALL;
-- NONSPLITPOLYS: a list of non-splitting polynomials, used in verbose mode
-- FIELD: the computed splitting field
local NONSPLITPOLYS;
local FIELD;
-- the variable used for field extension, default to `a` like for GF
local VARIABLE;
VARIABLE = "a";
-------------------------------------------------------------------------------
-- Takes a number n and a list
-- Produces the direct product of the list with itself n times
-- Same as lst**lst**lst**... but with flattened elements
--
pow = (n, lst) -> (
    if n < 0 then error "power of a list is not defined for n negative";
    result := {{}};
    for i in (0..<n) do result = toList \ flatten \ (result ** lst);
    return result;
);
-------------------------------------------------------------------------------
-- Utility function that shuffles two list u and v according to a list ind
-- of indices in 0 (for u) and 1 (for v)
--
shuffle = (ind, u, v) -> (
    (i,j) := (-1,-1);
    return for k in (0..<#ind) list if ind_k==0 then (i=i+1; u_i) else (j=j+1; v_j);
);
-------------------------------------------------------------------------------
-- Takes in a finite field and returns the list of its elements
--
listFiniteFieldElements = k -> (
    p := char k;
    if isFinitePrimeField k then return apply(p, i->i*1_k);
    -- else k is either a GF or constructed using `toField`
    d := degree baseRing k;
    if p^d < 2^29 then return (pow_d apply(p, i->i*1_k)) / (c -> sum apply(d, i -> c_i * k_0^i))
    else error "the field is too big";
);
-------------------------------------------------------------------------------
-- Class of a projective point, mainly used for pretty printing
--
ProjectivePoint = new Type of List;
net ProjectivePoint := (pt) -> net "("| horizontalJoin between(net" : ", net \ pt)|net ")";
texMath ProjectivePoint := (pt) -> "\\left("|concatenate between(":", texMath \ pt)|"\\right)";
undocumented{(net, ProjectivePoint),(texMath, ProjectivePoint)};
-------------------------------------------------------------------------------
ProjectiveVariety VisibleList := (X, x) -> (
    R := ring X;
    k := coefficientRing R;
    if not numgens R == #x then error "wrong number of coordinates";
    if all(x, zero) then error "homogeneous coordinates cannot be all zero";
    if not zero (map(k,ring ideal R,toList x)) ideal R then error "the point does not lie on the variety";
    return new ProjectivePoint from apply(homogCoord x, xi->xi*1_k);
);
ProjectivePoint == ProjectivePoint := (p1, p2) -> (
    if not #p1 == #p2 then error "wrong number of coordinates";
    p1 = homogCoord toList p1;
    p2 = homogCoord toList p2;
    return p1 == p2;
);
ring ProjectivePoint := pt -> (ring pt_0);
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-- This part handles various aspects of number fields
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-- Cache the following values for a number field for quick access
--
primitiveElement = (cacheValue symbol primitiveElement) (F -> (baseRing F)_0);
primitiveElementValues = (cacheValue symbol primitiveElementValues) (F -> roots minpoly F);
extensionDegree = (cacheValue symbol extensionDegree) (F -> degree baseRing F);
toBaseRing = (cacheValue symbol toBaseRing) (F -> (map(baseRing F, F, gens baseRing F)));
embed := (p,x) -> (map(CC, ring x, {p})) x;
embed = memoize embed;
-------------------------------------------------------------------------------
-- Get the matrix of an integral principal ideal in Hermite normal form
--
getIntegralMatrix := x -> (
    F := ring x;
    n := extensionDegree F;
    a := primitiveElement F;
    x = (toBaseRing F) x;
    return gens gb lift(sum apply(n, i->(multiTable F)_i*coefficient(a^i, x)), ZZ);
);
getIntegralMatrix = memoize getIntegralMatrix;
-------------------------------------------------------------------------------
-- multiplicative height function
--
ht := x -> (
    k := ring x_0;
    -- Over QQ is easy
    if k === QQ then return (lcm \\ denominator \ x) / (gcd \\ abs@@numerator \ x) * max(abs \ x);
    d := extensionDegree k;
    -- First get rid of the denominators
    denom := lcm flatten apply(x, xi->denominator \ last \ listForm((toBaseRing k) xi));
    x = apply(x, xi->xi*denom);
    -- Compute the norm of the gcd (which is an integral ideal instead of an integer!)
    finitePlaces := 1 / det gens gb fold((a,b)->(a|b), getIntegralMatrix \ x);
    place := (p, x) -> max(x / (xi -> abs embed_p xi)); -- this is less accurate than `baseChange` but a lot faster
    infinitePlaces := product for p in primitiveElementValues k list place_p x;
    return (finitePlaces * infinitePlaces)^(1/d);
);
ht = memoize ht;
-------------------------------------------------------------------------------
-- exported interface
--
globalHeight = method(Options=>{});
globalHeight List := opts -> x -> (
    if not (same(ring \ x) and (k := ring x_0; isField k and char k == 0))
    then error "height is only defined for a list of elements in a number field";
    return ht x;
);
-------------------------------------------------------------------------------
-- Compute a basis of the ring of integers
-- Zassenhaus' Round 2, Algorithm 6.1.8 in Cohen, A course in computational
-- algebraic number theory
--
round2Integers := F -> (
    n := extensionDegree F;
    T := minpoly F;
    x := gens ring T;
    D := lift(((-1)^(n*(n-1)/2) * resultant(T, diff(x_0, T), x_0)), ZZ);
    recover := map(F, baseRing F, gens F);
    F0 := F;
    F = baseRing F;
    a := F_0;
    result := apply(n, i->a^i);
    getCoeffs := x -> apply(n, i->coefficient(a^i, x));
    getAbsMatrix := x -> transpose matrix for w in x list getCoeffs w;
    getMatrix := x -> inverse getAbsMatrix result * getAbsMatrix x;
    for fct in factor abs D do if (last fct // 2 > 0) then (
        p := first fct;
        modp := (ZZ/p)(monoid[x]);
        T' := sub(T, modp);
        g := product for f in factor T' list first f;
        h := T'//g;
        f := sub((sub(g, ring T)*sub(h, ring T) - T) / p, modp);
        Z := gcd(f, g, h);
        m := (degree Z)_0;
        if m == 0 then continue;
        U := (map(F, modp, {a})) (T'//Z);
        M := gens gb (lift(getMatrix apply(result_{0..m-1}, w->w*U), ZZ) | p*id_(ZZ^n));
        result = first entries (1/p*matrix{result}*sub(M, F));
        if (last fct // 2) >= m+1 then while true do (
            q := p;
            while (q<n) do q = q*p;
            A := sub(getMatrix apply(result, w->w^q), ZZ/p);
            alpha := first entries (matrix{result} * sub(gens kernel A, F));
            if alpha == {} then (alpha = for w in result list p*w; M = inverse getAbsMatrix alpha;)
            else (
                M = getAbsMatrix alpha;
                r := rank M;
                for w in result do (
                    B := getAbsMatrix{p*w};
                    if rank (M|B) > r then (r=r+1; alpha=alpha|{p*w}; M=M|B;);
                    if r == n then break;
                );
                M = inverse M;
            );
            C := sub(transpose matrix apply(result, w->flatten entries (M * getAbsMatrix apply(alpha, a->w*a))), ZZ/p);
            kerC := first entries (matrix{result} * lift(gens kernel C, ZZ));
            if #kerC == 0 then break;
            M = gens gb (lift(getMatrix kerC, ZZ) | p*id_(ZZ^n));
            result = first entries (1/p*matrix{result}*sub(M, F));
        );
    );
    disc := lift(D * (product apply(n, i->coefficient(a^i, result_i)))^2, ZZ);
    return (recover \ result, disc);
);
round2Integers = memoize round2Integers;
-------------------------------------------------------------------------------
-- Multiplication table in the basis O_F
--
multiTable = (cacheValue symbol multiTable) (F -> (
    n := extensionDegree F;
    a := primitiveElement F;
    OF := (toBaseRing F) \ first round2Integers F;
    getCoeffs := s -> apply(n, i->coefficient(a^i, s));
    getAbsMatrix := s -> transpose matrix for w in s list getCoeffs w;
    M := inverse getAbsMatrix OF;
    return apply(n, i->M*getAbsMatrix apply(OF, x->x*a^i));
));
-------------------------------------------------------------------------------
-- Exported interface
--
integers = method();
integers Ring := F -> (
    if not (isField F and char F == 0) then error "expect a number field";
    return if F === QQ then {1_QQ} else first round2Integers F;
);
discriminant Ring := ZZ => o -> F -> (
    if not (isField F and char F == 0) then error "expect a number field";
    return if F === QQ then 1 else last round2Integers F;
);
-------------------------------------------------------------------------------
-- Compute the Hermite normal form of a fractional ideal
--
hermiteNormalForm = method();
hermiteNormalForm RingElement := x -> hermiteNormalForm {x};
hermiteNormalForm List := x -> (
    F := ring x_0;
    n := extensionDegree F;
    a := primitiveElement F;
    -- First get rid of the denominators
    denom := lcm flatten apply(x, xi->denominator \ last \ listForm((toBaseRing F) xi));
    x = apply(x, xi->xi*denom);
    M := gens gb fold((a,b)->(a|b), getIntegralMatrix \ x);
    return (1/denom)*sub(M, QQ);
);
-------------------------------------------------------------------------------
-- List elements in a number field with bounded height
-- Depends on Sage for number fields other than QQ
--
listNumberFieldElements = (k, B) -> (
    if k === QQ then return unique flatten for q in (1..floor B) list for p in (-floor B..floor B) list (p/q)*1_QQ;
    try return bddHeight(k, B) else notImplemented();
);
-------------------------------------------------------------------------------
-- The interface for Sage
-- Sage needs around 1 second to start up every time...
--
bddHeight = (k, B) -> (
    sage := findProgram("sage", "sage -v");
    p := minpoly k;
    R := ring p;
    d := (degree p)_0;
    R' := k(monoid[getSymbol "x"]);
    p = (map(R', R, gens R')) p;
    INPUT := temporaryFileName()|".py";
    F := openOut INPUT;
    F << "from sage.all import NumberField, PolynomialRing, QQ, RR\n"
      << "x = PolynomialRing(QQ, 'x').gen()\n"
      << "K = NumberField("|replace("\\^", "**", toString p)|", 'a')\n"
      << "pts = list(K.elements_of_bounded_height(bound=RR("|toString (B^d)|")))\n"
      << "print(','.join(str(pt.list()) for pt in pts))" << close;
    runResult := runProgram(sage, INPUT);
    assert zero runResult#"return value";
    coeffs := value("{"|runResult#"output"|"}");
    removeFile INPUT;
    if VERBOSE then print("-- Sage has finished and returned "|toString(#coeffs)|" elements");
    return for c in coeffs list sum apply(d, i->(c_i)*(k_0)^i);
);
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-- This part handles field extensions
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-- Takes in the `listForm` of a univariate polynomial
-- E.g. {({2}, 1), ({0}, 1)} -> 1+x^2
-- Returns a list of zeros and a list of non-splitting factors
-- E.g. 2-3*x+x^2 -> ({1, 2}, {})
--      1+x^2     -> ({}, {x^2+1})
--
getZeros := (coeffs) -> (
    k := ring coeffs_0_1;
    R := k(monoid[getSymbol VARIABLE]);
    getValue := c -> -coefficient(1_R, c) / coefficient(R_0, c);
    p := sum apply(coeffs, t -> t_1 * R_0^(t_0_0));
    factors := value \ first \ toList factor p;
    linearFactors := select(factors, c->(first degree c == 1));
    otherFactors := select(factors, c->(first degree c > 1));
    return (getValue \ linearFactors, otherFactors);
);
getZeros = memoize getZeros;
-------------------------------------------------------------------------------
-- The exported interface
--
zeros = method(Options => {});
zeros(RingElement) := opts -> p -> first getZeros listForm p;
zeros(Ring, RingElement) := opts -> (F, p) -> zeros baseChange_F p;
-------------------------------------------------------------------------------
-- Utility function that tests if F is prime
--
isPrimeField = F -> (F === QQ or isFinitePrimeField F or (instance(F, GaloisField) and F.degree == 1));
-------------------------------------------------------------------------------
-- Given a field F and a polynomial p with coefficients in the prime field k of
-- F, check if x is indeed a zero, or if x == FindOne, try find a zero
-- Returns the zero or an error will be raised
--
checkZero := (F, p, x) -> (
    k := coefficientRing ring p; -- k is the prime field
    assert isField F;
    if x === FindOne then (
        coeffs := apply(listForm p, t->(t_0, (map(F,k,{})) t_1));
        lifts := first getZeros coeffs;
        assert(#lifts > 0);
        return lifts_0;
    );
    assert(F === ring x);
    assert zero sum apply(listForm p, t->t_1*x^(t_0_0));
    return x;
);
checkZero = memoize checkZero;
-------------------------------------------------------------------------------
-- el is some element to perform the base change on
-- It can be an element in a base field F0 / a polynomial with coefficients in
-- F0 / an ideal with coefficients in F0
--
baseChange = method(Options => {PrimitiveElement=>FindOne});
baseChange(Ring, Number) := baseChange(Ring, RingElement) := baseChange(Ring, Ideal) := opts -> (F, el) -> (
    if not isField F then error "the coefficient ring is not a field";
    R := ring el;
    prim := opts.PrimitiveElement;
    if isField R then ( -- el is a field element
        if isPrimeField R then return (map(F,R,{})) el;
        try (prim = checkZero(F, minpoly R, prim)) else error "the primitive element is not valid";
        return (map(F,R,{prim})) el;
    ) else if isPolynomialRing R then ( -- el is a polynomial or an ideal
        F0 := coefficientRing R;
        if not isField F0 then error "the coefficient ring is not a field";
        R' := F(R.monoid); -- this way the monomial order will be the same for R'
        if isPrimeField F0 then return (map(R', R, gens R')) el;
        -- else do a base change for the coefficients
        try (prim = checkZero(F, minpoly F0, prim)) else error "the primitive element is not valid";
        return (map(R', R, gens R'|{prim})) el;
    ) else error "cannot perform the base change";
);
-------------------------------------------------------------------------------
-- Base change to a prime characteristic
--
baseChange(Number, Number) := baseChange(Number, RingElement) := baseChange(Number, Ideal) := opts -> (p, el) -> (
    R := ring el;
    if not isPrime p then error "expect a prime number";
    if char R != 0 then error "already in positive characteristic";
    F := if isField R then R else if isPolynomialRing R then coefficientRing R else error "cannot perform the base change";
    F = extField(ZZ/p, {minpoly F});
    return baseChange(F, el);
);
-------------------------------------------------------------------------------
-- Base change to CC
-- First cache the complex values of the primitive element
--
baseChangeCC := (x, p) -> (
    xAprox := embed_p x;
    return first sort(roots minpoly x, v -> abs(v-xAprox));
);
baseChangeCC = memoize baseChangeCC;
-------------------------------------------------------------------------------
baseChange(InexactFieldFamily, Number) := baseChange(InexactFieldFamily, RingElement) := opts -> (C, x) -> (
    assert(C === CC); -- `roots` only returns values in CC
    F := ring x;
    if not (isField F and char F == 0) then error "not a number field";
    if F === QQ then return toCC x;
    prim := opts.PrimitiveElement;
    if prim === FindOne then prim = first primitiveElementValues F
    else prim = first sort(primitiveElementValues F, v -> abs(v-prim));
    return baseChangeCC(x, prim);
);
-------------------------------------------------------------------------------
-- Rudimentary simplification of the minimal polynomial
-- Makes the polynomial monic and divides out multiple factors
-- E.g. 4*x^2+9 -> x^2+1
-- Will first try to use `polredbest` if PARI/GP is available
--
reducePoly = p -> (
    R := ring p;
    assert(isPolynomialRing R and #gens R==1 and coefficientRing R===QQ);
    try return polredbest p; -- try use `polredbest` first
    d := (degree p)_0;
    denom := lcm \\ denominator \ last \ listForm p;
    c := leadCoefficient p;
    p = denom^d * c^(d-1) * (map(R, R, {R_0/c/denom})) p;
    getPower := (d, n) -> (if n == 0 then 0
        else value apply(factor abs numerator n, f->((first f)^(last f // d))));
    c = gcd apply(drop(listForm p, 1), t -> getPower(d-t_0_0, t_1));
    p = c^(-d) * (map(R, R, {R_0*c})) p;
    return p;
);
-------------------------------------------------------------------------------
-- Allows PARISIZE to be set
--
PARISIZE = 8000000;
setPariSize = n -> (PARISIZE = n);
-------------------------------------------------------------------------------
-- The interface for using `polredbest` in GP
--
polredbest = p -> (
    gp := findProgram("gp", "gp --version");
    R := ring p;
    k := coefficientRing R;
    d := (degree p)_0;
    UID := temporaryFileName();
    INPUT := UID|".gp";
    OUTPUT := UID|"-output";
    F := openOut INPUT;
    F << "allocatemem("|toString PARISIZE|")\n"
      << "f=polredbest("|toString p|")\n"
      << "for(d=0,poldegree(f),write1(\""|OUTPUT|"\",polcoeff(f,d),\",\"))\n"
      << "quit()" << close;
    assert zero (runProgram(gp, "-q <"|INPUT))#"return value";
    coeffs := value("{"|get OUTPUT|"}");
    removeFile \ {INPUT, OUTPUT};
    return sum apply(d+1, i -> coeffs_i*R_0^i);
);
-------------------------------------------------------------------------------
-- Perform a field extension by adding a zero of an irreducible polynomial
-- The complicated part is when the base field is already an extension of type
-- k[z]/minpoly(z). For finite fields only the degree is computed; for infinite
-- fields we try to find the new primitive element using elements of form
-- x+c*y.
--
getField := (F, p) -> (
    if not (isPolynomialRing ring p and #gens ring p == 1) then error "expect a univariate polynomial";
    ch := char F;
    k := if ch == 0 then QQ else ZZ/ch; -- k is the prime field
    if not member(char ring p, {0, ch}) then error "the field has different characteristic than the polynomial";
    -- decompose p in F
    coeffs := apply(listForm p, t->(t_0, baseChange_F t_1));
    factors := last getZeros coeffs;
    if #factors == 0 then return F; -- p splits in F already
    p = factors_0; -- otherwise we will take one irreducible factor
    d := (degree p)_0;
    R := F(monoid[getSymbol VARIABLE]); -- F[x]
    x := gens R;
    p = (map(R, ring p, x)) p;
    local F';
    if isPrimeField F then F' = (if ch == 0 then toField(R / reducePoly p)
        else GF(ch^d, Variable=>getSymbol VARIABLE))
    else ( -- case where F itself is already some extension k[y]/minpoly(y)
        q := minpoly F;
        e := (degree q)_0;
        if ch > 0 then F' = GF(ch^(d*e), Variable=>getSymbol VARIABLE)
        else (
            -- The extension is supposed to have degree d*e over the prime field
            -- We will try to find a primitive element of form x+c*y
            F' = toField first flattenRing(R/ideal p, CoefficientRing=>k);
            kx := k(monoid[x]); -- k[x]
            c := 1;
            while true do (
                p = (kernel map(F', kx, {F'_0+c*F'_1}))_0;
                if (degree p)_0 == d*e then break
                else c = c + 1;
            );
            F' = toField(kx / reducePoly p);
        );
    );
    return F';
);
getField = memoize getField;
-------------------------------------------------------------------------------
-- When the polynomial is explicitly typed in by the user using `extField`, set
-- the value of the symbol to its image in the extension field, like `toField`
--
assignZero = (F, p) -> (
    var := toString baseName (ring p)_0;
    -- only set the value when it is visible to the user
    if ring value var === ring p then getSymbol var <- first zeros baseChange_F p;
);
-------------------------------------------------------------------------------
-- Interface for doing a field extension
-- Use `Variable` to specify the variable used for the primitive element,
-- default to `a` like for GF
-- Assign the new value in the extension field to the variable of p
--
extField = method(Options => {Variable=>"a"});
extField(Ring, RingElement) := opts -> (F, p) -> (
    VARIABLE = if opts.Variable === "a" then "a" else toString baseName opts.Variable;
    F = getField(F, p);
    assignZero(F, p); -- to mimic the behavior of `toField`, discutable
    return F;
);
-------------------------------------------------------------------------------
extField(RingElement) := opts -> p -> extField(coefficientRing ring p, p, opts);
-------------------------------------------------------------------------------
extField(Ring, List) := opts -> (F, ps) -> (
    VARIABLE = if opts.Variable === "a" then "a" else toString baseName opts.Variable;
    for p in ps do F = getField(F, p);
    return F;
);
-------------------------------------------------------------------------------
extField(List) := opts -> ps -> (
    F := coefficientRing ring first ps;
    return extField(F, ps, opts);
);
-------------------------------------------------------------------------------
-- Find the splitting field of a polynomial by repeatedly applying `getField`
--
splittingField = method(Options => {Variable=>"a"});
splittingField(Ring, RingElement) := opts -> (F, p) -> (
    VARIABLE = if opts.Variable === "a" then "a" else toString baseName opts.Variable;
    if not (isPolynomialRing ring p and #gens ring p == 1) then error "expect a univariate polynomial";
    ch := char ring p;
    k := if ch == 0 then QQ else ZZ/ch; -- k is the prime field
    if char F != ch then error "the field has different characteristic than the polynomial";
    -- decompose p in F
    coeffs := apply(listForm p, t->(t_0, baseChange_F t_1));
    factors := last getZeros coeffs;
    if #factors == 0 then return F; -- p splits in F already
    for q in factors do ( -- for each factor q, apply repeatedly `getField`
        d := if isPrimeField F then 1 else extensionDegree F;
        while true do (
            newF := getField(F, q);
            if newF === F or (extensionDegree newF // d == ((degree q)_0)!) then (F = newF; break;);
            F = newF;
        );
    );
    return F;
);
-------------------------------------------------------------------------------
splittingField(RingElement) := opts -> p -> splittingField(coefficientRing ring p, p, opts);
-------------------------------------------------------------------------------
-- characteristic and minimal polynomial over prime field
--
charpoly = method(Options=>{Variable=>"x"});
charpoly(Number) := charpoly(RingElement) := opts -> (x) -> (
    F := ring x;
    if not isField F then error "the coefficient ring is not a field";
    d := if isPrimeField F then 1 else degree baseRing F;
    p := minpoly(x, opts);
    p ^ (d//(degree p)_0)
);
-------------------------------------------------------------------------------
minpoly = method(Options=>{Variable=>"x"});
minpoly(Number) := minpoly(RingElement) := opts -> (x) -> (
    VARIABLE = if opts.Variable === "x" then "x" else toString baseName opts.Variable;
    F := ring x;
    if not isField F then error "the coefficient ring is not a field";
    if isPrimeField F then (
        S := F(monoid[getSymbol VARIABLE]);
        S_0 - x
    ) else (
        k := coefficientRing baseRing F;
        p := (kernel map(F, k(monoid[getSymbol VARIABLE]), {x}))_0;
        p // leadCoefficient p
    )
);
-------------------------------------------------------------------------------
minpoly(Ring) := opts -> (F) -> (
    if not isField F then error "expect a field";
    if isPrimeField F then return minpoly(1_F, opts)
    else return (ideal baseRing F)_0;
);
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-- Finally, codes for finding rational points
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
-- Takes in a linear ideal and returns the list of solutions
--
linearSolve := (I) -> (
    R := ring I;
    k := coefficientRing R;
    x := gens R;
    n := #x;
    system := first entries gens gb I;
    if system == {} then return pow_n ELS;
    -- solve a system of linear equations Ax=b
    A := matrix apply(system, f->apply(x, xi->coefficient(xi, f)));
    b := matrix apply(system, f->{-coefficient(1_R, f)});
    v := transpose solve(A, b);
    if #system == n then return {first entries v};
    M := transpose gens kernel A;
    getSolution := x -> first entries (v + matrix{toList x} * M);
    return getSolution \ pow_(n-#system) ELS;
);
linearSolve = memoize linearSolve;
-------------------------------------------------------------------------------
-- Utility function that normalizes a list of homogeneous coordinates so that
-- the first non-zero coordinate is 1
--
homogCoord = p -> (a := p_(position(p, x->(x!=0))); (x->x/a) \ p);
-------------------------------------------------------------------------------
-- The main function that carries out the enumeration of points
-- Takes in an ideal and enumerate the points coordinate by coordinate
-- Returns the list of points or the number of points if `Amount` is set to
-- true
--
findPoints = I -> (
    result := if AMOUNT then 0 else {};
    if I == 1 then return result; -- trivial case
    R := ring I;
    k := coefficientRing R;
    x := gens R;
    n := #x;
    if n == 0 then return if AMOUNT then 1 else {{}};
    -- Enforce the `Eliminate` monomial order
    if not member(Weights=>toList(n-1:1), (options R).MonomialOrder) then (
        R = newRing(R, MonomialOrder=>Eliminate(n-1));
        x = gens R;
        I = (map(R, ring I, x)) I;
    );
    I = ideal gens gb I;
    if BOUND === null then ( -- these tricks do not get along well with height
        -- use `linearSolve` if I is linear
        if all(I_*, c -> degree c == {1}) then (
            return if AMOUNT then (#ELS)^(dim I) else linearSolve I;
        );
        -- use `findProjPoints` if I is homogeneous
        if isHomogeneous I then (
            result = findProjPoints I;
            return if AMOUNT then 1+(#ELS-1)*result
            else {toList(n:0_k)} | join flatten((v->(for a in ELS_{1..<#ELS} list a*v)) \ result);
        );
    ) else if I == 0 then return pow_n ELS;
    possibleValues := ELS;
    -- In the `Eliminate` monomial order with weight (1,...,1,0), if the first
    -- polynomial in the Groebner basis contains x_(n-1) as its only variable,
    -- then we can solve it to obtain all the possible values for x_(n-1)
    elim := I_0;
    if support elim == {x_(n-1)} then (
        coeffs := apply(listForm elim, t->(t_0_{n-1}, t_1));
        possibleValues = first getZeros coeffs;
        if BOUND =!= null and not KEEPALL then possibleValues = select(possibleValues, q -> ht {q,1_k} <= BOUND);
        nonSplit := last getZeros coeffs;
        -- handles the verbose outputs
        if #nonSplit > 0 and not SPLIT and VERBOSE and BOUND === null then (
            if #NONSPLITPOLYS == 0 then print "-- the following polynomials do not split";
            kx := k(monoid[getSymbol "x"]);
            nonSplit = (p -> toExternalString ((map(kx, ring p, gens kx)) p)) \ nonSplit;
            for fct in nonSplit do (
                if not member(fct, NONSPLITPOLYS) then (
                    print ("   -- " | fct);
                    NONSPLITPOLYS = NONSPLITPOLYS | {fct};
                );
            );
        );
        if #nonSplit > 0 and SPLIT then (
            k' := extField(k, {first nonSplit});
            I = baseChange_k' I;
            if #zeros_FIELD first nonSplit == 0 then (
                FIELD = extField(FIELD, {first nonSplit});
                if VERBOSE then print(
                    if char k > 0 then "-- base change to " | toString describe FIELD
                    else "-- base change to the field " | toExternalString baseRing FIELD);
            );
            return findPoints I;
        );
        -- replace the eliminant by 0 so we don't need to reevaluate it
        I = ideal ({0_R} | drop(I_*, 1));
    );
    R' := k(monoid[x_{0..<n-1},MonomialOrder=>Eliminate(n-2)]);
    x = gens R';
    for v in possibleValues do (
        I' := (map(R', R, x|{v})) I;
        if AMOUNT then result = result + findPoints I'
        else result = result | (x->x|{v}) \ findPoints I';
    );
    return result;
);
-------------------------------------------------------------------------------
-- Takes in a homogeneous ideal and enumerate points in the projective space
--
findProjPoints = I -> (
    R := ring I;
    k := coefficientRing R;
    x := gens R;
    n := #x - 1;
    result := if AMOUNT then 0 else {};
    -- find all rational points of form (1,...), (0,1,...), (0,0,1,...) etc.
    for i in (0..n) do (
        R' := k(monoid[x_{0..<n-i},MonomialOrder=>Eliminate(n-i-1)]);
        x = gens R';
        v := toList(i:0_k) | {1_k};
        I' := (map(R', R, v|x)) I;
        if AMOUNT then result = result + findPoints I'
        else (
            newPoints := (x->v|x) \ findPoints I';
            if BOUND =!= null and not KEEPALL then newPoints = select(newPoints, x -> ht x <= BOUND);
            result = result | newPoints;
        );
    );
    return result;
);
-------------------------------------------------------------------------------
-- Main interface
-------------------------------------------------------------------------------
-- Takes in an ideal
-- Pre-processing: trim the ideal and reorder the variables
-- Enumeration of the points
-- Post-processing: reconstruct the points according to the original order
-- Returns the list of points or the number of points if `Amount` is set to
-- true
--
rationalPoints = method(Options => {Projective=>false, Amount=>false, Verbose=>false, Split=>false, Bound=>null, KeepAll=>false});
rationalPoints(Ideal) := opts -> I -> (
    (AMOUNT, PROJECTIVE, VERBOSE, SPLIT) = (opts.Amount, opts.Projective, opts.Verbose, opts.Split);
    (BOUND, KEEPALL) = (opts.Bound, opts.KeepAll);
    if VERBOSE then NONSPLITPOLYS = {};
    --
    R := ring I;
    if not isPolynomialRing R then error "expect an ideal of a polynomial ring";
    k := coefficientRing R;
    if not isField k then error "the coefficient ring is not a field";
    ch := char k;
    if ch > 0 then BOUND = null; -- disable BOUND
    if BOUND =!= null then AMOUNT = false; -- disable AMOUNT
    FIELD = k;
    n := numgens R;
    x := gens R;
    result := if AMOUNT then 0 else {};
    --
    if (not PROJECTIVE and dim I <= 0) or (PROJECTIVE and dim I <= 1) then (
        -- In dim 0 use only elimination, which will avoid listing field
        -- elements when the field is big
        ELS = {};
        result = if PROJECTIVE then findProjPoints I else findPoints I;
        if SPLIT and not opts.Amount then result = result / (p -> p / baseChange_FIELD);
    ) else ( -- in positive dimension
        if SPLIT then error "not implemented for positive dimensional ideals";
        if ch > 0 then ( -- finite field
            ELS = listFiniteFieldElements k;
            -- we first get rid of the unused variables; they will be added at the end
            supp := set flatten(for p in first entries gens I list support p);
            ind := apply(n, i -> if member(x_i, supp) then 1 else 0);
            unused := n - sum ind;
            if unused > 0 then (
                unusedR := k(monoid[x_{0..<unused}]);
                R = k(monoid[x_{0..<n-unused}]);
                x = gens R;
                I = (map(R, ring I, shuffle(ind, (unused:0), x))) I;
            );
            -- enumeration of points and post-processing
            if PROJECTIVE then ( -- projective case
                if not isHomogeneous I then error "I is not a homogeneous ideal";
                result = findProjPoints I;
                if unused > 0 then ( -- reconstruction
                    if AMOUNT then result = result*(#ELS)^unused+((#ELS)^unused-1)//(#ELS-1)
                    else (
                        -- shuffle and homogenize so that the first non-zero coordinate is 1
                        result = flatten table(pow_unused ELS, result, homogCoord @@ shuffle_ind);
                        -- extra points lying in a projective subspace
                        PP := findProjPoints(ideal 0_unusedR);
                        result = result | PP / (x->homogCoord shuffle(ind, x, ((n-unused):0_k)));
                    );
                );
            ) else ( -- affine case
                result = findPoints I;
                if unused > 0 then ( -- reconstruction
                    if AMOUNT then result = result*(#ELS)^unused
                    else result = flatten table(pow_unused ELS, result, shuffle_ind);
                );
            );
        ) else ( -- number field
            if BOUND === null then error "need to specify a bound for the height";
            ELS = listNumberFieldElements(k, BOUND);
            if PROJECTIVE then ( -- projective case
                if not isHomogeneous I then error "I is not a homogeneous ideal";
                result = findProjPoints I;
            ) else ( -- affine case
                result = findPoints I;
            );
        );
    );
    if opts.Amount and instance(result, List) then result = #result;
    if not opts.Amount and PROJECTIVE then result = result / (p->new ProjectivePoint from p);
    return result;
);
-------------------------------------------------------------------------------
rationalPoints AffineVariety := opts -> X -> rationalPoints(ideal X, opts, Projective=>false);
rationalPoints ProjectiveVariety := opts -> X -> rationalPoints(ideal X, opts, Projective=>true);
rationalPoints(Ring, Ideal) := opts -> (k, I) -> rationalPoints(baseChange_k I, opts, Split=>false);
rationalPoints(Ring, AffineVariety) := opts -> (k, X) -> rationalPoints(k, ideal X, opts, Projective=>false);
rationalPoints(Ring, ProjectiveVariety) := opts -> (k, X) -> rationalPoints(k, ideal X, opts, Projective=>true);
-------------------------------------------------------------------------------
beginDocumentation()
-------------------------------------------------------------------------------
doc ///
Key
 RationalPoints2
Headline
 Find the rational points on a variety
Description
 Text
  {\em RationalPoints2} is a package for enumerating rational points on a
  variety defined by an ideal of a polynomial ring. The main function is @TO
  rationalPoints@. Over a finite field it will list all the rational points.
  Over a number field it can find all points of a 0-dimensional ideal in its
  splitting field, and all points with bounded height for a positive
  dimensional ideal. The package also provides some utility functions related
  to field extensions, which allow the computation of rational points with
  prescribed coefficient field.
 Text
  An example over a finite field.
 Example
  ZZ/2[x,y,z]; rationalPoints(ideal(y^2*z+y*z^2-x^3-x*z^2), Projective=>true)
 Text
  An example over a number field: we enumerate the rational points on the unit
  circle with bounded height.
 Example
  QQ[x,y]; rationalPoints(ideal(x^2+y^2-1), Bound=>5)
 Text
  And an example of a 0-dimensional ideal.
 Example
  QQ[x,y]; I = ideal(x^2+1,y^2-2);
  rationalPoints(I, Verbose=>true)
  rationalPoints(I, Split=>true, Verbose=>true)
Caveat
 Currently the functionality of positive dimensional ideals over number fields
 other than {\tt QQ} depends on {\tt Sage} (an algorithm by Doyle--Krumm for
 enumerating field elements with bounded height).
Subnodes
 :Main function that enumerates the rational points
  rationalPoints
 :Utility functions for handling field extensions and number fields
  baseChange
  extField
  charpoly
  integers
  hermiteNormalForm
///
-------------------------------------------------------------------------------
doc ///
Key
 baseChange
 (baseChange, Ring, Number)
 (baseChange, Ring, RingElement)
 (baseChange, Ring, Ideal)
 (baseChange, Number, Number)
 (baseChange, Number, RingElement)
 (baseChange, Number, Ideal)
 (baseChange, InexactFieldFamily, Number)
 (baseChange, InexactFieldFamily, RingElement)
 PrimitiveElement
 [baseChange, PrimitiveElement]
Headline
 Perform base change for field extensions
Usage
 baseChange(F, -)
 baseChange(p, -)
Inputs
 F:
  the extension field
 p:
  a prime number
 :
  a polynomial or @ofClass Ideal@ with coefficients in the base field
 PrimitiveElement=>
  specify the image of the primitive element
Outputs
 :
  a polynomial or @ofClass Ideal@ with coefficients in the extension field
Description
 Text
  The base change requires finding an image for the primitive element of the
  base field, which can be specified using the option {\tt PrimitiveElement}.
  The default is {\tt FindOne}, in which case we solve the minimal polynomial
  of the base field in the extension field and take the first zero. Beware that
  it is well-defined only up to the action of the Galois group.
 Text
  In this example we first obtain an extension {\tt F} of the rational numbers
  by adding a cube root of 2, then we compute its image in the splitting field
  {\tt F'} of {\tt x^3-2}.
 Example
  F = toField(QQ[c]/(c^3-2));
  QQ[x]; F' = splittingField(x^3-2);
  c' = baseChange_F' c
  c'^3 == 2
  F[x,y]; I = ideal(x-c, y-c^2);
  baseChange_F' I
 Text
  We get a different embedding by specifying the image of the primitive
  element.
 Example
  F'[x]; c'' = last zeros(x^3-2);
  baseChange(F', I, PrimitiveElement => c'')
 Text
  If a prime number {\tt p} is input instead of a field, the element will be
  reduced to characteristic {\tt p}.
 Example
  QQ[x]; baseChange_7 (x^3-2)
  F[x]; baseChange_7 ideal(x-c)
 Text
  We can also base change an element of a number field to the complex numbers.
  Note that to avoid accumulated errors, it is better to do computations in the
  field and only make the conversion once in the end.
 Example
  baseChange_CC c
  oo^3
  baseChange_CC (c^3)
 Text
  For complex numbers we only need to provide an approximate value for {\tt
  PrimitiveElement}.
 Example
  baseChange(CC, c, PrimitiveElement => -.6+1.1*ii)
///
-------------------------------------------------------------------------------
doc ///
Key
 extField
 (extField, RingElement)
 (extField, Ring, RingElement)
 (extField, List)
 (extField, Ring, List)
 [extField, Variable]
 splittingField
 (splittingField, RingElement)
 (splittingField, Ring, RingElement)
 [splittingField, Variable]
 setPariSize
Headline
 Define field extensions
Usage
 extField p / extField {p, p', ...}
 extField(F, p) / extField(F, {p, p', ...})
 splittingField p
 splittingField(F, p)
Inputs
 p:
  a univariate polynomial or @ofClass List@ of them
 F:
  the base field; if not specified then it is taken to be the field of
  definition of {\tt p}
 Variable=>Symbol
  the variable to use for the primitive element of the extension field
Outputs
 F':
  an extension field over the field of definition of {\tt p}, or over the base
  field {\tt F}
Description
 Text
  Methods for defining field extensions. The method {\tt extField} accepts an
  irreducible polynomial {\tt p} and adds one of its zeros into the base field.
  The base field can itself be a simple extension represented by a primitive
  element, in which case a new primitive element is computed. Use the option
  {\tt Variable} to specify the symbol used for the primitive element (the
  default is {\tt a} like @TO GaloisField@). {\bf N.B. Any data previously
  stored in this symbol will be lost!} As is the case for {\tt GaloisField}.
 Text
  Here we first add {\tt i} into the rational numbers. Similar to @TO toField@,
  the variable used in the polynomial will be set to its image in the
  extension, but note that this is not necessarily the primitive element.
 Example
  QQ[x]; F = extField(4*x^2+1, Variable=>i); x
  baseRing F
  i^2 + 1 == 0
 Text
  Now we perform another field extension over {\tt F} by adding the square root
  of 2. A new primitive element {\tt a} will be computed.
 Example
  QQ[q]; F' = extField(F, q^2-2); q
  q^2 == 2
  baseRing F'
 Text
  For elements of the base field {\tt F}, we can use @TO baseChange@ to
  determine their images in the extension field {\tt F'}.
 Example
  j = baseChange_F' i
  j^2 + 1 == 0
 Text
  It is also possible to do multiple extensions at the same time. Note that
  {\tt x} is now set to {\tt i/2}, so in order to use it as a variable we need
  to use @TO "symbol"@ to reclaim it.
 Example
  QQ[symbol x]; F' = extField {4*x^2+1, x^2-2}; baseRing F'
 Text
  When passing @ofClass List@ of polynomials, the variable {\tt x} will not be
  set. This can also be used for one polynomial to avoid setting the value of
  {\tt x}.
 Example
  extField {4*x^2+1}; x
 Text
  Repeatedly applying {\tt extField} using the same polynomial will add all its
  zeros into the field. This can also be done using the method {\tt
  splittingField}. For example we may construct the splitting field of {\tt
  x^3-2} as follows.
 Example
  QQ[x]; p = x^3-2; F = extField {p}; #zeros_F p
  F' = extField {p, p}; #zeros_F' p
  F' = splittingField p; #zeros_F' p
 Text
  Over finite fields, the degree is first computed while the rest is handled by
  @TO GaloisField@.
 Example
  ZZ/11[x]; F = extField(x^4+1); describe F
  x
  x^4 + 1 == 0
Caveat
 If PARI/GP is available on the system, the program will try to use the
 function {\tt polredbest} to reduce the minimal polynomial. This will give an
 equivalent but usually much optimal result. Note that for large inputs, one
 could increase the pari stack size by using {\tt setPariSize n}.
///
-------------------------------------------------------------------------------
doc ///
Key
 integers
 (integers, Ring)
 (discriminant, Ring)
Headline
 Compute a basis for the integers of a number field
Usage
 integers F
 discriminant F
Inputs
 F:
  a number field
Outputs
 :
  return an integral basis of the algebraic integers of the field {\tt F} as
  @ofClass List@, and compute the discriminant
Description
 Text
  This is an implementation of the Zassenhaus' Round 2 algorithm, following the
  textbook {\em A course in computational algebraic number theory} by Henri
  Cohen (Algorithm 6.1.8). Given a number field $F=\mathbf Q(\theta)$ defined
  as a simple extension over the rational numbers by an integral primitive
  element $\theta$, it computes an integral basis of the algebraic integers,
  expressed as polynomials in $\theta$.
 Example
  integers QQ
  integers toField(QQ[i]/(i^2+1))
  integers toField(QQ[a]/(a^2+3))
 Text
  The discriminant is computed using the same algorithm.
 Example
  discriminant toField(QQ[a]/(a^2+3))
  discriminant toField(QQ[a]/(a^4-a+2))
///
-------------------------------------------------------------------------------
doc ///
Key
 hermiteNormalForm
 (hermiteNormalForm, RingElement)
 (hermiteNormalForm, List)
Headline
 Compute the Hermite normal form of a fractional ideal in a number field
Usage
 hermiteNormalForm a
 hermiteNormalForm {a, a', ...}
Inputs
 a:
  an element of a number field, or @ofClass List@ of them
Outputs
 :Matrix
  representing the Hermite normal form of the fractional ideal generated by
  {\tt a}
Description
 Text
  This method computes the Hermite normal form of a fractional ideal in terms
  of the integral basis given by @TO integers@.
 Example
  F = toField(QQ[q]/(q^2+3));
  hermiteNormalForm 1_F
  hermiteNormalForm(q/3)
///
-------------------------------------------------------------------------------
doc ///
Key
 charpoly
 (charpoly, Number)
 (charpoly, RingElement)
 [charpoly, Variable]
 minpoly
 (minpoly, Number)
 (minpoly, RingElement)
 (minpoly, Ring)
 [minpoly, Variable]
Headline
 Characteristic and minimal polynomials over the prime field
Usage
 charpoly a
 minpoly a
Inputs
 a:
  an element of a field
 Variable=>Symbol
  the variable to use for the polynomial
Outputs
 :
  the characteristic / minimal polynomial of {\tt a} over its prime field
Description
 Text
  Obtain the characteristic / minimal polynomial of an element over its prime
  field.
 Example
  QQ[x]; F = splittingField((x^2+1)*(x^2-2));
  minpoly a
  charpoly(a^2+1, Variable=>y)
  minpoly(a^2+1, Variable=>y)
  GF 81; minpoly(a+1)
 Text
  The method {\tt minpoly} can also be used on a field to recover the
  polynomial used in its definition.
 Example
  minpoly F
///
-------------------------------------------------------------------------------
doc ///
Key
 zeros
 (zeros, RingElement)
 (zeros, Ring, RingElement)
Headline
 List the zeros of a polynomial
Usage
 zeros p
 zeros(F, p)
Inputs
 p:
  a polynomial
 F:
  the coefficient field; if not specified then it is taken to be the field of
  definition of {\tt p}
Outputs
 :
  @ofClass List@ of zeros of {\tt p} in the coefficient field
Description
 Text
  Get @ofClass List@ of zeros of a polynomial using @TO factor@, baby version
  of @TO rationalPoints@.
 Example
  QQ[x]; p = (x-2)^2 * (x^2-2) * (x^3-x-1); zeros p
  F = toField(QQ[q]/(q^2-2)); zeros_F p
  F = splittingField p; #zeros_F p
 Text
  Note that when the degree is big, the expression of each zero in terms of a
  primitive element is usually complicated.
 Example
  last zeros_F p
  (map(F, ring p, {oo})) p
 Text
  Over finite fields.
 Example
  q = baseChange_13 p; zeros q
  F = splittingField q; describe F
  #zeros_F q
///
-------------------------------------------------------------------------------
doc ///
Key
 ProjectivePoint
 "ProjectivePoint == ProjectivePoint"
 (ring, ProjectivePoint)
Headline
 Class of a projective point
Description
 Text
  A subclass of @TO List@, mainly used for pretty printing. Can be constructed
  using @ofClass ProjectiveVariety@ plus @ofClass VisibleList@.
 Example
  P = Proj(QQ[x,y,z]);
  pt = P(1,2,3)
  pt == P{2,4,6}
 Text
  Use {\tt ring} to retrieve the coefficient field.
 Example
  ring pt
///
-------------------------------------------------------------------------------
doc ///
Key
 globalHeight
 (globalHeight, List)
Headline
 Multiplicative height function
Usage
 globalHeight pt
Inputs
 pt:
  @ofClass List@ representing homogeneous coordinates or @ofClass
  ProjectivePoint@
Outputs
 :
  the multiplicative height function
Description
 Text
  This method computes the multiplicative height function given by
  $(x_0:\dots:x_n)\mapsto \prod_{v}\max_i|x_i|_v ^{d_v/d}$. Over rational
  numbers the result will be given as a precise value.
 Example
  globalHeight {1/1,2/3,5/8}
  PP = Proj(QQ[x,y,z]); pt = PP(1/1,2/3,5/8);
  globalHeight pt
  F = toField(QQ[u]/(u^3-5));
  globalHeight {u,u^2/5,1_F}
///
-------------------------------------------------------------------------------
doc ///
Key
 rationalPoints
 (rationalPoints, Ideal)
 (rationalPoints, AffineVariety)
 (rationalPoints, ProjectiveVariety)
 (rationalPoints, Ring, Ideal)
 (rationalPoints, Ring, AffineVariety)
 (rationalPoints, Ring, ProjectiveVariety)
 Amount
 Bound
 KeepAll
 Projective
 Split
 Verbose
 [rationalPoints, Amount]
 [rationalPoints, Bound]
 [rationalPoints, KeepAll]
 [rationalPoints, Projective]
 [rationalPoints, Split]
 [rationalPoints, Verbose]
Headline
 Find the rational points on a variety
Usage
 rationalPoints I
 rationalPoints(F, I)
Inputs
 I:Ideal
  viewed as an affine variety; @ofClass AffineVariety@ / @ofClass
  ProjectiveVariety@ can also be used
 F:
  the coefficient field; if not specified then it is taken to be the field of
  definition of {\tt I}
 Amount=>Boolean
  whether to only compute the number of rational points.
 Projective=>Boolean
  whether to treat the ideal as a homogeneous ideal and consider the
  corresponding projective variety
 Split=>Boolean
  whether to compute all rational points over the splitting field, works only
  if {\tt I} is 0-dimensional (ignored if {\tt F} is specified)
 Verbose=>Boolean
  whether to print the polynomials computed that do not split in the
  coefficient field (if {\tt Split} is set to {\tt true}, it will instead print
  information when a field extension is performed)
 Bound=>Number
  over a number field, specify the upper bound of the multiplicative height for
  the points to enumerate
 KeepAll=>Boolean
  whether to keep the rational points found with height exceeding the bound
Outputs
 :
  Either @ofClass List@ of points, or the number of such points if {\tt Amount}
  is set to true. In the affine case each point is represented by @ofClass
  List@ of coordinates, and in the projective case each point is @ofClass
  ProjectivePoint@.
Description
 Text
  This method takes in an ideal {\tt I} and returns a list of rational points
  on the variety defined by the ideal. The algorithm uses brute force plus some
  elimination. Options are available for different use cases.
 Tree
  :Over finite fields
 Text
  Over a finite field ({\tt ZZ/p} or {\tt GF q}), the ideal can have arbitrary
  dimension. By default it will be regarded as an affine variety.
 Example
  R = ZZ/5[x,y,z]; I = ideal(x^3-y*z, x+y);
  rationalPoints I
  #rationalPoints Spec(R/I)
  #rationalPoints_(GF 25) I
 Text
  To consider projective varieties, one can use the option {\tt
  Projective=>true} or pass directly @ofClass ProjectiveVariety@. In this
  case, the rational points are given in homogeneous coordinates. The first
  non-zero coordinate will be normalized to 1.
 Text
  For example we can take the twisted cubic, which is a smooth rational curve
  so it has q+1 points over a finite field of q elements.
 Example
  ZZ/5[x,y,z,w]; I = ideal "xz-y2,xw-yz,yw-z2";
  rationalPoints(I, Projective => true)
  #rationalPoints variety I
  #rationalPoints_(GF 25) variety I
 Text
  If only the number of rational points is needed, set {\tt Amount} to {\tt
  true} can sometimes speed up the computation.
 Example
  ZZ/101[u_0..u_10]; f = sum toList(u_0..u_10); I = ideal f
  time rationalPoints(I, Amount => true)
 Tree
  :Over number fields
 Text
  Over a number field one can use the option {\tt Bound} to specify a maximal
  multiplicative height given by $(x_0:\dots:x_n)\mapsto \prod_{v}\max_i|x_i|_v
  ^{d_v/d}$ (this is also available as a method @TO globalHeight@).
 Example
  QQ[x,y,z]; I = homogenize(ideal(y^2-x*(x-1)*(x-2)*(x-5)*(x-6)), z);
  rationalPoints(variety I, Bound=>12)
  globalHeight \ oo
 Text
  During the enumeration, sometimes there are rational points found with height
  exceeding the bound. These are cast away by default. Use {\tt KeepAll=>true}
  to keep such points (this will also avoid the heavy task of computing the
  height function, especially over number fields other than {\tt QQ}).
 Example
  rationalPoints(variety I, Bound=>12, KeepAll=>true)
  globalHeight \ oo
 Text
  In the affine case, the height function is used on individual coordinates, so
  for example {\tt {2,1/2}} will be considered as of height 2 instead of 4.
 Example
  QQ[x,y]; rationalPoints(ideal(x-2), Bound=>2)
 Tree
  :Zero-dimensional ideals
 Text
  For a 0-dimensional ideal, all the rational points can be found, including
  those in the splitting field. Set the option {\tt Verbose} to {\tt true} will
  print the polynomials found not splitting in the coefficient field, and set
  {\tt Split} to {\tt true} will automatically make field extensions and
  compute all the points over the splitting field. The information on the
  splitting field can either be printed using {\tt Verbose=>true}, or be
  retrieved from any coordinate of the points. One can refer to @TO extField@
  and @TO baseChange@ for the methods used in handling field extensions.
 Example
  R = QQ[x,y]; I = ideal(x^2+y^2-1,x^3+y^3-1);
  rationalPoints(I, Verbose => true)
  rationalPoints(I, Verbose => true, Split => true)
  ring \ first oo
 Text
  Note that the computations over number fields of large degree are slow. In
  such cases it might be better to switch to a large finite field where the
  program runs much faster. One could also use numerical algorithms implemented
  in the packages @TO "EigenSolver"@ (which uses the same strategy of
  elimination) and @TO "NumericalAlgebraicGeometry"@ (which uses homotopy
  continuation).
 Text
  The next example computes the 31 nodes on a @HREF
  {"https://en.wikipedia.org/wiki/Togliatti_surface", "Togliatti surface"}@.
  Over characteristic 0 the splitting field has degree 8 over {\tt QQ},
  so the computation isn't so heavy.
 Example
  F = toField(QQ[q]/(q^4-10*q^2+20)); R = F[x,y,z,w];
  I = ideal "64(x-w)(x4-4x3w-10x2y2-4x2w2+16xw3-20xy2w+5y4+16w4-20y2w2)-5q(2z-qw)(4(x2+y2-z2)+(1+3(5-q2))w2)2";
  nodes = I + ideal jacobian I;
  time rationalPoints(variety nodes, Split=>true, Verbose=>true);
  #oo
 Text
  Still it runs a lot faster when reduced to a positive characteristic.
 Example
  nodes' = baseChange_32003 nodes;
  time #rationalPoints(variety nodes', Split=>true, Verbose=>true)
Subnodes
 zeros
 ProjectivePoint
 globalHeight
Caveat
 For a number field other than {\tt QQ}, the enumeration of elements with
 bounded height depends on an algorithm by Doyle--Krumm, which is currently
 only implemented in {\tt Sage}.
///
-------------------------------------------------------------------------------
TEST /// -- basic sanity tests
 -- rational points
    R = ZZ/13[t];
    assert(#rationalPoints ideal(0_R) == 13);
    assert(#rationalPoints variety ideal(0_R) == 1);
    assert(#rationalPoints ideal(1_R) == 0);
    assert(#rationalPoints variety ideal(1_R) == 0);
    assert(#rationalPoints ideal(t-2) == 1);
    assert(#rationalPoints ideal((t-2)^100) == 1);
    R = ZZ/5[x_1..x_4];
    I = ideal(x_2^2+x_1*x_2+1, x_1*x_2*x_3*x_4+1);
    assert(#rationalPoints I == 8);
    R = GF 9[m,n,j]; I = ideal(m+a, n*j);
    assert(rationalPoints(I,Amount => true) == 17);
    R = QQ[x]; I = ideal (x^2*(x+1)*(x^2-2));
    assert(rationalPoints(I,Amount => true) == 2);
 -- field extensions
    QQ[x]; F = extField(x^2+4, Variable=>i);
    assert(isField F and x^2+4==0 and i^2+1==0);
    QQ[symbol x]; assert(extField(F, {x^2+1}) === F);
    QQ[y]; F' = extField(F, y^2-2);
    assert(isField F' and y^2==2 and ((baseChange_F' i)^2+1==0));
    QQ[symbol y]; F = splittingField(y^3-2);
    assert(degree baseRing F == 6);
///
-------------------------------------------------------------------------------
TEST ///
 -- number fields
    assert(discriminant QQ == 1);
    assert(discriminant toField(QQ[x]/(x^2+1)) == -4);
    F = toField(QQ[x]/(x^2+3));
    assert(discriminant F == -3 and integers F == {1,(x+1)/2});
    assert(discriminant toField(QQ[x]/(x^4-2*x^3+x^2-5)) == -1975);
    assert(discriminant toField(QQ[x]/(x^4-x+2)) == 2021);
 -- height
    assert(globalHeight{2/1,1/2} == 4);
    F = toField(QQ[x]/(x^3-5));
    assert(abs(globalHeight{x,x^2/5,1_F} - 25^(1/3)) < 1e-10);
    F = toField(QQ[x]/(x^2+3));
    assert(abs(globalHeight{4*x+1,7_F} - 7^(1/2)) < 1e-10);
    R = QQ[x];
    assert(#rationalPoints(ideal 0_R, Bound=>5) == 39);
    R = QQ[x,y,z,w];
    I = ideal "xz-y2,xw-yz,yw-z2";
    assert(#rationalPoints(variety I, Bound=>2) == 4);
    assert(rationalPoints(variety ideal 0_R, Bound=>0.9, Amount=>true) == 0);
///
-------------------------------------------------------------------------------
TEST ///
    assert(#rationalPoints Grassmannian(1,3,CoefficientRing=>ZZ/2) == 36);
    R = ZZ/101[u_0..u_10]; I = ideal sum toList(u_0..u_10);
    assert(rationalPoints(I, Amount => true) == 101^10);
    R = ZZ/101[x,y,z]; I = homogenize(ideal(y^2-x^3-x^2), z);
    assert(#rationalPoints variety I == 101);
    R = QQ[x,y]; I = ideal(x^2+y^2-1,x^3+y^3-1);
    F = toField(QQ[q]/(q^2+2));
    assert(#rationalPoints_F I == 4);
    R = GF 4[x]; I = ideal(x^3 - a);
    assert all(rationalPoints_(GF 64) I / (p -> (x:=p_0; x^6+x^3+1)), zero);
///
-------------------------------------------------------------------------------
TEST ///
    R = QQ[x,y,z,w]; I = Fano_1 ideal(x^3+y^3+z^3+w^3);
    assert(#rationalPoints(variety I, Split=>true) == 27);
///
-------------------------------------------------------------------------------
TEST ///
    needsPackage "Points";
    k = ZZ/101; (m, n) = (20, 3); R = k[x_1..x_n];
    M = random(k^n, k^m);
    mul = toList(1..m) / (i->random(1, 3));
    I = ideal last affineFatPoints(M, mul, R);
    S1 = entries transpose M;
    S2 = rationalPoints I;
    assert(set S1 === set S2);
///
-------------------------------------------------------------------------------
TEST ///
    F = ZZ/32003; S = F[x_0..x_9];
    sigma0 = {(0,2,5),(1,3,6),(2,4,7),(3,0,8),(4,1,9),(0,9,7),(1,5,8),(2,6,9),(3,7,5),(4,8,6)};
    delta = (x,y,v) -> table(10,10,(i,j) -> if i==x and j==y then v else 0);
    skew = (i,j,k) -> sum(delta \ {(i,j,x_k),(j,k,x_i),(k,i,x_j),(j,i,-x_k),(k,j,-x_i),(i,k,-x_j)});
    I = trim pfaffians_6 matrix sum(skew \ sigma0);
    assert(#rationalPoints(variety I, Split=>true) == 55);
///
-------------------------------------------------------------------------------
TEST ///
    QQ[x,y,z,w]; I = ideal "xz-y2,xw-yz,yw-z2";
    assert(#rationalPoints(variety I, Bound=>8) == 8);
///
-------------------------------------------------------------------------------
TEST ///
    if findProgram("sage", "sage -v", RaiseError=>false) =!= null then (
        F = toField(QQ[a]/(a^2+1));
        X = Proj(F[x,y,z]/((4*a+3)*x-5*z));
        assert(#rationalPoints(X, Bound=>3+1e-10) == 30);
    ) else continue;
///
-------------------------------------------------------------------------------
endPackage "RationalPoints2"
-------------------------------------------------------------------------------
end
needsPackage "RationalPoints2"
installPackage("RationalPoints2", RemakeAllDocumentation=>true)
-------------------------------------------------------------------------------
-- Some examples too large to be included as tests
-------------------------------------------------------------------------------
k = ZZ/7823;
R = k[x,y,z,w];
I = ideal(4*x*z + 2*x*w + y^2 + 4*y*w + 7821*z^2 + 7820*w^2,
4*x^2 + 4*x*y + 7821*x*w + 7822*y^2 + 7821*y*w + 7821*z^2 + 7819*z*w + 7820*w^2);
elapsedTime << #rationalPoints variety I; -- 7824
-------------------------------------------------------------------------------
QQ[x]; F = extField {x^3-5};
elapsedTime << #rationalPoints(Proj(F[x,y,z]), Bound=>2.5, Verbose=>true); -- 2557
-------------------------------------------------------------------------------