SageMath code for sampling from affine p-adic manifolds

Main functions

[1]:
import numpy as np

set_verbose(-1)

def init(X):

    n = X.dimension()

    Aff = X.ambient_space()
    gens = list(Aff.gens())
    N = len(gens)

    return [n,gens,N]


def pAdicSolver(eqns, K, package = 'Sage'):

    PolyRing = eqns[0].parent()
    gens = list(PolyRing.gens())
    N = len(gens)

    if package == 'Sage':
        PolyRing = PolynomialRing(QQ, gens, order = 'lex')
        eqns = [PolyRing(eqn) for eqn in eqns]
        GB = ideal(eqns).groebner_basis()
    elif package == 'Magma':
        GB = magma.eval('_<'+','.join(str(e) for e in gens)+'> := PolynomialRing(Rationals(),'+str(N)+'); eqns := '+str(eqns)+'; I := Ideal(eqns); GB := GroebnerBasis(I); print GB;')
        GB = magma('GB').sage()
    else:
        raise TypeError('package must be either Sage or Magma.')

    if list(GB[-1].variables()) != [gens[-1]]:
        raise TypeError('Ideal does not admit a shape position basis, we can not solve the system.')
    for i in range(len(GB) - 1):
        if list(GB[i].variables()) != [gens[i],gens[-1]]:
            raise TypeError('Ideal does not admit a shape position basis, we can not solve the system.')

    last_poly = PolynomialRing(K,[gens[-1]])(GB[-1])
    roots = last_poly.roots(multiplicities = False)

    if len(roots) == 0:
        return []
    else:
        sols = []
        for rt in roots:
            sol = []
            for i in range(len(GB)-1):
                pol = PolynomialRing(QQ,[gens[i],gens[-1]])(GB[i]).change_ring(K)
                sol.append(gens[i] - pol(gens[i],rt))
            sol.append(rt)
            sols.append(sol)
        return sols


"""
def NewPrec(eqns, K, package = 'Sage'):

    sols = pAdicSolver(eqns, K, package)

    vals = []
    for eqn in eqns:
        for sol in sols:
            vals.append(K(eqn(sol)).valuation())
    min_val = min(vals)

    if min_val == infinity:
        new_prec = K.precision_cap()
    else:
        new_prec = min_val

    return new_prec
"""


def AffineTangentSpace(x):

    X = x.scheme()
    n, _, N = init(X)

    TxX = X.tangent_space(x)
    Jac = TxX.Jacobian_matrix().change_ring(X.base_ring()).transpose()
    ker = Jac.kernel()
    ker_basis = ker.basis()
    basis = [list(ker_basis[i]) for i in range(len(ker_basis))]
    W = Matrix(n,N,basis).transpose()
    _,Uinv,_  = W.smith_form()
    U = Uinv^(-1)
    W = U.submatrix(0,0,N,n)

    return W


def IsAnOPoint(v, K):
    min_val = min([K(v[i]).valuation() for i in range(len(v))])
    if min_val >= 0:
        return 1
    else:
        return 0


def AffineWeight(x, K):

    X = x.scheme()
    n, _, N = init(X)
    p = K.prime()

    W = AffineTangentSpace(x)

    val_x = min([K(x[i]).valuation() for i in range(N)])
    norm_x = p^(-val_x)

    diags = [1 for i in range(N-1)] + [p^(max(0,-val_x))]
    S_x = diagonal_matrix(diags)
    Wtilde = S_x*W
    D,_,_ = Wtilde.smith_form()

    diags_D = [p^(-K(D[i][i]).valuation()) for i in range(n)]
    Nr_x = prod(diags_D)

    return max(1,norm_x) * prod([1+p^(-i) for i in range(1,n+1)]) / Nr_x


def AffineOpointWeight(X, K):

    n, _, _ = init(X)
    p = K.prime()

    return prod([1+p^(-i) for i in range(1,n+1)])


def AffineSampleOPoints(X, K, package = 'Sage'):

    n, gens, N = init(X)
    p = K.prime()

    X_proj = X.projective_closure()
    deg = X_proj.degree()

    #print("Thing are set up, entering the while loop \n")

    while True:

        #print("trying\n")

        A = Matrix(QQ, n , N, lambda i, j: choice(range(p^(N+1))))
        b = Matrix(QQ, n , 1, lambda i, j: choice(range(p^(N+1))))

        eqns_planes = A*Matrix(N,gens) - b
        eqns_planes = [eqns_planes[i][0] for i in range(n)]

        eqns_X = list(X.defining_polynomials())
        eqns = eqns_X + eqns_planes

        X_Ab = pAdicSolver(eqns, K, package)

        if len(X_Ab) != 0:

            #new_prec = NewPrec(eqns, K, package)
            #K = Qp(p,new_prec)
            try:
                XK = X.change_ring(K)
                #X_Ab = [[K(X_Ab[i][j]) for j in range(N)] for i in range(len(X_Ab))]

                f_bar = sum([IsAnOPoint(sol,K)*AffineWeight(XK(sol),K) for sol in X_Ab])
                kappa = f_bar / (deg * AffineOpointWeight(X,K))

                #print(kappa, "\n")

                Z = np.random.choice([1,0],1,[kappa, 1 - kappa])[0]

                if Z == 1:

                    O_pt_inds = [i for i in range(len(X_Ab)) if IsAnOPoint(X_Ab[i],K)]
                    i = choice(O_pt_inds)

                    return X_Ab[i]
            except:
                continue

Examples

Elliptic curve

[7]:
Aff.<x,y> = AffineSpace(QQ, 2)

a = 2

b = 4

eqns_X = [y^2 - x^3 - a*x - b]

X = Aff.subscheme(eqns_X)

p, prec = 101, 50
K = Qp(p,prec)

k = K.residue_field()

X_modp = EllipticCurve(k,[a,b])


orders = []


M = 100 # number of sampled points from X

for i in range(M):


    point = AffineSampleOPoints(X,K)
    point_modp = [k(point[0]), k(point[1])]
    P = X_modp(point_modp[0], point_modp[1])


    orders.append(P.order())

orders = np.array(orders)

unique, counts = np.unique(orders, return_counts=True)

D = dict(zip(unique, counts))

D
[7]:
{3: 2,
 5: 5,
 6: 1,
 8: 5,
 10: 8,
 12: 4,
 15: 7,
 20: 5,
 24: 6,
 30: 5,
 40: 17,
 60: 10,
 120: 25}

\(SO(2, \mathcal{O})\)

[11]:
Aff.<a,b,c,d> = AffineSpace(QQ, 4)
eqns_X = [a^2 + b^2 - 1,
          c^2 + d^2 - 1,
          a*c + b*d,
          a*d - b*c - 1]
