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

refactor miller-rabin a bit 

still uses too many locals, but at least they're not symbols
db4
Doug Coleman 2008-05-10 13:06:40 -05:00
parent 4f1e524142
commit 9da8bed8f9
1 changed files with 21 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -3,21 +3,9 @@ math.functions math.ranges namespaces random sequences
hashtables sets ; hashtables sets ;
IN: math.miller-rabin IN: math.miller-rabin
SYMBOL: a : >even ( n -- int ) dup even? [ 1- ] unless ; foldable
SYMBOL: n : >odd ( n -- int ) dup even? [ 1+ ] when ; foldable
SYMBOL: r : next-odd ( m -- n ) dup even? [ 1+ ] [ 2 + ] if ;
SYMBOL: s
SYMBOL: count
SYMBOL: trials
: >even ( n -- int )
dup even? [ 1- ] unless ; foldable
: >odd ( n -- int )
dup even? [ 1+ ] when ; foldable
: next-odd ( m -- n )
dup even? [ 1+ ] [ 2 + ] if ;
TUPLE: positive-even-expected n ; TUPLE: positive-even-expected n ;
@ -28,34 +16,30 @@ TUPLE: positive-even-expected n ;
#! factor an integer into s * 2^r #! factor an integer into s * 2^r
0 swap (factor-2s) ; 0 swap (factor-2s) ;
:: (miller-rabin) ( n prime?! -- ? ) :: (miller-rabin) ( n trials -- ? )
n 1- factor-2s s set r set [let | r [ n 1- factor-2s drop ]
trials get [ s [ n 1- factor-2s nip ]
n 1- [1,b] random a set prime?! [ t ]
a get s get n ^mod 1 = [ a! [ 0 ]
0 count set count! [ 0 ] |
r get [ trials [
2^ s get * a get swap n ^mod n - -1 = [ n 1- [1,b] random a!
count [ 1+ ] change a s n ^mod 1 = [
r get + 0 count!
] when r [
2^ s * a swap n ^mod n - -1 =
[ count 1+ count! r + ] when
] each ] each
count get zero? [ count zero? [ f prime?! trials + ] when
f prime?! ] unless drop
trials get + ] each prime? ] ;
] when
] unless
drop
] each prime? ;
TUPLE: miller-rabin-bounds ;
: miller-rabin* ( n numtrials -- ? ) : miller-rabin* ( n numtrials -- ? )
over { over {
{ [ dup 1 <= ] [ 3drop f ] } { [ dup 1 <= ] [ 3drop f ] }
{ [ dup 2 = ] [ 3drop t ] } { [ dup 2 = ] [ 3drop t ] }
{ [ dup even? ] [ 3drop f ] } { [ dup even? ] [ 3drop f ] }
[ [ drop trials set t (miller-rabin) ] with-scope ] [ [ drop (miller-rabin) ] with-scope ]
} cond ; } cond ;
: miller-rabin ( n -- ? ) 10 miller-rabin* ; : miller-rabin ( n -- ? ) 10 miller-rabin* ;