stevenhao/gmp.jl

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

module GMP

export BigInt

type Hi
end

import Base: *, +, -, /, <, <<, >>, >>>, <=, ==, >, >=, ^, (~), (&), (|), ($),
             binomial, cmp, convert, div, divrem, factorial, fld, gcd, gcdx, lcm, mod,
             ndigits, promote_rule, rem, show, isqrt, string, isprime, powermod,
             sum, trailing_zeros, trailing_ones, count_ones, base, tryparse_internal,
             bin, oct, dec, hex, isequal, invmod, prevpow2, nextpow2, ndigits0z, widen, signed, unsafe_trunc, trunc

if Clong == Int32
    typealias ClongMax Union{Int8, Int16, Int32}
    typealias CulongMax Union{UInt8, UInt16, UInt32}
else
    typealias ClongMax Union{Int8, Int16, Int32, Int64}
    typealias CulongMax Union{UInt8, UInt16, UInt32, UInt64}
end
typealias CdoubleMax Union{Float16, Float32, Float64}

gmp_version() = VersionNumber(bytestring(unsafe_load(cglobal((:__gmp_version, :libgmp), Ptr{Cchar}))))
gmp_bits_per_limb() = Int(unsafe_load(cglobal((:__gmp_bits_per_limb, :libgmp), Cint)))

const GMP_VERSION = gmp_version()
const GMP_BITS_PER_LIMB = gmp_bits_per_limb()

# GMP's mp_limb_t is by default a typedef of `unsigned long`, but can also be configured to be either
# `unsigned int` or `unsigned long long int`. The correct unsigned type is here named Limb, and must
# be used whenever mp_limb_t is in the signature of ccall'ed GMP functions.
if GMP_BITS_PER_LIMB == 32
    typealias Limb UInt32
elseif GMP_BITS_PER_LIMB == 64
    typealias Limb UInt64
else
    error("GMP: cannot determine the type mp_limb_t (__gmp_bits_per_limb == $GMP_BITS_PER_LIMB)")
end


type BigInt <: Integer
    alloc::Cint
    size::Cint
    d::Tuple
    function BigInt()
        b = new(zero(Cint), zero(Cint), ())
        ccall((:__gmpz_init,:libgmp), Void, (Ptr{BigInt},), &b)
        # finalizer(b, _gmp_clear_func)
        return b
    end
end

_gmp_clear_func = C_NULL
_mpfr_clear_func = C_NULL

function __init__()
    try
        if gmp_version().major != GMP_VERSION.major || gmp_bits_per_limb() != GMP_BITS_PER_LIMB
            error(string("The dynamically loaded GMP library (version $(gmp_version()) with __gmp_bits_per_limb == $(gmp_bits_per_limb()))\n",
                         "does not correspond to the compile time version (version $GMP_VERSION with __gmp_bits_per_limb == $GMP_BITS_PER_LIMB).\n",
                         "Please rebuild Julia."))
        end

        global _gmp_clear_func = cglobal((:__gmpz_clear, :libgmp))
        global _mpfr_clear_func = cglobal((:mpfr_clear, :libmpfr))

        # JL_DLLEXPORT void *jl_gc_counted_malloc(size_t sz)
        function tuple_malloc(sz::Csize_t)
#            JL_DLLEXPORT jl_value_t *jl_new_struct_uninit(jl_datatype_t *type);
            ret = pointer_from_objref(ccall(:jl_new_struct_uninit, Any, (Any,), NTuple{Int64(sz), Int32}))
            # STATIC_INLINE void jl_gc_wb(void *parent, void *ptr)
            # how do we get the parent here?
            return ret
        end

        #JL_DLLEXPORT void *jl_gc_counted_realloc_with_old_size(void *p, size_t old,
        #                                               size_t sz)
        function tuple_realloc(p::Ptr{Void}, old::Csize_t, sz::Csize_t)
            ret = pointer_from_objref(ccall(:jl_new_struct_uninit, Any, (Any,), NTuple{Int64(sz), Int32}))
            return ret
        end

        #JL_DLLEXPORT void jl_gc_counted_free(void *p, size_t sz)
        function tuple_free(p::Ptr{Void}, sz::Csize_t)
            return # don't need to do anything; GC will handle freeing when it's time
        end

        ccall((:__gmp_set_memory_functions, :libgmp), Void,
              (Ptr{Void},Ptr{Void},Ptr{Void}),
              cfunction(tuple_malloc, Ptr{Void}, (Csize_t,)),
              cfunction(tuple_realloc, Ptr{Void}, (Ptr{Void}, Csize_t, Csize_t)),
              cfunction(tuple_free, Void, (Ptr{Void}, Csize_t)))
    catch ex
        Base.showerror_nostdio(ex,
            "WARNING: Error during initialization of module GMP")
    end
end

widen(::Type{Int128})  = BigInt
widen(::Type{UInt128}) = BigInt
widen(::Type{BigInt})  = BigInt

signed(x::BigInt) = x

convert(::Type{BigInt}, x::BigInt) = x

function tryparse_internal(::Type{BigInt}, s::AbstractString, startpos::Int, endpos::Int, base::Int, raise::Bool)
    _n = Nullable{BigInt}()

    # don't make a copy in the common case where we are parsing a whole bytestring
    bstr = startpos == start(s) && endpos == endof(s) ? bytestring(s) : bytestring(SubString(s,startpos,endpos))

    sgn, base, i = Base.parseint_preamble(true,base,bstr,start(bstr),endof(bstr))
    if i == 0
        raise && throw(ArgumentError("premature end of integer: $(repr(bstr))"))
        return _n
    end
    z = BigInt()
    if Base.containsnul(bstr)
        err = -1 # embedded NUL char (not handled correctly by GMP)
    else
        err = ccall((:__gmpz_set_str, :libgmp),
                    Int32, (Ptr{BigInt}, Ptr{UInt8}, Int32),
                    &z, pointer(bstr)+(i-start(bstr)), base)
    end
    if err != 0
        raise && throw(ArgumentError("invalid BigInt: $(repr(bstr))"))
        return _n
    end
    Nullable(sgn < 0 ? -z : z)
