|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r77353 - in sandbox/big_number: boost/multiprecision libs/multiprecision/doc libs/multiprecision/doc/html libs/multiprecision/doc/html/boost_multiprecision libs/multiprecision/example libs/multiprecision/test
From: john_at_[hidden]
Date: 2012-03-17 05:18:22
Author: johnmaddock
Date: 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
New Revision: 77353
URL: http://svn.boost.org/trac/boost/changeset/77353
Log:
Add simple version of the Miller Rabin test
Added:
sandbox/big_number/boost/multiprecision/miller_rabin.hpp (contents, props changed)
sandbox/big_number/libs/multiprecision/example/safe_prime.cpp (contents, props changed)
sandbox/big_number/libs/multiprecision/test/test_miller_rabin.cpp (contents, props changed)
Text files modified:
sandbox/big_number/boost/multiprecision/random.hpp | 2 ++
sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/ref.html | 6 +++---
sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/tut.html | 1 +
sandbox/big_number/libs/multiprecision/doc/html/index.html | 3 ++-
sandbox/big_number/libs/multiprecision/doc/multiprecision.qbk | 28 ++++++++++++++++++++++++++++
sandbox/big_number/libs/multiprecision/test/Jamfile.v2 | 8 ++++++++
6 files changed, 44 insertions(+), 4 deletions(-)
Added: sandbox/big_number/boost/multiprecision/miller_rabin.hpp
==============================================================================
--- (empty file)
+++ sandbox/big_number/boost/multiprecision/miller_rabin.hpp 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -0,0 +1,106 @@
+///////////////////////////////////////////////////////////////
+// Copyright 2012 John Maddock. Distributed under the Boost
+// Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
+
+#ifndef BOOST_MP_MR_HPP
+#define BOOST_MP_MR_HPP
+
+#include <boost/multiprecision/random.hpp>
+
+namespace boost{
+namespace multiprecision{
+//
+// Calculate (a^b)%c:
+//
+template <class Backend, bool ExpressionTemplates>
+void expmod(const mp_number<Backend, ExpressionTemplates>& a, mp_number<Backend, ExpressionTemplates> b, const mp_number<Backend, ExpressionTemplates>& c, mp_number<Backend, ExpressionTemplates>& result)
+{
+ typedef mp_number<Backend, ExpressionTemplates> number_type;
+ number_type x(1), y(a);
+ while(b > 0)
+ {
+ if(b & 1)
+ {
+ x = (x * y) % c;
+ }
+ y = (y * y) % c;
+ b /= 2;
+ }
+ result = x % c;
+}
+
+template <class Backend, bool ExpressionTemplates, class Engine>
+bool miller_rabin_test(const mp_number<Backend, ExpressionTemplates>& n, unsigned trials, Engine& gen)
+{
+ typedef mp_number<Backend, ExpressionTemplates> number_type;
+
+ static const unsigned small_factors[] = {
+ 3u, 5u, 7u, 11u, 13u, 17u, 19u, 23u, 29u, 31u, 37u, 41u, 43u, 47u, 53u };
+
+ if(n < 2)
+ return false;
+ if((n & 1) == 0)
+ return false;
+ //
+ // Sanity check for small factors:
+ //
+
+ number_type q = (n - 1) >> 1;
+ unsigned k = 1;
+ while((q & 1) == 0)
+ {
+ q >>= 1;
+ ++k;
+ }
+ // Declare our random number generator:
+ boost::random::uniform_int_distribution<number_type> dist(0, n);
+ number_type nm1 = n - 1;
+ //
+ // Execute the trials:
+ //
+ for(unsigned i = 0; i < trials; ++i)
+ {
+ number_type x = dist(gen);
+ number_type y;
+ expmod(x, q, n, y);
+ unsigned j = 0;
+ while(true)
+ {
+ if((y == nm1) || ((y == 1) && (j == 0)))
+ break;
+ if(y == 1)
+ return false; // test failed
+ if(++j == k)
+ return false; // failed
+ y = (y * y) % n;
+ }
+ }
+ return true; // Yeheh! probably prime.
+}
+
+template <class Backend, bool ExpressionTemplates>
+bool miller_rabin_test(const mp_number<Backend, ExpressionTemplates>& x, unsigned trials)
+{
+ static mt19937 gen;
+ return miller_rabin_test(x, trials, gen);
+}
+
+template <class tag, class Arg1, class Arg2, class Arg3, class Engine>
+bool miller_rabin_test(const detail::mp_exp<tag, Arg1, Arg2, Arg3> & n, unsigned trials, Engine& gen)
+{
+ typedef typename detail::mp_exp<tag, Arg1, Arg2, Arg3>::result_type number_type;
+ return miller_rabin_test(number_type(n), trials, gen);
+}
+
+template <class tag, class Arg1, class Arg2, class Arg3>
+bool miller_rabin_test(const detail::mp_exp<tag, Arg1, Arg2, Arg3> & n, unsigned trials)
+{
+ typedef typename detail::mp_exp<tag, Arg1, Arg2, Arg3>::result_type number_type;
+ return miller_rabin_test(number_type(n), trials);
+}
+
+}} // namespaces
+
+#endif
+
Modified: sandbox/big_number/boost/multiprecision/random.hpp
==============================================================================
--- sandbox/big_number/boost/multiprecision/random.hpp (original)
+++ sandbox/big_number/boost/multiprecision/random.hpp 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -1,4 +1,6 @@
///////////////////////////////////////////////////////////////
+// Copyright Jens Maurer 2006-1011
+// Copyright Steven Watanabe 2011
// Copyright 2012 John Maddock. Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
Modified: sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/ref.html
==============================================================================
--- sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/ref.html (original)
+++ sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/ref.html 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -6,12 +6,12 @@
<meta name="generator" content="DocBook XSL Stylesheets V1.76.1">
<link rel="home" href="../index.html" title="Chapter 1. Boost.Multiprecision">
<link rel="up" href="../index.html" title="Chapter 1. Boost.Multiprecision">
-<link rel="prev" href="tut/random.html" title="Generating Random Numbers">
+<link rel="prev" href="tut/primetest.html" title="Primality Testing">
<link rel="next" href="ref/mp_number.html" title="mp_number">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<div class="spirit-nav">
-<a accesskey="p" href="tut/random.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a><a accesskey="n" href="ref/mp_number.html"><img src="../images/next.png" alt="Next"></a>
+<a accesskey="p" href="tut/primetest.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a><a accesskey="n" href="ref/mp_number.html"><img src="../images/next.png" alt="Next"></a>
</div>
<div class="section boost_multiprecision_ref">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
@@ -32,7 +32,7 @@
</tr></table>
<hr>
<div class="spirit-nav">
-<a accesskey="p" href="tut/random.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a><a accesskey="n" href="ref/mp_number.html"><img src="../images/next.png" alt="Next"></a>
+<a accesskey="p" href="tut/primetest.html"><img src="../images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../images/home.png" alt="Home"></a><a accesskey="n" href="ref/mp_number.html"><img src="../images/next.png" alt="Next"></a>
</div>
</body>
</html>
Modified: sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/tut.html
==============================================================================
--- sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/tut.html (original)
+++ sandbox/big_number/libs/multiprecision/doc/html/boost_multiprecision/tut.html 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -22,6 +22,7 @@
<dt><span class="section">Floating Point Numbers</span></dt>
<dt><span class="section">Rational Number Types</span></dt>
<dt><span class="section">Generating Random Numbers</span></dt>
+<dt><span class="section">Primality Testing</span></dt>
</dl></div>
<p>
In order to use this library you need to make two choices: what kind of number
Modified: sandbox/big_number/libs/multiprecision/doc/html/index.html
==============================================================================
--- sandbox/big_number/libs/multiprecision/doc/html/index.html (original)
+++ sandbox/big_number/libs/multiprecision/doc/html/index.html 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -34,6 +34,7 @@
<dt><span class="section">Floating Point Numbers</span></dt>
<dt><span class="section">Rational Number Types</span></dt>
<dt><span class="section">Generating Random Numbers</span></dt>
+<dt><span class="section">Primality Testing</span></dt>
</dl></dd>
<dt><span class="section">Reference</span></dt>
<dd><dl>
@@ -57,7 +58,7 @@
</div>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
-<td align="left"><p><small>Last revised: March 15, 2012 at 18:34:52 GMT</small></p></td>
+<td align="left"><p><small>Last revised: March 17, 2012 at 09:04:10 GMT</small></p></td>
<td align="right"><div class="copyright-footer"></div></td>
</tr></table>
<hr>
Modified: sandbox/big_number/libs/multiprecision/doc/multiprecision.qbk
==============================================================================
--- sandbox/big_number/libs/multiprecision/doc/multiprecision.qbk (original)
+++ sandbox/big_number/libs/multiprecision/doc/multiprecision.qbk 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -24,6 +24,7 @@
[import ../example/tommath_snips.cpp]
[import ../example/cpp_int_snips.cpp]
[import ../example/random_snips.cpp]
+[import ../example/safe_prime.cpp]
[template mpfr[] [@http://www.mpfr.org MPFR]]
[template gmp[] [@http://gmplib.org GMP]]
@@ -698,6 +699,33 @@
[endsect]
+[section:primetest Primality Testing]
+
+The library implements a fairly basic (meaning unoptimized, and possibly slow) Miller-Rabin test for primality:
+
+ #include <boost/multiprecision/miller_rabin.hpp>
+
+ template <class Backend, bool ExpressionTemplates, class Engine>
+ bool miller_rabin_test(const mp_number<Backend, ExpressionTemplates>& n, unsigned trials, Engine& gen);
+
+ template <class Backend, bool ExpressionTemplates, class Engine>
+ bool miller_rabin_test(const mp_number<Backend, ExpressionTemplates>& n, unsigned trials);
+
+These functions perform a Miller-Rabin test for primality, if the result is `false` then /n/ is definitely composite,
+while if the result is `true` then /n/ is prime with probability 0.25^trials. The algorithm used is straight out of
+Knuth Vol 2, which recomends 25 trials for a pretty strong likelyhood that /n/ is prime.
+
+The third optional argument is for a Uniform Random Number Generator from Boost.Random. When not provided the `mt19937`
+generator is used. Note that when producing random primes then you should probably use a different random number generator
+to produce possible primes, than used internally for testing whether the value is prime. It also helps of course to seed the
+generators with some source of randomness.
+
+The following example searches for a prime `p` for which `(p-1)/2` is also probably prime:
+
+[safe_prime]
+
+[endsect]
+
[endsect]
[section:ref Reference]
Added: sandbox/big_number/libs/multiprecision/example/safe_prime.cpp
==============================================================================
--- (empty file)
+++ sandbox/big_number/libs/multiprecision/example/safe_prime.cpp 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -0,0 +1,45 @@
+///////////////////////////////////////////////////////////////
+// Copyright 2012 John Maddock. Distributed under the Boost
+// Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
+
+//[safe_prime
+
+#include <boost/multiprecision/gmp.hpp>
+#include <boost/multiprecision/miller_rabin.hpp>
+#include <iostream>
+#include <iomanip>
+
+int main()
+{
+ using namespace boost::random;
+ using namespace boost::multiprecision;
+
+ typedef mpz_int int_type;
+ mt11213b base_gen(clock());
+ independent_bits_engine<mt11213b, 256, int_type> gen(base_gen);
+ //
+ // We must use a different generator for the tests and number generation, otherwise
+ // we get false positives.
+ //
+ mt19937 gen2(clock());
+
+ for(unsigned i = 0; i < 100000; ++i)
+ {
+ int_type n = gen();
+ if(miller_rabin_test(n, 25, gen2))
+ {
+ // Value n is probably prime, see if (n-1)/2 is also prime:
+ std::cout << "We have a probable prime with value: " << std::hex << std::showbase << n << std::endl;
+ if(miller_rabin_test((n-1)/2, 25, gen2))
+ {
+ std::cout << "We have a safe prime with value: " << std::hex << std::showbase << n << std::endl;
+ return 0;
+ }
+ }
+ }
+ std::cout << "Ooops, no safe primes were found" << std::endl;
+ return 1;
+}
+
+//]
Modified: sandbox/big_number/libs/multiprecision/test/Jamfile.v2
==============================================================================
--- sandbox/big_number/libs/multiprecision/test/Jamfile.v2 (original)
+++ sandbox/big_number/libs/multiprecision/test/Jamfile.v2 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -602,6 +602,14 @@
release # otherwise runtime is too slow!!
;
+run test_miller_rabin.cpp gmp
+ : # command line
+ : # input files
+ : # requirements
+ [ check-target-builds ../config//has_gmp : : <build>no ]
+ release # otherwise runtime is too slow!!
+ ;
+
run test_rational_io.cpp $(TOMMATH)
: # command line
: # input files
Added: sandbox/big_number/libs/multiprecision/test/test_miller_rabin.cpp
==============================================================================
--- (empty file)
+++ sandbox/big_number/libs/multiprecision/test/test_miller_rabin.cpp 2012-03-17 05:18:20 EDT (Sat, 17 Mar 2012)
@@ -0,0 +1,47 @@
+///////////////////////////////////////////////////////////////
+// Copyright 2012 John Maddock. Distributed under the Boost
+// Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
+
+#include <boost/multiprecision/gmp.hpp>
+#include <boost/multiprecision/miller_rabin.hpp>
+#include <iostream>
+#include <iomanip>
+#include "test.hpp"
+
+int main()
+{
+ //
+ // Very simple test program to verify that the GMP's Miller-Rabin
+ // implementation and this one agree on whether some random numbers
+ // are prime or not. Of course these are probabilistic tests so there's
+ // no reason why they should actually agree - except the probability of
+ // disagreement for 25 trials is almost infinitely small.
+ //
+ using namespace boost::random;
+ using namespace boost::multiprecision;
+ independent_bits_engine<mt11213b, 256, mpz_int> gen;
+ //
+ // We must use a different generator for the tests and number generation, otherwise
+ // we get false positives. Further we use the same random number engine for the
+ // Miller Rabin test as GMP uses internally:
+ //
+ mt19937 gen2;
+
+ for(unsigned i = 0; i < 10000; ++i)
+ {
+ mpz_int n = gen();
+ bool is_prime_boost = miller_rabin_test(n, 25, gen2);
+ bool is_gmp_prime = mpz_probab_prime_p(n.backend().data(), 25);
+ if(is_prime_boost && is_gmp_prime)
+ {
+ std::cout << "We have a prime: " << std::hex << std::showbase << n << std::endl;
+ }
+ if(is_prime_boost != is_gmp_prime)
+ std::cout << std::hex << std::showbase << "n = " << n << std::endl;
+ BOOST_CHECK_EQUAL(is_prime_boost, is_gmp_prime);
+ }
+ return 0;
+}
+
+
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk