Skip to content

Commit

Permalink
work on Simon two-descent
Browse files Browse the repository at this point in the history
  • Loading branch information
@fchapoton
fchapoton committed Sep 19, 2024
1 parent a880b3f commit c999704
Show file tree
Hide file tree
Showing 3 changed files with 75 additions and 106 deletions.
9 changes: 0 additions & 9 deletions src/sage/schemes/elliptic_curves/ell_number_field.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,15 +250,6 @@ def simon_two_descent(self, verbose=0, lim1=2, lim3=4, limtriv=2,
listpoints = [[Mod(1/2*y + 3/2, y^2 + 7), Mod(-y - 2, y^2 + 7), 1]]
(1, 1, [(1/2*a + 3/2 : -a - 2 : 1)])
sage: v = E.simon_two_descent(verbose=2)
K = bnfinit(y^2 + 7);
a = Mod(y,K.pol);
bnfellrank(K, [0, 0, 0, 1, a], [[Mod(1/2*y + 3/2, y^2 + 7), Mod(-y - 2, y^2 + 7)]]);
...
v = [1, 1, [[Mod(1/2*y + 3/2, y^2 + 7), Mod(-y - 2, y^2 + 7)]]]
sage: v
(1, 1, [(1/2*a + 3/2 : -a - 2 : 1)])
A curve with 2-torsion::
sage: K.<a> = NumberField(x^2 + 7)
Expand Down
69 changes: 32 additions & 37 deletions src/sage/schemes/elliptic_curves/ell_rational_field.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1837,6 +1837,10 @@ def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=3,
...
DeprecationWarning: Use E.rank(algorithm="pari") instead, as this script has been ported over to pari.
See https://github.com/sagemath/sage/issues/35621 for details.
doctest:warning
...
DeprecationWarning: please use the 2-descent algorithm over QQ inside pari
See https://github.com/sagemath/sage/issues/38461 for details.
(0, 0, [])
sage: E = EllipticCurve('37a1')
sage: E.simon_two_descent()
Expand Down Expand Up @@ -1937,7 +1941,7 @@ def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=3,

return rank_low_bd, two_selmer_rank, pts

two_descent_simon = simon_two_descent
two_descent_simon = simon_two_descent # deprecated in #35621

def three_selmer_rank(self, algorithm='UseSUnits'):
r"""
Expand Down Expand Up @@ -2388,21 +2392,9 @@ def _compute_gens(self, proof,
TESTS::
sage: P = E.lift_x(611429153205013185025/9492121848205441)
sage: set(E.gens(use_database=False, algorithm='pari',pari_effort=4)) <= set([P+T for T
....: in E.torsion_points()] + [-P+T for T in E.torsion_points()])
True
sage: E = EllipticCurve([-157^2,0])
sage: E.gens(use_database=False, algorithm='pari')
Traceback (most recent call last):
...
RuntimeError: generators could not be determined. So far we found []. Hint: increase pari_effort.
sage: ge = E.gens(use_database=False, algorithm='pari',pari_effort=10)
sage: ge #random
[(-166136231668185267540804/2825630694251145858025 : 167661624456834335404812111469782006/150201095200135518108761470235125 : 1)]
sage: P = E.lift_x(-166136231668185267540804/2825630694251145858025)
sage: set(E.gens(use_database=False, algorithm='pari',pari_effort=4)) <= set([P+T for T
....: in E.torsion_points()] + [-P+T for T in E.torsion_points()])
sage: ge = set(E.gens(use_database=False, algorithm='pari',pari_effort=4))
sage: ge <= set([P+T for T in E.torsion_points()]
....: + [-P+T for T in E.torsion_points()])
True
"""
# If the optional extended database is installed and an
Expand Down Expand Up @@ -6209,10 +6201,11 @@ def point_preprocessing(free,tor):
e1,e2,e3 = ei
if r >= 1: #preprocessing of mw_base only necessary if rank > 0
mw_base = point_preprocessing(mw_base, tors_points)
#at most one point in E^{egg}
# at most one point in E^{egg}

elif disc < 0: # one real component => 1 root in RR (=: e3),
# 2 roots in C (e1,e2)
elif disc < 0:
# one real component => 1 root in RR (=: e3),
# 2 roots in C (e1,e2)
roots = pol.roots(C,multiplicities=False)
e3 = pol.roots(R,multiplicities=False)[0]
roots.remove(e3)
Expand Down Expand Up @@ -6695,16 +6688,16 @@ def test_with_T(R):
for T in tors_points:
test(R+T)

# For small rank and small H_q perform simple search
# For small rank and small H_q perform simple search
if r == 1 and N <= 10:
for P in multiples(mw_base[0],N+1):
test_with_T(P)
return xs

# explicit computation and testing linear combinations
# ni loops through all tuples (n_1,...,n_r) with |n_i| <= N
# stops when (0,0,...,0) is reached because after that, only inverse points of
# previously tested points would be tested
# explicit computation and testing linear combinations
# ni loops through all tuples (n_1,...,n_r) with |n_i| <= N
# stops when (0,0,...,0) is reached because after that, only inverse points of
# previously tested points would be tested

E0 = E(0)
ni = [-N for i in range(r)]
Expand Down Expand Up @@ -6857,13 +6850,14 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
prec *= 2
RR = RealField(prec)
ei = pol.roots(RR,multiplicities=False)
e1,e2,e3 = ei
elif disc < 0: # one real component => 1 root in RR (=: e3),
# 2 roots in C (e1,e2)
e1, e2, e3 = ei
elif disc < 0:
# one real component => 1 root in RR (=: e3),
# 2 roots in C (e1,e2)
roots = pol.roots(C,multiplicities=False)
e3 = pol.roots(R,multiplicities=False)[0]
roots.remove(e3)
e1,e2 = roots
e1, e2 = roots

len_tors = len(tors_points)
n = r + 1
Expand Down Expand Up @@ -6896,10 +6890,10 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
if verbose:
print('k1,k2,k3,k4', k1, k2, k3, k4)
sys.stdout.flush()
#H_q -> [PZGH]:N_0 (due to consistency to integral_points())
# H_q -> [PZGH]:N_0 (due to consistency to integral_points())
H_q = R(((k1/2+k2)/lamda).sqrt())

#computation of logs
# computation of logs
mw_base_log = [(pts.elliptic_logarithm().abs())*(len_tors/w1) for pts in mw_base]
mw_base_p_log = []
beta = []
Expand All @@ -6909,7 +6903,7 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
Np = E.Np(p)
cp = E.tamagawa_exponent(p)
mp_temp = Z(len_tors).lcm(cp*Np)
mp.append(mp_temp) #only necessary because of verbose below
mp.append(mp_temp) # only necessary because of verbose below
p_prec = 30+E.discriminant().valuation(p)
p_prec_ok = False
while not p_prec_ok:
Expand All @@ -6920,7 +6914,7 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
p_prec_ok = True
except ValueError:
p_prec *= 2
#reorder mw_base_p: last value has minimal valuation at p
# reorder mw_base_p: last value has minimal valuation at p
mw_base_p_log_val = [mw_base_p_log[tmp][i].valuation() for i in range(r)]
if verbose:
print("mw_base_p_log_val = ",mw_base_p_log_val)
Expand All @@ -6930,7 +6924,7 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
print("min_psi = ", min_psi)
mw_base_p_log[tmp].remove(min_psi)
mw_base_p_log[tmp].append(min_psi)
#beta needed for reduction at p later on
# beta needed for reduction at p later on
try:
beta.append([-mw_base_p_log[tmp][j]/min_psi for j in range(r)])
except ValueError:
Expand All @@ -6945,7 +6939,8 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
print('mw_base_p_log', mw_base_p_log)
sys.stdout.flush()

#constants in reduction (not needed to be computed every reduction step)
# constants in reduction
# (not needed to be computed every reduction step)
k5 = R((2*len_tors)/(3*w1))
k6 = R((k2/len_S).exp())
k7 = R(lamda/len_S)
Expand All @@ -6956,20 +6951,20 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):

