From 806e4b477d1996a3840c603bb9460139fee05f43 Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Fri, 30 Jun 2023 16:05:55 -0400
Subject: [PATCH] add rational_approx()

---
 math.scad            | 28 ++++++++++++++++++++++++++++
 tests/test_math.scad | 17 +++++++++++++++++
 2 files changed, 45 insertions(+)

diff --git a/math.scad b/math.scad
index 784e42e..889bab8 100644
--- a/math.scad
+++ b/math.scad
@@ -314,8 +314,36 @@ function lcm(a,b=[]) =
         assert(len(arglist)>0, "Invalid call to lcm with empty list(s)")
         _lcmlist(arglist);
 
+// Function rational_approx()
+// Usage:
+//   pq = rational_approx(x, maxq);
+// Description:
+//   Finds the best rational approximation p/q to the number x so that q<=maxq.  Returns
+//   the result as `[p,q]`.  If the input is zero, then returns `[0,1]`.  
+// Example:
+//   pq1 = rational_approx(PI,10);        // Returns: [22,7]
+//   pq2 = rational_approx(PI,10000);     // Returns: [355, 113]
+//   pq3 = rational_approx(221/323,500);  // Returns: [13,19]
+//   pq4 = rational_approx(0,50);         // Returns: [0,1]
+function rational_approx(x, maxq, cfrac=[], p, q) =
+  let(
+       next = floor(x),
+       fracpart = x-next,
+       cfrac = [each cfrac, next],
+       pq = _cfrac_to_pq(cfrac)
+  )
+    approx(fracpart,0) ? pq 
+  : pq[1]>maxq ? [p,q]
+  : rational_approx(1/fracpart,maxq,cfrac, pq[0], pq[1]);
 
 
+// Converts a continued fraction given as a list with leading integer term
+// into a fraction in the form p / q, returning [p,q]. 
+function _cfrac_to_pq(cfrac,p=0,q=1,ind) =
+    is_undef(ind) ? _cfrac_to_pq(cfrac,p,q,len(cfrac)-1)
+  : ind==0 ? [p+q*cfrac[0], q]
+  : _cfrac_to_pq(cfrac, q, cfrac[ind]*q+p, ind-1);
+
 
 // Section: Hyperbolic Trigonometry
 
diff --git a/tests/test_math.scad b/tests/test_math.scad
index a787165..64e59ae 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -505,6 +505,23 @@ module test_lcm() {
 }
 test_lcm();
 
+module test_rational_approx()
+{
+   pq1 = rational_approx(PI,10);       // Returns: [22,7]
+   pq2 = rational_approx(PI,10000);    // Returns: [355, 113]
+   pq3 = rational_approx(221/323,500); // Returns: [13,19]
+   pq4 = rational_approx(0,50);        // Returns: [0,1]
+   assert_equal(pq1,[22,7]);
+   assert_equal(pq2,[355,113]);
+   assert_equal(pq3,[13,19]);
+   assert_equal(pq4,[0,1]);
+   assert_equal(rational_approx(-PI,10),[-22,7]);
+   assert_equal(rational_approx(7,10), [7,1]);
+}
+test_rational_approx();
+
+
+
 
 
 module test_complex(){