Skip to content

Instantly share code, notes, and snippets.

@ScottPJones
Created June 29, 2016 03:19
Show Gist options
  • Save ScottPJones/37aef2198eb2bc37c5241887cb7117a2 to your computer and use it in GitHub Desktop.
Save ScottPJones/37aef2198eb2bc37c5241887cb7117a2 to your computer and use it in GitHub Desktop.
Latest (WIP) version of BigFlt (to replace BigFloat) module, with Flt type
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment