Skip to content

Primality is decided by Miller-Rabin and Lucas-Lehmer, and finishes #2232

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Oct 30, 2026
Merged

Primality is decided by Miller-Rabin and Lucas-Lehmer, and finishes #2232

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 23 additions & 5 deletions BREAKING-CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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`,
Expand Down Expand Up @@ -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

Expand Down
3 changes: 3 additions & 0 deletions Sources/.editorconfig
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
4 changes: 2 additions & 2 deletions Sources/AngouriMath/Core/Entity/Omni/Entity.Omni.Classes.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -761,8 +761,8 @@ public override bool TryContains(Entity entity, out bool contains)
}

/// <summary>
/// Whether the whole number is prime, or <see langword="null"/> where it does
/// not fit a machine word and trial division is not attempted.
/// Whether the whole number is prime, or <see langword="null"/> where it is not
/// decided: see <see cref="Functions.Primes.IsPrime"/>.
/// </summary>
internal static bool? IsDecidablyPrime(Integer whole)
=> Functions.Primes.IsPrime(whole.EInteger);
Expand Down
81 changes: 70 additions & 11 deletions Sources/AngouriMath/Functions/NumberTheory/Primes.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

/// <summary>
/// Whether the whole number is prime, or <see langword="null"/> 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 <c>PP</c> is left open rather than answered slowly or
/// wrongly. https://github.com/asc-community/AngouriMath/issues/1450
/// Whether the whole number is prime, or <see langword="null"/> where it is not decided. Miller–Rabin to the
/// first thirteen primes is exact below <c>3317044064679887385961981</c>; past it a witness still proves a
/// number composite, a Mersenne number <c>2^p - 1</c> 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
/// </summary>
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;
}

/// <summary>
/// The first prime at or after <paramref name="atLeast"/>, or <see langword="null"/>
/// where the search would leave the machine word or run past the steps allowed.
Expand Down
9 changes: 2 additions & 7 deletions Sources/AngouriMath/Numerics/NumberStatements.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,7 @@ namespace AngouriMath.Numerics
/// <summary>
/// Membership in the sets of numbers and divisibility, for the compilers, on the terms evaluation decides
/// them: no truth value at <c>NaN</c>, an infinity in <c>RR</c> and <c>CC</c> 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
/// </summary>
internal static class NumberStatements
Expand All @@ -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;

/// <summary>Whether <paramref name="z"/> is in the set of numbers <paramref name="set"/>, a <see cref="NumberSet"/>.</summary>
public static bool? In(NumericsComplex z, int set)
{
Expand All @@ -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())
};
}

Expand Down
29 changes: 14 additions & 15 deletions Sources/Tests/UnitTests/Core/Sets/MembershipInASetOperationTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,39 +14,38 @@ namespace AngouriMath.Tests.Core.Sets
{
/// <summary>
/// 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: <c>2^89 - 1</c>, a Mersenne prime, is past the machine word
/// and its membership in <c>PP</c> is left open. A union of <c>PP</c> with a set excluding it was
/// <c>False</c>. <a href="https://github.com/asc-community/AngouriMath/issues/2219">#2219</a>
/// connectives, where a side may leave it open: <c>10^30 + 57</c>, a prime past what the prime test proves, has
/// its membership in <c>PP</c> left open. A union of <c>PP</c> with a set excluding it was <c>False</c>. <a href="https://github.com/asc-community/AngouriMath/issues/2219">#2219</a>
/// </summary>
[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<Set.Inf>(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);
Expand Down
96 changes: 96 additions & 0 deletions Sources/Tests/UnitTests/Core/Sets/PrimalityTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
{
/// <summary>
/// Primality by Miller–Rabin to the first thirteen primes, exact below <c>3317044064679887385961981</c> (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
/// <c>2^63</c> never finished, the square of the divisor overflowing.
/// <a href="https://github.com/asc-community/AngouriMath/issues/2230">#2230</a>
/// </summary>
[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);
}

/// <summary>
/// Past the exact range, a number that passes every round is left open: <c>10^30 + 57</c> is prime, and
/// <c>3317044064679887385961981</c>, the least strong pseudoprime to all thirteen bases, is not.
/// </summary>
[Theory]
[InlineData("(10^30 + 57) in PP")]
[InlineData("3317044064679887385961981 in PP")]
public void LeftOpen(string input)
=> Assert.IsType<Set.Inf>(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}");
}
}
}
}
Loading