Skip to content

Use Mazur-bound small-multiple torsion check in rank search #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
asiryan wants to merge 4 commits into
base: codex/verify-and-implement-rank-calculation-algorithm
Choose a base branch
from
Open

Use Mazur-bound small-multiple torsion check in rank search #6

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
cb1f8b9
Optimize torsion filtering in rank search
asiryan Dec 19, 2025
6d8deb5
Cap two-isogeny Selmer enumeration
asiryan Dec 19, 2025
839d181
Add size guards for rank bound factorization
asiryan Dec 19, 2025
e190b3f
Avoid heavy factorization in rank bounds
asiryan Dec 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
110 changes: 100 additions & 10 deletions sources/EllipticCurveQ.Analysis.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,8 @@ private BigInteger ComputeConductor()

// Coarse bound (very safe but often loose): 2*ω(2*Δ_min) + 2.
BigInteger absMinDisc = BigInteger.Abs(MinimalDiscriminant);
int omega = InternalMath.FactorAbs(absMinDisc == 0 ? 0 : 2 * absMinDisc).Count;
BigInteger twoDelta = absMinDisc == 0 ? 0 : 2 * absMinDisc;
int omega = CountDistinctPrimeFactorsUpperBound(twoDelta, trialLimit: 10000);
int upper = 2 * omega + 2;

// Tighter bound when E has a rational 2-torsion point: 2-isogeny descent.
Expand Down Expand Up @@ -176,6 +177,8 @@ private static List<BigRational> RationalRootsOfTwoTorsionCubic(BigInteger b2, B
// Return rational roots (in lowest terms). Any rational root has denominator dividing 4.
var roots = new List<BigRational>();

if (BitLength(b6) > 256) return roots;

if (b6.IsZero)
roots.Add(BigRational.Zero);

Expand Down Expand Up @@ -287,7 +290,9 @@ private static int TwoAdicExponent(BigInteger n)
{
if (b.IsZero) return null;

var primes = InternalMath.FactorAbs(b).Keys.ToList();
if (BitLength(b) > 256) return null;

if (!TryFactorPrimesByTrialDivision(b, trialLimit: 10000, out var primes)) return null;
primes.Sort((x, y) => x.CompareTo(y));

// Selmer elements live in the squareclass group generated by -1 and primes dividing b.
Expand All @@ -297,9 +302,11 @@ private static int TwoAdicExponent(BigInteger n)
if (n > 30) return null;

BigInteger discPart = 2 * b * (a * a - 4 * b);
if (!TryFactorPrimesByTrialDivision(discPart, trialLimit: 10000, out var discPrimes)) return null;
var checkPrimes = new List<int>();
foreach (var p in InternalMath.FactorAbs(discPart).Keys)
for (int i = 0; i < discPrimes.Count; i++)
{
var p = discPrimes[i];
if (p < 2) continue;
if (p > int.MaxValue) return null;
checkPrimes.Add((int)p);
Expand All @@ -308,6 +315,10 @@ private static int TwoAdicExponent(BigInteger n)

int total = 1 << n;

// Guard against exponential blow-ups in the Selmer enumeration.
const int maxSelmerMasks = 1 << 18; // 262,144
if (total > maxSelmerMasks) return null;

// Gaussian elimination basis (bit i corresponds to generator i).
ulong[] basis = new ulong[n];
int rank = 0;
Expand Down Expand Up @@ -352,6 +363,76 @@ private static int TwoAdicExponent(BigInteger n)
return rank;
}

private static int BitLength(BigInteger n)
{
if (n.Sign < 0) n = BigInteger.Abs(n);
if (n.IsZero) return 0;
return n.ToByteArray().Length * 8;
}

private static int CountDistinctPrimeFactorsUpperBound(BigInteger n, int trialLimit)
{
if (n.Sign < 0) n = BigInteger.Abs(n);
if (n <= 1) return 0;

int count = 0;
BigInteger remaining = n;

if ((remaining & 1) == 0)
{
count++;
while ((remaining & 1) == 0) remaining >>= 1;
}

for (int p = 3; p <= trialLimit && remaining > 1; p += 2)
{
if (remaining % p != 0) continue;
count++;
while (remaining % p == 0) remaining /= p;
}

if (remaining > 1)
{
count += Math.Max(1, BitLength(remaining));
}

return count;
}

private static bool TryFactorPrimesByTrialDivision(BigInteger n, int trialLimit, out List<BigInteger> primes)
{
primes = new List<BigInteger>();
if (n.Sign < 0) n = BigInteger.Abs(n);
if (n <= 1) return true;

BigInteger remaining = n;
if ((remaining & 1) == 0)
{
primes.Add(2);
while ((remaining & 1) == 0) remaining >>= 1;
}

for (int p = 3; p <= trialLimit && remaining > 1; p += 2)
{
if (remaining % p != 0) continue;
primes.Add(p);
while (remaining % p == 0) remaining /= p;
}

if (remaining == 1) return true;
if (remaining == 0) return false;

if (remaining > 1)
{
// If the remainder is composite with large factors, skip the expensive path.
if (BitLength(remaining) > 32) return false;
primes.Add(remaining);
return true;
}

return true;
}

private static bool RealSolvableIsogenyQuartic(BigInteger a, BigInteger b, BigInteger d)
{
if (d.Sign > 0) return true;
Expand Down Expand Up @@ -556,12 +637,8 @@ private int RankLowerBoundFromSearch(int searchNumMax, int searchDenMax, int max
// the numerical linear algebra bounded.
//
// Torsion filter:
// Over Q, the torsion subgroup is classified (Mazur). In particular, every torsion point
// has order dividing lcm(1,2,3,4,5,6,7,8,9,10,12,16) = 5040, hence 5040*P = O iff P is torsion.
//
// This avoids computing the full torsion subgroup (your TorsionPoints cache), which can
// dominate runtime when rank bounds are requested repeatedly.
const int torsionKiller = 5040;
// Over Q, possible torsion orders are at most 16 (Mazur). Checking up to 16 avoids
// huge scalar multiplications while still discarding torsion points efficiently.

var basis = new List<EllipticCurvePoint>();
var heights = new List<double>();
Expand All @@ -572,7 +649,7 @@ private int RankLowerBoundFromSearch(int searchNumMax, int searchDenMax, int max
{
if (processed++ >= maxPoints) break;
if (p.IsInfinity) continue;
if (Multiply(p, torsionKiller).IsInfinity) continue;
if (IsTorsionByMazurBound(p)) continue;

// Approximate canonical height.
double hp = CanonicalHeightApprox(p, heightDoublings);
Expand Down Expand Up @@ -611,6 +688,19 @@ private int RankLowerBoundFromSearch(int searchNumMax, int searchDenMax, int max
return basis.Count;
}

private bool IsTorsionByMazurBound(EllipticCurvePoint p)
{
// Any torsion point over Q has order <= 16, so test multiples up to 16.
if (p.IsInfinity) return true;
var q = EllipticCurvePoint.Infinity;
for (int k = 1; k <= 16; k++)
{
q = Add(q, p);
if (q.IsInfinity) return true;
}
return false;
}

private static bool NumericallyIndependent(double hp, double[] v, List<double[]> gram)
{
// Given Gram matrix G of current basis and vector v = <P, basis>,
Expand Down