Gaussian ABC conjecture

GaussianABC.cs

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;
		}
	}
}