prime_numbers.hpp

Go to the documentation of this file.

/******************************************************************************
** Copyright © 2012 by J.M.McGuiness, coder@hussar.me.uk
**
** This library is free software; you can redistribute it and/or
** modify it under the terms of the GNU Lesser General Public
** License as published by the Free Software Foundation; either
** version 2.1 of the License, or (at your option) any later version.
**
** This library is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
** Lesser General Public License for more details.
**
** You should have received a copy of the GNU Lesser General Public
** License along with this library; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
 
#include "config.h"
 
#include <algorithm>
#include <bitset>
#include <vector>
 
namespace jmmcg { namespace LIBJMMCG_VER_NAMESPACE {
 
namespace dyn {
   typedef std::vector<unsigned long long> primes_colln;
 
   /// Find all of the prime numbers in the range [2...max_num) using the Sieve of Eratosthenes.
   /**
      Complexity: time: O(n(log(log(n))))
 
      \param   max_num  The largest number in the range up to which the primes should be found.
      \return  The primes within the range [2...max_num).
   */
   inline primes_colln
   sieve_of_eratosthenes(primes_colln::size_type max_num) noexcept(false) {
      bool is_prime[max_num];
      std::fill_n(is_prime, max_num, true);
      for (primes_colln::size_type i=2; i<std::sqrt(max_num); ++i) {
         if (is_prime[i]) {
            for (primes_colln::size_type j=i*i; j<max_num; j+=i) {
               is_prime[j]=false;
            }
         }
      }
      primes_colln::size_type num_primes=0;
      for (primes_colln::size_type i=2; i<max_num; ++i) {
         if (is_prime[i]) {
            ++num_primes;
         }
      }
      assert(num_primes<=max_num);
      primes_colln primes;
      primes.reserve(num_primes);
      for (primes_colln::size_type i=2; i<max_num; ++i) {
         if (is_prime[i]) {
            primes.push_back(i);
         }
      }
      return primes;
   }
}
 
namespace mpl {
   template<class E, E MN>
   class sieve_of_eratosthenes {
   public:
      typedef E element_type;
      typedef std::vector<element_type> primes_colln;
 
      static constexpr element_type max_num=MN;
 
   private:
      static constexpr std::bitset<max_num>
      filter_primes() noexcept(true) {
         std::bitset<max_num> is_prime;
         for (typename primes_colln::size_type i=2; i<std::sqrt(max_num); ++i) {
            if (!is_prime[i]) {
               for (std::size_t j=i*i; j<max_num; j+=i) {
                  is_prime[j]=true;
               }
            }
         }
         return is_prime;
      }
 
 
      static constexpr typename primes_colln::size_type
      count_primes(std::bitset<max_num> const &is_prime) noexcept(true) __attribute__((pure)) {
         typename primes_colln::size_type num_primes=0;
         for (typename primes_colln::size_type i=2; i<max_num; ++i) {
            if (is_prime[i]) {
               ++num_primes;
            }
         }
         return num_primes;
      }
 
   public:
      static primes_colln
      result() noexcept(false) {
         const std::bitset<max_num> is_prime(filter_primes());
         primes_colln primes;
         primes.reserve(count_primes(is_prime));
         for (typename primes_colln::size_type i=2; i<max_num; ++i) {
            if (!is_prime[i]) {
               primes.push_back(i);
            }
         }
         return primes;
      }
   };
}
 
} }