text:
const pyjlarraytype = pynew()

function pyjl_getaxisindex(x::AbstractUnitRange{<:Integer}, k::Py)
    if pyisslice(k)
        a = @pyconvert_and_del Union{Int,Nothing} k.start begin
            errset(pybuiltins.TypeError, "slice components must be integers")
            pythrow()
        end
        b = @pyconvert_and_del Union{Int,Nothing} k.step begin
            errset(pybuiltins.TypeError, "slice components must be integers")
            pythrow()
        end
        c = @pyconvert_and_del Union{Int,Nothing} k.stop begin
            errset(pybuiltins.TypeError, "slice components must be integers")
            pythrow()
        end
        # step defaults to 1
        b′ = b === nothing ? 1 : b
        if a === nothing && c === nothing
            # when neither is specified, start and stop default to the full range,
            # which is reversed when the step is negative
            if b′ > 0
                a′ = Int(first(x))
                c′ = Int(last(x))
            elseif b′ < 0
                a′ = Int(last(x))
                c′ = Int(first(x))
            else
                errset(pybuiltins.ValueError, "step must be non-zero")
                pythrow()
            end
        else
            # start defaults
            a′ = Int(a === nothing ? first(x) : a < 0 ? (last(x) + a + 1) : (first(x) + a))
            c′ = Int(
                c === nothing ? last(x) :
                c < 0 ? (last(x) + 1 + c - sign(b′)) : (first(x) + c - sign(b′)),
            )
        end
        r = StepRange{Int,Int}(a′, b′, c′)
        if checkbounds(Bool, x, r)
            return r
        else
            errset(pybuiltins.IndexError, "array index out of bounds")
            pythrow()
        end
    else
        j = @pyconvert Int k begin
            errset(pybuiltins.TypeError, "index must be slice or integer, got '$(pytype(k).__name__)'")
            pythrow()
        end
        r = Int(j < 0 ? (last(x) + j + 1) : (first(x) + j))
        if checkbounds(Bool, x, r)
            return r
        else
            errset(pybuiltins.IndexError, "array index out of bounds")
            pythrow()
        end
    end
end

function pyjl_getarrayindices(x::AbstractArray{T,N}, ks::Py) where {T,N}
    if pyistuple(ks)
        if pylen(ks) == N
            return ntuple(N) do i
                k = pytuple_getitem(ks, i-1)
                ans = pyjl_getaxisindex(axes(x, i), k)
                pydel!(k)
                return ans
            end
        else
            errset(pybuiltins.TypeError, "expecting $N indices, got $(pylen(ks))")
            pythrow()
        end
    elseif N == 1
        return (pyjl_getaxisindex(axes(x, 1), ks),)
    else
        errset(pybuiltins.TypeError, "expecting $N indices, got 1")
    end
end

function pyjlarray_getitem(x::AbstractArray{T,N}, k_::Py) where {T,N}
    k = pyjl_getarrayindices(x, k_)
    if k isa NTuple{N,Int}
        return Py(x[k...])
    else
        return Py(view(x, k...))
    end
end

function pyjlarray_setitem(x::AbstractArray{T,N}, k_::Py, v_::Py) where {T,N}
    k = pyjl_getarrayindices(x, k_)
    if k isa NTuple{N,Int}
        v = pyconvertarg(T, v_, "value")
        x[k...] = v
    else
        v = pyconvertarg(Any, v_, "value")
        x[k...] .= v
    end
    return Py(nothing)
end

function pyjlarray_delitem(x::AbstractArray{T,N}, k_::Py) where {T,N}
    if N == 1
        k = pyjl_getarrayindices(x, k_)
        deleteat!(x, k...)
    else
        errset(pybuiltins.TypeError, "can only delete from 1D arrays")
        pythrow()
    end
    return Py(nothing)
end
pyjl_handle_error_type(::typeof(pyjlarray_delitem), x, exc::MethodError) = exc.f === deleteat! ? pybuiltins.TypeError : PyNULL

function pyjlarray_reshape(x::AbstractArray, shape_::Py)
    shape = pyconvertarg(Union{Int,Vector{Int}}, shape_, "shape")
    return Py(reshape(x, shape...))
end

pyjlarray_isbufferabletype(::Type{T}) where {T} = T in (
    Int8,
    Int16,
    Int32,
    Int64,
    UInt8,
    UInt16,
    UInt32,
    UInt64,
    Float16,
    Float32,
    Float64,
    Complex{Float16},
    Complex{Float32},
    Complex{Float64},
    Bool,
    Ptr{Cvoid},
)
pyjlarray_isbufferabletype(::Type{T}) where {T<:Tuple} =
    isconcretetype(T) &&
    PythonCall.allocatedinline(T) &&
    all(pyjlarray_isbufferabletype, fieldtypes(T))
pyjlarray_isbufferabletype(::Type{NamedTuple{names,T}}) where {names,T} =
    pyjlarray_isbufferabletype(T)

function pyjlarray_buffer_info(x::AbstractArray{T,N}) where {T,N}
    if pyjlarray_isbufferabletype(T)
        C.PyBufferInfo{N}(
            ptr = Base.unsafe_convert(Ptr{T}, x),
            readonly = !Utils.ismutablearray(x),
            itemsize = sizeof(T),
            format = pybufferformat(T),
            shape = size(x),
            strides = strides(x) .* Base.aligned_sizeof(T),
        )
    else
        error("element type is not bufferable")
    end
end

const PYBUFFERFORMAT = IdDict{Type,String}()

pybufferformat(::Type{T}) where {T} = get!(PYBUFFERFORMAT, T) do
    T == Cchar ? "b" :
    T == Cuchar ? "B" :
    T == Cshort ? "h" :
    T == Cushort ? "H" :
    T == Cint ? "i" :
    T == Cuint ? "I" :
    T == Clong ? "l" :
    T == Culong ? "L" :
    T == Clonglong ? "q" :
    T == Culonglong ? "Q" :
    T == Float16 ? "e" :
    T == Cfloat ? "f" :
    T == Cdouble ? "d" :
    T == Complex{Float16} ? "Ze" :
    T == Complex{Cfloat} ? "Zf" :
    T == Complex{Cdouble} ? "Zd" :
    T == Bool ? "?" :
    T == Ptr{Cvoid} ? "P" :
    if isstructtype(T) && isconcretetype(T) && allocatedinline(T)
        n = fieldcount(T)
        flds = []
        for i = 1:n
            nm = fieldname(T, i)
            tp = fieldtype(T, i)
            push!(flds, string(pybufferformat(tp), nm isa Symbol ? ":$nm:" : ""))
            d =
                (i == n ? sizeof(T) : fieldoffset(T, i + 1)) -
                (fieldoffset(T, i) + sizeof(tp))
            @assert d ≥ 0
            d > 0 && push!(flds, "$(d)x")
        end
        string("T{", join(flds, " "), "}")
    else
        "$(sizeof(T))x"
    end
end

pyjlarray_isarrayabletype(::Type{T}) where {T} = T in (
    UInt8,
    Int8,
    UInt16,
    Int16,
    UInt32,
    Int32,
    UInt64,
    Int64,
    Bool,
    Float16,
    Float32,
    Float64,
    Complex{Float16},
    Complex{Float32},
    Complex{Float64},
)
pyjlarray_isarrayabletype(::Type{T}) where {T<:Tuple} =
    isconcretetype(T) &&
    Base.allocatedinline(T) &&
    all(pyjlarray_isarrayabletype, T.parameters)
pyjlarray_isarrayabletype(::Type{NamedTuple{names,types}}) where {names,types} =
    pyjlarray_isarrayabletype(types)

const PYTYPESTRDESCR = IdDict{Type,Tuple{String,Py}}()

pytypestrdescr(::Type{T}) where {T} = get!(PYTYPESTRDESCR, T) do
    c = Utils.islittleendian() ? '<' : '>'
    T == Bool ? ("$(c)b$(sizeof(Bool))", PyNULL) :
    T == Int8 ? ("$(c)i1", PyNULL) :
    T == UInt8 ? ("$(c)u1", PyNULL) :
    T == Int16 ? ("$(c)i2", PyNULL) :
    T == UInt16 ? ("$(c)u2", PyNULL) :
    T == Int32 ? ("$(c)i4", PyNULL) :
    T == UInt32 ? ("$(c)u4", PyNULL) :
    T == Int64 ? ("$(c)i8", PyNULL) :
    T == UInt64 ? ("$(c)u8", PyNULL) :
    T == Float16 ? ("$(c)f2", PyNULL) :
    T == Float32 ? ("$(c)f4", PyNULL) :
    T == Float64 ? ("$(c)f8", PyNULL) :
    T == Complex{Float16} ? ("$(c)c4", PyNULL) :
    T == Complex{Float32} ? ("$(c)c8", PyNULL) :
    T == Complex{Float64} ? ("$(c)c16", PyNULL) :
    if isstructtype(T) && isconcretetype(T) && Base.allocatedinline(T)
        n = fieldcount(T)
        flds = []
        for i = 1:n
            nm = fieldname(T, i)
            tp = fieldtype(T, i)
            ts, ds = pytypestrdescr(tp)
            isempty(ts) && return ("", PyNULL)
            push!(
                flds,
                (nm isa Integer ? "f$(nm-1)" : string(nm), ds === nothing ? ts : ds),
            )
            d = (i == n ? sizeof(T) : fieldoffset(T, i + 1)) - (fieldoffset(T, i) + sizeof(tp))
            @assert d ≥ 0
            d > 0 && push!(flds, ("", "|V$(d)"))
        end
        ("|$(sizeof(T))V", pylist(flds))
    else
        ("", PyNULL)
    end
end

pyjlarray_array__array(x::AbstractArray) = x isa Array ? Py(nothing) : pyjl(Array(x))
pyjlarray_array__pyobjectarray(x::AbstractArray) = pyjl(PyObjectArray(x))

function pyjlarray_array_interface(x::AbstractArray{T,N}) where {T,N}
    if pyjlarray_isarrayabletype(eltype(x))
        data = (UInt(Base.unsafe_convert(Ptr{T}, x)), !Utils.ismutablearray(x))
        typestr, descr = pytypestrdescr(T)
        if !isempty(typestr)
            d = pydict()
            d["shape"] = size(x)
            d["typestr"] = typestr
            d["data"] = data
            d["strides"] = strides(x) .* Base.aligned_sizeof(T)
            d["version"] = 3
            if !ispynull(descr)
                d["descr"] = descr
            end
            return d
        end
    end
    errset(pybuiltins.AttributeError, "__array_interface__")
    return pynew()
end
pyjl_handle_error_type(::typeof(pyjlarray_array_interface), x, exc) = pybuiltins.AttributeError

function init_jlwrap_array()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class ArrayValue(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        _jl_buffer_info = $(pyjl_methodnum(pyjlarray_buffer_info))
        @property
        def ndim(self):
            return self._jl_callmethod($(pyjl_methodnum(Py ∘ ndims)))
        @property
        def shape(self):
            return self._jl_callmethod($(pyjl_methodnum(Py ∘ size)))
        def copy(self):
            return self._jl_callmethod($(pyjl_methodnum(Py ∘ copy)))
        def reshape(self, shape):
            return self._jl_callmethod($(pyjl_methodnum(pyjlarray_reshape)), shape)
        def __getitem__(self, k):
            return self._jl_callmethod($(pyjl_methodnum(pyjlarray_getitem)), k)
        def __setitem__(self, k, v):
            self._jl_callmethod($(pyjl_methodnum(pyjlarray_setitem)), k, v)
        def __delitem__(self, k):
            self._jl_callmethod($(pyjl_methodnum(pyjlarray_delitem)), k)
        @property
        def __array_interface__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlarray_array_interface)))
        def __array__(self):
            # convert to an array-like object
            arr = self
            if not (hasattr(arr, "__array_interface__") or hasattr(arr, "__array_struct__")):
                # the first attempt collects into an Array
                arr = self._jl_callmethod($(pyjl_methodnum(pyjlarray_array__array)))
                if not (hasattr(arr, "__array_interface__") or hasattr(arr, "__array_struct__")):
                    # the second attempt collects into a PyObjectArray
                    arr = self._jl_callmethod($(pyjl_methodnum(pyjlarray_array__pyobjectarray)))
            # convert to a numpy array if numpy is available
            try:
                import numpy
                arr = numpy.array(arr)
            except ImportError:
                pass
            return arr
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlarraytype, jl.ArrayValue)
end

