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

   1 ! Copyright (c) 2008-2009 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 combinators.short-circuit math.bitwise
   5 math math.order ;
   6 IN: math.miller-rabin
   7
   8 : >odd ( n -- int ) 0 set-bit ; foldable
   9
  10 : >even ( n -- int ) 0 clear-bit ; foldable
  11
  12 : next-even ( m -- n ) >even 2 + ;
  13
  14 : next-odd ( m -- n ) dup even? [ 1 + ] [ 2 + ] if ;
  15
  16 <PRIVATE
  17
  18 :: (miller-rabin) ( n trials -- ? )
  19     n 1 - :> n-1
  20     n-1 factor-2s :> s :> r
  21     0 :> a!
  22     trials [
  23         drop
  24         2 n 2 - [a,b] random a!
  25         a s n ^mod 1 = [
  26             f
  27         ] [
  28             r iota [
  29                 2^ s * a swap n ^mod n - -1 =
  30             ] any? not
  31         ] if
  32     ] any? not ;
  33
  34 PRIVATE>
  35
  36 : miller-rabin* ( n numtrials -- ? )
  37     over {
  38         { [ dup 1 <= ] [ 3drop f ] }
  39         { [ dup 2 = ] [ 3drop t ] }
  40         { [ dup even? ] [ 3drop f ] }
  41         [ drop (miller-rabin) ]
  42     } cond ;
  43
  44 : miller-rabin ( n -- ? ) 10 miller-rabin* ;
  45
  46 ERROR: prime-range-error n ;
  47
  48 : next-prime ( n -- p )
  49     dup 1 < [ prime-range-error ] when
  50     dup 1 = [
  51         drop 2
  52     ] [
  53         next-odd dup miller-rabin [ next-prime ] unless
  54     ] if ;
  55
  56 : random-bits* ( numbits -- n )
  57     1 - [ random-bits ] keep set-bit ;
  58
  59 : random-prime ( numbits -- p )
  60     random-bits* next-prime ;
  61
  62 ERROR: no-relative-prime n ;
  63
  64 <PRIVATE
  65
  66 : (find-relative-prime) ( n guess -- p )
  67     over 1 <= [ over no-relative-prime ] when
  68     dup 1 <= [ drop 3 ] when
  69     2dup gcd nip 1 > [ 2 + (find-relative-prime) ] [ nip ] if ;
  70
  71 PRIVATE>
  72
  73 : find-relative-prime* ( n guess -- p )
  74     #! find a prime relative to n with initial guess
  75     >odd (find-relative-prime) ;
  76
  77 : find-relative-prime ( n -- p )
  78     dup random find-relative-prime* ;
  79
  80 ERROR: too-few-primes ;
  81
  82 : unique-primes ( numbits n -- seq )
  83     #! generate two primes
  84     swap
  85     dup 5 < [ too-few-primes ] when
  86     2dup [ random-prime ] curry replicate
  87     dup all-unique? [ 2nip ] [ drop unique-primes ] if ;
  88
  89 ! Safe primes are of the form p = 2q + 1, p,q are prime
  90 ! See http://en.wikipedia.org/wiki/Safe_prime
  91
  92 <PRIVATE
  93
  94 : safe-prime-candidate? ( n -- ? )
  95     1 + 6 divisor? ;
  96
  97 : next-safe-prime-candidate ( n -- candidate )
  98     next-prime dup safe-prime-candidate?
  99     [ next-safe-prime-candidate ] unless ;
 100
 101 PRIVATE>
 102
 103 : safe-prime? ( q -- ? )
 104     {
 105         [ 1 - 2 / dup integer? [ miller-rabin ] [ drop f ] if ]
 106         [ miller-rabin ]
 107     } 1&& ;
 108
 109 : next-safe-prime ( n -- q )
 110     next-safe-prime-candidate
 111     dup safe-prime? [ next-safe-prime ] unless ;
 112
 113 : random-safe-prime ( numbits -- p )
 114     random-bits* next-safe-prime ;