basis/math/miller-rabin/miller-rabin.factor

   1 ! Copyright (C) 2008 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: combinators io locals kernel math math.functions
   4 math.ranges namespaces random sequences hashtables sets ;
   5 IN: math.miller-rabin
   6
   7 : >even ( n -- int ) dup even? [ 1- ] unless ; foldable
   8 : >odd ( n -- int ) dup even? [ 1+ ] when ; foldable
   9 : next-odd ( m -- n ) dup even? [ 1+ ] [ 2 + ] if ;
  10
  11 TUPLE: positive-even-expected n ;
  12
  13 :: (miller-rabin) ( n trials -- ? )
  14     [let | r [ n 1- factor-2s drop ]
  15            s [ n 1- factor-2s nip ]
  16            prime?! [ t ]
  17            a! [ 0 ]
  18            count! [ 0 ] |
  19         trials [
  20             n 1- [1,b] random a!
  21             a s n ^mod 1 = [
  22                 0 count!
  23                 r [
  24                     2^ s * a swap n ^mod n - -1 =
  25                     [ count 1+ count! r + ] when
  26                 ] each
  27                 count zero? [ f prime?! trials + ] when
  28             ] unless drop
  29         ] each prime? ] ;
  30
  31 : miller-rabin* ( n numtrials -- ? )
  32     over {
  33         { [ dup 1 <= ] [ 3drop f ] }
  34         { [ dup 2 = ] [ 3drop t ] }
  35         { [ dup even? ] [ 3drop f ] }
  36         [ [ drop (miller-rabin) ] with-scope ]
  37     } cond ;
  38
  39 : miller-rabin ( n -- ? ) 10 miller-rabin* ;
  40
  41 : next-prime ( n -- p )
  42     next-odd dup miller-rabin [ next-prime ] unless ;
  43
  44 : random-prime ( numbits -- p )
  45     random-bits next-prime ;
  46
  47 ERROR: no-relative-prime n ;
  48
  49 : (find-relative-prime) ( n guess -- p )
  50     over 1 <= [ over no-relative-prime ] when
  51     dup 1 <= [ drop 3 ] when
  52     2dup gcd nip 1 > [ 2 + (find-relative-prime) ] [ nip ] if ;
  53
  54 : find-relative-prime* ( n guess -- p )
  55     #! find a prime relative to n with initial guess
  56     >odd (find-relative-prime) ;
  57
  58 : find-relative-prime ( n -- p )
  59     dup random find-relative-prime* ;
  60
  61 ERROR: too-few-primes ;
  62
  63 : unique-primes ( numbits n -- seq )
  64     #! generate two primes
  65     swap
  66     dup 5 < [ too-few-primes ] when
  67     2dup [ random-prime ] curry replicate
  68     dup all-unique? [ 2nip ] [ drop unique-primes ] if ;