X = Aff.subscheme(eqns_X)

p, prec = 7, 40
K = Qp(p,prec)

AffineSampleOPoints(X,K)
[11]:
[3*7^4 + 3*7^6 + 3*7^7 + 5*7^8 + 6*7^10 + 6*7^13 + 6*7^14 + 6*7^15 + 2*7^17 + 6*7^18 + 7^19 + 7^20 + 4*7^22 + 7^23 + 5*7^24 + 7^25 + 4*7^27 + 3*7^28 + 6*7^29 + 2*7^30 + 2*7^31 + 2*7^32 + 7^33 + 4*7^36 + 6*7^37 + 2*7^38 + 6*7^39 + O(7^40),
 6 + 6*7 + 6*7^2 + 6*7^3 + 6*7^4 + 6*7^5 + 6*7^6 + 6*7^7 + 4*7^9 + 5*7^10 + 6*7^11 + 6*7^12 + 4*7^13 + 3*7^15 + 2*7^16 + 6*7^17 + 5*7^19 + 7^20 + 5*7^21 + 3*7^22 + 7^23 + 7^24 + 6*7^26 + 4*7^27 + 6*7^28 + 7^29 + 6*7^30 + 3*7^31 + 3*7^32 + 3*7^33 + 7^34 + 4*7^36 + 3*7^37 + 3*7^38 + O(7^40),
 1 + 6*7^8 + 2*7^9 + 7^10 + 2*7^13 + 6*7^14 + 3*7^15 + 4*7^16 + 6*7^18 + 7^19 + 5*7^20 + 7^21 + 3*7^22 + 5*7^23 + 5*7^24 + 6*7^25 + 2*7^27 + 5*7^29 + 3*7^31 + 3*7^32 + 3*7^33 + 5*7^34 + 6*7^35 + 2*7^36 + 3*7^37 + 3*7^38 + 6*7^39 + O(7^40),
 3*7^4 + 3*7^6 + 3*7^7 + 5*7^8 + 6*7^10 + 6*7^13 + 6*7^14 + 6*7^15 + 2*7^17 + 6*7^18 + 7^19 + 7^20 + 4*7^22 + 7^23 + 5*7^24 + 7^25 + 4*7^27 + 3*7^28 + 6*7^29 + 2*7^30 + 2*7^31 + 2*7^32 + 7^33 + 4*7^36 + 6*7^37 + 2*7^38 + 6*7^39 + O(7^40)]

\(SO(3, \mathcal{O})\)

[12]:
Aff.<a,b,c, d,e,f, g,h,i> = AffineSpace(QQ, 9)

eqns_X = [a^2 + b^2 + c^2 - 1,
          d^2 + e^2 + f^2 - 1,
          g^2 + h^2 + i^2 - 1,
          a*d + b*e + c*f,
          a*g + b*h + c*i,
          d*g + e*h + f*i,
          a*(e*i - h*f) - b*(d*i - g*f) + c*(d*h - e*g) - 1]

X = Aff.subscheme(eqns_X)

p, prec = 7, 40
K = Qp(p,prec)

AffineSampleOPoints(X,K)
[12]:
[3*7 + 6*7^3 + 5*7^4 + 2*7^5 + 4*7^6 + 5*7^7 + 4*7^8 + 5*7^10 + 7^11 + 5*7^12 + 2*7^13 + 4*7^14 + 5*7^15 + 3*7^16 + 2*7^17 + 7^18 + 6*7^19 + 5*7^20 + 6*7^21 + 7^25 + 3*7^26 + 2*7^27 + 6*7^28 + 7^29 + 4*7^30 + 4*7^33 + 5*7^34 + 5*7^36 + 5*7^37 + 7^38 + 4*7^39 + O(7^40),
 2 + 3*7 + 7^2 + 7^3 + 3*7^4 + 2*7^5 + 7^6 + 3*7^7 + 6*7^8 + 4*7^9 + 4*7^10 + 2*7^11 + 7^12 + 2*7^13 + 4*7^14 + 4*7^15 + 6*7^16 + 5*7^17 + 5*7^18 + 5*7^19 + 4*7^20 + 2*7^21 + 5*7^23 + 4*7^24 + 4*7^28 + 5*7^29 + 2*7^30 + 4*7^32 + 5*7^33 + 5*7^34 + 7^35 + 6*7^36 + 5*7^37 + 3*7^38 + 7^39 + O(7^40),
 5 + 4*7 + 7^2 + 3*7^3 + 7^4 + 4*7^5 + 5*7^6 + 6*7^7 + 6*7^8 + 5*7^9 + 2*7^10 + 3*7^11 + 2*7^12 + 5*7^13 + 3*7^15 + 7^17 + 3*7^18 + 4*7^19 + 6*7^20 + 5*7^21 + 2*7^22 + 7^23 + 4*7^25 + 2*7^26 + 3*7^27 + 7^28 + 7^29 + 6*7^30 + 4*7^31 + 7^32 + 6*7^33 + 6*7^35 + 2*7^36 + 7^37 + 3*7^38 + 6*7^39 + O(7^40),
 1 + 4*7^2 + 4*7^3 + 5*7^4 + 7^5 + 4*7^6 + 3*7^7 + 5*7^8 + 2*7^9 + 3*7^10 + 3*7^11 + 6*7^12 + 7^15 + 7^17 + 2*7^18 + 3*7^19 + 2*7^20 + 6*7^21 + 5*7^22 + 7^23 + 5*7^24 + 3*7^26 + 2*7^27 + 2*7^28 + 2*7^29 + 3*7^31 + 3*7^33 + 5*7^34 + 3*7^36 + 5*7^37 + O(7^40),
 5*7 + 6*7^2 + 7^3 + 5*7^5 + 2*7^6 + 7^7 + 5*7^8 + 5*7^9 + 2*7^10 + 2*7^11 + 7^12 + 2*7^14 + 2*7^15 + 2*7^18 + 7^19 + 6*7^20 + 5*7^21 + 2*7^22 + 4*7^23 + 7^24 + 6*7^25 + 7^26 + 4*7^27 + 3*7^28 + 5*7^29 + 5*7^30 + 3*7^31 + 4*7^33 + 5*7^34 + 4*7^35 + 7^36 + 3*7^37 + 2*7^38 + 2*7^39 + O(7^40),
 3*7 + 4*7^2 + 4*7^3 + 5*7^4 + 5*7^6 + 4*7^8 + 6*7^9 + 4*7^10 + 5*7^11 + 7^12 + 6*7^13 + 2*7^14 + 7^15 + 7^16 + 2*7^17 + 4*7^18 + 5*7^19 + 3*7^20 + 5*7^21 + 6*7^22 + 5*7^23 + 6*7^24 + 4*7^25 + 5*7^26 + 4*7^27 + 5*7^28 + 7^29 + 6*7^31 + 5*7^33 + 2*7^34 + 2*7^35 + 5*7^36 + 2*7^37 + 7^38 + O(7^40),
 2*7 + 6*7^2 + 4*7^3 + 7^4 + 4*7^5 + 3*7^6 + 6*7^7 + 4*7^8 + 5*7^9 + 3*7^10 + 6*7^11 + 5*7^12 + 4*7^13 + 4*7^14 + 4*7^15 + 7^16 + 7^17 + 7^18 + 5*7^19 + 2*7^20 + 6*7^21 + 7^22 + 3*7^24 + 3*7^25 + 3*7^28 + 4*7^29 + 5*7^30 + 6*7^31 + 3*7^33 + 4*7^34 + 4*7^35 + 5*7^37 + 5*7^38 + 5*7^39 + O(7^40),
 5 + 4*7 + 5*7^2 + 6*7^4 + 6*7^6 + 7^7 + 4*7^8 + 4*7^9 + 2*7^10 + 4*7^11 + 5*7^12 + 6*7^13 + 2*7^14 + 3*7^15 + 3*7^16 + 5*7^17 + 3*7^18 + 7^19 + 4*7^20 + 5*7^21 + 4*7^22 + 3*7^23 + 2*7^24 + 6*7^25 + 3*7^26 + 6*7^27 + 5*7^28 + 4*7^29 + 5*7^30 + 7^31 + 3*7^33 + 6*7^34 + 4*7^35 + 6*7^36 + 7^38 + O(7^40),
 5 + 3*7 + 5*7^2 + 4*7^3 + 3*7^4 + 6*7^5 + 4*7^6 + 6*7^8 + 7^9 + 5*7^10 + 3*7^11 + 6*7^12 + 7^14 + 6*7^15 + 6*7^18 + 5*7^19 + 6*7^20 + 7^21 + 6*7^23 + 6*7^24 + 6*7^25 + 6*7^26 + 5*7^27 + 5*7^28 + 7^29 + 7^30 + 3*7^31 + 6*7^33 + 3*7^35 + 5*7^36 + 5*7^37 + 4*7^38 + 6*7^39 + O(7^40)]

