WIP on BigFloat replacement

BigFloats.jl

# This file is a part of Julia. License is MIT: http://julialang.org/license

module BigFloats

export
    FloatRef,
    setprecision,
    BigFlt, Flt, Flt128, Flt256, Flt512


export convert!, idiv!, div!, pow!, add!, sub!, mul!, fma!, neg!, sqrt!, ldexp!

import
    Base: (*), +, -, /, <, <=, ==, >, >=, ^, besselj, besselj0, besselj1, bessely,
        bessely0, bessely1, ceil, cmp, convert, copysign, div,
        exp, exp2, exponent, factorial, floor, fma, hypot, isinteger,
        isfinite, isinf, isnan, ldexp, log, log2, log10, max, min, mod, modf,
        nextfloat, prevfloat, promote_rule, rem, round, show,
        sum, sqrt, string, print, trunc, precision, exp10, expm1,
        gamma, lgamma, digamma, erf, erfc, zeta, eta, log1p, airyai,
        eps, signbit, sin, cos, tan, sec, csc, cot, acos, asin, atan,
        cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, atan2,
        cbrt, typemax, typemin, unsafe_trunc, realmin, realmax, rounding,
        setrounding, maxintfloat, widen, significand, frexp, tryparse

import Base.Rounding: rounding_raw, setrounding_raw

import Base.GMP: ClongMax, CulongMax, CdoubleMax, Limb, GMP_BITS_PER_LIMB

import Base.Math.lgamma_r

# Basic type and initialization definitions

type FloatRef <: AbstractFloat
    prec::Clong
    sign::Cint
    exp::Clong
    d::Ptr{Limb}
    # Not recommended for general use
    FloatRef(prec::Clong, sign::Cint, exp::Clong, d::Ptr{Void}) = new(prec, sign, exp, d)
end

const _RM = Vector{Cint}()
const _DP = Vector{Int}()

function _fill_fr(prec::Clong = Clong(256))
    _res = FloatRef(zero(Clong), zero(Cint), zero(Clong), C_NULL)
    ccall((:mpfr_init2,:libmpfr), Void, (Ref{FloatRef}, Clong), _res, prec)
    finalizer(_res, cglobal((:mpfr_clear, :libmpfr)))
    return _res
end
_blank_fr() = FloatRef(Clong(256), Cint(1), Clong(0), C_NULL)

immutable FloatRegs
    R::FloatRef
    X::FloatRef
    Y::FloatRef
    Z::FloatRef
    FloatRegs() = new(_fill_fr(),_blank_fr(),_blank_fr(),_blank_fr())
end

const _F = Vector{FloatRegs}(0)

@inline RM() = (@inbounds x = _RM[Base.Threads.threadid()]; x)
@inline DP() = (@inbounds x = _DP[Base.Threads.threadid()]; x)

FloatRef() = _fill_fr(DP())

immutable BigFlt <: AbstractFloat
    prec::Clong # store sign in high bit of prec
    exp::Clong
    d::Vector{Limb}
end

immutable Flt{N} <: AbstractFloat
    prec::Clong
    exp::Clong
    d::NTuple{N,Limb}
end

for i = (128, 256, 512)
    len = div(i+GMP_BITS_PER_LIMB-1, GMP_BITS_PER_LIMB)
    @eval typealias $(Symbol(:Flt,i)) Flt{$len}
end

typealias FixedFltTypes Union{BigFlt, Flt128, Flt256, Flt512}
typealias FltTypes Union{FloatRef, FixedFltTypes}
typealias NumTypes Union{CulongMax, ClongMax, CdoubleMax}

#widen(::Type{Float64}) = FloatRef
widen{T<:FltTypes}(::Type{T}) = T

function __init__()
    try
        # set exponent to full range by default
        set_emin!(get_emin_min())
        set_emax!(get_emax_max())
        N = Base.Threads.nthreads()
        resize!(_F, N)
        resize!(_RM, N)
        resize!(_DP, N)
        for i = 1:N
            _F[i] = FloatRegs()
            _RM[i] = 0
            _DP[i] = 256
        end
    catch ex
        Base.showerror_nostdio(ex,
            "WARNING: Error during initialization of module BigFloats")
    end
end

@inline function _copy!(dst::FloatRef, x::BigFlt)
    dst.prec = abs(x.prec)
    dst.sign = sign(x.prec)
    dst.exp  = x.exp
    dst.d = pointer(x.d)
    dst
end

@inline function BigFlt(z::FloatRef)
    len = div(z.prec+GMP_BITS_PER_LIMB-1, GMP_BITS_PER_LIMB)
    x = BigFlt(z.prec*z.sign, z.exp, Vector{Limb}(len))
    unsafe_copy!(pointer(x.d), z.d, len)
    x
end

@inline function _copy!{T<:Flt}(dst::FloatRef, x::T)
    dst.prec = abs(x.prec)
    dst.sign = sign(x.prec)
    dst.exp  = x.exp
    dst.d = pointer_from_objref(x.d)
    dst
end

function Flt(z::FloatRef)
    N = div(z.prec+GMP_BITS_PER_LIMB-1, GMP_BITS_PER_LIMB)
    Flt{N}(z.prec*z.sign, z.exp, unsafe_load(reinterpret(Ptr{NTuple{N,Limb}},z.d)))
end

# Wrap functions that take a BigFlt and return a BigFlt
@inline function _wrap1{T<:FixedFltTypes}(f::Function, x::T)
    local old = gc_enable(false)::Bool
    try
        @inbounds _f = _F[Base.Threads.threadid()]
        f(_f.R, _copy!(_f.X, x))
        T(_f.R)
    finally
        old && gc_enable(true)
    end
end

# Wrap functions that take BigFlt, something not BigFlt and return a BigFlt
@inline function _wrap2{T<:FixedFltTypes}(f::Function, x::T, y)
    local old = gc_enable(false)::Bool
    try
        @inbounds _f = _F[Base.Threads.threadid()]
        f(_f.R, _copy!(_f.X, x), y)
        T(_f.R)
    finally
        old && gc_enable(true)
    end
end