end

function convert(::Type{BigInt}, x::Union{Clong,Int32})
    z = BigInt()
    ccall((:__gmpz_set_si, :libgmp), Void, (Ptr{BigInt}, Clong), &z, x)
    return z
end
function convert(::Type{BigInt}, x::Union{Culong,UInt32})
    z = BigInt()
    ccall((:__gmpz_set_ui, :libgmp), Void, (Ptr{BigInt}, Culong), &z, x)
    return z
end

convert(::Type{BigInt}, x::Bool) = BigInt(UInt(x))


function unsafe_trunc(::Type{BigInt}, x::Union{Float32,Float64})
    z = BigInt()
    ccall((:__gmpz_set_d, :libgmp), Void, (Ptr{BigInt}, Cdouble), &z, x)
    return z
end

function convert(::Type{BigInt}, x::Union{Float32,Float64})
    isinteger(x) || throw(InexactError())
    unsafe_trunc(BigInt,x)
end

function trunc(::Type{BigInt}, x::Union{Float32,Float64})
    isfinite(x) || throw(InexactError())
    unsafe_trunc(BigInt,x)
end

convert(::Type{BigInt}, x::Float16) = BigInt(Float64(x))
convert(::Type{BigInt}, x::Float32) = BigInt(Float64(x))

function convert(::Type{BigInt}, x::Integer)
    if x < 0
        if typemin(Clong) <= x
            return BigInt(convert(Clong,x))
        end
        b = BigInt(0)
        shift = 0
        while x < -1
            b += BigInt(~UInt32(x&0xffffffff))<<shift
            x >>= 32
            shift += 32
        end
        return -b-1
    else
        if x <= typemax(Culong)
            return BigInt(convert(Culong,x))
        end
        b = BigInt(0)
        shift = 0
        while x > 0
            b += BigInt(UInt32(x&0xffffffff))<<shift
            x >>>= 32
            shift += 32
        end
        return b
    end
end


rem(x::BigInt, ::Type{Bool}) = ((x&1)!=0)
function rem{T<:Union{Unsigned,Signed}}(x::BigInt, ::Type{T})
    u = zero(T)
    for l = 1:min(abs(x.size), cld(sizeof(T),sizeof(Limb)))
        u += (unsafe_load(x.d,l)%T) << ((sizeof(Limb)<<3)*(l-1))
    end
    x.size < 0 ? -u : u
end

function convert{T<:Unsigned}(::Type{T}, x::BigInt)
    if sizeof(T) < sizeof(Limb)
        convert(T, convert(Limb,x))
    else
        0 <= x.size <= cld(sizeof(T),sizeof(Limb)) || throw(InexactError())
        x % T
    end
end

function convert{T<:Signed}(::Type{T}, x::BigInt)
    n = abs(x.size)
    if sizeof(T) < sizeof(Limb)
        SLimb = typeof(Signed(one(Limb)))
        convert(T, convert(SLimb, x))
    else
        0 <= n <= cld(sizeof(T),sizeof(Limb)) || throw(InexactError())
        y = x % T
        (x.size > 0) $ (y > 0) && throw(InexactError()) # catch overflow
        y
    end
end


function call(::Type{Float64}, n::BigInt, ::RoundingMode{:ToZero})
    ccall((:__gmpz_get_d, :libgmp), Float64, (Ptr{BigInt},), &n)
end

function call{T<:Union{Float16,Float32}}(::Type{T}, n::BigInt, ::RoundingMode{:ToZero})
    T(Float64(n,RoundToZero),RoundToZero)
end

function call{T<:CdoubleMax}(::Type{T}, n::BigInt, ::RoundingMode{:Down})
    x = T(n,RoundToZero)
    x > n ? prevfloat(x) : x
end
function call{T<:CdoubleMax}(::Type{T}, n::BigInt, ::RoundingMode{:Up})
    x = T(n,RoundToZero)
    x < n ? nextfloat(x) : x
end

function call{T<:CdoubleMax}(::Type{T}, n::BigInt, ::RoundingMode{:Nearest})
    x = T(n,RoundToZero)
    if maxintfloat(T) <= abs(x) < T(Inf)
        r = n-BigInt(x)
        h = eps(x)/2
        if iseven(reinterpret(Unsigned,x)) # check if last bit is odd/even
            if r < -h
                return prevfloat(x)
            elseif r > h
                return nextfloat(x)
            end
        else
            if r <= -h
                return prevfloat(x)
            elseif r >= h
                return nextfloat(x)
            end
        end
    end
    x
end

convert(::Type{Float64}, n::BigInt) = Float64(n,RoundNearest)
convert(::Type{Float32}, n::BigInt) = Float32(n,RoundNearest)
convert(::Type{Float16}, n::BigInt) = Float16(n,RoundNearest)

promote_rule{T<:Integer}(::Type{BigInt}, ::Type{T}) = BigInt

# Binary ops
for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
                 (:fld, :fdiv_q), (:div, :tdiv_q), (:mod, :fdiv_r), (:rem, :tdiv_r),
                 (:gcd, :gcd), (:lcm, :lcm),
                 (:&, :and), (:|, :ior), (:$, :xor))
    @eval begin
        function ($fJ)(x::BigInt, y::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y)
            return z
        end
    end
end

function invmod(x::BigInt, y::BigInt)
    z = BigInt()
    y = abs(y)
    if y == 1
        return big(0)
    end
    if (y==0 || ccall((:__gmpz_invert, :libgmp), Cint, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y) == 0)
        error("no inverse exists")
    end
    return z
end

