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