fmpz.jl

###############################################################################
#
#   ZZRingElem.jl : BigInts
#
###############################################################################

# Copyright (c) 2009-2014: Jeff Bezanson, Stefan Karpinski, Viral B. Shah,
# and other contributors:
#
# https://github.com/JuliaLang/julia/contributors
#
# Copyright (C) 2014, 2015 William Hart
# Copyright (C) 2015, Claus Fieker
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.


###############################################################################
#
#   Data type and parent methods
#
###############################################################################

parent_type(::Type{ZZRingElem}) = ZZRing

@doc raw"""
    parent(a::ZZRingElem)

Returns the unique Flint integer parent object `ZZ`.
"""
parent(a::ZZRingElem) = ZZ

elem_type(::Type{ZZRing}) = ZZRingElem

base_ring_type(::Type{ZZRing}) = typeof(Union{})

base_ring(a::ZZRing) = Union{}

is_domain_type(::Type{ZZRingElem}) = true

###############################################################################
#
#   Internal data
#
###############################################################################

data(a::ZZRingElem) = a.d
data(a::Ref{ZZRingElem}) = a[].d
data(a::Ptr{ZZRingElem}) = unsafe_load(reinterpret(Ptr{Int}, a))

###############################################################################
#
#   Ranges
#
###############################################################################

# ZZRingElem needs to be iterable for
# Base.throw_boundserror(something, ::ZZRingElem) to print correctly
Base.iterate(x::ZZRingElem) = (x, nothing)
Base.iterate(x::ZZRingElem, i::Any) = nothing

# Note that we cannot get a UnitRange as this is only legal for subtypes of Real.
# So, we use an AbstractUnitRange here mostly copied from `base/range.jl`.
# `StepRange`s on the other hand work out of the box thanks to duck typing.

# See FlintTypes.jl for the struct ZZRingElemUnitRange

fmpz_unitrange_last(start::ZZRingElem, stop::ZZRingElem) =
ifelse(stop >= start, stop, start - one(ZZRingElem))

Base.:(:)(a::ZZRingElem, b::ZZRingElem) = ZZRingElemUnitRange(a, b)

@inline function getindex(r::ZZRingElemUnitRange, i::ZZRingElem)
  val = r.start + (i - 1)
  @boundscheck _in_unit_range(r, val) || Base.throw_boundserror(r, i)
  val
end
_in_unit_range(r::ZZRingElemUnitRange, val::ZZRingElem) = r.start <= val <= r.stop

show(io::IO, r::ZZRingElemUnitRange) = print(io, repr(first(r)), ':', repr(last(r)))

in(x::IntegerUnion, r::ZZRingElemUnitRange) = first(r) <= x <= last(r)

mod(i::IntegerUnion, r::ZZRingElemUnitRange) = mod(i - first(r), length(r)) + first(r)

Base.:(:)(a::ZZRingElem, b::Integer) = (:)(promote(a, b)...)
Base.:(:)(a::Integer, b::ZZRingElem) = (:)(promote(a, b)...)

Base.:(:)(x::Int, y::ZZRingElem) = ZZRingElem(x):y
Base.:(:)(x::ZZRingElem, y::Int) = x:ZZRingElem(y)

# Construct StepRange{ZZRingElem, T} where +(::ZZRingElem, zero(::T)) must be defined
Base.:(:)(a::ZZRingElem, s, b::Integer) = ((a_, b_) = promote(a, b); a_:s:b_)
Base.:(:)(a::Integer, s, b::ZZRingElem) = ((a_, b_) = promote(a, b); a_:s:b_)

# `length` should return an Integer, so BigInt seems appropriate as ZZRingElem is not <: Integer
# this method is useful in particular to enable rand(ZZ(n):ZZ(m))
function Base.length(r::StepRange{ZZRingElem})
  n = div((last(r) - first(r)) + step(r), step(r))
  isempty(r) ? zero(BigInt) : BigInt(n)
end

function Base.in(x::IntegerUnion, r::AbstractRange{ZZRingElem})
  if isempty(r) || first(r) > x || x > last(r)
    return false
  end
  return mod(convert(ZZRingElem, x), step(r)) == mod(first(r), step(r))
end

function Base.getindex(a::StepRange{ZZRingElem}, i::ZZRingElem)
  res = first(a) + (i - 1) * Base.step(a)
  ok = false
  if step(a) > 0
    ok = res <= last(a) && res >= first(a)
  else
    ok = res >= last(a) && res <= first(a)
  end
  @boundscheck ((i > 0) && ok) || Base.throw_boundserror(a, i)
  return res
end

################################################################################
#
#   Hashing
#
################################################################################

# Similar to hash for BigInt found in julia/base

@inline function _fmpz_is_small(a::ZZRingElem)
  return __fmpz_is_small(a.d)
end

@inline function _fmpz_limbs(a::ZZRingElem)
  return __fmpz_limbs(a.d)
end

function hash_integer(a::ZZRingElem, h::UInt)
  return _hash_integer(a.d, h)
end

function hash(a::ZZRingElem, h::UInt)
  return hash_integer(a, h)
end

@inline function __fmpz_is_small(a::Int)
  return (unsigned(a) >> (Sys.WORD_SIZE - 2) != 1)
end

function __fmpz_limbs(a::Int)
  if __fmpz_is_small(a)
    return Cint(0)
  end
  b = unsafe_load(convert(Ptr{Cint}, unsigned(a)<<2), 2)
  return b
end

function _hash_integer(a::Int, h::UInt)
  s::Cint = __fmpz_limbs(a)
  s == 0 && return Base.hash_integer(a, h)
  # get the pointer after the first two Cint
  d = convert(Ptr{Ptr{UInt}}, unsigned(a) << 2) + 2*sizeof(Cint)
  p = unsafe_load(d)
  b = unsafe_load(p)
  h = xor(Base.hash_uint(xor(ifelse(s < 0, -b, b), h)), h)
  for k = 2:abs(s)
    h = xor(Base.hash_uint(xor(unsafe_load(p, k), h)), h)
  end
  return h
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function deepcopy_internal(a::ZZRingElem, dict::IdDict)
  z = ZZRingElem()
  set!(z, a)
  return z
end

Base.copy(a::ZZRingElem) = deepcopy(a)

characteristic(R::ZZRing) = 0

one(::ZZRing) = ZZRingElem(1)

zero(::ZZRing) = ZZRingElem(0)

one(::Type{ZZRingElem}) = ZZRingElem(1)

zero(::Type{ZZRingElem}) = ZZRingElem(0)

one(::ZZRingElem) = ZZRingElem(1)

zero(::ZZRingElem) = ZZRingElem(0)

@doc raw"""
    sign(a::ZZRingElem)

Return the sign of $a$, i.e. $+1$, $0$ or $-1$.
"""
sign(a::ZZRingElem) = ZZRingElem(ccall((:fmpz_sgn, libflint), Cint, (Ref{ZZRingElem},), a))

sign(::Type{Int}, a::ZZRingElem) = Int(ccall((:fmpz_sgn, libflint), Cint, (Ref{ZZRingElem},), a))

Base.signbit(a::ZZRingElem) = signbit(sign(Int, a))

is_negative(n::ZZRingElem) = sign(Int, n) < 0
is_positive(n::ZZRingElem) = sign(Int, n) > 0

@doc raw"""
    fits(::Type{Int}, a::ZZRingElem)

Return `true` if $a$ fits into an `Int`, otherwise return `false`.
"""
fits(::Type{Int}, a::ZZRingElem) = ccall((:fmpz_fits_si, libflint), Bool,
                                         (Ref{ZZRingElem},), a)

@doc raw"""
    fits(::Type{UInt}, a::ZZRingElem)

Return `true` if $a$ fits into a `UInt`, otherwise return `false`.
"""
@inline function fits(::Type{UInt}, a::ZZRingElem)
  a < 0 && return false
  return ccall((:fmpz_abs_fits_ui, libflint), Bool, (Ref{ZZRingElem},), a)
end

if Culong !== UInt
  function fits(::Type{Culong}, a::ZZRingElem)
    return 0 <= a && a <= UInt(typemax(Culong))
  end
end

@doc raw"""
    size(a::ZZRingElem)

Return the number of limbs required to store the absolute value of $a$.
"""
size(a::ZZRingElem) = Int(ccall((:fmpz_size, libflint), Cint, (Ref{ZZRingElem},), a))

is_unit(a::ZZRingElemOrPtr) = data(a) == 1 || data(a) == -1

is_zero(a::ZZRingElemOrPtr) = data(a) == 0

is_one(a::ZZRingElemOrPtr) = data(a) == 1

isinteger(::ZZRingElem) = true

isfinite(::ZZRingElem) = true

isinf(::ZZRingElem) = false

@doc raw"""
    denominator(a::ZZRingElem)

Return the denominator of $a$ thought of as a rational. Always returns $1$.
"""
function denominator(a::ZZRingElem)
  return ZZRingElem(1)
end

@doc raw"""
    numerator(a::ZZRingElem)

Return the numerator of $a$ thought of as a rational. Always returns $a$.
"""
function numerator(a::ZZRingElem)
  return a
end

isodd(a::ZZRingElem)  = ccall((:fmpz_is_odd,  libflint), Cint, (Ref{ZZRingElem},), a) % Bool
iseven(a::ZZRingElem) = ccall((:fmpz_is_even, libflint), Cint, (Ref{ZZRingElem},), a) % Bool

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

# ZZRingElem is allowed as a leaf, and the following code is needed by AA's api
expressify(x::ZZRingElem; context = nothing) = x

function AbstractAlgebra.get_syntactic_sign_abs(obj::ZZRingElem)
  return obj < 0 ? (-1, -obj) : (1, obj)
end

AbstractAlgebra.is_syntactic_one(x::ZZRingElem) = isone(x)

AbstractAlgebra.is_syntactic_zero(x::ZZRingElem) = iszero(x)

function AbstractAlgebra.print_obj(S::AbstractAlgebra.printer, mi::MIME,
    obj::ZZRingElem, left::Int, right::Int)
  AbstractAlgebra.print_integer_string(S, mi, string(obj), left, right)
end

string(x::ZZRingElem) = dec(x)

show(io::IO, x::ZZRingElem) = print(io, string(x))

function show(io::IO, a::ZZRing)
  # deliberately no @show_name or @show_special here as this is a singleton type
  if is_terse(io)
    print(pretty(io), LowercaseOff(), "ZZ")
  else
    print(io, "Integer ring")
  end
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(x::ZZRingElem) = x < 0 ? ZZRingElem(-1) : ZZRingElem(1)

###############################################################################
#
#   Unary operators and functions, e.g. -ZZRingElem(12), ~ZZRingElem(12)
#
###############################################################################

function -(x::ZZRingElem)
  z = ZZRingElem()
  neg!(z, x)
  return z
end

function ~(x::ZZRingElem)
  z = ZZRingElem()
  ccall((:fmpz_complement, libflint), Nothing, (Ref{ZZRingElem}, Ref{ZZRingElem}), z, x)
  return z
end

function abs(x::ZZRingElem)
  z = ZZRingElem()
  ccall((:fmpz_abs, libflint), Nothing, (Ref{ZZRingElem}, Ref{ZZRingElem}), z, x)
  return z
end

floor(x::ZZRingElem) = x
ceil(x::ZZRingElem) = x
trunc(x::ZZRingElem) = x
round(x::ZZRingElem) = x

floor(::Type{ZZRingElem}, x::ZZRingElem) = x
ceil(::Type{ZZRingElem}, x::ZZRingElem) = x
trunc(::Type{ZZRingElem}, x::ZZRingElem) = x
round(::Type{ZZRingElem}, x::ZZRingElem) = x

conj(x::ZZRingElem) = x

###############################################################################
#
#   Binary operators and functions
#
###############################################################################

# Metaprogram to define functions +, -, *, gcd, lcm,
#                                 &, |, $ (xor)

