BigInt to Float conversions are unnecessarily slow

[simonbyrne]

julia> using BenchmarkTools

julia> x = big"10"^300;

julia> @btime Float64($x)
  468.070 ns (6 allocations: 328 bytes)
1.0e300

julia> @btime Float64(BigFloat($x, precision=precision(BigFloat)))
  74.334 ns (2 allocations: 112 bytes)
1.0e300

This could probably be done even more efficiently by doing bit-shifting on the limbs.

[JeffBezanson]

julia> @btime Float64($x);
  287.156 ns (6 allocations: 328 bytes)

julia> @btime Float64($x, RoundToZero);
  8.383 ns (0 allocations: 0 bytes)

Yikes!

[narendrakpatel]

I am interested in solving this issue. I experimented with Float64(::UInt128) to create efficient Float64(::BigInt). But I was not able to implement RoundNearest rounding mode. Can you provide me with some leads? @simonbyrne @JeffBezanson

[simonbyrne]

Great!

It's a little bit tricky (which is why I suspect GMP hasn't implemented it), but one approach is:

0. Do some preprocessing: check if the integer is too big (so would be an Inf), or if it can be converted directly via a UInt64.
1. find the position of the highest-order bit
2. load the 54 highest-order bits into a UInt64 (Float64 has 53 bits of precision, we add an extra one for rounding).
3. adjust the last bit to do round-to-nearest, ties-to-even. What you would need here is check if the last two bits of the UInt64 are 01, and change them to 00 if all the remaining bits of the BigInt are 0.
4. Convert the UInt64 to a Float64 and adjust the exponent appropriately (can use ldexp for this).

[narendrakpatel]

I came up with this implementation

function proposed_Float64(x::BigInt)
           x == 0.0 && return 0.0
           if abs(x.size) > 16
               temp = Inf
           elseif abs(x.size) == 1
               temp = Float64(unsafe_load(x.d))
           else
               temp = Float64((UInt128(unsafe_load(x.d,abs(x.size))) << 64)
                                       +unsafe_load(x.d, (abs(x.size)-1)))
               temp *= 2.0^((abs(x.size)-2)*64)
           end
           signbit(x.size) && return -temp
           return temp
       end

By using current Float64(::BigInt) I noticed that it does not check for all the remaining bits of the BigInt. Thus this method is correct(assuming current method is correct).

julia> d = big"2"^300 + 64
2037035976334486086268445688409378161051468393665936250636140449354381299763336706183397440

julia> @btime proposed_Float64(d)
  26.645 ns (0 allocations: 0 bytes)
2.037035976334486e90

julia> @btime Float64(d)
  281.984 ns (7 allocations: 168 bytes)
2.037035976334486e90

Any suggestions? @simonbyrne

[simonbyrne]

There are a couple of issues:

1. on some platforms the limbs are only 32-bits, so you may need 3 of them to get the significant bits.
2. this can suffer from double-rounding:

julia> x = big(2)^256 + big(2)^(256-53) + 1
115792089237316208279075339080610112188966723394941384217081534270255812050945

julia> Float64(x)
1.1579208923731622e77

julia> proposed_Float64(x)
1.157920892373162e77

[narendrakpatel]

This implementation solves the double rounding problem. I still have to clean the code and implement for 32 bits. I just wanted to know, if this type of approach is okay or not? Also in 32 bit implementation should I just use if/else or is there any better way? @simonbyrne

function Float64(n::BigInt)
    n == 0 && return 0.0
    if Int == Int64
        if abs(x.size) > 16
            temp = Inf
        elseif abs(x.size) == 1
            temp = Float64(unsafe_load(x.d))
        else
            y1 = unsafe_load(x.d, x.size)
            n = 64 - leading_zeros(y1)
            if n > 53
                y = (y1 >> (n-54))
                y = (y + 1) >> 1
                y &= ~UInt64(trailing_zeros(x) == (n-54 + (x.size-1)*64))
                d = (n+1021) << 52
                temp = reinterpret(Float64, d+y)
                temp = ldexp(temp, (x.size - 1)*64)
            else
                y = y1 << (54 - n) + unsafe_load(x.d, x.size-1) >> (10 + n)
                y = (y + 1) >> 1
                y &= ~UInt64(trailing_zeros(x) == (10+n + (x.size-2)*64))
                d = (n+1021) << 52
                temp = reinterpret(Float64, d+y)
                temp = ldexp(temp, (x.size - 1)*64)
            end
        end
        signbit(x.size) && return -temp
        return temp
    end

    #TODO: implement for 32 bit machine
end

[simonbyrne]

Something along those lines was what I had in mind. Did you want to open a draft pull request? It might be easier to provide feedback on that.