# More efficient commutative operations
for (fJ, fC) in ((:+, :add), (:*, :mul), (:&, :and), (:|, :ior), (:$, :xor))
    @eval begin
        function ($fJ)(a::BigInt, b::BigInt, c::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
            return z
        end
        function ($fJ)(a::BigInt, b::BigInt, c::BigInt, d::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &d)
            return z
        end
        function ($fJ)(a::BigInt, b::BigInt, c::BigInt, d::BigInt, e::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &d)
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &e)
            return z
        end
    end
end

# Basic arithmetic without promotion
function +(x::BigInt, c::CulongMax)
    z = BigInt()
    ccall((:__gmpz_add_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end
+(c::CulongMax, x::BigInt) = x + c

function -(x::BigInt, c::CulongMax)
    z = BigInt()
    ccall((:__gmpz_sub_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end
function -(c::CulongMax, x::BigInt)
    z = BigInt()
    ccall((:__gmpz_ui_sub, :libgmp), Void, (Ptr{BigInt}, Culong, Ptr{BigInt}), &z, c, &x)
    return z
end
+(x::BigInt, c::ClongMax) = c < 0 ? -(x, -(c % Culong)) : x + convert(Culong, c)
+(c::ClongMax, x::BigInt) = c < 0 ? -(x, -(c % Culong)) : x + convert(Culong, c)
-(x::BigInt, c::ClongMax) = c < 0 ? +(x, -(c % Culong)) : -(x, convert(Culong, c))
-(c::ClongMax, x::BigInt) = c < 0 ? -(x + -(c % Culong)) : -(convert(Culong, c), x)

function *(x::BigInt, c::CulongMax)
    z = BigInt()
    ccall((:__gmpz_mul_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end
*(c::CulongMax, x::BigInt) = x * c
function *(x::BigInt, c::ClongMax)
    z = BigInt()
    ccall((:__gmpz_mul_si, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Clong), &z, &x, c)
    return z
end
*(c::ClongMax, x::BigInt) = x * c

# unary ops
for (fJ, fC) in ((:-, :neg), (:~, :com))
    @eval begin
        function ($fJ)(x::BigInt)
            z = BigInt()
            ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
            return z
        end
    end
end

function <<(x::BigInt, c::Int)
    c < 0 && throw(DomainError())
    c == 0 && return x
    z = BigInt()
    ccall((:__gmpz_mul_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end

function >>(x::BigInt, c::Int)
    c < 0 && throw(DomainError())
    c == 0 && return x
    z = BigInt()
    ccall((:__gmpz_fdiv_q_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
    return z
end

>>>(x::BigInt, c::Int) = x >> c

trailing_zeros(x::BigInt) = Int(ccall((:__gmpz_scan1, :libgmp), Culong, (Ptr{BigInt}, Culong), &x, 0))
trailing_ones(x::BigInt) = Int(ccall((:__gmpz_scan0, :libgmp), Culong, (Ptr{BigInt}, Culong), &x, 0))

count_ones(x::BigInt) = Int(ccall((:__gmpz_popcount, :libgmp), Culong, (Ptr{BigInt},), &x))

function divrem(x::BigInt, y::BigInt)
    z1 = BigInt()
    z2 = BigInt()
    ccall((:__gmpz_tdiv_qr, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z1, &z2, &x, &y)
    z1, z2
end

function cmp(x::BigInt, y::BigInt)
    ccall((:__gmpz_cmp, :libgmp), Int32, (Ptr{BigInt}, Ptr{BigInt}), &x, &y)
end
function cmp(x::BigInt, y::ClongMax)
    ccall((:__gmpz_cmp_si, :libgmp), Int32, (Ptr{BigInt}, Clong), &x, y)
end
function cmp(x::BigInt, y::CulongMax)
    ccall((:__gmpz_cmp_ui, :libgmp), Int32, (Ptr{BigInt}, Culong), &x, y)
end
cmp(x::BigInt, y::Integer) = cmp(x,big(y))
cmp(x::Integer, y::BigInt) = -cmp(y,x)

function cmp(x::BigInt, y::CdoubleMax)
    isnan(y) && throw(DomainError())
    ccall((:__gmpz_cmp_d, :libgmp), Int32, (Ptr{BigInt}, Cdouble), &x, y)
end
cmp(x::CdoubleMax, y::BigInt) = -cmp(y,x)

function isqrt(x::BigInt)
    z = BigInt()
    ccall((:__gmpz_sqrt, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
    return z
end

function ^(x::BigInt, y::Culong)
    z = BigInt()
    ccall((:__gmpz_pow_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, y)
    return z
end

function bigint_pow(x::BigInt, y::Integer)
    if y<0; throw(DomainError()); end
    if x== 1; return x; end
    if x==-1; return isodd(y) ? x : -x; end
    if y>typemax(Culong)
       x==0 && return x

       #At this point, x is not 1, 0 or -1 and it is not possible to use
       #gmpz_pow_ui to compute the answer. Note that the magnitude of the
       #answer is:
       #- at least 2^(2^32-1) ≈ 10^(1.3e9) (if Culong === UInt32).
       #- at least 2^(2^64-1) ≈ 10^(5.5e18) (if Culong === UInt64).
       #
       #Assume that the answer will definitely overflow.

       throw(OverflowError())
    end
    return x^convert(Culong, y)
end

^(x::BigInt , y::BigInt ) = bigint_pow(x, y)
^(x::BigInt , y::Bool   ) = y ? x : one(x)
^(x::BigInt , y::Integer) = bigint_pow(x, y)
^(x::Integer, y::BigInt ) = bigint_pow(BigInt(x), y)

function powermod(x::BigInt, p::BigInt, m::BigInt)
    p < 0 && throw(DomainError())
    r = BigInt()
    ccall((:__gmpz_powm, :libgmp), Void,
          (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
          &r, &x, &p, &m)
    return m < 0 && r > 0 ? r + m : r # choose sign conistent with mod(x^p, m)
end

powermod(x::Integer, p::Integer, m::BigInt) = powermod(big(x), big(p), m)

function gcdx(a::BigInt, b::BigInt)
    if b == 0 # shortcut this to ensure consistent results with gcdx(a,b)
        return a < 0 ? (-a,-one(BigInt),zero(BigInt)) : (a,one(BigInt),zero(BigInt))
    end
    g = BigInt()
    s = BigInt()
    t = BigInt()
    ccall((:__gmpz_gcdext, :libgmp), Void,
        (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
        &g, &s, &t, &a, &b)
    if t == 0
        # work around a difference in some versions of GMP
        if a == b
            return g, t, s
        elseif abs(a)==abs(b)
            return g, t, -s
        end
    end
    g, s, t
end

function sum(arr::AbstractArray{BigInt})
    n = BigInt(0)
    for i in arr
        ccall((:__gmpz_add, :libgmp), Void,
            (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
            &n, &n, &i)
    end
    return n
end

function factorial(x::BigInt)
    x.size < 0 && return BigInt(0)
    z = BigInt()
    ccall((:__gmpz_fac_ui, :libgmp), Void, (Ptr{BigInt}, Culong), &z, x)
    return z
end

function binomial(n::BigInt, k::UInt)
    z = BigInt()
    ccall((:__gmpz_bin_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &n, k)
    return z
end
binomial(n::BigInt, k::Integer) = k < 0 ? BigInt(0) : binomial(n, UInt(k))

==(x::BigInt, y::BigInt) = cmp(x,y) == 0
==(x::BigInt, i::Integer) = cmp(x,i) == 0
==(i::Integer, x::BigInt) = cmp(x,i) == 0
==(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) == 0
==(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) == 0

<=(x::BigInt, y::BigInt) = cmp(x,y) <= 0
<=(x::BigInt, i::Integer) = cmp(x,i) <= 0
<=(i::Integer, x::BigInt) = cmp(x,i) >= 0
<=(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) <= 0
<=(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) >= 0

<(x::BigInt, y::BigInt) = cmp(x,y) < 0
<(x::BigInt, i::Integer) = cmp(x,i) < 0
<(i::Integer, x::BigInt) = cmp(x,i) > 0
<(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) < 0
<(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) > 0

string(x::BigInt) = dec(x)
show(io::IO, x::BigInt) = print(io, string(x))

bin(n::BigInt) = base( 2, n)
oct(n::BigInt) = base( 8, n)
dec(n::BigInt) = base(10, n)
hex(n::BigInt) = base(16, n)

function base(b::Integer, n::BigInt)
    2 <= b <= 62 || throw(ArgumentError("base must be 2 ≤ base ≤ 62, got $b"))
    p = ccall((:__gmpz_get_str,:libgmp), Ptr{UInt8}, (Ptr{UInt8}, Cint, Ptr{BigInt}), C_NULL, b, &n)
    len = Int(ccall(:strlen, Csize_t, (Cstring,), p))
    ASCIIString(pointer_to_array(p,len,true))
end

function ndigits0z(x::BigInt, b::Integer=10)
    b < 2 && throw(DomainError())
    if ispow2(b)
        Int(ccall((:__gmpz_sizeinbase,:libgmp), Culong, (Ptr{BigInt}, Int32), &x, b))
    else
        # non-base 2 mpz_sizeinbase might return an answer 1 too big
        # use property that log(b, x) < ndigits(x, b) <= log(b, x) + 1
        n = Int(ccall((:__gmpz_sizeinbase,:libgmp), Culong, (Ptr{BigInt}, Int32), &x, 2))
        lb = log2(b) # assumed accurate to <1ulp (true for openlibm)
        q,r = divrem(n,lb)
        iq = Int(q)
        maxerr = q*eps(lb) # maximum error in remainder
        if r-1.0 < maxerr
            abs(x) >= big(b)^iq ? iq+1 : iq
        elseif lb-r < maxerr
            abs(x) >= big(b)^(iq+1) ? iq+2 : iq+1
        else
            iq+1
        end
    end
end
ndigits(x::BigInt, b::Integer=10) = x.size == 0 ? 1 : ndigits0z(x,b)

isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp), Cint, (Ptr{BigInt}, Cint), &x, reps) > 0

prevpow2(x::BigInt) = x.size < 0 ? -prevpow2(-x) : (x <= 2 ? x : one(BigInt) << (ndigits(x, 2)-1))
nextpow2(x::BigInt) = x.size < 0 ? -nextpow2(-x) : (x <= 2 ? x : one(BigInt) << ndigits(x-1, 2))

Base.checked_abs(x::BigInt) = abs(x)
Base.checked_neg(x::BigInt) = -x
Base.checked_add(a::BigInt, b::BigInt) = a + b
Base.checked_sub(a::BigInt, b::BigInt) = a - b
Base.checked_mul(a::BigInt, b::BigInt) = a * b
Base.checked_div(a::BigInt, b::BigInt) = div(a, b)
Base.checked_rem(a::BigInt, b::BigInt) = rem(a, b)
Base.checked_fld(a::BigInt, b::BigInt) = fld(a, b)
Base.checked_mod(a::BigInt, b::BigInt) = mod(a, b)
Base.checked_cld(a::BigInt, b::BigInt) = cld(a, b)

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

	module GMP

	export BigInt

	type Hi
	end

	import Base: *, +, -, /, <, <<, >>, >>>, <=, ==, >, >=, ^, (~), (&), (\|), ($),
	binomial, cmp, convert, div, divrem, factorial, fld, gcd, gcdx, lcm, mod,
	ndigits, promote_rule, rem, show, isqrt, string, isprime, powermod,
	sum, trailing_zeros, trailing_ones, count_ones, base, tryparse_internal,
	bin, oct, dec, hex, isequal, invmod, prevpow2, nextpow2, ndigits0z, widen, signed, unsafe_trunc, trunc

	if Clong == Int32
	typealias ClongMax Union{Int8, Int16, Int32}
	typealias CulongMax Union{UInt8, UInt16, UInt32}
	else
	typealias ClongMax Union{Int8, Int16, Int32, Int64}
	typealias CulongMax Union{UInt8, UInt16, UInt32, UInt64}
	end
	typealias CdoubleMax Union{Float16, Float32, Float64}

	gmp_version() = VersionNumber(bytestring(unsafe_load(cglobal((:__gmp_version, :libgmp), Ptr{Cchar}))))
	gmp_bits_per_limb() = Int(unsafe_load(cglobal((:__gmp_bits_per_limb, :libgmp), Cint)))

	const GMP_VERSION = gmp_version()
	const GMP_BITS_PER_LIMB = gmp_bits_per_limb()

	# GMP's mp_limb_t is by default a typedef of `unsigned long`, but can also be configured to be either
	# `unsigned int` or `unsigned long long int`. The correct unsigned type is here named Limb, and must
	# be used whenever mp_limb_t is in the signature of ccall'ed GMP functions.
	if GMP_BITS_PER_LIMB == 32
	typealias Limb UInt32
	elseif GMP_BITS_PER_LIMB == 64
	typealias Limb UInt64
	else
	error("GMP: cannot determine the type mp_limb_t (__gmp_bits_per_limb == $GMP_BITS_PER_LIMB)")
	end


	type BigInt <: Integer
	alloc::Cint
	size::Cint
	d::Tuple
	function BigInt()
	b = new(zero(Cint), zero(Cint), ())
	ccall((:__gmpz_init,:libgmp), Void, (Ptr{BigInt},), &b)
	# finalizer(b, _gmp_clear_func)
	return b
	end
	end

	_gmp_clear_func = C_NULL
	_mpfr_clear_func = C_NULL

	function __init__()
	try
	if gmp_version().major != GMP_VERSION.major \|\| gmp_bits_per_limb() != GMP_BITS_PER_LIMB
	error(string("The dynamically loaded GMP library (version $(gmp_version()) with __gmp_bits_per_limb == $(gmp_bits_per_limb()))\n",
	"does not correspond to the compile time version (version $GMP_VERSION with __gmp_bits_per_limb == $GMP_BITS_PER_LIMB).\n",
	"Please rebuild Julia."))
	end

	global _gmp_clear_func = cglobal((:__gmpz_clear, :libgmp))
	global _mpfr_clear_func = cglobal((:mpfr_clear, :libmpfr))

	# JL_DLLEXPORT void *jl_gc_counted_malloc(size_t sz)
	function tuple_malloc(sz::Csize_t)
	# JL_DLLEXPORT jl_value_t jl_new_struct_uninit(jl_datatype_t type);
	ret = pointer_from_objref(ccall(:jl_new_struct_uninit, Any, (Any,), NTuple{Int64(sz), Int32}))
	# STATIC_INLINE void jl_gc_wb(void parent, void ptr)
	# how do we get the parent here?
	return ret
	end

	#JL_DLLEXPORT void jl_gc_counted_realloc_with_old_size(void p, size_t old,
	# size_t sz)
	function tuple_realloc(p::Ptr{Void}, old::Csize_t, sz::Csize_t)
	ret = pointer_from_objref(ccall(:jl_new_struct_uninit, Any, (Any,), NTuple{Int64(sz), Int32}))
	return ret
	end

	#JL_DLLEXPORT void jl_gc_counted_free(void *p, size_t sz)
	function tuple_free(p::Ptr{Void}, sz::Csize_t)
	return # don't need to do anything; GC will handle freeing when it's time
	end

	ccall((:__gmp_set_memory_functions, :libgmp), Void,
	(Ptr{Void},Ptr{Void},Ptr{Void}),
	cfunction(tuple_malloc, Ptr{Void}, (Csize_t,)),
	cfunction(tuple_realloc, Ptr{Void}, (Ptr{Void}, Csize_t, Csize_t)),
	cfunction(tuple_free, Void, (Ptr{Void}, Csize_t)))
	catch ex
	Base.showerror_nostdio(ex,
	"WARNING: Error during initialization of module GMP")
	end
	end

	widen(::Type{Int128}) = BigInt
	widen(::Type{UInt128}) = BigInt
	widen(::Type{BigInt}) = BigInt

	signed(x::BigInt) = x

	convert(::Type{BigInt}, x::BigInt) = x

	function tryparse_internal(::Type{BigInt}, s::AbstractString, startpos::Int, endpos::Int, base::Int, raise::Bool)
	_n = Nullable{BigInt}()

	# don't make a copy in the common case where we are parsing a whole bytestring
	bstr = startpos == start(s) && endpos == endof(s) ? bytestring(s) : bytestring(SubString(s,startpos,endpos))

	sgn, base, i = Base.parseint_preamble(true,base,bstr,start(bstr),endof(bstr))
	if i == 0
	raise && throw(ArgumentError("premature end of integer: $(repr(bstr))"))
	return _n
	end
	z = BigInt()
	if Base.containsnul(bstr)
	err = -1 # embedded NUL char (not handled correctly by GMP)
	else
	err = ccall((:__gmpz_set_str, :libgmp),
	Int32, (Ptr{BigInt}, Ptr{UInt8}, Int32),
	&z, pointer(bstr)+(i-start(bstr)), base)
	end
	if err != 0
	raise && throw(ArgumentError("invalid BigInt: $(repr(bstr))"))
	return _n
	end
	Nullable(sgn < 0 ? -z : z)
	end

	function convert(::Type{BigInt}, x::Union{Clong,Int32})
	z = BigInt()
	ccall((:__gmpz_set_si, :libgmp), Void, (Ptr{BigInt}, Clong), &z, x)
	return z
	end
	function convert(::Type{BigInt}, x::Union{Culong,UInt32})
	z = BigInt()
	ccall((:__gmpz_set_ui, :libgmp), Void, (Ptr{BigInt}, Culong), &z, x)
	return z
	end

	convert(::Type{BigInt}, x::Bool) = BigInt(UInt(x))


	function unsafe_trunc(::Type{BigInt}, x::Union{Float32,Float64})
	z = BigInt()
	ccall((:__gmpz_set_d, :libgmp), Void, (Ptr{BigInt}, Cdouble), &z, x)
	return z
	end

	function convert(::Type{BigInt}, x::Union{Float32,Float64})
	isinteger(x) \|\| throw(InexactError())
	unsafe_trunc(BigInt,x)
	end

	function trunc(::Type{BigInt}, x::Union{Float32,Float64})
	isfinite(x) \|\| throw(InexactError())
	unsafe_trunc(BigInt,x)
	end

	convert(::Type{BigInt}, x::Float16) = BigInt(Float64(x))
	convert(::Type{BigInt}, x::Float32) = BigInt(Float64(x))

	function convert(::Type{BigInt}, x::Integer)
	if x < 0
	if typemin(Clong) <= x
	return BigInt(convert(Clong,x))
	end
	b = BigInt(0)
	shift = 0
	while x < -1
	b += BigInt(~UInt32(x&0xffffffff))<<shift
	x >>= 32
	shift += 32
	end
	return -b-1
	else
	if x <= typemax(Culong)
	return BigInt(convert(Culong,x))
	end
	b = BigInt(0)
	shift = 0
	while x > 0
	b += BigInt(UInt32(x&0xffffffff))<<shift
	x >>>= 32
	shift += 32
	end
	return b
	end
	end


	rem(x::BigInt, ::Type{Bool}) = ((x&1)!=0)
	function rem{T<:Union{Unsigned,Signed}}(x::BigInt, ::Type{T})
	u = zero(T)
	for l = 1:min(abs(x.size), cld(sizeof(T),sizeof(Limb)))
	u += (unsafe_load(x.d,l)%T) << ((sizeof(Limb)<<3)*(l-1))
	end
	x.size < 0 ? -u : u
	end

	function convert{T<:Unsigned}(::Type{T}, x::BigInt)
	if sizeof(T) < sizeof(Limb)
	convert(T, convert(Limb,x))
	else
	0 <= x.size <= cld(sizeof(T),sizeof(Limb)) \|\| throw(InexactError())
	x % T
	end
	end

	function convert{T<:Signed}(::Type{T}, x::BigInt)
	n = abs(x.size)
	if sizeof(T) < sizeof(Limb)
	SLimb = typeof(Signed(one(Limb)))
	convert(T, convert(SLimb, x))
	else
	0 <= n <= cld(sizeof(T),sizeof(Limb)) \|\| throw(InexactError())
	y = x % T
	(x.size > 0) $ (y > 0) && throw(InexactError()) # catch overflow
	y
	end
	end


	function call(::Type{Float64}, n::BigInt, ::RoundingMode{:ToZero})
	ccall((:__gmpz_get_d, :libgmp), Float64, (Ptr{BigInt},), &n)
	end

	function call{T<:Union{Float16,Float32}}(::Type{T}, n::BigInt, ::RoundingMode{:ToZero})
	T(Float64(n,RoundToZero),RoundToZero)
	end

	function call{T<:CdoubleMax}(::Type{T}, n::BigInt, ::RoundingMode{:Down})
	x = T(n,RoundToZero)
	x > n ? prevfloat(x) : x
	end
	function call{T<:CdoubleMax}(::Type{T}, n::BigInt, ::RoundingMode{:Up})
	x = T(n,RoundToZero)
	x < n ? nextfloat(x) : x
	end

	function call{T<:CdoubleMax}(::Type{T}, n::BigInt, ::RoundingMode{:Nearest})
	x = T(n,RoundToZero)
	if maxintfloat(T) <= abs(x) < T(Inf)
	r = n-BigInt(x)
	h = eps(x)/2
	if iseven(reinterpret(Unsigned,x)) # check if last bit is odd/even
	if r < -h
	return prevfloat(x)
	elseif r > h
	return nextfloat(x)
	end
	else
	if r <= -h
	return prevfloat(x)
	elseif r >= h
	return nextfloat(x)
	end
	end
	end
	x
	end

	convert(::Type{Float64}, n::BigInt) = Float64(n,RoundNearest)
	convert(::Type{Float32}, n::BigInt) = Float32(n,RoundNearest)
	convert(::Type{Float16}, n::BigInt) = Float16(n,RoundNearest)

	promote_rule{T<:Integer}(::Type{BigInt}, ::Type{T}) = BigInt

	# Binary ops
	for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
	(:fld, :fdiv_q), (:div, :tdiv_q), (:mod, :fdiv_r), (:rem, :tdiv_r),
	(:gcd, :gcd), (:lcm, :lcm),
	(:&, :and), (:\|, :ior), (:$, :xor))
	@eval begin
	function ($fJ)(x::BigInt, y::BigInt)
	z = BigInt()
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y)
	return z
	end
	end
	end

	function invmod(x::BigInt, y::BigInt)
	z = BigInt()
	y = abs(y)
	if y == 1
	return big(0)
	end
	if (y==0 \|\| ccall((:__gmpz_invert, :libgmp), Cint, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y) == 0)
	error("no inverse exists")
	end
	return z
	end

	# More efficient commutative operations
	for (fJ, fC) in ((:+, :add), (:*, :mul), (:&, :and), (:\|, :ior), (:$, :xor))
	@eval begin
	function ($fJ)(a::BigInt, b::BigInt, c::BigInt)
	z = BigInt()
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
	return z
	end
	function ($fJ)(a::BigInt, b::BigInt, c::BigInt, d::BigInt)
	z = BigInt()
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &d)
	return z
	end
	function ($fJ)(a::BigInt, b::BigInt, c::BigInt, d::BigInt, e::BigInt)
	z = BigInt()
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &a, &b)
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &c)
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &d)
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &z, &e)
	return z
	end
	end
	end

	# Basic arithmetic without promotion
	function +(x::BigInt, c::CulongMax)
	z = BigInt()
	ccall((:__gmpz_add_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
	return z
	end
	+(c::CulongMax, x::BigInt) = x + c

	function -(x::BigInt, c::CulongMax)
	z = BigInt()
	ccall((:__gmpz_sub_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
	return z
	end
	function -(c::CulongMax, x::BigInt)
	z = BigInt()
	ccall((:__gmpz_ui_sub, :libgmp), Void, (Ptr{BigInt}, Culong, Ptr{BigInt}), &z, c, &x)
	return z
	end
	+(x::BigInt, c::ClongMax) = c < 0 ? -(x, -(c % Culong)) : x + convert(Culong, c)
	+(c::ClongMax, x::BigInt) = c < 0 ? -(x, -(c % Culong)) : x + convert(Culong, c)
	-(x::BigInt, c::ClongMax) = c < 0 ? +(x, -(c % Culong)) : -(x, convert(Culong, c))
	-(c::ClongMax, x::BigInt) = c < 0 ? -(x + -(c % Culong)) : -(convert(Culong, c), x)

	function *(x::BigInt, c::CulongMax)
	z = BigInt()
	ccall((:__gmpz_mul_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
	return z
	end
	(c::CulongMax, x::BigInt) = x c
	function *(x::BigInt, c::ClongMax)
	z = BigInt()
	ccall((:__gmpz_mul_si, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Clong), &z, &x, c)
	return z
	end
	(c::ClongMax, x::BigInt) = x c

	# unary ops
	for (fJ, fC) in ((:-, :neg), (:~, :com))
	@eval begin
	function ($fJ)(x::BigInt)
	z = BigInt()
	ccall(($(string(:__gmpz_,fC)), :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
	return z
	end
	end
	end

	function <<(x::BigInt, c::Int)
	c < 0 && throw(DomainError())
	c == 0 && return x
	z = BigInt()
	ccall((:__gmpz_mul_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
	return z
	end

	function >>(x::BigInt, c::Int)
	c < 0 && throw(DomainError())
	c == 0 && return x
	z = BigInt()
	ccall((:__gmpz_fdiv_q_2exp, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, c)
	return z
	end

	>>>(x::BigInt, c::Int) = x >> c

	trailing_zeros(x::BigInt) = Int(ccall((:__gmpz_scan1, :libgmp), Culong, (Ptr{BigInt}, Culong), &x, 0))
	trailing_ones(x::BigInt) = Int(ccall((:__gmpz_scan0, :libgmp), Culong, (Ptr{BigInt}, Culong), &x, 0))

	count_ones(x::BigInt) = Int(ccall((:__gmpz_popcount, :libgmp), Culong, (Ptr{BigInt},), &x))

	function divrem(x::BigInt, y::BigInt)
	z1 = BigInt()
	z2 = BigInt()
	ccall((:__gmpz_tdiv_qr, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z1, &z2, &x, &y)
	z1, z2
	end

	function cmp(x::BigInt, y::BigInt)
	ccall((:__gmpz_cmp, :libgmp), Int32, (Ptr{BigInt}, Ptr{BigInt}), &x, &y)
	end
	function cmp(x::BigInt, y::ClongMax)
	ccall((:__gmpz_cmp_si, :libgmp), Int32, (Ptr{BigInt}, Clong), &x, y)
	end
	function cmp(x::BigInt, y::CulongMax)
	ccall((:__gmpz_cmp_ui, :libgmp), Int32, (Ptr{BigInt}, Culong), &x, y)
	end
	cmp(x::BigInt, y::Integer) = cmp(x,big(y))
	cmp(x::Integer, y::BigInt) = -cmp(y,x)

	function cmp(x::BigInt, y::CdoubleMax)
	isnan(y) && throw(DomainError())
	ccall((:__gmpz_cmp_d, :libgmp), Int32, (Ptr{BigInt}, Cdouble), &x, y)
	end
	cmp(x::CdoubleMax, y::BigInt) = -cmp(y,x)

	function isqrt(x::BigInt)
	z = BigInt()
	ccall((:__gmpz_sqrt, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}), &z, &x)
	return z
	end

	function ^(x::BigInt, y::Culong)
	z = BigInt()
	ccall((:__gmpz_pow_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &x, y)
	return z
	end

	function bigint_pow(x::BigInt, y::Integer)
	if y<0; throw(DomainError()); end
	if x== 1; return x; end
	if x==-1; return isodd(y) ? x : -x; end
	if y>typemax(Culong)
	x==0 && return x

	#At this point, x is not 1, 0 or -1 and it is not possible to use
	#gmpz_pow_ui to compute the answer. Note that the magnitude of the
	#answer is:
	#- at least 2^(2^32-1) ≈ 10^(1.3e9) (if Culong === UInt32).
	#- at least 2^(2^64-1) ≈ 10^(5.5e18) (if Culong === UInt64).
	#
	#Assume that the answer will definitely overflow.

	throw(OverflowError())
	end
	return x^convert(Culong, y)
	end

	^(x::BigInt , y::BigInt ) = bigint_pow(x, y)
	^(x::BigInt , y::Bool ) = y ? x : one(x)
	^(x::BigInt , y::Integer) = bigint_pow(x, y)
	^(x::Integer, y::BigInt ) = bigint_pow(BigInt(x), y)

	function powermod(x::BigInt, p::BigInt, m::BigInt)
	p < 0 && throw(DomainError())
	r = BigInt()
	ccall((:__gmpz_powm, :libgmp), Void,
	(Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
	&r, &x, &p, &m)
	return m < 0 && r > 0 ? r + m : r # choose sign conistent with mod(x^p, m)
	end

	powermod(x::Integer, p::Integer, m::BigInt) = powermod(big(x), big(p), m)

	function gcdx(a::BigInt, b::BigInt)
	if b == 0 # shortcut this to ensure consistent results with gcdx(a,b)
	return a < 0 ? (-a,-one(BigInt),zero(BigInt)) : (a,one(BigInt),zero(BigInt))
	end
	g = BigInt()
	s = BigInt()
	t = BigInt()
	ccall((:__gmpz_gcdext, :libgmp), Void,
	(Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
	&g, &s, &t, &a, &b)
	if t == 0
	# work around a difference in some versions of GMP
	if a == b
	return g, t, s
	elseif abs(a)==abs(b)
	return g, t, -s
	end
	end
	g, s, t
	end

	function sum(arr::AbstractArray{BigInt})
	n = BigInt(0)
	for i in arr
	ccall((:__gmpz_add, :libgmp), Void,
	(Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
	&n, &n, &i)
	end
	return n
	end

	function factorial(x::BigInt)
	x.size < 0 && return BigInt(0)
	z = BigInt()
	ccall((:__gmpz_fac_ui, :libgmp), Void, (Ptr{BigInt}, Culong), &z, x)
	return z
	end

	function binomial(n::BigInt, k::UInt)
	z = BigInt()
	ccall((:__gmpz_bin_ui, :libgmp), Void, (Ptr{BigInt}, Ptr{BigInt}, Culong), &z, &n, k)
	return z
	end
	binomial(n::BigInt, k::Integer) = k < 0 ? BigInt(0) : binomial(n, UInt(k))

	==(x::BigInt, y::BigInt) = cmp(x,y) == 0
	==(x::BigInt, i::Integer) = cmp(x,i) == 0
	==(i::Integer, x::BigInt) = cmp(x,i) == 0
	==(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) == 0
	==(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) == 0

	<=(x::BigInt, y::BigInt) = cmp(x,y) <= 0
	<=(x::BigInt, i::Integer) = cmp(x,i) <= 0
	<=(i::Integer, x::BigInt) = cmp(x,i) >= 0
	<=(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) <= 0
	<=(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) >= 0

	<(x::BigInt, y::BigInt) = cmp(x,y) < 0
	<(x::BigInt, i::Integer) = cmp(x,i) < 0
	<(i::Integer, x::BigInt) = cmp(x,i) > 0
	<(x::BigInt, f::CdoubleMax) = isnan(f) ? false : cmp(x,f) < 0
	<(f::CdoubleMax, x::BigInt) = isnan(f) ? false : cmp(x,f) > 0

	string(x::BigInt) = dec(x)
	show(io::IO, x::BigInt) = print(io, string(x))

	bin(n::BigInt) = base( 2, n)
	oct(n::BigInt) = base( 8, n)
	dec(n::BigInt) = base(10, n)
	hex(n::BigInt) = base(16, n)

	function base(b::Integer, n::BigInt)
	2 <= b <= 62 \|\| throw(ArgumentError("base must be 2 ≤ base ≤ 62, got $b"))
	p = ccall((:__gmpz_get_str,:libgmp), Ptr{UInt8}, (Ptr{UInt8}, Cint, Ptr{BigInt}), C_NULL, b, &n)
	len = Int(ccall(:strlen, Csize_t, (Cstring,), p))
	ASCIIString(pointer_to_array(p,len,true))
	end

	function ndigits0z(x::BigInt, b::Integer=10)
	b < 2 && throw(DomainError())
	if ispow2(b)
	Int(ccall((:__gmpz_sizeinbase,:libgmp), Culong, (Ptr{BigInt}, Int32), &x, b))
	else
	# non-base 2 mpz_sizeinbase might return an answer 1 too big
	# use property that log(b, x) < ndigits(x, b) <= log(b, x) + 1
	n = Int(ccall((:__gmpz_sizeinbase,:libgmp), Culong, (Ptr{BigInt}, Int32), &x, 2))
	lb = log2(b) # assumed accurate to <1ulp (true for openlibm)
	q,r = divrem(n,lb)
	iq = Int(q)
	maxerr = q*eps(lb) # maximum error in remainder
	if r-1.0 < maxerr
	abs(x) >= big(b)^iq ? iq+1 : iq
	elseif lb-r < maxerr
	abs(x) >= big(b)^(iq+1) ? iq+2 : iq+1
	else
	iq+1
	end
	end
	end
	ndigits(x::BigInt, b::Integer=10) = x.size == 0 ? 1 : ndigits0z(x,b)

	isprime(x::BigInt, reps=25) = ccall((:__gmpz_probab_prime_p,:libgmp), Cint, (Ptr{BigInt}, Cint), &x, reps) > 0

	prevpow2(x::BigInt) = x.size < 0 ? -prevpow2(-x) : (x <= 2 ? x : one(BigInt) << (ndigits(x, 2)-1))
	nextpow2(x::BigInt) = x.size < 0 ? -nextpow2(-x) : (x <= 2 ? x : one(BigInt) << ndigits(x-1, 2))

	Base.checked_abs(x::BigInt) = abs(x)
	Base.checked_neg(x::BigInt) = -x
	Base.checked_add(a::BigInt, b::BigInt) = a + b
	Base.checked_sub(a::BigInt, b::BigInt) = a - b
	Base.checked_mul(a::BigInt, b::BigInt) = a * b
	Base.checked_div(a::BigInt, b::BigInt) = div(a, b)
	Base.checked_rem(a::BigInt, b::BigInt) = rem(a, b)
	Base.checked_fld(a::BigInt, b::BigInt) = fld(a, b)
	Base.checked_mod(a::BigInt, b::BigInt) = mod(a, b)
	Base.checked_cld(a::BigInt, b::BigInt) = cld(a, b)

	end # module