Latest (WIP) version of BigFlt (to replace BigFloat) module, with Flt type

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!, prec!

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

"""
    A FloatRef is a mutable container of a MPFR arbitrary precision floating point number

    It is designed so that one can efficiently accumulate results, can specify the desired
    precision and rounding mode for the result (if applicable).
    It also uses memory assigned by Julia, so that the overhead of setting a finalizer is not
    required.
"""
type FloatRef <: AbstractFloat
    prec::Clong
    sign::Cint
    exp::Clong
    d::Ptr{Limb}
    vec::Vector{Limb}	# This isn't part of the MPFR type
    # Not recommended for general use
    FloatRef(prec::Clong, sign::Cint, exp::Clong, v::Vector{Limb}) =
        new(prec, sign, exp, pointer(v), v)
end

typealias FloatRefPnt Ptr{Void}

Base.unsafe_convert(::Type{FloatRefPnt}, x::FloatRef) = pointer_from_objref(x)
num_limbs(prec) = div(prec+GMP_BITS_PER_LIMB-1,GMP_BITS_PER_LIMB)
limbs_bits(n)   = n*GMP_BITS_PER_LIMB

const MaxPrec  = 512
const MaxLimbs = num_limbs(MaxPrec)
const _RM = Vector{Cint}()
const _DP = Vector{Int}()
const emptyLimbs = Vector{Limb}()

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

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

const FloatExpZer = typemin(Clong)+1
const FloatExpNaN = typemin(Clong)+2
const FloatExpInf = typemin(Clong)+3

_fill_fr(prec::Clong = Clong(MaxPrec)) =
    FloatRef(prec, Cint(1), FloatExpNaN, Vector{Limb}(num_limbs(prec)))

_blank_fr() = FloatRef(Clong(0), Cint(1), FloatExpNaN, emptyLimbs)

immutable FloatRegs
    # These are real FloatRefs, to hold results
    R::FloatRef
    S::FloatRef
    # These are filled in using a pointer to the Vector{Limb} of a BigFlt
    X::FloatRef
    Y::FloatRef
    Z::FloatRef
    # These are used for Flt{N} operations
    FX::FloatRef
    FY::FloatRef
    FZ::FloatRef
end

FloatRegs() =
    FloatRegs(_fill_fr(),_fill_fr(),
              _blank_fr(),_blank_fr(),_blank_fr(),
              _fill_fr(),_fill_fr(),_fill_fr())

const _F = Vector{FloatRegs}(0)

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

"""
Create a new FloatRef with the default precision
"""
FloatRef() = _fill_fr(DP())

for i = (128, 256, 512)
    @eval typealias $(Symbol(:Flt,i)) Flt{$(num_limbs(i))}
end
#const FltBitTypes = (Flt128, Flt256, Flt512)

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

#widen(::Type{Float64}) = Flt128
widen(::Type{Flt128}) = Flt256
widen(::Type{Flt256}) = Flt512

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

for ind = 1:3
    s = (:X,:Y,:Z)[ind]
    FX = Symbol('F',s)
    fS = Symbol("_copy",s,'!')
    @eval begin
        @inline function ($fS)(f::FloatRegs, x::BigFlt)
            r = f.$s
            r.prec = abs(x.prec)
            r.sign = sign(x.prec)
            r.exp  = x.exp
            r.vec = x.vec
            r.d = pointer(r.vec)
            r
        end
        @inline function ($fS){N}(f::FloatRegs, x::Flt{N})
            r = f.$FX
            r.prec = abs(x.prec)
            r.sign = sign(x.prec)
            r.exp  = x.exp
            unsafe_copy!(r.d, reinterpret(Ptr{Limb},pointer_from_objref(x.vec)), N)
            r
        end
    end
end

@inline function BigFlt(z::FloatRef)
    len = num_limbs(z.prec)
    x = BigFlt(z.prec*z.sign, z.exp, Vector{Limb}(len))
    unsafe_copy!(pointer(x.vec), z.d, len)
    x