\(O(2, \mathcal{O})\)

[14]:
Aff.<a,b,c,d> = AffineSpace(QQ, 4)
eqns_X = [a^2 + b^2 - 1,
          c^2 + d^2 - 1,
          a*c + b*d]
X = Aff.subscheme(eqns_X)

p, prec = 7, 40
K = Qp(p,prec)

AffineSampleOPoints(X,K)

[14]:
[2 + 6*7 + 6*7^3 + 6*7^4 + 7^5 + 5*7^6 + 7^7 + 7^8 + 2*7^9 + 3*7^11 + 4*7^12 + 5*7^13 + 2*7^14 + 5*7^15 + 7^16 + 4*7^17 + 3*7^18 + 3*7^20 + 7^21 + 2*7^22 + 6*7^24 + 5*7^25 + 2*7^26 + 4*7^27 + 3*7^28 + 6*7^29 + 2*7^30 + 5*7^31 + 2*7^32 + 3*7^33 + 3*7^35 + 6*7^36 + 4*7^37 + 3*7^38 + 7^39 + O(7^40),
 2 + 6*7 + 3*7^2 + 7^3 + 3*7^4 + 4*7^5 + 7^6 + 4*7^7 + 4*7^8 + 2*7^9 + 6*7^10 + 7^12 + 4*7^13 + 3*7^14 + 3*7^15 + 3*7^16 + 6*7^17 + 6*7^18 + 7^20 + 7^21 + 7^22 + 4*7^23 + 6*7^24 + 5*7^25 + 5*7^27 + 7^28 + 7^29 + 3*7^30 + 4*7^31 + 7^32 + 3*7^33 + 4*7^34 + 7^35 + 3*7^36 + 4*7^37 + 6*7^38 + 3*7^39 + O(7^40),
 5 + 3*7^2 + 5*7^3 + 3*7^4 + 2*7^5 + 5*7^6 + 2*7^7 + 2*7^8 + 4*7^9 + 6*7^11 + 5*7^12 + 2*7^13 + 3*7^14 + 3*7^15 + 3*7^16 + 6*7^19 + 5*7^20 + 5*7^21 + 5*7^22 + 2*7^23 + 7^25 + 6*7^26 + 7^27 + 5*7^28 + 5*7^29 + 3*7^30 + 2*7^31 + 5*7^32 + 3*7^33 + 2*7^34 + 5*7^35 + 3*7^36 + 2*7^37 + 3*7^39 + O(7^40),
 2 + 6*7 + 6*7^3 + 6*7^4 + 7^5 + 5*7^6 + 7^7 + 7^8 + 2*7^9 + 3*7^11 + 4*7^12 + 5*7^13 + 2*7^14 + 5*7^15 + 7^16 + 4*7^17 + 3*7^18 + 3*7^20 + 7^21 + 2*7^22 + 6*7^24 + 5*7^25 + 2*7^26 + 4*7^27 + 3*7^28 + 6*7^29 + 2*7^30 + 5*7^31 + 2*7^32 + 3*7^33 + 3*7^35 + 6*7^36 + 4*7^37 + 3*7^38 + 7^39 + O(7^40)]

Picard curves

[15]:
Aff.<x,y> = AffineSpace(QQ, 2)

eqns_X = [y^3 - (x^4 + x + 1)]


X = Aff.subscheme(eqns_X)

p, prec = 7, 100

K = Qp(p,prec)

