{{ message }}

Instantly share code, notes, and snippets.

# pjt33/GaussianABC.cs

Last active Apr 5, 2018
Gaussian ABC conjecture
 using System; using System.Collections.Generic; using System.Linq; using System.Text; using BigInteger = System.Numerics.BigInteger; namespace Sandbox { // https://math.stackexchange.com/questions/2642142/complex-number-abc-conjecture // https://stackoverflow.com/questions/2269810/whats-a-nice-method-to-factor-gaussian-integers class GaussianABC { internal static void Main() { var units = new GaussianInteger[] { GaussianInteger.One, GaussianInteger.I, -GaussianInteger.One, -GaussianInteger.I }; foreach (var p in GaussianPrimes) { if (p.Im > 0) continue; Console.WriteLine($"*{p}"); // Crude progress tracking foreach (var q in GaussianPrimes) { if (q == p) break; var ppow = GaussianInteger.One; for (int x = 1; x < 9; x++) { ppow *= p; var qpow = GaussianInteger.One; for (int y = 1; y < 9; y++) { qpow *= q; foreach (var u in units) { var a = ppow; var b = qpow * u; var c = a + b; // Q = ln |c| / \sum_{prime p | abc} ln |p| // But wlog we can rearrange and tweak units to maximise c var numerator = a.NormSquared.Max(b.NormSquared).Max(c.NormSquared).Log(); var denominatorSum = p.NormSquared.Log() + q.NormSquared.Log(); // Factor C var norm = c.NormSquared; var remaining = c; var factorisation = new StringBuilder(); foreach (var prime in IntegerPrimes) { if (norm == 1) break; if (prime * (BigInteger)prime > norm) break; // Save effort on low-quality sums if (denominatorSum >= numerator) { norm = 1; break; } // Sanity check: stick to small-ish cases if (prime > (1 << 20)) { numerator = 0; norm = 1; break; } if (norm % prime == 0) { foreach (var candidateFactor in GaussianPrimesFor(prime)) { if (remaining % candidateFactor == GaussianInteger.Zero) { denominatorSum += candidateFactor.NormSquared.Log(); int exp = 1; remaining /= candidateFactor; while (remaining % candidateFactor == GaussianInteger.Zero) { exp++; remaining /= candidateFactor; } factorisation.Append($"{candidateFactor}^{{{exp}}}"); } } norm /= prime; while (norm % prime == 0) norm /= prime; } } if (norm > 1) { foreach (var candidateFactor in GaussianPrimesFor(norm)) { if (c % candidateFactor == GaussianInteger.Zero) { denominatorSum += candidateFactor.NormSquared.Log(); int exp = 1; remaining /= candidateFactor; while (remaining % candidateFactor == GaussianInteger.Zero) { exp++; remaining /= candidateFactor; } factorisation.Append($"{candidateFactor}^{{{exp}}}"); } } } var Q = numerator / denominatorSum; if (Q > 1.2) { var maxNorm = a.NormSquared.Max(b.NormSquared).Max(c.NormSquared); if (c.NormSquared == maxNorm) Console.WriteLine($"Q={Q:F5}:& {p}^{x} + {u}*{q}^{y} &=& {remaining}*{factorisation}"); else if (a.NormSquared == maxNorm) Console.WriteLine($"Q={Q:F5}:& {remaining}*{factorisation} + {-u}*{q}^{y} &=& {p}^{x}"); else Console.WriteLine($"Q={Q:F5}:& {remaining}*{factorisation} - {p}^{x} &=& {u}*{q}^{y}"); } } } } } } } private static List _IntegerPrimes = new List { 2, 3, 5 }; private static IEnumerable IntegerPrimes { get { // This could be more efficient int len = _IntegerPrimes.Count; for (int i = 0; i < len; i++) yield return _IntegerPrimes[i]; for (var p = _IntegerPrimes[len - 1] + 2; true; p += 2) { // Test division bool isComposite = false; foreach (var q in _IntegerPrimes) { if (p % q == 0) { isComposite = true; break; } if (q * q > p) break; } if (isComposite) continue; yield return p; // Handle multi-iterator, albeit crappily if (p > _IntegerPrimes[_IntegerPrimes.Count - 1]) _IntegerPrimes.Add(p); } } } private static IEnumerable GaussianPrimes => IntegerPrimes.SelectMany(p => GaussianPrimesFor(p)); private static IEnumerable GaussianPrimesFor(BigInteger p) { // Assume p is an integer prime switch ((int)(p % 4)) { case 0: throw new ArgumentOutOfRangeException(nameof(p)); case 2: yield return new GaussianInteger(1, 1); break; case 3: yield return new GaussianInteger(p, 0); break; default: // Find k such that (k-i)(k+i) == 0 (mod p) var k = Range().Where(n => n.ModPow((p - 1) / 2, p) == p - 1).Select(n => n.ModPow((p - 1) / 4, p)).First(); var a = new GaussianInteger(p, 0); var b = new GaussianInteger(k, 1); while (b.Re != 0 || b.Im != 0) { var tmp = a % b; a = b; b = tmp; } yield return new GaussianInteger(a.Re.Abs(), a.Im.Abs()); yield return new GaussianInteger(a.Re.Abs(), -a.Im.Abs()); break; } } private static IEnumerable Range() { BigInteger a = 2; while (true) { yield return a; a++; } } struct GaussianInteger { public static GaussianInteger Zero => new GaussianInteger(0, 0); public static GaussianInteger One => new GaussianInteger(1, 0); public static GaussianInteger I => new GaussianInteger(0, 1); public GaussianInteger(BigInteger re, BigInteger im) { Re = re; Im = im; } public BigInteger Re { get; private set; } public BigInteger Im { get; private set; } public static bool operator ==(GaussianInteger a, GaussianInteger b) => a.Re == b.Re && a.Im == b.Im; public static bool operator !=(GaussianInteger a, GaussianInteger b) => !(a == b); public static GaussianInteger operator +(GaussianInteger a, GaussianInteger b) => new GaussianInteger(a.Re + b.Re, a.Im + b.Im); public static GaussianInteger operator -(GaussianInteger a, GaussianInteger b) => new GaussianInteger(a.Re - b.Re, a.Im - b.Im); public static GaussianInteger operator -(GaussianInteger a) => new GaussianInteger(-a.Re, -a.Im); public static GaussianInteger operator *(GaussianInteger a, GaussianInteger b) => new GaussianInteger(a.Re * b.Re - a.Im * b.Im, a.Re * b.Im + a.Im * b.Re); public static GaussianInteger operator %(GaussianInteger a, GaussianInteger b) => a - a / b * b; public static GaussianInteger operator /(GaussianInteger a, GaussianInteger b) { // (c + di)(c - di) = c^2 + d^2 // Therefore 1 / (c + di) = (c - di) / (c^2 + d^2) // Therefore (a + bi) / (c + di) = (a + bi)(c - di) / (c^2 + d^2) var n = a * b.Conjugate; var d = b.NormSquared; // Need to round each part to nearest return new GaussianInteger((n.Re + n.Re.Sign() * d / 2) / d, (n.Im + n.Im.Sign() * d / 2) / d); } public GaussianInteger Conjugate => new GaussianInteger(Re, -Im); public BigInteger NormSquared => Re * Re + Im * Im; public override int GetHashCode() => (int)(Re + 37 * Im); public override bool Equals(object obj) => obj is GaussianInteger && Equals((GaussianInteger)obj); public bool Equals(GaussianInteger z) => Re == z.Re && Im == z.Im; public override string ToString() { if (Im == 0) return $"{Re}"; if (Re == 0) { if (Im == 1) return$"i"; if (Im == -1) return $"-i"; return$"{Im}i"; } if (Im == 1) return $"({Re}+i)"; if (Im == -1) return$"({Re}-i)"; if (Im > 0) return $"({Re}+{Im}i)"; return$"({Re}{Im}i)"; } } } static class BigIntegerExt { public static BigInteger Max(this BigInteger a, BigInteger b) => a > b ? a : b; public static BigInteger Abs(this BigInteger a) => a > 0 ? a : -a; public static int Sign(this BigInteger a) => a == 0 ? 0 : a > 0 ? 1 : -1; public static double Log(this BigInteger a) => Math.Log((double)a); public static BigInteger ModPow(this BigInteger a, BigInteger pow, BigInteger modulus) { BigInteger accum = 1; BigInteger x = a; BigInteger exp = pow; while (exp > 0 && x != 1) { if ((exp & 1) == 1) accum = accum * x % modulus; x = x * x % modulus; exp >>= 1; } return accum; } } }