end

@inline function Flt(z::FloatRef)
    N = num_limbs(z.prec)
    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)
    @inbounds _f = _F[Base.Threads.threadid()]
    _f.R.prec = abs(x.prec)
    f(_f.R, _copyX!(_f, x))
    T(_f.R)
end

for t1 in (FloatRef, BigFlt, Flt, Flt128, Flt256, Flt512),
    st in (Int8,Int16,Int32,Int64,Bool,UInt8,UInt16,UInt32,UInt64)
    @eval promote_rule(::Type{$t1}, ::Type{$st}  ) = $t1
end

promote_rule(::Type{Flt128}, ::Type{Float16}) = Flt128
promote_rule(::Type{Flt128}, ::Type{Float32}) = Flt128
promote_rule(::Type{Flt128}, ::Type{Float64}) = Flt128
promote_rule(::Type{Flt256}, ::Type{Flt128})  = Flt256
promote_rule(::Type{Flt512}, ::Type{Flt256})  = Flt512

# Wrap functions that take BigFlt, something not BigFlt and return a BigFlt
@inline function _wrap2{T<:FixedFltTypes}(f::Function, x::T, y)
    @inbounds _f = _F[Base.Threads.threadid()]
    _f.R.prec = abs(x.prec)
    f(_f.R, _copyX!(_f, x), y)
    T(_f.R)
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)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    _f.R.prec = abs(x.prec)
    @inbounds f(_f.R, _copyX!(_f, x), _RM[id])
    T(_f.R)
end

@inline function _wrap_rm2{T<:FixedFltTypes}(f::Function, x::T, y::T)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    _f.R.prec = abs(x.prec)
    @inbounds f(_f.R, _copyX!(_f, x), _copyY!(_f, y), _RM[id])
    T(_f.R)
end

@inline function _wrap_rm3{T<:FixedFltTypes}(f::Function, x::T, y::T)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    _f.R.prec = abs(x.prec)
    @inbounds f(_f.R, _copyX!(_f, x), _copyY!(_f, y), _copyZ!(_f, z), _RM[id])
    T(_f.R)
end

# Wrap functions that take 1 or 2 BigFlts, and return something not a BigFlt
@inline _wrap_c1{T<:FixedFltTypes}(f::Function, x::T) =
    f(_copyX!(_F[Base.Threads.threadid()], x))
@inline function _wrap_c2{T<:FixedFltTypes}(f::Function, x::T, y::T)
    @inbounds _f = _F[Base.Threads.threadid()]
    f(_copyX!(_f, x), _copyY!(_f, y))
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, z); res)

# Wrap functions that take 1 or 2 BigFlt or Flts, optionally with a rounding mode, return a BigFlt
@inline function _wrap_res1{T<:Union{BigFlt,Flt}}(f::Function, res::FloatRef, x::T)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    @inbounds f(res, _copyX!(_f, x), _RM[id])
end
@inline function _wrap_res1{T<:Union{BigFlt,Flt}}(f::Function, res::FloatRef, x::T, rm::Cint)
    @inbounds _f = _F[Base.Threads.threadid()]
    @inbounds f(res, _copyX!(_f, x), rm)
end
@inline function _wrap_res2{T<:Union{BigFlt,Flt}}(f::Function, res::FloatRef, x::T, y::T)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    @inbounds f(res, _copyX!(_f, x), _RM[id])
end
@inline function _wrap_res2{T<:Union{BigFlt,Flt}}(f::Function,res::FloatRef,x::T,y::T,rm::Cint)
    @inbounds _f = _F[Base.Threads.threadid()]
    @inbounds f(res, _copyX!(_f, x), rm)
end

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,
                  (FloatRefPnt, ($fC), Cint), pointer_from_objref(res), x, rm)
        convert(::Type{FloatRef}, x::$fC) = (res = FloatRef(); convert!(res, x); res)
    end
end

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

