diff --git a/sources/EllipticCurveQ.Analysis.cs b/sources/EllipticCurveQ.Analysis.cs index 59b9cce..b8b0e76 100644 --- a/sources/EllipticCurveQ.Analysis.cs +++ b/sources/EllipticCurveQ.Analysis.cs @@ -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. @@ -176,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); @@ -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. @@ -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(); - 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); @@ -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; @@ -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 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; @@ -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(); var heights = new List(); @@ -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); @@ -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 gram) { // Given Gram matrix G of current basis and vector v = ,