break_cond = 0
M = MatrixSpace(Z,n)
#Reduction of initial bound
# Reduction of initial bound
if verbose:
print('initial bound', H_q)
sys.stdout.flush()

while break_cond < 0.9:
#reduction at infinity
# reduction at infinity
bound_list = []
c = R((H_q**n)*100)
m = copy(M.identity_matrix())
for i in range(r):
m[i, r] = R(c*mw_base_log[i]).round()
m[r,r] = max(Z(1), R(c*w1).round())
#LLL - implemented in sage - operates on rows not on columns
# LLL - implemented in sage - operates on rows not on columns
m_LLL = m.LLL()
m_gram = m_LLL.gram_schmidt()[0]
b1_norm = R(m_LLL.row(0).norm())
Expand Down
103 changes: 43 additions & 60 deletions src/sage/schemes/elliptic_curves/gp_simon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,33 +16,20 @@
#
# https://www.gnu.org/licenses/
# ****************************************************************************
from pathlib import Path

from sage.structure.parent_gens import localvars
from cypari2.handle_error import PariError

from sage.interfaces.gp import Gp
from sage.misc.sage_eval import sage_eval
from sage.env import SAGE_EXTCODE
from sage.libs.pari import pari
from sage.misc.randstate import current_randstate
from sage.misc.superseded import deprecation
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.structure.parent_gens import localvars


gp = None


def init():
"""
Function to initialize the gp process
"""
global gp
if gp is None:
import os
from sage.env import DOT_SAGE
logfile = os.path.join(DOT_SAGE, 'gp-simon.log')
gp = Gp(script_subdirectory='simon', logfile=logfile)
gp.read("ellQ.gp")
gp.read("ell.gp")
gp.read("qfsolve.gp")
gp.read("resultant3.gp")
simon_dir = Path(SAGE_EXTCODE) / 'pari' / 'simon'


