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
[ ]: