https://github.com/JuliaLang/julia
Tip revision: 26d7bb86c41a2343ad65e4da3e2bea2c7bc72e1c authored by Jameson Nash on 27 February 2016, 05:20:43 UTC
merge calls to identical methods (differing only by type signature) after doing union splitting
merge calls to identical methods (differing only by type signature) after doing union splitting
Tip revision: 26d7bb8
gmp.jl
# This file is a part of Julia. License is MIT: http://julialang.org/license
module GMP
export BigInt
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::Ptr{Limb}
function BigInt()
b = new(zero(Cint), zero(Cint), C_NULL)
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))
ccall((:__gmp_set_memory_functions, :libgmp), Void,
(Ptr{Void},Ptr{Void},Ptr{Void}),
cglobal(:jl_gc_counted_malloc),
cglobal(:jl_gc_counted_realloc_with_old_size),
cglobal(:jl_gc_counted_free))
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