function convert{N}(::Type{FloatRef}, x::Flt{N})
    r = _fill_fr(limbs_bits(N))
    r.sign = sign(x.prec)
    r.exp  = x.exp
    unsafe_copy!(r.d, pointer_from_objref(x.d), N)
    r
end

function convert(::Type{FloatRef}, x::BigFlt)
    r = _fill_fr(abs(x.prec))
    r.sign = sign(x.prec)
    r.exp  = x.exp
    unsafe_copy!(r.d, pointer(x.d), length(x.d))
    r
end

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)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    @inbounds _f.R.prec = _DP[id]
    convert!(_f.R, x, _RM[id])
    T(_f.R)
end

function convert{T<:FixedFltTypes}(::Type{T}, x::Float64)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    @inbounds _f.R.prec = _DP[id]
    convert!(_f.R, x, _RM[id])
    T(_f.R)
end

function convert{N}(::Type{Flt{N}}, x::FloatRef)
    limbs_bits(N) == x.prec &&
        return Flt{N}(x.prec*x.sign, x.exp, unsafe_load(reinterpret(Ptr{NTuple{N,Limb}},x.d)))
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
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,
                (FloatRefPnt, 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)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    err = ccall((:mpfr_set_str, :libmpfr), Cint,
                (FloatRefPnt, Cstring, Cint, Cint), _f.R, s, base, _RM[id])
    err == 0 ? Nullable(T(_f.R)) : Nullable{T}()
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,
          (FloatRefPnt, Cint), x, ri)
unsafe_cast(::Type{UInt64}, x::FloatRef, ri::Cint) =
    ccall((:__gmpfr_mpfr_get_uj,:libmpfr), Cuintmax_t,
          (FloatRefPnt, 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}, FloatRefPnt, Cint),
          _res, x, ri)
    _res
end

## BigFlt -> Integer
function unsafe_cast{T<:FixedFltTypes}(::Type{Int64}, x::T, ri::Cint)
    @inbounds _f = _F[Base.Threads.threadid()]
    ccall((:__gmpfr_mpfr_get_sj,:libmpfr), Cintmax_t,
          (Ptr{Void}, Cint), _copyX!(_f, x), ri)
end
function unsafe_cast{T<:FixedFltTypes}(::Type{UInt64}, x::T, ri::Cint)
    @inbounds _f = _F[Base.Threads.threadid()]
    ccall((:__gmpfr_mpfr_get_uj,:libmpfr), Cuintmax_t,
          (FloatRefPnt, Cint), _copyX!(_f, x), ri)
end
function unsafe_cast{T<:FixedFltTypes}(::Type{BigInt}, x::T, ri::Cint)
    @inbounds _f = _F[Base.Threads.threadid()]
    # actually safe, just keep naming consistent
    _res = BigInt()
    ccall((:mpfr_get_z, :libmpfr), Cint, (Ref{BigInt}, FloatRefPnt, Cint),
          _res, _copyX!(_f, x), ri)
    _res
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, (FloatRefPnt,Cint), x, RM())
convert(::Type{Float32}, x::FloatRef) =
    ccall((:mpfr_get_flt,:libmpfr), Float32, (FloatRefPnt,Cint), x, RM())

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

## BigFlt -> AbstractFloat
function convert{T<:FixedFltTypes}(::Type{Float64}, x::T)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    ccall((:mpfr_get_d,:libmpfr), Float64, (FloatRefPnt, Cint), _copyX!(_f, x), _RM[id])
end
function convert{T<:FixedFltTypes}(::Type{Float32}, x::T)
    id = Base.Threads.threadid()
    @inbounds _f = _F[id]
    ccall((:mpfr_get_flt,:libmpfr), Float32, (FloatRefPnt, Cint), _copyX!(_f, x), _RM[id])
end
(::Type{Float64}){T<:FixedFltTypes}(x::T, r::RoundingMode) =
    ccall((:mpfr_get_d,:libmpfr), Float64,
          (FloatRefPnt, Cint), _copyX!(_f, _F[Base.Threads.threadid()], x), to_mpfr(r))
