Skip to content

Instantly share code, notes, and snippets.

@pjt33

pjt33/GaussianABC.cs

Last active Apr 5, 2018
Embed
What would you like to do?
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<BigInteger> _IntegerPrimes = new List<BigInteger> { 2, 3, 5 };
private static IEnumerable<BigInteger> 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<GaussianInteger> GaussianPrimes => IntegerPrimes.SelectMany(p => GaussianPrimesFor(p));
private static IEnumerable<GaussianInteger> 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<BigInteger> 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;
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.