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 kernel locals math math.functions math.ranges
   4 random sequences sets ;
   5 IN: math.miller-rabin
   6
   7 <PRIVATE
   8
   9 : >odd ( n -- int ) dup even? [ 1+ ] when ; foldable
  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 PRIVATE>
  32
  33 : next-odd ( m -- n ) dup even? [ 1+ ] [ 2 + ] if ;
  34
  35 : miller-rabin* ( n numtrials -- ? )
  36     over {
  37         { [ dup 1 <= ] [ 3drop f ] }
  38         { [ dup 2 = ] [ 3drop t ] }
  39         { [ dup even? ] [ 3drop f ] }
  40         [ drop (miller-rabin) ]
  41     } cond ;
  42
  43 : miller-rabin ( n -- ? ) 10 miller-rabin* ;
  44
  45 : next-prime ( n -- p )
  46     next-odd dup miller-rabin [ next-prime ] unless ;
  47
  48 : random-prime ( numbits -- p )
  49     random-bits next-prime ;
  50
  51 ERROR: no-relative-prime n ;
  52
  53 <PRIVATE
  54
  55 : (find-relative-prime) ( n guess -- p )
  56     over 1 <= [ over no-relative-prime ] when
  57     dup 1 <= [ drop 3 ] when
  58     2dup gcd nip 1 > [ 2 + (find-relative-prime) ] [ nip ] if ;
  59
  60 PRIVATE>
  61
  62 : find-relative-prime* ( n guess -- p )
  63     #! find a prime relative to n with initial guess
  64     >odd (find-relative-prime) ;
  65
  66 : find-relative-prime ( n -- p )
  67     dup random find-relative-prime* ;
  68
  69 ERROR: too-few-primes ;
  70
  71 : unique-primes ( numbits n -- seq )
  72     #! generate two primes
  73     swap
  74     dup 5 < [ too-few-primes ] when
  75     2dup [ random-prime ] curry replicate
  76     dup all-unique? [ 2nip ] [ drop unique-primes ] if ;