forked from mpusz/mp-units
Chip Hogg
cfd9ddb675
This can mark a number as either "probably prime", or "definitely composite". The first parameter is the base, and the second is the number to test. Future PRs will build up the Strong Lucas test which complements this, and then form the Baillie-PSW test by composing the two. Helps #506.
126 lines
5.7 KiB
C++
126 lines
5.7 KiB
C++
// The MIT License (MIT)
|
|
//
|
|
// Copyright (c) 2018 Mateusz Pusz
|
|
//
|
|
// Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
// of this software and associated documentation files (the "Software"), to deal
|
|
// in the Software without restriction, including without limitation the rights
|
|
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
// copies of the Software, and to permit persons to whom the Software is
|
|
// furnished to do so, subject to the following conditions:
|
|
//
|
|
// The above copyright notice and this permission notice shall be included in all
|
|
// copies or substantial portions of the Software.
|
|
//
|
|
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
// SOFTWARE.
|
|
|
|
#include <mp-units/ext/prime.h>
|
|
#ifdef MP_UNITS_IMPORT_STD
|
|
import std;
|
|
#else
|
|
#include <array>
|
|
#include <cstddef>
|
|
#include <utility>
|
|
#endif
|
|
|
|
using namespace mp_units::detail;
|
|
|
|
namespace {
|
|
|
|
inline constexpr auto MAX_U64 = std::numeric_limits<std::uint64_t>::max();
|
|
|
|
template<std::size_t BasisSize, std::size_t... Is>
|
|
constexpr bool check_primes(std::index_sequence<Is...>)
|
|
{
|
|
return ((Is < 2 || wheel_factorizer<BasisSize>::is_prime(Is) == is_prime_by_trial_division(Is)) && ...);
|
|
}
|
|
|
|
static_assert(check_primes<2>(std::make_index_sequence<122>{}));
|
|
|
|
// This is the smallest number that can catch the bug where we use only _prime_ numbers in the first wheel, rather than
|
|
// numbers which are _coprime to the basis_.
|
|
//
|
|
// The basis for N = 4 is {2, 3, 5, 7}, so the wheel size is 210. 11 * 11 = 121 is within the first wheel. It is
|
|
// coprime with every element of the basis, but it is _not_ prime. If we keep only prime numbers, then we will neglect
|
|
// using numbers of the form (210 * n + 121) as trial divisors, which is a problem if any are prime. For n = 1, we have
|
|
// a divisor of (210 + 121 = 331), which happens to be prime but will not be used. Thus, (331 * 331 = 109561) is a
|
|
// composite number which could wrongly appear prime if we skip over 331.
|
|
static_assert(wheel_factorizer<4>::is_prime(109'561) == is_prime_by_trial_division(109'561));
|
|
|
|
static_assert(wheel_factorizer<1>::coprimes_in_first_wheel.size() == 1);
|
|
static_assert(wheel_factorizer<2>::coprimes_in_first_wheel.size() == 2);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel.size() == 8);
|
|
static_assert(wheel_factorizer<4>::coprimes_in_first_wheel.size() == 48);
|
|
static_assert(wheel_factorizer<5>::coprimes_in_first_wheel.size() == 480);
|
|
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[0] == 1);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[1] == 7);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[2] == 11);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[3] == 13);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[4] == 17);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[5] == 19);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[6] == 23);
|
|
static_assert(wheel_factorizer<3>::coprimes_in_first_wheel[7] == 29);
|
|
|
|
static_assert(!wheel_factorizer<1>::is_prime(0));
|
|
static_assert(!wheel_factorizer<1>::is_prime(1));
|
|
static_assert(wheel_factorizer<1>::is_prime(2));
|
|
|
|
static_assert(!wheel_factorizer<2>::is_prime(0));
|
|
static_assert(!wheel_factorizer<2>::is_prime(1));
|
|
static_assert(wheel_factorizer<2>::is_prime(2));
|
|
|
|
static_assert(!wheel_factorizer<3>::is_prime(0));
|
|
static_assert(!wheel_factorizer<3>::is_prime(1));
|
|
static_assert(wheel_factorizer<3>::is_prime(2));
|
|
|
|
// Modular arithmetic.
|
|
static_assert(add_mod(1u, 2u, 5u) == 3u);
|
|
static_assert(add_mod(4u, 4u, 5u) == 3u);
|
|
static_assert(add_mod(MAX_U64 - 1u, MAX_U64 - 2u, MAX_U64) == MAX_U64 - 3u);
|
|
|
|
static_assert(sub_mod(2u, 1u, 5u) == 1u);
|
|
static_assert(sub_mod(1u, 2u, 5u) == 4u);
|
|
static_assert(sub_mod(MAX_U64 - 2u, MAX_U64 - 1u, MAX_U64) == MAX_U64 - 1u);
|
|
static_assert(sub_mod(1u, MAX_U64 - 1u, MAX_U64) == 2u);
|
|
|
|
static_assert(mul_mod(6u, 7u, 10u) == 2u);
|
|
static_assert(mul_mod(13u, 11u, 50u) == 43u);
|
|
static_assert(mul_mod(MAX_U64 / 2u, 10u, MAX_U64) == MAX_U64 - 5u);
|
|
|
|
static_assert(half_mod_odd(0u, 11u) == 0u);
|
|
static_assert(half_mod_odd(10u, 11u) == 5u);
|
|
static_assert(half_mod_odd(1u, 11u) == 6u);
|
|
static_assert(half_mod_odd(9u, 11u) == 10u);
|
|
static_assert(half_mod_odd(MAX_U64 - 1u, MAX_U64) == (MAX_U64 - 1u) / 2u);
|
|
static_assert(half_mod_odd(MAX_U64 - 2u, MAX_U64) == MAX_U64 - 1u);
|
|
|
|
static_assert(pow_mod(5u, 8u, 9u) == ((5u * 5u * 5u * 5u) * (5u * 5u * 5u * 5u)) % 9u);
|
|
static_assert(pow_mod(2u, 64u, MAX_U64) == 1u);
|
|
|
|
// Miller-Rabin primality testing.
|
|
static_assert(miller_rabin_probable_prime(2u, 5u));
|
|
static_assert(miller_rabin_probable_prime(2u, 7u));
|
|
static_assert(!miller_rabin_probable_prime(2u, 9u));
|
|
static_assert(miller_rabin_probable_prime(2u, 11u));
|
|
|
|
static_assert(miller_rabin_probable_prime(2u, 2047u), "Known base 2 pseudoprime");
|
|
static_assert(miller_rabin_probable_prime(2u, 3277u), "Known base 2 pseudoprime");
|
|
|
|
static_assert(miller_rabin_probable_prime(3u, 121u), "Known base 3 pseudoprime");
|
|
static_assert(miller_rabin_probable_prime(3u, 703u), "Known base 3 pseudoprime");
|
|
|
|
static_assert(miller_rabin_probable_prime(2u, 225'653'407'801u), "Large known prime");
|
|
static_assert(miller_rabin_probable_prime(2u, 334'524'384'739u), "Large known prime");
|
|
static_assert(miller_rabin_probable_prime(2u, 9'007'199'254'740'881u), "Large known prime");
|
|
|
|
static_assert(miller_rabin_probable_prime(2u, 18'446'744'073'709'551'557u), "Largest 64-bit prime");
|
|
|
|
} // namespace
|