From 5a4270f77749221fbfcd70160bee4b8d9e2d4201 Mon Sep 17 00:00:00 2001
From: Doug Coleman <erg@jobim.local>
Date: Wed, 6 May 2009 00:54:14 -0500
Subject: [PATCH] fix miller-rabin, it's correct but a little ugly still.  bed
 time

---
 .../miller-rabin/miller-rabin-tests.factor    | 12 ++++-
 basis/math/miller-rabin/miller-rabin.factor   | 52 +++++++++++++++----
 2 files changed, 52 insertions(+), 12 deletions(-)

diff --git a/basis/math/miller-rabin/miller-rabin-tests.factor b/basis/math/miller-rabin/miller-rabin-tests.factor
index 5f1b9835e4..676c4bf20d 100644
--- a/basis/math/miller-rabin/miller-rabin-tests.factor
+++ b/basis/math/miller-rabin/miller-rabin-tests.factor
@@ -1,4 +1,4 @@
-USING: math.miller-rabin tools.test ;
+USING: math.miller-rabin tools.test kernel sequences ;
 IN: math.miller-rabin.tests
 
 [ f ] [ 473155932665450549999756893736999469773678960651272093993257221235459777950185377130233556540099119926369437865330559863 miller-rabin ] unit-test
@@ -8,4 +8,12 @@ IN: math.miller-rabin.tests
 [ t ] [ 37 miller-rabin ] unit-test
 [ 101 ] [ 100 next-prime ] unit-test
 [ t ] [ 2135623355842621559 miller-rabin ] unit-test
-[ 100000000000031 ] [ 100000000000000 next-prime ] unit-test
\ No newline at end of file
+[ 100000000000031 ] [ 100000000000000 next-prime ] unit-test
+
+[ 863 ] [ 862 next-safe-prime ] unit-test
+[ f ] [ 862 safe-prime? ] unit-test
+[ t ] [ 7 safe-prime? ] unit-test
+[ f ] [ 31 safe-prime? ] unit-test
+[ t ] [ 863 safe-prime? ] unit-test
+
+[ f ] [ 1000 [ drop 15 miller-rabin ] any? ] unit-test
diff --git a/basis/math/miller-rabin/miller-rabin.factor b/basis/math/miller-rabin/miller-rabin.factor
index 62d8ee4432..93d7f4c582 100755
--- a/basis/math/miller-rabin/miller-rabin.factor
+++ b/basis/math/miller-rabin/miller-rabin.factor
@@ -1,7 +1,7 @@
 ! 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 ;
+random sequences sets combinators.short-circuit ;
 IN: math.miller-rabin
 
 <PRIVATE
@@ -14,17 +14,16 @@ TUPLE: positive-even-expected n ;
     n 1 - :> n-1
     n-1 factor-2s :> s :> r
     0 :> a!
-
+    t :> prime?!
     trials [
-        drop
-        n-1 [1,b] random a!
+        n 1 - [1,b] random a!
         a s n ^mod 1 = [
-            f
-        ] [
-            r [ 2^ s * a swap n ^mod n - -1 = ] any?
-        ] if
-    ] any? ;
-    
+            r iota [
+                2^ s * a swap n ^mod n - -1 =
+            ] any? not [ f prime?! trials + ] when
+        ] unless drop
+    ] each prime? ;
+
 PRIVATE>
 
 : next-odd ( m -- n ) dup even? [ 1 + ] [ 2 + ] if ;
@@ -71,3 +70,36 @@ ERROR: too-few-primes ;
     dup 5 < [ too-few-primes ] when
     2dup [ random-prime ] curry replicate
     dup all-unique? [ 2nip ] [ drop unique-primes ] if ;
+
+! Safe primes are of the form p = 2q + 1, p,q are prime
+! See http://en.wikipedia.org/wiki/Safe_prime
+
+<PRIVATE
+
+: >safe-prime-form ( q -- p ) 2 * 1 + ;
+
+: safe-prime-candidate? ( n -- ? )
+    >safe-prime-form
+    1 + 6 divisor? ;
+
+: next-safe-prime-candidate ( n -- candidate )
+    1 - 2/
+    next-prime dup safe-prime-candidate?
+    [ next-safe-prime-candidate ] unless ;
+
+PRIVATE>
+
+: safe-prime? ( q -- ? )
+    {
+        [ 1 - 2 / dup integer? [ miller-rabin ] [ drop f ] if ]
+        [ miller-rabin ]
+    } 1&& ;
+
+: next-safe-prime ( n -- q )
+    next-safe-prime-candidate
+    dup >safe-prime-form
+    dup miller-rabin
+    [ nip ] [ drop next-safe-prime ] if ;
+
+: random-safe-prime ( numbits -- p )
+    random-bits next-safe-prime ;