From 2c8d3bad08cff72abc8ed879846f3149ad3cf07f Mon Sep 17 00:00:00 2001
From: Samuel Tardieu <sam@rfc1149.net>
Date: Mon, 2 Mar 2009 19:10:43 +0100
Subject: [PATCH] Continued fractions

---
 extra/math/continued-fractions/authors.txt    |  1 +
 .../continued-fractions-docs.factor           | 14 +++++++++++
 .../continued-fractions-tests.factor          | 21 ++++++++++++++++
 .../continued-fractions.factor                | 24 +++++++++++++++++++
 extra/math/continued-fractions/summary.txt    |  1 +
 5 files changed, 61 insertions(+)
 create mode 100644 extra/math/continued-fractions/authors.txt
 create mode 100644 extra/math/continued-fractions/continued-fractions-docs.factor
 create mode 100644 extra/math/continued-fractions/continued-fractions-tests.factor
 create mode 100644 extra/math/continued-fractions/continued-fractions.factor
 create mode 100644 extra/math/continued-fractions/summary.txt

diff --git a/extra/math/continued-fractions/authors.txt b/extra/math/continued-fractions/authors.txt
new file mode 100644
index 0000000000..f3b0233f74
--- /dev/null
+++ b/extra/math/continued-fractions/authors.txt
@@ -0,0 +1 @@
+Samuel Tardieu
diff --git a/extra/math/continued-fractions/continued-fractions-docs.factor b/extra/math/continued-fractions/continued-fractions-docs.factor
new file mode 100644
index 0000000000..667deb7f06
--- /dev/null
+++ b/extra/math/continued-fractions/continued-fractions-docs.factor
@@ -0,0 +1,14 @@
+USING: help.markup help.syntax ;
+IN: math.continued-fractions
+
+HELP: approx
+{ $values { "epsilon" "a positive floating point number representing the absolute acceptable error" } { "float" "a positive floating point number to approximate" } { "a/b" "a fractional number containing the approximation" } }
+{ $description "Give a rational approximation of " { $snippet "float" } " with a precision of " { $snippet "epsilon" } " using the smallest possible denominator." } ;
+
+HELP: >ratio
+{ $values { "seq" "a sequence representing a continued fraction" } { "a/b" "a fractional number" } }
+{ $description "Transform " { $snippet "seq" } " into its rational representation." } ;
+
+HELP: next-approx
+{ $values { "seq" "a mutable sequence" } }
+{ $description "Compute the next step in continued fraction calculation." } ;
diff --git a/extra/math/continued-fractions/continued-fractions-tests.factor b/extra/math/continued-fractions/continued-fractions-tests.factor
new file mode 100644
index 0000000000..d8fac0beb2
--- /dev/null
+++ b/extra/math/continued-fractions/continued-fractions-tests.factor
@@ -0,0 +1,21 @@
+USING: kernel math.constants math.continued-fractions tools.test ;
+
+[ V{ 2 2.0 } ] [ V{ 2.5 } dup next-approx ] unit-test
+[ V{ 2 2 } ] [ V{ 2.5 } dup next-approx dup next-approx ] unit-test
+
+[ 5/2 ] [ V{ 2 2.1 } >ratio ] unit-test
+[ 5/2 ] [ V{ 2 1.9 } >ratio ] unit-test
+[ 5/2 ] [ V{ 2 2.0 } >ratio ] unit-test
+[ 5/2 ] [ V{ 2 2 } >ratio ] unit-test
+
+[ 3 ] [ 1 pi approx ] unit-test
+[ 22/7 ] [ 0.1 pi approx ] unit-test
+[ 355/113 ] [ 0.00001 pi approx ] unit-test
+
+[ 2 ] [ 1 2 approx ] unit-test
+[ 2 ] [ 0.1 2 approx ] unit-test
+[ 2 ] [ 0.00001 2 approx ] unit-test
+
+[ 3 ] [ 1 2.5 approx ] unit-test
+[ 5/2 ] [ 0.1 2.5 approx ] unit-test
+[ 5/2 ] [ 0.0001 2.5 approx ] unit-test
diff --git a/extra/math/continued-fractions/continued-fractions.factor b/extra/math/continued-fractions/continued-fractions.factor
new file mode 100644
index 0000000000..26454a3e90
--- /dev/null
+++ b/extra/math/continued-fractions/continued-fractions.factor
@@ -0,0 +1,24 @@
+! Copyright (C) 2009 Samuel Tardieu.
+! See http://factorcode.org/license.txt for BSD license.
+USING: kernel math math.functions sequences vectors ;
+IN: math.continued-fractions
+
+<PRIVATE
+
+: split-float ( f -- d i ) dup >integer [ - ] keep ;
+
+: closest ( seq -- newseq ) unclip-last round >integer suffix ;
+
+PRIVATE>
+
+: next-approx ( seq -- )
+    dup [ pop split-float ] [ push ] bi
+    dup zero? [ 2drop ] [ recip swap push ] if ;
+
+: >ratio ( seq -- a/b )
+    closest reverse unclip-slice [ swap recip + ] reduce ;
+
+: approx ( epsilon float -- a/b )
+    dup 1vector
+    [ 3dup >ratio - abs < ] [ dup next-approx ] while
+    2nip >ratio ;
diff --git a/extra/math/continued-fractions/summary.txt b/extra/math/continued-fractions/summary.txt
new file mode 100644
index 0000000000..e8b2f66654
--- /dev/null
+++ b/extra/math/continued-fractions/summary.txt
@@ -0,0 +1 @@
+Continued fractions