# Wrap functions that take 1 to 3 BigFlts with a rounding mode, and return a BigFlt
@inline function _wrap_rm1{T<:FixedFltTypes}(f::Function, x::T)
    local old = gc_enable(false)::Bool
    try
        id = Base.Threads.threadid()
        @inbounds _f = _F[id]
        @inbounds f(_f.R, _copy!(_f.X, x), _RM[id])
        T(_f.R)
    finally
        old && gc_enable(true)
    end
end

@inline function _wrap_rm2{T<:FixedFltTypes}(f::Function, x::T, y::T)
    local old = gc_enable(false)
    try
        id = Base.Threads.threadid()
        @inbounds _f = _F[id]
        @inbounds f(_f.R, _copy!(_f.X, x), _copy!(_f.Y, y), _RM[id])
        T(_f.R)
    finally
        old && gc_enable(true)
    end
end

@inline function _wrap_rm3{T<:FixedFltTypes}(f::Function, x::T, y::T)
    local old = gc_enable(false)
    try
        id = Base.Threads.threadid()
        @inbounds _f = _F[id]
        @inbounds f(_f.R, _copy!(_f.X, x), _copy!(_f.Y, y), _copy!(_f.Z, z), _RM[id])
        T(_f.R)
    finally
        old && gc_enable(true)
    end
end

# Wrap functions that take 1 or 2 BigFlts, and return something not a BigFlt
@inline function _wrap_c1{T<:FixedFltTypes}(f::Function, x::T)
    local old = gc_enable(false)::Bool
    try
        @inbounds _f = _F[Base.Threads.threadid()]
        f(_copy!(_f.X, x))
    finally
        old && gc_enable(true)
    end
end

@inline function _wrap_c2{T<:FixedFltTypes}(f::Function, x::T, y::T)
    local old = gc_enable(false)
    try
        @inbounds _f = _F[Base.Threads.threadid()]
        f(_copy!(_f.X, x), _copy!(_f.Y, y))
    finally
        old && gc_enable(true)
    end
end

@inline _wrap1(f::Function, x::FloatRef)    = (res = FloatRef(); f(res, x); res)
@inline _wrap2(f::Function, x::FloatRef, y) = (res = FloatRef(); f(res, x, y); res)
@inline _wrap_c1(f::Function, x::FloatRef)  = f(x)
@inline _wrap_c2(f::Function, x::FloatRef, y::FloatRef) = f(x, y)
@inline _wrap_rm1(f::Function, x::FloatRef) =
    (res = FloatRef(); f(res, x); res)
@inline _wrap_rm2(f::Function, x::FloatRef, y::FloatRef) =
    (res = FloatRef(); f(res, x, y); res)
@inline _wrap_rm3(f::Function, x::FloatRef, y::FloatRef, z::FloatRef) =
    (res = FloatRef(); f(res, x, y); res)

convert{T<:FltTypes}(::Type{T}, x::T) = x

# convert to FloatRef
for (fJ, fC) in ((:si,:Clong), (:ui,:Culong), (:d,:Float64))
    @eval begin
        convert!(res::FloatRef, x::($fC), rm::Cint = RM()) =
            ccall(($(string(:mpfr_set_,fJ)), :libmpfr), Cint,
                  (Ref{FloatRef}, ($fC), Cint), res, x, rm)
    end
end

convert!(res::FloatRef, x::BigInt, rm::Cint = RM()) =
    ccall((:mpfr_set_z, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{BigInt}, Cint), res, x, rm)

convert{T<:FltTypes}(::Type{T}, x::Integer) = T(BigInt(x))

convert(::Type{FloatRef}, x::BigInt) = (res = FloatRef(); convert!(res, x); res)
function convert{T<:FixedFltTypes}(::Type{T}, x::BigInt)
    local old = gc_enable(false)::Bool
    try
        id = Base.Threads.threadid()
        @inbounds _f = _F[id]
        convert!(_f.R, x, _RM[id])
        T(_f.R)
    finally
        old && gc_enable(true)
    end
end

convert{T<:FltTypes}(::Type{T}, x::Union{Bool,Int8,Int16,Int32}) = T(convert(Clong,x))
convert{T<:FltTypes}(::Type{T}, x::Union{UInt8,UInt16,UInt32}) = T(convert(Culong,x))
convert{T<:FltTypes}(::Type{T}, x::Union{Float16,Float32}) = T(Float64(x))
# Note: this could be optimized for BigFlt
convert{T<:FltTypes}(::Type{T}, x::Rational) = T(num(x)) / T(den(x))

function tryparse!(res::FloatRef, s::AbstractString, base::Int=0, rm::Cint = RM())
    err = ccall((:mpfr_set_str, :libmpfr), Cint,
                (Ref{FloatRef}, Cstring, Cint, Cint), res, s, base, rm)
    err == 0 ? Nullable(res) : Nullable{FloatRef}()
end

tryparse(::Type{FloatRef}, s::AbstractString) = tryparse!(FloatRef(), s)
tryparse(::Type{FloatRef}, s::AbstractString, base::Int) = tryparse!(FloatRef(), s, base)

# Make this threadsafe and protect from GC
function tryparse{T<:FixedFltTypes}(::Type{T}, s::AbstractString, base::Int=0)
    local old = gc_enable(false)::Bool
    try
        id = Base.Threads.threadid()
        @inbounds _f = _F[id]
        err = ccall((:mpfr_set_str, :libmpfr), Cint,
                    (Ref{FloatRef}, Cstring, Cint, Cint), _f.R, s, base, _RM[id])
        err == 0 ? Nullable(T(_f.R)) : Nullable{T}()
    finally
        old && gc_enable(true)
    end
end

convert(::Type{Rational}, x::FloatRef) = convert(Rational{BigInt}, x)
#convert(::Type{AbstractFloat}, x::BigInt) = FloatRef(x)

## FloatRef -> Integer
unsafe_cast(::Type{Int64}, x::FloatRef, ri::Cint) =
    ccall((:__gmpfr_mpfr_get_sj,:libmpfr), Cintmax_t,
          (Ref{FloatRef}, Cint), x, ri)
unsafe_cast(::Type{UInt64}, x::FloatRef, ri::Cint) =
    ccall((:__gmpfr_mpfr_get_uj,:libmpfr), Cuintmax_t,
          (Ref{FloatRef}, Cint), x, ri)