AffineSampleOPoints(X,K)
[15]:
[5 + 4*7 + 6*7^2 + 7^3 + 6*7^4 + 3*7^5 + 2*7^6 + 4*7^7 + 7^8 + 3*7^9 + 3*7^10 + 4*7^12 + 6*7^13 + 5*7^14 + 4*7^16 + 6*7^17 + 6*7^18 + 6*7^20 + 4*7^21 + 5*7^23 + 4*7^24 + 5*7^25 + 2*7^26 + 3*7^27 + 2*7^28 + 6*7^29 + 6*7^30 + 5*7^31 + 7^32 + 6*7^33 + 5*7^34 + 4*7^35 + 4*7^36 + 6*7^37 + 3*7^38 + 7^39 + 6*7^40 + 3*7^41 + 7^42 + 7^43 + 7^44 + 5*7^45 + 6*7^46 + 4*7^47 + 5*7^48 + 2*7^49 + 4*7^52 + 2*7^53 + 7^55 + 7^56 + 5*7^57 + 7^58 + 4*7^59 + 4*7^62 + 7^63 + 7^64 + 4*7^65 + 3*7^66 + 3*7^67 + 4*7^69 + 2*7^70 + 7^72 + 7^73 + 5*7^74 + 2*7^75 + 2*7^76 + 2*7^77 + 2*7^79 + 3*7^81 + 5*7^82 + 2*7^83 + 6*7^84 + 5*7^85 + 5*7^86 + 4*7^87 + 5*7^88 + 7^89 + 6*7^90 + 2*7^91 + 7^92 + 3*7^93 + 5*7^94 + 5*7^95 + 2*7^96 + 7^97 + 5*7^98 + O(7^100),
 4 + 7 + 6*7^2 + 4*7^3 + 7^4 + 3*7^5 + 6*7^7 + 6*7^9 + 5*7^10 + 4*7^11 + 7^12 + 5*7^13 + 5*7^14 + 4*7^15 + 7^16 + 5*7^17 + 4*7^18 + 2*7^19 + 5*7^20 + 7^21 + 6*7^22 + 3*7^23 + 3*7^24 + 3*7^25 + 6*7^26 + 4*7^27 + 7^29 + 4*7^30 + 7^31 + 3*7^32 + 7^33 + 4*7^34 + 2*7^35 + 2*7^36 + 3*7^37 + 3*7^38 + 5*7^39 + 4*7^40 + 6*7^41 + 3*7^42 + 4*7^43 + 2*7^47 + 5*7^48 + 7^49 + 4*7^50 + 2*7^51 + 5*7^52 + 7^53 + 4*7^54 + 7^55 + 6*7^56 + 3*7^57 + 6*7^58 + 3*7^59 + 3*7^61 + 6*7^62 + 6*7^63 + 5*7^64 + 7^66 + 5*7^67 + 6*7^68 + 5*7^69 + 4*7^70 + 4*7^72 + 6*7^73 + 2*7^74 + 3*7^75 + 2*7^76 + 7^78 + 2*7^79 + 3*7^80 + 4*7^81 + 3*7^82 + 7^83 + 7^84 + 3*7^86 + 3*7^87 + 7^88 + 7^89 + 2*7^90 + 5*7^91 + 4*7^93 + 4*7^94 + 3*7^97 + 5*7^98 + 5*7^99 + O(7^100)]

The symplitic group Sp_1 ( or SL_2 )

[16]:
Aff.<a,b,c,d> = AffineSpace(QQ,4)

eqns_X = [a*d - b*c - 1]

X = Aff.subscheme(eqns_X)

p , prec = 7, 100

K = Qp(p,prec)

AffineSampleOPoints(X,K)
[16]:
[2 + 4*7 + 2*7^4 + 3*7^5 + 6*7^6 + 3*7^7 + 3*7^8 + 4*7^9 + 4*7^10 + 2*7^11 + 7^12 + 6*7^14 + 6*7^15 + 6*7^16 + 7^18 + 6*7^19 + 4*7^20 + 3*7^21 + 7^22 + 7^23 + 5*7^24 + 2*7^25 + 7^26 + 3*7^27 + 7^28 + 5*7^29 + 5*7^30 + 2*7^31 + 3*7^34 + 2*7^35 + 4*7^36 + 5*7^38 + 2*7^39 + 4*7^40 + 3*7^43 + 4*7^45 + 3*7^46 + 3*7^47 + 7^48 + 2*7^49 + 2*7^50 + 2*7^51 + 7^53 + 2*7^54 + 4*7^55 + 6*7^56 + 2*7^58 + 4*7^59 + 6*7^60 + 4*7^61 + 5*7^62 + 7^63 + 2*7^64 + 3*7^65 + 3*7^66 + 2*7^67 + 5*7^68 + 3*7^69 + 5*7^70 + 4*7^71 + 4*7^72 + 7^74 + 6*7^75 + 7^77 + 4*7^78 + 2*7^79 + 3*7^80 + 6*7^81 + 4*7^83 + 7^84 + 7^87 + 6*7^88 + 2*7^89 + 6*7^90 + 6*7^91 + 5*7^92 + 7^93 + 4*7^95 + 2*7^96 + 3*7^97 + 2*7^98 + 6*7^99 + O(7^100),
 1 + 6*7 + 6*7^2 + 4*7^3 + 3*7^4 + 4*7^5 + 2*7^6 + 6*7^7 + 2*7^9 + 5*7^10 + 4*7^11 + 2*7^12 + 6*7^13 + 4*7^14 + 6*7^15 + 7^16 + 3*7^17 + 3*7^18 + 4*7^19 + 6*7^20 + 3*7^21 + 3*7^22 + 7^24 + 2*7^26 + 4*7^28 + 5*7^29 + 2*7^30 + 7^31 + 7^32 + 3*7^33 + 5*7^34 + 6*7^35 + 2*7^36 + 6*7^37 + 7^38 + 3*7^39 + 7^40 + 4*7^41 + 3*7^42 + 7^43 + 3*7^44 + 2*7^45 + 2*7^48 + 6*7^49 + 6*7^50 + 2*7^51 + 2*7^52 + 2*7^54 + 5*7^55 + 5*7^56 + 4*7^57 + 3*7^58 + 6*7^59 + 4*7^60 + 4*7^61 + 4*7^62 + 5*7^63 + 3*7^64 + 4*7^65 + 5*7^67 + 3*7^68 + 6*7^69 + 7^70 + 7^71 + 4*7^72 + 2*7^73 + 5*7^74 + 7^75 + 6*7^76 + 5*7^77 + 5*7^78 + 5*7^80 + 2*7^81 + 6*7^84 + 4*7^86 + 7^87 + 6*7^88 + 7^89 + 2*7^90 + 4*7^91 + 5*7^92 + 6*7^93 + 6*7^94 + 6*7^95 + 6*7^97 + 3*7^98 + 2*7^99 + O(7^100),
 2 + 7 + 3*7^3 + 3*7^4 + 7^5 + 4*7^6 + 2*7^7 + 6*7^9 + 7^10 + 6*7^11 + 5*7^13 + 3*7^14 + 7^15 + 7^16 + 6*7^17 + 7^18 + 4*7^19 + 5*7^21 + 7^22 + 6*7^23 + 4*7^24 + 7^25 + 2*7^26 + 2*7^27 + 3*7^28 + 4*7^29 + 3*7^30 + 6*7^31 + 4*7^32 + 7^33 + 4*7^34 + 3*7^35 + 6*7^36 + 3*7^37 + 2*7^38 + 6*7^39 + 2*7^40 + 5*7^41 + 5*7^43 + 4*7^44 + 5*7^45 + 3*7^46 + 3*7^47 + 5*7^48 + 4*7^49 + 7^50 + 2*7^51 + 4*7^52 + 7^53 + 2*7^54 + 4*7^55 + 6*7^57 + 5*7^58 + 3*7^59 + 5*7^60 + 5*7^61 + 6*7^62 + 6*7^65 + 3*7^66 + 6*7^67 + 2*7^68 + 7^70 + 4*7^71 + 3*7^72 + 7^73 + 5*7^74 + 3*7^75 + 3*7^77 + 7^78 + 6*7^79 + 3*7^80 + 4*7^81 + 2*7^82 + 3*7^83 + 7^84 + 5*7^85 + 4*7^86 + 3*7^87 + 3*7^88 + 5*7^89 + 5*7^91 + 7^92 + 4*7^93 + 6*7^94 + 2*7^95 + 7^96 + 2*7^97 + 6*7^98 + 4*7^99 + O(7^100),
 5 + 3*7 + 2*7^2 + 7^3 + 2*7^4 + 4*7^6 + 5*7^7 + 2*7^8 + 2*7^9 + 2*7^11 + 5*7^13 + 2*7^14 + 6*7^16 + 5*7^17 + 6*7^18 + 2*7^19 + 3*7^20 + 2*7^22 + 7^23 + 4*7^24 + 5*7^25 + 5*7^27 + 5*7^28 + 3*7^30 + 6*7^31 + 2*7^32 + 5*7^34 + 5*7^35 + 6*7^36 + 3*7^37 + 3*7^38 + 7^39 + 2*7^40 + 6*7^41 + 3*7^43 + 2*7^44 + 7^45 + 7^46 + 5*7^47 + 4*7^48 + 2*7^49 + 4*7^51 + 6*7^52 + 6*7^53 + 5*7^55 + 2*7^56 + 6*7^57 + 2*7^58 + 3*7^59 + 7^60 + 7^61 + 2*7^62 + 2*7^63 + 7^64 + 7^66 + 4*7^68 + 6*7^69 + 4*7^70 + 3*7^71 + 6*7^74 + 4*7^75 + 7^76 + 6*7^77 + 7^78 + 7^80 + 4*7^82 + 5*7^83 + 6*7^84 + 3*7^85 + 3*7^86 + 6*7^88 + 2*7^89 + 6*7^91 + 7^93 + 3*7^94 + 3*7^95 + 2*7^96 + 5*7^98 + 7^99 + O(7^100)]