pyjl(v::AbstractArray) = pyjl(pyjlarraytype, v)


text:
# run all serial tests in a single session

args_orig = copy(ARGS)

cd("./advection")
include(joinpath(pwd(), "runtests.jl"))

cd("../euler")
resize!(ARGS, length(args_orig))  # reset to default
copy!(ARGS, args_orig)
include(joinpath(pwd(), "runtests.jl"))

cd("../navier_stokes")
resize!(ARGS, length(args_orig))  # reset to default
copy!(ARGS, args_orig)
include(joinpath(pwd(), "runtests.jl"))

cd("../simpleODE")
resize!(ARGS, length(args_orig))  # reset to default
copy!(ARGS, args_orig)
include(joinpath(pwd(), "runtests.jl"))

cd("../elliptic")
resize!(ARGS, length(args_orig))
copy!(ARGS, args_orig)
include(joinpath(pwd(), "runtests.jl"))

cd("../")

ext: jl


text:
################################################################################
#
#  NfRel/NfRel.jl : Relative number field extensions
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2017 Claus Fieker, Tommy Hofmann
#
################################################################################

export absolute_field

add_assert_scope(:NfRel)
################################################################################
#
#  Copy
#
################################################################################

function Base.deepcopy_internal(a::NfRelElem{T}, dict::IdDict) where T
  z = NfRelElem{T}(Base.deepcopy_internal(data(a), dict))
  z.parent = parent(a)
  return z
end

################################################################################
#
#  Comply with ring interface
#
################################################################################

needs_parentheses(::NfRelElem) = true

isnegative(x::NfRelElem) = isnegative(data(x))

show_minus_one(::Type{NfRelElem{T}}) where {T} = true

function iszero(a::NfRelElem)
  reduce!(a)
  return iszero(data(a))
end

function isone(a::NfRelElem)
  reduce!(a)
  return isone(data(a))
end

zero(K::NfRel) = K(zero(parent(K.pol)))

one(K::NfRel) = K(one(parent(K.pol)))

function zero!(a::NfRelElem)
  zero!(a.data)
  return a
end

################################################################################
#
#  Promotion
#
################################################################################

AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{fmpz}) where {T <: NumFieldElem} = NfRelElem{T}

AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{fmpq}) where {T <: NumFieldElem} = NfRelElem{T}

AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{T}) where {T} = NfRelElem{T}

AbstractAlgebra.promote_rule(::Type{T}, ::Type{NfRelElem{T}}) where {T} = NfRelElem{T}

AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{NfRelElem{T}}) where T <: NumFieldElem = NfRelElem{T}

function AbstractAlgebra.promote_rule(::Type{NfRelElem{T}}, ::Type{U}) where {T <: NumFieldElem, U <: NumFieldElem}
  AbstractAlgebra.promote_rule(T, U) == T ? NfRelElem{T} : Union{}
end

function AbstractAlgebra.promote_rule(::Type{U}, ::Type{NfRelElem{T}}) where {T <: NumFieldElem, U <: NumFieldElem}
  AbstractAlgebra.promote_rule(T, U) == T ? NfRelElem{T} : Union{}
end

################################################################################
#
#  Order types
#
################################################################################

order_type(K::NfRel{T}) where {T} = NfRelOrd{T, fractional_ideal_type(order_type(base_field(K))), NfRelElem{T}}

order_type(::Type{NfRel{T}}) where {T} = NfRelOrd{T, fractional_ideal_type(order_type(parent_type(T))), NfRelElem{T}}

################################################################################
#
#  Field access
#
################################################################################

@inline base_field(a::NfRel{T}) where {T} = a.base_ring::parent_type(T)

@inline data(a::NfRelElem) = a.data

@inline parent(a::NfRelElem{T}) where {T} = a.parent::NfRel{T}

@inline issimple(a::NfRel) = true

################################################################################
#
#  Coefficient setting and getting
#
################################################################################

@inline coeff(a::NfRelElem{T}, i::Int) where {T} = coeff(a.data, i)

@inline setcoeff!(a::NfRelElem{T}, i::Int, c::T) where {T} = setcoeff!(a.data, i, c)

# copy does not do anything (so far), this is only for compatibility with coeffs(::AbsAlgAssElem)
function coeffs(a::NfRelElem{T}; copy::Bool = false) where {T}
  return T[coeff(a, i) for i = 0:degree(parent(a)) - 1 ]
end

################################################################################
#
#  Constructors
#
################################################################################
#TODO: I don't understand those
(K::NfRel)(x::NfRelElem) = K(base_field(K)(x))

(K::NfRel)(x::nf_elem) = K(base_field(K)(x))

(K::NfRel{T})(x::NfRelElem{T}) where {T} = K(x.data)

function (K::NfRel{NfRelElem{T}})(x::NfRelElem{T}) where {T}
  if parent(x) == base_field(K)
    return K(parent(K.pol)(x))
  end
  return force_coerce_throwing(K, x)
end

function (K::NfRel{nf_elem})(x::nf_elem)
  if parent(x) == base_field(K)
    return K(parent(K.pol)(x))
  end
  return force_coerce_throwing(K, x)
end

################################################################################
#
#  Degree
#
################################################################################

@inline degree(L::Hecke.NfRel) = degree(L.pol)

################################################################################
#
#  Reduction
#
################################################################################

function reduce!(a::NfRelElem)
  a.data = mod(a.data, parent(a).pol)
  return a
end

################################################################################
#
#  mod(::NfRelElem, ::fmpz) as in the absolute case
#
################################################################################

function mod(a::NfRelElem{T}, p::fmpz) where T
  K = parent(a)
  b = data(a)
  coeffs = Vector{T}(undef, degree(K)+1)
  for i = 0:degree(K)
    coeffs[i+1] = mod(coeff(b, i), p)
  end
  Kx = parent(b)
  return K(Kx(coeffs))
end

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, a::NfRel)
  print(io, "Relative number field over with defining polynomial ", a.pol)
  print(io, "over \n", a.base_ring)
end

function AbstractAlgebra.expressify(a::NfRelElem; context = nothing)
  return AbstractAlgebra.expressify(data(a), var(a.parent), context = context)
end

function Base.show(io::IO, a::NfRelElem)
  print(io, AbstractAlgebra.obj_to_string(a, context = io))
end

################################################################################
#
#  Constructors and parent object overloading
#
################################################################################

function NumberField(f::PolyElem{T}, s::String;
                     cached::Bool = false, check::Bool = true)  where {T <: NumFieldElem}
  S = Symbol(s)
  check && !isirreducible(f) && throw(error("Polynomial must be irreducible"))
  K = NfRel{T}(f, S, cached)
  return K, K(gen(parent(f)))
end

function NumberField(f::PolyElem{<: NumFieldElem}; cached::Bool = false, check::Bool = true)
  return NumberField(f, "_\$", cached = cached, check = check)
end

#Conversion to absolute non simple
function number_field(::Type{AnticNumberField}, L::NfRel{nf_elem})
  @assert degree(base_field(L)) == 1
  pol = isunivariate(map_coeffs(FlintQQ, L.pol))[2]
  return number_field(pol, check = false)
end

function (K::NfRel{T})(a::Generic.Poly{T}) where T
  z = NfRelElem{T}(mod(a, K.pol))
  z.parent = K
  return z
end

function (K::NfRel{T})(a::T) where T
  parent(a) != base_ring(parent(K.pol)) == error("Cannot coerce")
  z = NfRelElem{T}(parent(K.pol)(a))
  z.parent = K
  return z
end

(K::NfRel)(a::Integer) = K(parent(K.pol)(a))

(K::NfRel)(a::Rational{T}) where {T <: Integer} = K(parent(K.pol)(a))

(K::NfRel)(a::fmpz) = K(parent(K.pol)(a))

(K::NfRel)(a::fmpq) = K(parent(K.pol)(a))

(K::NfRel)() = zero(K)

gen(K::NfRel) = K(gen(parent(K.pol)))

################################################################################
#
#  Unary operators
#
################################################################################

function Base.:(-)(a::NfRelElem)
  return parent(a)(-data(a))
end

################################################################################
#
#  Binary operators
#
################################################################################

function Base.:(+)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  parent(a) == parent(b) || force_op(+, a, b)::NfRelElem{T}
  return parent(a)(data(a) + data(b))
end

function Base.:(-)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  parent(a) == parent(b) || force_op(-, a, b)::NfRelElem{T}
  return parent(a)(data(a) - data(b))
end

function Base.:(*)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  parent(a) == parent(b) || force_op(*, a, b)::NfRelElem{T}
  return parent(a)(data(a) * data(b))
end