function unsafe_cast(::Type{BigInt}, x::FloatRef, ri::Cint)
    # actually safe, just keep naming consistent
    _res = BigInt()
    ccall((:mpfr_get_z, :libmpfr), Cint, (Ref{BigInt}, Ref{FloatRef}, Cint),
          _res, x, ri)
    _res
end

## BigFlt -> Integer
function unsafe_cast{T<:FixedFltTypes}(::Type{Int64}, x::T, ri::Cint)
    local old = gc_enable(false)::Bool
    try
        @inbounds _f = _F[Base.Threads.threadid()]
        ccall((:__gmpfr_mpfr_get_sj,:libmpfr), Cintmax_t,
              (Ref{FloatRef}, Cint), _copy!(_f.X, x), ri)
    finally
        old && gc_enable(true)
    end
end
function unsafe_cast{T<:FixedFltTypes}(::Type{UInt64}, x::T, ri::Cint)
    local old = gc_enable(false)::Bool
    try
        @inbounds _f = _F[Base.Threads.threadid()]
        ccall((:__gmpfr_mpfr_get_uj,:libmpfr), Cuintmax_t,
              (Ref{FloatRef}, Cint), _copy!(_f.X, x), ri)
    finally
        old && gc_enable(true)
    end
end
function unsafe_cast{T<:FixedFltTypes}(::Type{BigInt}, x::T, ri::Cint)
    local old = gc_enable(false)::Bool
    try
        @inbounds _f = _F[Base.Threads.threadid()]
        # actually safe, just keep naming consistent
        _res = BigInt()
        ccall((:mpfr_get_z, :libmpfr), Cint, (Ref{BigInt}, Ref{FloatRef}, Cint),
              _res, _copy!(_f.X, x), ri)
        _res
    finally
        old && gc_enable(true)
    end
end

unsafe_cast{T<:Signed}(::Type{T}, x::FltTypes, ri::Cint)   = unsafe_cast(Int64, x, ri) % T
unsafe_cast{T<:Unsigned}(::Type{T}, x::FltTypes, ri::Cint) = unsafe_cast(UInt64, x, ri) % T
unsafe_cast(::Type{Int128}, x::FltTypes, ri::Cint) = Int128(unsafe_cast(BigInt,x,ri))
unsafe_cast(::Type{UInt128}, x::FltTypes, ri::Cint) = UInt128(unsafe_cast(BigInt,x,ri))
unsafe_cast{T<:Integer}(::Type{T}, x::FltTypes, r::RoundingMode) = unsafe_cast(T,x,to_mpfr(r))
unsafe_trunc{T<:Integer}(::Type{T}, x::FltTypes) = unsafe_cast(T,x,RoundToZero)

function trunc{T<:Union{Signed,Unsigned}}(::Type{T}, x::FltTypes)
    (typemin(T) <= x <= typemax(T)) || throw(InexactError())
    unsafe_cast(T,x,RoundToZero)
end
function floor{T<:Union{Signed,Unsigned}}(::Type{T}, x::FltTypes)
    (typemin(T) <= x <= typemax(T)) || throw(InexactError())
    unsafe_cast(T,x,RoundDown)
end
function ceil{T<:Union{Signed,Unsigned}}(::Type{T}, x::FltTypes)
    (typemin(T) <= x <= typemax(T)) || throw(InexactError())
    unsafe_cast(T,x,RoundUp)
end

function round{T<:Union{Signed,Unsigned}}(::Type{T}, x::FltTypes)
    (typemin(T) <= x <= typemax(T)) || throw(InexactError())
    unsafe_cast(T,x,RM())
end

trunc(::Type{BigInt}, x::FltTypes) = unsafe_cast(BigInt, x, RoundToZero)
floor(::Type{BigInt}, x::FltTypes) = unsafe_cast(BigInt, x, RoundDown)
ceil(::Type{BigInt}, x::FltTypes) = unsafe_cast(BigInt, x, RoundUp)
round(::Type{BigInt}, x::FltTypes) = unsafe_cast(BigInt, x, RM())

# convert/round/trunc/floor/ceil(Integer, x) should return a BigInt
trunc(::Type{Integer}, x::FltTypes) = trunc(BigInt, x)
floor(::Type{Integer}, x::FltTypes) = floor(BigInt, x)
ceil(::Type{Integer}, x::FltTypes) = ceil(BigInt, x)
round(::Type{Integer}, x::FltTypes) = round(BigInt, x)

convert(::Type{Bool}, x::FltTypes) = x==0 ? false : x==1 ? true : throw(InexactError())
function convert(::Type{BigInt},x::FltTypes)
    isinteger(x) || throw(InexactError())
    trunc(BigInt,x)
end

function convert{T<:Integer}(::Type{T},x::FltTypes)
    isinteger(x) || throw(InexactError())
    trunc(T,x)
end

## FloatRef -> AbstractFloat
convert(::Type{Float64}, x::FloatRef) =
    ccall((:mpfr_get_d,:libmpfr), Float64, (Ref{FloatRef},Cint), x, RM())
convert(::Type{Float32}, x::FloatRef) =
    ccall((:mpfr_get_flt,:libmpfr), Float32, (Ref{FloatRef},Cint), x, RM())

(::Type{Float64})(x::FloatRef, r::RoundingMode) =
    ccall((:mpfr_get_d,:libmpfr), Float64, (Ref{FloatRef},Cint), x, to_mpfr(r))
(::Type{Float32})(x::FloatRef, r::RoundingMode) =
    ccall((:mpfr_get_flt,:libmpfr), Float32, (Ref{FloatRef},Cint), x, to_mpfr(r))

## BigFlt -> AbstractFloat
function convert{T<:FixedFltTypes}(::Type{Float64}, x::T)
    local old = gc_enable(false)::Bool
    try
        id = Base.Threads.threadid()
        @inbounds _f = _F[id]
        ccall((:mpfr_get_d,:libmpfr), Float64, (Ref{FloatRef}, Cint), _copy!(_f.X, x), _RM[id])
    finally
        old && gc_enable(true)
    end
end
function convert{T<:FixedFltTypes}(::Type{Float32}, x::T)
    local old = gc_enable(false)::Bool
    try
        id = Base.Threads.threadid()
        @inbounds _f = _F[id]
        ccall((:mpfr_get_flt,:libmpfr), Float32, (Ref{FloatRef}, Cint), _copy!(_f.X, x), _RM[id])
    finally
        old && gc_enable(true)
    end
