Port limit_denominator from cpython, used while encoding phases

quizx/src/json/phase.rs

@@ -18,6 +18,7 @@
 
 use super::VertexPhase;
 use crate::graph::VType;
+use crate::scalar::utils::limit_denominator;
 
 use num::{FromPrimitive, One, Rational64, Zero};
 
@@ -37,14 +38,11 @@ impl VertexPhase {
             return Self("0".to_string());
         }
 
-        #[allow(clippy::if_same_then_else)]
-        let simstr = if *phase.denom() > 256 {
-            // TODO: This should approximate the phase to a Rational64 number with a small denominator.
-            //       This is not currently implemented.
-            //       See https://docs.python.org/3/library/fractions.html#fractions.Fraction.limit_denominator.
-            ""
+        let (phase, simstr) = if *phase.denom() > 256 {
+            let phase = limit_denominator(phase, 256);
+            (phase, "~")
         } else {
-            ""
+            (phase, "")
         };
 
         let numer = if *phase.numer() == 1 {

quizx/src/scalar.rs

@@ -14,6 +14,8 @@
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
+pub mod utils;
+
 use approx::AbsDiffEq;
 use num::complex::Complex;
 use num::rational::Rational64;

quizx/src/scalar/utils.rs

@@ -0,0 +1,80 @@
+//! Utility functions for scalar values.
+
+use num::rational::Ratio;
+use num::Integer;
+
+/// Approximate a fraction to a Rational64 number with a small denominator.
+///
+/// Mimics python's
+/// [https://docs.python.org/3/library/fractions.html#fractions.Fraction.limit_denominator].
+///
+/// # Panics
+///
+/// Panics if `max_denom` is 0.
+///
+/// # Example
+/// ```
+/// # use quizx::scalar::utils::limit_denominator;
+/// # use num::FromPrimitive;
+/// # use num::rational::Rational64;
+/// assert_eq!(
+///     limit_denominator(Rational64::from_f64(3.141592653589793).unwrap(), 10),
+///     Rational64::new(22, 7)
+/// );
+/// assert_eq!(
+///     limit_denominator(Rational64::from_f64(3.141592653589793).unwrap(), 100),
+///     Rational64::new(311, 99)
+/// );
+/// assert_eq!(
+///     limit_denominator(Rational64::new(4321, 8765), 10000),
+///     Rational64::new(4321, 8765)
+/// );
+/// ```
+pub fn limit_denominator<T>(fraction: Ratio<T>, max_denom: T) -> Ratio<T>
+where
+    T: Clone + Integer,
+{
+    // Using cpython's limit_denominator as a reference.
+    //
+    // [https://github.com/python/cpython/blob/bfc57d43d8766120ba0c8f3f6d7b2ac681a81d8a/Lib/fractions.py#L357]
+    if max_denom <= T::one() {
+        panic!("max_denom must be greater than 1");
+    }
+
+    let mut numer = fraction.numer().clone();
+    let mut denom = fraction.denom().clone();
+    if denom <= max_denom {
+        return fraction;
+    }
+
+    // Upper and lower bounds for the approximation.
+    let mut numer_0 = T::zero();
+    let mut denom_0 = T::one();
+    let mut numer_1 = T::one();
+    let mut denom_1 = T::zero();
+    loop {
+        let a = numer.div_floor(&denom);
+        let new_numer = denom_0.clone() + a.clone() * denom_1.clone();
+        if new_numer > max_denom {
+            break;
+        }
+        let new_denom = numer_0 + a.clone() * numer_1.clone();
+        numer_0 = numer_1;
+        denom_0 = denom_1;
+        numer_1 = new_denom;
+        denom_1 = new_numer;
+
+        let tmp_denom = numer - a * denom.clone();
+        numer = denom;
+        denom = tmp_denom;
+    }
+    let k = (max_denom - denom_0.clone()).div_floor(&denom_1);
+
+    if (T::one() + T::one()) * denom * (denom_0.clone() + k.clone() * denom_1.clone())
+        <= *fraction.denom()
+    {
+        Ratio::new_raw(numer_1, denom_1)
+    } else {
+        Ratio::new_raw(numer_0 + k.clone() * numer_1, denom_0 + k * denom_1)
+    }
+}