math.extras: adding Stein's algorithm

diff --git a/extra/math/extras/extras-tests.factor b/extra/math/extras/extras-tests.factor
index 88240eb03ae9c6f13c06c4daf6e31d4906f44693..786062df3107e1e3b5902f57b6f01f933ed1b55e 100644 (file)
--- a/extra/math/extras/extras-tests.factor
+++ b/extra/math/extras/extras-tests.factor
@@ -151,3 +151,7 @@ tools.test ;
 { 3 } [ 12 integer-sqrt ] unit-test
 { 4 } [ 16 integer-sqrt ] unit-test
 { 44 } [ 2019 integer-sqrt ] unit-test
+
+{ 1 } [ 11 13 stein ] unit-test
+{ 2 } [ 14 52 stein ] unit-test
+{ 7 } [ 14 7 stein ] unit-test

diff --git a/extra/math/extras/extras.factor b/extra/math/extras/extras.factor
index 0c3c9c9ce66355cd82dfa521ce670beb97dc62a3..97c0ce096be560b17f20b5f4ec4e7c425d2dcbdc 100644 (file)
--- a/extra/math/extras/extras.factor
+++ b/extra/math/extras/extras.factor
@@ -358,3 +358,27 @@ PRIVATE>
         ] each
         a a sq m > [ 1 - ] when
     ] if-zero ;
+
+<PRIVATE
+
+: reduce-evens ( value u v -- value' u' v' )
+    [ 2dup [ even? ] both? ] [ [ 2 * ] [ 2/ ] [ 2/ ] tri* ] while ;
+
+: reduce-odds ( value u v -- value' u' v' )
+    [
+        [ [ dup even? ] [ 2/ ] while ] bi@
+        2dup <=> {
+            { +eq+ [ over '[ _ * ] 2dip f ] }
+            { +lt+ [ swap [ - ] keep t ] }
+            { +gt+ [ [ - ] keep t ] }
+        } case
+    ] loop ;
+
+PRIVATE>
+
+: stein ( u v -- w )
+    2dup [ zero? ] both? [ "gcd for zeros is undefined" throw ] when
+    [ dup 0 < [ neg ] when ] bi@
+    [ 1 ] 2dip reduce-evens reduce-odds 2drop ;
+
+