(::Type{Float32}){T<:FixedFltTypes}(x::T, r::RoundingMode) =
    ccall((:mpfr_get_flt,:libmpfr), Float32,
          (FloatRefPnt, Cint), _copyX!(_f, _F[Base.Threads.threadid()], x), to_mpfr(r))

# 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,
                  (FloatRefPnt, FloatRefPnt, FloatRefPnt, 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,fS)),:libmpfr), Cint,
                      (FloatRefPnt, FloatRefPnt, $RT, Cint),
                      res, x, y, rm)
            ($fJ)(res::FloatRef, c::$T, x::FloatRef) = ($fJ)(res, x, c)
            ($fJ)(res::FloatRef, c::$T, x::FloatRef, rm) = ($fJ)(res, x, c, rm)
        end
    end
end

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

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

for (fC, TX, TY, RX, RY) in ((:sub,    FloatRef, FloatRef,   FloatRefPnt, FloatRefPnt),
                             (:sub_ui, FloatRef, CulongMax,  FloatRefPnt, Culong),
                             (:ui_sub, CulongMax, FloatRef,  Culong, FloatRefPnt),
                             (:sub_si, FloatRef, ClongMax,   FloatRefPnt, Clong),
                             (:si_sub, ClongMax, FloatRef,   Clong, FloatRefPnt),
                             (:sub_d,  FloatRef, CdoubleMax, FloatRefPnt, Cdouble),
                             (:d_sub,  CdoubleMax, FloatRef, Cdouble, FloatRefPnt),
                             (:sub_z,  FloatRef, BigInt,     FloatRefPnt, Ref{BigInt}))
    @eval begin
        function sub!(res::FloatRef, x::$TX, y::$TY, rm::Cint = RM())
            ccall(($(string(:mpfr_,fC)),:libmpfr), Cint,
                  (FloatRefPnt, $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)
            @inbounds _f = _F[id]
            @inbounds rm = _RM[id]
            $fC(_f.R, _copyX!(_f, a), _copyY!(_f, b), rm)
            $fC(_f.R, _f.R, _copyX!(_f, c), rm)
            T(_f.R)
        end
        function ($fJ){T<:FixedFltTypes}(a::T, b::T, c::T, d::T)
            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, _copyX!(_f, c), rm)
            $fC(_f.R, _f.R, _copyX!(_f, d), rm)
            T(_f.R)
        end
        function ($fJ)(a::BigFlt, b::BigFlt, c::BigFlt, d::BigFlt, e::BigFlt)
            id = Base.Threads.threadid()
            @inbounds _f = _F[id]
            @inbounds rm = _RM[id]
            $fC(_f.R, _copyX!(_f, a), _copyY!(_f, b), rm)
            $fC(_f.R, _f.R, _copyX!(_f, c), rm)
            $fC(_f.R, _f.R, _copyX!(_f, d), rm)
            $fC(_f.R, _f.R, _copyX!(_f, e), rm)
            T(_f.R)
        end
    end
end

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

fma!(res::FloatRef, x::FloatRef, y::FloatRef, z::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_fma, :libmpfr), Cint,
          (FloatRefPnt, FloatRefPnt, FloatRefPnt, FloatRefPnt, Cint),
          res, x, y, z, rm)

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

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

for (fC, T, RT) in ((Symbol(""), FloatRef, FloatRefPnt),
                    (:_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,
                  (FloatRefPnt, FloatRefPnt, $RT, Cint),
                  res, x, y, rm)
    end
end

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,
                  (FloatRefPnt, FloatRefPnt, Cint), res, x, rm)
    end
end

# return log(2)
big_ln2!(res::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_const_log2, :libmpfr), Cint, (FloatRefPnt, 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,
          (FloatRefPnt, FloatRefPnt, Clong, Cint), res, x, n, rm)