end
function (::Type{Float64}){T<:FixedFltTypes}(x::T, r::RoundingMode)
    local old = gc_enable(false)::Bool
    try
        ccall((:mpfr_get_d,:libmpfr), Float64,
              (Ref{FloatRef}, Cint), _copy!(_F[Base.Threads.threadid()].X, x), to_mpfr(r))
    finally
        old && gc_enable(true)
    end
end
function (::Type{Float32}){T<:FixedFltTypes}(x::T, r::RoundingMode)
    local old = gc_enable(false)::Bool
    try
        ccall((:mpfr_get_flt,:libmpfr), Float32,
              (Ref{FloatRef}, Cint), _copy!(_F[Base.Threads.threadid()].X, x), to_mpfr(r))
    finally
        old && gc_enable(true)
    end
end


# TODO: avoid double rounding
convert(::Type{Float16}, x::FltTypes) = convert(Float16, convert(Float32, x))
# TODO: avoid double rounding
(::Type{Float16})(x::FltTypes, r::RoundingMode) =
    convert(Float16, Float32(x, r))

promote_rule{T<:Real, S<:FltTypes}(::Type{S}, ::Type{T}) = S
promote_rule{T<:AbstractFloat, S<:FltTypes}(::Type{S},::Type{T}) = S

#This conflicts with BigFloat definition for now
#promote_rule{T<:AbstractFloat}(::Type{BigInt},::Type{T}) = FloatRef

function convert(::Type{Rational{BigInt}}, x::FltTypes)
    isnan(x) && return zero(BigInt)//zero(BigInt)
    isinf(x) && return copysign(one(BigInt),x)//zero(BigInt)
    x == 0   && return zero(BigInt) // one(BigInt)
    s = max(precision(x) - exponent(x), 0)
    BigInt(ldexp(x,s)) // (BigInt(1) << s)
end

 # Basic arithmetic without promotion
for (fJ,fC) in ((:add!, :add), (:mul!, :mul))
    @eval begin
        ($fJ)(res::FloatRef, x::FloatRef, y::FloatRef, rm::Cint = RM()) =
            ccall(($(string(:mpfr_,fC)),:libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint),
                  res, x, y, rm)
    end
    for (fS, T, RT) in ((:_ui, CulongMax, Culong),
                        (:_si, ClongMax, Clong),
                        (:_d,  CdoubleMax, Cdouble),
                        (:_z,  BigInt, Ref{BigInt}))
        @eval begin
            ($fJ)(res::FloatRef, x::FloatRef, y::$T, rm::Cint = RM()) =
                ccall(($(string(:mpfr_,fC)),:libmpfr), Cint,
                      (Ref{FloatRef}, Ref{FloatRef}, $RT, Cint),
                      res, x, y, rm)
            ($fJ)(c::$T, x::FloatRef) = ($fJ)(x,c)
        end
    end
end

sub!(res::FloatRef, c::BigInt, x::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_z_sub, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{BigInt}, Ref{FloatRef}, Cint),
          res, c, x, rm)

fma!(res::FloatRef, x::FloatRef, y::FloatRef, z::FloatRef, rm::Cint = RM()) =
    ccall(("mpfr_fma",:libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint),
          res, x, y, z, rm)

# div
for (fC, TX, TY, RX, RY) in ((:div,    FloatRef, FloatRef,   Ref{FloatRef}, Ref{FloatRef}),
                             (:div_ui, FloatRef, CulongMax,  Ref{FloatRef}, Culong),
                             (:ui_div, CulongMax, FloatRef,  Culong, Ref{FloatRef}),
                             (:div_si, FloatRef, ClongMax,   Ref{FloatRef}, Clong),
                             (:si_div, ClongMax, FloatRef,   Clong, Ref{FloatRef}),
                             (:div_d,  FloatRef, CdoubleMax, Ref{FloatRef}, Cdouble),
                             (:d_div,  CdoubleMax, FloatRef, Cdouble, Ref{FloatRef}),
                             (:div_z,  FloatRef, BigInt,     Ref{FloatRef}, Ref{BigInt}))
    @eval begin
        function idiv!(res::FloatRef, x::$TX, y::$TY)
            ccall(($(string(:mpfr_,fC)),:libmpfr), Cint,
                  (Ref{FloatRef}, $RX, $RY, Cint),
                  res, x, y, to_mpfr(RoundToZero))
            ccall((:mpfr_trunc, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}), res, res)
        end
    end
end

# /
for (fC, TX, TY, RX, RY) in ((:div,    FloatRef, FloatRef,   Ref{FloatRef}, Ref{FloatRef}),
                             (:div_ui, FloatRef, CulongMax,  Ref{FloatRef}, Culong),
                             (:ui_div, CulongMax, FloatRef,  Culong, Ref{FloatRef}),
                             (:div_si, FloatRef, ClongMax,   Ref{FloatRef}, Clong),
                             (:si_div, ClongMax, FloatRef,   Clong, Ref{FloatRef}),
                             (:div_d,  FloatRef, CdoubleMax, Ref{FloatRef}, Cdouble),
                             (:d_div,  CdoubleMax, FloatRef, Cdouble, Ref{FloatRef}),
                             (:div_z,  FloatRef, BigInt,     Ref{FloatRef}, Ref{BigInt}))
    @eval begin
        function div!(res::FloatRef, x::$TX, y::$TY, rm::Cint = RM())
            ccall(($(string(:mpfr_,fC)),:libmpfr), Cint,
                  (Ref{FloatRef}, $RX, $RY, Cint),
                  res, x, y, rm)
        end
    end
end

for (fC, TX, TY, RX, RY) in ((:sub,    FloatRef, FloatRef,   Ref{FloatRef}, Ref{FloatRef}),
                             (:sub_ui, FloatRef, CulongMax,  Ref{FloatRef}, Culong),
                             (:ui_sub, CulongMax, FloatRef,  Culong, Ref{FloatRef}),
                             (:sub_si, FloatRef, ClongMax,   Ref{FloatRef}, Clong),
                             (:si_sub, ClongMax, FloatRef,   Clong, Ref{FloatRef}),
                             (:sub_d,  FloatRef, CdoubleMax, Ref{FloatRef}, Cdouble),
                             (:d_sub,  CdoubleMax, FloatRef, Cdouble, Ref{FloatRef}),
                             (:sub_z,  FloatRef, BigInt,     Ref{FloatRef}, Ref{BigInt}))
    @eval begin
        function sub!(res::FloatRef, x::$TX, y::$TY, rm::Cint = RM())
            ccall(($(string(:mpfr_,fC)),:libmpfr), Cint,
                  (Ref{FloatRef}, $RX, $RY, Cint),
                  res, x, y, rm)
        end
    end
