Last active
December 28, 2015 15:19
-
-
Save vmx/7521098 to your computer and use it in GitHub Desktop.
Rationalize implementation for Rust based on the Clisp one.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
extern mod extra; | |
use std::cast; | |
use std::num::{One, Zero}; | |
use extra::rational::BigRational; | |
use extra::rational::Ratio; | |
use extra::bigint::{BigInt, BigUint, Sign, Plus, Minus}; | |
// Ported from | |
// http://www.javadocexamples.com/java_source/org/armedbear/lisp/FloatFunctions.java.html | |
// (2013-11-17) | |
fn integer_decode_float(f: f64) -> (u64, i16, i8) { | |
let bits: u64 = unsafe { | |
cast::transmute(f) | |
}; | |
let sign: i8 = if bits >> 63 == 0 { 1 } else { -1 }; | |
let mut exponent: i16 = ((bits >> 52) & 0x7ff) as i16; | |
let mantissa = if exponent == 0 { | |
(bits & 0xfffffffffffff) << 1 | |
} else { | |
(bits & 0xfffffffffffff) | 0x10000000000000 | |
}; | |
// Exponent bias + mantissa shift | |
exponent -= 1023 + 52; | |
(mantissa, exponent, sign) | |
} | |
// Ported from | |
// http://www.ginac.de/CLN/cln.git/?p=cln.git;a=blob;f=src/real/misc/cl_R_rationalize.cc;h=835c4b97fc3a98427898f7e0952094816a8750dc;hb=HEAD | |
// http://sourceforge.net/p/sbcl/sbcl/ci/master/tree/src/code/float.lisp#l850 | |
// (2013-11-17) | |
fn rationalize(f: f64) -> BigRational { | |
let (mantissa, exponent, sign) = integer_decode_float(f); | |
println!("{} {} {}", sign, exponent, mantissa); | |
if mantissa == 0 || exponent >= 0 { | |
let mut numer: BigUint = FromPrimitive::from_u64(mantissa).unwrap(); | |
numer = numer << (exponent as uint); | |
let bigintSign: Sign = match sign { | |
-1 => Minus, | |
1 | _ => Plus | |
}; | |
return Ratio::from_integer(BigInt::from_biguint(bigintSign, 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)); | |
println!("a: {}", a.to_str()) | |
println!("b: {}", b.to_str()) | |
let mut p_iminus1: BigInt = Zero::zero(); | |
let mut p_i: BigInt = One::one(); | |
let mut q_iminus1: BigInt = One::one(); | |
let mut q_i: BigInt = Zero::zero(); | |
let mut c; | |
loop { | |
println!(""); | |
println!("a: {}", a.to_str()); | |
println!("b: {}", b.to_str()); | |
c = a.ceil(); | |
if c >= b { | |
let k = c.to_integer() - One::one(); | |
let p_iplus1 = k * p_i + p_iminus1; | |
p_iminus1 = p_i; | |
p_i = p_iplus1; | |
let q_iplus1 = k * q_i + q_iminus1; | |
q_iminus1 = q_i; | |
q_i = q_iplus1; | |
let tmp = (b - Ratio::from_integer(k.clone())).recip(); | |
b = (a - Ratio::from_integer(k.clone())).recip(); | |
a = tmp; | |
} else { | |
break; | |
} | |
}; | |
let p_last: BigInt = c.to_integer() * p_i + p_iminus1; | |
let q_last: BigInt = c.to_integer() * q_i + q_iminus1; | |
if sign == 1 { | |
Ratio::new(p_last, q_last) | |
} else { | |
Ratio::new(-p_last, q_last) | |
} | |
} | |
fn assert_eq(given: f64, expected: (&str, &str)) { | |
let ratio: BigRational = rationalize(given); | |
let (numer, denom) = expected; | |
assert_eq!(ratio, Ratio::new( | |
FromStr::from_str(numer).unwrap(), | |
FromStr::from_str(denom).unwrap() | |
)); | |
} | |
fn main() { | |
//let f = 3.14159265359f64; | |
//let ratio: BigRational = rationalize(f); | |
//println!("{}", ratio.to_str()); | |
assert_eq(3.14159265359f64, ("226883371", "72219220")); | |
assert_eq(2f64.pow(&100.), ("1267650600228229401496703205376", "1")); | |
assert_eq(-2f64.pow(&100.), ("-1267650600228229401496703205376", "1")); | |
assert_eq(1.0 / 2f64.pow(&100.), ("1", "1267650600228229260759214850049")); | |
} |
(Without getting permission from the copyright holders, that is)
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
You can't use the port of integer_decode_float in the stdlib, its license is incompatible. Same with rational, unfortunately.