The “sphere”

[17]:
Aff.<x,y,z> = AffineSpace(QQ,3)

eqns_X = [x^2 + y^2 + z^2 - 1]

X = Aff.subscheme(eqns_X)

p , prec = 7, 100

K = Qp(p,prec)

AffineSampleOPoints(X,K)
[17]:
[4 + 2*7 + 7^2 + 3*7^3 + 7^4 + 3*7^5 + 5*7^6 + 2*7^7 + 7^8 + 2*7^9 + 2*7^10 + 2*7^11 + 5*7^12 + 6*7^14 + 7^15 + 3*7^16 + 4*7^17 + 3*7^18 + 4*7^19 + 2*7^20 + 5*7^21 + 2*7^22 + 4*7^23 + 2*7^24 + 6*7^25 + 2*7^26 + 5*7^27 + 7^28 + 7^29 + 4*7^30 + 7^31 + 3*7^32 + 7^33 + 5*7^34 + 6*7^35 + 5*7^36 + 5*7^39 + 2*7^40 + 7^41 + 5*7^42 + 3*7^43 + 5*7^44 + 2*7^45 + 7^46 + 7^47 + 2*7^48 + 2*7^49 + 2*7^50 + 4*7^51 + 2*7^52 + 6*7^53 + 5*7^54 + 5*7^55 + 7^56 + 5*7^57 + 3*7^58 + 4*7^59 + 3*7^60 + 7^61 + 6*7^62 + 2*7^63 + 6*7^64 + 5*7^65 + 3*7^66 + 6*7^67 + 4*7^68 + 6*7^69 + 5*7^70 + 7^71 + 4*7^72 + 6*7^73 + 7^74 + 5*7^75 + 4*7^76 + 4*7^77 + 2*7^78 + 7^79 + 2*7^80 + 4*7^81 + 5*7^82 + 5*7^83 + 2*7^84 + 3*7^86 + 6*7^87 + 3*7^88 + 7^89 + 7^90 + 7^91 + 4*7^92 + 5*7^93 + 7^94 + 6*7^95 + 3*7^96 + 2*7^97 + 5*7^99 + O(7^100),
 4 + 3*7 + 3*7^2 + 5*7^3 + 7^4 + 3*7^5 + 7^6 + 5*7^7 + 4*7^8 + 7^9 + 3*7^11 + 6*7^12 + 7^13 + 4*7^14 + 6*7^15 + 3*7^16 + 3*7^17 + 3*7^18 + 6*7^19 + 6*7^20 + 3*7^21 + 4*7^22 + 5*7^23 + 5*7^24 + 6*7^25 + 4*7^26 + 2*7^27 + 5*7^28 + 5*7^29 + 7^30 + 3*7^32 + 6*7^33 + 5*7^34 + 3*7^35 + 2*7^36 + 3*7^37 + 3*7^38 + 6*7^39 + 2*7^40 + 4*7^41 + 3*7^43 + 6*7^44 + 6*7^45 + 6*7^46 + 3*7^47 + 2*7^48 + 7^49 + 5*7^50 + 7^51 + 2*7^52 + 2*7^53 + 6*7^54 + 6*7^55 + 5*7^56 + 7^57 + 7^58 + 6*7^60 + 6*7^61 + 4*7^62 + 4*7^63 + 3*7^64 + 7^66 + 7^67 + 4*7^68 + 3*7^69 + 4*7^70 + 6*7^72 + 4*7^73 + 2*7^74 + 4*7^75 + 6*7^76 + 3*7^77 + 4*7^78 + 3*7^79 + 7^80 + 3*7^82 + 6*7^84 + 7^86 + 6*7^87 + 3*7^88 + 5*7^89 + 3*7^90 + 4*7^91 + 4*7^92 + 2*7^93 + 5*7^94 + 3*7^95 + 6*7^96 + 4*7^98 + 5*7^99 + O(7^100),
 2 + 7 + 6*7^2 + 6*7^3 + 3*7^4 + 6*7^5 + 6*7^6 + 2*7^7 + 2*7^8 + 3*7^9 + 5*7^11 + 7^12 + 4*7^13 + 4*7^14 + 4*7^15 + 7^16 + 6*7^17 + 2*7^18 + 2*7^19 + 4*7^20 + 4*7^21 + 3*7^22 + 4*7^23 + 6*7^25 + 3*7^26 + 6*7^27 + 4*7^28 + 6*7^29 + 3*7^30 + 4*7^31 + 7^32 + 2*7^33 + 2*7^35 + 3*7^36 + 6*7^37 + 5*7^39 + 6*7^40 + 7^41 + 6*7^42 + 7^44 + 5*7^46 + 5*7^47 + 7^49 + 5*7^51 + 5*7^53 + 7^54 + 6*7^55 + 2*7^58 + 2*7^60 + 3*7^61 + 7^62 + 5*7^63 + 7^64 + 6*7^65 + 5*7^66 + 2*7^68 + 6*7^69 + 3*7^71 + 5*7^72 + 7^73 + 3*7^74 + 3*7^75 + 3*7^76 + 2*7^77 + 6*7^78 + 3*7^79 + 7^80 + 7^81 + 6*7^82 + 4*7^83 + 3*7^84 + 4*7^85 + 3*7^87 + 6*7^88 + 4*7^90 + 4*7^91 + 2*7^92 + 5*7^93 + 4*7^94 + 3*7^95 + 5*7^96 + 4*7^97 + 4*7^99 + O(7^100)]