end
    
# More efficient commutative operations
for (fJ, fC) in ((:+, :add!), (:*, :mul!))
    @eval begin
        function ($fJ)(a::FloatRef, b::FloatRef, c::FloatRef)
            rm = RM()
            res = FloatRef()
            $fC(res, a, b, rm)
            $fC(res, res, c, rm)
            return res
        end
        function ($fJ)(a::FloatRef, b::FloatRef, c::FloatRef, d::FloatRef)
            rm = RM()
            res = FloatRef()
            $fC(res, a, b, rm)
            $fC(res, res, c, rm)
            $fC(res, res, d, rm)
            return res
        end
        function ($fJ)(a::FloatRef, b::FloatRef, c::FloatRef, d::FloatRef, e::FloatRef)
            rm = RM()
            res = FloatRef()
            $fC(res, a, b, rm)
            $fC(res, res, c, rm)
            $fC(res, res, d, rm)
            $fC(res, res, e, rm)
            return res
        end
        function ($fJ){T<:FixedFltTypes}(a::T, b::T, c::T)
            local old = gc_enable(false)::Bool
            try
                id = Base.Threads.threadid()
                @inbounds _f = _F[id]
                @inbounds rm = _RM[id]
                $fC(_f.R, _copy!(_f.x, a), _copy!(_f.y, b), rm)
                $fC(_f.R, _f.R, _copy!(_f.x, c), rm)
                T(_f.R)
            finally
                old && gc_enable(true)
            end
        end
        function ($fJ){T<:FixedFltTypes}(a::T, b::T, c::T, d::T)
            local old = gc_enable(false)::Bool
            try
                id = Base.Threads.threadid()
                @inbounds _f = _F[id]
                @inbounds rm = _RM[id]
                $fC(_f.R, _copy!(_f.x, a), _copy!(_f.y, b), rm)
                $fC(_f.R, _f.R, _copy!(_f.x, c), rm)
                $fC(_f.R, _f.R, _copy!(_f.x, d), rm)
                T(_f.R)
            finally
                old && gc_enable(true)
            end
        end
        function ($fJ)(a::BigFlt, b::BigFlt, c::BigFlt, d::BigFlt, e::BigFlt)
            try
                id = Base.Threads.threadid()
                @inbounds _f = _F[id]
                @inbounds rm = _RM[id]
                $fC(_f.R, _copy!(_f.x, a), _copy!(_f.y, b), rm)
                $fC(_f.R, _f.R, _copy!(_f.x, c), rm)
                $fC(_f.R, _f.R, _copy!(_f.x, d), rm)
                $fC(_f.R, _f.R, _copy!(_f.x, e), rm)
                T(_f.R)
            finally
                old && gc_enable(true)
            end
        end
    end
end

neg!(res::FloatRef, x::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_neg, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Cint),
          res, x, rm)

function sqrt!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    isnan(x) && return x
    ccall((:mpfr_sqrt, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Cint),
          res, x, rm)
    isnan(res) && throw(DomainError())
end

-(x::FltTypes) = _wrap_rm1(neg!, x)
sqrt(x::FltTypes) = isnan(x) ? x : _wrap_rm1(sqrt!, x)

#sqrt(x::BigInt) = sqrt(FloatRef(x))

for (fC, T, RT) in ((Symbol(""), FloatRef, Ref{FloatRef}),
                    (:_ui, CulongMax, Culong),
                    (:_si, ClongMax, Clong),
                    (:_d,  CdoubleMax, Cdouble),
                    (:_z,  BigInt, Ref{BigInt}))
    @eval begin
        pow!(res::FloatRef, x::FloatRef, y::$T, rm::Cint = RM()) =
            ccall(($(string(:mpfr_pow,fC)),:libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}, $RT, Cint),
                  res, x, y, rm)
        ^(x::FltTypes, y::$T) = _wrap_rm2(pow!, x, y)
    end
end

^(x::FltTypes, y::Integer)  = typemin(Clong)  <= y <= typemax(Clong)  ? x^Clong(y)  : x^BigInt(y)
^(x::FltTypes, y::Unsigned) = typemin(Culong) <= y <= typemax(Culong) ? x^Culong(y) : x^BigInt(y)

for f in (:exp, :exp2, :exp10, :expm1, :digamma, :erf, :erfc, :zeta,
          :cosh,:sinh,:tanh,:sech,:csch,:coth, :cbrt)
    fe = Symbol(f,'!')
    @eval begin
        export $fe
        ($fe)(res::FloatRef, x::FloatRef, rm::Cint = RM()) =
            ccall(($(string(:mpfr_,f)), :libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, rm)
    end
end

# return log(2)
big_ln2!(res::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_const_log2, :libmpfr), Cint, (Ref{FloatRef}, Cint), res, rm)
big_ln2() = (res = FloatRef(); big_ln2!(res); res)

function eta(x::FloatRef)
    x == 1 ? big_ln2() : -zeta(x) * expm1(big_ln2()*(1-x))
end

ldexp!(res::FloatRef, x::FloatRef, n::Clong, rm::Cint = RM()) =
    ccall((:mpfr_mul_2si, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Clong, Cint), res, x, n, rm)
ldexp!(res::FloatRef, x::FloatRef, n::Culong, rm::Cint = RM()) =
    ccall((:mpfr_mul_2ui, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Culong, Cint), res, x, n, rm)

ldexp(x::FloatRef, n::Union{Clong, Culong}) = (res = FloatRef(); ldexp!(res, x, n); res)

ldexp(x::FltTypes, n::ClongMax) = ldexp(x, convert(Clong, n))
ldexp(x::FltTypes, n::CulongMax) = ldexp(x, convert(Culong, n))
ldexp(x::FloatRef, n::Integer) = x*exp2(FloatRef(n))
#ldexp(x::BigFlt, n::Integer) = x*exp2(FloatRef(n))

