steve
/
factor
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

fix miller-rabin, it's correct but a little ugly still. bed time

db4
Doug Coleman 2009-05-06 00:54:14 -05:00
parent d168f76ab0
commit 5a4270f777
2 changed files with 52 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
basis/math/miller-rabin/miller-rabin-tests.factor
View File

@ -1,4 +1,4 @@
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
@ -9,3 +9,11 @@ IN: math.miller-rabin.tests
[ 101 ] [ 100 next-prime ] unit-test
[ t ] [ 2135623355842621559 miller-rabin ] unit-test
[ 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

50
basis/math/miller-rabin/miller-rabin.factor
View File

@ -1,7 +1,7 @@
! 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
@ -14,16 +14,15 @@ TUPLE: positive-even-expected n ;
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>
@ -71,3 +70,36 @@ ERROR: too-few-primes ;
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 ;