function divexact(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  b == 0 && error("Element not invertible")
  parent(a) == parent(b) || force_op(divexact, a, b)::NfRelElem{T}
  return a*inv(b)
end

Base.:(//)(a::NfRelElem{T}, b::NfRelElem{T}) where {T} = divexact(a, b)

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

function Base.inv(a::NfRelElem)
  a == 0 && throw(error("Element not invertible"))
  g, s, _ = gcdx(data(a), parent(a).pol)
  @assert g == 1
  return parent(a)(s)
end

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

Base.:(^)(a::NfRelElem, n::UInt) = a^Int(n)

function Base.:(^)(a::NfRelElem, n::Int)
  K = parent(a)
  if iszero(a)
    return zero(K)
  end

  if n == 0
    return one(K)
  end

  if n < 0 && iszero(a)
    error("Element is not invertible")
  end

  return K(powmod(data(a), n, K.pol))
end

function Base.:(^)(a::NfRelElem, b::fmpz)
  if fits(Int, b)
    return a^Int(b)
  end
  if b < 0
    return inv(a)^(-b)
  elseif b == 0
    return parent(a)(1)
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  else
    #mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

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

function Base.:(==)(a::NfRelElem{T}, b::NfRelElem{T}) where T
  reduce!(a)
  reduce!(b)
  parent(a) == parent(b) || force_op(==, a, b)::Bool
  return data(a) == data(b)
end

if !isdefined(Nemo, :promote_rule1)
  function Base.:(==)(a::NfRelElem{T}, b::Union{Integer, Rational}) where T
    return a == parent(a)(b)
  end
end

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

function Base.:(*)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(data(a) * b)
end

Base.:(*)(a::T, b::NfRelElem{T}) where {T <: NumFieldElem} = b * a

function Base.:(+)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(data(a) + b)
end

Base.:(+)(a::T, b::NfRelElem{T}) where {T <: NumFieldElem} = b + a

function Base.:(-)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(data(a) - b)
end

function Base.:(-)(a::T, b::NfRelElem{T}) where {T <: NumFieldElem}
  return parent(b)(a - data(b))
end

function divexact(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(divexact(data(a), b))
end

Base.:(//)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem} = divexact(a, b)

for F in [fmpz, fmpq, Int]
  @eval begin
    function Base.:(*)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(data(a) * b)
    end

    Base.:(*)(a::$F, b::NfRelElem{T}) where {T <: NumFieldElem} = b * a

    function Base.:(+)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(data(a) + b)
    end

    Base.:(+)(a::$F, b::NfRelElem{T}) where {T <: NumFieldElem} = b + a

    function Base.:(-)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(data(a) - b)
    end

    function Base.:(-)(a::$F, b::NfRelElem{T}) where {T <: NumFieldElem}
      return parent(b)(a - data(b))
    end

    function divexact(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(divexact(data(a), b))
    end

    Base.:(//)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem} = divexact(a, b)
  end
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function mul!(c::NfRelElem{T}, a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  mul!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

function mul!(c::NfRelElem{T}, a::NfRelElem{T}, b::T) where {T}
  mul!(c.data, a.data, b)
  c = reduce!(c)
  return c
end

function addeq!(b::NfRelElem{T}, a::NfRelElem{T}) where {T}
  addeq!(b.data, a.data)
  b = reduce!(b)
  return b
end

function add!(c::NfRelElem{T}, a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  add!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

################################################################################
#
#  Isomorphism to absolute field (AnticNumberField)
#
################################################################################

#@doc Markdown.doc"""
#    absolute_field(K::NfRel{nf_elem}, cached::Bool = false) -> AnticNumberField, Map, Map
#Given an extension $K/k/Q$, find an isomorphic extension of $Q$.
#"""
function absolute_field(K::NfRel{nf_elem}; cached::Bool = false, simplify::Bool = false)
  if simplify
    return simplified_absolute_field(K, cached = cached)
  end
  Ka, a, b, c = _absolute_field(K, cached)
  h1 = hom(Ka, K, c, inverse = (a, b))
  h2 = hom(base_field(K), Ka, a, check = false)
  embed(h1)
  embed(MapFromFunc(x->preimage(h1, x), K, Ka))
  embed(h2)
  return Ka, h1, h2
end

#@doc Markdown.doc"""
#    absolute_field(K::NfRel{NfRelElem}, cached::Bool = false) -> NfRel, Map, Map
#Given an extension $E/K/k$, find an isomorphic extension of $k$.
#In a tower, only the top-most steps are collapsed.
#"""
function absolute_field(K::NfRel{NfRelElem{T}}, cached::Bool = false) where T
  Ka, a, b, c = _absolute_field(K)
  h1 = hom(Ka, K, c, inverse = (a, b))
  h2 = hom(base_field(K), Ka, a, check = false)
  embed(h1)
  embed(MapFromFunc(x->preimage(h1, x), K, Ka))
  embed(h2)
  return Ka, h1, h2
end

#Trager: p4, Algebraic Factoring and Rational Function Integration
function _absolute_field(K::NfRel, cached::Bool = false)
  f = K.pol
  kx = parent(f)
  k = base_ring(kx)
  Qx = parent(k.pol)

  l = 0
  g = f
  N = norm(g)

  while true
    @assert degree(N) == degree(g) * degree(k)

    if !isconstant(N) && issquarefree(N)
      break
    end

    l += 1

    g = compose(f, gen(kx) - l*gen(k))
    N = norm(g)
  end

  Ka, gKa = NumberField(N, "x", cached = cached, check = false)

  KaT, T = PolynomialRing(Ka, "T", cached = false)

  # map Ka -> K: gen(Ka) -> gen(K)+ k gen(k)

  # gen(k) -> Root(gcd(g, poly(k)))  #gcd should be linear:
  # g in kx = (Q[a])[x]. Want to map x -> gen(Ka), a -> T

  gg = zero(KaT)
  for i=degree(g):-1:0
    auxp = change_base_ring(Ka, Qx(coeff(g, i)), parent = KaT)
    gg = gg*gKa
    add!(gg, gg, auxp)
    #gg = gg*gKa + auxp
  end

  q = gcd(gg, change_base_ring(Ka, k.pol, parent = KaT))
  @assert degree(q) == 1
  al = -trailing_coefficient(q)//lead(q)
  be = gKa - l*al
  ga = gen(K) + l*gen(k)

  #al -> gen(k) in Ka
  #be -> gen(K) in Ka
  #ga -> gen(Ka) in K
  return Ka, al, be, ga
end

function simplified_absolute_field(K::NfRel{nf_elem}; cached::Bool = false)
  Kabs, mKabs, mk = absolute_field(K, false)
  OK = maximal_order(K)
  new_basis = Vector{nf_elem}(undef, degree(Kabs))
  B = pseudo_basis(OK)
  ideals = Dict{NfOrdIdl, Vector{nf_elem}}()
  for i = 1:length(B)
    I = B[i][2].num
    if !haskey(ideals, I)
      bas = lll_basis(I)
      ideals[I] = nf_elem[mKabs\(K(x)) for x in bas]
    end
  end
  ind = 1
  for i = 1:degree(K)
    I = B[i][2]
    bI = ideals[I.num]
    el = mKabs\(B[i][1])
    for j = 1:length(bI)
      new_basis[ind] = divexact(el*bI[j], I.den)
      ind += 1
    end
  end
  O = NfOrd(new_basis)
  O.disc = absolute_discriminant(OK)
  O.ismaximal = 1
  OLLL = lll(O)
  _set_maximal_order_of_nf(Kabs, OLLL)
  Ks, mKs = simplify(Kabs)
  mKsi = inv(mKs)
  return Ks, mKs*mKabs, mk*mKsi
end

function check_parent(a, b)
  return a==b
end

function hash(a::Hecke.NfRelElem{nf_elem}, b::UInt)
  return hash(a.data, b)
end

################################################################################
#
#  Coercion and in function
#
################################################################################

function (K::AnticNumberField)(a::NfRelElem{nf_elem})
  K != base_field(parent(a)) && return force_coerce_throwing(K, a)
  for i in 2:degree(parent(a))
    @assert iszero(coeff(a, i - 1))
  end
  return coeff(a, 0)
end

function in(a::NfRelElem{nf_elem}, K::AnticNumberField)
  L = parent(a)
  @assert base_field(L) == K
  for i in 2:degree(parent(a))
    if !iszero(coeff(a, i - 1))
      return false
    end
  end
  return true
end

################################################################################
#
#  Basis & representation matrix
#
################################################################################

function basis_matrix(v::Vector{<: NfRelElem})
  K = parent(v[1])
  k = base_field(K)
  z = zero_matrix(k, length(v), degree(K))
  for i in 1:length(v)
    for j in 0:(degree(K) - 1)
      z[i, j + 1] = coeff(v[i], j)
    end
  end
  return z
end

function elem_to_mat_row!(M::Generic.Mat{T}, i::Int, a::NfRelElem{T}) where T
  for c = 1:ncols(M)
    M[i, c] = deepcopy(coeff(a, c - 1))
  end
  return nothing
end

function elem_from_mat_row(L::NfRel{T}, M::Generic.Mat{T}, i::Int) where T
  t = L(1)
  a = L()
  for c = 1:ncols(M)
    a += M[i, c]*t
    mul!(t, t, gen(L))
  end
  return a
end

function representation_matrix(a::NfRelElem)
  L = a.parent
  n = degree(L)
  M = zero_matrix(base_field(L), n, n)
  t = gen(L)
  b = deepcopy(a)
  for i = 1:n-1
    elem_to_mat_row!(M, i, b)
    mul!(b, b, t)
  end
  elem_to_mat_row!(M, n, b)
  return M
end

function norm(a::NfRelElem{nf_elem}, new::Bool = !true)
  if new && ismonic(parent(a).pol) #should be much faster - eventually
    return resultant_mod(a.data, parent(a).pol)
  end
  M = representation_matrix(a)
  return det(M)
end

function norm(a::NfRelElem, new::Bool = true)
  if new && ismonic(parent(a).pol)
    return resultant(a.data, parent(a).pol)
  end
  M = representation_matrix(a)
  return det(M)
end

################################################################################
#
#  Trace
#
################################################################################

function assure_trace_basis(K::NfRel)
  if isdefined(K, :trace_basis)
    return nothing
  end
  F = base_field(K)
  trace_basis = Vector{elem_type(F)}(undef, degree(K))
  trace_basis[1] = F(degree(K))
  a = gen(K)
  for i = 2:degree(K)
    #We can do better, probably...
    M = representation_matrix(a)
    trace_basis[i] = tr(M)
    a *= gen(K)
  end
  K.trace_basis = trace_basis
  return nothing
end

function tr(a::NfRelElem)
  K = parent(a)
  assure_trace_basis(K)
  t = coeff(a, 0)*K.trace_basis[1]
  for i = 2:degree(K)
    c = coeff(a, i-1)
    if !iszero(c)
      t += c*K.trace_basis[i]
    end
  end
  return t
end

function tr(a::NfRelElem{nf_elem})
  K = parent(a)
  assure_trace_basis(K)
  t = coeff(a, 0)*K.trace_basis[1]
  for i = 2:degree(K)
    c = coeff(a, i-1)
    if !iszero(c)
      add!(t, t, c*K.trace_basis[i])
    end
  end
  return t
end

################################################################################
#
#  Minimal and characteristic polynomial
#
################################################################################

function _poly_norm_to(f, k::T) where {T}
  if base_ring(f) isa T
    @assert base_ring(f) == k
    return f
  else
    return _poly_norm_to(norm(f), k)
  end
end

#TODO: investigate charpoly/ minpoly from power_sums, aka tr(a^i) and
#      Newton identities
#TODO: cache traces of powers of the generator on the field, then
#      the trace does not need the matrix

@doc Markdown.doc"""
    charpoly(a::NfRelElem) -> PolyElem

Given an element $a$ in an extension $L/K$, this function returns the
characteristic polynomial of $a$ over $K$.
"""
function charpoly(a::NfRelElem)
  M = representation_matrix(a)
  R = PolynomialRing(base_field(parent(a)), cached = false)[1]
  return charpoly(R, M)
end

@doc Markdown.doc"""
    minpoly(a::NfRelElem) -> PolyElem

Given an element $a$ in an extension $L/K$, this function returns the minimal
polynomial of $a$ over $K$.
"""
function minpoly(a::NfRelElem{S}) where {S}
  M = representation_matrix(a)
  R = PolynomialRing(base_field(parent(a)), cached = false)[1]
  return minpoly(R, M, false)::Generic.Poly{S}
end

function charpoly(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})
  f = charpoly(a)
  return _poly_norm_to(f, k)
end

function absolute_charpoly(a::NfRelElem)
  return charpoly(a, FlintQQ)
end

function minpoly(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})

  if parent(a) == k
    R, y  = PolynomialRing(k, cached = false)
    return y - a
  end

  f = minpoly(a)
  while base_ring(f) != k && !(base_ring(f) isa FlintRationalField && k isa FlintRationalField)
    f = norm(f)
    g = gcd(f, derivative(f))
    if !isone(g)
      f = divexact(f, g)
    end
  end
  return f
end

function absolute_minpoly(a::NfRelElem)
  return minpoly(a, FlintQQ)
end

norm(a::fmpq_poly) = a

#

#

function (R::Generic.PolyRing{nf_elem})(a::NfRelElem{nf_elem})
  if base_ring(R)==base_field(parent(a))
    return a.data
  end
  error("wrong ring")
end

function (R::Generic.PolyRing{NfRelElem{T}})(a::NfRelElem{NfRelElem{T}}) where T
  if base_ring(R)==base_field(parent(a))
    return a.data
  end
  error("wrong ring")
end

################################################################################
#
#  issubfield and isisomorphic
#
################################################################################

function issubfield(K::NfRel, L::NfRel)
  @assert base_field(K) == base_field(L)
  f = K.pol
  g = L.pol
  if mod(degree(g), degree(f)) != 0
    return false, hom(K, L, zero(L), check = false)
  end
  Lx, x = PolynomialRing(L, "x", cached = false)
  fL = Lx()
  for i = 0:degree(f)
    setcoeff!(fL, i, L(coeff(f, i)))
  end
  fac = factor(fL)
  for (a, b) in fac
    if degree(a) == 1
      c1 = coeff(a, 0)
      c2 = coeff(a, 1)
      h = parent(K.pol)(-c1*inv(c2))
      return true, hom(K, L, h(gen(L)), check = false)
    end
  end
  return false, hom(K, L, zero(L), check = false)
end

function isisomorphic(K::NfRel, L::NfRel)
  @assert base_field(K) == base_field(L)
  f = K.pol
  g = L.pol
  if degree(f) != degree(g)
    return false, hom(K, L, zero(L), check = false)
  end
  return issubfield(K, L)
end

################################################################################
#
#  Normal basis
#
################################################################################

# Mostly the same as in the absolute case
function normal_basis(L::NfRel{nf_elem}, check::Bool = false)
  O = EquationOrder(L)
  K = base_field(L)
  OK = base_ring(O)
  d = discriminant(O)

  for p in PrimeIdealsSet(OK, degree(L), -1, indexdivisors = false, ramified = false)
    if valuation(d, p) != 0
      continue
    end

    # Check if p is totally split
    F, mF = ResidueField(OK, p)
    mmF = extend(mF, K)
    Ft, t = PolynomialRing(F, "t", cached = false)
    ft = map_coeffs(mmF, L.pol, parent = Ft)
    pt = powmod(t, order(F), ft)

    if degree(gcd(ft, pt - t)) == degree(ft)
      # Lift an idempotent of O/pO
      immF = pseudo_inv(mmF)
      fac = factor(ft)
      gt = divexact(ft, first(keys(fac.fac)))
      g = fq_poly_to_nf_elem_poly(parent(L.pol), immF, gt)
      return L(g)
    end
  end
end

################################################################################
#
#  Linear disjointness
#
################################################################################

function islinearly_disjoint(K1::NfRel, K2::NfRel)
  if base_field(K1) != base_field(K2)
    throw(error("Number fields must have the same base field"))
  end

  if gcd(degree(K1), degree(K2)) == 1
    return true
  end

  f = change_base_ring(K2, defining_polynomial(K1))
  return isirreducible(f)
end

################################################################################
#
#  Random elements
#
################################################################################

RandomExtensions.maketype(L::NfRel, B) = elem_type(L)

function rand(rng::AbstractRNG, sp::SamplerTrivial{<:Make2{<:NfRelElem,<:NfRel,<:UnitRange}})
  L, B = sp[][1:end]
  k = base_field(L)
  pb = basis(L)
  z = zero(L)
  for i = 1:length(pb)
    t = rand(rng, k, B)
    z += t*pb[i]
  end
  return z
end

rand(L::NfRel, B::UnitRange{Int}) = rand(GLOBAL_RNG, L, B)
rand(rng::AbstractRNG, L::NfRel, B::UnitRange{Int}) = rand(rng, make(L, B))

################################################################################
#
#  Find Kummer generator
#
################################################################################

@doc Markdown.doc"""
    kummer_generator(K::NfRel{nf_elem}) -> nf_elem

Given an extension $K/k$ which is a cyclic Kummer extension of degree $n$, returns an element $a\in k$
such that $K = k(\sqrt[n]{a})$. Throws an error if the extension is not a cyclic Kummer extension.
"""
function kummer_generator(K::NfRel{nf_elem})
  n = degree(K)
  k = base_field(K)
  tuo, gen_tu = _torsion_units_gen(k)
  if !divisible(tuo, n)
    error("Not a Kummer extension!")
  end
  zeta = gen_tu^divexact(tuo, n)
  roots = powers(zeta, n-1)
  auts = automorphisms(K)
  if length(auts) != n
    error("Not a Kummer extension!")
  end
  #TODO: Not very serious, fix this.
  gens = small_generating_set(auts, *)
  if length(gens) > 1
    error("Not a Kummer extension!")
  end
  gen_aut = gens[1]
  a = zero(K)
  stable = 10
  inter = 1
  first = true
  while iszero(a)
    stable -= 1
    if iszero(stable)
      inter += 1
      stable = 10
    end
    if first
      first = false
      b = gen(K)
    else
      b = rand(K, -inter:inter)
    end
    a = b
    new_b = b
    for i = 1:n-1
      new_b = gen_aut(new_b)
      a += roots[i+1]*new_b
    end

  end
  res = k(a^n)
  #We even reduce the support....
  res1 = reduce_mod_powers(res, n)
  return res1
end

################################################################################
#
#  Relative extension
#
################################################################################

#TODO: Put some more thought in it.
@doc Markdown.doc"""
    relative_extension(K::AnticNumberField, k::AnticNumberField) -> NfRel{nf_elem}

Given two fields $K\supset k$, it returns $K$ as a relative
extension of $k$ and an isomorphism between it and $K$.
"""
function relative_extension(m::NfToNfMor)
  k = domain(m)
  K = codomain(m)
  lf = factor(K.pol, k)
  rel_deg = divexact(degree(K), degree(k))
  pols = [f for (f, v) in lf if degree(f) == rel_deg]
  p = pols[1]
  if length(pols) > 1
    i = 2
    while !iszero(map_coeffs(m, p)(gen(K)))
      p = pols[i]
      i += 1
    end
  end
  L, b = number_field(p, cached = false, check = false)
  mp = hom(K, L, b, inverse = (image_primitive_element(m), gen(K)))
  return L, mp
end

function relative_extension(K::AnticNumberField, k::AnticNumberField)
  fl, mp = issubfield(k, K)
  if !fl
    error("Not a subfield!")
  end
  return relative_extension(mp)
end

################################################################################
#
#  Signature
#
################################################################################

function signature(L::NfRel)
  c = get_special(L, :signature)
  if c isa Tuple{Int, Int}
    return c::Tuple{Int, Int}
  end
  K = base_field(L)
  rlp = real_places(K)
  rL = 0
  for P in rlp
    rL += number_real_roots(defining_polynomial(L), P)
  end
  @assert mod(absolute_degree(L) - rL, 2) == 0
  r, s = rL, div(absolute_degree(L) - rL, 2)
  set_special(L, :signature => (r, s))
  return r, s
end

function istotally_real(L::NfRel)
  r, s = signature(L)
  return s == 0
end

function istotally_complex(L::NfRel)
  r, s = signature(L)
  return r == 0
end


text:
@testset "translate.jl" begin

    @testset "scale=1" begin

        map = Map_Translate()

        @test parameter_size(map) == 1

        θ = Float32[3]
        X_hat = Int[1,2,3]
        X = eval_map(map,X_hat,θ)

        @test eltype(X) == Float32
        @test X ≈ Float32[1,2,3] .+ (first(θ)-1)

    end

    @testset "scale=10" begin

        map = Map_Translate(10.0)

        @test parameter_size(map) == 1

        θ = Float32[3]
        X_hat = Int[1,2,3]
        X = eval_map(map,X_hat,θ)

        @test eltype(X) == Float64
        @test X ≈ Float64[1,2,3] .+ 10 * (first(θ)-1)

    end

end


text:
#####
##### LPCM API types/functions/stubs
#####

const LPCM_FORMAT_REGISTRY = Any[]

"""
    Onda.register_lpcm_format!(create_constructor)

Register an `AbstractLPCMFormat` constructor so that it can automatically be used when
[`format`](@ref) is called. Authors of new `AbstractLPCMFormat` subtypes should call this
function for their subtype.

`create_constructor` should be a unary function that accepts a single `file_format::AbstractString`
argument, and return either a matching `AbstractLPCMFormat` constructor or `nothing`. Any returned
`AbstractLPCMFormat` constructor `f` should be of the form `f(info::SamplesInfo; kwargs...)::AbstractLPCMFormat`.

Note that if `Onda.register_lpcm_format!` is called in a downstream package, it must be called
within the `__init__` function of the package's top-level module to ensure that the function
is always invoked when the module is loaded (not just during precompilation). For details,
see https://docs.julialang.org/en/v1/manual/modules/#Module-initialization-and-precompilation.
"""
register_lpcm_format!(create_constructor) = push!(LPCM_FORMAT_REGISTRY, create_constructor)

"""
    format(file_format::AbstractString, info::SamplesInfo; kwargs...)

Return `f(info; kwargs...)` where `f` constructs the `AbstractLPCMFormat` instance that
corresponds to `file_format`. `f` is determined by matching `file_format` to a suitable
format constuctor registered via [`register_lpcm_format!`](@ref).

See also: [`deserialize_lpcm`](@ref), [`serialize_lpcm`](@ref)
"""
function format(file_format::AbstractString, info::SamplesInfo; kwargs...)
    for create_constructor in LPCM_FORMAT_REGISTRY
        f = create_constructor(file_format)
        f === nothing && continue
        return f(info; kwargs...)
    end
    throw(ArgumentError("unrecognized file_format: \"$file_format\""))
end

"""
    AbstractLPCMFormat

A type whose subtypes represents byte/stream formats that can be (de)serialized
to/from Onda's standard interleaved LPCM representation.

All subtypes of the form `F<:AbstractLPCMFormat` must call [`Onda.register_lpcm_format!`](@ref)
and define an appropriate [`file_format_string`](@ref) method.

See also:

- [`format`](@ref)
- [`deserialize_lpcm`](@ref)
- [`deserialize_lpcm_callback`](@ref)
- [`serialize_lpcm`](@ref)
- [`LPCMFormat`](@ref)
- [`LPCMZstFormat`](@ref)
- [`AbstractLPCMStream`](@ref)
"""
abstract type AbstractLPCMFormat end

"""
    AbstractLPCMStream

A type that represents an LPCM (de)serialization stream.

See also:

- [`deserializing_lpcm_stream`](@ref)
- [`serializing_lpcm_stream`](@ref)
- [`finalize_lpcm_stream`](@ref)
"""
abstract type AbstractLPCMStream end

"""
    deserialize_lpcm_callback(format::AbstractLPCMFormat, samples_offset, samples_count)

Return `(callback, required_byte_offset, required_byte_count)` where `callback` accepts the
byte block specified by `required_byte_offset` and `required_byte_count` and returns the
samples specified by `samples_offset` and `samples_count`.

As a fallback, this function returns `(callback, missing, missing)`, where `callback`
requires all available bytes. `AbstractLPCMFormat` subtypes that support partial/block-based
deserialization (e.g. the basic `LPCMFormat`) can overload this function to only request
exactly the byte range that is required for the sample range requested by the caller.

This allows callers to handle the byte block retrieval themselves while keeping
Onda's LPCM Serialization API agnostic to the caller's storage layer of choice.
"""
function deserialize_lpcm_callback(format::AbstractLPCMFormat, samples_offset, samples_count)
    callback = bytes -> deserialize_lpcm(format, bytes, samples_offset, samples_count)
    return callback, missing, missing
end

"""
    deserializing_lpcm_stream(format::AbstractLPCMFormat, io)

Return a `stream::AbstractLPCMStream` that wraps `io` to enable direct LPCM
deserialization from `io` via [`deserialize_lpcm`](@ref).

Note that `stream` must be finalized after usage via [`finalize_lpcm_stream`](@ref).
Until `stream` is finalized, `io` should be considered to be part of the internal
state of `stream` and should not be directly interacted with by other processes.
"""
function deserializing_lpcm_stream end

"""
    serializing_lpcm_stream(format::AbstractLPCMFormat, io)

Return a `stream::AbstractLPCMStream` that wraps `io` to enable direct LPCM
serialization to `io` via [`serialize_lpcm`](@ref).

Note that `stream` must be finalized after usage via [`finalize_lpcm_stream`](@ref).
Until `stream` is finalized, `io` should be considered to be part of the internal
state of `stream` and should not be directly interacted with by other processes.
"""
function serializing_lpcm_stream end

"""
    finalize_lpcm_stream(stream::AbstractLPCMStream)::Bool

Finalize `stream`, returning `true` if the underlying I/O object used to construct
`stream` is still open and usable. Otherwise, return `false` to indicate that
underlying I/O object was closed as result of finalization.
"""
function finalize_lpcm_stream end

"""
    deserialize_lpcm(format::AbstractLPCMFormat, bytes,
                     samples_offset::Integer=0,
                     samples_count::Integer=typemax(Int))
    deserialize_lpcm(stream::AbstractLPCMStream,
                     samples_offset::Integer=0,
                     samples_count::Integer=typemax(Int))

Return a channels-by-timesteps `AbstractMatrix` of interleaved LPCM-encoded
sample data by deserializing the provided `bytes` in the given `format`, or
from the given `stream` constructed by [`deserializing_lpcm_stream`](@ref).

Note that this operation may be performed in a zero-copy manner such that the
returned sample matrix directly aliases `bytes`.

The returned segment is at most `sample_offset` samples offset from the start of
`stream`/`bytes` and contains at most `sample_count` samples. This ensures that
overrun behavior is generally similar to the behavior of `Base.skip(io, n)` and
`Base.read(io, n)`.

This function is the inverse of the corresponding [`serialize_lpcm`](@ref) method, i.e.:

```
serialize_lpcm(format, deserialize_lpcm(format, bytes)) == bytes
```
"""
function deserialize_lpcm end

"""
    serialize_lpcm(format::AbstractLPCMFormat, samples::AbstractMatrix)
    serialize_lpcm(stream::AbstractLPCMStream, samples::AbstractMatrix)

Return the `AbstractVector{UInt8}` of bytes that results from serializing `samples`
to the given `format` (or serialize those bytes directly to `stream`) where `samples`
is a channels-by-timesteps matrix of interleaved LPCM-encoded sample data.

Note that this operation may be performed in a zero-copy manner such that the
returned `AbstractVector{UInt8}` directly aliases `samples`.

This function is the inverse of the corresponding [`deserialize_lpcm`](@ref)
method, i.e.:

```
deserialize_lpcm(format, serialize_lpcm(format, samples)) == samples
```
"""
function serialize_lpcm end

"""
    file_format_string(format::AbstractLPCMFormat)

Return the `String` representation of `format` to be written to the `file_format` field of a `*.signals` file.
"""
function file_format_string end

#####
##### read_lpcm/write_lpcm
#####

read_lpcm(path, format::AbstractLPCMFormat) = deserialize_lpcm(format, read(path))

function read_lpcm(path, format::AbstractLPCMFormat, sample_offset, sample_count)
    deserialize_requested_samples,
    required_byte_offset,
    required_byte_count = deserialize_lpcm_callback(format,
                                                    sample_offset,
                                                    sample_count)
    bytes = read_byte_range(path, required_byte_offset, required_byte_count)
    return deserialize_requested_samples(bytes)
end

write_lpcm(path, format::AbstractLPCMFormat, data) = write_full_path(path, serialize_lpcm(format, data))

#####
##### `LPCMFormat`
#####

"""
    LPCMFormat(channel_count::Int, sample_type::Type)
    LPCMFormat(info::SamplesInfo)

Return a `LPCMFormat<:AbstractLPCMFormat` instance corresponding to Onda's default
interleaved LPCM format assumed for sample data files with the "lpcm"
extension.

`channel_count` corresponds to `length(info.channels)`, while `sample_type`
corresponds to `sample_type(info)`

Note that bytes (de)serialized to/from this format are little-endian (per the
Onda specification).
"""
struct LPCMFormat{S<:LPCM_SAMPLE_TYPE_UNION} <: AbstractLPCMFormat
    channel_count::Int
    sample_type::Type{S}
end

LPCMFormat(info::SamplesInfo) = LPCMFormat(length(info.channels), sample_type(info))

register_lpcm_format!(file_format -> file_format == "lpcm" ? LPCMFormat : nothing)

file_format_string(::LPCMFormat) = "lpcm"

function _validate_lpcm_samples(format::LPCMFormat{S}, samples::AbstractMatrix) where {S}
    if format.channel_count != size(samples, 1)
        throw(ArgumentError("""
                            `samples` row count ($(size(samples, 1))) does not
                            match expected channel count ($(format.channel_count))
                            """))
    elseif !(eltype(samples) <: S)
        throw(ArgumentError("""
                            `samples` eltype ($(eltype(samples))) does not
                            match expected eltype ($S)
                            """))
    end
    return nothing
end

_bytes_per_sample(format::LPCMFormat{S}) where {S} = sizeof(S) * format.channel_count

struct LPCMStream{S<:LPCM_SAMPLE_TYPE_UNION,I} <: AbstractLPCMStream
    format::LPCMFormat{S}
    io::I
end

deserializing_lpcm_stream(format::LPCMFormat, io) = LPCMStream(format, io)

serializing_lpcm_stream(format::LPCMFormat, io) = LPCMStream(format, io)

finalize_lpcm_stream(::LPCMStream) = true

function deserialize_lpcm(format::LPCMFormat{S}, bytes, sample_offset::Integer=0,
                          sample_count::Integer=typemax(Int)) where {S}
    sample_interpretation = reinterpret(S, bytes)
    sample_start = min((format.channel_count * sample_offset) + 1, length(sample_interpretation))
    sample_end = format.channel_count * (sample_offset + sample_count)
    sample_end = sample_end >= 0 ? sample_end : typemax(Int) # handle overflow
    sample_end = min(sample_end, length(sample_interpretation))
    sample_view = view(sample_interpretation, sample_start:sample_end)
    timestep_count = min(Int(length(sample_view) / format.channel_count), sample_count)
    return reshape(sample_view, (format.channel_count, timestep_count))
end

function deserialize_lpcm_callback(format::LPCMFormat{S}, samples_offset, samples_count) where {S}
    callback = bytes -> deserialize_lpcm(format, bytes)
    bytes_per_sample = _bytes_per_sample(format)
    return callback, samples_offset * bytes_per_sample, samples_count * bytes_per_sample
end

function deserialize_lpcm(stream::LPCMStream, sample_offset::Integer=0,
                          sample_count::Integer=typemax(Int))
    bytes_per_sample = _bytes_per_sample(stream.format)
    jump(stream.io, bytes_per_sample * sample_offset)
    byte_count = bytes_per_sample * sample_count
    byte_count = byte_count >= 0 ? byte_count : typemax(Int) # handle overflow
    return deserialize_lpcm(stream.format, read(stream.io, byte_count))
end

function serialize_lpcm(format::LPCMFormat, samples::AbstractMatrix)
    _validate_lpcm_samples(format, samples)
    samples isa Matrix && return reinterpret(UInt8, vec(samples))
    io = IOBuffer()
    write(io, samples)
    return resize!(io.data, io.size)
end

function serialize_lpcm(stream::LPCMStream, samples::AbstractMatrix)
    _validate_lpcm_samples(stream.format, samples)
    return write(stream.io, samples)
end

#####
##### `LPCMZstFormat`
#####

"""
    LPCMZstFormat(lpcm::LPCMFormat; level=3)
    LPCMZstFormat(info::SamplesInfo; level=3)

Return a `LPCMZstFormat<:AbstractLPCMFormat` instance that corresponds to Onda's
default interleaved LPCM format compressed by `zstd`. This format is assumed
for sample data files with the "lpcm.zst" extension.

The `level` keyword argument sets the same compression level parameter as the
corresponding flag documented by the `zstd` command line utility.

See https://facebook.github.io/zstd/ for details about `zstd`.
"""
struct LPCMZstFormat{S} <: AbstractLPCMFormat
    lpcm::LPCMFormat{S}
    level::Int
    LPCMZstFormat(lpcm::LPCMFormat{S}; level=3) where {S} = new{S}(lpcm, level)
end

LPCMZstFormat(info::SamplesInfo; kwargs...) = LPCMZstFormat(LPCMFormat(info); kwargs...)

register_lpcm_format!(file_format -> file_format == "lpcm.zst" ? LPCMZstFormat : nothing)

file_format_string(::LPCMZstFormat) = "lpcm.zst"

function deserialize_lpcm(format::LPCMZstFormat, bytes, args...)
    decompressed_bytes = zstd_decompress(unsafe_vec_uint8(bytes))
    return deserialize_lpcm(format.lpcm, decompressed_bytes, args...)
end

function serialize_lpcm(format::LPCMZstFormat, samples::AbstractMatrix)
    decompressed_bytes = unsafe_vec_uint8(serialize_lpcm(format.lpcm, samples))
    return zstd_compress(decompressed_bytes, format.level)
end

struct LPCMZstStream{L<:LPCMStream} <: AbstractLPCMStream
    stream::L
end

function deserializing_lpcm_stream(format::LPCMZstFormat, io)
    stream = LPCMStream(format.lpcm, ZstdDecompressorStream(io))
    return LPCMZstStream(stream)
end

function serializing_lpcm_stream(format::LPCMZstFormat, io)
    stream = LPCMStream(format.lpcm, ZstdCompressorStream(io; level=format.level))
    return LPCMZstStream(stream)
end

function finalize_lpcm_stream(stream::LPCMZstStream)
    if stream.stream.io isa ZstdCompressorStream
        # write `TranscodingStreams.TOKEN_END` and change the `ZstdCompressorStream`'s
        # mode to `:close`, which flushes any remaining buffered data and finalizes the
        # underlying codec to free its resources without closing the underlying I/O object.
        write(stream.stream.io, TranscodingStreams.TOKEN_END)
        TranscodingStreams.changemode!(stream.stream.io, :close)
        return true
    else
        close(stream.stream.io)
        return false
    end
end

deserialize_lpcm(stream::LPCMZstStream, args...) = deserialize_lpcm(stream.stream, args...)

serialize_lpcm(stream::LPCMZstStream, args...) = serialize_lpcm(stream.stream, args...)


text:
using Test
using PETSc

@testset "options tests" begin
    kw_opts = (
        ksp_monitor = nothing,
        ksp_view = true,
        false_opt = false,
        da_grid_x = 100,
        da_grid_y = 100,
        pc_type = "mg",
        pc_mg_levels = 1,
        mg_levels_0_pc_type = "ilu",
    )
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)

        opts = PETSc.Options(petsclib; kw_opts...)

        # Check that all the keys got added
        for (key, val) in pairs(kw_opts)
            key = string(key)
            if val === false
                @test_throws KeyError opts[key]
            else
                val = isnothing(val) || val === true ? "" : string(val)
                @test val == opts[key]
            end
        end

        # Check that we throw when a bad key is asked for
        @test_throws KeyError opts["bad key"]

        # try to insert some new keys
        opts["new_opt"] = 1
        opts["nothing_opt"] = nothing
        @test "" == opts["nothing_opt"]
        @test "1" == opts["new_opt"]

        # Check that viewer is working
        _stdout = stdout
        (rd, wr) = redirect_stdout()
        @show opts

        @test readline(rd) == "opts = #PETSc Option Table entries:"
        @test readline(rd) == "-da_grid_x 100"
        @test readline(rd) == "-da_grid_y 100"
        @test readline(rd) == "-ksp_monitor"
        @test readline(rd) == "-ksp_view"
        @test readline(rd) == "-mg_levels_0_pc_type ilu"
        @test readline(rd) == "-new_opt 1"
        @test readline(rd) == "-nothing_opt"
        @test readline(rd) == "-pc_mg_levels 1"
        @test readline(rd) == "-pc_type mg"
        @test readline(rd) == "#End of PETSc Option Table entries"
        @test readline(rd) == ""

        glo_opts = PETSc.GlobalOptions(petsclib)
        show(stdout, "text/plain", glo_opts)
        @test readline(rd) == "#No PETSc Option Table entries"

        # Try to set some options and check that they are set
        PETSc.with(opts) do
            show(stdout, "text/plain", glo_opts)
        end
        @test readline(rd) == "#PETSc Option Table entries:"
        @test readline(rd) == "-da_grid_x 100"
        @test readline(rd) == "-da_grid_y 100"
        @test readline(rd) == "-ksp_monitor"
        @test readline(rd) == "-ksp_view"
        @test readline(rd) == "-mg_levels_0_pc_type ilu"
        @test readline(rd) == "-new_opt 1"
        @test readline(rd) == "-nothing_opt"
        @test readline(rd) == "-pc_mg_levels 1"
        @test readline(rd) == "-pc_type mg"
        @test readline(rd) == "#End of PETSc Option Table entries"

        show(stdout, "text/plain", glo_opts)
        @test readline(rd) == "#No PETSc Option Table entries"

        redirect_stdout(_stdout)

        PETSc.finalize(petsclib)
    end
end

@testset "parse_options tests" begin
    @test begin
        julia = joinpath(Sys.BINDIR, Base.julia_exename())
        run(`$(julia) --startup-file=no --project -e "using PETSc
                                 using Test
                                 opts = PETSc.parse_options(ARGS)
                                 @test length(opts) == 4
                                 @test opts.ksp_monitor === nothing
                                 @test opts.ksp_view === nothing
                                 @test opts.pc_type === \"mg\"
                                 @test opts.da_grid_x === \"100\"
                                 @test_throws ErrorException opts.bad_opt
                                 "
                                 --
                                 -ksp_monitor
                                 -da_grid_x 100
                                 -ksp_view
                                 -pc_type mg`)
        true
    end
end


text:
@testset "Connectivity" begin
    g6 = smallgraph(:house)
    gx = PathGraph(4)
    add_vertices!(gx,10)
    add_edge!(gx,5,6)
    add_edge!(gx,6,7)
    add_edge!(gx,8,9)
    add_edge!(gx,10,9)

    for g in testgraphs(gx)
      @test @inferred(!is_connected(g))
      cc = @inferred(connected_components(g))
      label = zeros(eltype(g), nv(g))
      @inferred(LightGraphs.connected_components!(label, g))
      @test label[1:10] == [1,1,1,1,5,5,5,8,8,8]
      import LightGraphs: components, components_dict
      cclab = @inferred(components_dict(label))
      @test cclab[1] == [1,2,3,4]
      @test cclab[5] == [5,6,7]
      @test cclab[8] == [8,9,10]
      @test length(cc) >= 3 && sort(cc[3]) == [8,9,10]
    end
    for g in testgraphs(g6)
      @test @inferred(is_connected(g))
    end

    g10 = DiGraph(4)
    add_edge!(g10,1,3)
    add_edge!(g10,2,4)
    for g in testdigraphs(g10)
      @test @inferred(is_bipartite(g))
    end
    add_edge!(g10,1,4)
    for g in testdigraphs(g10)
      @test @inferred(is_bipartite(g))
    end

    g10 = DiGraph(20)
    for g in testdigraphs(g10)
      for m=1:50
        i = rand(1:10)
        j = rand(11:20)
        if rand() < 0.5
            i, j = j, i
        end
        if !has_edge(g, i, j)
            add_edge!(g, i, j)
            @test @inferred(is_bipartite(g))
        end
      end
    end

    # graph from https://en.wikipedia.org/wiki/Strongly_connected_component
    h = DiGraph(8)
    add_edge!(h,1,2); add_edge!(h,2,3); add_edge!(h,2,5);
    add_edge!(h,2,6); add_edge!(h,3,4); add_edge!(h,3,7);
    add_edge!(h,4,3); add_edge!(h,4,8); add_edge!(h,5,1);
    add_edge!(h,5,6); add_edge!(h,6,7); add_edge!(h,7,6);
    add_edge!(h,8,4); add_edge!(h,8,7)
    for g in testdigraphs(h)
      @test @inferred(is_weakly_connected(g))
      scc = @inferred(strongly_connected_components(g))
      wcc = @inferred(weakly_connected_components(g))

      @test length(scc) == 3 && sort(scc[3]) == [1,2,5]
      @test length(wcc) == 1 && length(wcc[1]) == nv(g)

    end

    function scc_ok(graph)
      """Check that all SCC really are strongly connected"""
      scc = @inferred(strongly_connected_components(graph))
      scc_as_subgraphs = map(i -> graph[i], scc)
      return all(is_strongly_connected, scc_as_subgraphs)
    end

    # the two graphs below are isomorphic (exchange 2 <--> 4)
    h = DiGraph(4);  add_edge!(h, 1, 4); add_edge!(h, 4, 2); add_edge!(h, 2, 3); add_edge!(h, 1, 3);
    for g in testdigraphs(h)
      @test scc_ok(g)
    end

    h2 = DiGraph(4); add_edge!(h2, 1, 2); add_edge!(h2, 2, 4); add_edge!(h2, 4, 3); add_edge!(h2, 1, 3);
    for g in testdigraphs(h2)
      @test scc_ok(g)
    end

    h = DiGraph(6)
    add_edge!(h,1,3); add_edge!(h,3,4); add_edge!(h,4,2); add_edge!(h,2,1)
    add_edge!(h,3,5); add_edge!(h,5,6); add_edge!(h,6,4)
    for g in testdigraphs(h)
      scc = @inferred(strongly_connected_components(g))
      @test length(scc) == 1 && sort(scc[1]) == [1:6;]
    end
    # tests from Graphs.jl
    h = DiGraph(4)
    add_edge!(h,1,2); add_edge!(h,2,3); add_edge!(h,3,1); add_edge!(h,4,1)
    for g in testdigraphs(h)
      scc = @inferred(strongly_connected_components(g))
      @test length(scc) == 2 && sort(scc[1]) == [1:3;] && sort(scc[2]) == [4]
    end

    h = DiGraph(12)
    add_edge!(h,1,2); add_edge!(h,2,3); add_edge!(h,2,4); add_edge!(h,2,5);
    add_edge!(h,3,6); add_edge!(h,4,5); add_edge!(h,4,7); add_edge!(h,5,2);
    add_edge!(h,5,6); add_edge!(h,5,7); add_edge!(h,6,3); add_edge!(h,6,8);
    add_edge!(h,7,8); add_edge!(h,7,10); add_edge!(h,8,7); add_edge!(h,9,7);
    add_edge!(h,10,9); add_edge!(h,10,11); add_edge!(h,11,12); add_edge!(h,12,10)

    for g in testdigraphs(h)
      scc = @inferred(strongly_connected_components(g))
      @test length(scc) == 4
      @test sort(scc[1]) == [7,8,9,10,11,12]
      @test sort(scc[2]) == [3, 6]
      @test sort(scc[3]) == [2, 4, 5]
      @test scc[4] == [1]
    end
    # Test examples with self-loops from
    # Graph-Theoretic Analysis of Finite Markov Chains by J.P. Jarvis & D. R. Shier

    # figure 1 example
    fig1 = spzeros(5,5)
    fig1[[3,4,9,10,11,13,18,19,22,24]] = [.5,.4,.1,1.,1.,.2,.3,.2,1.,.3]
    fig1 = DiGraph(fig1)
    scc_fig1 = Vector[[2,5],[1,3,4]]

    # figure 2 example
    fig2 = spzeros(5,5)
    fig2[[3, 10, 11, 13, 14, 17, 18, 19, 22]] = 1
    fig2 = DiGraph(fig2)

    # figure 3 example
    fig3 = spzeros(8,8)
    fig3[[1,7,9,13,14,15,18,20,23,27,28,31,33,34,37,45,46,49,57,63,64]] = 1
    fig3 = DiGraph(fig3)
    scc_fig3 = Vector[[3,4],[2,5,6],[8],[1,7]]
    fig3_cond = DiGraph(4);
    add_edge!(fig3_cond,4,3); add_edge!(fig3_cond,2,1)
    add_edge!(fig3_cond,4,1); add_edge!(fig3_cond,4,2)

    # construct a n-number edge ring graph (period = n)
    n = 10
    n_ring_m = spdiagm(ones(n-1),1,n,n); n_ring_m[end,1] = 1
    n_ring = DiGraph(n_ring_m)

    n_ring_shortcut = copy(n_ring); add_edge!(n_ring_shortcut,1,4)

    # figure 8 example
    fig8 = spzeros(6,6)
    fig8[[2,10,13,21,24,27,35]] = 1
    fig8 = DiGraph(fig8)

    @test Set(@inferred(strongly_connected_components(fig1))) == Set(scc_fig1)
    @test Set(@inferred(strongly_connected_components(fig3))) == Set(scc_fig3)

    @test @inferred(period(n_ring)) == n
    @test @inferred(period(n_ring_shortcut)) == 2

    @test @inferred(condensation(fig3)) == fig3_cond

    @test @inferred(attracting_components(fig1)) == Vector[[2,5]]
    @test @inferred(attracting_components(fig3)) == Vector[[3,4],[8]]

    g10 = StarGraph(10)
    for g in testgraphs(g10)
      @test @inferred(neighborhood(g, 1 , 0)) == [1]
      @test length(@inferred(neighborhood(g, 1, 1))) == 10
      @test length(@inferred(neighborhood(g, 2, 1))) == 2
      @test length(@inferred(neighborhood(g, 1, 2))) == 10
      @test length(@inferred(neighborhood(g, 2, 2))) == 10
    end
    g10 = StarDiGraph(10)
    for g in testdigraphs(g10)
      @test @inferred(neighborhood(g10, 1 , 0, dir=:out)) == [1]
      @test length(@inferred(neighborhood(g, 1, 1, dir=:out))) == 10
      @test length(@inferred(neighborhood(g, 2, 1, dir=:out))) == 1
      @test length(@inferred(neighborhood(g, 1, 2, dir=:out))) == 10
      @test length(@inferred(neighborhood(g, 2, 2, dir=:out))) == 1
      @test @inferred(neighborhood(g, 1 , 0, dir=:in)) == [1]
      @test length(@inferred(neighborhood(g, 1, 1, dir=:in))) == 1
      @test length(@inferred(neighborhood(g, 2, 1, dir=:in))) == 2
      @test length(@inferred(neighborhood(g, 1, 2, dir=:in))) == 1
      @test length(@inferred(neighborhood(g, 2, 2, dir=:in))) == 2
    end
    @test @inferred(!isgraphical([1,1,1]))
    @test @inferred(isgraphical([2,2,2]))
    @test @inferred(isgraphical(fill(3,10)))
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
# This file was generated, do not modify it. # hide
mutable struct CompositeModel2 <: DeterministicNetwork
    std_model::Standardizer
    box_model::UnivariateBoxCoxTransformer
    ridge_model::RidgeRegressor
end

function MLJ.fit(m::CompositeModel2, verbosity::Int, X, y)
    Xs = source(X)
    ys = source(y, kind=:target)
    W = transform(machine(m.std_model, Xs), Xs)
    box = machine(m.box_model, ys)
    z = transform(box, ys)
    ẑ = predict(machine(m.ridge_model, W, z), W)
    ŷ = inverse_transform(box, ẑ)
    fit!(ŷ, verbosity=0)
    return fitresults(ŷ)
end

mdl = CompositeModel2(Standardizer(), UnivariateBoxCoxTransformer(),
                      RidgeRegressor(lambda=0.1))
cm = machine(mdl, X, y)
res = evaluate!(cm, resampling=Holdout(fraction_train=0.8), measure=rms)
round(res.measurement[1], sigdigits=3)


text:
gamma(x::Float64) = nan_dom_err(ccall((:tgamma,libm),  Float64, (Float64,), x), x)
gamma(x::Float32) = nan_dom_err(ccall((:tgammaf,libm),  Float32, (Float32,), x), x)
gamma(x::Real) = gamma(float(x))
@vectorize_1arg Number gamma

function lgamma_r(x::Float64)
    signp = Array(Int32, 1)
    y = ccall((:lgamma_r,libm),  Float64, (Float64, Ptr{Int32}), x, signp)
    return y, signp[1]
end
function lgamma_r(x::Float32)
    signp = Array(Int32, 1)
    y = ccall((:lgammaf_r,libm),  Float32, (Float32, Ptr{Int32}), x, signp)
    return y, signp[1]
end
lgamma_r(x::Real) = lgamma_r(float(x))

lfact(x::Real) = (x<=1 ? zero(float(x)) : lgamma(x+one(x)))
@vectorize_1arg Number lfact

const clg_coeff = [76.18009172947146,
                   -86.50532032941677,
                   24.01409824083091,
                   -1.231739572450155,
                   0.1208650973866179e-2,
                   -0.5395239384953e-5]

function clgamma_lanczos(z)
    const sqrt2pi = 2.5066282746310005

    y = x = z
    temp = x + 5.5
    zz = log(temp)
    zz = zz * (x+0.5)
    temp -= zz
    ser = complex(1.000000000190015, 0)
    for j=1:6
        y += 1.0
        zz = clg_coeff[j]/y
        ser += zz
    end
    zz = sqrt2pi*ser / x
    return log(zz) - temp
end

function lgamma(z::Complex)
    if real(z) <= 0.5
        a = clgamma_lanczos(1-z)
        b = log(sinpi(z))
        const logpi = 1.14472988584940017
        z = logpi - b - a
    else
        z = clgamma_lanczos(z)
    end
    complex(real(z), angle_restrict_symm(imag(z)))
end

gamma(z::Complex) = exp(lgamma(z))

# Bernoulli numbers B_{2k}, using tabulated numerators and denominators from
# the online encyclopedia of integer sequences.  (They actually have data
# up to k=249, but we stop here at k=20.)  Used for generating the polygamma
# (Stirling series) coefficients below.
#   const A000367 = map(BigInt, split("1,1,-1,1,-1,5,-691,7,-3617,43867,-174611,854513,-236364091,8553103,-23749461029,8615841276005,-7709321041217,2577687858367,-26315271553053477373,2929993913841559,-261082718496449122051", ","))
#   const A002445 = [1,6,30,42,30,66,2730,6,510,798,330,138,2730,6,870,14322,510,6,1919190,6,13530]
#   const bernoulli = A000367 .// A002445 # even-index Bernoulli numbers

function digamma(z::Union(Float64,Complex{Float64}))
    # Based on eq. (12), without looking at the accompanying source
    # code, of: K. S. Kölbig, "Programs for computing the logarithm of
    # the gamma function, and the digamma function, for complex
    # argument," Computer Phys. Commun.  vol. 4, pp. 221–226 (1972).
    x = real(z)
    if x <= 0 # reflection formula
        ψ = -π * cot(π*z)
        z = 1 - z
        x = real(z)
    else
        ψ = zero(z)
    end
    if x < 7
        # shift using recurrence formula
        n = 7 - ifloor(x)
        for ν = 1:n-1
            ψ -= inv(z + ν)
        end
        ψ -= inv(z)
        z += n
    end
    t = inv(z)
    ψ += log(z) - 0.5*t
    t *= t # 1/z^2
    # the coefficients here are float64(bernoulli[2:9] .// (2*(1:8)))
    ψ -= t * @evalpoly(t,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686)
end

function trigamma(z::Union(Float64,Complex{Float64}))
    # via the derivative of the Kölbig digamma formulation
    x = real(z)
    if x <= 0 # reflection formula
        return (π * csc(π*z))^2 - trigamma(1 - z)
    end
    ψ = zero(z)
    if x < 8
        # shift using recurrence formula
        n = 8 - ifloor(x)
        ψ += inv(z)^2
        for ν = 1:n-1
            ψ += inv(z + ν)^2
        end
        z += n
    end
    t = inv(z)
    w = t * t # 1/z^2
    ψ += t + 0.5*w
    # the coefficients here are float64(bernoulli[2:9])
    ψ += t*w * @evalpoly(w,0.16666666666666666,-0.03333333333333333,0.023809523809523808,-0.03333333333333333,0.07575757575757576,-0.2531135531135531,1.1666666666666667,-7.092156862745098)
end

signflip(m::Number, z) = (-1+0im)^m * z
signflip(m::Integer, z) = iseven(m) ? z : -z

# (-1)^m d^m/dz^m cot(z) = p_m(cot z), where p_m is a polynomial
# that satisfies the recurrence p_{m+1}(x) = p_m′(x) * (1 + x^2).
# Note that p_m(x) has only even powers of x if m is odd, and
# only odd powers of x if m is even.  Therefore, we write p_m(x)
# as p_m(x) = q_m(x^2) m! for m odd and p_m(x) = x q_m(x^2) m! if m is even.
# Hence the polynomials q_m(y) satisfy the recurrence:
#     m odd:  q_{m+1}(y) = 2 q_m′(y) * (1+y) / (m+1)
#    m even:  q_{m+1}(y) = [q_m(y) + 2 y q_m′(y)] * (1+y) / (m+1)
# This function computes the coefficients of the polynomial q_m(y),
# returning an array of the coefficients of 1, y, y^2, ...,
function cotderiv_q(m::Int)
    m < 0 && throw(ArgumentError("$m < 0 not allowed"))
    m == 0 && return [1.0]
    m == 1 && return [1.0, 1.0]
    q₋ = cotderiv_q(m-1)
    d = length(q₋) - 1 # degree of q₋
    if isodd(m-1)
        q = Array(Float64, length(q₋))
        q[end] = d * q₋[end] * 2/m
        for i = 1:length(q)-1
            q[i] = ((i-1)*q₋[i] + i*q₋[i+1]) * 2/m
        end
    else # iseven(m-1)
        q = Array(Float64, length(q₋) + 1)
        q[1] = q₋[1] / m
        q[end] = (1 + 2d) * q₋[end] / m
        for i = 2:length(q)-1
            q[i] = ((1 + 2(i-1))*q₋[i] + (1 + 2(i-2))*q₋[i-1]) / m
        end
    end
    return q
end

# precompute a table of cot derivative polynomials
const cotderiv_Q = [cotderiv_q(m) for m in 1:100]

# Evaluate (-1)^m d^m/dz^m [π cot(πz)] / m!.  If m is small, we
# use the explicit derivative formula (a polynomial in cot(πz));
# if m is large, we use the derivative of Euler's harmonic series:
#          π cot(πz) = ∑ inv(z + n)
function cotderiv(m::Integer, z)
    isinf(imag(z)) && return zero(z)
    if m <= 0
        m == 0 && return π * cot(π*z)
        throw(DomainError())
    end
    if m <= length(cotderiv_Q)
        q = cotderiv_Q[m]
        x = cot(π*z)
        y = x*x
        s = q[1] + q[2] * y
        t = y
        # evaluate q(y) using Horner (TODO: Knuth for complex z?)
        for i = 3:length(q)
            t *= y
            s += q[i] * t
        end
        return π^(m+1) * (isodd(m) ? s : x*s)
    else # m is large, series derivative should converge quickly
        p = m+1
        z -= round(real(z))
        s = inv(z^p)
        n = 1
        sₒ = zero(s)
        while s != sₒ
            sₒ = s
            a = (z+n)^p
            b = (z-n)^p
            s += (a + b) / (a * b)
            n += 1
        end
        return s
    end
end

# Helper macro for polygamma(m, z):
#   Evaluate p[1]*c[1] + x*p[2]*c[2] + x^2*p[3]*c[3] + ...
#   where c[1] = m + 1
#         c[k] = c[k-1] * (2k+m-1)*(2k+m-2) / ((2k-1)*(2k-2)) = c[k-1] * d[k]
#         i.e. d[k] = c[k]/c[k-1] = (2k+m-1)*(2k+m-2) / ((2k-1)*(2k-2))
#   by a modified version of Horner's rule:
#      c[1] * (p[1] + d[2]*x * (p[2] + d[3]*x * (p[3] + ...))).
# The entries of p must be literal constants and there must be > 1 of them.
macro pg_horner(x, m, p...)
    k = length(p)
    me = esc(m)
    xe = esc(x)
    ex = :(($me + $(2k-1)) * ($me + $(2k-2)) * $(p[end]/((2k-1)*(2k-2))))
    for k = length(p)-1:-1:2
        cdiv = 1 / ((2k-1)*(2k-2))
        ex = :(($cdiv * ($me + $(2k-1)) * ($me + $(2k-2))) *
               ($(p[k]) + $xe * $ex))
    end
    :(($me + 1) * ($(p[1]) + $xe * $ex))
end

# compute inv(oftype(x, y)) efficiently, choosing the correct branch cut
inv_oftype(x::Complex, y::Complex) = oftype(x, inv(y))
function inv_oftype(x::Complex, y::Real)
    yi = inv(y) # using real arithmetic for efficiency
    oftype(x, Complex(yi, -zero(yi))) # get correct sign of zero!
end
inv_oftype(x::Real, y::Real) = oftype(x, inv(y))

# Hurwitz zeta function, which is related to polygamma
# (at least for integer m > 0 and real(z) > 0) by:
#    polygamma(m, z) = (-1)^(m+1) * gamma(m+1) * zeta(m+1, z).
# Our algorithm for the polygamma is just the m-th derivative
# of our digamma approximation, and this derivative process yields
# a function of the form
#          (-1)^(m) * gamma(m+1) * (something)
# So identifying the (something) with the -zeta function, we get
# the zeta function for free and might as well export it, especially
# since this is a common generalization of the Riemann zeta function
# (which Julia already exports).
function zeta(s::Union(Int,Float64,Complex{Float64}),
              z::Union(Float64,Complex{Float64}))
    ζ = zero(promote_type(typeof(s), typeof(z)))

    # like sqrt, require complex inputs to get complex outputs
    !isa(s,Integer) && isa(ζ, Real) && z < 0 && throw(DomainError())

    z == 1 && return oftype(ζ, zeta(s))
    s == 2 && return oftype(ζ, trigamma(z))

    x = real(z)

    # annoying s = Inf or NaN cases:
    if !isfinite(s)
        (isnan(s) || isnan(z)) && return (s*z)^2 # type-stable NaN+Nan*im
        if real(s) == Inf
            z == 1 && return one(ζ)
            if x > 1 || (x >= 0.5 ? abs(z) > 1 : abs(z - round(x)) > 1)
                return zero(ζ) # distance to poles is > 1
            end
            x > 0 && imag(z) == 0 && imag(s) == 0 && return oftype(ζ, Inf)
        end
        throw(DomainError()) # nothing clever to return
    end

    # We need a different algorithm for the real(s) < 1 domain
    real(s) < 1 && throw(ArgumentError("order $s < 1 is not implemented (issue #7228)"))

    m = s - 1

    # Algorithm is just the m-th derivative of digamma formula above,
    # with a modified cutoff of the final asymptotic expansion.

    # Note: we multiply by -(-1)^m m! in polygamma below, so this factor is
    #       pulled out of all of our derivatives.

    isnan(x) && return oftype(ζ, imag(z)==0 && isa(s,Int) ? x : Complex(x,x))

    cutoff = 7 + real(m) + imag(m) # TODO: this cutoff is too conservative?
    if x < cutoff
        # shift using recurrence formula
        xf = floor(x)
        if x <= 0 && xf == z
            if isa(s, Int)
                iseven(s) && return oftype(ζ, Inf)
                x == 0 && return oftype(ζ, inv(x))
            end
            throw(DomainError()) # or return NaN?
        end
        nx = int(xf)
        n = iceil(cutoff - nx)
        ζ += inv_oftype(ζ, z)^s
        for ν = -nx:-1:1
            ζₒ= ζ
            ζ += inv_oftype(ζ, z + ν)^s
            ζ == ζₒ && break # prevent long loop for large -x > 0
                             # FIXME: still slow for small m, large Im(z)
        end
        for ν = max(1,1-nx):n-1
            ζₒ= ζ
            ζ += inv_oftype(ζ, z + ν)^s
            ζ == ζₒ && break # prevent long loop for large m
        end
        z += n
    end

    t = inv(z)
    w = isa(t, Real) ? conj(oftype(ζ, t))^m : oftype(ζ, t)^m
    ζ += w * (inv(m) + 0.5*t)

    t *= t # 1/z^2
    ζ += w*t * @pg_horner(t,m,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686,3.0539543302701198)

    return ζ
end

function polygamma(m::Integer, z::Union(Float64,Complex{Float64}))

    m == 0 && return digamma(z)
    m == 1 && return trigamma(z)

    # In principle, we could support non-integer m here, but the
    # extension to complex m seems to be non-unique, the obvious
    # extension (just plugging in a complex m below) does not seem to
    # be the one used by Mathematica or Maple, and sources do not
    # agree on what the "right" extension is (e.g. Mathematica & Maple
    # disagree).   So, at least for now, we require integer m.

    # real(m) < 0 is valid, but I don't think our asymptotic expansion
    # here works in this case.  m < 0 polygamma is called a
    # "negapolygamma" function in the literature, and there are
    # multiple possible definitions arising from different integration
    # constants. We throw a DomainError() since the definition is unclear.
    real(m) < 0 && throw(DomainError())

    s = m+1
    if real(z) <= 0 # reflection formula
        (zeta(s, 1-z) + signflip(m, cotderiv(m,z))) * (-gamma(s))
    else
        signflip(m, zeta(s,z) * (-gamma(s)))
    end
end

# If we really cared about single precision, we could make a faster
# Float32 version by truncating the Stirling series at a smaller cutoff.
for (f,T) in ((float32,Float32), (float16,Float16))
    @eval begin
        zeta(s::Integer, z::Union($T,Complex{$T})) =
            $f(zeta(int(s), float64(z)))
        zeta(s::Union(Float64,Complex{Float64}), z::Union($T,Complex{$T})) =
            zeta(s, float64(z))
        zeta(s::Number, z::Union($T,Complex{$T})) =
            $f(zeta(float64(s), float64(z)))
        polygamma(m::Integer, z::Union($T,Complex{$T})) =
            $f(polygamma(int(m), float64(z)))
        digamma(z::Union($T,Complex{$T})) =
            $f(digamma(float64(z)))
        trigamma(z::Union($T,Complex{$T})) =
            $f(trigamma(float64(z)))
    end
end

zeta(s::Integer, z::Number) = zeta(int(s), float64(z))
zeta(s::Number, z::Number) = zeta(float64(s), float64(z))
for f in (:digamma, :trigamma)
    @eval begin
        $f(z::Number) = $f(float64(z))
        @vectorize_1arg Number $f
    end
end
polygamma(m::Integer, z::Number) = polygamma(m, float64(z))
@vectorize_2arg Number polygamma
@vectorize_2arg Number zeta

# Inverse digamma function:
# Implementation of fixed point algorithm described in
#  "Estimating a Dirichlet distribution" by Thomas P. Minka, 2000
function invdigamma(y::Float64)
    # Closed form initial estimates
    if y >= -2.22
        x_old = exp(y) + 0.5
        x_new = x_old
    else
        x_old = -1.0 / (y - digamma(1.0))
        x_new = x_old
    end

    # Fixed point algorithm
    delta = Inf
    iteration = 0
    while delta > 1e-12 && iteration < 25
        iteration += 1
        x_new = x_old - (digamma(x_old) - y) / trigamma(x_old)
        delta = abs(x_new - x_old)
        x_old = x_new
    end

    return x_new
end
invdigamma(x::Float32) = float32(invdigamma(float64(x)))
invdigamma(x::Real) = invdigamma(float64(x))
@vectorize_1arg Real invdigamma

function beta(x::Number, w::Number)
    yx, sx = lgamma_r(x)
    yw, sw = lgamma_r(w)
    yxw, sxw = lgamma_r(x+w)
    return copysign(exp(yx + yw - yxw), sx*sw*sxw)
end
lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)
@vectorize_2arg Number beta
@vectorize_2arg Number lbeta

# Riemann zeta function; algorithm is based on specializing the Hurwitz
# zeta function above for z==1.
function zeta(s::Union(Float64,Complex{Float64}))
    # blows up to ±Inf, but get correct sign of imaginary zero
    s == 1 && return NaN + zero(s) * imag(s)

    if !isfinite(s) # annoying NaN and Inf cases
        isnan(s) && return imag(s) == 0 ? s : s*s
        if isfinite(imag(s))
            real(s) > 0 && return 1.0 - zero(s)*imag(s)
            imag(s) == 0 && return NaN + zero(s)
        end
        return NaN*zero(s) # NaN + NaN*im
    elseif real(s) < 0.5
        if abs(real(s)) + abs(imag(s)) < 1e-3 # Taylor series for small |s|
            return @evalpoly(s, -0.5,
                             -0.918938533204672741780329736405617639861,
                             -1.0031782279542924256050500133649802190,
                             -1.00078519447704240796017680222772921424,
                             -0.9998792995005711649578008136558752359121)
        end
        return zeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * (2π)^s / π
    end

    m = s - 1

    # shift using recurrence formula:
    #   n is a semi-empirical cutoff for the Stirling series, based
    #   on the error term ~ (|m|/n)^18 / n^real(m)
    n = iceil(6 + 0.7*abs(imag(s-1))^inv(1 + real(m)*0.05))
    ζ = one(s)
    for ν = 2:n
        ζₒ= ζ
        ζ += inv(ν)^s
        ζ == ζₒ && break # prevent long loop for large m
    end
    z = 1 + n
    t = inv(z)
    w = t^m
    ζ += w * (inv(m) + 0.5*t)

    t *= t # 1/z^2
    ζ += w*t * @pg_horner(t,m,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686,3.0539543302701198)

    return ζ
end

zeta(x::Integer) = zeta(float64(x))
zeta(x::Real)    = oftype(float(x),zeta(float64(x)))
zeta(z::Complex) = oftype(float(z),zeta(complex128(z)))
@vectorize_1arg Number zeta

function eta(z::Union(Float64,Complex{Float64}))
    δz = 1 - z
    if abs(real(δz)) + abs(imag(δz)) < 7e-3 # Taylor expand around z==1
        return 0.6931471805599453094172321214581765 *
               @evalpoly(δz,
                         1.0,
                         -0.23064207462156020589789602935331414700440,
                         -0.047156357547388879740146103148112380421254,
                         -0.002263576552598880778433550956278702759143568,
                         0.001081837223249910136105931217561387128141157)
    else
        return -zeta(z) * expm1(0.6931471805599453094172321214581765*δz)
    end
end
eta(x::Integer) = eta(float64(x))
eta(x::Real)    = oftype(float(x),eta(float64(x)))
eta(z::Complex) = oftype(float(z),eta(complex128(z)))
@vectorize_1arg Number eta

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
"""
    get_proportion(for_prop_arr::AbstractVector{Float64})::AbstractVector{Float64}
Compute the proportional values for each element in the 1-Dim array.
returns a AbstractVector{Float64}
# Example
```
julia> get_proportion([1.0,2.0,3.0,4.0,5.0,34.0000034,2423,5656.98988])
=[0.00012301651427098885, 0.0002460330285419777, 0.00036904954281296654, 0.0004920660570839554, 0.0006150825713549442, 0.004182561903469769, 0.298069014078606, 0.6959031763038594]
```
"""
function get_proportion(for_prop_arr::AbstractVector{Float64})::AbstractVector{Float64}
    return for_prop_arr ./ sum(for_prop_arr)
end

"""
    get_proportion_round(for_prop_arr::AbstractVector{Float64};round_digits::Int64)::AbstractVector{Float64}
Compute the proportional values for each element in the array and round to a specified decimal digits.
If the sum of all the elements is not 1 , then the adjustment is done in the 1st largest element
returns a AbstractVector{Float64}
# Example
```
julia> get_proportion_round([1.33,1.33,1.23,1.1111],round_digits = 2)
=[0.26, 0.27, 0.25, 0.22]
```
"""
function get_proportion_round(for_prop_arr::AbstractVector{Float64};round_digits::Int64 = 2)::AbstractVector{Float64}
    temparr=  round.(for_prop_arr ./ sum(for_prop_arr), digits=round_digits, base = 10)
    diff = 1 - sum(temparr) #delta naming
    if diff != 0.0
        max = argmax(temparr)
        temparr[max] =  round(temparr[max] + diff, digits=round_digits, base = 10)
    end
    return temparr
end

"""
    get_proportion_round_add1251(for_prop_arr::AbstractVector{Float64};round_digits::Int64)::AbstractVector{Float64}
Compute the proportional values for each element in the array and round to a specified decimal digits.
If the sum of all the elements is not 1 , then add 1251 to the array and try looking for proportional correctness
after 2 itterations , the adjustment is done in the 1st largest element
returns a AbstractVector{Float64}
# Example
```
julia> get_proportion_round_add1251([1.33,1.33,1.23,1.1111],round_digits = 2)
=[0.26, 0.27, 0.25, 0.22]
```
"""
function get_proportion_round_add1251(for_prop_arr::AbstractVector{Float64};round_digits::Int64 = 2)::AbstractVector{Float64}
    temparr=  round.(for_prop_arr ./ sum(for_prop_arr), digits=round_digits, base = 10)
    diff = 1 - sum(temparr)
    if diff != 0.0
        check = 1
        while check < 2
            sqarr= (for_prop_arr) .+1251
            temparr=  round.(sqarr ./ sum(sqarr), digits=round_digits, base = 10)
            diff = 1 - sum(temparr)
            if diff != 0.00
                check = check + 1
                if check == 2
                    max = argmax(temparr)
                    temparr[max] = round(temparr[max]+diff,digits=round_digits, base = 10)
                end
            else
                break
            end
        end
    end
    return temparr
end

max_stars_repo_name: hillelawaskar/Proportions

max_issues_repo_name: hillelawaskar/Proportions