# 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'):

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

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


[ ]: