v1.4.1

README.md

@@ -42,7 +42,7 @@ namespace cppr {
 ```
 
 It returns true if the input value is a prime number; otherwise, it returns false.
-If you want to reduce the size of the executable file, use this function instead of `cppr::IsPrime` because `cppr::IsPrime` uses a 40KB table for performance optimization.
+If you want to reduce the size of the executable file, use this function instead of `cppr::IsPrime` because `cppr::IsPrime` uses a 640KB table for performance optimization.
 
 #### example
 
@@ -188,8 +188,42 @@ Workflow: [bench.yml](https://github.com/sortA0329/libcpprime/actions/workflows/
     />
 </p>
 
+## Algorithm
+
+The core algorithm of `libcpprime` is the Miller–Rabin primality test. Although traditionally a randomized algorithm, it has been proven that any integer below 2^64 can be deterministically tested using seven fixed bases [1].
+
+Using a hash table optimizes this process further. For composite numbers that act as pseudoprimes to a specific base (such as base 2), their hash values are calculated to assign them into separate buckets. Bases are then selected for each bucket to correctly test every number mapped to it [2]. Finding these optimal bases took under a minute, aided by GPU acceleration and precomputed base-2 pseudoprimes.
+
+Thanks to this strategy, `cppr::IsPrime` requires only 2 bases, while `cppr::IsPrimeCompact` uses 5. For numbers smaller than 2^32, Bradley Berg's algorithm is applied, which relies on a single base [3].
+
+Beyond core algorithmic choices, micro-optimizations play a critical role. Because CPU division takes longer than multiplication, Montgomery reduction is employed. For numbers below 2^32, accepting both `MR(x)` and `MR(x) + n` (where `MR` denotes the Montgomery representation) eliminates unpredictable conditional branches. For values smaller than 2^21, fast modulo reduction using Lemire's method is applied [4] [5].
+
+Furthermore, tests against multiple bases are executed concurrently within a single loop. This design maximizes instruction-level parallelism (ILP) and minimizes conditional branching. Although this increases execution time for composite numbers that would otherwise be filtered out early in sequential testing, `libcpprime` prioritizes optimizing worst-case performance.
+
+Additional micro-optimizations include:
+
+- Trial division by small primes
+- Precomputed lookup flags for small numbers below 2^21 (`cppr::IsPrime`) or 2^10 (`cppr::IsPrimeCompact`)
+- Skipping the second step of the Miller–Rabin test (repeated squaring) when n ≡ 3 mod 4
+- Newton's method for computing `-n^-1 mod R` [6]
+- GCD calculation with products of small primes using Stein's algorithm (binary GCD) instead of division [7]
+- Standard division for numbers between 2^21 and 2^32
+- Inline assembly (GCC/Clang) and `_udiv128` (MSVC) for 128-bit modulo arithmetic, avoiding overhead from compiler built-in functions for `unsigned __int128`
+- `libdivide` for 128-bit modulo arithmetic when hardware/intrinsic support is unavailable
+- Compiler hints (e.g., `__builtin_assume`) providing value ranges for improved code generation
+
+[1] https://miller-rabin.appspot.com/
+[2] https://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html
+[3] https://www.techneon.com/download/is.prime.32.base.data
+[4] https://lemire.me/blog/2016/06/27/a-fast-alternative-to-the-modulo-reduction/
+[5] https://en.algorithmica.org/hpc/arithmetic/division/
+[6] https://rsk0315.hatenablog.com/entry/2022/11/27/060616
+[7] https://lpha-z.hatenablog.com/entry/2020/05/31/231500
+
 ## Releases
 
+- 2026/10/04 v1.4.1
+  - Improve performance of `cppr::IsPrime`, which resulted in a 600KB increase in table size
 - 2026/10/01 v1.4.0
   - Rename `cppr::IsPrimeNoTable` to `cppr::IsPrimeCompact`
   - Improve performance of `cppr::IsPrimeCompact`
@@ -247,11 +281,3 @@ Workflow: [bench.yml](https://github.com/sortA0329/libcpprime/actions/workflows/
   - Add `cppr::IsPrime` with a table
 - 2024/12/18 v1.0.0
   - Add `cppr::IsPrime`
-
-## References
-
-- https://miller-rabin.appspot.com/
-- https://zenn.dev/mizar/articles/791698ea860581
-- https://www.techneon.com/download/is.prime.32.base.data
-- https://www.techneon.com/download/is.prime.64.base.data
-- https://lemire.me/blog/2016/06/27/a-fast-alternative-to-the-modulo-reduction/

include/libcpprime/IsPrime.hpp

@@ -6,24 +6,6 @@
  * SPDX-License-Identifier: MIT
  *
  **/
-/**
- *
- * The algorithm in this library is based on Bradley Berg's method.
- * See this page for more information:
- * https://www.techneon.com/download/is.prime.64.base.data
- *
- * Copyright 2018 Bradley Berg   < (My last name) @ t e c h n e o n . c o m >
- *
- * Permission to use, copy, modify, and distribute this software for any
- * purpose with or without fee is hereby granted, provided that the above
- * copyright notice and this permission notice appear in all copies.
- *
- * This algorithm is deliberately unpatented. The license above also
- * lets you even freely use it in commercial code.
- *
- * Primality testing using a hash table of bases originated with Steven Worley.
- *
- **/
 
 #ifndef CPPR_INTERNAL_INCLUDED_IS_PRIME
 #define CPPR_INTERNAL_INCLUDED_IS_PRIME
@@ -38,23 +20,23 @@ namespace cppr {
 
 namespace internal {
 
-constexpr std::uint64_t FlagTable17[1024] = {
-#include "internal/IsPrimeTable17.txt"
+constexpr std::uint64_t FlagTable21[16384] = {
+#include "internal/IsPrimeTable21.txt"
 };
-CPPR_INTERNAL_CONSTEXPR_INLINE bool IsPrime17(const std::uint64_t n) noexcept { return n == 2 || (n % 2 == 1 && (FlagTable17[n / 128] & (1ull << (n % 128 / 2)))); }
+CPPR_INTERNAL_CONSTEXPR_INLINE bool IsPrime21(const std::uint64_t n) noexcept { return n == 2 || (n % 2 == 1 && (FlagTable21[n / 128] & (1ull << (n % 128 / 2)))); }
 
-constexpr std::uint16_t Bases64[16384] = {
+constexpr std::uint16_t Bases64[262144] = {
 #include "internal/IsPrimeBases64.txt"
 };
-CPPR_INTERNAL_CONSTEXPR_INLINE std::uint16_t GetBase(std::uint64_t x) noexcept { return Bases64[(0xad625b89u * static_cast<std::uint32_t>(x)) >> 18]; }
-CPPR_INTERNAL_CONSTEXPR_INLINE bool IsPrime49(const std::uint64_t x) noexcept {
-    const MontgomeryModint64Impl<false> mint(x);
+template <bool Strict>
+CPPR_INTERNAL_CONSTEXPR_INLINE bool IsPrime64(const std::uint64_t x) noexcept {
+    const MontgomeryModint64Impl<Strict> mint(x);
     const std::int32_t S = CountrZero(x - 1);
     const std::uint64_t D = (x - 1) >> S;
     const auto one = mint.one();
     const auto mone = mint.mone();
     auto c = mint.raw(2);
-    auto d = mint.raw(GetBase(x));
+    auto d = mint.raw(Bases64[(0xad625b89u * static_cast<std::uint32_t>(x)) >> 14]);
     auto a = c;
     auto b = d;
     if (D != 1) {
@@ -87,73 +69,18 @@ CPPR_INTERNAL_CONSTEXPR_INLINE bool IsPrime49(const std::uint64_t x) noexcept {
     }
     return res1 && res2;
 }
-template <bool Strict>
-CPPR_INTERNAL_CONSTEXPR_INLINE bool IsPrime64(const std::uint64_t x) noexcept {
-    const MontgomeryModint64Impl<Strict> mint(x);
-    const std::int32_t S = CountrZero(x - 1);
-    const std::uint64_t D = (x - 1) >> S;
-    const auto one = mint.one();
-    const auto mone = mint.mone();
-    const std::uint32_t base = GetBase(x);
-    const std::uint64_t base_mask = static_cast<std::uint64_t>(15ull | (135ull << 8) | (13ull << 16) | (60ull << 24) | (15ull << 32) | (117ull << 40) | (65ull << 48) | (29ull << 56));
-    auto d = mint.raw(2);
-    auto e = mint.raw(base);
-    auto f = mint.raw((base_mask >> (8 * (base >> 13))) & 0xff);
-    auto a = d;
-    auto b = e;
-    auto c = f;
-    if (D != 1) {
-        d = mint.mul(d, d);
-        e = mint.mul(e, e);
-        f = mint.mul(f, f);
-        std::uint64_t ex = D >> 1;
-        while (ex != 1) {
-            const auto g = mint.mul(d, d);
-            const auto h = mint.mul(e, e);
-            const auto i = mint.mul(f, f);
-            if (ex & 1) {
-                a = mint.mul(a, d);
-                b = mint.mul(b, e);
-                c = mint.mul(c, f);
-            }
-            d = g;
-            e = h;
-            f = i;
-            ex >>= 1;
-        }
-    }
-    a = mint.mul(a, d);
-    b = mint.mul(b, e);
-    c = mint.mul(c, f);
-    bool res1 = mint.same(a, one) || mint.same(a, mone);
-    bool res2 = mint.same(b, one) || mint.same(b, mone);
-    bool res3 = mint.same(c, one) || mint.same(c, mone);
-    if (x % 4 == 1 && !(res1 && res2 && res3)) {
-        for (std::int32_t i = 0; i != S - 1; ++i) {
-            a = mint.mul(a, a);
-            b = mint.mul(b, b);
-            c = mint.mul(c, c);
-            res1 |= mint.same(a, mone);
-            res2 |= mint.same(b, mone);
-            res3 |= mint.same(c, mone);
-        }
-    }
-    return res1 && res2 && res3;
-}
 
 }  // namespace internal
 
 CPPR_INTERNAL_CONSTEXPR bool IsPrime(std::uint64_t n) noexcept {
-    if (n < 131072) {
-        return internal::IsPrime17(n);
+    if (n < (1ull << 21)) {
+        return internal::IsPrime21(n);
     } else if (n <= 0xffffffff) {
         if (internal::TrialDivision32(static_cast<std::uint32_t>(n))) return false;
         return internal::IsPrime32(static_cast<std::uint32_t>(n));
     } else {
         if (internal::TrialDivision64(n)) return false;
-        if (n < (std::uint64_t(1) << 49)) {
-            return internal::IsPrime49(n);
-        } else if (n < (std::uint64_t(1) << 62)) {
+        if (n < (std::uint64_t(1) << 62)) {
             return internal::IsPrime64<false>(n);
         } else {
             return internal::IsPrime64<true>(n);

include/libcpprime/IsPrimeCompact.hpp

@@ -156,7 +156,7 @@ CPPR_INTERNAL_CONSTEXPR_INLINE bool IsPrime64Compact(const std::uint64_t x) noex
 }  // namespace internal
 
 CPPR_INTERNAL_CONSTEXPR bool IsPrimeCompact(std::uint64_t n) noexcept {
-    if (n < 1024) {
+    if (n < (1ull << 10)) {
         return internal::IsPrime10(n);
     } else if (n <= 0xffffffff) {
         if (internal::TrialDivision32(static_cast<std::uint32_t>(n))) return false;