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

start cleanup of miller-rabin

db4
Doug Coleman 2008-01-12 20:07:49 -10:00
parent bdb160dd6c
commit 304aa98fe7
1 changed files with 19 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
extra/math/miller-rabin/miller-rabin.factor
View File

@ -18,9 +18,7 @@ SYMBOL: trials
: next-odd ( m -- n ) : next-odd ( m -- n )
dup even? [ 1+ ] [ 2 + ] if ; dup even? [ 1+ ] [ 2 + ] if ;
: random-bits ( m -- n ) : random-bits ( m -- n ) 2^ random ; foldable
#! Top bit is always set
2^ [ random ] keep -1 shift bitor ; foldable
: (factor-2s) ( s n -- s n ) : (factor-2s) ( s n -- s n )
dup even? [ -1 shift >r 1+ r> (factor-2s) ] when ; dup even? [ -1 shift >r 1+ r> (factor-2s) ] when ;
@ -32,29 +30,24 @@ SYMBOL: trials
] unless 0 swap (factor-2s) ; ] unless 0 swap (factor-2s) ;
:: (miller-rabin) | n prime?! | :: (miller-rabin) | n prime?! |
n dup 1 = over even? or [ n 1- factor-2s s set r set
drop f trials get [
] [ n 1- [1,b] random a set
1- factor-2s s set r set a get s get n ^mod 1 = [
trials get [ 0 count set
n 1- [1,b] random a set r get [
a get s get n ^mod 1 = [ 2^ s get * a get swap n ^mod n - -1 = [
0 count set count [ 1+ ] change
r get [ r get +
2^ s get * a get swap n ^mod n - -1 = [
count [ 1+ ] change
r get +
] when
] each
count get zero? [
f prime?!
trials get +
] when ] when
] unless ] each
drop count get zero? [
] each f prime?!
prime? trials get +
] if ; ] when
] unless
drop
] each prime? ;
TUPLE: miller-rabin-bounds ; TUPLE: miller-rabin-bounds ;
@ -62,6 +55,7 @@ TUPLE: miller-rabin-bounds ;
over { over {
{ [ dup 1 <= ] [ 3drop f ] } { [ dup 1 <= ] [ 3drop f ] }
{ [ dup 2 = ] [ 3drop t ] } { [ dup 2 = ] [ 3drop t ] }
{ [ dup even? ] [ 3drop f ] }
{ [ t ] [ [ drop trials set t (miller-rabin) ] with-scope ] } { [ t ] [ [ drop trials set t (miller-rabin) ] with-scope ] }
} cond ; } cond ;