Factorise past the longs, and decide IsPrime by Miller-Rabin

BREAKING-CHANGES.md

@@ -472,6 +472,23 @@ arcsine of that is `NaN`. Of `arcsin(k + c)` for `k` from 1 to 40 and `c` in `pi
 | `"arcsin(6 + pi)".ToEntity().Evaled` | `NaN` | `1.5707963267948966 - 2.902976737696685i`, to 15 digits |
 | `"arccos(6 + pi)".ToEntity().Evaled` | `NaN` | `2.902976737696685i`, to 15 digits |
 
+### Factorisation, the divisor count, the totient and `IsPrime` answer past 2^63
+
+**Not silent: what threw or did not return now answers.** The three trial-divided over `long`, and `IsPrime` was
+whether the divisor count was 2: past 2^63 each threw `OverflowException`, `phi(7^30)` out of `Simplify` among them,
+and the largest prime below 2^63 did not finish. They divide in `EInteger` now, as far as 65536, and ask
+Miller–Rabin and Lucas–Lehmer of what is left, splitting a composite rest by Pollard's rho. Where a factor is neither
+proven prime nor split, `Factorize`, `CountDivisors`, `Phi` and `IsPrime` throw `NotSufficientlySupportedException`,
+and `phi(n)` in an expression is left as written. https://github.com/asc-community/AngouriMath/issues/2341,
+https://github.com/asc-community/AngouriMath/issues/2342
+
+| Input | Was (2.6.0) | Now |
+|---|---|---|
+| `MathS.NumberTheory.Factorize(2^64)` | `OverflowException` | `(2, 64)` |
+| `"phi(7^30)".Simplify()` | `OverflowException` | `19319434534879078361025642` |
+| `(2^89 - 1).IsPrime` | `OverflowException` | `true` |
+| `MathS.NumberTheory.Factorize(9223372036854775783)` | does not return | `(9223372036854775783, 1)` |
+
 ### An exact rational is zero, positive or negative as it is
 
 **Silent: 0 divided by a rational smaller than about 1e-199 is 0, where it was NaN.** Whether an exact rational was zero,

Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Integer.Definition.cs

@@ -89,7 +89,16 @@ public static Integer Create(EInteger value)
                 /// <summary>
                 /// Detemine whether integer is prime or not.
                 /// </summary>
-                public bool IsPrime => CountDivisors() == 2;
+                /// <remarks>
+                /// By Miller–Rabin and Lucas–Lehmer, as membership in <c>PP</c> is. It was whether the number had two
+                /// divisors, by trial division over <c>long</c>, so <c>2^89 - 1</c> threw <c>OverflowException</c> and the
+                /// largest prime below 2^63 did not finish. Where primality is not decided -- past
+                /// <c>3317044064679887385961981</c>, for a number that is not a Mersenne number and passes every test --
+                /// it throws <see cref="Core.Exceptions.NotSufficientlySupportedException"/>.
+                /// https://github.com/asc-community/AngouriMath/issues/2341
+                /// </remarks>
+                public bool IsPrime => Functions.Primes.IsPrime(EInteger)
+                    ?? throw new Core.Exceptions.NotSufficientlySupportedException($"Whether {EInteger} is prime is not decided");
 
                 /// <summary>
                 /// Deconstructs as record

Sources/AngouriMath/Core/Transformations/Matching/MatchedRules.cs

@@ -584,7 +584,7 @@ internal static class MatchedRules
                 "eulers-totient-of-a-prime-power",
                 MatchPattern.Node<Phif>(
                     MatchPattern.Node<Powf>(
-                        MatchPattern.Any<Integer>("p", p => p.IsPrime),
+                        MatchPattern.Any<Integer>("p", p => Functions.Primes.IsPrime(p.EInteger) is true),
                         MatchPattern.Any("k"))),
                 // (Integer) on the prime and not on the exponent, and the difference is the
                 // point: p - 1 is arithmetic that folds to a number, k - 1 is a tree because k

Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialFactoring.cs

@@ -425,9 +425,9 @@ private static IEnumerable<EInteger> Divisors(EInteger n)
             {
                 var extended = new List<EInteger>(divisors);
                 var multiplier = EInteger.One;
-                for (long i = 0; i < power; i++)
+                for (var i = EInteger.Zero; i < power; i += 1)
                 {
-                    multiplier = multiplier.Multiply(EInteger.FromInt64(prime));
+                    multiplier = multiplier.Multiply(prime);
                     foreach (var divisor in divisors)
                         extended.Add(divisor.Multiply(multiplier));
                 }

Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs

@@ -306,8 +306,8 @@ public static System.Numerics.Complex Round(System.Numerics.Complex a)
         public static System.Numerics.Complex Phi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
         public static double Phi(double a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
         public static float Phi(float a) => (float)AngouriMath.Numerics.SpecialFunctions.Phi(a);
-        public static long Phi(long a) => a.Phi();
-        public static int Phi(int a) => (int)((long)a).Phi();
+        public static long Phi(long a) => PeterO.Numbers.EInteger.FromInt64(a).Phi().ToInt64Checked();
+        public static int Phi(int a) => PeterO.Numbers.EInteger.FromInt32(a).Phi().ToInt32Checked();
         public static BigInteger Phi(BigInteger a) => PeterO.Numbers.EInteger.FromBytes(a.ToByteArray(), littleEndian: true).Phi().ToBigInteger();
 
         // The special functions, in double precision as the kernels define them, and the factorial,

Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt

@@ -225,8 +225,8 @@ namespace AngouriMath.Core.Compilation.IntoLinq
         public static System.Numerics.Complex Phi(System.Numerics.Complex a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
         public static double Phi(double a) => AngouriMath.Numerics.SpecialFunctions.Phi(a);
         public static float Phi(float a) => (float)AngouriMath.Numerics.SpecialFunctions.Phi(a);
-        public static long Phi(long a) => a.Phi();
-        public static int Phi(int a) => (int)((long)a).Phi();
+        public static long Phi(long a) => PeterO.Numbers.EInteger.FromInt64(a).Phi().ToInt64Checked();
+        public static int Phi(int a) => PeterO.Numbers.EInteger.FromInt32(a).Phi().ToInt32Checked();
         public static BigInteger Phi(BigInteger a) => PeterO.Numbers.EInteger.FromBytes(a.ToByteArray(), littleEndian: true).Phi().ToBigInteger();
 
         // The special functions, in double precision as the kernels define them, and the factorial,

Sources/AngouriMath/Functions/Evaluation/Evaluation.Discrete/Evaluation.Discrete.Classes.cs

@@ -786,7 +786,9 @@ protected override Entity InnerSimplify(bool isExact)
                 => ExpandOnOneArgument(Argument,
                     a => a switch
                     {
-                        Integer integer => integer.Phi(),
+                        // Left as written where the argument cannot be factorised, rather than throwing out of a
+                        // simplification. https://github.com/asc-community/AngouriMath/issues/2342
+                        Integer integer => integer.EInteger.TryPhi() is { } totient ? Integer.Create(totient) : null,
                         Number n => MathS.NaN,
                         _ => null
                     },

Sources/AngouriMath/Functions/InternalAMExtensions.cs

@@ -321,106 +321,151 @@ public static EInteger Lcm(this EInteger bigintFirst, EInteger bigintSecond) =>
         public static EInteger Combinations(this EInteger n, EInteger k) =>
             n.Factorial() / ((n - k).Factorial() * k.Factorial());
 
+        /// <summary>
+        /// Euler's totient of <paramref name="n"/>, 0 for an <paramref name="n"/> that is not positive, or
+        /// <see langword="null"/> where <paramref name="n"/> cannot be factorised (<see cref="TryFactorize"/>).
+        /// <a href="https://en.wikipedia.org/wiki/Euler%27s_totient_function"/>
+        /// </summary>
+        internal static EInteger? TryPhi(this EInteger n)
+        {
+            if (n.Sign <= 0)
+                return EInteger.Zero;
+            if (n.TryFactorize() is not { } factors)
+                return null;
+            var result = EInteger.One;
+            foreach (var (prime, power) in factors)
+                result *= prime.Pow(power - 1) * (prime - 1);
+            return result;
+        }
+
         /// <summary>
         /// Computes Euler phi function
         /// <a href="https://en.wikipedia.org/wiki/Euler%27s_totient_function"/>
         /// </summary>
         /// If integer x is non-positive, the result will be 0
         public static EInteger Phi(this EInteger n)
-        {
-            if (n <= 0)
-                return 0;
+            => n.TryPhi() ?? throw Unfactorised(n);
 
-            var result = n.ToInt64Checked();
-            var original = result;
+        /// <summary>
+        /// Factorization of integer
+        /// </summary>
+        public static IEnumerable<(EInteger prime, EInteger power)> Factorize(this EInteger n)
+            => n.TryFactorize() ?? throw Unfactorised(n);
 
-            for (long i = 2; i * i <= result; i++) {
-                if (original % i == 0)
-                {
-                    while (original % i == 0) original /= i;
-                    result -= result / i;
-                }
-            }
+        private static NotSufficientlySupportedException Unfactorised(EInteger n)
+            => new($"{n} has a factor that is neither proven prime nor split, and is not factorised");
 
-            if (original > 1)
-                result -= result / original;
+        /// <summary>How far trial division goes before the rest is asked whether it is prime.</summary>
+        private const int TrialDivisionLimit = 1 << 16;
+
+        /// <summary>How many steps Pollard's rho takes, for each of its constants, before it gives up on a split.</summary>
+        private const int RhoSteps = 1 << 18;
 
-            return result;
-        }
         /// <summary>
-        /// Computes Euler phi function (for long)
-        /// <a href="https://en.wikipedia.org/wiki/Euler%27s_totient_function"/>
+        /// The prime factorisation of <paramref name="n"/>, the primes in increasing order, empty for an
+        /// <paramref name="n"/> below 2, or <see langword="null"/> where a factor is neither proven prime nor split.
         /// </summary>
-        /// If integer x is non-positive, the result will be 0
-        public static long Phi(this long n)
+        /// <remarks>
+        /// This was trial division over <c>long</c>, and for every <paramref name="n"/> past 2^63 it threw
+        /// <c>OverflowException</c>, though the first division finds a factor like 2 or 7 at any size: <c>phi(7^30)</c>
+        /// threw from <c>Simplify</c>. Its loop ran while the square of the divisor was at most the running totient,
+        /// not what was left to factor, so <c>2^62</c> took 1.5e9 divisions after the first had divided every 2 out;
+        /// and past 3037000499 that square wrapped negative, so the largest prime below 2^63 never finished. Now trial
+        /// division in <see cref="EInteger"/> to <see cref="TrialDivisionLimit"/>, bounded by what is left, then
+        /// <see cref="Primes.IsPrime"/>, Miller–Rabin and Lucas–Lehmer, for the rest, and Pollard's rho to split a
+        /// rest that is composite. https://github.com/asc-community/AngouriMath/issues/2342
+        /// </remarks>
+        internal static List<(EInteger Prime, EInteger Power)>? TryFactorize(this EInteger n)
         {
-            if (n <= 0)
-                return 0;
-
-            var result = n;
-            var original = result;
-
-            for (long i = 2; i * i <= result; i++)
+            var byPrime = new SortedDictionary<EInteger, EInteger>();
+            if (n.CompareTo(EInteger.FromInt32(2)) < 0)
+                return new();
+            var rest = n;
+            void Add(EInteger prime, EInteger power)
+                => byPrime[prime] = (byPrime.TryGetValue(prime, out var already) ? already : EInteger.Zero) + power;
+            for (var divisor = 2; divisor <= TrialDivisionLimit; divisor += divisor == 2 ? 1 : 2)
             {
-                if (original % i == 0)
+                var asInteger = EInteger.FromInt32(divisor);
+                if (asInteger * asInteger > rest)
+                    break;
+                var power = EInteger.Zero;
+                while (rest.Remainder(asInteger).IsZero)
                 {
-                    while (original % i == 0) original /= i;
-                    result -= result / i;
+                    rest /= asInteger;
+                    power += 1;
                 }
+                if (!power.IsZero)
+                    Add(asInteger, power);
             }
-
-            if (original > 1)
-                result -= result / original;
-
-            return result;
+            var pending = new Stack<EInteger>();
+            if (!rest.Equals(EInteger.One))
+                pending.Push(rest);
+            while (pending.Count > 0)
+            {
+                Core.Multithreading.MultithreadingFunctional.ExitIfCancelled();
+                var factor = pending.Pop();
+                // Past the trial division a factor whose square is past it has no divisor it did not try.
+                var limit = EInteger.FromInt32(TrialDivisionLimit);
+                var isPrime = factor.CompareTo(limit * limit) <= 0 ? true : Primes.IsPrime(factor);
+                if (isPrime is true)
+                {
+                    Add(factor, EInteger.One);
+                    continue;
+                }
+                if (isPrime is null || PollardRho(factor) is not { } split)
+                    return null;
+                pending.Push(split);
+                pending.Push(factor / split);
+            }
+            var factors = new List<(EInteger, EInteger)>(byPrime.Count);
+            foreach (var pair in byPrime)
+                factors.Add((pair.Key, pair.Value));
+            return factors;
         }
 
-        /// <summary>
-        /// Factorization of integer
-        /// </summary>
-        public static IEnumerable<(long prime, long power)> Factorize(this EInteger n)
+        /// <summary>A proper divisor of the odd composite <paramref name="n"/>, or <see langword="null"/>.</summary>
+        private static EInteger? PollardRho(EInteger n)
         {
-            var result = n.ToInt64Checked();
-            var original = result;
-
-            for (long i = 2; i * i <= result; i++)
+            for (var constant = 1; constant <= 4; constant++)
             {
-                if (original % i == 0)
+                var c = EInteger.FromInt32(constant);
+                EInteger Step(EInteger x) => (x * x + c).Remainder(n);
+                var (slow, fast, divisor) = (EInteger.FromInt32(2), EInteger.FromInt32(2), EInteger.One);
+                for (var step = 0; step < RhoSteps && divisor.Equals(EInteger.One); step++)
                 {
-                    long power = 0;
-                    while (original % i == 0)
-                    {
-                        original /= i;
-                        power++;
-                    }
-
-                    result -= result / i;
-                    yield return (i, power);
+                    slow = Step(slow);
+                    fast = Step(Step(fast));
+                    divisor = (slow - fast).Abs().Gcd(n);
                 }
+                if (!divisor.Equals(EInteger.One) && !divisor.Equals(n))
+                    return divisor;
             }
-
-            if (original > 1)
-                yield return (original, 1L);
+            return null;
         }
 
         /// <summary>
-        /// Count of all divisors of an integer
+        /// The number of divisors of <paramref name="n"/>, 0 for an <paramref name="n"/> that is not positive, or
+        /// <see langword="null"/> where it cannot be factorised.
         /// </summary>
-        /// If integer x is non-positive, the result will be 0
-        public static EInteger CountDivisors(this EInteger n)
+        internal static EInteger? TryCountDivisors(this EInteger n)
         {
-            if (n <= 0)
-                return 0;
-
-            EInteger result = 1;
-            foreach ((var prime, var power) in Factorize(n))
-            {
+            if (n.Sign <= 0)
+                return EInteger.Zero;
+            if (n.TryFactorize() is not { } factors)
+                return null;
+            var result = EInteger.One;
+            foreach (var (_, power) in factors)
                 result *= power + 1;
-            }
-
             return result;
         }
 
+        /// <summary>
+        /// Count of all divisors of an integer
+        /// </summary>
+        /// If integer x is non-positive, the result will be 0
+        public static EInteger CountDivisors(this EInteger n)
+            => n.TryCountDivisors() ?? throw Unfactorised(n);
+
         public static bool GreaterThan(this EDecimal bigDecimalOne, EDecimal bigDecimalTwo) =>
             bigDecimalOne.CompareTo(bigDecimalTwo) > 0;
         public static bool GreaterThanOrEquals(this EDecimal bigDecimalOne, EDecimal bigDecimalTwo) =>

Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.NumberTheory.cs

@@ -14,7 +14,7 @@ internal static partial class Patterns
         [AddressableRules]
         internal static Entity PhiFunctionRules(Entity x) => x switch
         {
-            Phif(Powf(Integer prime, var variable)) when prime.IsPrime => new Powf(prime, variable - 1) * (prime - 1), 
+            Phif(Powf(Integer prime, var variable)) when Primes.IsPrime(prime.EInteger) is true => new Powf(prime, variable - 1) * (prime - 1), 
             _ => x
         };
     }

Sources/AngouriMath/Numerics/SpecialFunctions.cs

@@ -503,7 +503,8 @@ private static NumericsComplex LogSine(NumericsComplex w)
         public static NumericsComplex Phi(NumericsComplex z) => z.Imaginary == 0 ? Phi(z.Real) : NaN;
 
         public static double Phi(double x)
-            => x == Math.Floor(x) && x >= long.MinValue && x < long.MaxValue ? ((long)x).Phi() : double.NaN;
+            => x == Math.Floor(x) && x >= long.MinValue && x < long.MaxValue && PeterO.Numbers.EInteger.FromInt64((long)x).TryPhi() is { } totient
+                ? PeterO.Numbers.EDecimal.FromEInteger(totient).ToDouble() : double.NaN;
 
         /// <summary>
         /// Half of <paramref name="value"/>, part by part: dividing a complex number by 2 multiplies it