本文介绍了一个由Christian Stigen Larsen实现的米勒-拉宾素性测试算法，该算法是一种概率性的素数判定方法。文章包含了完整的C语言实现代码，详细解释了算法原理及其实现细节。
https://github.com/cslarsen/miller-rabin/blob/master/miller-rabin.cpp
/*
* The Miller-Rabin primality test
*
* Written by Christian Stigen Larsen, 2012-01-10
* http://csl.sublevel3.org
*
* Distributed under the modified BSD license
*
* NOTE: I implemented this probabilistic algorithm purely as a recreational
* challenge. The code has big room for improvements, but it does work
* as advertised.
*/
#include <stdlib.h>
// rand
#include <stdint.h>
// uint64_t
#include "miller-rabin.h"
/*
* Which PRNG function to use; libc rand() by default
*/
static
int
(
*
rand_func
)(
void
)
=
rand
;
/*
* Maximum number that rand_func can return.
*/
static
int
rand_max
=
RAND_MAX
;
/*
* Fast calculation of `a^x mod n´ by using right-to-left
* binary modular exponentiation.
*
* This algorithm is taken from Bruce Schneier's book
* APPLIED CRYPTOGRAPHY.
*
* See http://en.wikipedia.org/wiki/Modular_exponentiation
*/
static
uint64_t
pow_mod
(
uint64_t
a
,
uint64_t
x
,
uint64_t
n
)
{
/*
* Note that this code is sensitive to overflowing for testing
* of large prime numbers. The `a*r´ and `a*a´ operations can
* overflow. One easy way of solving this is to use 128-bit
* precision for calculating a*b % n, since the mod operator
* should always get us back to 64bits again.
*
* You can either use GCC's built-in __int128_t or use
*
* typedef unsigned int uint128_t __attribute__((mode(TI)));
*
* to create a 128-bit datatype.
*/
uint64_t
r
=
1
;
while
(
x
)
{
if
(
(
x
&
1
)
==
1
)
//r = (__int128_t)a*r % n; // Slow
r
=
a
*
r
%
n
;
x
>>=
1
;
//a = (__int128_t)a*a % n; // SLow
a
=
a
*
a
%
n
;
}
return
r
;
}
/*
* Return an integer between a and b.
*
* Note that we use rand() here, meaning that all its pathological cases
* will apply here as well --- i.e., it's slow and not very random --- but
* should suffice.
*
*/
static
uint64_t
rand_between
(
uint64_t
a
,
uint64_t
b
)
{
// Assume rand_func() is 32 bits
uint64_t
r
=
(
static_cast
<
uint64_t
>
(
rand_func
())
<<
32
)
|
rand_func
();
return
a
+
(
uint64_t
)((
double
)(
b
-
a
+
1
)
*
r
/
(
UINT64_MAX
+
1.0
));
}
/*
* The Miller-Rabin probabilistic primality test.
*
* Returns true if ``n´´ is PROBABLY prime, false if it's composite.
* The parameter ``k´´ is the accuracy.
*
* The running time should be somewhere around O(k log_3 n).
*
*/
bool
isprime
(
uint64_t
n
,
int
k
)
{
// Must have ODD n greater than THREE
if
(
n
==
2
||
n
==
3
)
return
true
;
if
(
n
<=
1
||
!
(
n
&
1
)
)
return
false
;
// Write n-1 as d*2^s by factoring powers of 2 from n-1
int
s
=
0
;
for
(
uint64_t
m
=
n
-
1
;
!
(
m
&
1
);
++
s
,
m
>>=
1
)
;
// loop
uint64_t
d
=
(
n
-
1
)
/
(
1
<<
s
);
for
(
int
i
=
0
;
i
<
k
;
++
i
)
{
uint64_t
a
=
rand_between
(
2
,
n
-
2
);
uint64_t
x
=
pow_mod
(
a
,
d
,
n
);
if
(
x
==
1
||
x
==
n
-
1
)
continue
;
for
(
int
r
=
1
;
r
<=
s
-
1
;
++
r
)
{
x
=
pow_mod
(
x
,
2
,
n
);
if
(
x
==
1
)
return
false
;
if
(
x
==
n
-
1
)
goto
LOOP
;
}
return
false
;
LOOP:
continue
;
}
// n is *probably* prime
return
true
;
}
/*
* Set which rand function to use.
*
* If passed a NULL parameter, it will revert back to the default libc
* rand().
*/
void
setrand
(
int
(
*
pf
)(
void
),
const
int
rmax
)
{
if
(
pf
!=
NULL
)
{
rand_func
=
pf
;
rand_max
=
rmax
;
}
else
{
rand_func
=
rand
;
rand_max
=
RAND_MAX
;
}
}
扫一扫
1318
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
