floating_bar/
r64_t.rs

1use core::fmt;
2use core::cmp::Ordering;
3use core::ops::*;
4use core::str::FromStr;
5
6use gcd::Gcd;
7
8use super::{ParseRatioErr, r32};
9
10/// The 64-bit floating bar type.
11#[allow(non_camel_case_types)]
12#[derive(Clone, Copy, Eq, Default)]
13pub struct r64(u64);
14
15const DSIZE_SIZE: u32 = 6;
16const FRACTION_SIZE: u64 = 58;
17
18const FRACTION_FIELD: u64 = (1 << FRACTION_SIZE) - 1;
19
20impl r64 {
21	// PRIVATE methods
22	#[inline]
23	fn get_frac_size(n: i128, d: u128) -> u64 {
24		let dsize = 128 - d.leading_zeros() - 1;
25		let nsize = if n >= 0 {
26			128 - n.leading_zeros() + 1
27		} else {
28			128 - n.leading_ones() + 1
29		};
30		
31		(nsize + dsize) as u64
32	}
33	
34	// PUBLIC API
35	
36	/// Creates a rational number without checking the values.
37	/// 
38	/// # Safety
39	/// 
40	/// The values must fit in the fraction field.
41	#[inline]
42	pub const unsafe fn new_unchecked(numer: i64, denom: u64) -> r64 {
43		let denom_size = 64 - denom.leading_zeros() - 1;
44		let denom_mask = (1 << denom_size as u64) - 1;
45		let numer_mask = FRACTION_FIELD & !denom_mask;
46		
47		r64(
48			(denom_size as u64) << FRACTION_SIZE |
49			((numer << denom_size) as u64) & numer_mask |
50			denom & denom_mask
51		)
52	}
53	
54	/// Creates a rational number if the given values both fit in the fraction
55	/// field.
56	pub fn new(mut numer: i64, mut denom: u64) -> Option<r64> {
57		let gcd = numer.unsigned_abs().gcd(denom);
58		numer /= gcd as i64;
59		denom /= gcd;
60		
61		let denom_size = 64 - denom.leading_zeros() - 1;
62		let numer_size = if numer >= 0 {
63			64 - numer.leading_zeros() + 1
64		} else {
65			64 - numer.leading_ones() + 1
66		};
67		
68		if numer_size + denom_size > FRACTION_SIZE as u32 {
69			return None;
70		}
71		
72		// SAFETY: we just checked if the values fit.
73		unsafe {
74			Some(r64::new_unchecked(numer, denom))
75		}
76	}
77	
78	/// Calculates the approximate square root of the value.
79	/// 
80	/// **Warning**: This method can give a number that overflows easily, so
81	/// use it with caution, and discard it as soon as you're done with it.
82	pub fn sqrt(self) -> r64 {
83		// shh...
84		let f: f64 = self.into();
85		r64::from(f.sqrt())
86	}
87	
88	/// Calculates the approximate cube root of the value.
89	/// 
90	/// **Warning**: This method can give a number that overflows easily, so
91	/// use it with caution, and discard it as soon as you're done with it.
92	pub fn cbrt(self) -> r64 {
93		// shh...
94		let f: f64 = self.into();
95		r64::from(f.cbrt())
96	}
97	
98	/// Checked rational addition. Computes `self + rhs`, returning `None` if
99	/// overflow occurred.
100	pub fn checked_add(self, rhs: r64) -> Option<r64> {
101		if self.is_nan() || rhs.is_nan() { return Some(r64::NAN) }
102		// self = a/b, other = c/d
103		
104		// num = ad + bc
105		let mut num =
106			self.numer() as i128 * rhs.denom() as i128
107			+ self.denom() as i128 * rhs.numer() as i128;
108		// den = bd
109		let mut den = self.denom() as u128 * rhs.denom() as u128;
110		
111		let mut size = r64::get_frac_size(num, den);
112		
113		if size > FRACTION_SIZE {
114			let gcd = num.unsigned_abs().gcd(den);
115			num /= gcd as i128;
116			den /= gcd;
117			size = r64::get_frac_size(num, den);
118		}
119		
120		if size <= FRACTION_SIZE {
121			unsafe { Some(r64::new_unchecked(num as i64, den as u64)) }
122		} else {
123			None
124		}
125	}
126	
127	/// Checked rational multiplication. Computes `self * rhs`, returning `None`
128	/// if overflow occurred.
129	pub fn checked_mul(self, rhs: r64) -> Option<r64> {
130		match (self.is_nan(), rhs.is_nan()) {
131			(true, false) if rhs.numer() == 0 => return Some(r64(0)),
132			(false, true) if self.numer() == 0 => return Some(r64(0)),
133			(false, false) => {}
134			_ => return Some(r64::NAN),
135		}
136		
137		// a/b * c/d = ac/bd
138		let mut n = self.numer() as i128 * rhs.numer() as i128;
139		let mut d = self.denom() as u128 * rhs.denom() as u128;
140		
141		let mut size = r64::get_frac_size(n, d);
142		
143		if size > FRACTION_SIZE {
144			let gcd = n.unsigned_abs().gcd(d);
145			n /= gcd as i128;
146			d /= gcd;
147			size = r64::get_frac_size(n, d);
148		}
149		
150		if size <= FRACTION_SIZE {
151			unsafe { Some(r64::new_unchecked(n as i64, d as u64)) }
152		} else {
153			None
154		}
155	}
156	
157	/// Checked conversion to `i64`.
158	///
159	/// Returns `i64` value if denominator is 1. Otherwise, returns `None`.
160	#[inline]
161	pub fn to_i64(self) -> Option<i64> {
162		let norm = self.normalize();
163		if norm.denom_size() == 0 {
164			Some(norm.numer())
165		} else {
166			None
167		}
168	}
169}
170
171crate::impl_ratio_type! { r64 u64 i64 NonZeroU64 }
172
173impl From<u16> for r64 {
174	#[inline]
175	fn from(v: u16) -> Self { r64(v as u64) }
176}
177
178impl From<i16> for r64 {
179	fn from(v: i16) -> Self {
180		unsafe { r64::new_unchecked(v as i64, 1) }
181	}
182}
183
184impl From<u32> for r64 {
185	#[inline]
186	fn from(v: u32) -> Self { r64(v as u64) }
187}
188
189impl From<i32> for r64 {
190	fn from(v: i32) -> Self {
191		unsafe { r64::new_unchecked(v as i64, 1) }
192	}
193}
194
195impl From<f32> for r64 {
196	fn from(f: f32) -> Self { r64::from(f as f64) }
197}
198
199impl From<f64> for r64 {
200	/// Based on [John D. Cook's Best Rational Approximation post](https://www.johndcook.com/blog/2010/10/20/best-rational-approximation/)
201	fn from(mut f: f64) -> Self {
202		// why 29? bc it's fraction_size / 2 + 1
203		// div by 2 is to have enough space for both numer and denom.
204		// plus 1 is to count implicit bit bc numer and denom can both have 29
205		// bits of precision here.
206		const N: u64 = (1 << 29) - 1; // 2^29 - 1 = 536870911
207		let is_neg = f < 0.0;
208		
209		if f.is_nan() || f.is_infinite() {
210			return r64::NAN;
211		}
212		
213		if is_neg { f = f.abs() }
214		
215		let (mut a, mut b) = (0, 1); // lower
216		let (mut c, mut d) = (1, 0); // upper
217		let mut is_mediant = false;
218		
219		// while neither denominator is too big,
220		while b <= N && d <= N {
221			let mediant = (a + c) as f64 / (b + d) as f64;
222			
223			if f == mediant {
224				is_mediant = true;
225				break;
226			} else if f > mediant {
227				a += c;
228				b += d;
229			} else {
230				c += a;
231				d += b;
232			}
233		}
234		
235		let result = if is_mediant {
236			// if N can contain sum of both denoms,
237			if b + d <= N { (a + c, b + d) } // use sum of numers & sum of denoms
238			// else if upper bound denom is bigger than lower bound denom,
239			else if d > b { (c, d) } // use upper bound
240			else          { (a, b) } // else, use lower bound
241		} else {
242			// if lower bound denom is too big,
243			if b > N { (c, d) } // use upper bound
244			else     { (a, b) } // else, lower bound
245		};
246		
247		// SAFETY: values were given a maximum in which they could fit.
248		unsafe {
249			if is_neg {
250				r64::new_unchecked(-result.0, result.1)
251			} else {
252				r64::new_unchecked(result.0, result.1)
253			}
254		}
255	}
256}
257
258impl From<r32> for r64 {
259	fn from(v: r32) -> Self {
260		unsafe { r64::new_unchecked(v.numer() as i64, v.denom() as u64) }
261	}
262}
263
264impl From<r64> for f32 {
265	fn from(r: r64) -> f32 {
266		if r.is_nan() {
267			f32::NAN
268		} else {
269			(r.numer() as f32) / (r.denom() as f32)
270		}
271	}
272}
273
274impl From<r64> for f64 {
275	fn from(r: r64) -> f64 {
276		if r.is_nan() {
277			f64::NAN
278		} else {
279			(r.numer() as f64) / (r.denom() as f64)
280		}
281	}
282}
283
284impl PartialOrd for r64 {
285	fn partial_cmp(&self, other: &r64) -> Option<Ordering> {
286		// both are nan or both are zero
287		if self.is_nan() && other.is_nan() {
288			return Some(Ordering::Equal);
289		}
290		
291		// only one of them is nan
292		if self.is_nan() || other.is_nan() {
293			return None;
294		}
295		
296		// compare signs
297		self.is_positive()
298		.partial_cmp(&other.is_positive())
299		.map(|c| c.then(
300			// compare numbers
301			// a/b = c/d <=> ad = bc
302			// when a, b, c, and d are all > 0
303			(self.numer() as i128 * other.denom() as i128)
304			.cmp(&(self.denom() as i128 * other.numer() as i128))
305		))
306	}
307}
308
309
310#[cfg(test)]
311mod tests {
312	#[cfg(feature = "bench")]
313	extern crate test;
314	
315	use super::*;
316	use super::super::r64;
317	
318	#[test]
319	fn normalize() {
320		assert_eq!(r64!(4/2).normalize(), r64!(2));
321		assert_eq!(r64!(-4/2).normalize(), r64!(-2));
322	}
323	
324	#[test]
325	fn neg() {
326		assert_eq!(-r64(0), r64(0));
327		assert_eq!(-r64(1), r64!(-1));
328		assert_eq!(-r64!(-1), r64(1));
329	}
330	
331	#[test]
332	fn signum() {
333		assert_eq!(r64::NAN.signum(), r64::NAN);
334		assert_eq!(r64(0).signum(), r64(0));
335		assert_eq!(r64(1).signum(), r64(1));
336		assert_eq!(r64(2).signum(), r64(1));
337		assert_eq!(r64!(-1).signum(), r64!(-1));
338		assert_eq!(r64!(-2).signum(), r64!(-1));
339	}
340
341	#[test]
342	fn pow() {
343		assert_eq!(r64::NAN.pow(0), r64(1));
344		
345		assert_eq!(r64(0).pow(0),   r64(1));
346		assert_eq!(r64(1).pow(1),   r64(1));
347		
348		assert_eq!(r64(3).pow(2),    r64(9));
349		assert_eq!(r64(3).pow(-2),   r64!(1/9));
350		assert_eq!(r64!(-3).pow(2),  r64(9));
351		assert_eq!(r64!(-3).pow(-2), r64!(1/9));
352		
353		assert_eq!(r64(2).pow(3),     r64(8));
354		assert_eq!(r64(2).pow(-3),    r64!(1/8));
355		assert_eq!(r64!(1/2).pow(3),  r64!(1/8));
356		assert_eq!(r64!(1/2).pow(-3), r64(8));
357		
358		assert_eq!(r64!(-2/1).pow(3),  r64!(-8/1));
359		assert_eq!(r64!(-2/1).pow(-3), r64!(-1/8));
360		assert_eq!(r64!(-1/2).pow(3),  r64!(-1/8));
361		assert_eq!(r64!(-1/2).pow(-3), r64!(-8/1));
362	}
363
364	#[test]
365	fn checked_pow() {
366		assert_eq!(r64::NAN.checked_pow(0), Some(r64(1)));
367		assert_eq!(r64(3).checked_pow(60), None);
368	}
369
370	#[test]
371	#[cfg(feature = "roots")]
372	fn checked_sqrt() {
373		assert_eq!(r64(0).checked_sqrt(), Some(r64(0)));
374		assert_eq!(r64(1).checked_sqrt(), Some(r64(1)));
375		assert_eq!(r64(2).checked_sqrt(), None);
376		assert_eq!(r64(4).checked_sqrt(), Some(r64(2)));
377	}
378	
379	#[test]
380	fn fract() {
381		assert_eq!(r64(5).fract(),     r64(0));
382		assert_eq!(r64!(3/2).fract(),  r64!( 1/2));
383		assert_eq!(r64!(-3/2).fract(), r64!(-1/2));
384	}
385	
386	#[test]
387	fn trunc() {
388		assert_eq!(r64(5).trunc(),     r64(5));
389		assert_eq!(r64!( 1/2).trunc(), r64(0));
390		assert_eq!(r64!(-1/2).trunc(), r64(0));
391		assert_eq!(r64!( 3/2).trunc(), r64(1));
392		assert_eq!(r64!(-3/2).trunc(), r64!(-1));
393	}
394
395	#[test]
396	fn floor() {
397		assert_eq!(r64!( 3/2).floor(), r64(1));
398		assert_eq!(r64!( 2/1).floor(), r64(2));
399		assert_eq!(r64!(-3/2).floor(), r64!(-2));
400		assert_eq!(r64!(-2/1).floor(), r64!(-2));
401	}
402
403	#[test]
404	fn ceil() {
405		assert_eq!(r64!( 3/2).ceil(), r64(2));
406		assert_eq!(r64!( 2/1).ceil(), r64(2));
407		assert_eq!(r64!(-3/2).ceil(), r64!(-1));
408		assert_eq!(r64!(-2/1).ceil(), r64!(-2));
409	}
410
411	#[test]
412	fn round() {
413		assert_eq!(r64(1).round(),     r64(1));
414		assert_eq!(r64!(-1).round(),   r64!(-1));
415		assert_eq!(r64!( 3/2).round(), r64(2));
416		assert_eq!(r64!(-3/2).round(), r64!(-2));
417	}
418	
419	#[test]
420	fn min() {
421		assert_eq!(r64::NAN.min(r64::NAN), r64::NAN);
422		assert_eq!(r64::NAN.min(r64(0)),   r64(0));
423		assert_eq!(r64(0).min(r64::NAN),   r64(0));
424		assert_eq!(r64(0).min(r64(1)),     r64(0));
425	}
426	
427	#[test]
428	fn max() {
429		assert_eq!(r64::NAN.max(r64::NAN), r64::NAN);
430		assert_eq!(r64::NAN.max(r64(0)),   r64(0));
431		assert_eq!(r64(0).max(r64::NAN),   r64(0));
432		assert_eq!(r64(0).max(r64(1)),     r64(1));
433	}
434	
435	#[test]
436	fn recip() {
437		assert_eq!(r64(5).recip(),    r64!(1/5));
438		assert_eq!(r64!(5/2).recip(), r64!(2/5));
439		assert_eq!(r64(1).recip(),    r64(1));
440	}
441	
442	#[test]
443	fn cmp() {
444		assert!(r64(0) == r64(0));
445		
446		assert!(r64(0) < r64(1));
447		assert!(r64(2) < r64(3));
448		assert!(r64(0) > -r64(1));
449		assert!(r64(2) > -r64(3));
450		
451		// TODO add more assertions here
452	}
453	
454	#[test]
455	fn mul() {
456		assert_eq!(r64(0) * r64(0), r64(0));
457		
458		assert_eq!(r64(0) * r64(1), r64(0));
459		assert_eq!(r64(1) * r64(0), r64(0));
460		assert_eq!(r64(1) * r64(1), r64(1));
461		
462		assert_eq!(-r64(1) *  r64(1), -r64(1));
463		assert_eq!( r64(1) * -r64(1), -r64(1));
464		assert_eq!(-r64(1) * -r64(1),  r64(1));
465		
466		assert_eq!(r64(1) * r64(2), r64(2));
467		assert_eq!(r64(2) * r64(2), r64(4));
468		
469		assert_eq!(
470			r64!(1/2) * r64!(1/2), r64!(1/4)
471		);
472		assert_eq!(
473			r64!(-1/2) * r64!(1/2), r64!(-1/4)
474		);
475		assert_eq!(
476			r64!(2/3) * r64!(2/3), r64!(4/9)
477		);
478		assert_eq!(
479			r64!(3/2) * r64!(2/3), r64(1)
480		);
481	}
482	
483	#[test] #[should_panic]
484	fn mul_invalid() {
485		let _ = r64(1 << FRACTION_SIZE - 1) * r64(1 << FRACTION_SIZE - 1);
486	}
487	
488	#[test]
489	fn div() {
490		assert_eq!(r64(0) / r64(1), r64(0));
491		assert_eq!(r64(0) / r64(2), r64(0));
492		assert_eq!(r64(1) / r64(1), r64(1));
493		
494		assert_eq!(-r64(1) / r64(1), -r64(1));
495		assert_eq!(r64(1) / -r64(1), -r64(1));
496		assert_eq!(-r64(1) / -r64(1), r64(1));
497		
498		assert_eq!(r64(1) / r64(2), r64!(1/2));
499		assert_eq!(r64(2) / r64(1), r64(2));
500		assert_eq!(r64(2) / r64(2), r64(1));
501	}
502
503	#[test]
504	fn rem() {
505		assert_eq!(r64(5) % r64(2), r64(1));
506		assert_eq!(r64(6) % r64(2), r64(0));
507		assert_eq!(r64(8) % r64!(3/2), r64!(1/2));
508		
509		assert_eq!(-r64(5) %  r64(2),  r64(1));
510		assert_eq!( r64(5) % -r64(2), -r64(1));
511		assert_eq!(-r64(5) % -r64(2), -r64(1));
512	}
513	
514	#[test]
515	fn add() {
516		assert_eq!(r64(0) + r64(0), r64(0));
517		assert_eq!(-r64(0) + r64(0), r64(0));
518		
519		assert_eq!( r64(1) +  r64(1),  r64(2));
520		assert_eq!( r64(1) + -r64(1),  r64(0));
521		assert_eq!(-r64(1) +  r64(1),  r64(0));
522		assert_eq!(-r64(1) + -r64(1), -r64(2));
523		
524		assert_eq!(r64(2) + r64(2), r64(4));
525		assert_eq!(
526			r64!(1/2) + r64!(3/4), r64!(5/4)
527		);
528		assert_eq!(
529			r64!(1/2) + r64!(-3/4), r64!(-1/4)
530		);
531		assert_eq!(
532			r64!(-1/2) + r64!(3/4), r64!(1/4)
533		);
534	}
535	
536	#[test] #[should_panic]
537	fn add_invalid() {
538		let _ = r64(1 << FRACTION_SIZE - 1) + r64(1 << FRACTION_SIZE - 1);
539	}
540	
541	#[test]
542	fn from_str() {
543		assert_eq!("0".parse::<r64>().unwrap(),   r64(0));
544		assert_eq!("1".parse::<r64>().unwrap(),   r64(1));
545		assert_eq!("+1".parse::<r64>().unwrap(),  r64(1));
546		assert_eq!("-1".parse::<r64>().unwrap(),  r64!(-1));
547		assert_eq!("1/1".parse::<r64>().unwrap(), r64(1));
548	}
549	
550	#[test] #[should_panic]
551	fn from_str_fail() {
552		"1/-1".parse::<r64>().unwrap();
553		"/1".parse::<r64>().unwrap();
554		"1/".parse::<r64>().unwrap();
555		"1/0".parse::<r64>().unwrap();
556	}
557	
558	#[test]
559	fn from_f32() {
560		//assert_eq!(r64::from(std::f32::consts::E), r64(2850325) / r64(1048576));
561		//assert_eq!(r64::from(std::f32::consts::TAU), r64(13176795) / r64(2097152));
562	}
563	/*
564	#[test]
565	fn from_f64() {
566		assert_eq!(r64::from(0.0), r64(0));
567		assert_eq!(r64::from(1.0), r64(1));
568		assert_eq!(r64::from(-1.0), -r64(1));
569		assert_eq!(r64::from(0.2), r64(1) / r64(5));
570		assert_eq!(r64::from(1.0 - 0.7), r64(3) / r64(10));
571		//assert_eq!(r64::from(std::f64::consts::E), r64(268876667) / r64(98914198));
572		//assert_eq!(r64::from(std::f64::consts::TAU), r64(411557987) / r64(65501488));
573	}
574	*/
575}