Add float to Ratio conversion

[vmx]

Add a function to create a Ratio out of a float
(f32 or f64).

I'm a newbie, hence this commit is meant as a first step to get something similar committed. I originally tried to port the as_integer_ratio() function from CPython, but then found out that the IronPython's implementation is even simpler.

So please let me know what needs to be done to make that a committable change.

[bill-myers]

It seems this is nice to have but needs to work very differently:

• It should extract the exponent from the float and definitely not use a loop
• It should handle positive exponents
• It should handle denormals
• It should handle infinities and NaNs and thus return Option<Ratio<BigInt>> or maybe ExtendedNumber<Ratio<BigInt>> with enum ExtendedNumber<T> {NaN, Infinity(bool), Finite(T)}.
• It should call Ratio::new_raw after normalizing internally by calling trailing_zeros() on the mantissa instead of using gcd in Ratio::new() since that's faster when the gcd is always a power of 2

Also, it might make sense to preserve the sign of floating point zeros by making Ratio::reduce no longer force the denominator to be positive, but I'm not sure if this causes issues elsewhere.

[brson]

Returning Option to handle infinity an NaN is probably the least we should do. @bill-myers I don't know maths but are your other points mostly about performance? Is it otherwise correct?

[huonw]

Some googling suggests that an algorithm like http://www.math.niu.edu/~rusin/known-math/95/rationalize might be better. It uses continued fractions, which would imply that it's going to be close to the closest rational to the input float for a given upper bound on the denominator; I would suggest that allowing specifying an optional upper bound would be good too, either via max_denom: Option<&BigInt> or a separate method.

[vmx]

I'm sorry for not making any progress on this, I was just too busy with work. I'll try to get something done next week. Thanks @huonw for the link, I'll have a look.

[vmx]

I did a bit of research last week about how other programming languages deal with float to ratio conversion.

I let them ran with 0.1 and 3.14159265359:

The results were:
Python, Ruby, Java, Haskell, Go:

0.1: 3602879701896397/36028797018963968
3.14159265359: 3537118876014453/1125899906842624

Lisp (SBCL), GMP (through Python binding):

0.1: 1/10
3.14159265359: 226883371/72219220

Note that Lisp has the distinction between rational and rationalize. The difference
is the performance. rational returns the same result as the other programming languages,
where rationalize is slower but returns the same as GMP.

I currently learn towards following the Lisp route and implementing a rational as
well as a rationalize function for float.

I'm sorry that I don't remember the name, but I'd like to thank the person that gave me the idea of looking into the Lisp implementation (it was at the last Rust meetup in Mountain View).

[catamorphism]

@vmx What's the status of this -- are you stuck on anything?

[vmx]

@catamorphism Thanks for asking. I'm not really stuck. I've implemented the algorithm from Lisp in Python to get a better understanding. Next step is porting it to Rust. I'd expect to have a version by the end of the week (rebased on top of the newest master of course).

[vmx]

I've create a Gist with an implementation of rationalize. It seems to work. Tests and proper integration with libextra is still missing, but I thought I post my progress. Once done I'll add it properly to this pull request.

[huonw]

That looks cool, although the shifts seems to be prone to overflow, e.g. changing let f = <pi> to let f = let f = 2f64.pow(&100.); gives 0/1, presumably because mantissa << exponent overflows (this is actually undefined behaviour for LLVM too); I think it should do something like from_i64(sign as i64 * mantissa as i64) << exponent (or whatever casting on needs to do to get that to typecheck). Similarly for the 1 << (1 - exponent) ones below.

Also, integer_decode_float seems to depend on IEEE754 floats; I don't know if Rust (and/or LLVM) mandates this. (It would be very reasonable to.)

[vmx]

I haven't found the time to integrate it into Rust yet, but I updated my Gist according to the comments from @huonw.

According to the Rust manual are f32 and f64 IEEE 754-2008 floats, so integer_decode_float can be used.

[vmx]

I did a force push to my branch as the first version wasn't really going anywhere (if you know a way to preserve the old commits of a pull request and adding a completely unrelated branch, let me know).

I dared to have the current commit based on a slightly outdated Rust version and I also only compiled/tested it with make check-stage2-extra NO_REBUILD=1 only as I still consider it WIP and expect major complaints.

I would move integer_decode_float() to libstd f32/f64, though I'm not sure if that's the way to go, please let me know.

[vmx]

Currently the code doesn't handle Infinity and NaN. That's what I'm going to fix next. Anything else?

[vmx]

