diff --git a/src/libextra/num/rational.rs b/src/libextra/num/rational.rs index 1aa03350305fc..89ea5d84b7d5e 100644 --- a/src/libextra/num/rational.rs +++ b/src/libextra/num/rational.rs @@ -14,7 +14,7 @@ use std::cmp; use std::from_str::FromStr; use std::num::{Zero,One,ToStrRadix,FromStrRadix,Round}; -use super::bigint::BigInt; +use super::bigint::{BigInt, BigUint, Sign, Plus, Minus}; /// Represents the ratio between 2 numbers. #[deriving(Clone)] @@ -107,6 +107,66 @@ impl } } +impl Ratio { + // Ported from + // http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/real/misc/cl_R_rationalize.cc;hb=HEAD + // http://sourceforge.net/p/sbcl/sbcl/ci/master/tree/src/code/float.lisp#l850 + // (2013-11-17) + /// Converts a float into a rational number + fn rationalize(f: T) -> Option { + if !f.is_finite() { + return None; + } + let (mantissa, exponent, sign) = f.integer_decode(); + if mantissa == 0 || exponent >= 0 { + let mut numer: BigUint = FromPrimitive::from_u64(mantissa).unwrap(); + numer = numer << (exponent as uint); + let bigint_sign: Sign = if sign == 1 { Plus } else { Minus }; + return Some(Ratio::from_integer(BigInt::from_biguint(bigint_sign, numer))); + } + + let one: BigInt = One::one(); + let mut a: BigRational = Ratio::new( + FromPrimitive::from_u64(2 * mantissa - 1).unwrap(), + one << ((1 - exponent) as uint)); + let mut b: BigRational = Ratio::new( + FromPrimitive::from_u64(2 * mantissa + 1).unwrap(), + one << ((1 - exponent) as uint)); + + let mut p0: BigInt = Zero::zero(); + let mut p1: BigInt = One::one(); + let mut q0: BigInt = One::one(); + let mut q1: BigInt = Zero::zero(); + let mut c; + + loop { + c = a.ceil(); + if c >= b { + let k = c.to_integer() - One::one(); + + let p2 = k * p1 + p0; + p0 = p1; + p1 = p2; + let q2 = k * q1 + q0; + q0 = q1; + q1 = q2; + + let tmp = (b - Ratio::from_integer(k.clone())).recip(); + b = (a - Ratio::from_integer(k)).recip(); + a = tmp; + } else { + break; + } + }; + + let p_last: BigInt = c.to_integer() * p1 + p0; + let q_last: BigInt = c.to_integer() * q1 + q0; + let p = if sign == 1 { p_last } else { -p_last }; + Some(Ratio::new(p, q_last)) + } +} + + /* Comparisons */ // comparing a/b and c/d is the same as comparing a*d and b*c, so we @@ -621,4 +681,43 @@ mod test { test(s); } } + + #[test] + fn test_rationalize() { + fn test(given: T, expected: (&str, &str)) { + let ratio: BigRational = Ratio::rationalize(given).unwrap(); + let (numer, denom) = expected; + assert_eq!(ratio, Ratio::new( + FromStr::from_str(numer).unwrap(), + FromStr::from_str(denom).unwrap())); + } + + // f32 + test(3.14159265359f32, ("93343", "29712")); + test(2f32.pow(&100.), ("1267650600228229401496703205376", "1")); + test(-2f32.pow(&100.), ("-1267650600228229401496703205376", "1")); + test(1.0 / 2f32.pow(&100.), ("1", "1267650524670370179181738721296")); + test(684729.48391f32, ("1369459", "2")); + test(-8573.5918555f32, ("-420106", "49")); + + // f64 + test(3.14159265359f64, ("226883371", "72219220")); + test(2f64.pow(&100.), ("1267650600228229401496703205376", "1")); + test(-2f64.pow(&100.), ("-1267650600228229401496703205376", "1")); + test(1.0 / 2f64.pow(&100.), ("1", "1267650600228229260759214850049")); + test(684729.48391f64, ("68472948391", "100000")); + test(-8573.5918555f64, ("-13743853556", "1603045")); + } + + #[test] + fn test_rationalize_fail() { + use std::{f32, f64}; + + assert_eq!(Ratio::rationalize(f32::NAN), None); + assert_eq!(Ratio::rationalize(f32::INFINITY), None); + assert_eq!(Ratio::rationalize(f32::NEG_INFINITY), None); + assert_eq!(Ratio::rationalize(f64::NAN), None); + assert_eq!(Ratio::rationalize(f64::INFINITY), None); + assert_eq!(Ratio::rationalize(f64::NEG_INFINITY), None); + } }