for (fJ,fC) in ((:airyai,   :ai),
                (:besselj0, :j0),
                (:besselj1, :j1),
                (:besselj,  :jn))
    fe = Symbol(fJ,'!')
    @eval begin
        export $fe
        ($fe)(res::FloatRef, x::FloatRef, rm::Cint = RM()) =
            ccall((string(:mpfr_,$fC), :libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, rm)
        ($fJ)(x::FltTypes) = _wrap_rm1($fe, x)
    end
end

for (fJ,fC) in ((:bessely0, :y0),
                (:bessely1, :y1),
                (:bessely,  :yn))
    fe = Symbol(fJ,'!')
    @eval begin
        export $fe
        function ($fe)(res::FloatRef, x::FloatRef, rm::Cint = RM())
            x < 0 && throw(DomainError())
            ccall((string(:mpfr_,$fC), :libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, rm)
        end
        ($fJ)(x::FltTypes) = _wrap_rm1($fe, x)
    end
end

function factorial!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    if x < 0 || !isinteger(x)
        throw(DomainError())
    end
    ui = convert(Culong, x)
    ccall((:mpfr_fac_ui, :libmpfr), Cint,
        (Ref{FloatRef}, Culong, Cint), res, ui, rm)
end
factorial(x::FltTypes) = _wrap_rm1(factorial!, x)

hypot!(res::FloatRef, x::FloatRef, y::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_hypot, :libmpfr), Cint,
        (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, y, rm)
hypot{T<:FltTypes}(x::T, y::T) = _wrap_rm2(hypot!, x, y)

for f in (:log, :log2, :log10)
    fe = Symbol(f,'!')
    @eval begin
        function ($fe)(res::FloatRef, x::FloatRef, rm::Cint = RM())
            x < 0 && throw(DomainError())
            ccall(($(string(:mpfr_,f)), :libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, rm)
        end
        $f(x::FltTypes) = _wrap_rm1($fe, x)
    end
end

function log1p!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    x < -1 && throw(DomainError())
    ccall((:mpfr_log1p, :libmpfr), Cint,
        (Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, rm)
end

max!(res::FloatRef, x::FloatRef, y::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_max, :libmpfr), Cint,
        (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, y, rm)

min!(res::FloatRef, x::FloatRef, y::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_min, :libmpfr), Cint,
        (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, y, rm)

function modf(x::FloatRef)
    isinf(x) && return (FloatRef(NaN), x)
    zint = FloatRef()
    zfloat = FloatRef()
    ccall((:mpfr_modf, :libmpfr), Cint,
        (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), zint, zfloat, x, RM())
    return (zfloat, zint)
end

function modf(x::BigFlt)
    isinf(x) && return (BigFlt(NaN), x)
    zint = FloatRef()
    zfloat = FloatRef()
    ccall((:mpfr_modf, :libmpfr), Cint,
        (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), zint, zfloat, x, RM())
    return (BigFlt(zfloat), BigFlt(zint))
end

rem!(res::FloatRef, x::FloatRef, y::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_fmod, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, y, rm)

function sum!(res::FloatRef, arr::AbstractArray{FloatRef})
    for val in arr
        ccall((:mpfr_add, :libmpfr), Cint,
            (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint),
            res, res, val, 0)
    end
end
                                              
# Functions for which NaN results are converted to DomainError, following Base
for f in (:sin,:cos,:tan,:sec,:csc,:acos,:asin,:atan,:acosh,:asinh,:atanh,:gamma)
    fe = Symbol(f,'!')
    @eval begin
        ($fe)(res::FloatRef, x::FloatRef, rm::Cint = RM()) =
            ccall(($(string(:mpfr_,f)), :libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, rm)
        function ($f)(x::FloatRef)
            isnan(x) && return x
            res = FloatRef()
            ($fe)(res, x)
            isnan(res) && throw(DomainError())
            res
        end
        function ($f)(x::BigFlt)
            isnan(x) && return x
            res = _wrap_rm1($fe, x)
            isnan(res) && throw(DomainError())
            res
        end
    end
end

# log of absolute value of gamma function
const lgamma_signp = Array{Cint}(1)
lgamma!(res::FloatRef, x::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_lgamma,:libmpfr), Cint,
          (Ref{FloatRef}, Ptr{Cint}, Ref{FloatRef}, Cint),
          res, lgamma_signp, x, rm)

lgamma_r(x::FloatRef) = (lgamma(x), lgamma_signp[1])

atan2!(res::FloatRef, y::FloatRef, x::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_atan2, :libmpfr), Cint,
        (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), res, y, x, rm)

# Utility functions
==(x::FloatRef, y::FloatRef) =
    ccall((:mpfr_equal_p, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}), x, y) != 0
<=(x::FloatRef, y::FloatRef) =
    ccall((:mpfr_lessequal_p, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}), x, y) != 0
>=(x::FloatRef, y::FloatRef) =
    ccall((:mpfr_greaterequal_p, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}), x, y) != 0
<(x::FloatRef, y::FloatRef) =
    ccall((:mpfr_less_p, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}), x, y) != 0
>(x::FloatRef, y::FloatRef) =
    ccall((:mpfr_greater_p, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}), x, y) != 0

function cmp(x::FloatRef, y::BigInt)
    isnan(x) && throw(DomainError())
    ccall((:mpfr_cmp_z, :libmpfr), Cint, (Ref{FloatRef}, Ref{BigInt}), x, y)
end
function cmp(x::FloatRef, y::ClongMax)
    isnan(x) && throw(DomainError())
    ccall((:mpfr_cmp_si, :libmpfr), Cint, (Ref{FloatRef}, Clong), x, y)
end
function cmp(x::FloatRef, y::CulongMax)
    isnan(x) && throw(DomainError())
    ccall((:mpfr_cmp_ui, :libmpfr), Cint, (Ref{FloatRef}, Culong), x, y)
end

cmp(x::FltTypes, y::Integer) = cmp(x,big(y))
cmp(x::Integer, y::FltTypes) = -cmp(y,x)

function cmp(x::FloatRef, y::CdoubleMax)
    (isnan(x) || isnan(y)) && throw(DomainError())
    ccall((:mpfr_cmp_d, :libmpfr), Cint, (Ref{FloatRef}, Cdouble), x, y)
