diff --git a/src/tests/util/numerics/CMakeLists.txt b/src/tests/util/numerics/CMakeLists.txt index c60b11ec11..57ca5f7b1c 100644 --- a/src/tests/util/numerics/CMakeLists.txt +++ b/src/tests/util/numerics/CMakeLists.txt @@ -19,3 +19,6 @@ add_test(float ${CMAKE_CURRENT_BINARY_DIR}/float) add_executable(xnumeral xnumeral.cpp) target_link_libraries(xnumeral ${EXTRA_LIBS}) add_test(xnumeral ${CMAKE_CURRENT_BINARY_DIR}/xnumeral) +add_executable(primes primes.cpp) +target_link_libraries(primes ${EXTRA_LIBS}) +add_test(primes ${CMAKE_CURRENT_BINARY_DIR}/primes) diff --git a/src/tests/util/numerics/primes.cpp b/src/tests/util/numerics/primes.cpp new file mode 100644 index 0000000000..a2ea9a1099 --- /dev/null +++ b/src/tests/util/numerics/primes.cpp @@ -0,0 +1,31 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Leonardo de Moura +*/ +#include +#include "util/test.h" +#include "util/numerics/primes.h" +using namespace lean; + +#define NUM_SMALL_PRIMES 168 +static uint64 small_primes[NUM_SMALL_PRIMES] = { + 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997 +}; + +static void tst1() { + prime_iterator it; + for (unsigned i = 0; i < NUM_SMALL_PRIMES; i++) { + uint64 p = it.next(); + lean_assert_eq(p, small_primes[i]); + std::cout << p << " "; + } + std::cout << "\n"; +} + +int main() { + tst1(); + return has_violations() ? 1 : 0; +} + diff --git a/src/util/numerics/CMakeLists.txt b/src/util/numerics/CMakeLists.txt index 335deb96aa..48979cb8ba 100644 --- a/src/util/numerics/CMakeLists.txt +++ b/src/util/numerics/CMakeLists.txt @@ -1,2 +1,2 @@ -add_library(numerics gmp_init.cpp mpz.cpp mpq.cpp mpbq.cpp mpfp.cpp float.cpp double.cpp numeric_traits.cpp) +add_library(numerics gmp_init.cpp mpz.cpp mpq.cpp mpbq.cpp mpfp.cpp float.cpp double.cpp numeric_traits.cpp primes.cpp) target_link_libraries(numerics ${LEAN_LIBS} ${EXTRA_LIBS}) diff --git a/src/util/numerics/primes.cpp b/src/util/numerics/primes.cpp new file mode 100644 index 0000000000..6b623e921f --- /dev/null +++ b/src/util/numerics/primes.cpp @@ -0,0 +1,109 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Leonardo de Moura +*/ +#include +#include +#include "util/uint64.h" +#include "util/debug.h" +#include "util/exception.h" +#include "util/numerics/primes.h" + +#ifndef LEAN_PRIME_LIST_MAX_SIZE +#define LEAN_PRIME_LIST_MAX_SIZE 1<<20 +#endif + +namespace lean { +class prime_generator { + std::vector m_primes; + void process_next_k_numbers(uint64 k) { + std::vector todo; + uint64 begin = m_primes.back() + 2; + uint64 end = begin + k; + for (uint64 i = begin; i < end; i += 2) { + todo.push_back(i); + } + unsigned j = 1; + lean_assert(m_primes[j] == 3); + while (!todo.empty()) { + unsigned sz = m_primes.size(); + for (; j < sz; j++) { + uint64 p = m_primes[j]; + unsigned todo_sz = todo.size(); + unsigned k1 = 0; + unsigned k2 = 0; + for (; k1 < todo_sz; k1++) { + if (todo[k1] % p == 0) + continue; + todo[k2] = todo[k1]; + k2++; + } + todo.resize(k2); + if (k2 == 0) + return; + if (p > (todo[k2-1] / p) + 1) { + // all numbers in todo are primes + for (unsigned k1 = 0; k1 < k2; k1++) { + m_primes.push_back(todo[k1]); + } + return; + } + } + uint64 p = m_primes.back(); + p = p*p; + unsigned todo_sz = todo.size(); + unsigned k1 = 0; + for (k1 = 0; k1 < todo_sz; k1++) { + if (todo[k1] < p) { + m_primes.push_back(todo[k1]); + } + break; + } + unsigned k2 = 0; + for (; k1 < todo_sz; k1++, k2++) { + todo[k2] = todo[k1]; + } + todo.resize(k2); + } + } + +public: + prime_generator() { + m_primes.push_back(2); + m_primes.push_back(3); + process_next_k_numbers(128); + } + + uint64 operator()(unsigned idx) { + if (idx < m_primes.size()) + return m_primes[idx]; + if (idx > LEAN_PRIME_LIST_MAX_SIZE) + throw exception("prime generator capacity exceeded"); + process_next_k_numbers(1024); + if (idx < m_primes.size()) + return m_primes[idx]; + while (idx <= m_primes.size()) + process_next_k_numbers(1024*16); + return m_primes[idx]; + } +}; + +static prime_generator g_prime_generator; +static std::mutex g_prime_generator_mutex; + + +prime_iterator::prime_iterator(): + m_idx(0) { +} + +uint64 prime_iterator::next() { + unsigned idx = m_idx; + m_idx++; + { + std::lock_guard guard(g_prime_generator_mutex); + return g_prime_generator(idx); + } +} +} diff --git a/src/util/numerics/primes.h b/src/util/numerics/primes.h new file mode 100644 index 0000000000..7cb33fd1ba --- /dev/null +++ b/src/util/numerics/primes.h @@ -0,0 +1,19 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Leonardo de Moura +*/ +#include "util/uint64.h" + +namespace lean { +/** + \brief Prime number iterator. It can be used to enumerate the first LEAN_PRIME_LIST_MAX_SIZE primes. +*/ +class prime_iterator { + unsigned m_idx; +public: + prime_iterator(); + uint64 next(); +}; +} diff --git a/src/util/uint64.h b/src/util/uint64.h new file mode 100644 index 0000000000..e6679b3615 --- /dev/null +++ b/src/util/uint64.h @@ -0,0 +1,10 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Leonardo de Moura +*/ +namespace lean { +typedef unsigned long long uint64; +static_assert(sizeof(uint64) == 8, "unexpected uint64 size"); +}