Document miller-rabin, more unit tests for some corner cases

diff --git a/basis/math/miller-rabin/miller-rabin-docs.factor b/basis/math/miller-rabin/miller-rabin-docs.factor
new file mode 100644 (file)
index 0000000..4aa318f
--- /dev/null
+++ b/basis/math/miller-rabin/miller-rabin-docs.factor
@@ -0,0 +1,100 @@
+! Copyright (C) 2009 Doug Coleman.
+! See http://factorcode.org/license.txt for BSD license.
+USING: help.markup help.syntax kernel sequences math ;
+IN: math.miller-rabin
+
+HELP: find-relative-prime
+{ $values
+    { "n" integer }
+    { "p" integer }
+}
+{ $description "Returns a number that is relatively prime to " { $snippet "n" } "." } ;
+
+HELP: find-relative-prime*
+{ $values
+    { "n" integer } { "guess" integer }
+    { "p" integer }
+}
+{ $description "Returns a number that is relatively prime to " { $snippet "n" } ", starting by trying " { $snippet "guess" } "." } ;
+
+HELP: miller-rabin
+{ $values
+    { "n" integer }
+    { "?" "a boolean" }
+}
+{ $description "Returns true if the number is a prime. Calls " { $link miller-rabin* } " with a default of 10 Miller-Rabin tests." } ;
+
+{ miller-rabin miller-rabin* } related-words
+
+HELP: miller-rabin*
+{ $values
+    { "n" integer } { "numtrials" integer }
+    { "?" "a boolean" }
+}
+{ $description "Performs " { $snippet "numtrials" } " trials of the Miller-Rabin probabilistic primality test algorithm and returns true if prime." } ;
+
+HELP: next-prime
+{ $values
+    { "n" integer }
+    { "p" integer }
+}
+{ $description "Tests consecutive numbers for primality with " { $link miller-rabin } " and returns the next prime." } ;
+
+HELP: next-safe-prime
+{ $values
+    { "n" integer }
+    { "q" integer }
+}
+{ $description "Tests consecutive numbers and returns the next safe prime. A safe prime is desirable in cryptography applications such as Diffie-Hellman and SRP6." } ;
+
+HELP: random-bits*
+{ $values
+    { "numbits" integer }
+    { "n" integer }
+}
+{ $description "Returns an integer exactly " { $snippet "numbits" } " in length, with the topmost bit set to one." } ;
+
+HELP: random-prime
+{ $values
+    { "numbits" integer }
+    { "p" integer }
+}
+{ $description "Returns a prime number exactly " { $snippet "numbits" } " bits in length, with the topmost bit set to one." } ;
+
+HELP: random-safe-prime
+{ $values
+    { "numbits" integer }
+    { "p" integer }
+}
+{ $description "Returns a safe prime number " { $snippet "numbits" } " bits in length, with the topmost bit set to one." } ;
+
+HELP: safe-prime?
+{ $values
+    { "q" integer }
+    { "?" "a boolean" }
+}
+{ $description "Tests whether the number is a safe prime. A safe prime " { $snippet "p" } " must be prime, as must " { $snippet "(p - 1) / 2" } "." } ;
+
+HELP: unique-primes
+{ $values
+    { "numbits" integer } { "n" integer }
+    { "seq" sequence }
+}
+{ $description "Generates a sequence of " { $snippet "n" } " unique prime numbers with exactly " { $snippet "numbits" } " bits." } ;
+
+ARTICLE: "math.miller-rabin" "Miller-Rabin probabilistic primality test"
+"The " { $vocab-link "math.miller-rabin" } " vocabulary implements the Miller-Rabin probabilistic primality test and utility words that use it in order to generate random prime numbers." $nl
+"The Miller-Rabin probabilistic primality test:"
+{ $subsection miller-rabin }
+{ $subsection miller-rabin* }
+"Generating relative prime numbers:"
+{ $subsection find-relative-prime }
+{ $subsection find-relative-prime* }
+"Generating prime numbers:"
+{ $subsection next-prime }
+{ $subsection random-prime }
+"Generating safe prime numbers:"
+{ $subsection next-safe-prime }
+{ $subsection random-safe-prime } ;
+
+ABOUT: "math.miller-rabin"

diff --git a/basis/math/miller-rabin/miller-rabin-tests.factor b/basis/math/miller-rabin/miller-rabin-tests.factor
index 676c4bf20d93c52d57a43f170cee7c488a3fb1e8..9981064ec076dbaaa6916e83473dfa489548fcfc 100644 (file)
--- a/basis/math/miller-rabin/miller-rabin-tests.factor
+++ b/basis/math/miller-rabin/miller-rabin-tests.factor
@@ -1,4 +1,5 @@
-USING: math.miller-rabin tools.test kernel sequences ;
+USING: math.miller-rabin tools.test kernel sequences
+math.miller-rabin.private math ;
 IN: math.miller-rabin.tests
 
 [ f ] [ 473155932665450549999756893736999469773678960651272093993257221235459777950185377130233556540099119926369437865330559863 miller-rabin ] unit-test