def simon_two_descent(E, verbose=0, lim1=None, lim3=None, limtriv=None,
Expand All @@ -59,6 +46,9 @@ def simon_two_descent(E, verbose=0, lim1=None, lim3=None, limtriv=None,
sage: import sage.schemes.elliptic_curves.gp_simon
sage: E = EllipticCurve('389a1')
sage: sage.schemes.elliptic_curves.gp_simon.simon_two_descent(E)
doctest:warning...:
DeprecationWarning: please use the 2-descent algorithm over QQ inside pari
See https://github.com/sagemath/sage/issues/38461 for details.
(2, 2, [(5/4 : 5/8 : 1), (-3/4 : 7/8 : 1)])
TESTS::
Expand Down Expand Up @@ -94,75 +84,68 @@ def simon_two_descent(E, verbose=0, lim1=None, lim3=None, limtriv=None,
sage: E.simon_two_descent() # long time
(0, 2, [])
"""
init()
pari.read(simon_dir / "ellQ.gp")
pari.read(simon_dir / "ell.gp")
pari.read(simon_dir / "qfsolve.gp")
pari.read(simon_dir / "resultant3.gp")

current_randstate().set_seed_gp(gp)
current_randstate().set_seed_pari()

K = E.base_ring()
K_orig = K
E_orig = E

# The following is to correct the bug at #5204: the gp script
# fails when K is a number field whose generator is called 'x'.
# It also deals with relative number fields.
E_orig = E

if K is not QQ:
K = K_orig.absolute_field('a')
y = K.gen()
from_K, to_K = K.structure()
E = E_orig.change_ring(to_K)
known_points = [P.change_ring(to_K) for P in known_points]
# Simon's program requires that this name be y.
with localvars(K.polynomial().parent(), 'y'):
gp.eval("K = bnfinit(%s);" % K.polynomial())
if verbose >= 2:
print("K = bnfinit(%s);" % K.polynomial())
gp.eval("%s = Mod(y,K.pol);" % K.gen())
if verbose >= 2:
print("%s = Mod(y,K.pol);" % K.gen())
# Simon's program requires that this name be y.
K_pari = pari.bnfinit(K.polynomial())
known_points = [P.change_ring(to_K) for P in known_points]
else:
deprecation(38461, "please use the 2-descent algorithm over QQ inside pari")
from_K = lambda x: x
to_K = lambda x: x

# The block below mimics the defaults in Simon's scripts, and needs to be changed
# when these are updated.
# The block below mimics the defaults in Simon's scripts.
# They need to be changed when these are updated.
if K is QQ:
cmd = 'ellQ_ellrank(%s, %s);' % (list(E.ainvs()), [P.__pari__() for P in known_points])
over_QQ = True
if lim1 is None:
lim1 = 5
if lim3 is None:
lim3 = 50
if limtriv is None:
limtriv = 3
else:
cmd = 'bnfellrank(K, %s, %s);' % (list(E.ainvs()), [P.__pari__() for P in known_points])
over_QQ = False
if lim1 is None:
lim1 = 2
if lim3 is None:
lim3 = 4
if limtriv is None:
limtriv = 2

gp('DEBUGLEVEL_ell=%s; LIM1=%s; LIM3=%s; LIMTRIV=%s; MAXPROB=%s; LIMBIGPRIME=%s;' % (
verbose, lim1, lim3, limtriv, maxprob, limbigprime))

if verbose >= 2:
print(cmd)
s = gp.eval('ans=%s;' % cmd)
if s.find(" *** ") != -1:
raise RuntimeError("\n%s\nAn error occurred while running Simon's 2-descent program" % s)
if verbose > 0:
print(s)
v = gp.eval('ans')
if v == 'ans': # then the call to ellQ_ellrank() or bnfellrank() failed
raise RuntimeError("An error occurred while running Simon's 2-descent program")
if verbose >= 2:
print("v = %s" % v)

# pari represents field elements as Mod(poly, defining-poly)
# so this function will return the respective elements of K
def _gp_mod(*args):
return args[0]
ans = sage_eval(v, {'Mod': _gp_mod, 'y': K.gen(0)})
lower = ZZ(ans[0])
upper = ZZ(ans[1])
points = [E_orig([from_K(c) for c in list(P)]) for P in ans[2]]
pari('DEBUGLEVEL_ell=%s; LIM1=%s; LIM3=%s; LIMTRIV=%s; MAXPROB=%s; LIMBIGPRIME=%s;' % (
verbose, lim1, lim3, limtriv, maxprob, limbigprime))

try:
if over_QQ:
ans = pari("ellQ_ellrank")(E, known_points)
else:
ans = pari("bnfellrank")(K_pari, E, known_points)
except PariError as err:
raise RuntimeError("an error occurred while running Simon's 2-descent program") from err

loc = {} if over_QQ else {'y': y}
lower, upper, pts = ans.sage(locals=loc)
lower = ZZ(lower)
upper = ZZ(upper)
points = [E_orig([from_K(c) for c in P]) for P in pts]
points = [P for P in points if P.has_infinite_order()]
return lower, upper, points

0 comments on commit c999704

Please sign in to comment.