From d947f85aaebfea0fcec770786443e9825eaaef52 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 30 Oct 2026 08:11:14 +0000 Subject: [PATCH] Primality is decided by Miller-Rabin and Lucas-Lehmer, and finishes Primes.IsPrime was trial division, decided only for a number that fits a long, and it had three defects: - It hung. The loop was divisor * divisor <= value, and past divisor = 3037000499 the square overflows long and wraps negative. So for a prime between 3037000499^2 and 2^63 the loop ran until the divisor reached the number itself. 9223372036854775783 in PP, the largest prime below 2^63, did not finish on 2.6.0 or on master. - It was slow: about 1.5e9 divisions for a prime near 2^63, and 1.1 s for 2^61 - 1. - It left every number past 2^63 open, the Mersenne prime 2^89 - 1 and the composite 2^67 - 1 alike. @Happypig375 pointed to that on #2221. It is now: - Miller-Rabin to the first thirteen primes as bases. That is exact below psi13 = 3317044064679887385961981, the least strong pseudoprime to all thirteen (Sorenson and Webster, "Strong pseudoprimes to twelve prime bases", Math. Comp. 86 (2017)). The first twelve bases already suffice below 2^64. - Lucas-Lehmer for a Mersenne number 2^p - 1, after checking that p is prime, for p up to 5000. - Past psi13 and up to 4096 bits, a failed round still proves a number composite. A number that passes every round is left open, since passing is not a proof. 9223372036854775783 in PP is True in 81 ms, 2^89 - 1 is prime, and 10^30 + 1 is composite. 10^30 + 57, a prime past psi13, stays open, as does psi13 itself. The compiled PP membership calls the same test for any whole double. It was left open past 2^63. #2221's MembershipInASetOperationTest, and its BREAKING-CHANGES entry, used 2^89 - 1 as a number PP leaves open. Lucas-Lehmer now decides it, so both use 10^30 + 57 instead, with an interval that holds it, [0; 10^31], and the entry's values re-measured on 2.6.0 and on this build. PrimalityTest has 16 decided cases under a 30-second cap, among them the largest prime below 2^63, seven Mersenne numbers from 2^61 - 1 to 2^4423 - 1, the least strong pseudoprime to the first twelve bases, Carmichael numbers and base-2 strong pseudoprimes. It has two cases that stay open, every n below 200000 against a sieve, and 2000 random odd n below 10^12 against trial division. Fixes #2230. Co-Authored-By: Claude Opus 5.5 (1M context) --- BREAKING-CHANGES.md | 28 +++++- Sources/.editorconfig | 3 + .../Core/Entity/Omni/Entity.Omni.Classes.cs | 4 +- .../Functions/NumberTheory/Primes.cs | 81 +++++++++++++--- .../AngouriMath/Numerics/NumberStatements.cs | 9 +- .../Sets/MembershipInASetOperationTest.cs | 29 +++--- .../UnitTests/Core/Sets/PrimalityTest.cs | 96 +++++++++++++++++++ 7 files changed, 210 insertions(+), 40 deletions(-) create mode 100644 Sources/Tests/UnitTests/Core/Sets/PrimalityTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index b05ac8721..c7f83978a 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -121,6 +121,24 @@ carried it once combined: a second derivative of `1/sqrt(x^2 + 13)` stated `not | `"(a provided b > 0) + (c provided b > 0)".ToEntity().InnerSimplified` | `a + c provided b > 0 and b > 0` | `a + c provided b > 0` | | `MathS.Equations("sqrt(y) - x", "y - a").Solve("x", "y")` | `[[sqrt(a), (sqrt(a)) ^ 2], [-sqrt(a), (-sqrt(a)) ^ 2]]` | `[[sqrt(a) provided C(sqrt(a)), (sqrt(a)) ^ 2 provided C(sqrt(a))], [-sqrt(a) provided C(-sqrt(a)), (-sqrt(a)) ^ 2 provided C(-sqrt(a))]]`, `C(v)` being `v = 0 or -1/2 * pi < arg(v) and arg(v) <= 1/2 * pi` | +### Primality is decided by Miller–Rabin and Lucas–Lehmer, and finishes + +**Silent: a membership in `PP` left open past the machine word is decided where a test proves it, and one that did +not finish does.** Primality was trial division, exact only below `2^63`. For a prime past `3037000499^2`, just below +`2^63`, the square of the divisor overflowed and the division never finished, and a prime near `2^63` took about +`1.5e9` divisions. It is now Miller–Rabin to the first thirteen primes as bases, exact below +`3317044064679887385961981` (Sorenson and Webster, "Strong pseudoprimes to twelve prime bases", Math. Comp. 86 +(2017)), and Lucas–Lehmer for a Mersenne number `2^p - 1`. Past those a failed round still proves a number composite, +and a number that passes every round is left open, since passing is not a proof +([#2230](https://github.com/asc-community/AngouriMath/issues/2230)). + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"9223372036854775783 in PP".ToEntity().Evaled` | does not finish; stopped after 60 s | `True`, in 81 ms | +| `"(2^89 - 1) in PP".ToEntity().Evaled` | `618970019642690137449562111 in PP` | `True` | +| `"(10^30 + 1) in PP".ToEntity().Evaled` | `1000000000000000000000000000001 in PP` | `False` | +| `"(10^30 + 57) in PP".ToEntity().Evaled` | `1000000000000000000000000000057 in PP` | `1000000000000000000000000000057 in PP` | + ### A system with a product equation is solved branch by branch **Silent: solutions that were missing are answered.** `EquationSystem.Solve` eliminated `x * y = 0` for `y` as `0`, @@ -603,16 +621,16 @@ around a principal `n`-th root of an `n`-th power, for a whole `n >= 2`, establi **Silent: a membership answered `False` is left open, where it was wrong; and some left open are decided.** Membership in a union, an intersection and a set difference is decided from each side's by Kleene's connectives. A union was -`False` where one side left membership open and the other excluded the element: `2^89 - 1` is a Mersenne prime past -the machine word, `PP` leaves its membership open, and it was outside `PP \/ [0; 1]`. An intersection and a difference +`False` where one side left membership open and the other excluded the element: `10^30 + 57` is a prime past what the +prime test proves, `PP` leaves its membership open, and it was outside `PP \/ [0; 1]`. An intersection and a difference were decided only where both sides were, and are now decided where one side settles it ([#2219](https://github.com/asc-community/AngouriMath/issues/2219)). | Input | Was (2.6.0) | Now | |---|---|---| -| `"(2^89 - 1) in (PP \/ [0; 1])".ToEntity().Evaled` | `False` | `618970019642690137449562111 in PP \/ [0; 1]` | -| `"(2^89 - 1) in ([0; 1] \ PP)".ToEntity().Evaled` | `618970019642690137449562111 in [0; 1] \ PP` | `False` | -| `"(2^89 - 1) in (PP \ [0; 10^30])".ToEntity().Evaled` | `618970019642690137449562111 in PP \ [0; 1000000000000000000000000000000]` | `False` | +| `"(10^30 + 57) in (PP \/ [0; 1])".ToEntity().Evaled` | `False` | `1000000000000000000000000000057 in PP \/ [0; 1]` | +| `"(10^30 + 57) in ([0; 1] \ PP)".ToEntity().Evaled` | `1000000000000000000000000000057 in [0; 1] \ PP` | `False` | +| `"(10^30 + 57) in (PP \ [0; 10^31])".ToEntity().Evaled` | `1000000000000000000000000000057 in PP \ [0; 10000000000000000000000000000000]` | `False` | ### A system is not answered with a copy of one of its rows diff --git a/Sources/.editorconfig b/Sources/.editorconfig index bea86d119..cddc2c8b3 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -478,6 +478,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Core/Sets/MembershipInASetOperationTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n +[Tests/UnitTests/Core/Sets/PrimalityTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + [Tests/UnitTests/Calculus/InverseHyperbolicSecantIntegralTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Core/Entity/Omni/Entity.Omni.Classes.cs b/Sources/AngouriMath/Core/Entity/Omni/Entity.Omni.Classes.cs index d5b79ca72..c9a41d0e2 100644 --- a/Sources/AngouriMath/Core/Entity/Omni/Entity.Omni.Classes.cs +++ b/Sources/AngouriMath/Core/Entity/Omni/Entity.Omni.Classes.cs @@ -761,8 +761,8 @@ public override bool TryContains(Entity entity, out bool contains) } /// - /// Whether the whole number is prime, or where it does - /// not fit a machine word and trial division is not attempted. + /// Whether the whole number is prime, or where it is not + /// decided: see . /// internal static bool? IsDecidablyPrime(Integer whole) => Functions.Primes.IsPrime(whole.EInteger); diff --git a/Sources/AngouriMath/Functions/NumberTheory/Primes.cs b/Sources/AngouriMath/Functions/NumberTheory/Primes.cs index 2aaa7046c..dcafde233 100644 --- a/Sources/AngouriMath/Functions/NumberTheory/Primes.cs +++ b/Sources/AngouriMath/Functions/NumberTheory/Primes.cs @@ -105,29 +105,88 @@ static void AddPrimeI() } } + // The first thirteen primes. Miller–Rabin to them as bases is exact below psi13, the least strong + // pseudoprime to all thirteen: Sorenson and Webster, "Strong pseudoprimes to twelve prime bases", + // Math. Comp. 86 (2017). The first twelve already suffice below 2^64. + [ConstantField] private static readonly int[] witnesses = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41 }; + [ConstantField] private static readonly EInteger psi13 = EInteger.FromString("3317044064679887385961981"); + + // Past these sizes neither test is attempted, and primality is left open: Lucas–Lehmer at 2^5000 - 1 + // takes about 5000 squarings of 5000 bits, and Miller–Rabin at 4096 bits thirteen powers of that size. + private const long LargestMersenneExponent = 5000; + private const long LargestTestedBits = 4096; + /// - /// Whether the whole number is prime, or where it is not decided: - /// trial division, exact for a number that fits a machine word, and nothing claimed past - /// that, so that a membership in PP is left open rather than answered slowly or - /// wrongly. https://github.com/asc-community/AngouriMath/issues/1450 + /// Whether the whole number is prime, or where it is not decided. Miller–Rabin to the + /// first thirteen primes is exact below 3317044064679887385961981; past it a witness still proves a + /// number composite, a Mersenne number 2^p - 1 is decided by Lucas–Lehmer, and a number that passes is + /// left open, since passing is not a proof. This was trial division, exact only below 2^63, about 1.5e9 + /// divisions for a prime near it, and for one past 3037000499^2 never finishing, the square of the divisor + /// overflowing. https://github.com/asc-community/AngouriMath/issues/2230 /// internal static bool? IsPrime(EInteger n) { if (n.CompareTo(EInteger.FromInt32(2)) < 0) return false; - if (!n.CanFitInInt64()) - return null; - var value = n.ToInt64Checked(); - if (value < 4) + foreach (var witness in witnesses) + { + if (n.CompareTo(EInteger.FromInt32(witness)) == 0) + return true; + if (n.Remainder(EInteger.FromInt32(witness)).IsZero) + return false; + } + // No factor up to 41 and below 43^2: prime. + if (n.CompareTo(EInteger.FromInt32(43 * 43)) < 0) return true; - if (value % 2 == 0) + var bits = n.GetUnsignedBitLengthAsInt64(); + // 2^p - 1, whose successor is a power of two. + if (n.Add(EInteger.One).And(n).IsZero && bits <= LargestMersenneExponent) + return IsPrime(EInteger.FromInt64(bits)) is true && LucasLehmer(n, bits); + if (bits > LargestTestedBits) + return null; + if (!IsStrongProbablePrime(n)) return false; - for (long divisor = 3; divisor * divisor <= value; divisor += 2) - if (value % divisor == 0) + return n.CompareTo(psi13) < 0 ? true : null; + } + + // Whether n passes Miller–Rabin to every witness, n odd and past them. + private static bool IsStrongProbablePrime(EInteger n) + { + var less = n.Subtract(EInteger.One); + var odd = less; + var twos = 0; + while (odd.IsEven) + { + odd = odd.ShiftRight(1); + twos++; + } + foreach (var witness in witnesses) + { + var x = EInteger.FromInt32(witness).ModPow(odd, n); + if (x.Equals(EInteger.One) || x.Equals(less)) + continue; + var reached = false; + for (var i = 1; i < twos && !reached; i++) + { + x = x.Multiply(x).Remainder(n); + reached = x.Equals(less); + } + if (!reached) return false; + } return true; } + // 2^p - 1, p an odd prime, is prime exactly when s reaches 0 modulo it after p - 2 steps of s -> s^2 - 2 + // from 4. + private static bool LucasLehmer(EInteger mersenne, long exponent) + { + var s = EInteger.FromInt32(4); + for (long i = 0; i < exponent - 2; i++) + s = s.Multiply(s).Subtract(EInteger.FromInt32(2)).Mod(mersenne); + return s.IsZero; + } + /// /// The first prime at or after , or /// where the search would leave the machine word or run past the steps allowed. diff --git a/Sources/AngouriMath/Numerics/NumberStatements.cs b/Sources/AngouriMath/Numerics/NumberStatements.cs index ba79475b4..03aa2c528 100644 --- a/Sources/AngouriMath/Numerics/NumberStatements.cs +++ b/Sources/AngouriMath/Numerics/NumberStatements.cs @@ -15,8 +15,7 @@ namespace AngouriMath.Numerics /// /// Membership in the sets of numbers and divisibility, for the compilers, on the terms evaluation decides /// them: no truth value at NaN, an infinity in RR and CC as evaluation has it and in - /// none of the sets of whole numbers, and primality left open past the machine word, where the interpreter - /// leaves it. Divisibility is a statement about whole numbers, and has no truth value at any other. + /// none of the sets of whole numbers, and primality by the interpreter's own test. Divisibility is a statement about whole numbers, and has no truth value at any other. /// https://github.com/asc-community/AngouriMath/issues/2195 /// internal static class NumberStatements @@ -40,10 +39,6 @@ internal enum NumberSet { Complexes, Reals, Integers, PositiveIntegers, NonNegat _ => null }; - // 2^63: a whole double below it is a long, and the trial division the interpreter decides primality by - // takes it. - private const double PastTheWord = 9223372036854775808.0; - /// Whether is in the set of numbers , a . public static bool? In(NumericsComplex z, int set) { @@ -65,7 +60,7 @@ internal enum NumberSet { Complexes, Reals, Integers, PositiveIntegers, NonNegat NumberSet.Integers => true, NumberSet.PositiveIntegers => z.Real > 0, NumberSet.NonNegativeIntegers => z.Real >= 0, - _ => z.Real < 2 ? false : z.Real >= PastTheWord ? null : Primes.IsPrime(EInteger.FromInt64((long)z.Real)) + _ => z.Real < 2 ? false : Primes.IsPrime(EFloat.FromDouble(z.Real).ToEInteger()) }; } diff --git a/Sources/Tests/UnitTests/Core/Sets/MembershipInASetOperationTest.cs b/Sources/Tests/UnitTests/Core/Sets/MembershipInASetOperationTest.cs index bba9830ab..dbf1f248f 100644 --- a/Sources/Tests/UnitTests/Core/Sets/MembershipInASetOperationTest.cs +++ b/Sources/Tests/UnitTests/Core/Sets/MembershipInASetOperationTest.cs @@ -14,39 +14,38 @@ namespace AngouriMath.Tests.Core.Sets { /// /// Membership in a union, an intersection and a difference of sets is decided from each side's by Kleene's - /// connectives, where a side may leave it open: 2^89 - 1, a Mersenne prime, is past the machine word - /// and its membership in PP is left open. A union of PP with a set excluding it was - /// False. #2219 + /// connectives, where a side may leave it open: 10^30 + 57, a prime past what the prime test proves, has + /// its membership in PP left open. A union of PP with a set excluding it was False. #2219 /// [Trait("Area", "Core")] public sealed class MembershipInASetOperationTest { [Theory] // One side open, the other excluding it: open, in either order. - [InlineData("(2^89 - 1) in (PP \\/ [0; 1])")] - [InlineData("(2^89 - 1) in ([0; 1] \\/ PP)")] - [InlineData("(2^89 - 1) in (PP \\/ (-oo; 0))")] - [InlineData("(2^89 - 1) in (PP \\/ ({2} \\/ [5; 7]))")] + [InlineData("(10^30 + 57) in (PP \\/ [0; 1])")] + [InlineData("(10^30 + 57) in ([0; 1] \\/ PP)")] + [InlineData("(10^30 + 57) in (PP \\/ (-oo; 0))")] + [InlineData("(10^30 + 57) in (PP \\/ ({2} \\/ [5; 7]))")] // One side open, the other containing it: open for an intersection, and for a difference whose right side is open. - [InlineData("(2^89 - 1) in (PP /\\ [0; 10^30])")] - [InlineData("(2^89 - 1) in ([0; 10^30] \\ PP)")] - [InlineData("(2^89 - 1) in (PP \\ [0; 1])")] + [InlineData("(10^30 + 57) in (PP /\\ [0; 10^31])")] + [InlineData("(10^30 + 57) in ([0; 10^31] \\ PP)")] + [InlineData("(10^30 + 57) in (PP \\ [0; 1])")] public void LeftOpen(string input) => Assert.IsType(input.ToEntity().Evaled); [Theory] // A union holds where either side does. - [InlineData("(2^89 - 1) in (PP \\/ [0; 10^30])", "True")] + [InlineData("(10^30 + 57) in (PP \\/ [0; 10^31])", "True")] [InlineData("2 in (PP \\/ [5; 7])", "True")] [InlineData("6 in (PP \\/ [5; 7])", "True")] [InlineData("4 in (PP \\/ [5; 7])", "False")] // An intersection fails where either side does. - [InlineData("(2^89 - 1) in (PP /\\ (-oo; 0))", "False")] - [InlineData("(2^89 - 1) in (PP /\\ [0; 1])", "False")] + [InlineData("(10^30 + 57) in (PP /\\ (-oo; 0))", "False")] + [InlineData("(10^30 + 57) in (PP /\\ [0; 1])", "False")] [InlineData("7 in (PP /\\ [5; 10^30])", "True")] // A difference fails where its left side does, or its right side holds. - [InlineData("(2^89 - 1) in ([0; 1] \\ PP)", "False")] - [InlineData("(2^89 - 1) in (PP \\ [0; 10^30])", "False")] + [InlineData("(10^30 + 57) in ([0; 1] \\ PP)", "False")] + [InlineData("(10^30 + 57) in (PP \\ [0; 10^31])", "False")] [InlineData("7 in (PP \\ [0; 5])", "True")] public void Decided(string input, string expected) => Assert.Equal(expected.ToEntity(), input.ToEntity().Evaled); diff --git a/Sources/Tests/UnitTests/Core/Sets/PrimalityTest.cs b/Sources/Tests/UnitTests/Core/Sets/PrimalityTest.cs new file mode 100644 index 000000000..09b264ef6 --- /dev/null +++ b/Sources/Tests/UnitTests/Core/Sets/PrimalityTest.cs @@ -0,0 +1,96 @@ +// +// 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; +using AngouriMath.Extensions; +using PeterO.Numbers; +using Xunit; +using static AngouriMath.Entity; + +namespace AngouriMath.Tests.Core.Sets +{ + /// + /// Primality by Miller–Rabin to the first thirteen primes, exact below 3317044064679887385961981 (Sorenson and + /// Webster, Math. Comp. 86 (2017)), and by Lucas–Lehmer for a Mersenne number; past both a witness still proves a + /// number composite and a number that passes is left open. It was trial division, which for a prime just below + /// 2^63 never finished, the square of the divisor overflowing. + /// #2230 + /// + [Trait("Area", "Core")] + public sealed class PrimalityTest + { + [Theory] + // The largest prime below 2^63, where trial division did not finish. + [InlineData("9223372036854775783 in PP", "True")] + [InlineData("9223372036854775781 in PP", "False")] + [InlineData("(2^61 - 1) in PP", "True")] + // Mersenne numbers, by Lucas–Lehmer: 2^89 - 1 and 2^127 - 1 are prime, 2^67 - 1 and 2^101 - 1 are not. + [InlineData("(2^89 - 1) in PP", "True")] + [InlineData("(2^127 - 1) in PP", "True")] + [InlineData("(2^67 - 1) in PP", "False")] + [InlineData("(2^101 - 1) in PP", "False")] + // 2^4423 - 1 is prime, and 2^4421 - 1, with a prime exponent, is not. + [InlineData("(2^4423 - 1) in PP", "True")] + [InlineData("(2^4421 - 1) in PP", "False")] + // The least strong pseudoprime to the first twelve prime bases, which 41 exposes. + [InlineData("318665857834031151167461 in PP", "False")] + // Carmichael numbers and strong pseudoprimes to base 2. + [InlineData("561 in PP", "False")] + [InlineData("41041 in PP", "False")] + [InlineData("2047 in PP", "False")] + [InlineData("3215031751 in PP", "False")] + // Past the exact range a witness still proves a number composite. + [InlineData("(10^30 + 1) in PP", "False")] + [InlineData("(10^30 + 3) in PP", "False")] + public void Decided(string input, string expected) + { + var task = Task.Run(() => input.ToEntity().Evaled); + Assert.True(task.Wait(TimeSpan.FromSeconds(30)), $"{input} did not finish"); + Assert.Equal(expected.ToEntity(), task.Result); + } + + /// + /// Past the exact range, a number that passes every round is left open: 10^30 + 57 is prime, and + /// 3317044064679887385961981, the least strong pseudoprime to all thirteen bases, is not. + /// + [Theory] + [InlineData("(10^30 + 57) in PP")] + [InlineData("3317044064679887385961981 in PP")] + public void LeftOpen(string input) + => Assert.IsType(input.ToEntity().Evaled); + + [Fact] + public void BelowTwoHundredThousandItIsTheSieve() + { + const int limit = 200000; + var composite = new bool[limit]; + for (var i = 2; i * i < limit; i++) + if (!composite[i]) + for (var j = i * i; j < limit; j += i) + composite[j] = true; + for (var n = 0; n < limit; n++) + Assert.True(AngouriMath.Functions.Primes.IsPrime(EInteger.FromInt32(n)) == (n >= 2 && !composite[n]), $"{n}"); + } + + [Fact] + public void BelowTenToTheTwelveItIsTrialDivision() + { + var random = new Random(2230); + for (var t = 0; t < 2000; t++) + { + var n = (long)(random.NextDouble() * 1e12) | 1; + var prime = n >= 2; + for (long d = 3; d * d <= n && prime; d += 2) + prime = n % d != 0; + Assert.True(AngouriMath.Functions.Primes.IsPrime(EInteger.FromInt64(n)) == prime, $"{n}"); + } + } + } +}