From cb1f8b94be17b5b6f902a2d2d764002c16e17aad Mon Sep 17 00:00:00 2001 From: Valery Asiryan Date: Fri, 19 Dec 2025 18:08:30 +0300 Subject: [PATCH 1/4] Optimize torsion filtering in rank search --- sources/EllipticCurveQ.Analysis.cs | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/sources/EllipticCurveQ.Analysis.cs b/sources/EllipticCurveQ.Analysis.cs index 59b9cce..fa5fa8c 100644 --- a/sources/EllipticCurveQ.Analysis.cs +++ b/sources/EllipticCurveQ.Analysis.cs @@ -556,12 +556,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(); var heights = new List(); @@ -572,7 +568,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); @@ -611,6 +607,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 gram) { // Given Gram matrix G of current basis and vector v = , From 6d8deb50394f8d4b51b17974217fc0bc60d1c084 Mon Sep 17 00:00:00 2001 From: Valery Asiryan Date: Fri, 19 Dec 2025 18:12:45 +0300 Subject: [PATCH 2/4] Cap two-isogeny Selmer enumeration --- sources/EllipticCurveQ.Analysis.cs | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sources/EllipticCurveQ.Analysis.cs b/sources/EllipticCurveQ.Analysis.cs index fa5fa8c..efec076 100644 --- a/sources/EllipticCurveQ.Analysis.cs +++ b/sources/EllipticCurveQ.Analysis.cs @@ -308,6 +308,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; From 839d181ab0dd01adabe14b77eacdcc8c6e97a696 Mon Sep 17 00:00:00 2001 From: Valery Asiryan Date: Fri, 19 Dec 2025 18:15:26 +0300 Subject: [PATCH 3/4] Add size guards for rank bound factorization --- sources/EllipticCurveQ.Analysis.cs | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/sources/EllipticCurveQ.Analysis.cs b/sources/EllipticCurveQ.Analysis.cs index efec076..db55d7e 100644 --- a/sources/EllipticCurveQ.Analysis.cs +++ b/sources/EllipticCurveQ.Analysis.cs @@ -106,7 +106,18 @@ 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; + + // Avoid expensive full factorization for large discriminants: use ω(n) ≤ log2(n). + if (BitLength(twoDelta) > 256) + { + omega = BitLength(twoDelta); + } + else + { + omega = InternalMath.FactorAbs(twoDelta).Count; + } int upper = 2 * omega + 2; // Tighter bound when E has a rational 2-torsion point: 2-isogeny descent. @@ -287,6 +298,8 @@ private static int TwoAdicExponent(BigInteger n) { if (b.IsZero) return null; + if (BitLength(b) > 256) return null; + var primes = InternalMath.FactorAbs(b).Keys.ToList(); primes.Sort((x, y) => x.CompareTo(y)); @@ -297,6 +310,7 @@ private static int TwoAdicExponent(BigInteger n) if (n > 30) return null; BigInteger discPart = 2 * b * (a * a - 4 * b); + if (BitLength(discPart) > 256) return null; var checkPrimes = new List(); foreach (var p in InternalMath.FactorAbs(discPart).Keys) { @@ -356,6 +370,13 @@ 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 bool RealSolvableIsogenyQuartic(BigInteger a, BigInteger b, BigInteger d) { if (d.Sign > 0) return true; From e190b3f052de1d6e6d0c138b2fbe5988c24829b1 Mon Sep 17 00:00:00 2001 From: Valery Asiryan Date: Fri, 19 Dec 2025 18:45:04 +0300 Subject: [PATCH 4/4] Avoid heavy factorization in rank bounds --- sources/EllipticCurveQ.Analysis.cs | 84 +++++++++++++++++++++++++----- 1 file changed, 70 insertions(+), 14 deletions(-) diff --git a/sources/EllipticCurveQ.Analysis.cs b/sources/EllipticCurveQ.Analysis.cs index db55d7e..b8b0e76 100644 --- a/sources/EllipticCurveQ.Analysis.cs +++ b/sources/EllipticCurveQ.Analysis.cs @@ -107,17 +107,7 @@ private BigInteger ComputeConductor() // Coarse bound (very safe but often loose): 2*ω(2*Δ_min) + 2. BigInteger absMinDisc = BigInteger.Abs(MinimalDiscriminant); BigInteger twoDelta = absMinDisc == 0 ? 0 : 2 * absMinDisc; - int omega; - - // Avoid expensive full factorization for large discriminants: use ω(n) ≤ log2(n). - if (BitLength(twoDelta) > 256) - { - omega = BitLength(twoDelta); - } - else - { - omega = InternalMath.FactorAbs(twoDelta).Count; - } + int omega = CountDistinctPrimeFactorsUpperBound(twoDelta, trialLimit: 10000); int upper = 2 * omega + 2; // Tighter bound when E has a rational 2-torsion point: 2-isogeny descent. @@ -187,6 +177,8 @@ private static List RationalRootsOfTwoTorsionCubic(BigInteger b2, B // Return rational roots (in lowest terms). Any rational root has denominator dividing 4. var roots = new List(); + if (BitLength(b6) > 256) return roots; + if (b6.IsZero) roots.Add(BigRational.Zero); @@ -300,7 +292,7 @@ private static int TwoAdicExponent(BigInteger n) if (BitLength(b) > 256) return null; - var primes = InternalMath.FactorAbs(b).Keys.ToList(); + 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. @@ -310,10 +302,11 @@ private static int TwoAdicExponent(BigInteger n) if (n > 30) return null; BigInteger discPart = 2 * b * (a * a - 4 * b); - if (BitLength(discPart) > 256) return null; + if (!TryFactorPrimesByTrialDivision(discPart, trialLimit: 10000, out var discPrimes)) return null; var checkPrimes = new List(); - 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); @@ -377,6 +370,69 @@ private static int BitLength(BigInteger n) 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 primes) + { + primes = new List(); + 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;