end
cmp(x::CdoubleMax, y::FltTypes) = -cmp(y,x)

==(x::FltTypes, y::Integer)   = !isnan(x) && cmp(x,y) == 0
==(x::Integer, y::FltTypes)   = y == x
==(x::FltTypes, y::CdoubleMax) = !isnan(x) && !isnan(y) && cmp(x,y) == 0
==(x::CdoubleMax, y::FltTypes) = y == x

<(x::FltTypes, y::Integer)   = !isnan(x) && cmp(x,y) < 0
<(x::Integer, y::FltTypes)   = !isnan(y) && cmp(y,x) > 0
<(x::FltTypes, y::CdoubleMax) = !isnan(x) && !isnan(y) && cmp(x,y) < 0
<(x::CdoubleMax, y::FltTypes) = !isnan(x) && !isnan(y) && cmp(y,x) > 0

<=(x::FltTypes, y::Integer)   = !isnan(x) && cmp(x,y) <= 0
<=(x::Integer, y::FltTypes)   = !isnan(y) && cmp(y,x) >= 0
<=(x::FltTypes, y::CdoubleMax) = !isnan(x) && !isnan(y) && cmp(x,y) <= 0
<=(x::CdoubleMax, y::FltTypes) = !isnan(x) && !isnan(y) && cmp(y,x) >= 0

signbit(x::FloatRef) = ccall((:mpfr_signbit, :libmpfr), Cint, (Ref{FloatRef},), x) != 0
signbit(x::FltTypes) = _wrap1(signbit, x)

# precision of an object of type FloatRef
precision(x::FloatRef) = ccall((:mpfr_get_prec, :libmpfr), Clong, (Ref{FloatRef},), x)
precision(x::BigFlt) = abs(x.prec)
precision{T<:Flt}(x::T) = abs(x.prec)

precision(::Type{FloatRef}) = DP()  # precision of the type FloatRef itself

"""
    setprecision([T=FloatRef,] precision::Int)

Set the precision (in bits) to be used for `T` arithmetic.
"""
function setprecision{T<:FltTypes}(::Type{T}, precision::Int)
    precision < 2 && throw(DomainError())
    @inbounds _DP[Base.Threads.threadid()] = precision
    precision
end

setprecision(precision::Int) = setprecision(FloatRef, precision)

maxintfloat{T<:FltTypes}(x::T) = T(2)^precision(x)
maxintfloat{T<:FltTypes}(::Type{T}) = T(2)^precision(T)

to_mpfr(::RoundingMode{:Nearest}) = Cint(0)
to_mpfr(::RoundingMode{:ToZero}) = Cint(1)
to_mpfr(::RoundingMode{:Up}) = Cint(2)
to_mpfr(::RoundingMode{:Down}) = Cint(3)
to_mpfr(::RoundingMode{:FromZero}) = Cint(4)

function from_mpfr(c::Integer)
    if c == 0
        return RoundNearest
    elseif c == 1
        return RoundToZero
    elseif c == 2
        return RoundUp
    elseif c == 3
        return RoundDown
    elseif c == 4
        return RoundFromZero
    else
        throw(ArgumentError("invalid MPFR rounding mode code: $c"))
    end
    RoundingMode(c)
end

rounding_raw{T<:FltTypes}(::Type{T}) = RM()
setrounding_raw{T<:FltTypes}(::Type{T},i::Integer) = (_RM[Base.Threads.threadid()] = i)

rounding{T<:FltTypes}(::Type{T}) = from_mpfr(rounding_raw(FloatRef))
setrounding{T<:FltTypes}(::Type{T},r::RoundingMode) = setrounding_raw(FloatRef,to_mpfr(r))

copysign!(res::FloatRef, x::FloatRef, y::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_copysign, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, y, rm)

function exponent(x::FloatRef)
    if x == 0 || !isfinite(x)
        throw(DomainError())
    end
    # The '- 1' is to make it work as Base.exponent
    return ccall((:mpfr_get_exp, :libmpfr), Clong, (Ref{FloatRef},), x) - 1
end

function frexp(x::FloatRef)
    z = FloatRef()
    c = Clong[0]
    ccall((:mpfr_frexp, :libmpfr), Cint,
          (Ptr{Clong}, Ref{FloatRef}, Ref{FloatRef}, Cint), c, z, x, RM())
    return (z, c[1])
end

function significand!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    c = Clong[0]
    ccall((:mpfr_frexp, :libmpfr), Cint,
          (Ptr{Clong}, Ref{FloatRef}, Ref{FloatRef}, Cint), c, res, x, rm)
    # Double the significand to make it work as Base.significand
    ccall((:mpfr_mul_si, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Clong, Cint), res, res, 2, rm)
end

isinteger(x::FloatRef) = ccall((:mpfr_integer_p, :libmpfr), Cint, (Ref{FloatRef},), x) != 0

for f in (:ceil, :floor, :trunc)
    fe = Symbol(f,'!')
    @eval begin
        function ($fe)(res::FloatRef, x::FloatRef)
            ccall(($(string(:mpfr_,f)), :libmpfr), Cint,
                  (Ref{FloatRef}, Ref{FloatRef}), res, x)
        end
    end
end

round!(res::FloatRef, x::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_rint, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}, Cint), res, x, rm)

round!(res::FloatRef, x::FloatRef, ::RoundingMode{:NearestTiesAway}) =
    ccall((:mpfr_round, :libmpfr), Cint,
          (Ref{FloatRef}, Ref{FloatRef}), res, x)

isinf(x::FloatRef) = ccall((:mpfr_inf_p, :libmpfr), Cint, (Ref{FloatRef},), x) != 0
isnan(x::FloatRef) = ccall((:mpfr_nan_p, :libmpfr), Cint, (Ref{FloatRef},), x) != 0
isfinite(x::FloatRef) = !isinf(x) && !isnan(x)

for fC in (:exponent, :isinteger, :isinf, :isnan, :isfinite)
    @eval ($fC)(x::BigFlt) = _wrap1($fC, x)
end

function nextfloat!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    ccall((:mpfr_set, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}, Cint),
          res, x, rm)
    ccall((:mpfr_nextabove, :libmpfr), Cint, (Ref{FloatRef},), res) != 0
