ezno_checker/utilities/
modulo_class.rs

1type BetterF64 = ordered_float::NotNan<f64>;
2
3/// x ≡ *class* [mod *modulo*]
4#[derive(Debug, Clone, Copy, PartialEq)]
5pub struct ModuloClass {
6	pub modulo: BetterF64,
7	pub offset: BetterF64,
8}
9
10// TODO more operations
11impl ModuloClass {
12	#[must_use]
13	pub fn new(modulo: BetterF64, offset: BetterF64) -> Self {
14		debug_assert!(modulo != 0.);
15		if modulo > 0f64.try_into().unwrap() {
16			Self { offset: offset % modulo, modulo }
17		} else {
18			// TODO think this is correct. [1]_{-3} = [2]_{3}
19			let modulo = -modulo;
20			Self { offset: modulo - (offset % modulo), modulo }
21		}
22	}
23
24	#[must_use]
25	pub fn contains(self, value: BetterF64) -> bool {
26		// Note -0. = 0.
27		(value - self.offset) % self.modulo == 0.
28	}
29
30	/// WIP
31	#[must_use]
32	pub fn disjoint(self, other: Self) -> bool {
33		if let Ok(gcd) = gcd_of_float(self.modulo, other.modulo) {
34			crate::utilities::notify!("{:?}", gcd);
35			(self.offset - other.offset) % gcd != 0.
36		} else {
37			crate::utilities::notify!("Here");
38			true
39		}
40	}
41
42	#[must_use]
43	pub fn intersection(self, _other: Self) -> Option<Self> {
44		todo!()
45	}
46
47	#[must_use]
48	pub fn get_cover(self, _other: Self) -> Option<Self> {
49		todo!()
50	}
51
52	#[must_use]
53	pub fn offset(self, offset: BetterF64) -> Self {
54		// TODO temp fix
55		if self.is_default() {
56			self
57		} else {
58			Self::new(self.modulo, self.offset + offset)
59		}
60	}
61
62	#[must_use]
63	pub fn multiply(self, multiple: BetterF64) -> Self {
64		// TODO temp fix
65		if self.is_default() {
66			self
67		} else {
68			Self::new(self.modulo * multiple, self.offset)
69		}
70	}
71
72	#[must_use]
73	pub fn negate(self) -> Self {
74		// TODO temp fix
75		if self.is_default() {
76			self
77		} else {
78			Self::new(self.modulo, self.modulo - self.offset)
79		}
80	}
81
82	#[must_use]
83	pub fn is_default(self) -> bool {
84		self.modulo == f64::EPSILON
85	}
86}
87
88/// Using Farey algoirthm
89/// TODO is there a faster way implemntation
90/// Note that numerator and denominator are coprime
91fn try_get_numerator_denominator(input: BetterF64) -> Result<(i32, i32), ()> {
92	const STEPS: usize = 50;
93	const MARGIN: f64 = 1e-4;
94
95	let integer_part = input.trunc() as i32;
96	let fractional_part = input.fract();
97	if fractional_part == 0. {
98		return Ok((integer_part, 1));
99	}
100
101	let (mut a, mut b, mut c, mut d) = (0, 1, 1, 1);
102
103	for _ in 0..STEPS {
104		let mediant_float = (f64::from(a) + f64::from(b)) / (f64::from(c) + f64::from(d));
105		if (fractional_part - mediant_float).abs() < MARGIN {
106			let numerator = a + b + integer_part * (c + d);
107			let denominator = c + d;
108			return Ok((numerator, denominator));
109		} else if fractional_part > mediant_float {
110			a += b;
111			c += d;
112		} else {
113			b += a;
114			d += c;
115		}
116	}
117
118	Err(())
119}
120
121fn gcd_of_float(n1: BetterF64, n2: BetterF64) -> Result<BetterF64, ()> {
122	fn gcd(mut n1: i32, mut n2: i32) -> i32 {
123		while n2 != 0 {
124			let t = n2;
125			n2 = n1 % n2;
126			n1 = t;
127		}
128		n1
129	}
130
131	/// n1*n2 = gcd(n1, n2)*lcm(n1, n2)
132	fn lcm(n1: i32, n2: i32) -> i32 {
133		(n1 * n2) / gcd(n1, n2)
134	}
135
136	let (a, b) = try_get_numerator_denominator(n1)?;
137	let (c, d) = try_get_numerator_denominator(n2)?;
138
139	// gcd(a / b, c / d) = gcd(a, c) / lcm(b, d)
140	Ok(BetterF64::new(f64::from(gcd(a, c)) / f64::from(lcm(b, d))).unwrap())
141}
142
143// hmmm
144impl Default for ModuloClass {
145	fn default() -> Self {
146		Self { modulo: f64::EPSILON.try_into().unwrap(), offset: BetterF64::new(0.).unwrap() }
147	}
148}
149
150#[cfg(test)]
151mod tests {
152	use super::*;
153
154	// TODO test negatives etc
155	#[test]
156	fn gcd() {
157		assert_eq!(
158			gcd_of_float(BetterF64::new(1. / 3.).unwrap(), BetterF64::new(3. / 2.).unwrap()),
159			Ok(BetterF64::new(0.16666666666666666).unwrap())
160		);
161	}
162}
ezno_checker/utilities/modulo_class.rs

ezno_checker/utilities/
modulo_class.rs