ldexp!(res::FloatRef, x::FloatRef, n::Culong, rm::Cint = RM()) =
    ccall((:mpfr_mul_2ui, :libmpfr), Cint,
          (FloatRefPnt, FloatRefPnt, 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,
                  (FloatRefPnt, FloatRefPnt, 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,
                  (FloatRefPnt, FloatRefPnt, 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,
          (FloatRefPnt, Culong, Cint), res, ui, rm)
end

hypot!(res::FloatRef, x::FloatRef, y::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_hypot, :libmpfr), Cint,
          (FloatRefPnt, FloatRefPnt, FloatRefPnt, Cint), res, x, y, rm)
function log1p!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    x < -1 && throw(DomainError())
    ccall((:mpfr_log1p, :libmpfr), Cint,
          (FloatRefPnt, FloatRefPnt, Cint), res, x, rm)
end

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

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

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,
                  (FloatRefPnt, FloatRefPnt, Cint), res, x, rm)
        end
    end
end

function modf(x::FloatRef)
    isinf(x) && return (FloatRef(NaN), x)
    zint = FloatRef()
    zfloat = FloatRef()
    ccall((:mpfr_modf, :libmpfr), Cint,
          (FloatRefPnt, FloatRefPnt, FloatRefPnt, 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,
          (FloatRefPnt, FloatRefPnt, FloatRefPnt, 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,
          (FloatRefPnt, FloatRefPnt, FloatRefPnt, Cint), res, x, y, rm) # 

function sum!(res::FloatRef, arr::AbstractArray{FloatRef})
    for val in arr
        ccall((:mpfr_add, :libmpfr), Cint,
              (FloatRefPnt, FloatRefPnt, FloatRefPnt, 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,
                  (FloatRefPnt, FloatRefPnt, 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

# TODO: this needs to be changed for thread-safety!
const lgamma_signp = Array{Cint}(1)
lgamma!(res::FloatRef, x::FloatRef, rm::Cint = RM()) =
    ccall((:mpfr_lgamma,:libmpfr), Cint,
          (FloatRefPnt, Ptr{Cint}, FloatRefPnt, 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,
          (FloatRefPnt, FloatRefPnt, FloatRefPnt, Cint), res, y, x, rm)

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

for fC in (:(==), :(<=), :(>=), :(<), :(>))
    @eval ($fC){T<:Union{BigFlt,Flt}}(x::T, y::T) = _wrap_c2($fC, x, y)
end

function cmp(x::FloatRef, y::BigInt)
    isnan(x) && throw(DomainError())
    ccall((:mpfr_cmp_z, :libmpfr), Cint, (FloatRefPnt, Ref{BigInt}), x, y)
end
function cmp(x::FloatRef, y::ClongMax)
    isnan(x) && throw(DomainError())
    ccall((:mpfr_cmp_si, :libmpfr), Cint, (FloatRefPnt, Clong), x, y)
end
function cmp(x::FloatRef, y::CulongMax)
    isnan(x) && throw(DomainError())
    ccall((:mpfr_cmp_ui, :libmpfr), Cint, (FloatRefPnt, Culong), x, y)
end
function cmp(x::FloatRef, y::CdoubleMax)
    (isnan(x) || isnan(y)) && throw(DomainError())
    ccall((:mpfr_cmp_d, :libmpfr), Cint, (FloatRefPnt, Cdouble), x, y)
end
signbit(x::FloatRef) = ccall((:mpfr_signbit, :libmpfr), Cint, (FloatRefPnt,), x) != 0

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

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::FltTypes) = _wrap1(signbit, x)

for fC in (:add!, :mul!, :sub!, :div!, :idiv!, :rem!)
    @eval begin
        ($fC){T<:Union{BigFlt,Flt}}(res::FloatRef, x::T, y::T) =
            _wrap_res2($fC, res, x, y)
        ($fC){T<:Union{BigFlt,Flt}}(res::FloatRef, x::T, y::T, rm::Cint) =
            _wrap_res2($fC, res, x, y, rm)
    end
end
    
# precision of an object of type FloatRef
precision(x::FloatRef) = ccall((:mpfr_get_prec, :libmpfr), Clong, (FloatRefPnt,), 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)

precision!(::Type{FloatRef}, precision::Int) = _fill_fr(Clong(precision))

function precision!(x::FloatRef, precision::Int)
    x.prec = sign(x.prec)*precision
end

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,
          (FloatRefPnt, FloatRefPnt, FloatRefPnt, 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, (FloatRefPnt,), x) - 1
end

function frexp(x::FloatRef)
    z = FloatRef()
    c = Clong[0]
    ccall((:mpfr_frexp, :libmpfr), Cint,
          (Ptr{Clong}, FloatRefPnt, FloatRefPnt, 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}, FloatRefPnt, FloatRefPnt, Cint), c, res, x, rm)
    # Double the significand to make it work as Base.significand
    ccall((:mpfr_mul_si, :libmpfr), Cint,
          (FloatRefPnt, FloatRefPnt, Clong, Cint), res, res, 2, rm)
end

isinteger(x::FloatRef) = ccall((:mpfr_integer_p, :libmpfr), Cint, (FloatRefPnt,), 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,
                  (FloatRefPnt, FloatRefPnt), res, x)
        end
    end
end

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

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

isinf(x::FloatRef) = ccall((:mpfr_inf_p, :libmpfr), Cint, (FloatRefPnt,), x) != 0
isnan(x::FloatRef) = ccall((:mpfr_nan_p, :libmpfr), Cint, (FloatRefPnt,), 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, (FloatRefPnt, FloatRefPnt, Cint),
          res, x, rm)
    ccall((:mpfr_nextabove, :libmpfr), Cint, (FloatRefPnt,), res) != 0
end

function prevfloat!(res::FloatRef, x::FloatRef, rm::Cint = RM())
    ccall((:mpfr_set, :libmpfr), Cint, (FloatRefPnt, FloatRefPnt, Cint),
          res, x, rm)
    ccall((:mpfr_nextbelow, :libmpfr), Cint, (FloatRefPnt,), 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 begin
        typemax(::Type{$fT}) = ($fT)(Inf)
        typemin(::Type{$fT}) = ($fT)(-Inf)
        nextfloat(x::$fT, rm::Cint = RM()) = _wrap2(nextfloat!, x, rm)
        prevfloat(x::$fT, rm::Cint = RM()) = _wrap2(prevfloat!, x, rm)
        round(x::$fT, ::RoundingMode{:NearestTiesAway}) =
            _wrap2(round!, x, RoundingMode{:NearestTiesAway})
        eps(::Type{$fT})     = nextfloat($fT(1)) - $fT(1)
        realmin(::Type{$fT}) = nextfloat(zero($fT))
        realmax(::Type{$fT}) = prevfloat($fT(Inf))
        (/){T<:Union{$fT, NumTypes}}(x::$fT, y::T) = _wrap2(div!, x, y)
        div{T<:Union{$fT, NumTypes}}(x::$fT, y::T) = _wrap2(idiv!, x, y)
    end
    for (fJ,fC) in ((:+, :add!), (:*, :mul!), (:-, :sub!), (:^, :pow!))
        @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

    # Single argument wrapping
    for f in (:exp, :exp2, :exp10, :expm1, :digamma, :erf, :erfc, :zeta,
              :cosh,:sinh,:tanh,:sech,:csch,:coth, :cbrt, :factorial,
              :log, :log2, :log10, :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, :hypot)
        @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

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

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

^(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)

"""
    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)
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}, FloatRefPnt...), 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}, FloatRefPnt...), buf, "%.$(k)Re", x)
    elseif lng > k + 8
        buf = Array{UInt8}(lng + 1)
        lng = ccall((:mpfr_snprintf,:libmpfr), Cint,
                    (Ptr{UInt8}, Culong, Ptr{UInt8}, FloatRefPnt...), 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(x.prec, x.sign, x.exp, copy(x.vec))
    stackdict[x] = z
    z
end

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

end #module