libs/multiprecision/test/test_miller_rabin.cpp

   1 ///////////////////////////////////////////////////////////////
   2 //  Copyright 2012 John Maddock. Distributed under the Boost
   3 //  Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
   5
   6 #ifdef _MSC_VER
   7 #define _SCL_SECURE_NO_WARNINGS
   8 #endif
   9
  10 #include <boost/multiprecision/gmp.hpp>
  11 #include <boost/multiprecision/cpp_int.hpp>
  12 #include <boost/multiprecision/miller_rabin.hpp>
  13 #include <boost/math/special_functions/prime.hpp>
  14 #include <iostream>
  15 #include <iomanip>
  16 #include "test.hpp"
  17
  18 template <class I>
  19 void test()
  20 {
  21    //
  22    // Very simple test program to verify that the GMP's Miller-Rabin
  23    // implementation and this one agree on whether some random numbers
  24    // are prime or not.  Of course these are probabilistic tests so there's
  25    // no reason why they should actually agree - except the probability of
  26    // disagreement for 25 trials is almost infinitely small.
  27    //
  28    using namespace boost::random;
  29    using namespace boost::multiprecision;
  30
  31    typedef I test_type;
  32
  33    static const unsigned test_bits =
  34        std::numeric_limits<test_type>::digits && (std::numeric_limits<test_type>::digits <= 256)
  35            ? std::numeric_limits<test_type>::digits
  36            : 128;
  37
  38    independent_bits_engine<mt11213b, test_bits, test_type> gen;
  39    //
  40    // We must use a different generator for the tests and number generation, otherwise
  41    // we get false positives.  Further we use the same random number engine for the
  42    // Miller Rabin test as GMP uses internally:
  43    //
  44    mt19937 gen2;
  45
  46    //
  47    // Begin by testing the primes in our table as all these should return true:
  48    //
  49    for (unsigned i = 1; i < boost::math::max_prime; ++i)
  50    {
  51       BOOST_TEST(miller_rabin_test(test_type(boost::math::prime(i)), 25, gen));
  52       BOOST_TEST(mpz_probab_prime_p(mpz_int(boost::math::prime(i)).backend().data(), 25));
  53    }
  54    //
  55    // Now test some random values and compare GMP's native routine with ours.
  56    //
  57    for (unsigned i = 0; i < 10000; ++i)
  58    {
  59       test_type n              = gen();
  60       bool      is_prime_boost = miller_rabin_test(n, 25, gen2);
  61       bool      is_gmp_prime   = mpz_probab_prime_p(mpz_int(n).backend().data(), 25) ? true : false;
  62       if (is_prime_boost && is_gmp_prime)
  63       {
  64          std::cout << "We have a prime: " << std::hex << std::showbase << n << std::endl;
  65       }
  66       if (is_prime_boost != is_gmp_prime)
  67          std::cout << std::hex << std::showbase << "n = " << n << std::endl;
  68       BOOST_CHECK_EQUAL(is_prime_boost, is_gmp_prime);
  69    }
  70 }
  71
  72 int main()
  73 {
  74    using namespace boost::multiprecision;
  75
  76    test<mpz_int>();
  77    test<number<gmp_int, et_off> >();
  78    test<boost::uint64_t>();
  79    test<boost::uint32_t>();
  80
  81    test<cpp_int>();
  82    test<number<cpp_int_backend<64, 64, unsigned_magnitude, checked, void>, et_off> >();
  83    test<checked_uint128_t>();
  84    test<checked_uint1024_t>();
  85
  86    return boost::report_errors();
  87 }