The symplitic group Sp_2 (too slow for Groebner basis with sage)

[36]:
Aff.<x11,x12,x13,x14, x21,x22,x23,x24, x31,x32,x33,x34, x41,x42,x43,x44 > = AffineSpace(QQ,16)

eqns_X =  [x11*x33 - x13*x31 + x21*x43 - x23*x41 - 1,
           x12*x34 - x14*x32 + x22*x44 - x24*x42 - 1,
           x11*x32 - x12*x31 + x21*x42 - x22*x41,
           x12*x33 - x13*x32 + x22*x43 - x23*x42,
           x11*x34 - x14*x31 + x21*x44 - x24*x41,
           x13*x34 - x14*x33 + x23*x44 - x24*x43]

X = Aff.subscheme(eqns_X)

p , prec = 7, 100

K = Qp(p,prec)

[36]:
Closed subscheme of Affine Space of dimension 16 over Rational Field defined by:
  -x12*x31 + x11*x32 - x22*x41 + x21*x42,
  -x13*x31 + x11*x33 - x23*x41 + x21*x43 - 1,
  -x13*x32 + x12*x33 - x23*x42 + x22*x43,
  -x14*x31 + x11*x34 - x24*x41 + x21*x44,
  -x14*x32 + x12*x34 - x24*x42 + x22*x44 - 1,
  -x14*x33 + x13*x34 - x24*x43 + x23*x44

The modular curve X_1(30)

[ ]:
"""
X, par1, par2, E: as in https://math.mit.edu/~drew/X1_optcurves.html
M : number of tries
"""

def TamagawaNumbers(X, par1, par2, E, p, prec, M):

    X = Aff.subscheme(X)
    K = Qp(p,prec)

    _ = magma.eval(' _<x,y> := PolynomialRing(Rationals(),2); par1_fc := '+str(par1)+'; par2_fc := '+str(par2)+'; E_fc := '+str(E)+'; ')
    _ = magma.eval(' p := '+str(p)+'; prec := '+str(prec)+'; K := pAdicField(p,prec); ')

    TamagawaNumbersList = []

    print("All is set up. Begining computations")

    for i in range(M):

        try:
            print("Sampling a point ...")
            pt = AffineSampleOPoints(X,K)

            _ = magma.eval(' x0 := K!('+str(QQ(pt[0]))+'); y0 := K!('+str(QQ(pt[1]))+'); ')
            _ = magma.eval(' E := [Evaluate(elt,[x0,y0]) : elt in E_fc]; ')
            _ = magma.eval(' E := EllipticCurve(E); ')
        except:
            continue

        _ = magma.eval(' _,_,_,cp,_,_ := Explode(LocalInformation(E)); ')
        cp = magma('cp').sage()
        TamagawaNumbersList.append(cp)

    return TamagawaNumbersList

####################################################################################

Aff.<x,y> = AffineSpace(QQ, 2)

N = 30
X = y^6 \
    + (x^6 - 5*x^5 + 6*x^4 + 3*x^3 - 6*x^2 + 7*x + 3)*y^5 \
    + (x^7 - 3*x^6 - 13*x^5 + 44*x^4 - 18*x^3 + x^2 + 18*x + 3)*y^4 \
    + (x^8 - 3*x^7 - 13*x^6 + 27*x^5 + 46*x^4 - 32*x^3 + 21*x^2 + 15*x + 1)*y^3 \
    + 2*x*(x^7 - 8*x^6 + 9*x^5 + 20*x^4 + 6*x^3 - 6*x^2 + 9*x + 2)*y^2 \
    - 4*x^2*(2*x^5 - 7*x^4 - 3*x^3 - 1)*y \
    + 8*x^6
q = y+1
t = 4*(y+1)*(x+y)/(x*y^3 - 4*x*y - 4*x - 3*y^3 - 6*y^2 - 4*y)
E = [0,t^2-2*q*t-2,0,-(t^2-1)*(q*t+1)^2,0]

p = 31
prec = 250
M = 130

TamNumsList = TamagawaNumbers(X, q, t, E, p, prec, M)
print('#########################################################################################')
print('  Prime:',p, '- Precision:',prec, '- Level:',N, '----- Number of tries:', M, '- Number of fails:', M - len(TamNumsList))
print('#########################################################################################')
for TamNum in set(TamNumsList):
    print('Tamagawa number:', TamNum,'- Number of times it appears:', TamNumsList.count(TamNum))

Hilbert Modular Surfaces