for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
                 (:&, :and), (:|, :or), (:xor, :xor))
  @eval begin
    function ($fJ)(x::ZZRingElem, y::ZZRingElem)
      z = ZZRingElem()
      ccall(($(string(:fmpz_, fC)), libflint), Nothing,
            (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
      return z
    end
  end
end

# Metaprogram to define functions fdiv, cdiv, tdiv, div

for (fJ, fC) in ((:fdiv, :fdiv_q), (:cdiv, :cdiv_q), (:tdiv, :tdiv_q),
                 (:div, :fdiv_q))
  @eval begin
    function ($fJ)(x::ZZRingElem, y::ZZRingElem)
      iszero(y) && throw(DivideError())
      z = ZZRingElem()
      ccall(($(string(:fmpz_, fC)), libflint), Nothing,
            (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
      return z
    end
  end
end

# N.B. we do not export the internal definition of div
# which agrees with the internal definition of AbstractAlgebra
# Here we set Base.div to a version that agrees with Base
function Base.div(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z = ZZRingElem()
  ccall((:fmpz_tdiv_q, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return z
end

function divexact(x::ZZRingElem, y::ZZRingElem; check::Bool=true)
  iszero(y) && throw(DivideError())
  z = ZZRingElem()
  if check
    r = ZZRingElem()
    ccall((:fmpz_tdiv_qr, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, r, x, y)
    r != 0 && throw(ArgumentError("Not an exact division"))
  else
    ccall((:fmpz_divexact, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  end
  return z
end

function divides(x::ZZRingElem, y::ZZRingElem)
  z = ZZRingElem()
  res = ccall((:fmpz_divides, libflint), Bool,
              (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return res, z
end

divides(x::ZZRingElem, y::Integer) = divides(x, ZZRingElem(y))

@doc raw"""
    is_divisible_by(x::ZZRingElem, y::ZZRingElem)

Return `true` if $x$ is divisible by $y$, otherwise return `false`.
"""
function is_divisible_by(x::ZZRingElem, y::ZZRingElem)
  if iszero(x)
    return true
  elseif iszero(y)
    return false
  elseif iseven(y) && isodd(x)
    return false
  elseif nbits(y) > nbits(x)
    return false
  else
    flag, q = divides(x, y)
    return flag
  end
end

@doc raw"""
    is_divisible_by(x::ZZRingElem, y::ZZRingElem)

Return `true` if $x$ is divisible by $y$, otherwise return `false`.
"""
function is_divisible_by(x::ZZRingElem, y::Integer)
  if iszero(x)
    return true
  elseif iszero(y)
    return false
  elseif iseven(y) && isodd(x)
    return false
  elseif ndigits(y, base=2) > nbits(x)
    return false
  else
    r = mod(x, y)
    return r == 0
  end
end

function is_divisible_by(x::Integer, y::ZZRingElem)
  return is_divisible_by(ZZRingElem(x), y)
end

function rem(x::ZZRingElem, c::ZZRingElem)
  iszero(c) && throw(DivideError())
  q = ZZRingElem()
  r = ZZRingElem()
  ccall((:fmpz_tdiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), q, r, x, c)
  return r
end

function rem(a::ZZRingElem, b::UInt)
  return ccall((:fmpz_fdiv_ui, libflint), UInt, (Ref{ZZRingElem}, UInt), a, b)
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function +(x::ZZRingElem, c::Int)
  z = ZZRingElem()
  add!(z, x, c)
  return z
end

+(c::Int, x::ZZRingElem) = x + c

function -(x::ZZRingElem, c::Int)
  z = ZZRingElem()
  sub!(z, x, c)
  return z
end

function -(c::Int, x::ZZRingElem)
  z = ZZRingElem()
  sub!(z, x, c)
  neg!(z, z)
  return z
end

function *(x::ZZRingElem, c::Int)
  z = ZZRingElem()
  mul!(z, x, c)
  return z
end

*(c::Int, x::ZZRingElem) = x * c

+(a::ZZRingElem, b::Integer) = a + ZZRingElem(b)

+(a::Integer, b::ZZRingElem) = ZZRingElem(a) + b

-(a::ZZRingElem, b::Integer) = a - ZZRingElem(b)

-(a::Integer, b::ZZRingElem) = ZZRingElem(a) - b

*(a::ZZRingElem, b::Integer) = a*ZZRingElem(b)

*(a::Integer, b::ZZRingElem) = ZZRingElem(a)*b

*(a::ZZRingElem, b::AbstractFloat) = BigInt(a) * b
*(a::AbstractFloat, b::ZZRingElem) = a * BigInt(b)

/(a::AbstractFloat, b::ZZRingElem) = a / BigInt(b)

###############################################################################
#
#   Rounding
#
###############################################################################

for sym in (:trunc, :round, :ceil, :floor)
  @eval begin
    # support `trunc(ZZRingElem, 1.23)` etc. for arbitrary reals
    Base.$sym(::Type{ZZRingElem}, a::Real) = ZZRingElem(Base.$sym(BigInt, a))
    Base.$sym(::Type{ZZRingElem}, a::Rational) = ZZRingElem(Base.$sym(BigInt, a))
    Base.$sym(::Type{ZZRingElem}, a::Rational{T}) where T = ZZRingElem(Base.$sym(BigInt, a))
    Base.$sym(::Type{ZZRingElem}, a::Rational{Bool}) = ZZRingElem(Base.$sym(BigInt, a))

    # for integers we don't need to round in between
    Base.$sym(::Type{ZZRingElem}, a::Integer) = ZZRingElem(a)

    # support `trunc(ZZRingElem, m)` etc. where m is a matrix of reals
    function Base.$sym(::Type{ZZMatrix}, a::Matrix{<:Real})
      s = Base.size(a)
      m = zero_matrix(ZZ, s[1], s[2])
      for i = 1:s[1], j = 1:s[2]
        m[i, j] = Base.$sym(ZZRingElem, a[i, j])
      end
      return m
    end

    # rounding QQFieldElem to integer via ZZRingElem
    function Base.$sym(::Type{T}, a::QQFieldElem) where T <: Integer
      return T(Base.$sym(ZZRingElem, a))
    end
  end
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

divexact(x::ZZRingElem, y::Integer; check::Bool=true) = divexact(x, ZZRingElem(y); check=check)

divexact(x::Integer, y::ZZRingElem; check::Bool=true) = divexact(ZZRingElem(x), y; check=check)

###############################################################################
#
#   Ad hoc division
#
###############################################################################

function tdivpow2(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Exponent must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_tdiv_q_2exp, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

function fdivpow2(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Exponent must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_fdiv_q_2exp, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

function fmodpow2(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Exponent must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_fdiv_r_2exp, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

function cdivpow2(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Exponent must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_cdiv_q_2exp, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

function tdiv(x::ZZRingElem, c::Int)
  c == 0 && throw(DivideError())
  z = ZZRingElem()
  ccall((:fmpz_tdiv_q_si, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

function fdiv(x::ZZRingElem, c::Int)
  c == 0 && throw(DivideError())
  z = ZZRingElem()
  ccall((:fmpz_fdiv_q_si, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

function cdiv(x::ZZRingElem, c::Int)
  c == 0 && throw(DivideError())
  z = ZZRingElem()
  ccall((:fmpz_cdiv_q_si, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

rem(x::Integer, y::ZZRingElem) = rem(ZZRingElem(x), y)

rem(x::ZZRingElem, y::Integer) = rem(x, ZZRingElem(y))

mod(x::Integer, y::ZZRingElem) = mod(ZZRingElem(x), y)

@doc raw"""
    mod(x::ZZRingElem, y::Integer)

Return the remainder after division of $x$ by $y$. The remainder will be
closer to zero than $y$ and have the same sign, or it will be zero.
"""
mod(x::ZZRingElem, y::Integer) = mod(x, ZZRingElem(y))

div(x::Integer, y::ZZRingElem) = div(ZZRingElem(x), y)

# Note Base.div is different to Nemo.div
Base.div(x::Integer, y::ZZRingElem) = Base.div(ZZRingElem(x), y)

div(x::ZZRingElem, y::Integer) = div(x, ZZRingElem(y))

# Note Base.div is different to Nemo.div
Base.div(x::ZZRingElem, y::Integer) = Base.div(x, ZZRingElem(y))

divrem(x::ZZRingElem, y::Integer) = divrem(x, ZZRingElem(y))

divrem(x::Integer, y::ZZRingElem) = divrem(ZZRingElem(x), y)

Base.divrem(x::ZZRingElem, y::Int) = (Base.div(x, y), Base.rem(x, y))

###############################################################################
#
#   Division with remainder
#
###############################################################################

function divrem(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z1 = ZZRingElem()
  z2 = ZZRingElem()
  ccall((:fmpz_fdiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z1, z2, x, y)
  z1, z2
end

# N.B. Base.divrem differs from Nemo.divrem
function Base.divrem(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z1 = ZZRingElem()
  z2 = ZZRingElem()
  ccall((:fmpz_tdiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z1, z2, x, y)
  z1, z2
end

function tdivrem(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z1 = ZZRingElem()
  z2 = ZZRingElem()
  ccall((:fmpz_tdiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z1, z2, x, y)
  z1, z2
end

function fdivrem(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z1 = ZZRingElem()
  z2 = ZZRingElem()
  ccall((:fmpz_fdiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z1, z2, x, y)
  z1, z2
end

function cdivrem(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z1 = ZZRingElem()
  z2 = ZZRingElem()
  ccall((:fmpz_cdiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z1, z2, x, y)
  z1, z2
end

function ntdivrem(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z1 = ZZRingElem()
  z2 = ZZRingElem()
  ccall((:fmpz_ndiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z1, z2, x, y)
  return z1, z2
end

function nfdivrem(a::ZZRingElem, b::ZZRingElem)
  q, r = tdivrem(a, b)
  c = ccall((:fmpz_cmp2abs, libflint), Cint, (Ref{ZZRingElem}, Ref{ZZRingElem}), b, r)
  if c <= 0
    if sign(Int, b) != sign(Int, r)
      sub!(q, q, UInt(1))
      add!(r, r, b)
    elseif c < 0
      add!(q, q, UInt(1))
      sub!(r, r, b)
    end
  end
  return (q, r)
end

function ncdivrem(a::ZZRingElem, b::ZZRingElem)
  q, r = tdivrem(a, b)
  c = ccall((:fmpz_cmp2abs, libflint), Cint, (Ref{ZZRingElem}, Ref{ZZRingElem}), b, r)
  if c <= 0
    if sign(Int, b) == sign(Int, r)
      add!(q, q, UInt(1))
      sub!(r, r, b)
    elseif c < 0
      sub!(q, q, UInt(1))
      add!(r, r, b)
    end
  end
  return (q, r)
end

function ndivrem(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  z1 = ZZRingElem()
  z2 = ZZRingElem()
  ccall((:fmpz_ndiv_qr, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z1, z2, x, y)
  z1, z2
end

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(x::ZZRingElem)
  if isone(x)
    return ZZRingElem(1)
  elseif x == -1
    return ZZRingElem(-1)
  end
  iszero(x) && throw(DivideError())
  throw(ArgumentError("not a unit"))
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(x::ZZRingElem, y::Union{Int, UInt, ZZRingElem})
  if isone(x) || iszero(y)
    one(x)
  elseif x == -1
    isodd(y) ? deepcopy(x) : one(x)
  elseif y < 0
    throw(DomainError(y, "Exponent must be non-negative"))
  elseif isone(y)
    deepcopy(x)
  else
    z = ZZRingElem()
    ccall((:fmpz_pow_ui, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, UInt), z, x, UInt(y))
    z
  end
end

###############################################################################
#
#  Generic powering by an ZZRingElem
#
###############################################################################

^(a::T, n::IntegerUnion) where {T<:RingElem} = _generic_power(a, n)

function _generic_power(a, n::IntegerUnion)
  fits(Int, n) && return a^Int(n)
  if is_negative(n)
    a = inv(a)
    n = -n
  end
  r = one(parent(a))
  for b = bits(n)
    r = mul!(r, r, r)
    if b
      r = mul!(r, r, a)
    end
  end
  return r
end

###############################################################################
#
#   Comparison
#
###############################################################################

function cmp(x::ZZRingElem, y::ZZRingElem)
  Int(ccall((:fmpz_cmp, libflint), Cint,
            (Ref{ZZRingElem}, Ref{ZZRingElem}), x, y))
end

==(x::ZZRingElem, y::ZZRingElem) = cmp(x,y) == 0

<=(x::ZZRingElem, y::ZZRingElem) = cmp(x,y) <= 0

<(x::ZZRingElem, y::ZZRingElem) = cmp(x,y) < 0

function cmpabs(x::ZZRingElem, y::ZZRingElem)
  Int(ccall((:fmpz_cmpabs, libflint), Cint,
            (Ref{ZZRingElem}, Ref{ZZRingElem}), x, y))
end

isless(x::ZZRingElem, y::ZZRingElem) = x < y

isless(x::ZZRingElem, y::Integer) = x < ZZRingElem(y)

isless(x::Integer, y::ZZRingElem) = ZZRingElem(x) < y

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

function cmp(x::ZZRingElem, y::Int)
  Int(ccall((:fmpz_cmp_si, libflint), Cint, (Ref{ZZRingElem}, Int), x, y))
end

==(x::ZZRingElem, y::Int) = cmp(x,y) == 0

<=(x::ZZRingElem, y::Int) = cmp(x,y) <= 0

<(x::ZZRingElem, y::Int) = cmp(x,y) < 0

==(x::Int, y::ZZRingElem) = cmp(y,x) == 0

<=(x::Int, y::ZZRingElem) = cmp(y,x) >= 0

<(x::Int, y::ZZRingElem) = cmp(y,x) > 0

function cmp(x::ZZRingElem, y::UInt)
  Int(ccall((:fmpz_cmp_ui, libflint), Cint, (Ref{ZZRingElem}, UInt), x, y))
end

==(x::ZZRingElem, y::UInt) = cmp(x,y) == 0

<=(x::ZZRingElem, y::UInt) = cmp(x,y) <= 0

<(x::ZZRingElem, y::UInt) = cmp(x,y) < 0

==(x::UInt, y::ZZRingElem) = cmp(y,x) == 0

<=(x::UInt, y::ZZRingElem) = cmp(y,x) >= 0

<(x::UInt, y::ZZRingElem) = cmp(y,x) > 0


==(x::ZZRingElem, y::Rational) = isinteger(y) && x == numerator(y)

<=(x::ZZRingElem, y::Rational) = x*denominator(y) <= numerator(y)

<(x::ZZRingElem, y::Rational) = x*denominator(y) < numerator(y)

==(x::Rational, y::ZZRingElem) = isinteger(x) && numerator(x) == y

<=(x::Rational, y::ZZRingElem) = numerator(x) <= y*denominator(x)

<(x::Rational, y::ZZRingElem) = numerator(x) < y*denominator(x)

function cmp(a::BigFloat, b::ZZRingElem)
  if _fmpz_is_small(b)
    return Int(ccall((:mpfr_cmp_si, :libmpfr), Cint, (Ref{BigFloat}, Int), a, b.d))
  end
  return Int(ccall((:mpfr_cmp_z, :libmpfr), Cint, (Ref{BigFloat}, UInt), a, unsigned(b.d) << 2))
end

==(x::ZZRingElem, y::BigFloat) = cmp(y, x) == 0

isless(x::ZZRingElem, y::BigFloat) = cmp(y, x) > 0

==(x::BigFloat, y::ZZRingElem) = cmp(x, y) == 0

isless(x::BigFloat, y::ZZRingElem) = cmp(x, y) < 0

isless(a::Float64, b::ZZRingElem) = isless(a, BigFloat(b))
isless(a::ZZRingElem, b::Float64) = isless(BigFloat(a), b)

###############################################################################
#
#   Shifting
#
###############################################################################

@doc raw"""
    <<(x::ZZRingElem, c::Int)

Return $2^cx$ where $c \geq 0$.
"""
function <<(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Exponent must be non-negative"))
  c == 0 && return x
  z = ZZRingElem()
  ccall((:fmpz_mul_2exp, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

@doc raw"""
    >>(x::ZZRingElem, c::Int)

Return $x/2^c$, discarding any remainder, where $c \geq 0$.
"""
function >>(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Exponent must be non-negative"))
  c == 0 && return x
  z = ZZRingElem()
  ccall((:fmpz_fdiv_q_2exp, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, c)
  return z
end

###############################################################################
#
#   Modular arithmetic
#
###############################################################################

function mod!(r::ZZRingElem, x::ZZRingElem, y::ZZRingElem)
  ccall((:fmpz_fdiv_r, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), r, x, y)
  return r
end

function mod(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  r = ZZRingElem()
  return mod!(r, x, y)
end

function mod(x::ZZRingElem, c::UInt)
  c == 0 && throw(DivideError())
  ccall((:fmpz_fdiv_ui, libflint), UInt, (Ref{ZZRingElem}, UInt), x, c)
end

function mod_sym(a::ZZRingElem, b::ZZRingElem)
  return mod_sym!(deepcopy(a), b)
end

function mod_sym!(a::ZZRingElem, b::ZZRingElem)
  mod!(a, a, b)
  if (b > 0 && a > div(b, 2)) || (b < 0 && a < div(b, 2))
    sub!(a, a, b)
  end
  return a
end


@doc raw"""
    powermod(x::ZZRingElem, p::ZZRingElem, m::ZZRingElem)

Return $x^p (\mod m)$. The remainder will be in the range $[0, m)$
"""
function powermod(x::ZZRingElem, p::ZZRingElem, m::ZZRingElem)
  m <= 0 && throw(DomainError(m, "Exponent must be non-negative"))
  if p < 0
    x = invmod(x, m)
    p = -p
  end
  r = ZZRingElem()
  ccall((:fmpz_powm, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}),
        r, x, p, m)
  return r
end

@doc raw"""
    powermod(x::ZZRingElem, p::Int, m::ZZRingElem)

Return $x^p (\mod m)$. The remainder will be in the range $[0, m)$
"""
function powermod(x::ZZRingElem, p::Int, m::ZZRingElem)
  m <= 0 && throw(DomainError(m, "Exponent must be non-negative"))
  if p < 0
    x = invmod(x, m)
    p = -p
  end
  r = ZZRingElem()
  ccall((:fmpz_powm_ui, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int, Ref{ZZRingElem}),
        r, x, p, m)
  return r
end

#square-and-multiply algorithm to compute f^e mod g
function powermod(f::T, e::ZZRingElem, g::T) where {T}
  #small exponent -> use powermod
  if nbits(e) <= 63
    return powermod(f, Int(e), g)
  else
    #go through binary representation of exponent and multiply with res
    #or (res and f)
    res = parent(f)(1)
    for b = bits(e)
      res = mod(res^2, g)
      if b
        res = mod(res * f, g)
      end
    end
    return res
  end
end

@doc raw"""
    invmod(x::ZZRingElem, m::ZZRingElem)

Return $x^{-1} (\mod m)$. The remainder will be in the range $[0, m)$
"""
function invmod(x::ZZRingElem, m::ZZRingElem)
  m <= 0 && throw(DomainError(m, "Modulus must be non-negative"))
  z = ZZRingElem()
  if isone(m)
    return ZZRingElem(0)
  end
  if ccall((:fmpz_invmod, libflint), Cint,
           (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, m) == 0
    error("Impossible inverse in invmod")
  end
  return z
end

@doc raw"""
    sqrtmod(x::ZZRingElem, m::ZZRingElem)

Return a square root of $x (\mod m)$ if one exists. The remainder will be in
the range $[0, m)$. We require that $m$ is prime, otherwise the algorithm may
not terminate.

# Examples

```jldoctest
julia> sqrtmod(ZZ(12), ZZ(13))
5
```
"""
function sqrtmod(x::ZZRingElem, m::ZZRingElem)
  m <= 0 && throw(DomainError(m, "Modulus must be non-negative"))
  z = ZZRingElem()
  if (ccall((:fmpz_sqrtmod, libflint), Cint,
            (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, m) == 0)
    error("no square root exists or modulus is not prime")
  end
  return z
end

function _normalize_crt(r::ZZRingElem, m::ZZRingElem, signed)
  s = sign(Int, m)
  if s > 0
    return signed ? nfdivrem(r, m)[2] : fdivrem(r, m)[2]
  elseif s < 0
    return signed ? ncdivrem(r, m)[2] : cdivrem(r, m)[2]
  else
    return r
  end
end

function _normalize_crt_with_lcm(r::ZZRingElem, m::ZZRingElem, signed)
  s = sign(Int, m)
  if s == 0
    return (r, m)
  elseif s < 0
    m = -m
  end
  return (signed ? nfdivrem(r, m)[2] : fdivrem(r, m)[2], m)
end

@doc raw"""
    crt(r1::ZZRingElem, m1::ZZRingElem, r2::ZZRingElem, m2::ZZRingElem, signed=false; check::Bool=true)
    crt(r1::ZZRingElem, m1::ZZRingElem, r2::Union{Int, UInt}, m2::Union{Int, UInt}, signed=false; check::Bool=true)
    crt(r::Vector{ZZRingElem}, m::Vector{ZZRingElem}, signed=false; check::Bool=true)
    crt_with_lcm(r1::ZZRingElem, m1::ZZRingElem, r2::ZZRingElem, m2::ZZRingElem, signed=false; check::Bool=true)
    crt_with_lcm(r1::ZZRingElem, m1::ZZRingElem, r2::Union{Int, UInt}, m2::Union{Int, UInt}, signed=false; check::Bool=true)
    crt_with_lcm(r::Vector{ZZRingElem}, m::Vector{ZZRingElem}, signed=false; check::Bool=true)

As per the AbstractAlgebra `crt` interface, with the following option.
If `signed = true`, the solution is the range $(-m/2, m/2]$, otherwise it is in
the range $[0,m)$, where $m$ is the least common multiple of the moduli.

# Examples

```jldoctest
julia> crt(ZZ(5), ZZ(13), ZZ(7), ZZ(37), true)
44

julia> crt(ZZ(5), ZZ(13), 7, 37, true)
44
```
"""
function crt(r1::ZZRingElem, m1::ZZRingElem, r2::ZZRingElem, m2::ZZRingElem, signed=false; check::Bool=true)
  r, m = AbstractAlgebra._crt_with_lcm_stub(r1, m1, r2, m2; check=check)
  return _normalize_crt(r, m, signed)
end

function crt(r::Vector{ZZRingElem}, m::Vector{ZZRingElem}, signed=false; check::Bool=true)
  r, m = AbstractAlgebra._crt_with_lcm_stub(r, m; check=check)
  return _normalize_crt(r, m, signed)
end

function crt_with_lcm(r1::ZZRingElem, m1::ZZRingElem, r2::ZZRingElem, m2::ZZRingElem, signed=false; check::Bool=true)
  r, m = AbstractAlgebra._crt_with_lcm_stub(r1, m1, r2, m2; check=check)
  return _normalize_crt_with_lcm(r, m, signed)
end

function crt_with_lcm(r::Vector{ZZRingElem}, m::Vector{ZZRingElem}, signed=false; check::Bool=true)
  r, m = AbstractAlgebra._crt_with_lcm_stub(r, m; check=check)
  return _normalize_crt_with_lcm(r, m, signed)
end

# requires a < b
function _gcdinv(a::UInt, b::UInt)
  s = Ref{UInt}()
  g = ccall((:n_gcdinv, libflint), UInt,
            (Ptr{UInt}, UInt, UInt),
            s, a, b)
  return g, s[]
end

function submod(a::UInt, b::UInt, n::UInt)
  return a >= b ? a - b : a - b + n
end

function _crt_with_lcm(r1::ZZRingElem, m1::ZZRingElem, r2::UInt, m2::UInt; check::Bool=true)
  if iszero(m2)
    check && !is_divisible_by(r2 - r1, m1) && error("no crt solution")
    return (ZZRingElem(r2), ZZRingElem(m2))
  end
  if r2 >= m2
    r2 = mod(r2, m2)
  end
  if iszero(m1)
    check && mod(r1, m2) != r2 && error("no crt solution")
    return (r1, m1)
  end
  g, s = _gcdinv(mod(m1, m2), m2)
  diff = submod(r2, mod(r1, m2), m2)
  if isone(g)
    return (r1 + mulmod(diff, s, m2)*m1, m1*m2)
  else
    m2og = divexact(m2, g; check=false)
    diff = divexact(diff, g; check=check)
    return (r1 + mulmod(diff, s, m2og)*m1, m1*m2og)
  end
end

function _crt_with_lcm(r1::ZZRingElem, m1::ZZRingElem, r2::Union{Int, UInt},
    m2::Union{Int, UInt}; check::Bool=true)
  if iszero(m2)
    check && !is_divisible_by(r2 - r1, m1) && error("no crt solution")
    return (ZZRingElem(r2), ZZRingElem(m2))
  end
  m2 = abs(m2)%UInt
  r2 = r2 < 0 ? mod(r2, m2) : UInt(r2)
  return _crt_with_lcm(r1, m1, r2::UInt, m2; check=check)
end

function crt(r1::ZZRingElem, m1::ZZRingElem, r2::Union{Int, UInt},
    m2::Union{Int, UInt}, signed = false; check::Bool=true)
  r, m = _crt_with_lcm(r1, m1, r2, m2; check=check)
  return _normalize_crt(r, m, signed)
end

function crt_with_lcm(r1::ZZRingElem, m1::ZZRingElem, r2::Union{Int, UInt},
    m2::Union{Int, UInt}, signed = false; check::Bool=true)
  r, m = _crt_with_lcm(r1, m1, r2, m2; check=check)
  return _normalize_crt_with_lcm(r, m, signed)
end

###############################################################################
#
#   Integer logarithm
#
###############################################################################

@doc raw"""
    flog(x::ZZRingElem, c::ZZRingElem)
    flog(x::ZZRingElem, c::Int)

Return the floor of the logarithm of $x$ to base $c$.

# Examples

```jldoctest
julia> flog(ZZ(12), ZZ(2))
3

julia> flog(ZZ(12), 3)
2

```
"""
function flog(x::ZZRingElem, c::ZZRingElem)
  c <= 0 && throw(DomainError(c, "Base must be non-negative"))
  x <= 0 && throw(DomainError(x, "Argument must be non-negative"))
  return ccall((:fmpz_flog, libflint), Int,
               (Ref{ZZRingElem}, Ref{ZZRingElem}), x, c)
end

function flog(x::ZZRingElem, c::Int)
  c <= 0 && throw(DomainError(c, "Base must be non-negative"))
  return ccall((:fmpz_flog_ui, libflint), Int,
               (Ref{ZZRingElem}, Int), x, c)
end

@doc raw"""
    clog(x::ZZRingElem, c::ZZRingElem)
    clog(x::ZZRingElem, c::Int)

Return the ceiling of the logarithm of $x$ to base $c$.

# Examples

```jldoctest
julia> clog(ZZ(12), ZZ(2))
4

julia> clog(ZZ(12), 3)
3

```
"""
function clog(x::ZZRingElem, c::ZZRingElem)
  c <= 0 && throw(DomainError(c, "Base must be non-negative"))
  x <= 0 && throw(DomainError(x, "Argument must be non-negative"))
  return ccall((:fmpz_clog, libflint), Int,
               (Ref{ZZRingElem}, Ref{ZZRingElem}), x, c)
end

function clog(x::ZZRingElem, c::Int)
  c <= 0 && throw(DomainError(c, "Base must be non-negative"))
  return ccall((:fmpz_clog_ui, libflint), Int,
               (Ref{ZZRingElem}, Int), x, c)
end

###############################################################################
#
#   Natural logarithm
#
###############################################################################

log(a::ZZRingElem) = log(BigInt(a))
log(a::ZZRingElem, b::ZZRingElem) = log(b) / log(a)

###############################################################################
#
#   GCD and LCM
#
###############################################################################

@doc raw"""
    gcd(x::ZZRingElem, y::ZZRingElem, z::ZZRingElem...)

Return the greatest common divisor of $(x, y, ...)$. The returned result will
always be non-negative and will be zero iff all inputs are zero.
"""
function gcd(x::ZZRingElem, y::ZZRingElem, z::ZZRingElem...)
  d = ZZRingElem()
  ccall((:fmpz_gcd, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), d, x, y)
  length(z) == 0 && return d

  for ix in 1:length(z)
    ccall((:fmpz_gcd, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), d, d, z[ix])
  end
  return d
end

@doc raw"""
    gcd(x::Vector{ZZRingElem})

Return the greatest common divisor of the elements of $x$. The returned
result will always be non-negative and will be zero iff all elements of $x$
are zero.
"""
function gcd(x::Vector{ZZRingElem})
  if length(x) == 0
    error("Array must not be empty")
  elseif length(x) == 1
    return x[1]
  end

  z = ZZRingElem()
  ccall((:fmpz_gcd, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x[1], x[2])

  for i in 3:length(x)
    ccall((:fmpz_gcd, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, z, x[i])
    if isone(z)
      return z
    end
  end

  return z
end

@doc raw"""
    lcm(x::ZZRingElem, y::ZZRingElem, z::ZZRingElem...)

Return the least common multiple of $(x, y, ...)$. The returned result will
always be non-negative and will be zero if any input is zero.
"""
function lcm(x::ZZRingElem, y::ZZRingElem, z::ZZRingElem...)
  m = ZZRingElem()
  ccall((:fmpz_lcm, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), m, x, y)
  length(z) == 0 && return m

  for ix in 1:length(z)
    ccall((:fmpz_lcm, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), m, m, z[ix])
  end
  return m
end

@doc raw"""
    lcm(x::Vector{ZZRingElem})

Return the least common multiple of the elements of $x$. The returned result
will always be non-negative and will be zero iff the elements of $x$ are zero.
"""
function lcm(x::Vector{ZZRingElem})
  if length(x) == 0
    error("Array must not be empty")
  elseif length(x) == 1
    return x[1]
  end

  z = ZZRingElem()
  ccall((:fmpz_lcm, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x[1], x[2])

  for i in 3:length(x)
    ccall((:fmpz_lcm, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, z, x[i])
  end

  return z
end

gcd(a::ZZRingElem, b::Integer) = gcd(a, ZZRingElem(b))

gcd(a::Integer, b::ZZRingElem) = gcd(ZZRingElem(a), b)

lcm(a::ZZRingElem, b::Integer) = lcm(a, ZZRingElem(b))

lcm(a::Integer, b::ZZRingElem) = lcm(ZZRingElem(a), b)

###############################################################################
#
#   Extended GCD
#
###############################################################################

@doc raw"""
    gcdx(a::ZZRingElem, b::ZZRingElem)

Return a tuple $g, s, t$ such that $g$ is the greatest common divisor of $a$
and $b$ and integers $s$ and $t$ such that $g = as + bt$.
"""
function gcdx(a::ZZRingElem, b::ZZRingElem)
  @static if VERSION < v"1.12.0-DEV.410"
    # Just to conform with Julia's definition
    a == b == 0 && return zero(ZZ), one(ZZ), zero(ZZ)
  end

  d = ZZ()
  x = ZZ()
  y = ZZ()
  ccall((:fmpz_xgcd_canonical_bezout, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), d, x, y, a, b)
  return d, x, y
end

@doc raw"""
    gcdinv(a::ZZRingElem, b::ZZRingElem)

Return a tuple $g, s$ where $g$ is the greatest common divisor of $a$ and
$b$ and where $s$ is the inverse of $a$ modulo $b$ if $g = 1$. This function
can be used to detect impossible inverses, i.e. where $a$ and $b$ are not
coprime, and to yield the common factor of $a$ and $b$ if they are not
coprime. We require $b \geq a \geq 0$.
"""
function gcdinv(a::ZZRingElem, b::ZZRingElem)
  a < 0 && throw(DomainError(a, "First argument must be non-negative"))
  b < a && throw(DomainError((a, b), "First argument must be smaller than second argument"))
  g = ZZRingElem()
  s = ZZRingElem()
  ccall((:fmpz_gcdinv, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}),
        g, s, a, b)
  return g, s
end

gcdx(a::ZZRingElem, b::Integer) = gcdx(a, ZZRingElem(b))

gcdx(a::Integer, b::ZZRingElem) = gcdx(ZZRingElem(a), b)

gcdinv(a::ZZRingElem, b::Integer) = gcdinv(a, ZZRingElem(b))

gcdinv(a::Integer, b::ZZRingElem) = gcdinv(ZZRingElem(a), b)

###############################################################################
#
#   Roots
#
###############################################################################

const sqrt_moduli = [3, 5, 7, 8]
const sqrt_residues = [[0, 1], [0, 1, 4], [0, 1, 2, 4], [0, 1, 4]]

@doc raw"""
    isqrt(x::ZZRingElem)

Return the floor of the square root of $x$.

# Examples

```jldoctest
julia> isqrt(ZZ(13))
3

```
"""
function isqrt(x::ZZRingElem)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_sqrt, libflint), Nothing, (Ref{ZZRingElem}, Ref{ZZRingElem}), z, x)
  return z
end

@doc raw"""
    isqrtrem(x::ZZRingElem)

Return a tuple $s, r$ consisting of the floor $s$ of the square root of $x$
and the remainder $r$, i.e. such that $x = s^2 + r$. We require $x \geq 0$.

# Examples

```jldoctest
julia> isqrtrem(ZZ(13))
(3, 4)

```
"""
function isqrtrem(x::ZZRingElem)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  s = ZZRingElem()
  r = ZZRingElem()
  ccall((:fmpz_sqrtrem, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), s, r, x)
  return s, r
end

function Base.sqrt(x::ZZRingElem; check=true)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  if check
    for i = 1:length(sqrt_moduli)
      res = mod(x, sqrt_moduli[i])
      !(res in sqrt_residues[i]) && error("Not a square")
    end
    s = ZZRingElem()
    r = ZZRingElem()
    ccall((:fmpz_sqrtrem, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), s, r, x)
    !iszero(r) && error("Not a square")
  else
    s = ZZRingElem()
    ccall((:fmpz_sqrt, libflint), Nothing, (Ref{ZZRingElem}, Ref{ZZRingElem}), s, x)
  end
  return s
end

is_square(x::ZZRingElem) = Bool(ccall((:fmpz_is_square, libflint), Cint,
                                      (Ref{ZZRingElem},), x))

function is_square_with_sqrt(x::ZZRingElem)
  if x < 0
    return false, zero(ZZRingElem)
  end
  for i = 1:length(sqrt_moduli)
    res = mod(x, sqrt_moduli[i])
    if !(res in sqrt_residues[i])
      return false, zero(ZZRingElem)
    end
  end
  s = ZZRingElem()
  r = ZZRingElem()
  ccall((:fmpz_sqrtrem, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), s, r, x)
  if !iszero(r)
    return false, zero(ZZRingElem)
  end
  return true, s
end

@doc raw"""
    root(x::ZZRingElem, n::Int; check::Bool=true)

Return the $n$-the root of $x$. We require $n > 0$ and that
$x \geq 0$ if $n$ is even. By default the function tests whether the input was
a perfect $n$-th power and if not raises an exception. If `check=false` this
check is omitted.

# Examples

```jldoctest
julia> root(ZZ(27), 3; check=true)
3
```
"""
function root(x::ZZRingElem, n::Int; check::Bool=true)
  x < 0 && iseven(n) && throw(DomainError((x, n), "Argument `x` must be positive if exponent `n` is even"))
  n <= 0 && throw(DomainError(n, "Exponent must be positive"))
  z = ZZRingElem()
  res = ccall((:fmpz_root, libflint), Bool,
              (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, n)
  check && !res && error("Not a perfect n-th power (n = $n)")
  return z
end

@doc raw"""
    iroot(x::ZZRingElem, n::Int)

Return the integer truncation of the $n$-the root of $x$ (round towards zero).
We require $n > 0$ and that $x \geq 0$ if $n$ is even.

# Examples

```jldoctest
julia> iroot(ZZ(13), 3)
2
```
"""
function iroot(x::ZZRingElem, n::Int)
  x < 0 && iseven(n) && throw(DomainError((x, n), "Argument `x` must be positive if exponent `n` is even"))
  n <= 0 && throw(DomainError(n, "Exponent must be positive"))
  z = ZZRingElem()
  ccall((:fmpz_root, libflint), Bool,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, n)
  return z
end

#TODO (Hard): Implement this properly.
@doc raw"""
    is_squarefree(n::Union{Int, ZZRingElem}) -> Bool

Return `true` if $n$ is squarefree, `false` otherwise.
"""
function is_squarefree(n::Union{Int, ZZRingElem})
  iszero(n) && return false
  is_unit(n) && return true
  e, b = is_perfect_power_with_data(n)
  if e > 1
    return false
  end
  return isone(maximum(values(factor(n).fac); init = 1))
end
################################################################################
#
#   Factor trial range
#
################################################################################

function _factor_trial_range(N::ZZRingElem, start::Int = 0, np::Int = 10^5)
  F = fmpz_factor()
  ccall((:fmpz_factor_trial_range, libflint), Nothing,
        (Ref{fmpz_factor}, Ref{ZZRingElem}, UInt, UInt), F, N, start, np)
  res = Dict{ZZRingElem, Int}()
  for i in 1:F.num
    z = ZZRingElem()
    ccall((:fmpz_factor_get_fmpz, libflint), Nothing,
          (Ref{ZZRingElem}, Ref{fmpz_factor}, Int), z, F, i - 1)
    res[z] = unsafe_load(F.exp, i)
  end
  return res, canonical_unit(N)
end

###############################################################################
#
#   Number theoretic/combinatorial
#
###############################################################################

@doc raw"""
    divisors(a::Union{Int, ZZRingElem})

Return the positive divisors of $a$ in an array, not necessarily in ascending
order. We require $a \neq 0$.
"""
function divisors end

function divisors(a::ZZRingElem)
  iszero(a) && throw(DomainError("Argument must be non-zero"))

  divs = ZZRingElem[one(ZZ)]
  isone(a) && return divs

  for (p,e) in factor(a)
    ndivs = copy(divs)
    for i = 1:e
      map!(d -> p*d, ndivs, ndivs)
      append!(divs, ndivs)
    end
  end

  return divs
end

divisors(a::Int) = Int.(divisors(ZZ(a)))

@doc raw"""
    prime_divisors(a::ZZRingElem)

Return the prime divisors of $a$ in an array. We require $a \neq 0$.
"""
function prime_divisors(a::ZZRingElem)
  iszero(a) && throw(DomainError("Argument must be non-zero"))
  ZZRingElem[p for (p, e) in factor(a)]
end

@doc raw"""
    prime_divisors(a::Int)

Return the prime divisors of $a$ in an array. We require $a \neq 0$.
"""
prime_divisors(a::Int) = Int.(prime_divisors(ZZ(a)))

is_prime(x::UInt) = Bool(ccall((:n_is_prime, libflint), Cint, (UInt,), x))

@doc raw"""
    is_prime(x::ZZRingElem)
    is_prime(x::Int)

Return `true` if $x$ is a prime number, otherwise return `false`.

# Examples

```jldoctest
julia> is_prime(ZZ(13))
true
```
"""
function is_prime(x::ZZRingElem)
  !is_probable_prime(x) && return false
  return Bool(ccall((:fmpz_is_prime, libflint), Cint, (Ref{ZZRingElem},), x))
end

function is_prime(n::Int)
  if n < 0
    return false
  end
  return is_prime(n % UInt)
end

@doc raw"""
    is_probable_prime(x::ZZRingElem)

Return `true` if $x$ is very probably a prime number, otherwise return
`false`. No counterexamples are known to this test, but it is conjectured
that infinitely many exist.
"""
is_probable_prime(x::ZZRingElem) = Bool(ccall((:fmpz_is_probabprime, libflint), Cint,
                                              (Ref{ZZRingElem},), x))

@doc raw"""
    next_prime(x::ZZRingElem, proved = true)

Return the smallest prime strictly greater than $x$.
If a second argument of `false` is specified, the return is only probably prime.
"""
function next_prime(x::ZZRingElem, proved::Bool = true)
  z = ZZRingElem()
  ccall((:fmpz_nextprime, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Cint),
        z, x, proved)
  return z
end

function next_prime(x::UInt, proved::Bool = true)
  if (Base.GMP.BITS_PER_LIMB == 64 && x >= 0xffffffffffffffc5) ||
    (Base.GMP.BITS_PER_LIMB == 32 && x >= 0xfffffffb)
    error("No larger single-limb prime exists")
  end
  return ccall((:n_nextprime, libflint), UInt,
               (UInt, Cint),
               x, proved)
end

function next_prime(x::Int, proved::Bool = true)
  return x < 2 ? 2 : Int(next_prime(x % UInt, proved))
end

function remove!(a::ZZRingElem, b::ZZRingElem)
  v = ccall((:fmpz_remove, libflint), Clong, (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), a, a, b)
  return v, a
end

function remove(x::ZZRingElem, y::ZZRingElem)
  iszero(y) && throw(DivideError())
  y <= 1 && error("Factor <= 1")
  z = ZZRingElem()
  num = ccall((:fmpz_remove, libflint), Int,
              (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return num, z
end

remove(x::ZZRingElem, y::Integer) = remove(x, ZZRingElem(y))

remove(x::Integer, y::ZZRingElem) = remove(ZZRingElem(x), y)

function remove(a::UInt, b::UInt)
  b <= 1 && error("Factor <= 1")
  a == 0 && error("Not yet implemented")
  q = Ref(a)
  binv = ccall((:n_precompute_inverse, libflint), Float64, (UInt,), b)
  v = ccall((:n_remove2_precomp, libflint), Cint,
            (Ptr{UInt}, UInt, Float64),
            q, b, binv)
  return (Int(v), q[])
end

function remove(a::Int, b::Int)
  b <= 1 && error("Factor <= 1")
  v, q = remove(abs(a)%UInt, b%UInt)
  return (v, a < 0 ? -q%Int : q%Int)
end

function remove(a::BigInt, b::BigInt)
  b <= 1 && error("Factor <= 1")
  a == 0 && error("Not yet implemented")
  q = BigInt()
  v =  ccall((:__gmpz_remove, :libgmp), Culong,
             (Ref{BigInt}, Ref{BigInt}, Ref{BigInt}),
             q, a, b)
  return (Int(v), q)
end

function remove(x::Integer, y::Integer)
  v, q = remove(ZZRingElem(x), ZZRingElem(y))
  return (v, convert(promote_type(typeof(x), typeof(y)), q))
end

@doc raw"""
    valuation(x::ZZRingElem, y::ZZRingElem)

Return the largest $n$ such that $y^n$ divides $x$.
"""
function valuation(x::ZZRingElem, y::ZZRingElem)
  iszero(x) && error("Not yet implemented")
  n, _ = remove(x, y)
  return n
end

valuation(x::ZZRingElem, y::Integer) = valuation(x, ZZRingElem(y))

valuation(x::Integer, y::ZZRingElem) = valuation(ZZRingElem(x), y)

valuation(x::Integer, y::Integer) = valuation(ZZRingElem(x), ZZRingElem(y))

@doc raw"""
    divisor_lenstra(n::ZZRingElem, r::ZZRingElem, m::ZZRingElem)

If $n$ has a factor which lies in the residue class $r (\mod m)$ for
$0 < r < m < n$, this function returns such a factor. Otherwise it returns
$0$. This is only efficient if $m$ is at least the cube root of $n$. We
require gcd$(r, m) = 1$ and this condition is not checked.
"""
function divisor_lenstra(n::ZZRingElem, r::ZZRingElem, m::ZZRingElem)
  r <= 0 && throw(DomainError(r, "Residue class must be non-negative"))
  m <= r && throw(DomainError((m, r), "Modulus must be bigger than residue class"))
  n <= m && throw(DomainError((n, m), "Argument must be bigger than modulus"))
  z = ZZRingElem()
  if !Bool(ccall((:fmpz_divisor_in_residue_class_lenstra, libflint),
                 Cint, (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, n, r, m))
    z = 0
  end
  return z
end

@doc raw"""
    factorial(x::ZZRingElem)

Return the factorial of $x$, i.e. $x! = 1.2.3\ldots x$. We require
$x \geq 0$.

# Examples

```jldoctest
julia> factorial(ZZ(100))
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
```
"""
function factorial(x::ZZRingElem)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_fac_ui, libflint), Nothing, (Ref{ZZRingElem}, UInt), z, UInt(x))
  return z
end

@doc raw"""
    rising_factorial(x::ZZRingElem, n::Int)

Return the rising factorial of $x$, i.e. $x(x + 1)(x + 2)\ldots (x + n - 1)$.
If $n < 0$ we throw a `DomainError()`.
"""
function rising_factorial(x::ZZRingElem, n::Int)
  n < 0 && throw(DomainError(n, "Argument must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_rfac_ui, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, UInt), z, x, UInt(n))
  return z
end

@doc raw"""
    rising_factorial(x::ZZRingElem, n::ZZRingElem)

Return the rising factorial of $x$, i.e. $x(x + 1)(x + 2)\cdots (x + n - 1)$.
If $n < 0$ we throw a `DomainError()`.
"""
rising_factorial(x::ZZRingElem, n::ZZRingElem) = rising_factorial(x, Int(n))

@doc raw"""
    primorial(x::Int)

Return the primorial of $x$, i.e. the product of all primes less than or
equal to $x$. If $x < 0$ we throw a `DomainError()`.
"""
function primorial(x::Int)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  # Up to 28 is OK for Int32, up to 52 is OK for Int64, up to 100 is OK for Int128; beyond is too large
  if (Int == Int32 && x > 28) || (Int == Int64 && x > 52) || (Int == Int128 && x > 100)
    throw(OverflowError("primorial(::Int)"))
  end
  z = ZZRingElem()
  ccall((:fmpz_primorial, libflint), Nothing,
        (Ref{ZZRingElem}, UInt), z, UInt(x))
  return Int(z)
end

@doc raw"""
    primorial(x::ZZRingElem)

Return the primorial of $x$, i.e. the product of all primes less than or
equal to $x$. If $x < 0$ we throw a `DomainError()`.
"""
function primorial(x::ZZRingElem)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_primorial, libflint), Nothing,
        (Ref{ZZRingElem}, UInt), z, UInt(x))
  return z
end

@doc raw"""
    fibonacci(x::Int)

Return the $x$-th Fibonacci number $F_x$. We define $F_1 = 1$, $F_2 = 1$ and
$F_{i + 1} = F_i + F_{i - 1}$ for all integers $i$.
"""
function fibonacci(x::Int)
  # Up to 46 is OK for Int32; up to 92 is OK for Int64; up to 184 is OK for Int128; beyond is too large
  if (Int == Int32 && abs(x) > 46) || (Int == Int64 && abs(x) > 92) || (Int == Int128 && abs(x) > 184)
    throw(OverflowError("fibonacci(::Int)"))
  end
  z = ZZRingElem()
  ccall((:fmpz_fib_ui, libflint), Nothing,
        (Ref{ZZRingElem}, UInt), z, UInt(abs(x)))
  return x < 0 ? (iseven(x) ? -Int(z) : Int(z)) : Int(z)
end

@doc raw"""
    fibonacci(x::ZZRingElem)

Return the $x$-th Fibonacci number $F_x$. We define $F_1 = 1$, $F_2 = 1$ and
$F_{i + 1} = F_i + F_{i - 1}$ for all integers $i$.
"""
function fibonacci(x::ZZRingElem)
  z = ZZRingElem()
  ccall((:fmpz_fib_ui, libflint), Nothing,
        (Ref{ZZRingElem}, UInt), z, UInt(abs(x)))
  return x < 0 ? (iseven(x) ? -z : z) : z
end

@doc raw"""
    bell(x::Int)

Return the Bell number $B_x$.
"""
function bell(x::Int)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  z = ZZRingElem()
  ccall((:arith_bell_number, libflint), Nothing,
        (Ref{ZZRingElem}, UInt), z, UInt(x))
  return Int(z)
end

@doc raw"""
    bell(x::ZZRingElem)

Return the Bell number $B_x$.
"""
function bell(x::ZZRingElem)
  x < 0 && throw(DomainError(x, "Argument must be non-negative"))
  z = ZZRingElem()
  ccall((:arith_bell_number, libflint), Nothing,
        (Ref{ZZRingElem}, UInt), z, UInt(x))
  return z
end

# fmpz_bin_uiui doesn't always work on UInt input as it just wraps mpz_bin_uiui,
# which has silly gnu problems on windows
# TODO: fib_ui, pow_ui, fac_ui ditto
function _binomial!(z::ZZRingElem, n::Culong, k::Culong)
  ccall((:fmpz_bin_uiui, libflint), Nothing,
        (Ref{ZZRingElem}, UInt, UInt), z, UInt(n), UInt(k))
  return z
end

function _binomial(n::ZZRingElem, k::ZZRingElem)
  @assert k >= 0
  z = ZZRingElem(1)
  if fits(Culong, n) && fits(Culong, k)
    _binomial!(z, Culong(n), Culong(k))
  elseif fits(UInt, k)
    for K in UInt(1):UInt(k)
      mul!(z, z, n - (K - 1))
      divexact!(z, z, K)
    end
  else
    # if called with n >= k
    error("result of binomial($n, $k) probably is too large")
  end
  return z
end


@doc raw"""
    binomial(n::ZZRingElem, k::ZZRingElem)

Return the binomial coefficient $\frac{n (n-1) \cdots (n-k+1)}{k!}$.
If $k < 0$ we return $0$, and the identity
`binomial(n, k) == binomial(n - 1, k - 1) + binomial(n - 1, k)` always holds
for integers `n` and `k`.
"""
function binomial(n::ZZRingElem, k::ZZRingElem)
  ksgn = cmp(k, 0)
  ksgn > 0 || return ZZRingElem(ksgn == 0)
  # k > 0 now
  negz = false
  if n < 0
    n = k - n - 1
    negz = isodd(k)
  end
  if n < k
    z = ZZRingElem(0)
  elseif 2*k <= n
    z = _binomial(n, k)
  else
    z = _binomial(n, n - k)
  end
  return negz ? neg!(z) : z
end

@doc raw"""
    binomial(n::UInt, k::UInt, ::ZZRing)

Return the binomial coefficient $\frac{n!}{(n - k)!k!}$ as an `ZZRingElem`.
"""
function binomial(n::UInt, k::UInt, ::ZZRing)
  z = ZZRingElem()
  ccall((:fmpz_bin_uiui, libflint), Nothing,
        (Ref{ZZRingElem}, UInt, UInt), z, n, k)
  return z
end

@doc raw"""
    moebius_mu(x::ZZRingElem)

Return the Moebius mu function of $x$ as an `Int`. The value
returned is either $-1$, $0$ or $1$. If $x \leq 0$ we throw a `DomainError()`.
"""
function moebius_mu(x::ZZRingElem)
  x <= 0 && throw(DomainError(x, "Argument must be positive"))
  return Int(ccall((:fmpz_moebius_mu, libflint), Cint,
                   (Ref{ZZRingElem},), x))
end

@doc raw"""
    moebius_mu(x::Int)

Return the Moebius mu function of $x$ as an `Int`. The value
returned is either $-1$, $0$ or $1$. If $x \leq 0$ we throw a `DomainError()`.
"""
moebius_mu(x::Int) = moebius_mu(ZZRingElem(x))

@doc raw"""
    jacobi_symbol(x::ZZRingElem, y::ZZRingElem)

Return the value of the Jacobi symbol $\left(\frac{x}{y}\right)$. The modulus
$y$ must be odd and positive, otherwise a `DomainError` is thrown.
"""
function jacobi_symbol(x::ZZRingElem, y::ZZRingElem)
  (y <= 0 || iseven(y)) && throw(DomainError(y, "Modulus must be odd and positive"))
  if x < 0 || x >= y
    x = mod(x, y)
  end
  return Int(ccall((:fmpz_jacobi, libflint), Cint,
                   (Ref{ZZRingElem}, Ref{ZZRingElem}), x, y))
end

jacobi_symbol(x::IntegerUnion, y::IntegerUnion) = jacobi_symbol(ZZ(x), ZZ(y))

@doc raw"""
    jacobi_symbol(x::Int, y::Int)

Return the value of the Jacobi symbol $\left(\frac{x}{y}\right)$. The modulus
$y$ must be odd and positive, otherwise a `DomainError` is thrown.
"""
function jacobi_symbol(x::Int, y::Int)
  (y <= 0 || mod(y, 2) == 0) && throw(DomainError(y, "Modulus must be odd and positive"))
  if x < 0 || x >= y
    x = mod(x, y)
  end
  return Int(ccall((:n_jacobi, libflint), Cint, (Int, UInt), x, UInt(y)))
end

@doc raw"""
    kronecker_symbol(x::ZZRingElem, y::ZZRingElem)
    kronecker_symbol(x::Int, y::Int)

Return the value of the Kronecker symbol $\left(\frac{x}{y}\right)$.
The definition is as per Henri Cohen's book, "A Course in Computational
Algebraic Number Theory", Definition 1.4.8.
"""
function kronecker_symbol(x::Int, y::Int)
  return Int(ccall((:z_kronecker, libflint), Cint,
                   (Int, Int), x, y))
end

function kronecker_symbol(x::ZZRingElem, y::ZZRingElem)
  return Int(ccall((:fmpz_kronecker, libflint), Cint,
                   (Ref{ZZRingElem}, Ref{ZZRingElem}), x, y))
end

@doc raw"""
    divisor_sigma(x::ZZRingElem, y::Int)
    divisor_sigma(x::ZZRingElem, y::ZZRingElem)
    divisor_sigma(x::Int, y::Int)

Return the value of the sigma function, i.e. $\sum_{0 < d \;| x} d^y$. If
$x \leq 0$ or $y < 0$ we throw a `DomainError()`.

# Examples

```jldoctest
julia> divisor_sigma(ZZ(32), 10)
1127000493261825

julia> divisor_sigma(ZZ(32), ZZ(10))
1127000493261825

julia> divisor_sigma(32, 10)
1127000493261825
```
"""
function divisor_sigma(x::ZZRingElem, y::Int)
  x <= 0 && throw(DomainError(x, "Argument must be positive"))
  y < 0 && throw(DomainError(y, "Power must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_divisor_sigma, libflint), Nothing,
        (Ref{ZZRingElem}, UInt, Ref{ZZRingElem}), z, UInt(y), x)
  return z
end

divisor_sigma(x::ZZRingElem, y::ZZRingElem) = divisor_sigma(x, Int(y))
divisor_sigma(x::Int, y::Int) = Int(divisor_sigma(ZZRingElem(x), y))

@doc raw"""
    euler_phi(x::ZZRingElem)
    euler_phi(x::Int)

Return the value of the Euler phi function at $x$, i.e. the number of
positive integers up to $x$ (inclusive) that are coprime with $x$. An
exception is raised if $x \leq 0$.

# Examples

```jldoctest
julia> euler_phi(ZZ(12480))
3072

julia> euler_phi(12480)
3072
```
"""
function euler_phi(x::ZZRingElem)
  x <= 0 && throw(DomainError(x, "Argument must be positive"))
  z = ZZRingElem()
  ccall((:fmpz_euler_phi, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}), z, x)
  return z
end

euler_phi(x::Int) = Int(euler_phi(ZZRingElem(x)))

function euler_phi(x::Fac{ZZRingElem})
  return prod((p - 1) * p^(v - 1) for (p, v) in x)
end

@doc raw"""
    number_of_partitions(x::Int)
    number_of_partitions(x::ZZRingElem)

Return the number of partitions of $x$.

# Examples

```jldoctest
julia> number_of_partitions(100)
190569292

julia> number_of_partitions(ZZ(1000))
24061467864032622473692149727991
```
"""
function number_of_partitions(x::Int)
  if x < 0
    return 0
  end
  z = ZZRingElem()
  ccall((:partitions_fmpz_ui, libflint), Nothing,
        (Ref{ZZRingElem}, UInt), z, x)
  return Int(z)
end

function number_of_partitions(x::ZZRingElem)
  z = ZZRingElem()
  if x < 0
    return z
  end
  ccall((:partitions_fmpz_fmpz, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, 0)
  return z
end

###############################################################################
#
#   Number bases/digits
#
###############################################################################

@doc raw"""
    bin(n::ZZRingElem)

Return $n$ as a binary string.

# Examples

```jldoctest
julia> bin(ZZ(12))
"1100"
```
"""
bin(n::ZZRingElem) = base(n, 2)

@doc raw"""
    oct(n::ZZRingElem)

Return $n$ as a octal string.

# Examples

```jldoctest
julia> oct(ZZ(12))
"14"
```
"""
oct(n::ZZRingElem) = base(n, 8)

@doc raw"""
    dec(n::ZZRingElem)

Return $n$ as a decimal string.

# Examples

```jldoctest
julia> dec(ZZ(12))
"12"
```
"""
dec(n::ZZRingElem) = base(n, 10)

@doc raw"""
    hex(n::ZZRingElem) = base(n, 16)

Return $n$ as a hexadecimal string.

# Examples

```jldoctest
julia> hex(ZZ(12))
"c"
```
"""
hex(n::ZZRingElem) = base(n, 16)

@doc raw"""
    base(n::ZZRingElem, b::Integer)

Return $n$ as a string in base $b$. We require $2 \leq b \leq 62$.

# Examples

```jldoctest
julia> base(ZZ(12), 13)
"c"
```
"""
function base(n::ZZRingElem, b::Integer)
  2 <= b <= 62 || error("invalid base: $b")
  p = ccall((:fmpz_get_str,libflint), Ptr{UInt8},
            (Ptr{UInt8}, Cint, Ref{ZZRingElem}), C_NULL, b, n)
  s = unsafe_string(p)
  ccall((:flint_free, libflint), Nothing, (Ptr{UInt8},), p)
  return s
end

@doc raw"""
    number_of_digits(x::ZZRingElem, b::Integer)

Return the number of digits of $x$ in the base $b$ (default is $b = 10$).

# Examples

```jldoctest
julia> number_of_digits(ZZ(12), 3)
3
```
"""
function number_of_digits(x::ZZRingElem, b::Integer)
  number_of_digits(x, base=b)::Int
end

function number_of_digits(a::ZZRingElem; base::Integer = 10, pad::Integer = 1)
  iszero(a) && return max(pad, 1)
  return max(pad, 1+flog(abs(a), ZZRingElem(abs(base))))
end

Base.digits(n::ZZRingElem; base::Integer = 10, pad::Integer = 1) =
digits(typeof(base), n, base = base, pad = pad)

function Base.digits(T::Type{<:Integer}, n::ZZRingElem; base::Integer = 10, pad::Integer = 1)
  digits!(zeros(T, ndigits(n, base=base, pad=pad)), n, base=base)
end

function Base.digits!(a::AbstractVector{T}, n::ZZRingElem; base::Integer = 10) where T<:Integer
  2 <= base || throw(DomainError(base, "base must be ≥ 2"))
  Base.hastypemax(T) && abs(base) - 1 > typemax(T) &&
  throw(ArgumentError("type $T too small for base $base"))
  isempty(a) && return a

  if nbits(n)/ndigits(base, base = 2) > 100
    c = div(div(nbits(n), ndigits(base, base = 2)), 2)
    nn = ZZRingElem(base)^c
    q, r = divrem(n, nn)

    digits!(view(a, 1:c), r, base = base)
    digits!(view(a, c+1:length(a)), q, base = base)
    return a
  end

  for i in eachindex(a)
    n, r = Base.divrem(n, base)
    a[i] = r
  end
  return a
end

@doc raw"""
    nbits(x::ZZRingElem)

Return the number of binary bits of $x$. We return zero if $x = 0$.

# Examples

```jldoctest
julia> nbits(ZZ(12))
4
```
"""
nbits(x::ZZRingElem) = iszero(x) ? 0 : Int(ccall((:fmpz_bits, libflint), Clong,
                                                 (Ref{ZZRingElem},), x))  

@doc raw"""
    nbits(a::Integer) -> Int

Return the number of bits necessary to represent $a$.
"""
function nbits(a::Integer)
  return ndigits(a, base=2)
end

###############################################################################
#
#   Bit fiddling
#
###############################################################################

@doc raw"""
    popcount(x::ZZRingElem)

Return the number of ones in the binary representation of $x$.

# Examples

```jldoctest
julia> popcount(ZZ(12))
2
```
"""
popcount(x::ZZRingElem) = Int(ccall((:fmpz_popcnt, libflint), UInt,
                                    (Ref{ZZRingElem},), x))

@doc raw"""
    prevpow2(x::ZZRingElem)

Return the previous power of $2$ up to including $x$.
"""
prevpow2(x::ZZRingElem) = x < 0 ? -prevpow2(-x) :
(x <= 2 ? x : one(ZZ) << (ndigits(x, 2) - 1))

@doc raw"""
    nextpow2(x::ZZRingElem)

Return the next power of $2$ that is at least $x$.

# Examples

```jldoctest
julia> nextpow2(ZZ(12))
16
```
"""
nextpow2(x::ZZRingElem) = x < 0 ? -nextpow2(-x) :
(x <= 2 ? x : one(ZZ) << ndigits(x - 1, 2))

@doc raw"""
    trailing_zeros(x::ZZRingElem)

Return the number of trailing zeros in the binary representation of $x$.
"""
trailing_zeros(x::ZZRingElem) = ccall((:fmpz_val2, libflint), Int,
                                      (Ref{ZZRingElem},), x)

###############################################################################
#
#   Bitwise operations (unsafe)
#
###############################################################################

@doc raw"""
    clrbit!(x::ZZRingElem, c::Int)

Clear bit $c$ of $x$, where the least significant bit is the $0$-th bit. Note
that this function modifies its input in-place.

# Examples

```jldoctest
julia> a = ZZ(12)
12

julia> clrbit!(a, 3)

julia> a
4
```
"""
function clrbit!(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Second argument must be non-negative"))
  ccall((:fmpz_clrbit, libflint), Nothing, (Ref{ZZRingElem}, UInt), x, c)
end

@doc raw"""
    setbit!(x::ZZRingElem, c::Int)

Set bit $c$ of $x$, where the least significant bit is the $0$-th bit. Note
that this function modifies its input in-place.

# Examples

```jldoctest
julia> a = ZZ(12)
12

julia> setbit!(a, 0)

julia> a
13
```
"""
function setbit!(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Second argument must be non-negative"))
  ccall((:fmpz_setbit, libflint), Nothing, (Ref{ZZRingElem}, UInt), x, c)
end

@doc raw"""
    combit!(x::ZZRingElem, c::Int)

Complement bit $c$ of $x$, where the least significant bit is the $0$-th bit.
Note that this function modifies its input in-place.

# Examples

```jldoctest
julia> a = ZZ(12)
12

julia> combit!(a, 2)

julia> a
8
```
"""
function combit!(x::ZZRingElem, c::Int)
  c < 0 && throw(DomainError(c, "Second argument must be non-negative"))
  ccall((:fmpz_combit, libflint), Nothing, (Ref{ZZRingElem}, UInt), x, c)
end

@doc raw"""
    tstbit(x::ZZRingElem, c::Int)

Return bit $i$ of x (numbered from 0) as `true` for 1 or `false` for 0.

# Examples

```jldoctest
julia> a = ZZ(12)
12

julia> tstbit(a, 0)
false

julia> tstbit(a, 2)
true
```
"""
function tstbit(x::ZZRingElem, c::Int)
  return c >= 0 && Bool(ccall((:fmpz_tstbit, libflint), Cint,
                              (Ref{ZZRingElem}, UInt),
                              x, c))
end

###############################################################################
#
#   Unsafe operators
#
###############################################################################

function zero!(z::ZZRingElemOrPtr)
  @ccall libflint.fmpz_zero(z::Ref{ZZRingElem})::Nothing
  return z
end

function one!(z::ZZRingElemOrPtr)
  set!(z, UInt(1))
end

function set!(z::ZZRingElemOrPtr, a::ZZRingElemOrPtr)
  @ccall libflint.fmpz_set(z::Ref{ZZRingElem}, a::Ref{ZZRingElem})::Nothing
  return z
end

function set!(z::ZZRingElemOrPtr, a::Int)
  @ccall libflint.fmpz_set_si(z::Ref{ZZRingElem}, a::Int)::Nothing
  return z
end

function set!(z::ZZRingElemOrPtr, a::UInt)
  @ccall libflint.fmpz_set_ui(z::Ref{ZZRingElem}, a::UInt)::Nothing
  return z
end

function swap!(a::ZZRingElem, b::ZZRingElem)
  ccall((:fmpz_swap, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}),
        a, b)
end

function add!(z::ZZRingElem, x::ZZRingElem, y::ZZRingElem)
  ccall((:fmpz_add, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return z
end

function add!(z::ZZRingElem, x::ZZRingElem, y::Int)
  ccall((:fmpz_add_si, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, y)
  return z
end

function add!(z::ZZRingElem, a::ZZRingElem, b::UInt)
  ccall((:fmpz_add_ui, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, UInt),
        z, a, b)
  return z
end

function add!(a::ZZRingElem, b::ZZRingElem, c::Ptr{Int})
  ccall((:fmpz_add, libflint), Nothing, (Ref{ZZRingElem}, Ref{ZZRingElem}, Ptr{Int}), a, b, c)
  return a
end

add!(z::ZZRingElem, a::ZZRingElem, b::Integer) = add!(z, a, ZZRingElem(b))
add!(z::ZZRingElem, x::Int, y::ZZRingElem) = add!(z, y, x)

function neg!(z::ZZRingElemOrPtr, a::ZZRingElemOrPtr)
  ccall((:fmpz_neg, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}),
        z, a)
  return z
end

neg!(a::ZZRingElemOrPtr) = neg!(a, a)

function sub!(z::ZZRingElemOrPtr, a::ZZRingElemOrPtr, b::ZZRingElemOrPtr)
  ccall((:fmpz_sub, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}),
        z, a, b)
  return z
end

function sub!(z::ZZRingElemOrPtr, a::ZZRingElemOrPtr, b::Int)
  ccall((:fmpz_sub_si, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int),
        z, a, b)
  return z
end

function sub!(z::ZZRingElemOrPtr, a::ZZRingElemOrPtr, b::UInt)
  ccall((:fmpz_sub_ui, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, UInt),
        z, a, b)
  return z
end

function sub!(z::ZZRingElemOrPtr, a::ZZRingElemOrPtr, b::Integer)
  return sub!(z, a, ZZRingElem(b))
end

function sub!(z::ZZRingElemOrPtr, b::Integer, a::ZZRingElemOrPtr)
  sub!(z, a, b)
  return neg!(z, z)
end


function mul!(z::ZZRingElemOrPtr, x::ZZRingElemOrPtr, y::ZZRingElemOrPtr)
  ccall((:fmpz_mul, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return z
end

function mul!(z::ZZRingElem, x::ZZRingElem, y::Int)
  ccall((:fmpz_mul_si, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, y)
  return z
end

function mul!(z::ZZRingElem, a::ZZRingElem, b::UInt)
  ccall((:fmpz_mul_ui, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, UInt),
        z, a, b)
  return z
end

mul!(z::ZZRingElem, a::ZZRingElem, b::Integer) = mul!(z, a, ZZRingElem(b))

mul!(z::ZZRingElem, x::Int, y::ZZRingElem) = mul!(z, y, x)

function mul!(a::ZZRingElem, b::ZZRingElem, c::Ptr{Int})
  ccall((:fmpz_mul, libflint), Nothing, (Ref{ZZRingElem}, Ref{ZZRingElem}, Ptr{Int}), a, b, c)
  return a
end

function addmul!(z::ZZRingElem, x::ZZRingElem, y::ZZRingElem)
  ccall((:fmpz_addmul, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return z
end

addmul!(z::ZZRingElem, x::ZZRingElem, y::ZZRingElem, ::ZZRingElem) = addmul!(z, x, y)

function addmul!(z::ZZRingElem, x::ZZRingElem, y::Int)
  ccall((:fmpz_addmul_si, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), z, x, y)
  return z
end

addmul!(z::ZZRingElem, x::ZZRingElem, y::Int, ::ZZRingElem) = addmul!(z, x, y)

function submul!(z::ZZRingElem, a::ZZRingElem, b::ZZRingElem)
  ccall((:fmpz_submul, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}),
        z, a, b)
  return z
end

@doc raw"""
    fmma!(r::ZZRingElem, a::ZZRingElem, b::ZZRingElem, c::ZZRingElem, d::ZZRingElem)

Return $r = a b + c d$, changing $r$ in-place.
"""
function fmma!(r::ZZRingElem, a::ZZRingElem, b::ZZRingElem, c::ZZRingElem, d::ZZRingElem)
  ccall((:fmpz_fmma, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), r, a, b, c, d)
  return r
end

@doc raw"""
    fmms!(r::ZZRingElem, a::ZZRingElem, b::ZZRingElem, c::ZZRingElem, d::ZZRingElem)

Return $r = a b - c d$, changing $r$ in-place.
"""
function fmms!(r::ZZRingElem, a::ZZRingElem, b::ZZRingElem, c::ZZRingElem, d::ZZRingElem)
  ccall((:fmpz_fmms, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), r, a, b, c, d)
  return r
end


function divexact!(z::ZZRingElem, a::ZZRingElem, b::UInt)
  ccall((:fmpz_divexact_ui, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, UInt),
        z, a, b)
  return z
end

function divexact!(z::ZZRingElem, a::ZZRingElem, b::ZZRingElem)
  ccall((:fmpz_divexact, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}),
        z, a, b)
  return z
end

divexact!(z::ZZRingElem, b::ZZRingElem) = divexact!(z, z, b)

function pow!(z::ZZRingElem, a::ZZRingElem, b::Union{Int, UInt})
  ccall((:fmpz_pow_ui, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, UInt),
        z, a, UInt(b))
  return z
end

function lcm!(z::ZZRingElem, x::ZZRingElem, y::ZZRingElem)
  ccall((:fmpz_lcm, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return z
end

function gcd!(z::ZZRingElem, x::ZZRingElem, y::ZZRingElem)
  ccall((:fmpz_gcd, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{ZZRingElem}, Ref{ZZRingElem}), z, x, y)
  return z
end

###############################################################################
#
#   Parent object overloads
#
###############################################################################

(::ZZRing)() = ZZRingElem()

(::ZZRing)(a::Integer) = ZZRingElem(a)

(::ZZRing)(a::AbstractString) = ZZRingElem(a)

(::ZZRing)(a::ZZRingElem) = a

(::ZZRing)(a::Float64) = ZZRingElem(a)

(::ZZRing)(a::Float32) = ZZRingElem(Float64(a))

(::ZZRing)(a::Float16) = ZZRingElem(Float64(a))

(::ZZRing)(a::BigFloat) = ZZRingElem(BigInt(a))

@doc raw"""
    (ZZ::ZZRing)(x)

Coerce `x` into an element of $\mathbb Z$. Note that `ZZ(x)` is equivalent to [`ZZRingElem(x)`](@ref).

# Examples

```jldoctest
julia> ZZ(2)
2

julia> ZZ(2)^100
1267650600228229401496703205376
```
"""
(::ZZRing)(x)

###############################################################################
#
#   String parser
#
###############################################################################

function parse(::Type{ZZRingElem}, s::AbstractString, base::Int = 10)
  if s[1] == '-'
    sgn = -1
    s2 = string(SubString(s, 2))
  else
    sgn = 1
    s2 = string(s)
  end
  z = ZZRingElem()
  err = ccall((:fmpz_set_str, libflint),
              Int32, (Ref{ZZRingElem}, Ptr{UInt8}, Int32),
              z, s2, base)
  err == 0 || error("Invalid big integer: $(repr(s))")
  return sgn < 0 ? -z : z
end

###############################################################################
#
#   Random generation
#
###############################################################################

RandomExtensions.maketype(R::ZZRing, _) = ZZRingElem

# define rand(make(ZZ, n:m))
rand(rng::AbstractRNG, sp::SamplerTrivial{<:Make2{ZZRingElem,ZZRing}}) =
sp[][1](rand(rng, sp[][2]))

rand(rng::AbstractRNG, R::ZZRing, n::AbstractArray) = R(rand(rng, n))

rand(R::ZZRing, n::AbstractArray) = rand(Random.GLOBAL_RNG, R, n)

@doc raw"""
    rand_bits(::ZZRing, b::Int)

Return a random signed integer whose absolute value has $b$ bits.
"""
function rand_bits(::ZZRing, b::Int)
  b >= 0 || throw(DomainError(b, "Bit count must be non-negative"))
  z = ZZRingElem()
  ccall((:fmpz_randbits, libflint), Nothing,(Ref{ZZRingElem}, Ref{rand_ctx}, Int),
        z, _flint_rand_states[Threads.threadid()], b)
  return z
end

@doc raw"""
    rand_bits_prime(::ZZRing, n::Int, proved::Bool=true)

Return a random prime number with the given number of bits. If only a
probable prime is required, one can pass `proved=false`.
"""
function rand_bits_prime(::ZZRing, n::Int, proved::Bool = true)
  n < 2 && throw(DomainError(n, "No primes with that many bits"))
  z = ZZRingElem()
  ccall((:fmpz_randprime, libflint), Nothing,
        (Ref{ZZRingElem}, Ref{rand_ctx}, Int, Cint),
        z, _flint_rand_states[Threads.threadid()], n, Cint(proved))
  return z
end

# rand in a range
# this mirrors the implementation for BigInt in the Random module

using Base.GMP: Limb, MPZ

# "views" a non-small ZZRingElem as a readonly BigInt
# the caller must call GC.@preserve appropriately
function _as_bigint(z::ZZRingElem)
  @assert !_fmpz_is_small(z)
  unsafe_load(Ptr{BigInt}(z.d << 2))
end


Random.Sampler(::Type{<:AbstractRNG}, r::StepRange{ZZRingElem, ZZRingElem}, ::Random.Repetition) =
SamplerFmpz(r)

struct SamplerFmpz <: Random.Sampler{ZZRingElem}
  a::ZZRingElem           # first
  m::BigInt         # range length - 1
  nlimbs::Int       # number of limbs in generated BigInt's (z ∈ [0, m])
  nlimbsmax::Int    # max number of limbs for z+a
  mask::Limb        # applied to the highest limb

  ## diverges from Random.SamplerBigInt:
  step::ZZRingElem
end

function SamplerFmpz(r::StepRange{ZZRingElem, ZZRingElem})
  r1 = first(r)
  r2 = last(r)
  s = step(r)

  if isone(s)
    m = BigInt(r2 - r1)
  else
    m = length(r)::BigInt # type assertion in case length is changed to return an ZZRingElem
    MPZ.sub_ui!(m, 1)
  end

  m < 0 && throw(ArgumentError("range must be non-empty"))
  nd = ndigits(m, base=2)
  nlimbs, highbits = divrem(nd, 8*sizeof(Limb))
  highbits > 0 && (nlimbs += 1)
  mask = highbits == 0 ? ~zero(Limb) : one(Limb)<<highbits - one(Limb)
  GC.@preserve r1 r2 begin
    a1 = _fmpz_is_small(r1) ? 1 : _as_bigint(r1).size
    a2 = _fmpz_is_small(r2) ? 1 : _as_bigint(r2).size
  end
  nlimbsmax = max(nlimbs, abs(a1), abs(a2))
  return SamplerFmpz(r1, m, nlimbs, nlimbsmax, mask, s)
end

function rand(rng::AbstractRNG, sp::SamplerFmpz)
  z = ZZRingElem()
  # this make sure z is backed up by an mpz_t object:
  ccall((:fmpz_init2, libflint), Nothing, (Ref{ZZRingElem}, UInt), z, sp.nlimbsmax)
  @assert !_fmpz_is_small(z)
  GC.@preserve z begin
    x = _as_bigint(z)
    limbs = Random.UnsafeView(x.d, sp.nlimbs)
    while true
      rand!(rng, limbs)
      limbs[end] &= sp.mask
      MPZ.mpn_cmp(x, sp.m, sp.nlimbs) <= 0 && break
    end
    # adjust x.size (normally done by mpz_limbs_finish, in GMP version >= 6)
    sz = sp.nlimbs
    while sz > 0
      limbs[sz] != 0 && break
      sz -= 1
    end
    # write sz in the .size field of the mpz object
    unsafe_store!(Ptr{Cint}(z.d << 2) + sizeof(Cint), sz)
  end
  if !isone(sp.step)
    mul!(z, z, sp.step)
  end
  add!(z, z, sp.a)
end

#TODO
# need to be mapped onto proper Flint primitives
# flints needs a proper interface to randomness - I think
# currently one simply cannot use it at all
#
# should be tied(?) to the Julia rng stuff?
# similarly, all the derived rand functions should probably also do this
#
# inspired by/copied from a former BigInt implementation from the stdlib in
# `Random/src/generation.jl`
#

function rand(rng::AbstractRNG, a::ZZRingElemUnitRange)
  m = Base.last(a) - Base.first(a)
  m < 0 && error("range empty")
  nd = ndigits(m, 2)
  nl, high = Base.divrem(nd, 8 * sizeof(Base.GMP.Limb))
  if high > 0
    mask = m >> (nl * 8 * sizeof(Base.GMP.Limb))
  end
  s = ZZRingElem(0)
  c = (8 * sizeof(Base.GMP.Limb))
  while true
    s = ZZRingElem(0)
    for i = 1:nl
      ccall((:fmpz_mul_2exp, libflint), Nothing,
            (Ref{ZZRingElem}, Ref{ZZRingElem}, Int), s, s, c)
      add!(s, s, rand(rng, Base.GMP.Limb))
    end
    if high > 0
      s = s << high
      s += rand(rng, 0:Base.GMP.Limb(mask))
    end
    if s <= m
      break
    end
  end
  return s + first(a)
end

struct RangeGeneratorfmpz# <: Base.Random.RangeGenerator
  a::StepRange{ZZRingElem,ZZRingElem}
end

function Random.RangeGenerator(r::StepRange{ZZRingElem,ZZRingElem})
  m = last(r) - first(r)
  m < 0 && throw(ArgumentError("range must be non-empty"))
  return RangeGeneratorfmpz(r)
end

function rand(rng::AbstractRNG, g::RangeGeneratorfmpz)
  return rand(rng, g.a)
end

function rand!(A::Vector{ZZRingElem}, v::StepRange{ZZRingElem,ZZRingElem})
  for i in 1:length(A)
    A[i] = rand(v)
  end
  return A
end

###############################################################################
#
#   Constructors
#
###############################################################################

ZZRingElem(s::AbstractString) = parse(ZZRingElem, s)

ZZRingElem(z::Integer) = ZZRingElem(BigInt(z))

ZZRingElem(z::Float16) = ZZRingElem(Float64(z))

ZZRingElem(z::Float32) = ZZRingElem(Float64(z))

ZZRingElem(z::BigFloat) = ZZRingElem(BigInt(z))

###############################################################################
#
#   Conversions and promotions
#
###############################################################################

convert(::Type{ZZRingElem}, a::Integer) = ZZRingElem(a)

function (::Type{BigInt})(a::ZZRingElem)
  r = BigInt()
  ccall((:fmpz_get_mpz, libflint), Nothing, (Ref{BigInt}, Ref{ZZRingElem}), r, a)
  return r
end

function (::Type{Int})(a::ZZRingElem)
  (a > typemax(Int) || a < typemin(Int)) && throw(InexactError(:convert, Int, a))
  return ccall((:fmpz_get_si, libflint), Int, (Ref{ZZRingElem},), a)
end

function (::Type{UInt})(a::ZZRingElem)
  (a > typemax(UInt) || a < 0) && throw(InexactError(:convert, UInt, a))
  return ccall((:fmpz_get_ui, libflint), UInt, (Ref{ZZRingElem}, ), a)
end

if Culong !== UInt
  function (::Type{Culong})(a::ZZRingElem)
    fits(Culong, a) || throw(InexactError(:convert, Culong, a))
    return ccall((:fmpz_get_ui, libflint), UInt, (Ref{ZZRingElem}, ), a)%Culong
  end
end

(::Type{T})(a::ZZRingElem) where T <: Union{Int8, Int16, Int32} = T(Int(a))

(::Type{T})(a::ZZRingElem) where T <: Union{UInt8, UInt16, UInt32} = T(UInt(a))

(::Type{Int128})(a::ZZRingElem) = Int128(BigInt(a))

Integer(a::ZZRingElem) = BigInt(a)

convert(::Type{T}, a::ZZRingElem) where T <: Integer = T(a)

function (::Type{Float64})(n::ZZRingElem)
  # rounds to zero
  ccall((:fmpz_get_d, libflint), Float64, (Ref{ZZRingElem},), n)
end

convert(::Type{Float64}, n::ZZRingElem) = Float64(n)

(::Type{Float32})(n::ZZRingElem) = Float32(Float64(n))

convert(::Type{Float32}, n::ZZRingElem) = Float32(n)

(::Type{Float16})(n::ZZRingElem) = Float16(Float64(n))

convert(::Type{Float16}, n::ZZRingElem) = Float16(n)

(::Type{BigFloat})(n::ZZRingElem) = BigFloat(BigInt(n))

convert(::Type{BigFloat}, n::ZZRingElem) = BigFloat(n)

Base.promote_rule(::Type{ZZRingElem}, ::Type{T}) where {T <: Integer} = ZZRingElem

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

###############################################################################
#
#  Perfect power detection
#
###############################################################################

# 1, 0, -1 are perfect powers
# ex is not guaranteed to be maximal
function _is_perfect_power(a::ZZRingElem)
  rt = ZZRingElem()
  ex = ccall((:fmpz_is_perfect_power, libflint), Int, (Ref{ZZRingElem}, Ref{ZZRingElem}), rt, a)
  return rt, ex
end

@doc raw"""
    is_perfect_power(a::IntegerUnion)

Return whether $a$ is a perfect power, that is, whether $a = m^r$ for some
integer $m$ and $r > 1$.
"""
function is_perfect_power(a::ZZRingElem)
  _, ex = _is_perfect_power(a)
  return ex > 0
end

is_perfect_power(a::Integer) = is_perfect_power(ZZRingElem(a))

#compare to Oscar/examples/PerfectPowers.jl which is, for large input,
#far superior over gmp/ fmpz_is_perfect_power

@doc raw"""
    is_perfect_power_with_data(a::ZZRingElem) -> Int, ZZRingElem
    is_perfect_power_with_data(a::Integer) -> Int, Integer

Return $e$, $r$ such that $a = r^e$ with $e$ maximal. Note: $1 = 1^0$.
"""
function is_perfect_power_with_data(a::ZZRingElem)
  if iszero(a)
    error("must not be zero")
  end
  if isone(a)
    return 0, a
  end
  if a < 0
    e, r = is_perfect_power_with_data(-a)
    if isone(e)
      return 1, a
    end
    v, s = iszero(e) ? (0, 0) : remove(e, 2)
    return s, -r^(2^v)
  end
  rt = ZZRingElem()
  e = 1
  while true
    rt, ex = _is_perfect_power(a)
    if ex == 1 || ex == 0
      return e, a
    end
    e *= ex
    a = rt
  end
end

function is_perfect_power_with_data(a::Integer)
  e, r = is_perfect_power_with_data(ZZRingElem(a))
  return e, typeof(a)(r)
end

@doc raw"""
    is_power(a::ZZRingElem, n::Int) -> Bool, ZZRingElem
    is_power(a::QQFieldElem, n::Int) -> Bool, QQFieldElem
    is_power(a::Integer, n::Int) -> Bool, Integer

Return `true` and the root if $a$ is an $n$-th power.
"""
function is_power(a::ZZRingElem, n::Int)
  if a < 0 && is_even(n)
    return false, a
  end
  b = iroot(a, n)
  return b^n == a, b
end

@doc raw"""
    is_prime_power(q::IntegerUnion) -> Bool

Returns whether $q$ is a prime power.
"""
is_prime_power(::IntegerUnion)

function is_prime_power(q::ZZRingElem)
  iszero(q) && return false
  e, a = is_perfect_power_with_data(q)
  return is_prime(a)
end

is_prime_power(q::Integer) = is_prime_power(ZZRingElem(q))

@doc raw"""
    is_prime_power_with_data(q::IntegerUnion) -> Bool, ZZRingElem, Int

Returns a flag indicating whether $q$ is a prime power and integers $e, p$ such
that $q = p^e$. If $q$ is a prime power, than $p$ is a prime.
"""
is_prime_power_with_data(::IntegerUnion)

function is_prime_power_with_data(q::ZZRingElem)
  iszero(q) && return false, 1, q
  e, a = is_perfect_power_with_data(q)
  return is_prime(a), e, a
end

function is_prime_power_with_data(q::Integer)
  e, a = is_perfect_power_with_data(ZZRingElem(q))
  return is_prime(a), e, typeof(q)(a)
end

###############################################################################
#
#   Convenience methods for arithmetics (since `ZZRingElem` is not a `Number` type)
#
###############################################################################

//(v::Vector{ZZRingElem}, x::ZZRingElem) = v .// x
*(x::ZZRingElem, v::Vector{ZZRingElem}) = x .* v
*(v::Vector{ZZRingElem}, x::ZZRingElem) = v .* x

###############################################################################
#
#   Bit iteration
#
###############################################################################

function bits end

module BitsMod

using ..Nemo

import Base: ^
import Base: eltype
import Base: getindex
import Base: iterate
import Base: length
import Base: show

export bits
export Limbs

const hb = UInt(1) << 63

struct Limbs
  a::ZZRingElem
  len::Int
  b::Ptr{UInt}
  function Limbs(a::ZZRingElem; MSW::Bool=true)
    if Nemo._fmpz_is_small(a)
      return new(a, 0, convert(Ptr{UInt}, 0))
    end
    z = convert(Ptr{Cint}, unsigned(a.d) << 2)
    len = unsafe_load(z, 2)
    d = convert(Ptr{Ptr{UInt}}, unsigned(a.d) << 2) + 2 * sizeof(Cint)
    p = unsafe_load(d)
    if !MSW
      new(a, -len, p)
    else
      new(a, len, p)
    end
  end
end

function show(io::IO, L::Limbs)
  print(io, "limb-access for: ", L.a)
end

@inline function getindex(L::Limbs, i::Int)
  if L.len == 0
    return UInt(abs(L.a.d)) #error???
  end
  @boundscheck @assert i <= abs(L.len)
  return unsafe_load(L.b, i)
end

function iterate(L::Limbs)
  is_zero(L.a) && return nothing # nbits(ZZ()) == 0
  L.len < 0 && return L[1], 1

  return L[L.len], L.len
end

function iterate(L::Limbs, i::Int)
  if L.len < 0
    i > -L.len && return nothing
    return L[i+1], i + 1
  end
  i == 0 && return nothing
  return L[i-1], i - 1
end

function length(L::Limbs)
  is_zero(L.a) && return 0 # nbits(ZZ()) == 0
  return L.len + 1
end

eltype(L::Limbs) = UInt

#=
#from https://github.com/JuliaLang/julia/issues/11592
#compiles directly down to the ror/rol in assembly
for T in Base.BitInteger_types
mask = UInt8(sizeof(T) << 3 - 1)
@eval begin
ror(x::$T, k::Integer) = (x >>> ($mask & k)) | (x <<  ($mask & -k))
rol(x::$T, k::Integer) = (x <<  ($mask & k)) | (x >>> ($mask & -k))
end
end
=#

struct BitsFmpz
  L::Limbs

  function BitsFmpz(b::ZZRingElem)
    return new(Limbs(b))
  end
end

function iterate(B::BitsFmpz)
  L = B.L
  is_zero(L.a) && return nothing
  a = L[L.len]
  b = UInt(1) << (nbits(a) - 1)
  return true, (b, L.len)
end

@inline function iterate(B::BitsFmpz, s::Tuple{UInt,Int})
  b = s[1] >> 1
  if b == 0
    l = s[2] - 1
    if l < 1
      return nothing
    end
    b = hb
    a = B.L[l]
    return a & b != 0, (b, l)
  end
  return B.L[s[2]] & b != 0, (b, s[2])
end

function show(io::IO, B::BitsFmpz)
  print(io, "bit iterator for: ", B.L.a)
end

length(B::BitsFmpz) = nbits(B.L.a)

eltype(B::BitsFmpz) = Bool

Nemo.bits(a::ZZRingElem) = BitsFmpz(a)
#= wrong order, thus disabled

function getindex(B::BitsFmpz, i::Int)
return ccall((:fmpz_tstbit, libflint), Int, (Ref{ZZRingElem}, Int), B.L.a, i) != 0
end
=#
end

using .BitsMod