Parsing of hex-float and bit-float strings

binary-parse.jl

# converts bit-float and hex-float strings to floating point numbers
# assumes round-to-nearest

import Base: significand_mask, significand_bits, exponent_bias, tryparse

for B in (2,16)
    @eval function tryparse{T<:Float64}(::Type{T}, s::String, ::Type{Val{$B}})
        f_neg = false
        r_ind = false # radix pt. indicator
        r = 0 # no. of places after radix pt.
        B_bits = $B == 16 ? 4 : 1 # bits per base unit
        
        u = reinterpret(Unsigned, zero(T)) # significand
        u_max = (significand_mask(Float64)+1)<<1

        i = start(s)
        done(s,i) && return Nullable{T}()
        c,i = next(s,i)
        
        # 0: ignore preceding white-space
        while isspace(c)
            done(s,i) && return Nullable{T}()
            c,i = next(s,i)
        end
        
        # 1: sign (optional)
        if c=='+' || c=='-'
            f_neg = (c=='-')
            done(s,i) && return Nullable{T}()
            c,i = next(s,i)
        end

        # 2: remove "0x"/"0b" prefix (optional)
        if c == '0'
            done(s,i) && return Nullable(zero(T))
            cc, ii = next(s,i)
            if $B==16 ? (cc=='x' || cc=='X') : (cc=='b' || cc=='B')
                done(s,ii) && return Nullable{T}()
                c,i = next(s,ii)
            end
        end

        # require at least 1 digit before exponent
        if c == '.'
            done(s,i) && return Nullable(zero(T))
            cc, ii = next(s,i)
        else
            cc = c
        end        
        if $B == 16
            '0' <= cc <= '9' || 'a' <= cc <= 'f' || 'A' <= cc <= 'F' || return Nullable{T}()
        else
            cc == '0' || cc == '1' || return Nullable{T}()
        end
        
        # 3: significand digits
        while true
            if c == '.'
                r_ind && return Nullable{T}() # more than 1 radix point
                r_ind = true
            elseif c == 'p' || c == 'P'
                done(s,i) && return Nullable{T}()
                c,i = next(s,i)
                break
            else
                if $B == 16
                    if '0' <= c <= '9'
                        n = c - '0'
                    elseif 'a' <= c <= 'f'
                        n = c - 'a' + 10
                    elseif 'A' <= c <= 'F'
                        n = c - 'A' + 10
                    else
                        return Nullable{T}()
                    end
                else
                    if c == '0'
                        n = 0
                    elseif c == '1'
                        n = 1
                    else
                        return Nullable{T}()
                    end
                end                
                
                if u < u_max
                    u = (u << B_bits) | n
                    r += r_ind
                else
                    # excess bits: use "round-to-odd"
                    u |= (n!=0)
                    r -= !r_ind
                end
            end
            if done(s,i)
                p = 0 # no exponent
                @goto finalise
            end
            c,i = next(s,i)
        end

        # 4: exponent decimal digits
        # 4a: exponent sign
        p_neg = false
        if c=='+' || c=='-'
            p_neg = (c=='-')
            done(s,i) && return Nullable{T}()
            c,i = next(s,i)
        end

        # 4b: exponent digits    
        '0' <= c <= '9'|| return Nullable{T}() # at least 1 digit
        p = 0
        while true
            if '0' <= c <= '9'
                if p < 214748364 # prevent overflow of extremely large exponents
                    n = c-'0'
                    p = p*10 + n
                end
            else
                break
            end
            done(s,i) && @goto exponent_sign
            c,i = next(s,i)
        end
        
        # 4c: discard remaining whitespace
        while true
            isspace(c) || return Nullable{T}()
            done(s,i) && break
            c,i = next(s,i)
        end

        # 4d: set exponent sign
        @label exponent_sign
        if p_neg
            p = -p
        end

        # 5: create result
        @label finalise
        k = p - B_bits*r
        while k < -exponent_bias(T)-significand_bits(T) && u >= u_max
            # avoid double rounding of subnormals
            # keep (at most) 1 extra bit beyond max precision
            #  - if u >= u_max, result is normal: T(u) rounded, ldexp exact
            #  - if u < u_max, result is subnormal: T(u) exact, ldexp rounded
            u = (u>>1) | (u&1)
            k -= 1
        end
        f = ldexp(T(u), k)
        if f_neg
            f = -f
        end
        Nullable(f)
    end
end