diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 2072155fb..3ffe3d0b9 100644 --- a/BREAKING-CHANGES.md +++ b/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, diff --git a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Integer.Definition.cs b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Integer.Definition.cs index e3e4c6ccd..45e7f481f 100644 --- a/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Integer.Definition.cs +++ b/Sources/AngouriMath/Core/Entity/Continuous/Entity.Continuous.Integer.Definition.cs @@ -89,7 +89,16 @@ public static Integer Create(EInteger value) /// /// Detemine whether integer is prime or not. /// - public bool IsPrime => CountDivisors() == 2; + /// + /// By Miller–Rabin and Lucas–Lehmer, as membership in PP is. It was whether the number had two + /// divisors, by trial division over long, so 2^89 - 1 threw OverflowException and the + /// largest prime below 2^63 did not finish. Where primality is not decided -- past + /// 3317044064679887385961981, for a number that is not a Mersenne number and passes every test -- + /// it throws . + /// https://github.com/asc-community/AngouriMath/issues/2341 + /// + public bool IsPrime => Functions.Primes.IsPrime(EInteger) + ?? throw new Core.Exceptions.NotSufficientlySupportedException($"Whether {EInteger} is prime is not decided"); /// /// Deconstructs as record diff --git a/Sources/AngouriMath/Core/Transformations/Matching/MatchedRules.cs b/Sources/AngouriMath/Core/Transformations/Matching/MatchedRules.cs index f6821fe78..f60717747 100644 --- a/Sources/AngouriMath/Core/Transformations/Matching/MatchedRules.cs +++ b/Sources/AngouriMath/Core/Transformations/Matching/MatchedRules.cs @@ -584,7 +584,7 @@ internal static class MatchedRules "eulers-totient-of-a-prime-power", MatchPattern.Node( MatchPattern.Node( - MatchPattern.Any("p", p => p.IsPrime), + MatchPattern.Any("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 diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialFactoring.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialFactoring.cs index bd4cb905b..0f872448f 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialFactoring.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/PolynomialFactoring.cs @@ -425,9 +425,9 @@ private static IEnumerable Divisors(EInteger n) { var extended = new List(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)); } diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs index dfee97c14..db44183e8 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.cs +++ b/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, diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt b/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt index 1cb7e022d..5fc8d061f 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/MathAllMethods.tt +++ b/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, diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Discrete/Evaluation.Discrete.Classes.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Discrete/Evaluation.Discrete.Classes.cs index 4934a9bd7..05ffd31a5 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Discrete/Evaluation.Discrete.Classes.cs +++ b/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 }, diff --git a/Sources/AngouriMath/Functions/InternalAMExtensions.cs b/Sources/AngouriMath/Functions/InternalAMExtensions.cs index 9813b0343..9c4948d54 100644 --- a/Sources/AngouriMath/Functions/InternalAMExtensions.cs +++ b/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()); + /// + /// Euler's totient of , 0 for an that is not positive, or + /// where cannot be factorised (). + /// + /// + 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; + } + /// /// Computes Euler phi function /// /// /// 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; + /// + /// Factorization of integer + /// + 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; + /// How far trial division goes before the rest is asked whether it is prime. + private const int TrialDivisionLimit = 1 << 16; + + /// How many steps Pollard's rho takes, for each of its constants, before it gives up on a split. + private const int RhoSteps = 1 << 18; - return result; - } /// - /// Computes Euler phi function (for long) - /// + /// The prime factorisation of , the primes in increasing order, empty for an + /// below 2, or where a factor is neither proven prime nor split. /// - /// If integer x is non-positive, the result will be 0 - public static long Phi(this long n) + /// + /// This was trial division over long, and for every past 2^63 it threw + /// OverflowException, though the first division finds a factor like 2 or 7 at any size: phi(7^30) + /// threw from Simplify. Its loop ran while the square of the divisor was at most the running totient, + /// not what was left to factor, so 2^62 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 to , bounded by what is left, then + /// , 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 + /// + 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(); + 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(); + 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; } - /// - /// Factorization of integer - /// - public static IEnumerable<(long prime, long power)> Factorize(this EInteger n) + /// A proper divisor of the odd composite , or . + 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; } /// - /// Count of all divisors of an integer + /// The number of divisors of , 0 for an that is not positive, or + /// where it cannot be factorised. /// - /// 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; } + /// + /// Count of all divisors of an integer + /// + /// 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) => diff --git a/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.NumberTheory.cs b/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.NumberTheory.cs index 9a3ef866b..7a250a00d 100644 --- a/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.NumberTheory.cs +++ b/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 }; } diff --git a/Sources/AngouriMath/Numerics/SpecialFunctions.cs b/Sources/AngouriMath/Numerics/SpecialFunctions.cs index 623ceb230..dac4fb95a 100644 --- a/Sources/AngouriMath/Numerics/SpecialFunctions.cs +++ b/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; /// /// Half of , part by part: dividing a complex number by 2 multiplies it diff --git a/Sources/Tests/UnitTests/Common/AFactorisationPastTheLongsTest.cs b/Sources/Tests/UnitTests/Common/AFactorisationPastTheLongsTest.cs new file mode 100644 index 000000000..465bb856a --- /dev/null +++ b/Sources/Tests/UnitTests/Common/AFactorisationPastTheLongsTest.cs @@ -0,0 +1,90 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using System.Linq; +using System.Threading.Tasks; +using AngouriMath.Extensions; +using PeterO.Numbers; +using Xunit; +using static AngouriMath.Entity.Number; + +namespace AngouriMath.Tests.Common +{ + /// + /// The factorisation, the divisor count and the totient trial-divided over long: past 2^63 each threw + /// OverflowException, so phi(7^30) threw out of Simplify; the loop was bounded by the running + /// totient rather than what was left, so 2^62 took 1.5e9 divisions; and past 3037000499 the divisor's square + /// wrapped negative, so the largest prime below 2^63 never finished. IsPrime went through the same loop. + /// https://github.com/asc-community/AngouriMath/issues/2341, + /// https://github.com/asc-community/AngouriMath/issues/2342 + /// + [Trait("Area", "Common")] + public sealed class AFactorisationPastTheLongsTest + { + [Theory] + [InlineData("18446744073709551616", "2^64")] + [InlineData("4611686018427387904", "2^62")] + [InlineData("9223372036854775783", "9223372036854775783^1")] + [InlineData("22539340290692258087863249", "7^30")] + [InlineData("618970019642690137449562111", "618970019642690137449562111^1")] + [InlineData("1000000016000000063", "1000000007^1 * 1000000009^1")] + [InlineData("360", "2^3 * 3^2 * 5^1")] + public void IsFactorised(string number, string expected) + { + var factors = Within(() => string.Join(" * ", MathS.NumberTheory.Factorize(Integer.Create(EInteger.FromString(number))) + .Select(f => $"{f.prime}^{f.power}"))); + Assert.Equal(expected, factors); + } + + [Theory] + [InlineData("18446744073709551616", "65")] + [InlineData("4611686018427387904", "63")] + [InlineData("9223372036854775783", "2")] + [InlineData("360", "24")] + public void HasItsDivisorsCounted(string number, string expected) + => Assert.Equal(Integer.Create(EInteger.FromString(expected)), + Within(() => MathS.NumberTheory.CountDivisors(Integer.Create(EInteger.FromString(number))))); + + [Theory] + [InlineData("phi(2^64)", "9223372036854775808")] + [InlineData("phi(2^62)", "2305843009213693952")] + [InlineData("phi(9223372036854775783)", "9223372036854775782")] + [InlineData("phi(7^30)", "19319434534879078361025642")] + [InlineData("phi(360)", "96")] + public void HasItsTotient(string input, string expected) + { + Assert.Equal(expected.ToEntity(), Within(() => input.ToEntity().Evaled)); + Assert.Equal(expected.ToEntity(), Within(() => input.ToEntity().Simplify())); + } + + [Theory] + [InlineData("2147483647", true)] + [InlineData("2305843009213693951", true)] + [InlineData("618970019642690137449562111", true)] + [InlineData("9223372036854775783", true)] + [InlineData("9223372036854775781", false)] + [InlineData("1", false)] + [InlineData("-7", false)] + public void IsDecidedPrime(string number, bool expected) + => Assert.Equal(expected, Within(() => Integer.Create(EInteger.FromString(number)).IsPrime)); + + [Fact] + public void APowerOfALargePrimeHasItsTotient() + { + Assert.Equal("(2 ^ 89 - 1) ^ (k - 1) * (2 ^ 89 - 2)".ToEntity().Simplify(), + Within(() => "phi((2^89 - 1)^k)".ToEntity().Simplify())); + } + + private static T Within(Func evaluation) + { + var task = Task.Factory.StartNew(evaluation, TaskCreationOptions.LongRunning); + Assert.True(task.Wait(TimeSpan.FromSeconds(30)), "did not finish in thirty seconds"); + return task.Result; + } + } +}