[11]:
def IgusaClebschInvariants(D,K,cord1,cord2): # Elkies-Kumar

    if D == 5:
        g, h = K(cord1), K(cord2)
        A1 = 1
        A = -3*g^2/4
        B1 = -(g+1/4)
        B = -(g^3/4 + h)
        B2 = -h^2
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 8:
        r, s = K(cord1), K(cord2)
        A1 =  2*r*s^2
        A = -(9*r*s+4*r^2+4*r+1)/3
        B1 = r*s^2*(3*s+8*r-2)/3
        B = -(54*r^2*s+81*r*s-16*r^3-24*r^2-12*r-2)/27
        B2 = r^2
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 12:
        e, f = K(cord1), K(cord2)
        A1 = (f-1)
        A = -(f^4+15*e*f+9*e)/3
        B1 = (2*f^3-2*f^2-3*f+3*e+3)/3
        B = (2*f^6-63*e*f^3-81*e*f^2-54*e^2)/27
        B2 = e^3
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 13:
        g, h = K(cord1), K(cord2)
        A1 = -1
        A = (48*g*h-4*g^2+4*g-1)/3
        B1 = (3*h^2-6*h+4*g+1)/3
        B =  -2*(108*g^2*h^2+180*g^2*h-144*g*h+8*g^3-12*g^2+6*g-1)/27
        B2 = 16*g^4*h^2
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 17:
        g, h = K(cord1), K(cord2)
        A1 = (2*g+1)^2
        A =  -(16*h^2-8*g^2*h+20*g*h+64*h+g^4+4*g^3+6*g^2+4*g+1)/3
        B1 = -(48*h^2-16*g^2*h-40*g*h-16*h+4*g^4-12*g^3-23*g^2-12*g-2)/3
        B = -2*(64*h^3-48*g^2*h^2-312*g*h^2+978*h^2+12*g^4*h-6*g^3*h+72*g^2*h+210*g*h \
            +120*h-g^6-6*g^5-15*g^4-20*g^3-15*g^2-6*g-1)/27
        B2 =  -64*h^3
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 21:
        r, s = K(cord1), K(cord2)
        A1 = -1
        A = -(s-r)^2*(s^2+(24*r^3+30*r^2-26*r-30)*s-15*r^4-30*r^3+7*r^2+30*r+9)/3
        B1 = (r^4-2*s*r^3+(s^2-5/3)*r^2+10/3*s*r+(-2/3*s^2+2*s+1))
        B = 2*(s-r)^4*(s^2 + (-72*r^3-63*r^2+70*r+63)*s \
            -27*r^6-189*r^5-63*r^4+441*r^3+280*r^2-252*r-189)/27
        B2 = (r-1)^6 * (r+1)^4 * (s-r)^6
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 24:
        a, d = K(cord1), K(cord2)
        A1 = 2*a^2*(d-1)^2*(d^2-a-1)
        A = (3*a*d^2 + 10*a*d + (-3*a^2 + 5*a - 1/3))
        B1 = (2*a^2/3)*(-3*d^6 + 6*d^5 + (6*a+2)*d^4 - (6*a+10)*d^3 \
            - (9*a^2+5*a-3)*d^2 + (6*a^2+4*a+4)*d + (3*a^3+3*a^2+a-2))
        B = ((-11*a^2+3*a)*d^2 + (-14*a^2+14/3*a)*d + (-5*a^2+7/3*a+2/27))
        B2 = (a/2)
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 28:
        f, g = K(cord1), K(cord2)
        A1 = -4*(f+1)^2*(g-f-2)*(g-f)^3*(g+f)^4
        A = -(25*g^4+30*g^3+(-50*f^2-18*f+8)*g^2 \
            -30*f^2*g+25*f^4+18*f^3-8*f^2+1)/3
        B1 = (g-f)^3*(g+f)^4*(3*g^5-3*f*g^4-6*g^4+10*f^2*g^3+20*f*g^3+4*g^3-10*f^3*g^2-40*f^2*g^2 \
             -44*f*g^2-8*g^2+35*f^4*g+172*f^3*g+288*f^2*g+200*f*g+52*g-35*f^5 \
             -146*f^4-248*f^3-200*f^2-68*f-8)/3
        B = -(196*g^6+504*g^5-372*f^2*g^4-216*f*g^4+201*g^4-1008*f^2*g^3+126*g^3 \
            +156*f^4*g^2+432*f^3*g^2-402*f^2*g^2-162*f*g^2-24*g^2+504*f^4*g \
            -126*f^2*g+20*f^6-216*f^5+201*f^4+162*f^3+24*f^2-2)/27
        B2 = 1
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    elif D == 29:
        f, g = K(cord1), K(cord2)
        A1 = 8*f^2*g^2*(g-4*f)^2;
        A = -(9*g^2+72*f*g-54*g+16*f^2-8*f+1)/3
        B1 = -2*g^2*(3*g^4-36*f*g^3+48*f^4*g^2-16*f^3*g^2+148*f^2*g^2-64*f^4*g \
             -224*f^3*g+512*f^5+64*f^4)/3
        B = -2*(108*f^2*g^2-162*f*g^2+405*g^2-648*f^2*g+648*f*g-135*g+64*f^3 \
            -48*f^2+12*f-1)/27
        B2 = -1/2
        I2,I4,I6,I10 = [-24*B1/A1, -12*A, 96*(A/A1)*B1-36*B, -4*A1*B2]

    else:
        return "This discriminant is not implemented"

    return [I2,I4,I6,I10]

def LinearForms(IguCleInv): # Helminck

    I2,I4,I6,I10 = IguCleInv

    J2  = I2/8
    J4  = (4*J2^2 - I4)/96
    J6  = (8*J2^3 - 160*J2*J4 - I6)/576
    J8  = (J2*J6 - J4^2)/4
    J10 = I10/4096 # Igusa Invariants (LMFDB: https://www.lmfdb.org/knowledge/show/g2c.igusa_invariants)

    I2  = J2/12
    I4  = J2^2 - 24*J4
    I6  = J6
    I8  = J8
    I12 = -8*J4^3 + 9*J2*J4*J6 - 27*J6^2 - J2^2*J8 # Extended Igusa Invariants

    w11 = 5*J2.valuation()  - 1*J10.valuation()
    w12 = 5*J4.valuation()  - 2*J10.valuation()
    w13 = 5*J6.valuation()  - 3*J10.valuation()
    w14 = 5*J8.valuation()  - 4*J10.valuation()
    w15 = 0 # 5*J10.valuation() - 5*J10.valuation()
    w1  = [w11,w12,w13,w14,w15]

    w21 = 6*J2.valuation()  - 1*I12.valuation()
    w22 = 6*J4.valuation()  - 2*I12.valuation()
    w23 = 6*J6.valuation()  - 3*I12.valuation()
    w24 = 6*J8.valuation()  - 4*I12.valuation()
    w25 = 6*J10.valuation() - 5*I12.valuation()
    w2  = [w21,w22,w23,w24,w25]

    w31  = 2*J2.valuation()  - 1*I4.valuation()
    w32  = 2*J4.valuation()  - 2*I4.valuation()
    w33  = 2*J6.valuation()  - 3*I4.valuation()
    w34  = 2*J8.valuation()  - 4*I4.valuation()
    w35  = 2*J10.valuation() - 5*I4.valuation()
    w3x  = I12.valuation()   - 3*I4.valuation()
    w3y1 = J4.valuation()    - I4.valuation()
    w3y2 = 2*J6.valuation()  - 3*I4.valuation()
    w3   = [w31,w32,w33,w34,w35,w3x,w3y1,w3y2]

    w41 = w31
    w42 = w32
    w43 = w33
    w44 = w34
    w45 = w35
    w4  = [w41,w42,w43,w44,w45]

    w51 = 3*I4.valuation() - J10.valuation() - I2.valuation()
    w52 = I12.valuation()  - J10.valuation() - I2.valuation()
    w5  = [w51,w52]

    w61 = 3*I4.valuation() - I12.valuation()
    w62 = J10.valuation()  + I2.valuation() - I12.valuation()
    w6  = [w61,w62]

    w71 = -w61
    w72 = J10.valuation() + I2.valuation() - 3*I4.valuation()
    w7  = [w71,w72]

    w2c1 = I4.valuation()  - 2*I2.valuation()
    w2c2 = J10.valuation() - 5*I2.valuation()
    w2c3 = I12.valuation() - 6*I2.valuation()
    w2c  = [w2c1,w2c2,w2c3]

    return w1, w2, w3, w4, w5, w6, w7, w2c

