;
end function;
procedure sortin(~seq, val)
// seq is sorted decreasingly
for j := 1 to #seq do
if val gt seq[j] then
Insert(~seq, j, val);
Prune(~seq);
return;
end if;
end for;
end procedure;
procedure sortinreplace(~seq, val, old)
// seq is sorted decreasingly
for j := 1 to #seq do
if val gt seq[j] then
Insert(~seq, j, val);
for k := j+1 to #seq do
if old eq seq[k] then
Remove(~seq, k);
return;
end if;
end for;
Prune(~seq);
return;
end if;
end for;
end procedure;
// update GI structure with information in order to reach target eps
procedure computeGI(~GI, eps, ~short, testfun)
// we assume target is not met with input GI
r := GI`rank;
pl := GI`pl;
maxp := pl[#pl];
J := GI`J;
C := Curve(J);
f := HyperellipticPolynomials(C);
seenp := {entry`p : entry in GI`info};
short := 0; // no "shortcut" proof with a single prime known (yet)
t := 1.0;
n := GI`last;
repeat
// find next prime to be used
deepflag := false;
flag := false;
repeat
n +:= 1;
pp, p, k := IsPrimePower(n);
if pp then
if k gt 1 then
flag := p in seenp and p le maxp and testfun(p, k);
if flag then deepflag := true; end if;
else
if not p in GI`badp then
Jp := BaseChange(J, GF(p));
fact, rem := TrialDivision(#Jp, maxp);
flag := IsEmpty(rem) and forall{a : a in fact | a[1] le maxp};
elif p gt 2 and p le maxp
and ((fl1 and fl2) where fl1, fl2 := IsRegularAtp(C, p))
then
fact, rem := TrialDivision(BadOrder(f, p), maxp);
flag := IsEmpty(rem) and forall{a : a in fact | a[1] le maxp};
else
flag := p lt 10;
end if;
end if;
end if;
until flag;
// compute info
if deepflag then
pos := Position([entry`p : entry in GI`info], p);
assert pos gt 0;
// increase level
Isold := GI`info[pos]`invs;
DeepenLevel(~GI`info[pos], ~GI);
entry := GI`info[pos];
if IsEmpty(entry`imC) then
vprintf GroupInfo, 1: " Contradiction at p = %o, level = %o!\n",
p, entry`deep`level;
short := p^entry`deep`level;
return;
end if;
Is := entry`invs;
// sort into maxfactors
for j := 1 to #pl do
for i := 1 to #Isold do
sortinreplace(~GI`maxfactors[j], Is[i,j], Isold[i,j]);
end for;
for i := #Isold + 1 to #Is do
sortin(~GI`maxfactors[j], Is[i,j]);
end for;
end for;
GI`last := n;
else
deg3r := Dred(f, GI`deg3, p);
if p in GI`badp then
if IsOdd(p) and (fl2 where _, fl2 := IsRegularAtp(C, p)) then
entry := InfoAtBadRegularp(C, GI`bas, deg3r, pl, p);
else
if not assigned GI`hp then GI`hp := HeightPairingMatrix(GI`bas); end if;
entry := InfoAtBadp(C, GI`bas, GI`deg3, pl, p, GI`hp);
end if;
else
entry := InfoAtp(Jp, GI`bas, deg3r, pl, p);
end if;
Append(~GI`info, entry);
if IsEmpty(entry`imC) then
vprintf GroupInfo, 1: " Contradiction at p = %o!\n", p;
short := p;
return;
end if;
Is := entry`invs;
// sort into maxfactors
for j := 1 to #pl do
for I in Is do
sortin(~GI`maxfactors[j], I[j]);
end for;
end for;
Include(~seenp, p);
GI`last := p;
end if;
if IsVerbose("GroupInfo", 3) then
printf " maxfactors:\n";
for k := 1 to morefactors+1 do
printf " ";
printfact([GI`maxfactors[j, r+k] : j in [1..#pl]], pl);
printf "\n";
end for;
end if;
// compute current target value
expsize(~GI, GI`bestfactor, ~t); // new value for old bestfactor
bf0 := GI`bestfactor;
vprintf User3: "<%o", n;
for k := 0 to morefactors do
// look at B given by invariants of product group
bf := [GI`maxfactors[j, r+1+k] : j in [1..#pl]];
expsize(~GI, bf, ~t1);
vprintf User3: ", %o", ChangePrecision(t1,4);
if k eq 0 then t0 := t1; end if; // save for later
if t1 lt t then
t := t1;
bf0 := bf;
end if;
end for;
vprintf User3: ">,\n";
GI`bestfactor := bf0;
if IsVerbose("GroupInfo", 2) then
printf " bestfactor: ";
printfact(GI`bestfactor, pl);
printf "\n";
end if;
vprintf GroupInfo, 1: " ==> New target value is %o [%o]\n",
ChangePrecision(t, 4), ChangePrecision(t0, 6);
until t le eps;
end procedure;
// findQSeq computes a `good' sequence of primes
// GI: group info structure
// eps: upper bound for expected size of set of candidates
procedure findQSeq(~GI, eps, ~result, ~Bvec)
pl := GI`pl;
max := GI`bestfactor;
// optimze according to minimal expected size at B_i
// if B_0 = 1, B_{i+1} = B_i*p_i
// B in factored form
agendas := [[<1, [0 : i in [1..#pl]], 1.0>]]; // entries
agendav := [1.0];
reached := {[0 : i in [1..#pl]]}; // lists the values of B that have been seen so far
while not IsEmpty(agendav) do
v, pos := Min(agendav); // try to extend cheapest path so far
l := agendas[pos];
vprintf MWSieve, 2: " findQSeq: l = %o,\n v = %o\n",
[l[i,1] : i in [2..#l]], v;
Remove(~agendas, pos);
Remove(~agendav, pos);
B := l[#l,2];
s := l[#l,3];
if s lt eps then
// reached success criterion
vprintf MWSieve, 1: " Multipliers: %o\n", [l[i,1] : i in [2..#l]];
vprintf MWSieve, 1: " Predicted sizes: %o\n", [Round(l[i,3]) : i in [2..#l]];
result := [ : i in [2..#l]];
Bvec := B;
return;
end if;
// find possible extensions
for j := 1 to #pl do
if B[j] lt max[j] then
Bp := B;
Bp[j] +:= 1;
if Bp notin reached then
expsize(~GI, Bp, ~size);
new := ;
// and update agenda
Append(~agendas, Append(l, new));
Append(~agendav, size);
Include(~reached, Bp);
end if;
end if;
end for;
end while;
result := [];
Bvec := [];
end procedure;
function convertToOldGI(GI, Bvec)
resGI := [];
pl := GI`pl;
for i := 1 to #GI`info do
relsize(~GI`info[i], pl, Bvec, ~s);
if s lt 1.0 then
entry := GI`info[i];
Append(~resGI, [* entry`p, entry`imbas,
ChangeUniverse(entry`imC, Universe(entry`imbas)) *]);
end if;
end for;
return resGI;
end function;
// =========================================================================
// Taken over from previous version:
// Given a set of tuples mod B, return all tuples mod B*p lifting them
// and compatible with GI. Return "primes" used in the process as second
// value.
function LiftInformation(S, GI, B, p)
// S = {[a1,...,ar], ...}
// new version (2008-04-08): split lifting into p-steps
// preparation: find relevant entries in GI
vprintf MWSieve, 2: "LiftInformation: B = %o, p = %o, #S = %o\n", B, p, #S;
r := #GI[1,2];
vpB := Valuation(B, p);
Bp := B*p;
// Select relevant entries from GI
tests := [e : e in GI | Valuation(Exponent(Universe(e[2])), p) gt vpB];
// Reduce to the information needed and remove entries that give no
// restrictions
tests := [[* e[1], [q(b) : b in e[2]], imC *]
: e in tests
| #imC lt #G1
where imC := {q(c) : c in e[3]}
where G1, q := quo< G | [Bp*g : g in Generators(G)]>
where G := Universe(e[2])];
// Add the homomorphism and their kernels to the data
MW := FreeAbelianGroup(r);
tests := [Append(Append(e, h), Kernel(h))
where h := hom Universe(e[2]) | e[2]>
: e in tests];
// Find a filtration
// B*MW = L_0 \supset L_1 \supset ... \supset L_m = B*p*MW
// such that the expected sizes of the sets of candidates
// S_j = {s in MW/L_j : p_ij(h_i(s)) in p_ij(S_i) for all i}
// where h_i : MW --> G_i and p_ij : G_i -->> G_i/h_i(L_j),
// are as small as possible.
// Do it in the greedy way.
// The current subgroup
BMW := sub;
// The target subgroup
BpMW := sub;
// The sequence of subgroups, from BMW down to BpMW (initialize)
steps := [BMW];
Lcurr := BMW;
testsleft := tests;
vprintf MWSieve, 2: " starting computation of steps...\n";
while not IsEmpty(testsleft) do
vprintf MWSieve, 3: " Size of testleft: %o\n", #testsleft;
Lset := {@ Parent(BMW) | @}; // indexed set of subgroups considered
Lval := []; // parallel sequence of "values" of these subgroups
for e in testsleft do
Lnew := Lcurr meet e[5]; // e[5] = ker(h : MW -> G_e)
G := Universe(e[2]);
h := e[4];
G1, q := quo;
pd := #G/#G1;
val := 1.0*#e[3]/#{q(c) : c in e[3]}/pd;
vprintf MWSieve, 3: " q = %o, dim = %o, val = %o -- ",
e[1], Valuation(pd, p), ChangePrecision(pd*val, 4);
// try to find Lnew in Lset
pos := Position(Lset, Lnew);
if pos eq 0 then
// not yet in: add it and initialize the value
vprintf MWSieve, 3: "new subgroup\n";
Include(~Lset, Lnew);
Append(~Lval, pd*val);
else
// already there: update the value
Lval[pos] *:= val;
vprintf MWSieve, 3: "known subgroup, new value = %o\n",
ChangePrecision(Lval[pos], 4);
end if;
end for; // e in testsleft
// Find best new subgroup
m, pos := Min(Lval);
vprintf MWSieve, 2: " expected relative growth: %o, dimension: %o\n",
ChangePrecision(1.0*m, 5),
Valuation(#quo, p);
// and put it into the list
Append(~steps, Lset[pos]);
Lcurr := Lset[pos];
// Remove entries that no longer can give information
testsleft := [e : e in tests | #e[4](Lcurr) gt 1];
end while;
if Lcurr ne BpMW then
vprintf MWSieve, 2: " 'Information dimension' is only %o\n",
r - Valuation(#quo, p);
end if;
vprintf MWSieve, 2: " starting the lifting...\n";
vprintf MWSieve, 2: " intermediate sizes:";
maxsize := 0;
// now lift S successively
primes := {};
Q1, q1 := quo;
S1 := {q1(MW!s) : s in S};
for j := 2 to #steps do
L := steps[j];
QL, qL := quo;
// first determine relevant elements of tests
tests1 := [[* e[1], [q(b) : b in e[2]], {q(c) : c in e[3]},
hom QG | [q(h(g @@ qL)) : g in OrderedGenerators(QL)]>
*]
where QG, q := quo
: e in tests
| GL ne h(steps[j-1])
where GL := h(L)
where h := hom G | e[2]>
where G := Universe(e[2])];
// now lift
hL := hom Q1 | [q1(l @@ qL) : l in OrderedGenerators(QL)]>;
KhL := Set(Kernel(hL));
cands := {};
// for each possible lift of an element of S1 to QL,
// check if it "survives" all the conditions in tests1
for a in S1 do
lifta := qL(a @@ q1); // some lift of a to QL
for k in KhL do
// all lifts are obtained by adding elements from the kernel
c := lifta + k;
if forall(e){e : e in tests1 | e[4](c) in e[3]} then
// we have to keep this candidate
Include(~cands, c);
else
// we can discard it.
// Note in primes that we have used this entry.
Include(~primes, e[1]);
end if;
end for; // k
end for; // a
vprintf MWSieve, 2: " %o", #cands;
maxsize := Max(maxsize, #cands);
S1 := cands;
Q1 := QL;
q1 := qL;
end for; // j
// lift all the way and convert back to sequences
if Lcurr ne BpMW then
last, qlast := quo;
hlast := hom Q1 | [q1(l @@ qlast) : l in OrderedGenerators(last)]>;
Klast := Set(Kernel(hlast));
cands := {(a @@ hlast) + k : k in Klast, a in S1};
maxsize := Max(maxsize, #cands);
vprintf MWSieve, 2: " %o", #cands;
end if;
S := {Eltseq(c) : c in cands};
if IsVerbose("MWSieve", 2) then
if not IsEmpty(primes) then
printf "\n entries used: %o", primes;
end if;
printf "\n size of conditions mod %o is %o\n", Bp, #S;
end if;
return S, primes, maxsize;
end function;
function BMCompute(GI, run)
if not IsEmpty(run) and Type(run[1]) eq RngIntElt then
run := [ : r in run];
end if;
maxsize := 1;
// initialize
B := 1;
r := #GI[1,2];
cands := {[0 : i in [1..r]]};
primes := {};
s := #run;
for j := 1 to #run do
p := run[j,1];
cands, pr, ms := LiftInformation(cands, GI, B, p);
primes join:= pr;
maxsize := Max(maxsize, ms);
B *:= p;
if IsVerbose("MWSieve", 2) then
if run[j,2] lt 0 then
printf "BMCompute: B = %o, size = %o\n", B, #cands;
else
printf "BMCompute: B = %o, size = %o (predicted: %o)\n",
B, #cands, Round(run[j,2]);
end if;
end if;
if IsEmpty(cands) then
if IsVerbose("MWSieve", 1) then
printf "SUCCESS with B = %o = %o\n", B, Factorization(B);
printf "Successive prime factors: %o\n", [run[i,1] : i in [1..j]];
printf "Primes used: %o\n\n", primes;
end if;
s := j;
break;
end if;
end for;
return [Integers() | run[i,1] : i in [1..s]], cands, primes, maxsize;
end function;
// =========================================================================
procedure CheckPoints1(gens, ls, deg3, ~hp, ~info, ~flag, ~point)
// Checks if ls[1]*gens[1] + ... gives
// a point on J that is of the form pt + can - deg3.
// Returns true and the point pt on the curve or false.
// hp is the height pairing matrix of gens
// info[1] = [..., p, ...]
// info[2] = [* ..., f mod p, ... *]
// info[3] = [* ..., [gens mod p], ...*]
// info[4] = [* ..., deg3 mod p, ... *]
J := Universe(gens);
prec := 0.0;
p := 101;
if not assigned info then
info := <[Integers()|], [* *], [* *], [* *]>;
end if;
if not assigned hp then
hp := HeightPairingMatrix(gens);
end if;
// a rough height bound
ht1 := (v*hp, v) + 5.0 where v := Vector(RealField(), ls);
f := PolynomialRing(Integers())!HyperellipticPolynomials(Curve(J));
disc := Discriminant(f);
while prec lt ht1 do
if not IsDivisibleBy(disc, p) then
pos := Position(info[1], p);
if pos eq 0 then
Jp := BaseChange(J, GF(p));
gensp := [Jp!g : g in gens];
d3p := Dred(f, deg3, p);
fp := PolynomialRing(GF(p))!f;
Append(~info[1], p);
Append(~info[2], fp);
Append(~info[3], gensp);
Append(~info[4], d3p);
else
fp := info[2][pos];
gensp := info[3][pos];
d3p := info[4][pos];
end if;
ptJ := TDK(fp, d3p, ls, gensp);
if ptJ[4] ne 0 then
// not on curve
flag := false;
return;
end if;
prec +:= Log(p);
end if;
p := NextPrime(p);
end while;
// Now there should be a point; find it!
pt := TDK(f, deg3, ls, gens);
if pt[4] ne 0 then
// not on curve
flag := false;
else
C := Curve(J);
flag := true;
if pt[3] eq 0 then
point := Rep(PointsAtInfinity(C));
else
point := Rep(Points(C, -pt[2]/pt[3]));
end if;
end if;
return;
end procedure;
// =========================================================================
// Main function
// Do a Mordell-Weil sieve computation in order to show that the curve of J
// (must be of genus 2, defined over Q and of the form y^2 = f(x)) has no
// rational points. bas is a sequence of generators of the Mordell-Weil group,
// deg3 specifies an embedding of the curve into J:
// deg3 is a polynomial b(x) such that f(x) - b(x)^2 has a factor of odd degree.
// SmoothBound influences the selection of primes.
//
// Returns:
// 1. true or false, according to whether there is a rational point or not
// 2. a point on the curve (if true) or _ (if false)
// 3. the "q sequence" of successive prime factors of N
// 4. the set of prime (powers) used for local information
// 5. timing and space information
intrinsic HasPointMWSieve(J::JacHyp, bas::[JacHypPt], deg3::RngUPolElt
: SmoothBound := 200, testfun := func,
eps := 0.01, eps1 := 0.1, CheckSmall := true)
-> BoolElt, PtHyp, SeqEnum, SetEnum, SeqEnum
{Do a Mordell-Weil sieve computation in order to find out if the curve of J
(must be of genus 2, defined over Q and of the form y^2 = f(x)) has any
rational points.
bas is a sequence of generators of the Mordell-Weil group,
deg3 is a polynomial b(x) such that f(x) - b(x)^2 has a factor of odd degree.
SmoothBound, testfun, eps and eps1 are technical parameters. If CheckSmall
is true, the function searches for small points on the curve first.
The first return value is true or false. If true, the second reeturn value
is a rational point on the curve. Further return values have a more technical
meaning.}
t0 := Cputime();
t1 := 0.0;
t2 := 0.0;
t3 := 0.0;
C := Curve(J);
require Genus(C) eq 2: "J must be a genus 2 Jacobian.";
f, h := HyperellipticPolynomials(C);
require h eq 0: "Curve must be of the form y^2 = f(x).";
fd3 := Factorization(f - deg3^2);
degs := [Degree(a[1]) : a in fd3];
if degs[1] eq 1 then
// found a rational point
return true, Rep(Points(C, Roots(fd3[1,1])[1,1])), _, _, _;
end if;
if CheckSmall then
// check for small points
pts := Points(C : Bound := 1000);
if not IsEmpty(pts) then
return true, Rep(pts), _, _, _;
end if;
end if;
pos := Position(degs, 3);
require pos gt 0: "deg3 must give rise to a rational divisor of degree 3.";
deg3 := ;
vprintf MWSieve, 1: "MWSieve: f = %o, rank = %o\n", f, #bas;
// construct GI object
GI := initGI(J, bas, deg3, SmoothBound);
deg3 := GI`deg3;;
repeat
tt := Cputime();
vprintf MWSieve, 1: " Computing GroupInfo...\n";
computeGI(~GI, eps, ~short, testfun);
t1a := Cputime(tt); t1 +:= t1a;
vprintf MWSieve, 1: " Time: %o s\n", t1a;
if short ne 0 then
vprintf MWSieve, 1: " Contradiction at p = %o\n", short;
vprintf MWSieve, 1: " Total time: %o s\n", Cputime(t0);
return false, _, [], {short}, [Cputime(t0), t1, 0.0, 0.0, 0];
end if;
tt := Cputime();
vprintf MWSieve, 1: " Computing q sequence...\n";
findQSeq(~GI, eps1, ~qseq, ~Bvec);
t2a := Cputime(tt); t2 +:= t2a;
vprintf MWSieve, 1: " Time: %o s\n", t2a;
tt := Cputime();
vprintf MWSieve, 1: " Starting sieve computation...\n";
run1, cands, primes, maxsize := BMCompute(convertToOldGI(GI, Bvec), qseq);
t3a := Cputime(tt); t3 +:= t3a;
vprintf MWSieve, 1: " Size of largest intermediate candidate set: %o\n", maxsize;
vprintf MWSieve, 1: " Time: %o s\n", t3a;
if IsEmpty(cands) then
t := Cputime(t0);
vprintf MWSieve, 1: " Total time: %o s\n", t;
return false, _, run1, primes, [t, t1, t2, t3, maxsize];
end if;
// now try to find a point
B := &*[(GI`pl[i])^Bvec[i] : i in [1..#GI`pl]];
vprintf MWSieve, 1: " No contradiction yet. Checking for points...\n";
for c in cands do
// find shortest representative of candidate
ls := [2*a gt B select a-B else a : a in c];
CheckPoints1(bas, ls, deg3, ~hp, ~info, ~flag, ~point);
if flag then
vprintf MWSieve, 1: " Point %o found!\n", point;
t := Cputime(t0);
vprintf MWSieve, 1: " Total time: %o s\n", t;
return true, point, _, _, [t, t1, t2, t3, maxsize];
end if;
end for;
// reduce eps and eps1 so that we will use more information in the next
// pass through the loop
eps *:= 0.1;
eps1 *:= 0.1;
until false;
end intrinsic;