@@ -6,6 +7,9 @@ IN: math.miller-rabin.tests
 [ t ] [ 3 miller-rabin ] unit-test
 [ f ] [ 36 miller-rabin ] unit-test
 [ t ] [ 37 miller-rabin ] unit-test
+[ 2 ] [ 1 next-prime ] unit-test
+[ 3 ] [ 2 next-prime ] unit-test
+[ 5 ] [ 3 next-prime ] unit-test
 [ 101 ] [ 100 next-prime ] unit-test
 [ t ] [ 2135623355842621559 miller-rabin ] unit-test
 [ 100000000000031 ] [ 100000000000000 next-prime ] unit-test
@@ -14,6 +18,12 @@ IN: math.miller-rabin.tests
 [ f ] [ 862 safe-prime? ] unit-test
 [ t ] [ 7 safe-prime? ] unit-test
 [ f ] [ 31 safe-prime? ] unit-test
+[ t ] [ 47 safe-prime-candidate? ] unit-test
+[ t ] [ 47 safe-prime? ] unit-test
 [ t ] [ 863 safe-prime? ] unit-test
 
 [ f ] [ 1000 [ drop 15 miller-rabin ] any? ] unit-test
+
+[ 47 ] [ 31 next-safe-prime ] unit-test
+[ 49 ] [ 50 random-prime log2 ] unit-test
+[ 49 ] [ 50 random-bits* log2 ] unit-test

diff --git a/basis/math/miller-rabin/miller-rabin.factor b/basis/math/miller-rabin/miller-rabin.factor
index 5e999aa9561959d0a24dfd5b4766c6c5813b72e2..9fd604a00378c3a1bc9d9f4facf8647ba6542ac9 100755 (executable)
--- a/basis/math/miller-rabin/miller-rabin.factor
+++ b/basis/math/miller-rabin/miller-rabin.factor
@@ -1,15 +1,20 @@
 ! Copyright (C) 2008 Doug Coleman.
 ! See http://factorcode.org/license.txt for BSD license.
 USING: combinators kernel locals math math.functions math.ranges
-random sequences sets combinators.short-circuit math.bitwise ;
+random sequences sets combinators.short-circuit math.bitwise
+math math.order ;
 IN: math.miller-rabin
 
 <PRIVATE
 
-: >odd ( n -- int ) dup even? [ 1 + ] when ; foldable
+: >odd ( n -- int ) 0 set-bit ; foldable
 
 : >even ( n -- int ) 0 clear-bit ; foldable
 
+: next-even ( m -- n ) >even 2 + ;
+
+: next-odd ( m -- n ) dup even? [ 1 + ] [ 2 + ] if ;
+
 TUPLE: positive-even-expected n ;
 
 :: (miller-rabin) ( n trials -- ? )
@@ -18,7 +23,7 @@ TUPLE: positive-even-expected n ;
     0 :> a!
     trials [
         drop
-        n 1 - [1,b] random a!
+        2 n 2 - [a,b] random a!
         a s n ^mod 1 = [
             f
         ] [
@@ -30,8 +35,6 @@ TUPLE: positive-even-expected n ;
 
 PRIVATE>
 
-: next-odd ( m -- n ) dup even? [ 1 + ] [ 2 + ] if ;
-
 : miller-rabin* ( n numtrials -- ? )
     over {
         { [ dup 1 <= ] [ 3drop f ] }
@@ -42,11 +45,21 @@ PRIVATE>
 
 : miller-rabin ( n -- ? ) 10 miller-rabin* ;
 
+ERROR: prime-range-error n ;
+
 : next-prime ( n -- p )
-    next-odd dup miller-rabin [ next-prime ] unless ;
+    dup 1 < [ prime-range-error ] when
+    dup 1 = [
+        drop 2
+    ] [
+        next-odd dup miller-rabin [ next-prime ] unless
+    ] if ;
+
+: random-bits* ( numbits -- n )
+    1 - [ random-bits ] keep set-bit ;
 
 : random-prime ( numbits -- p )
-    random-bits next-prime ;
+    random-bits* next-prime ;
 
 ERROR: no-relative-prime n ;
 
@@ -80,10 +93,7 @@ ERROR: too-few-primes ;
 
 <PRIVATE
 
-: >safe-prime-form ( q -- p ) 2 * 1 + ;
-
 : safe-prime-candidate? ( n -- ? )
-    >safe-prime-form
     1 + 6 divisor? ;
 
 : next-safe-prime-candidate ( n -- candidate )
@@ -99,14 +109,8 @@ PRIVATE>
     } 1&& ;
 
 : next-safe-prime ( n -- q )
-    1 - >even 2 /
     next-safe-prime-candidate
-    dup >safe-prime-form
-    dup miller-rabin
-    [ nip ] [ drop next-safe-prime ] if ;
-
-: random-bits* ( numbits -- n )
-    [ random-bits ] keep set-bit ;
+    dup safe-prime? [ next-safe-prime ] unless ;
 
 : random-safe-prime ( numbits -- p )
-    1- random-bits* next-safe-prime ;
+    random-bits* next-safe-prime ;