Another update of my commit. It now used the already merged Float::integer_decode(). I also cleaned up the variable names a bit.

I've left the comments where the algorithm comes from as I haven't received any permission to use it yet (And I think that's needed as the code is non-trivial). Other than that I think it's ready to have a closer review.

[huonw]

Convention would be bigint_sign.

[huonw]

A power of two should definitely be an even number, and 2**-100 is representable exactly as a float, so there shouldn't be any loses. (Specifically, this should be an exact flip of the 2f64.pow(&100.) case.)

[huonw]

Oh, also, I get different numbers in Haskell (and Python agrees from a incomplete check):

> -- f32
> let nums = [3.14159265359, 2.0^100, -2.0 ^ 100, 1.0 / (2.0^100), 684729.48391, -8573.5918555] :: [Float]
> mapM_ (\n -> putStrLn $ show n ++ "  \t" ++ show (toRational n)) nums
3.1415927       13176795 % 4194304
1.2676506e30    1267650600228229401496703205376 % 1
-1.2676506e30   (-1267650600228229401496703205376) % 1
7.888609e-31    1 % 1267650600228229401496703205376
684729.5        1369459 % 2
-8573.592       (-4389679) % 512

> -- f64
> let nums = [3.14159265359, 2.0^100, -2.0 ^ 100, 1.0 / (2.0^100), 684729.48391, -8573.5918555] :: [Double]
> mapM_ (\n -> putStrLn $ show n ++ "\t" ++ show (toRational n)) nums
3.14159265359           3537118876014453 % 1125899906842624
1.2676506002282294e30   1267650600228229401496703205376 % 1
-1.2676506002282294e30  (-1267650600228229401496703205376) % 1
7.888609052210118e-31   1 % 1267650600228229401496703205376
684729.48391            367611342500051 % 536870912
-8573.5918555           (-4713381968463931) % 549755813888

(The format is <float> <numerator> % <denominator> .)

[vmx]

Python and Haskell have a different implementation. They do the simple way. The algorithm I use tries to find the smallest possible value that still matches the float. Try Lisp (Clisp or SBCL), that's what I've used to verify my numbers. Let me write a small script to show the output.

[huonw]

Ah, I see (sorry for my clisp inexperience):

> (map nil #'(lambda (n) (format t "~S ~S~%" n (rationalize n))) `(3.14159265359d0  684729.48391d0 -8573.5918555d0 ,(expt 2.0d0 -100)))
3.14159265359d0 226883371/72219220
684729.48391d0 68472948391/100000
-8573.5918555d0 -13743853556/1603045
7.888609052210118d-31 1/1267650600228229260759214850049
NIL

I think it's a little peculiar that the exact 2.0**-100 isn't exact for f64 (it appears to be exactly one epsilon off); I guess we don't have much choice?

[vmx]

Yes, it's one epsilon off. When you look at the bit representation:

>>> import struct
>>> [struct.unpack('<Q', struct.pack('<d', f))[0] for f in (2 ** -100, 1.0/1267650600228229260759214850049, 1.0/1267650600228229401496703205376)]
[4156822456062967808, 4156822456062967809, 4156822456062967808]

We have a choice in fixing the algorithm. Not sure how easy that would be as I've only ported it and not really thought through it. I'll give it a try though.

[huonw]

Assuming can you get the license stuff OK-d, I'm happy to r+ as is, and we can investigate a perfectly precise algorithm later, if necessary.

[vmx]

I've sent a mail on Monday, still waiting for a reply :) Though I'll take the chance to have a closer look at the algorithm. I'm not sure if it makes sense to have it go in with a known bug.

[vmx]

There are three reasons for merging #11125 instead of this one:

• We found a bug in the algorithms when using 2 ** -100
• I couldn't get hold of the author of the algorithm to get permission to us it in Rust
• Implement rational() function for floats #11125 is a way simpler implementation. It doesn't return as small fractions as the algorithm above, but it probably also isn't really needed. The reason to use this function is to be able to do precise calculations on numbers, so you don't really care what the ratio looks like as long it is correct (I gained this insight thanks to @bobbl).

[huonw]

I like that reasoning. (If we want, we can add simplification methods to Rational, say self.restrict_denom(n) and self.close_approx(rat), which would (respectively) find the closest rational where the denominator is no larger than n, and the/a rational with smallest denominator that is within rat of self (or just one of them).)

[alexcrichton]

Closing in favor of #11125, because it sounds like that one should be merged over this one (but feel free to reopen if I'm wrong)