From 24263299fa0d59680ef571eb8a3970cf64b3089a Mon Sep 17 00:00:00 2001
From: John Benediktsson <mrjbq7@gmail.com>
Date: Sat, 13 Dec 2014 16:52:28 -0800
Subject: [PATCH] math.primes.solovay-strassen: adding Solovay-Strassen
 primality test.

---
 .../math/primes/solovay-strassen/authors.txt  |  1 +
 .../solovay-strassen-tests.factor             |  7 +++++
 .../solovay-strassen/solovay-strassen.factor  | 30 +++++++++++++++++++
 3 files changed, 38 insertions(+)
 create mode 100644 extra/math/primes/solovay-strassen/authors.txt
 create mode 100644 extra/math/primes/solovay-strassen/solovay-strassen-tests.factor
 create mode 100644 extra/math/primes/solovay-strassen/solovay-strassen.factor

diff --git a/extra/math/primes/solovay-strassen/authors.txt b/extra/math/primes/solovay-strassen/authors.txt
new file mode 100644
index 0000000000..e091bb8164
--- /dev/null
+++ b/extra/math/primes/solovay-strassen/authors.txt
@@ -0,0 +1 @@
+John Benediktsson
diff --git a/extra/math/primes/solovay-strassen/solovay-strassen-tests.factor b/extra/math/primes/solovay-strassen/solovay-strassen-tests.factor
new file mode 100644
index 0000000000..359e01fa5c
--- /dev/null
+++ b/extra/math/primes/solovay-strassen/solovay-strassen-tests.factor
@@ -0,0 +1,7 @@
+USING: kernel math.primes sequences tools.test ;
+IN: math.primes.solovay-strassen
+
+{ t } [
+    100,000 iota [ solovay-strassen ] filter
+    100,000 primes-upto =
+] unit-test
diff --git a/extra/math/primes/solovay-strassen/solovay-strassen.factor b/extra/math/primes/solovay-strassen/solovay-strassen.factor
new file mode 100644
index 0000000000..12a940f2f0
--- /dev/null
+++ b/extra/math/primes/solovay-strassen/solovay-strassen.factor
@@ -0,0 +1,30 @@
+! Copyright (C) 2014 John Benediktsson
+! See http://factorcode.org/license.txt for BSD license
+
+USING: combinators kernel locals math math.extras
+math.extras.private math.functions math.ranges random sequences ;
+
+IN: math.primes.solovay-strassen
+
+<PRIVATE
+
+:: (solovay-strassen) ( n numtrials -- ? )
+    numtrials iota [
+        drop
+        n 1 - [1,b) random :> a
+        a n fast-gcd 1 > [ t ] [
+            a n jacobi n mod'
+            a n 1 - 2 /i n ^mod = not
+        ] if
+    ] any? not ;
+
+PRIVATE>
+
+: solovay-strassen* ( n numtrials -- ? )
+    {
+        { [ over 1 <= ] [ 2drop f ] }
+        { [ over even? ] [ drop 2 = ] }
+        [ (solovay-strassen) ]
+    } cond ;
+
+: solovay-strassen ( n -- ? ) 32 solovay-strassen* ;