frac_approx/
lib.rs

1use std::ops::{Add, AddAssign, Mul, Neg};
2
3pub struct FracApproxIter<I, F> {
4    lower: Fraction<I>,
5    upper: Fraction<I>,
6    bound: Option<I>,
7    pred: F,
8}
9
10impl<I: SbInt, F: Fn(I, I) -> bool> FracApproxIter<I, F> {
11    fn new(pred: F, bound: Option<I>) -> Self {
12        let (lower, upper) = if !I::SIGNED {
13            (I::lowest_frac(), I::highest_frac())
14        } else if pred(I::ZERO, I::ONE) {
15            (I::zero_frac(), I::highest_frac())
16        } else {
17            (I::lowest_frac(), I::zero_frac())
18        };
19        let (lower, upper) = (lower.into(), upper.into());
20        Self { lower, upper, bound, pred }
21    }
22
23    fn iter_return(&mut self, q: I, tf: bool) -> (ApproxFrac<I>, I) {
24        (ApproxFrac::from(((self.lower, self.upper), tf)), q)
25    }
26}
27
28impl<I: SbInt, F: Fn(I, I) -> bool> Iterator for FracApproxIter<I, F> {
29    type Item = (ApproxFrac<I>, I);
30    fn next(&mut self) -> Option<Self::Item> {
31        let Self { lower, upper, bound, pred } = self;
32        let pred = |fr: Fraction<I>| {
33            let (lower, upper) = fr.into();
34            pred(lower, upper)
35        };
36        let mediant = *lower + *upper;
37
38        if mediant.is_deeper(*bound) {
39            return None;
40        }
41
42        let tf = pred(mediant);
43        let (from, to) = if tf { (*lower, *upper) } else { (*upper, *lower) };
44
45        let mut lo = I::ONE;
46        let mut hi = lo + I::ONE;
47        while pred(from + to * hi) == tf {
48            lo += lo;
49            hi += hi;
50            if (from + to * lo).is_deeper(*bound) {
51                let bound =
52                    if let Some(bound) = *bound { bound } else { continue };
53                let steps = bound.steps(from.into_inner(), to.into_inner());
54                if steps == I::ZERO {
55                    return None;
56                }
57                let front = from + to * steps;
58                *(if tf { lower } else { upper }) = front;
59                return Some(self.iter_return(steps, tf));
60            }
61        }
62
63        while lo.lt1(hi) {
64            let mid = lo.avg(hi);
65            let tmp = from + to * mid;
66            let cur = !tmp.is_deeper(*bound) && pred(tmp) == tf;
67            *(if cur { &mut lo } else { &mut hi }) = mid;
68        }
69
70        *(if tf { lower } else { upper }) = from + to * lo;
71        Some(self.iter_return(lo, tf))
72    }
73}
74
75pub enum ApproxFrac<I> {
76    Lower((I, I)),
77    Upper((I, I)),
78}
79use ApproxFrac::{Lower, Upper};
80
81impl<I: std::fmt::Display> std::fmt::Debug for ApproxFrac<I> {
82    fn fmt(&self, fmt: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
83        let (name, (num, den)) = match self {
84            Lower(x) => ("Lower", x),
85            Upper(x) => ("Upper", x),
86        };
87        fmt.debug_tuple(name).field(&format!("{}/{}", num, den)).finish()
88    }
89}
90
91impl<I> ApproxFrac<I> {
92    pub fn into_inner(self) -> (I, I) {
93        match self {
94            Lower(x) | Upper(x) => x,
95        }
96    }
97}
98
99impl<I: SbInt> From<((Fraction<I>, Fraction<I>), bool)> for ApproxFrac<I> {
100    fn from(((lower, upper), tf): ((Fraction<I>, Fraction<I>), bool)) -> Self {
101        if tf { Lower(lower.into()) } else { Upper(upper.into()) }
102    }
103}
104
105pub trait FracApprox<I, F> {
106    fn frac_approx_iter(self) -> FracApproxIter<I, F>;
107    fn frac_approx_iter_bound(self, bound: I) -> FracApproxIter<I, F>;
108    fn frac_approx_bound(self, bound: I) -> [(I, I); 2];
109}
110
111impl<I: SbInt, F: Fn(I, I) -> bool> FracApprox<I, F> for F {
112    fn frac_approx_iter(self) -> FracApproxIter<I, F> {
113        FracApproxIter::new(self, None)
114    }
115    fn frac_approx_iter_bound(self, bound: I) -> FracApproxIter<I, F> {
116        FracApproxIter::new(self, Some(bound))
117    }
118    fn frac_approx_bound(self, bound: I) -> [(I, I); 2] {
119        let mut lo = I::lowest_frac();
120        let mut hi = I::highest_frac();
121        for x in FracApproxIter::new(self, Some(bound)) {
122            match x {
123                (Lower(f), _) => lo = f,
124                (Upper(f), _) => hi = f,
125            }
126        }
127        [lo, hi]
128    }
129}
130
131pub trait SbInt:
132    Copy
133    + Eq
134    + PartialOrd<Self>
135    + AddAssign<Self>
136    + Add<Self, Output = Self>
137    + Mul<Self, Output = Self>
138    + std::fmt::Display
139    + std::fmt::Debug
140{
141    const ZERO: Self;
142    const ONE: Self;
143    const SIGNED: bool;
144    fn lt1(self, other: Self) -> bool;
145    fn avg(self, other: Self) -> Self;
146    fn abs(self) -> Self;
147    fn neg(self) -> Self;
148    fn steps(self, from: (Self, Self), to: (Self, Self)) -> Self;
149    fn lowest_frac() -> (Self, Self);
150    fn highest_frac() -> (Self, Self);
151    fn zero_frac() -> (Self, Self);
152}
153
154#[derive(Clone, Copy, Eq, PartialEq, std::fmt::Debug)]
155struct Fraction<I>(I, I);
156
157impl<I: SbInt> From<(I, I)> for Fraction<I> {
158    fn from((num, den): (I, I)) -> Self { Self(num, den) }
159}
160
161impl<I: SbInt> From<Fraction<I>> for (I, I) {
162    fn from(frac: Fraction<I>) -> Self { (frac.0, frac.1) }
163}
164
165impl<I: SbInt> Neg for Fraction<I> {
166    type Output = Self;
167    fn neg(self) -> Self { self.neg() }
168}
169
170impl<I: SbInt> Mul<I> for Fraction<I> {
171    type Output = Self;
172    fn mul(self, a: I) -> Self { Self(self.0 * a, self.1 * a) }
173}
174
175impl<I: SbInt> Add<Fraction<I>> for Fraction<I> {
176    type Output = Self;
177    fn add(self, other: Self) -> Self {
178        // mediant
179        Self(self.0 + other.0, self.1 + other.1)
180    }
181}
182
183impl<I: SbInt> Fraction<I> {
184    fn is_deeper(self, bound: Option<I>) -> bool {
185        bound.into_iter().any(|b| self.1 > b)
186    }
187    fn neg(self) -> Self { Self(self.0.neg(), self.1) }
188    fn into_inner(self) -> (I, I) { (self.0, self.1) }
189}
190
191macro_rules! impl_uint {
192    ( $($ty:ty)* ) => { $(
193        impl SbInt for $ty {
194            const ZERO: $ty = 0;
195            const ONE: $ty = 1;
196            const SIGNED: bool = false;
197            fn lt1(self, other: Self) -> bool { self + 1 < other }
198            fn avg(self, other: Self) -> Self {
199                self + (other - self) / 2
200            }
201            fn abs(self) -> Self { self }
202            fn neg(self) -> Self { self.wrapping_neg() } // not to be called
203            fn steps(self, from: (Self, Self), to: (Self, Self)) -> Self {
204                if to.0 == 0 {
205                    (self - from.1) / to.1
206                } else if to.1 == 0 {
207                    (self - from.0) / to.0
208                } else {
209                    ((self - from.0) / to.0).min((self - from.1) / to.1)
210                }
211            }
212            fn lowest_frac() -> (Self, Self) { (0, 1) }
213            fn highest_frac() -> (Self, Self) { (1, 0) }
214            fn zero_frac() -> (Self, Self) { (0, 1) }
215        }
216    )* }
217}
218
219impl_uint! { u8 u16 u32 u64 u128 usize }
220
221macro_rules! impl_int {
222    ( $($ty:ty)* ) => { $(
223        impl SbInt for $ty {
224            const ZERO: $ty = 0;
225            const ONE: $ty = 1;
226            const SIGNED: bool = true;
227            fn lt1(self, other: Self) -> bool { self + 1 < other }
228            fn avg(self, other: Self) -> Self {
229                self + (other - self) / 2
230            }
231            fn abs(self) -> Self { self.abs() }
232            fn neg(self) -> Self { -self}
233            fn steps(self, from: (Self, Self), to: (Self, Self)) -> Self {
234                if to.0 == 0 {
235                    (self - from.1) / to.1
236                } else if to.1 == 0 {
237                    (self - from.0) / to.0
238                } else {
239                    ((self - from.0) / to.0).min((self - from.1) / to.1)
240                }
241            }
242            fn lowest_frac() -> (Self, Self) { (-1, 0) }
243            fn highest_frac() -> (Self, Self) { (1, 0) }
244            fn zero_frac() -> (Self, Self) { (0, 1) }
245        }
246    )* }
247}
248
249impl_int! { i8 i16 i32 i64 i128 isize }
250
251#[test]
252fn sanity_check() {
253    let sqrt2 = |p: i128, q| p * p <= 2 * q * q;
254    for ap in sqrt2.frac_approx_iter().take(30) {
255        eprintln!("{ap:?}");
256    }
257
258    eprintln!("{:?}", sqrt2.frac_approx_bound(10000));
259
260    let frac = |p: i128, q| 113 * p <= 355 * q;
261    for ap in frac.frac_approx_iter_bound(113) {
262        eprintln!("{ap:?}");
263    }
264
265    const PI_INT: u128 = 3;
266    const PI_FRAC: &str = "141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647";
267    let pi = |p: u128, q| {
268        // p/q <= A/B
269        if p / q != PI_INT {
270            return p / q < PI_INT;
271        }
272        let mut r = p % q * 10;
273        for d in PI_FRAC.bytes().map(|b| (b - b'0') as u128) {
274            if r / q != d {
275                return r / q < d;
276            }
277            r = r % q * 10;
278        }
279        true
280    };
281    for ap in pi.frac_approx_iter_bound(10_u128.pow(18)) {
282        eprintln!("{ap:?}");
283    }
284}
frac_approx/lib.rs

frac_approx/
lib.rs
