-USING: math.miller-rabin tools.test ;
+USING: math.miller-rabin tools.test kernel sequences ;
IN: math.miller-rabin.tests
[ f ] [ 473155932665450549999756893736999469773678960651272093993257221235459777950185377130233556540099119926369437865330559863 miller-rabin ] unit-test
[ t ] [ 37 miller-rabin ] unit-test
[ 101 ] [ 100 next-prime ] unit-test
[ t ] [ 2135623355842621559 miller-rabin ] unit-test
-[ 100000000000031 ] [ 100000000000000 next-prime ] unit-test
\ No newline at end of file
+[ 100000000000031 ] [ 100000000000000 next-prime ] unit-test
+
+[ 863 ] [ 862 next-safe-prime ] unit-test
+[ f ] [ 862 safe-prime? ] unit-test
+[ t ] [ 7 safe-prime? ] unit-test
+[ f ] [ 31 safe-prime? ] unit-test
+[ t ] [ 863 safe-prime? ] unit-test
+
+[ f ] [ 1000 [ drop 15 miller-rabin ] any? ] unit-test
! Copyright (C) 2008 Doug Coleman.
! See http://factorcode.org/license.txt for BSD license.
USING: combinators kernel locals math math.functions math.ranges
-random sequences sets ;
+random sequences sets combinators.short-circuit ;
IN: math.miller-rabin
<PRIVATE
n 1 - :> n-1
n-1 factor-2s :> s :> r
0 :> a!
-
+ t :> prime?!
trials [
- drop
- n-1 [1,b] random a!
+ n 1 - [1,b] random a!
a s n ^mod 1 = [
- f
- ] [
- r [ 2^ s * a swap n ^mod n - -1 = ] any?
- ] if
- ] any? ;
-
+ r iota [
+ 2^ s * a swap n ^mod n - -1 =
+ ] any? not [ f prime?! trials + ] when
+ ] unless drop
+ ] each prime? ;
+
PRIVATE>
: next-odd ( m -- n ) dup even? [ 1 + ] [ 2 + ] if ;
dup 5 < [ too-few-primes ] when
2dup [ random-prime ] curry replicate
dup all-unique? [ 2nip ] [ drop unique-primes ] if ;
+
+! Safe primes are of the form p = 2q + 1, p,q are prime
+! See http://en.wikipedia.org/wiki/Safe_prime
+
+<PRIVATE
+
+: >safe-prime-form ( q -- p ) 2 * 1 + ;
+
+: safe-prime-candidate? ( n -- ? )
+ >safe-prime-form
+ 1 + 6 divisor? ;
+
+: next-safe-prime-candidate ( n -- candidate )
+ 1 - 2/
+ next-prime dup safe-prime-candidate?
+ [ next-safe-prime-candidate ] unless ;
+
+PRIVATE>
+
+: safe-prime? ( q -- ? )
+ {
+ [ 1 - 2 / dup integer? [ miller-rabin ] [ drop f ] if ]
+ [ miller-rabin ]
+ } 1&& ;
+
+: next-safe-prime ( n -- q )
+ next-safe-prime-candidate
+ dup >safe-prime-form
+ dup miller-rabin
+ [ nip ] [ drop next-safe-prime ] if ;
+
+: random-safe-prime ( numbits -- p )
+ random-bits next-safe-prime ;