end

function prevfloat!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    ccall((:mpfr_set, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}, Cint),
          res, x, rm)
    ccall((:mpfr_nextbelow, :libmpfr), Cint, (Ref{FloatRef},), res) != 0
end

# Non-mutating forms that return a new floating value
sum(arr::AbstractArray{FloatRef}) = (res = FloatRef(); sum!(res, arr); res)

for fT in (:FloatRef, :BigFlt, :Flt)
    @eval round(x::$fT, ::RoundingMode{:NearestTiesAway}) =
        _wrap2(round!, x, RoundingMode{:NearestTiesAway})
    @eval eps(::Type{$fT})     = nextfloat($fT(1)) - $fT(1)
    @eval realmin(::Type{$fT}) = nextfloat(zero($fT))
    @eval realmax(::Type{$fT}) = prevfloat($fT(Inf))

    for (fJ,fC) in ((:+, :add!), (:*, :mul!), (:-, :sub!))
        @eval ($fJ){T<:Union{$fT, NumTypes}}(x::$fT, y::T) = _wrap_rm2($fC, x, y)
    end
    for (fJ,fC) in ((:+, :add!), (:*, :mul!))
        @eval ($fJ)(x::NumTypes, y::$fT) = _wrap_rm2($fC, y, x)
    end
    @eval (/){T<:Union{$fT, NumTypes}}(x::$fT, y::T) = _wrap2(div!, x, y)
    @eval div{T<:Union{$fT, NumTypes}}(x::$fT, y::T) = _wrap2(idiv!, x, y)
    # Single argument wrapping
    for f in (:exp, :exp2, :exp10, :expm1, :digamma, :erf, :erfc, :zeta,
              :cosh,:sinh,:tanh,:sech,:csch,:coth, :cbrt,
              :log1p, :lgamma, :significand, :round)
        @eval $f(x::$fT) = _wrap_rm1($(Symbol(f, '!')), x)
    end
    # Two argument wrapping
    for f in (:max, :min, :rem, :atan2, :copysign)
        @eval $f(x::$fT, y::$fT) = _wrap_rm2($(Symbol(f, '!')), x, y)
    end
    # Three argument wrapping
    @eval fma(x::$fT, y::$fT, z::$fT) = _wrap_rm3(fma!, x, y, z)
end

"""
    setprecision(f::Function, [T=FloatRef,] precision::Integer)

Change the `T` arithmetic precision (in bits) for the duration of `f`.
It is logically equivalent to:

    old = precision(FloatRef)
    setprecision(FloatRef, precision)
    f()
    setprecision(FloatRef, old)

Often used as `setprecision(T, precision) do ... end`
"""
function setprecision{T}(f::Function, ::Type{T}, prec::Integer)
    old_prec = precision(T)
    setprecision(T, prec)
    try
        return f()
    finally
        setprecision(T, old_prec)
    end
end

setprecision(f::Function, precision::Integer) = setprecision(f, FloatRef, precision)

for fT in (:FloatRef, :BigFlt, :Flt)
    @eval typemax(::Type{$fT}) = ($fT)(Inf)
    @eval typemin(::Type{$fT}) = ($fT)(-Inf)
end

function string(x::FloatRef)
    # In general, the number of decimal places needed to read back the number exactly
    # is, excluding the most significant, ceil(log(10, 2^precision(x)))
    k = ceil(Cint, precision(x) * 0.3010299956639812)
    lng = k + Cint(8) # Add space for the sign, the most significand digit, the dot and the exponent
    buf = Array{UInt8}(lng + 1)
    # format strings are guaranteed to contain no NUL, so we don't use Cstring
    lng = ccall((:mpfr_snprintf,:libmpfr), Cint,
                (Ptr{UInt8}, Culong, Ptr{UInt8}, Ref{FloatRef}...), buf, lng + 1, "%.Re", x)
    if lng < k + 5 # print at least k decimal places
        lng = ccall((:mpfr_sprintf,:libmpfr), Cint,
                    (Ptr{UInt8}, Ptr{UInt8}, Ref{FloatRef}...), buf, "%.$(k)Re", x)
    elseif lng > k + 8
        buf = Array{UInt8}(lng + 1)
        lng = ccall((:mpfr_snprintf,:libmpfr), Cint,
                    (Ptr{UInt8}, Culong, Ptr{UInt8}, Ref{FloatRef}...), buf, lng + 1, "%.Re", x)
    end
    n = (1 <= x < 10 || -10 < x <= -1 || x == 0) ? lng - 4 : lng
    return String(buf[1:n])
end
string(x::FixedFltTypes) = _wrap_c1(string, x)

print(io::IO, b::FltTypes) = print(io, string(b))
show(io::IO, b::FltTypes) = print(io, string(b))

# get/set exponent min/max
get_emax() = ccall((:mpfr_get_emax, :libmpfr), Clong, ())
get_emax_min() = ccall((:mpfr_get_emax_min, :libmpfr), Clong, ())
get_emax_max() = ccall((:mpfr_get_emax_max, :libmpfr), Clong, ())

get_emin() = ccall((:mpfr_get_emin, :libmpfr), Clong, ())
get_emin_min() = ccall((:mpfr_get_emin_min, :libmpfr), Clong, ())
get_emin_max() = ccall((:mpfr_get_emin_max, :libmpfr), Clong, ())

set_emax!(x) = ccall((:mpfr_set_emax, :libmpfr), Void, (Clong,), x)
set_emin!(x) = ccall((:mpfr_set_emin, :libmpfr), Void, (Clong,), x)

function Base.deepcopy_internal(x::FloatRef, stackdict::ObjectIdDict)
    haskey(stackdict, x) && return stackdict[x]
    z = FloatRef(zero(Clong), zero(Cint), zero(Clong), C_NULL)
    ccall((:mpfr_init2,:libmpfr), Void, (Ref{FloatRef}, Clong), z, x.prec)
    finalizer(z, cglobal((:mpfr_clear, :libmpfr)))
    ccall((:mpfr_set, :libmpfr), Cint, (Ref{FloatRef}, Ref{FloatRef}, Cint),
          z, x, RM())
    stackdict[x] = z
    z
end

Base.deepcopy_internal(x::Union{Flt,BigFlt}, stackdict::ObjectIdDict) = x

end #module