def ReductionType(D,K,cord1,cord2):

    w1, w2, w3, w4, w5, w6, w7, w2c = LinearForms(IgusaClebschInvariants(D,K,cord1,cord2))

    w11,w12,w13,w14,w15               = w1
    w21,w22,w23,w24,w25               = w2
    w31,w32,w33,w34,w35,w3x,w3y1,w3y2 = w3
    w41,w42,w43,w44,w45               = w4
    w51,w52                                                   = w5
    w61,w62                                                   = w6
    w71,w72                                                   = w7
    w2c1,w2c2,w2c3                                    = w2c

    condI      = w11 >= 0 and w12 >= 0 and w13 >= 0 and w14 >= 0 and w15 >= 0
    condII     = w21 >= 0 and w22 >= 0 and w23 >= 0 and w24 >= 0 and w25 > 0
    condIII    = w31 >= 0 and w32 >= 0 and w33 >= 0 and w34 >= 0 and w35 > 0 \
                 and w3x >= 0 and (w3y1 == 0 or w3y2 == 0)
    condIV     = w42 > 0 and w43 > 0 and w44 > 0 and w45 > 0
    condV      = w2c1 > 0 and w2c2 > 0 and w2c3 > 0 and w51 >= 0 and w52 >= 0
    condVI     = w2c1 > 0 and w2c2 > 0 and w2c3 > 0 and w61 >= 0 and w62 > 0
    condVII    = w2c1 > 0 and w2c2 > 0 and w2c3 > 0 and w71 > 0 and w72 > 0
    conditions = [condI,condII,condIII,condIV,condV,condVI,condVII]

    if conditions.count(True) == 0:
        return 'None of the conditions are satisfied, something is wrong'
    elif conditions.count(True) > 1:
        return 'Conditions are not exclusive, something is wrong'

    elif condI:
        return 'I'
    elif condII:
        return 'II'
    elif condIII:
         return 'III'
    elif condIV:
        return 'IV'
    elif condV:
        return 'V'
    elif condVI:
        return 'VI'
    elif condVII:
        return 'VII'

def ReductionTypes(D, p, prec, M):

    if D == 5:
        Aff.<g,h,z> = AffineSpace(QQ, 3)
        X = z^2 - 2*(6250*h^2-4500*g^2*h-1350*g*h-108*h-972*g^5-324*g^4-27*g^3)

    elif D == 8:
        Aff.<r,s,z> = AffineSpace(QQ, 3)
        X = z^2 - 2*(16*r*s^2+32*r^2*s-40*r*s-s+16*r^3+24*r^2+12*r+2)

    elif D == 12:
        Aff.<e,f,z> = AffineSpace(QQ, 3)
        X = z^2 - (f-1)*(f+1)*(f^6-f^4-18*e*f^2+27*e^2+16*e)

    elif D == 13:
        Aff.<g,h,z> = AffineSpace(QQ, 3)
        X = z^2 - (108*g*h^3-27*g^2*h^2-468*g*h^2+4*h^2+656*g^2*h+568*g*h-16*h-128*g^3+192*g^2-96*g+16)

    elif D == 17:
        Aff.<g,h,z> = AffineSpace(QQ, 3)
        X = z^2 + (256*h^3-192*g^2*h^2-464*g*h^2-185*h^2+48*g^4*h-236*g^3*h-346*g^2*h \
                - 144*g*h-18*h-4*g^6-20*g^5-41*g^4-44*g^3-26*g^2-8*g-1)

    elif D == 21:
        Aff.<r,s,z> = AffineSpace(QQ, 3)
        X = z^2 - 16*s^4  -  8*r*(27*r^2-23) * s^3  +  (621*r^4-954*r^2+349)*s^2 \
                - 18*(r^3-r)*(33*r^2-29) * s  +  (r^2-1) * (189*(r^4-r^2)+16)

    elif D == 24:
        Aff.<a,d,z> = AffineSpace(QQ, 3)
        X = z^2 - (d^2 - a - 1)*(16*a*d^4-8*a^2*d^2-20*a*d^2+d^2+a^3-3*a^2+3*a-1)

    elif D == 28:
        Aff.<f,g,z> = AffineSpace(QQ, 3)
        X = z^2 + (g-f-2)*(g+f+2)*(8*g^4+92*f^2*g^2+180*f*g^2+71*g^2-100*f^4-180*f^3-71*f^2+4*f+4)

    elif D == 29:
        Aff.<f,g,z> = AffineSpace(QQ, 3)
        X = z^2 + (g^4-6*f*g^3-11*g^3+27*f^4*g^2-18*f^3*g^2+5*f^2*g^2+102*f*g^2-g^2 \
                - 288*f^4*g+200*f^3*g-280*f^2*g+8*f*g+1024*f^5-768*f^4+192*f^3-16*f^2)
    else:
        return "This discriminant is not implemented"

    X = Aff.subscheme(X)
    K = Qp(p,prec)

    ReductionTypesList = []

    for i in range(M):

        if i % 10000 == 0:
            print(str(i)+"\n")

        try:
            cord1,cord2,_ = AffineSampleOPoints(X,K)
            ReductionTypesList.append(ReductionType(D,K,cord1,cord2))
        except:
            continue

    return ReductionTypesList

####################################################################################

# D must be one of 5, 8, 12, 13, 17, 21, 24, 28, 29. The list can be extended easily.

D = 5 # discriminant
p = 5 # prime
prec = 1000 # precision
M = 20 # number of tries

RedTypesList = ReductionTypes(D, p, prec, M)

print('###################################################################################################')
print('  Discriminant:',D, '- Prime:',p, '- Precision:',prec, '----- Number of tries:', M, '- Number of fails:', M - len(RedTypesList))
print('###################################################################################################')
for RedType in set(RedTypesList):
    print('Reduction Type:', RedType,'- Number of times it appears:', RedTypesList.count(RedType))
print('\n')



0

###################################################################################################
  Discriminant: 5 - Prime: 5 - Precision: 1000 ----- Number of tries: 20 - Number of fails: 0
###################################################################################################
Reduction Type: V - Number of times it appears: 1
Reduction Type: I - Number of times it appears: 17
Reduction Type: III - Number of times it appears: 2


[ ]: