text:
export polynomial, polynomial!, polynomialtype, terms, nterms, coefficients, monomials
export coefficienttype, monomialtype
export mindegree, maxdegree, extdegree, effective_variables
export leadingterm, leadingcoefficient, leadingmonomial
export removeleadingterm, removemonomials, monic
export mapcoefficients, mapcoefficients!, mapcoefficients_to!

LinearAlgebra.norm(p::AbstractPolynomialLike, r::Int=2) = LinearAlgebra.norm(coefficients(p), r)

changecoefficienttype(::Type{TT}, ::Type{T}) where {TT<:AbstractTermLike, T} = termtype(TT, T)
changecoefficienttype(::Type{PT}, ::Type{T}) where {PT<:AbstractPolynomial, T} = polynomialtype(PT, T)

changecoefficienttype(p::PT, ::Type{T}) where {PT<:APL, T} = convert(changecoefficienttype(PT, T), p)

abstract type ListState end
abstract type UnsortedState <: ListState end
struct MessyState <: UnsortedState end
# No duplicates or zeros
struct UniqState <: UnsortedState end
sortstate(::MessyState) = SortedState()
sortstate(::UniqState) = SortedUniqState()
struct SortedState <: ListState end
struct SortedUniqState <: ListState end

"""
    polynomial(p::AbstractPolynomialLike)

Converts `p` to a value with polynomial type.

    polynomial(p::AbstractPolynomialLike, ::Type{T}) where T

Converts `p` to a value with polynomial type with coefficient type `T`.

    polynomial(a::AbstractVector, mv::AbstractVector{<:AbstractMonomialLike})

Creates a polynomial equal to `dot(a, mv)`.

    polynomial(terms::AbstractVector{<:AbstractTerm}, s::ListState=MessyState())

Creates a polynomial equal to `sum(terms)` where `terms` are guaranteed to be in state `s`.

    polynomial(f::Function, mv::AbstractVector{<:AbstractMonomialLike})

Creates a polynomial equal to `sum(f(i) * mv[i] for i in 1:length(mv))`.

### Examples

Calling `polynomial([2, 4, 1], [x, x^2*y, x*y])` should return ``4x^2y + xy + 2x``.
"""
polynomial(p::AbstractPolynomial) = p
polynomial(p::APL{T}, ::Type{T}) where T = polynomial(terms(p))
polynomial(p::APL{T}) where T = polynomial(p, T)
function polynomial(Q::AbstractMatrix, mv::AbstractVector)
    LinearAlgebra.dot(mv, Q * mv)
end
function polynomial(Q::AbstractMatrix, mv::AbstractVector, ::Type{T}) where T
    polynomial(polynomial(Q, mv), T)
end

polynomial(f::Function, mv::AbstractVector{<:AbstractMonomialLike}) = polynomial!([term(f(i), mv[i]) for i in 1:length(mv)])

function polynomial(a::AbstractVector, x::AbstractVector, s::ListState=MessyState())
    # If `x` is e.g. `[v, 1]` then it will contains terms that are convertible to monomials.
    return polynomial([term(α, convert(monomialtype(m), m)) for (α, m) in zip(a, x)], s)
end

polynomial(ts::AbstractVector, s::ListState=MessyState()) = sum(ts)
polynomial!(ts::AbstractVector, s::ListState=MessyState()) = sum(ts)

polynomial!(ts::AbstractVector{<:AbstractTerm}, s::SortedUniqState) = polynomial(coefficient.(ts), monomial.(ts), s)

function polynomial!(ts::AbstractVector{<:AbstractTerm}, s::SortedState)
    polynomial!(uniqterms!(ts), SortedUniqState())
end
function polynomial!(ts::AbstractVector{<:AbstractTerm}, s::UnsortedState=MessyState())
    polynomial!(sort!(ts, lt=(>)), sortstate(s))
end

_collect(v::Vector) = v
_collect(v::AbstractVector) = collect(v)
polynomial(ts::AbstractVector{<:AbstractTerm}, args::Vararg{ListState, N}) where {N} = polynomial!(MA.mutable_copy(_collect(ts)), args...)

"""
    polynomialtype(p::AbstractPolynomialLike)

Returns the type that `p` would have if it was converted into a polynomial.

    polynomialtype(::Type{PT}) where PT<:AbstractPolynomialLike

Returns the same as `polynomialtype(::PT)`.

    polynomialtype(p::AbstractPolynomialLike, ::Type{T}) where T

Returns the type that `p` would have if it was converted into a polynomial of coefficient type `T`.

    polynomialtype(::Type{PT}, ::Type{T}) where {PT<:AbstractPolynomialLike, T}

Returns the same as `polynomialtype(::PT, ::Type{T})`.
"""
polynomialtype(::Union{P, Type{P}}) where P <: APL = Base.promote_op(polynomial, P)
polynomialtype(::Union{P, Type{P}}) where P <: AbstractPolynomial = P
polynomialtype(::Union{M, Type{M}}) where M<:AbstractMonomialLike = polynomialtype(termtype(M))
polynomialtype(::Union{M, Type{M}}, ::Type{T}) where {M<:AbstractMonomialLike, T} = polynomialtype(termtype(M, T))
polynomialtype(::Union{P, Type{P}}, ::Type{T}) where {P <: APL, T} = polynomialtype(polynomialtype(P), T)
polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}) where PT <: APL = polynomialtype(PT)
polynomialtype(::Union{AbstractVector{PT}, Type{<:AbstractVector{PT}}}, ::Type{T}) where {PT <: APL, T} = polynomialtype(PT, T)

function uniqterms!(ts::AbstractVector{<: AbstractTerm})
    i = firstindex(ts)
    for j in Iterators.drop(eachindex(ts), 1)
        if !iszero(ts[j])
            if monomial(ts[i]) == monomial(ts[j])
                ts[i] = term(MA.add!!(coefficient(ts[i]), coefficient(ts[j])), monomial(ts[i]))
            else
                if !iszero(ts[i])
                    i += 1
                end
                ts[i] = MA.copy_if_mutable(ts[j])
            end
        end
    end
    if i < length(ts)
        if iszero(ts[i])
            i -= 1
        end
        resize!(ts, i)
    end
    ts
end

"""
    terms(p::AbstractPolynomialLike)

Returns an iterator over the nonzero terms of the polynomial `p` sorted in the decreasing monomial order.

### Examples

Calling `terms` on ``4x^2y + xy + 2x`` should return an iterator of ``[4x^2y, xy, 2x]``.
"""
terms(t::AbstractTermLike) = iszero(t) ? termtype(t)[] : [term(t)]
terms(p::AbstractPolynomialLike) = terms(polynomial(p))

"""
    nterms(p::AbstractPolynomialLike)

Returns the number of nonzero terms in `p`, i.e. `length(terms(p))`.

### Examples

Calling `nterms` on ``4x^2y + xy + 2x`` should return 3.
"""
function nterms end

nterms(t::AbstractTermLike) = iszero(t) ? 0 : 1
nterms(p::AbstractPolynomialLike) = length(terms(p))

"""
    coefficients(p::AbstractPolynomialLike)

Returns an iterator over the coefficients of `p` of the nonzero terms of the polynomial sorted in the decreasing monomial order.

    coefficients(p::AbstractPolynomialLike, X::AbstractVector)

Returns an iterator over the coefficients of the monomials of `X` in `p` where `X` is a monomial vector not necessarily sorted but with no duplicate entry.

### Examples

Calling `coefficients` on ``4x^2y + xy + 2x`` should return an iterator of ``[4, 1, 2]``.
Calling `coefficients(4x^2*y + x*y + 2x + 3, [x, 1, x*y, y])` should return an iterator of ``[2, 3, 1, 0]``.
"""
coefficients(p::APL) = coefficient.(terms(p))
function coefficients(p::APL{T}, X::AbstractVector) where T
    σ, mv = sortmonovec(X)
    @assert length(mv) == length(X) # no duplicate in X
    c = zeros(T, length(mv))
    i = 1
    for t in terms(p)
        m = monomial(t)
        while i <= length(mv) && mv[i] > m
            c[σ[i]] = zero(T)
            i += 1
        end
        if i <= length(mv) && mv[i] == m
            c[σ[i]] = coefficient(t)
            i += 1
        end
    end
    c
end

"""
    monomials(p::AbstractPolynomialLike)

Returns an iterator over the monomials of `p` of the nonzero terms of the polynomial sorted in the decreasing order.

    monomials(vars::Tuple, degs::AbstractVector{Int}, filter::Function = m -> true)

Builds the vector of all the monovec `m` with variables `vars` such that the degree `degree(m)` is in `degs` and `filter(m)` is `true`.

### Examples

Calling `monomials` on ``4x^2y + xy + 2x`` should return an iterator of ``[x^2y, xy, x]``.

Calling `monomials((x, y), [1, 3], m -> degree(m, y) != 1)` should return `[x^3, x*y^2, y^3, x]` where `x^2*y` and `y` have been excluded by the filter.
"""
monomials(p::APL) = monovec(monomial.(terms(p)))

function isconstant(p::APL)
    n = nterms(p)
    return iszero(n) || (isone(n) && isconstant(first(terms(p))))
end

#$(SIGNATURES)
"""
    mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Returns the minimal total degree of the monomials of `p`, i.e. `minimum(degree, terms(p))`.

    mindegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)

Returns the minimal degree of the monomials of `p` in the variable `v`, i.e. `minimum(degree.(terms(p), v))`.

### Examples
Calling `mindegree` on on ``4x^2y + xy + 2x`` should return 1, `mindegree(4x^2y + xy + 2x, x)` should return 1 and  `mindegree(4x^2y + xy + 2x, y)` should return 0.
"""
function mindegree(X::AbstractVector{<:AbstractTermLike}, args...)
    isempty(X) ? 0 : minimum(t -> degree(t, args...), X)
end
function mindegree(p::AbstractPolynomialLike, args...)
    mindegree(terms(p), args...)
end

#$(SIGNATURES)
"""
    maxdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Returns the maximal total degree of the monomials of `p`, i.e. `maximum(degree, terms(p))`.

    maxdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)

Returns the maximal degree of the monomials of `p` in the variable `v`, i.e. `maximum(degree.(terms(p), v))`.

### Examples
Calling `maxdegree` on on ``4x^2y + xy + 2x`` should return 3, `maxdegree(4x^2y + xy + 2x, x)` should return 2 and  `maxdegree(4x^2y + xy + 2x, y)` should return 1.
"""
function maxdegree(X::AbstractVector{<:AbstractTermLike}, args::Vararg{Any,N}) where {N}
    return mapreduce(t -> degree(t, args...), max, X, init=0)
end
function maxdegree(p::AbstractPolynomialLike, args::Vararg{Any,N}) where {N}
    maxdegree(terms(p), args...)
end

#$(SIGNATURES)
"""
    extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Returns the extremal total degrees of the monomials of `p`, i.e. `(mindegree(p), maxdegree(p))`.

    extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)

Returns the extremal degrees of the monomials of `p` in the variable `v`, i.e. `(mindegree(p, v), maxdegree(p, v))`.

### Examples
Calling `extdegree` on on ``4x^2y + xy + 2x`` should return `(1, 3)`, `extdegree(4x^2y + xy + 2x, x)` should return `(1, 2)` and  `maxdegree(4x^2y + xy + 2x, y)` should return `(0, 1)`.
"""
function extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, args...)
    (mindegree(p, args...), maxdegree(p, args...))
end

"""
    effective_variables(p::AbstractPolynomialLike)

Return a vector of `eltype` `variable_union_type(p)` (see [`variable_union_type`](@ref)),
containing all the variables that has nonzero degree in at least one term.
That is, return all the variables `v` such that `maxdegree(p, v)` is not zero.
The returned vector is sorted in decreasing order.
"""
function effective_variables(p::AbstractPolynomialLike)
    VT = variable_union_type(p)
    return VT[v for v in variables(p) if !iszero(maxdegree(p, v))]
end

"""
    leadingterm(p::AbstractPolynomialLike)

Returns the coefficient of the leading term, i.e. `first(terms(p))`.

### Examples

Calling `leadingterm` on ``4x^2y + xy + 2x`` should return ``4x^2y``.
"""
function leadingterm(p::AbstractPolynomialLike)
    if iszero(p)
        zeroterm(p)
    else
        first(terms(p))
    end
end
leadingterm(t::AbstractTermLike) = term(t)

#$(SIGNATURES)
"""
    leadingcoefficient(p::AbstractPolynomialLike)

Returns the coefficient of the leading term of `p`, i.e. `coefficient(leadingterm(p))`.

### Examples

Calling `leadingcoefficient` on ``4x^2y + xy + 2x`` should return ``4`` and calling it on ``0`` should return ``0``.
"""
function leadingcoefficient(p::AbstractPolynomialLike)
    coefficient(leadingterm(p))
end

#$(SIGNATURES)
"""
    leadingmonomial(p::AbstractPolynomialLike)

Returns the monomial of the leading term of `p`, i.e. `monomial(leadingterm(p))` or `first(monomials(p))`.

### Examples

Calling `leadingmonomial` on ``4x^2y + xy + 2x`` should return ``x^2y``.
"""
function leadingmonomial(p::AbstractPolynomialLike)
    # first(monomials(p)) would be more efficient for DynamicPolynomials but
    # monomial(leadingterm(p)) is more efficient for TypedPolynomials and is better if p is a term
    monomial(leadingterm(p))
end

#$(SIGNATURES)
"""
    removeleadingterm(p::AbstractPolynomialLike)

Returns a polynomial with the leading term removed in the polynomial `p`.

### Examples

Calling `removeleadingterm` on ``4x^2y + xy + 2x`` should return ``xy + 2x``.
"""
function removeleadingterm(p::AbstractPolynomialLike)
    # Iterators.drop returns an Interators.Drop which is not an AbstractVector
    polynomial(terms(p)[2:end], SortedUniqState())
end

#$(SIGNATURES)
"""

Returns a polynomial with the terms having their monomial in the monomial vector `mv` removed in the polynomial `p`.

### Examples

Calling `removemonomials(4x^2*y + x*y + 2x, [x*y])` should return ``4x^2*y + 2x``.
"""
function removemonomials(p::AbstractPolynomialLike, mv::AbstractVector{MT}) where {MT <: AbstractMonomialLike}
    smv = monovec(mv) # Make sure it is sorted
    i = 1
    q = zero(p)
    for t in terms(p)
        m = monomial(t)
        while i <= length(smv) && smv[i] > m
            i += 1
        end
        if i > length(smv) || smv[i] != m
            q += t
        end
    end
    q
end

"""
    monic(p::AbstractPolynomialLike)

Returns `p / leadingcoefficient(p)` where the leading coefficient of the returned polynomials is made sure to be exactly one to avoid rounding error.
"""
function monic(p::APL)
    α = leadingcoefficient(p)
    polynomial!(_divtoone.(terms(p), α))
end
monic(m::AbstractMonomialLike) = m
monic(t::AbstractTermLike{T}) where T = term(one(T), monomial(t))

function _divtoone(t::AbstractTermLike{T}, α::S) where {T, S}
    U = Base.promote_op(/, T, S)
    β = coefficient(t)
    if β == α
        term(one(U), monomial(t))
    else
        term((β / α), monomial(t))
    end
end

# TODO deprecate
mapcoefficientsnz(f::Function, p::APL) = mapcoefficients(f, p, nonzero = true)
mapcoefficientsnz_to!(output::APL, f::Function, p::APL) = mapcoefficients_to!(output, f, p, nonzero = true)

"""
    mapcoefficients(f::Function, p::AbstractPolynomialLike, nonzero = false)

Returns a polynomial with the same monomials as `p` but each coefficient `α` is replaced by `f(α)`.
The function may return zero in which case the term is dropped.
If the function is known to never returns zero for a nonzero input, `nonzero`
can be set to `true` to get a small speedup.

See also [`mapcoefficients!`](@ref) and [`mapcoefficients_to!`](@ref).

### Examples

Calling `mapcoefficients(α -> mod(3α, 6), 2x*y + 3x + 1)` should return `3x + 3`.
"""
function mapcoefficients end
function mapcoefficients(f::Function, p::AbstractPolynomialLike; nonzero = false) # Not used by either TypedPolynomials or DynamicPolynomials but used by CustomPoly in tests. FIXME Remove in a breaking release
    # Invariant: p has only nonzero coefficient
    # therefore f(α) will be nonzero for every coefficient α of p
    # hence we can use Uniq
    polynomial!(mapcoefficients.(f, terms(p)), nonzero ? SortedUniqState() : SortedState())
end
function mapcoefficients(f::Function, t::AbstractTermLike; nonzero = false)
    return term(f(coefficient(t)), monomial(t))
end

"""
    mapcoefficients!(f::Function, p::AbstractPolynomialLike, nonzero = false)

Mutate `p` by replacing each coefficient `α` by `f(α)`.
The function may return zero in which case the term is dropped.
If the function is known to never returns zero for a nonzero input, `nonzero`
can be set to `true` to get a small speedup.
The function returns `p`, which is identically equal to the second argument.

See also [`mapcoefficients`](@ref) and [`mapcoefficients_to!`](@ref).

### Examples

Let `p = 2x*y + 3x + 1`, after `mapcoefficients!(α -> mod(3α, 6), p)`, `p` is
equal to `3x + 3`.
"""
function mapcoefficients! end

"""
    mapcoefficients_to!(output::AbstractPolynomialLike, f::Function, p::AbstractPolynomialLike, nonzero = false)

Mutate `output` by replacing each coefficient `α` of `p` by `f(α)`.
The function may return zero in which case the term is dropped.
If the function is known to never returns zero for a nonzero input, `nonzero`
can be set to `true` to get a small speedup.
The function returns `output`, which is identically equal to the first argument.

See also [`mapcoefficients!`](@ref) and [`mapcoefficients`](@ref).
"""
function mapcoefficients_to! end

"""
    deg_num_leading_terms(p::AbstractPolynomialLike, var)

Return `deg, num` where `deg = maxdegree(p, var)` and `num` is the number of
terms `t` such that `degree(t, var) == deg`.
"""
function deg_num_leading_terms(p::AbstractPolynomialLike, var)
    deg = 0
    num = 0
    for mono in monomials(p)
        d = degree(mono, var)
        if d > deg
            deg = d
            num = 1
        elseif d == deg
            num += 1
        end
    end
    return deg, num
end

Base.ndims(::Union{Type{<:AbstractPolynomialLike}, AbstractPolynomialLike}) = 0
Base.broadcastable(p::AbstractPolynomialLike) = Ref(p)


text:
# core definition of the AbstractGPUArray type

export AbstractGPUArray

"""
    AbstractGPUArray{T, N} <: DenseArray{T, N}

Supertype for `N`-dimensional GPU arrays (or array-like types) with elements of type `T`.
Instances of this type are expected to live on the host, see [`AbstractDeviceArray`](@ref)
for device-side objects.
"""
abstract type AbstractGPUArray{T, N} <: DenseArray{T, N} end

const AbstractGPUVector{T} = AbstractGPUArray{T, 1}
const AbstractGPUMatrix{T} = AbstractGPUArray{T, 2}
const AbstractGPUVecOrMat{T} = Union{AbstractGPUArray{T, 1}, AbstractGPUArray{T, 2}}

device(::AbstractGPUDevice) = error("Not implemented") # COV_EXCL_LINE
backend(::Type{<:AbstractGPUDevice}) = error("Not implemented") # COV_EXCL_LINE

# convenience aliases for working with wrapped arrays

const WrappedGPUArray{T,N} = WrappedArray{T,N,AbstractGPUArray,AbstractGPUArray{T,N}}

const AbstractOrWrappedGPUArray{T,N} =
    Union{AbstractGPUArray{T,N},
          WrappedArray{T,N,AbstractGPUArray,AbstractGPUArray{T,N}}}

# input/output

## serialization

using Serialization: AbstractSerializer, serialize_type

function Serialization.serialize(s::AbstractSerializer, t::T) where T <: AbstractGPUArray
    serialize_type(s, T)
    serialize(s, Array(t))
end
function Serialization.deserialize(s::AbstractSerializer, ::Type{T}) where T <: AbstractGPUArray
    A = deserialize(s)
    T(A)
end

## convert to CPU (keeping wrapper type)

Adapt.adapt_storage(::Type{<:Array}, xs::AbstractArray) = convert(Array, xs)
convert_to_cpu(xs) = adapt(Array, xs)

## showing

# display
Base.print_array(io::IO, X::AbstractOrWrappedGPUArray) =
    Base.print_array(io, convert_to_cpu(X))

# show
Base._show_nonempty(io::IO, X::AbstractOrWrappedGPUArray, prefix::String) =
    Base._show_nonempty(io, convert_to_cpu(X), prefix)
Base._show_empty(io::IO, X::AbstractOrWrappedGPUArray) =
    Base._show_empty(io, convert_to_cpu(X))
Base.show_vector(io::IO, v::AbstractOrWrappedGPUArray, args...) =
    Base.show_vector(io, convert_to_cpu(v), args...)

## collect to CPU (discarding wrapper type)

collect_to_cpu(xs::AbstractArray) = collect(convert_to_cpu(xs))
Base.collect(X::AbstractOrWrappedGPUArray) = collect_to_cpu(X)

# memory copying

## basic linear copies of identically-typed memory

# expects the GPU array type to have linear `copyto!` methods (i.e. accepting an integer
# offset and length) from and to CPU arrays and between GPU arrays.

for (D, S) in ((AbstractOrWrappedGPUArray, Array),
               (Array, AbstractOrWrappedGPUArray),
               (AbstractOrWrappedGPUArray, AbstractOrWrappedGPUArray))
    @eval begin
        function Base.copyto!(dest::$D{<:Any, N}, rdest::NTuple{N, UnitRange},
                              src::$S{<:Any, N}, ssrc::NTuple{N, UnitRange}) where {N}
            drange = CartesianIndices(rdest)
            srange = CartesianIndices(ssrc)
            copyto!(dest, drange, src, srange)
        end

        function Base.copyto!(dest::$D, d_range::CartesianIndices{1},
                              src::$S, s_range::CartesianIndices{1})
            len = length(d_range)
            if length(s_range) != len
                throw(ArgumentError("Copy range needs same length. Found: dest: $len, src: $(length(s_range))"))
            end
            len == 0 && return dest
            d_offset = first(d_range)[1]
            s_offset = first(s_range)[1]
            copyto!(dest, d_offset, src, s_offset, len)
        end

        function Base.copyto!(dest::$D, src::$S)
            len = length(src)
            len == 0 && return dest
            copyto!(dest, 1, src, 1, len)
        end
    end
end

# kernel-based variant for copying between wrapped GPU arrays

function linear_copy_kernel!(ctx::AbstractKernelContext, dest, dstart, src, sstart, n)
    i = linear_index(ctx)-1
    if i < n
        @inbounds dest[dstart+i] = src[sstart+i]
    end
    return
end

function Base.copyto!(dest::AbstractOrWrappedGPUArray, dstart::Integer,
                      src::AbstractOrWrappedGPUArray, sstart::Integer, n::Integer)
    n == 0 && return dest
    n < 0 && throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative")))
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    (checkbounds(Bool, destinds, dstart) && checkbounds(Bool, destinds, dstart+n-1)) || throw(BoundsError(dest, dstart:dstart+n-1))
    (checkbounds(Bool, srcinds, sstart)  && checkbounds(Bool, srcinds, sstart+n-1))  || throw(BoundsError(src,  sstart:sstart+n-1))

    gpu_call(linear_copy_kernel!,
             dest, dstart, src, sstart, n;
             total_threads=n)
    return dest
end

# variants that materialize the GPU wrapper before copying from or to the CPU

function Base.copyto!(dest::Array, dstart::Integer,
                      src::WrappedGPUArray, sstart::Integer, n::Integer)
    temp = similar(src, n)
    copyto!(temp, 1, src, sstart, n)
    copyto!(dest, dstart, temp, 1, n)
    return dest
end

function Base.copyto!(dest::WrappedGPUArray, dstart::Integer,
                      src::Array, sstart::Integer, n::Integer)
    temp = similar(dest, n)
    copyto!(temp, 1, src, sstart, n)
    copyto!(dest, dstart, temp, 1, n)
    return dest
end

## generalized blocks of heterogeneous memory

function cartesian_copy_kernel!(ctx::AbstractKernelContext, dest, dest_offsets, src, src_offsets, shape, length)
    i = linear_index(ctx)
    if i <= length
        idx = CartesianIndices(shape)[i]
        @inbounds dest[idx + dest_offsets] = src[idx + src_offsets]
    end
    return
end

function Base.copyto!(dest::AbstractOrWrappedGPUArray{<:Any, N}, destcrange::CartesianIndices{N},
                      src::AbstractOrWrappedGPUArray{<:Any, N}, srccrange::CartesianIndices{N}) where {N}
    shape = size(destcrange)
    if shape != size(srccrange)
        throw(DimensionMismatch("Ranges don't match their size. Found: $shape, $(size(srccrange))"))
    end
    len = length(destcrange)
    len == 0 && return dest

    dest_offsets = first(destcrange) - oneunit(CartesianIndex{N})
    src_offsets = first(srccrange) - oneunit(CartesianIndex{N})
    gpu_call(cartesian_copy_kernel!,
             dest, dest_offsets, src, src_offsets, shape, len;
             total_threads=len)
    dest
end

# XXX: these generalizations between non-linear CPU and GPU memory are very inefficient,
#       because it first materializes as linear memory.
# TODO: loop over linear sections of memory and perform linear copies

# NOTE: typed with Array because of ambiguities

function Base.copyto!(dest::AbstractGPUArray{T, N}, destcrange::CartesianIndices{N},
                      src::Array{T, N}, srccrange::CartesianIndices{N}) where {T, N}
    src_gpu = typeof(dest)(map(idx-> src[idx], srccrange))
    nrange = CartesianIndices(size(src_gpu))
    copyto!(dest, destcrange, src_gpu, nrange)
    dest
end

function Base.copyto!(dest::Array{T, N}, destcrange::CartesianIndices{N},
                      src::AbstractGPUArray{T, N}, srccrange::CartesianIndices{N}) where {T, N}
    dest_gpu = similar(src, size(destcrange))
    nrange = CartesianIndices(size(dest_gpu))
    copyto!(dest_gpu, nrange, src, srccrange)
    copyto!(dest, destcrange, Array(dest_gpu), nrange)
    dest
end

## other

Base.copy(x::AbstractGPUArray) = error("Not implemented") # COV_EXCL_LINE

Base.deepcopy(x::AbstractGPUArray) = copy(x)

# reinterpret

#=
copied from julia base/array.jl
Copyright (c) 2009-2016: Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors:

https://github.com/JuliaLang/julia/contributors

Permission is hereby granted, free of charge, to any person obtaining a copie 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.
=#

import Base.reinterpret

"""
    unsafe_reinterpret(T, a, dims)

Reinterpret the array `a` to have a new element type `T` and size `dims`.
"""
function unsafe_reinterpret end

function reinterpret(::Type{T}, a::AbstractGPUArray{S,1}) where T where S
    nel = (length(a)*sizeof(S)) ÷ sizeof(T)
    # TODO: maybe check that remainder is zero?
    return reinterpret(T, a, (nel,))
end

function reinterpret(::Type{T}, a::AbstractGPUArray{S}) where T where S
    if sizeof(S) != sizeof(T)
        throw(ArgumentError("result shape not specified"))
    end
    reinterpret(T, a, size(a))
end

function reinterpret(::Type{T}, a::AbstractGPUArray{S}, dims::NTuple{N, Integer}) where T where S where N
    if !isbitstype(T)
        throw(ArgumentError("cannot reinterpret Array{$(S)} to ::Type{Array{$(T)}}, type $(T) is not a bits type"))
    end
    if !isbitstype(S)
        throw(ArgumentError("cannot reinterpret Array{$(S)} to ::Type{Array{$(T)}}, type $(S) is not a bits type"))
    end
    nel = div(length(a)*sizeof(S),sizeof(T))
    if prod(dims) != nel
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(nel)"))
    end
    unsafe_reinterpret(T, a, dims)
end

function Base._reshape(A::AbstractGPUArray{T}, dims::Dims) where T
    n = length(A)
    prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
    return unsafe_reinterpret(T, A, dims)
end
#ambig
function Base._reshape(A::AbstractGPUArray{T, 1}, dims::Tuple{Integer}) where T
    n = Base._length(A)
    prod(dims) == n || throw(DimensionMismatch("parent has $n elements, which is incompatible with size $dims"))
    return unsafe_reinterpret(T, A, dims)
end

# filtering

# TODO: filter!

# revert of JuliaLang/julia#31929
Base.filter(f, As::AbstractGPUArray) = As[map(f, As)::AbstractGPUArray{Bool}]

ext: jl


text:
using Test
using QuantumControlBase

@testset "simple getcontrols" begin

    H = (nothing, (nothing, nothing))
    @test_throws TypeError getcontrols(H)

    ϵ₁ = t -> t
    ϵ₂ = t -> 0.0
    H = (nothing, (nothing, ϵ₁), (nothing, ϵ₂))
    @test getcontrols(H) == (ϵ₁, ϵ₂)

    u₁ = [0.1, 1.0]
    u₂ = [0.1, 2.0]
    H = (nothing, (nothing, u₁), (nothing, u₂))
    @test getcontrols(H) == (u₁, u₂)

end


text:
module CPUOptimizedConvNet

using LoopVectorization, PaddedMatrices
using PaddedMatrices.VectorizationBase: maybestaticsize, maybestaticfirst, maybestaticlast, staticm1, staticp1

export clippedimageview, zero_edges!

include("utils.jl")
# include("optimize.jl")
include("convolution.jl")
include("maxpool.jl")
include("logitcrossentropy.jl")

end


text:
using ReduxImGui
using Documenter

makedocs(;
    modules=[ReduxImGui],
    authors="Yupei Qi <qiyupei@gmail.com>",
    repo="https://github.com/Gnimuc/ReduxImGui.jl/blob/{commit}{path}#L{line}",
    sitename="ReduxImGui.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Gnimuc.github.io/ReduxImGui.jl",
        assets=String[],
    ),
    pages=[
       "Introduction" => "index.md",
       "API Reference" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/Gnimuc/ReduxImGui.jl.git",
)


text:
struct WideTreeAgnosticSysIDAgent
    widetree::WideTreeEnv
    model_class::Array{WideTreeEnv}
end

function WideTreeAgnosticSysIDAgent(widetree::WideTreeEnv, rng::MersenneTwister)
    model_class::Array{WideTreeEnv} = []
    # Add bad model
    push!(model_class, WideTreeEnv(widetree.n_leaves, "bad", rng))
    # Add several good models
    for _ in 1:5
        push!(model_class, WideTreeEnv(widetree.n_leaves, "good", rng))
    end
    WideTreeAgnosticSysIDAgent(widetree, model_class)
end

function run(
    agent::WideTreeAgnosticSysIDAgent,
    rng::MersenneTwister;
    debug = false,
    value_aware = true,
)
    # Choose an initial model
    model = agent.model_class[1]
    # Compute policy and values
    policy, values = optimal_policy_and_values(model, rng)
    # Initialize losses
    losses = [0.0 for _ in 1:length(agent.model_class)]

    # Start iterations
    n = 100  # number of iterations
    m = 10  # number of rollouts using policy
    p = 10  # number of rollouts using random policy with random start

    for _ in 1:n
        all_transitions::Array{Tuple{Int64, Int64, Int64, Int64}} = []
        performance = 0.
        # Do rollouts in the real world using policy to collect data
        for _ in 1:m
            transitions, total_return = simulate_episode(
                agent.widetree, policy, rng,
            )
            performance += total_return
            all_transitions = vcat(all_transitions, transitions)
        end
        println("Cost to go ", performance/m)
        # Do rollouts in the real world using random policy to collect data
        for _ in 1:p
            random_policy, _ = random_policy_and_values(
                model, rng,
            )
            transitions, _ = simulate_episode(
                agent.widetree, random_policy, rng,
            )
            all_transitions = vcat(all_transitions, transitions)
        end
        # Collected data, now update model
        for j in 1:length(agent.model_class)
            if value_aware
                losses[j] += compute_model_advantage_loss(
                    agent, all_transitions, values, agent.model_class[j], rng,
                ) / (m + p)
            else
                losses[j] += compute_classification_loss(
                    agent, all_transitions, agent.model_class[j], rng,
                ) / (m + p)
            end
        end
        # Find best model
        # TODO: Doing FTL now, but need to do better
        model = agent.model_class[argmin(losses)]
        println(losses)
        println("Selected ", argmin(losses))
        # Compute policy and values
        policy, values = optimal_policy_and_values(model, rng)
    end
end

function simulate_episode(
    widetree::WideTreeEnv,
    policy::Array{Int64},
    rng::MersenneTwister,
)
    state = init(widetree)
    transitions::Array{Tuple{Int64, Int64, Int64, Int64}} = []
    total_return = 0
    while !checkTerminal(widetree, state)
        action = policy[state]
        next_state, cost = step(widetree, state, action, rng)
        push!(transitions, (state, action, cost, next_state))
        total_return += cost
        state = next_state
    end

    transitions, total_return
end

function compute_model_advantage_loss(
    agent::WideTreeAgnosticSysIDAgent,
    transitions::Array{Tuple{Int64, Int64, Int64, Int64}},
    values::Array{Int64},
    model::WideTreeEnv,
    rng::MersenneTwister,
)
    loss = 0.
    for transition in transitions
        predicted_state, _ = step(model, transition[1], transition[2], rng)
        loss += abs(values[transition[4]] - values[predicted_state])
    end
    loss
end

function compute_classification_loss(
    agent::WideTreeAgnosticSysIDAgent,
    transitions::Array{Tuple{Int64, Int64, Int64, Int64}},
    model::WideTreeEnv,
    rng::MersenneTwister,
)
    loss = 0.
    for transition in transitions
        predicted_state, _ = step(model, transition[1], transition[2], rng)
        loss += 1 - Float64(predicted_state == transition[4])
    end
    loss
end


text:
if basename(pwd()) == "aoc"
    cd("2016/15")
end

function loaddiscs(filename::AbstractString)
    [
        begin
            a, b, c = parse.(Int,
                [match.(
                    r"Disc #(\d+) has (\d+) positions; at time=0, it is at position (\d+)\.",
                    line)...]
            )
            [(a + c) % b, b]
        end
        for line in eachline(filename)
    ]
end

part1(filename::AbstractString) = capsuledrop(loaddiscs(filename))
function capsuledrop(discs)
    for starttime in Iterators.countfrom(0)
        if all(iszero, (starttime + a) % b for (a, b) in discs)
            return starttime
        end
    end
end

@assert part1("example.txt") == 5
part1("input.txt")

part2(filename::AbstractString) = capsuledrop(push!(loaddiscs(filename), [7, 11]))
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using Images, ImageFeatures, CoordinateTransformations, ImageDraw

img1 = Gray.(load("cat-3417184_640.jpg"))
img2 = Gray.(load("cat-3417184_640_watermarked.jpg"))

# rotate img2 around the center and resize it
rot = recenter(RotMatrix(5pi/6), [size(img2)...] .÷ 2)
tform = rot ∘ Translation(-50, -40)
img2 = warp(img2, tform, indices(img2))
img2 = imresize(img2, Int.(trunc.(size(img2) .* 0.7)))

keypoints_1 = Keypoints(fastcorners(img1, 12, 0.35));
keypoints_2 = Keypoints(fastcorners(img2, 12, 0.35));

freak_params = FREAK()
desc_1, ret_keypoints_1 = create_descriptor(img1, keypoints_1, freak_params);
desc_2, ret_keypoints_2 = create_descriptor(img2, keypoints_2, freak_params);
matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.2)

# create
img3 = zeros(size(img1, 1), size(img2, 2))
img3[1:size(img2, 1), 1:size(img2, 2)] = img2

grid = Gray.(hcat(img1, img3))
offset = CartesianIndex(0, size(img1, 2))
map(m -> draw!(grid, LineSegment(m[1], m[2] + offset)), matches)

imshow(grid)


text:
using {{{PKG}}}
using Documenter

makedocs(;
    modules=[{{{PKG}}}],
    authors="{{{AUTHORS}}}",
    repo="https://{{{REPO}}}/blob/{commit}{path}#L{line}",
    sitename="{{{PKG}}}.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
{{#CANONICAL}}
        canonical="{{{CANONICAL}}}",
{{/CANONICAL}}
        assets={{^HAS_ASSETS}}String{{/HAS_ASSETS}}[{{^HAS_ASSETS}}],{{/HAS_ASSETS}}
{{#ASSETS}}
            "assets/{{{.}}}",
{{/ASSETS}}
{{#HAS_ASSETS}}
        ],
{{/HAS_ASSETS}}
    ),
    pages=[
        "Home" => "index.md",
    ],
{{#MAKEDOCS_KWARGS}}
    {{{first}}}={{{second}}},
{{/MAKEDOCS_KWARGS}}
)
{{#HAS_DEPLOY}}

deploydocs(;
    repo="{{{REPO}}}",
)
{{/HAS_DEPLOY}}


text:
module Division
    export thirt
    const RESTS = [1, 10, 9, 12, 3, 4]
    function thirt(n::Int64)::Int64
        xi = n
        while true
            xi,xo = sum( p -> *(p...), zip( digits(xi), Iterators.cycle(RESTS) ) ), xi
            if xi ≥ xo break end
        end
        xi
    end
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using ClangCompiler
using ClangCompiler.LLVM

julia_include_dir = joinpath(Sys.BINDIR, "..", "include", "julia") |> normpath

src = joinpath(@__DIR__, "embedding.cpp")

args = get_compiler_args(; version=v"7.1.0")
push!(args, "-std=c++14")
Sys.isapple() && push!(args, "-stdlib=libc++")
push!(args, "-I$julia_include_dir")

# create JIT, generate LLVM IR, and call function
jit = LLJIT(;tm=JITTargetMachine())
irgen = IRGenerator(src, args)
cc = CXCompiler(irgen, jit)
link_process_symbols(cc)
compile(cc)

addr = lookup(jit, "main")
@eval main() = ccall($(pointer(addr)), Cint, ())

main()

# clean up
dispose(cc)


text:
# Helper types for implementing `push`

struct SingletonVector{T} <: AbstractVector{T}
    value::T
end

Base.size(::SingletonVector) = (1,)

function Base.getindex(v::SingletonVector, i::Integer)
    @boundscheck i == 1 || throw(BoundsError(v, i))
    return v.value
end

# Define table interface as a `SingletonVector{<:NamedTuple}`:
Tables.istable(::Type{<:SingletonVector{<:NamedTuple{names}}}) where {names} =
    @isdefined(names)
Tables.rowaccess(::Type{<:SingletonVector{<:NamedTuple{names}}}) where {names} =
    @isdefined(names)
Tables.columnaccess(::Type{<:SingletonVector{<:NamedTuple{names}}}) where {names} =
    @isdefined(names)

Tables.rows(x::SingletonVector{<:NamedTuple}) = [x.value]
Tables.columns(x::SingletonVector{<:NamedTuple{names}}) where {names} =
    NamedTuple{names}(map(x -> [x], Tuple(x.value)))

struct SingletonDict{K,V} <: AbstractDict{K,V}
    key::K
    value::V
end

Base.iterate(d::SingletonDict) = (d.key => d.value, nothing)
Base.iterate(d::SingletonDict, ::Nothing) = nothing

function Base.getindex(d::SingletonDict{K}, key::K) where {K}
    @boundscheck d.key == key || throw(BoundsError(d, key))
    return d.value
end

Base.length(::SingletonDict) = 1


text:
"""
A [`Plugin`](@ref) that only adds a [`Badge`](@ref) to the [`Readme`](@ref) file.

Concrete subtypes only need to implement a [`badges`](@ref) method.
"""
abstract type BadgePlugin <: Plugin end

"""
    BlueStyleBadge()

Adds a [`BlueStyle`](https://github.com/invenia/BlueStyle) badge to the [`Readme`](@ref) file.
"""
struct BlueStyleBadge <: BadgePlugin end

function badges(::BlueStyleBadge)
    return Badge(
        "Code Style: Blue",
        "https://img.shields.io/badge/code%20style-blue-4495d1.svg",
        "https://github.com/invenia/BlueStyle",
    )
end

"""
    ColPracBadge()

Adds a [`ColPrac`](https://github.com/SciML/ColPrac) badge to the [`Readme`](@ref) file.
"""
struct ColPracBadge <: BadgePlugin end

function badges(::ColPracBadge)
    return Badge(
        "ColPrac: Contributor's Guide on Collaborative Practices for Community Packages",
        "https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet",
        "https://github.com/SciML/ColPrac",
    )
end

"""
    PkgEvalBadge()

Adds a [`PkgEval` badge](https://github.com/JuliaCI/NanosoldierReports#pkgeval-reports) to the [`Readme`](@ref) file.
"""
struct PkgEvalBadge <: BadgePlugin end

function badges(::PkgEvalBadge)
    return Badge(
        "PkgEval",
        "https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/{{{PKG1}}}/{{{PKG}}}.svg",
        "https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/report.html"
    )
end

function view(::PkgEvalBadge, t::Template, pkg::AbstractString)
    return Dict("PKG1" => first(pkg), "PKG" => pkg)
end


text:
module Underscore

"""
    Underscore.scan_ast(e, fv)

    is a utility function that scans an Expr e for a free variable fv, which
    must be a symbol.  Returns true if fv is a leaf anywhere in the ast.
"""
scan_ast(e::Any, fv::Symbol) = false
scan_ast(e::Symbol, fv::Symbol) = e == fv
scan_ast(e::Expr, fv::Symbol) = mapreduce((x) -> scan_ast(x, fv), |, false, e.args)

# package code goes here
"""
    \@_

gives a simple way of remapping multi-argument functions in pipelines.  Normally,
a julia pipeline requires a naked function with one argument to be the right
operand of the pipeline operator; this allows you to easily call multi-argument
functions in the pipeline operator by syntatically calling the remaining
parameters.  \@_ will create a lambda encapsulating the remaining parameters and
injecting the result as the first parameter.

#Examples

```jldocstring
julia> "hello" |> @_ println(" world")
hello world
```

you may also redirect in the case where you cannot intercept the first value.

```jldocstring
julia> [1, 2, 3] |> @_ reduce(+, 0, \$1)
6
```

finally, you can also pass multi-valued function results (aka tuples) into the
recieving function, as you might expect.

```jldocstring
julia> (1, 2) |> @_ sum(3, \$1, \$2)
6
```

"""
macro _(expr)
    #first, make sure that expr is a call.
    (expr.head == :call) || throw(ArgumentError("@_ must precede a function call."))
    #find out how many dollar sign arguments we're going to use.
    positionalparams = Set{Int64}()
    for param in expr.args
        if (param isa Expr) && (param.head == :$) && (length(param.args) == 1) && (param.args[1] isa Int64)
            push!(positionalparams, param.args[1])
        end
    end

    #find the tuple length.  This is the longest parameter in the tuple
    tuple_length = length(positionalparams) == 0 ? 1 : maximum(positionalparams)

    #select the correct variable to use as the free variable in the lambda.
    #this can get tricky since it's possible the lambda'd function call contains
    #other free variables that *look like* the one we want to use.  We'll have to
    #traverse the tree to make sure this isn't the case.  Allocate one extra
    #free variable to use as a splatting variable, if needed.
    test_index = 1
    free_variables = Symbol[]
    for idx = 1:(tuple_length + 1)
        while true
            free_var = Symbol(:x, test_index)
            scan_ast(expr, free_var) || break;    #break if we can find a safe value.
            test_index += 1
            (test_index > 100) && throw(ArgumentError("it seems like you have an awful lot of free x variables"))
        end
        push!(free_variables, Symbol(:x, test_index))
        test_index += 1
    end

    if (length(positionalparams) == 0)
        #prepend the free variable into the ast of the expression.
        insert!(expr.args, 2, free_variables[1])
    else
        #do a splicing operation to put in the correct free variables.
        for idx = 1:length(expr.args)
            param = expr.args[idx]
            if (param isa Expr) && (param.head == :$) && (length(param.args) == 1) && (param.args[1] isa Int64)
                substitution = param.args[1]
                splice!(expr.args, idx:idx, [free_variables[substitution]])
            end
        end
    end

    #create a shell, empty lambda into which we will put our free variables, except
    #the tailing one.
    generated_lambda = :(()-> $expr)
    append!(generated_lambda.args[1].args, free_variables[1:end-1])

    if tuple_length <= 1
        #return the function and escape.
        esc(generated_lambda)
    else
        #splat a passed tuple.
        t = free_variables[end]
        esc(:(($t)->($generated_lambda)($t...)))
    end
end

#this doesn't collide with anything in the Base namespace, so go ahead and
#export it for convenience.
export @_

end # module


text:
# Use
#    @warnpcfail precompile(args...)
# if you want to be warned when a precompile directive fails
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol, Any},AbstractMatrix{T} where T})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol, Any},AbstractVector{T} where T,AbstractVector{T} where T,Function})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol, Any},Function,Number,Number})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol, Any},GroupBy,Any})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol, Any},Type{SliceIt},Any,Any,Any})
    Base.precompile(Tuple{typeof(RecipesBase.apply_recipe),AbstractDict{Symbol, Any},Vector{Function},Number,Number})
    Base.precompile(Tuple{typeof(_apply_type_recipe),Any,AbstractArray,Any})
    Base.precompile(Tuple{typeof(_apply_type_recipe),Any,Surface,Any})
    Base.precompile(Tuple{typeof(_compute_xyz),Vector{Float64},Function,Nothing,Bool})
    Base.precompile(Tuple{typeof(_extract_group_attributes),Vector{String},Vector{Float64}})
    Base.precompile(Tuple{typeof(_map_funcs),Function,StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}})
    Base.precompile(Tuple{typeof(_process_ribbon),Tuple{LinRange{Float64}, LinRange{Float64}},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_process_seriesrecipe),Any,Any})
    Base.precompile(Tuple{typeof(_process_seriesrecipes!),Any,Any})
    Base.precompile(Tuple{typeof(_scaled_adapted_grid),Function,Symbol,Symbol,Float64,Irrational{:π}})
    Base.precompile(Tuple{typeof(_series_data_vector),Int64,Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Matrix{Float64},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Surface{Matrix{Int64}},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{AbstractVector{Float64}},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{Function},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{Int64},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{Real},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{String},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{Union{Missing, Int64}},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{Vector{Float64}},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(_series_data_vector),Vector{Vector{T} where T},Dict{Symbol, Any}})
    Base.precompile(Tuple{typeof(filter_data!),Dict{Symbol, Any},Vector{Int64}})
    Base.precompile(Tuple{typeof(recipe_pipeline!),Any,Any,Any})
    Base.precompile(Tuple{typeof(unzip),Vector{Tuple{Float64, Float64}}})
    Base.precompile(Tuple{typeof(unzip),Vector{Tuple{Int64, Int64}}})
    Base.precompile(Tuple{typeof(unzip),Vector{Tuple{Int64, Real}}})
    Base.precompile(Tuple{typeof(unzip),Vector{Tuple{Vector{Float64}, Vector{Float64}}}})
    let fbody = try __lookup_kwbody__(which(_extract_group_attributes, (Vector{String},Vector{Float64},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Function,typeof(_extract_group_attributes),Vector{String},Vector{Float64},))
        end
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Carmichael 3 strong pseudoprimes

function carmichael{T<:Integer}(pmax::T)
    0 < pmax || throw(DomainError())
    car = T[]
    for p in primes(pmax)
        for h₃ in 2:(p-1)
            m = (p - 1)*(h₃ + p)
            pmh = mod(-p^2, h₃)
            for Δ in 1:(h₃+p-1)
                m%Δ==0 && Δ%h₃==pmh || continue
                q = div(m, Δ) + 1
                isprime(q) || continue
                r = div((p*q - 1), h₃) + 1
                isprime(r) && mod(q*r, (p-1))==1 || continue
                append!(car, [p, q, r])
            end
        end
    end
    reshape(car, 3, div(length(car), 3))
end

hi = 61
car = carmichael(hi)

curp = 0
tcnt = 0
print("Carmichael 3 (p\u00d7q\u00d7r) Pseudoprimes, up to p = ", hi, ":")
for j in sortperm(1:size(car)[2], by=x->(car[1,x], car[2,x], car[3,x]))
    p, q, r = car[:,j]
    c = prod(car[:,j])
    if p != curp
        curp = p
        print(@sprintf("\n\np = %d\n  ", p))
        tcnt = 0
    end
    if tcnt == 4
        print("\n  ")
        tcnt = 1
    else
        tcnt += 1
    end
    print(@sprintf("p\u00d7%d\u00d7%d = %d  ", q, r, c))
end
println("\n\n", size(car)[2], " results in total.")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Interpolations, Test, ForwardDiff

@testset "Issues" begin
    @testset "issue 34" begin
        A = rand(1:20, 100, 100)

        # In #34, this incantation throws
        itp = interpolate(A, BSpline(Quadratic(Flat(OnCell()))))
        # Sanity check that not only don't throw, but actually interpolate
        for i in 1:size(A,1), j in 1:size(A,2)
            @test itp[i,j] ≈ A[i,j]
        end
    end

    @testset "issue 129" begin
        xy = [3z+2y for z in range(0.0,stop=1,length=10),y in 0:0.2:1]
        itp = interpolate((1:10,1:6),xy,(Gridded(Linear()),NoInterp()))
        @test itp(3.3,1:6) == [itp(3.3,i) for i = 1:6]
        sitp = scale(itp ,range(0.0,stop=1.0,length=10),1:6)
        @test sitp(0.8,1:6) == [sitp(0.8,i) for i = 1:6]
    end

    @testset "issue 151" begin
        V = zeros(10,10,10,10)
        interpV = interpolate(V, BSpline(Cubic(Line(OnGrid()))))
        @test ndims(interpV) == 4
    end

    @testset "issue 158" begin
        A_x = 1.:2.:40.
        A = [log(x) for x in A_x]
        itp = interpolate(A, BSpline(Cubic(Line(OnGrid()))))
        sitp = scale(itp, A_x)
        @test_throws BoundsError sitp(42.0)
        @test_throws BoundsError sitp([3.0,42.0])
    end

    @testset "issue 165" begin
        V0 = Array{Float64, 2}(undef, 5, 6)
        yGrid = range(-1.0,stop= 1.0,length= 5)
        bGrid = range(2.0, stop=3.0, length=6)
        iV = scale(interpolate(V0, BSpline(Cubic(Line(OnCell())))), yGrid, bGrid)
        @test iV isa Interpolations.AbstractInterpolation
    end

    @testset "issue 191" begin
        function eff(p, thing)
            out = p[1] *  p[2]
            return out
        end
        function GradientWrapper(func, p, thing)
            f(x) = func(x, thing)
            return ForwardDiff.gradient(f,p)
        end
        io = IOBuffer()
        println(io, "Start.")
        n = 40
        z = [i * j for i in 1:n, j in 1:n]
        itp = interpolate(z, BSpline(Cubic(Natural(OnGrid()))))
        p = [5.0, 6.0]
        test = [itp, itp]
        println(io, "Using type $(typeof(test))")
        println(io, "Calling the derivative directly:")
        f(x) = eff(x,test)
        println(io, ForwardDiff.gradient(f, p))
        println(io, "Calling the wrapper:")
        println(io, GradientWrapper(eff, p, test))
        test = Array{Any}(undef, 2)
        test[1] = itp
        test[2] = itp
        println(io, "Using type $(typeof(test))")
        println(io, "Calling the derivative directly:")
        g(x) = eff(x, test)
        println(io, ForwardDiff.gradient(g, p))
        println(io, "Calling the wrapper:")
        println(io, GradientWrapper(eff, p, test))
        test = Array{AbstractInterpolation}(undef, 2)
        test[1] = itp
        test[2] = itp
        println(io, "Using type $(typeof(test))")
        println(io, "Calling the derivative directly:")
        h(x) = eff(x, test)
        println(io, ForwardDiff.gradient(h, p))
        println(io, "Calling the wrapper -- this fails:")
        println(io, GradientWrapper(eff, p, test))
        println(io, "Done.")
        str = String(take!(io))
        @test endswith(str, "Done.\n")
    end

    @testset "issue 200" begin
        grid = range(0, stop=2.0, length=10)
        vals = [-exp(-x) for x in grid]
        itp = interpolate(vals, BSpline(Cubic(Line(OnGrid()))))
        sitp = scale(itp, grid)
        s = sitp(0.0)
        @test sitp([0.0, 0.0]) == [s, s]
    end

    @testset "issue 202" begin
        xs = 1:10
        itp = scale(interpolate(zeros(10), BSpline(Cubic(Periodic(OnGrid())))), xs)
        extp = extrapolate(itp, Periodic())
        @test Interpolations.gradient(extp, 2.) == [0.0]
        @test_throws ErrorException Interpolations.gradient(extp, 2.0, 2.0)

        # The issue title says 3d, so let's try a 3d case
        A = (0.2.*(1:3)) .* (0.3.*(1:3))' .* reshape(0.4.*(1:3), 1, 1, 3)
        itp = interpolate(A, BSpline(Linear()))
        etp = extrapolate(itp, Line())
        @test Interpolations.gradient(etp, 2, 2, 2) ≈ [0.2*0.6*0.8, 0.4*0.3*0.8, 0.4*0.6*0.4]
        @test Interpolations.gradient(etp, 10, 2, 2) ≈ [0.2*0.6*0.8, 0.6*0.3*0.8, 0.6*0.6*0.4]
        @test Interpolations.gradient(etp, 2, 10, 2) ≈ [0.2*0.9*0.8, 0.4*0.3*0.8, 0.4*0.9*0.4]
        @test Interpolations.gradient(etp, 2, 2, 10) ≈ [0.2*0.6*1.2, 0.4*0.3*1.2, 0.4*0.6*0.4]
    end

    @testset "issue 213" begin
        function build_itp()
            A_x = 1.0:2.0:40.0
            A = A_x.^2
            knots = (A_x,)
            interpolate(knots, A, Gridded(Linear()))
        end
        function itp_test(x::AbstractVector)
            itp = build_itp()
            return itp(x...)
        end

        j = ForwardDiff.gradient(itp_test, [2.0])
        @test j ≈ Interpolations.gradient(build_itp(), 2.0)
    end

    @testset "issue 335" begin
        function build_itp()
            A_x = 1:5
            A = A_x
            knots = (A_x,)
            LinearInterpolation(knots, A, extrapolation_bc=Line())
        end
        function itp_test(x::AbstractVector)
            itp = build_itp()
            return itp(x...)
        end

        @test ForwardDiff.gradient(itp_test, [1.0]) ≈ Interpolations.gradient(build_itp(), 1.0)
        @test ForwardDiff.gradient(itp_test, [5.0]) ≈ Interpolations.gradient(build_itp(), 5.0)
    end
    @testset "issue 430" begin
        mat1 = rand(Float64, 64, 64)
        itp = LinearInterpolation((1:64, 1:64,), mat1)
        function itp_test(x::AbstractVector)
            return itp(x...)
        end
        @test ForwardDiff.gradient(itp_test, [1.0, 1.0]) ≈ Interpolations.gradient(itp, 1.0, 1.0)
        @test ForwardDiff.gradient(itp_test, [3.0, 2.0]) ≈ Interpolations.gradient(itp, 3.0, 2.0)
    end

end


text:
# test resolutuion of #66

using BEAST
using Test
using CompScienceMeshes
using LinearAlgebra

function hassemble(operator::BEAST.AbstractOperator,
    test_functions,
    trial_functions)

    blkasm = BEAST.blockassembler(operator, test_functions, trial_functions)

    function assembler(Z, tdata, sdata)
        store(v,m,n) = (Z[m,n] += v)
        blkasm(tdata,sdata,store)
    end

    mat = zeros(scalartype(operator),
                numfunctions(test_functions),
                numfunctions(trial_functions))

    assembler(mat, 1:numfunctions(test_functions), 1:numfunctions(trial_functions))
    return mat
end

c = 3e8
μ = 4π * 1e-7
ε = 1/(μ*c^2)
f = 1e8
λ = c/f
k = 2π/λ
ω = k*c
η = sqrt(μ/ε)

a = 1
Γ = CompScienceMeshes.meshcuboid(a,a,a,0.2)

𝓣 = Maxwell3D.singlelayer(wavenumber=k)
𝓚 = Maxwell3D.doublelayer(wavenumber=k)

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

println("Number of RWG functions: ", numfunctions(X))

T_blockassembler = hassemble(𝓣, X, X)
T_standardassembler = assemble(𝓣, X, X)

@test norm(T_blockassembler - T_standardassembler)/norm(T_standardassembler) ≈ 0.0 atol=1e-14

T_bc_blockassembler = hassemble(𝓣, Y, Y)
T_bc_standardassembler = assemble(𝓣, Y, Y)

@test norm(T_bc_blockassembler - T_bc_standardassembler)/norm(T_bc_standardassembler) ≈ 0.0 atol=1e-14

K_mix_blockassembler = hassemble(𝓚,Y,X)
K_mix_standardassembler = assemble(𝓚,Y,X)

T_mix_blockassembler = hassemble(𝓣, Y, X)
T_mix_standardassembler = assemble(𝓣, Y, X)

@test norm(K_mix_blockassembler - K_mix_standardassembler)/norm(K_mix_standardassembler) ≈ 0.0 atol=1e-14
@test norm(T_mix_blockassembler - T_mix_standardassembler)/norm(T_mix_standardassembler) ≈ 0.0 atol=1e-14


text:
# Autogenerated wrapper script for Catalyst_jll for x86_64-w64-mingw32-cxx03
export catalyst_replay_exe, libcatalyst

JLLWrappers.@generate_wrapper_header("Catalyst")
JLLWrappers.@declare_library_product(libcatalyst, "libcatalyst.dll")
JLLWrappers.@declare_executable_product(catalyst_replay_exe)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcatalyst,
        "bin\\libcatalyst.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        catalyst_replay_exe,
        "bin\\catalyst_replay.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
Chance Sampling Counterfactual Regret Minimization
=#

struct CSCFRSolver{K,G,I} <: AbstractCFRSolver{K,G,I}
    I::Dict{K, I}
    game::G
end

"""
    `CSCFRSolver(game::Game{H,K}; debug::Bool=false)`

Instantiate chance sampling CFR solver with some `game`.

If `debug=true`, record history of strategies over training period, allowing
for training history of individual information states to be plotted with
`Plots.plot(is::DebugInfoState)`

"""
function CSCFRSolver(game::Game{H,K}; debug::Bool=false) where {H,K}
    if debug
        return CSCFRSolver(Dict{K, DebugInfoState}(), game)
    else
        return CSCFRSolver(Dict{K, InfoState}(), game)
    end
end

function CFR(solver::CSCFRSolver, h, i, t, π_1, π_2)
    game = solver.game
    if isterminal(game, h)
        return utility(game, i, h)
    elseif iszero(player(game, h)) # chance player
        a = chance_action(game, h)
        h′ = next_hist(game,h,a)
        return CFR(solver, h′, i, t, π_1, π_2)
    end
    I = infoset(solver, h)
    A = actions(game, h)

    v_σ = 0.0
    v_σ_Ia = I._tmp_σ

    for (k,a) in enumerate(A)
        h′ = next_hist(game, h, a)
        if player(game, h) === 1
            v_σ_Ia[k] = CFR(solver, h′, i, t, I.σ[k]*π_1, π_2)
        else
            v_σ_Ia[k] = CFR(solver, h′, i, t, π_1, I.σ[k]*π_2)
        end
        v_σ += I.σ[k]*v_σ_Ia[k]
    end

    if player(game, h) == i
        π_i = i == 1 ? π_1 : π_2
        π_ni = i == 1 ? π_2 : π_1
        for (k,a) in enumerate(A)
            I.r[k] += π_ni*(v_σ_Ia[k] - v_σ)
            I.s[k] += π_i*I.σ[k]
        end
    end

    return v_σ
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
import Base.+
import Base.-
import Base.*

@generated function +{R1,R2}( a1::BInt{R1}, a2::BInt{R2} )
        l1,u1 = get_range_bounds( R1 )
        l2,u2 = get_range_bounds( R2 )
        RangeOut = make_range( l1+l2, u1+u2 )
        Tout = default_int_type( RangeOut )
        # convert the args to Tout ensure no overflow
        return :(
                UncheckedBInt{ $RangeOut, $Tout }( $Tout(a1.n) + $Tout(a2.n) )
                )

end

@generated function -{R1,R2}( a1::BInt{R1}, a2::BInt{R2} )
        l1,u1 = get_range_bounds( R1 )
        l2,u2 = get_range_bounds( R2 )
        RangeOut = make_range( l1-u2, u1-l2 )
        Tout = default_int_type( RangeOut )

        @dd @show RangeOut Tout
        # convert the args to Tout ensure no overflow
        return :(
                UncheckedBInt{ $RangeOut, $Tout }( $Tout(a1.n) - $Tout(a2.n) )
                )

end

@generated function *{R1,R2}( a1::BInt{R1}, a2::BInt{R2} )
        l1,u1 = get_range_bounds( R1 )
        l2,u2 = get_range_bounds( R2 )

        corners = l1*l2, l1*u2, u1*l2, u1*u2
        lower = min( corners... )
        upper = max( corners... )

        @dd @show l1 u1 l2 u2 corners lower upper

        RangeOut = make_range( lower, upper )
        Tout = default_int_type( RangeOut )

        @dd @show RangeOut Tout
        # convert the args to Tout ensure no overflow
        return :(
                UncheckedBInt{ $RangeOut, $Tout }( $Tout(a1.n) * $Tout(a2.n) )
                )
end

import Base.<<
@generated function <<{R1,R2}( a1::BInt{R1}, a2::BInt{R2} )
        l1,u1 = get_range_bounds( R1 )
        l2,u2 = get_range_bounds( R2 )

        lower = signbit(u2) ? l1 >> -u2 : l1 << u2
        upper = signbit(u2) ? u1 >> -u2 : u1 << u2

        RangeOut = make_range( lower, upper )
        Tout = default_int_type( RangeOut )

        # be careful with the type used for the actual <<
        # operation.  a2 could always be negative
        # and the result type be smaller than the arg types
        # but need to do the shift with larger type
        # also note that 100 << -1 is 0 in julia....
        # 100 << -1 is 50 in Julia 0.5

        if l2 >= 0
                # always positive shift
                # working type is Tout
                return :(
                        UncheckedBInt{ $RangeOut, $Tout }( $Tout(a1.n) << a2.n )
                        )
        elseif u2 <= 0
                # always negative shift
                # working type is the input type
                # avoid corner case of RHS being Int8(-128) and then
                # being negated
                return :(
                        UncheckedBInt{ $RangeOut, $Tout }( $Tout( a1.n >> -Int(a2.n) ) )
                        )
        else
                # mixed sign shift
                # working type is Tout
                return :(
                        UncheckedBInt{ $RangeOut, $Tout }(
                                signbit( a2.n ) ? $Tout(a1.n) >> -Int(a2.n)
                                        : $Tout(a1.n) << a2.n
                                )
                        )
        end
end
function <<{R1}( a1::BInt{R1}, a2::Unsigned )
        return default_int_type(R1)(a1.n) << a2
end

@generated function Base.abs{R}( x::BInt{R} )
        l,u = get_range_bounds( R )

        lower = (l <= 0 && u >= 0) ? 0 : min( abs(l), abs(u) )
        upper = max( abs(l), abs(u) )

        RangeOut = make_range( lower, upper )
        Tout = default_int_type( RangeOut )

        return :( UncheckedBInt{ $RangeOut, $Tout }( Base.abs( $Tout(x) ) ) )
end

# abs2(x) just returns x*x - but the bounds can be tighter than just letting the
# default implementation do it
@generated function Base.abs2{R}( x::BInt{R} )
        l,u = get_range_bounds( R )

        lower = (l <= 0 && u >= 0) ? 0 : min( abs(l), abs(u) )
        upper1 = max( abs(l), abs(u) )
        upper = upper1 * upper1

        RangeOut = make_range( lower, upper )
        Tout = default_int_type( RangeOut )

        return :( UncheckedBInt{ $RangeOut, $Tout }( $Tout(x) * $Tout(x) ) )
end

# there are corner cases where this fails
# need a function that gives the union type of the input and
# output range, and uses this for the working type
# and probably as the result type too.

# consider x with range 0:( typemax(Int64)+1 )
# result range is then typemin(Int64):0
# so using just the output range would suggest that Int64 was suitable
# when it is not, because converting ( typemax(Int64)+1 ) from an Int128
# to a Int64 throws an exception

@generated function -{R}( x::BInt{R} )
        l,u = get_range_bounds( R )

        lower = -u
        upper = -l

        RangeOut = make_range( lower, upper )
        Tout = default_int_type( range_union( R, RangeOut ) )

        return :( UncheckedBInt{ $RangeOut, $Tout }( -( $Tout(x) ) ) )
end

import Base.~
@generated function ~{R}( x::BInt{R} )
        l,u = get_range_bounds( R )

        lower = ~u
        upper = ~l

        RangeOut = make_range( lower, upper )
        Tout = default_int_type( range_union( R, RangeOut ) )

        return :( UncheckedBInt{ $RangeOut, $Tout }( ~( $Tout(x) ) ) )
end

Base.rand{R}(rng::AbstractRNG, ::Type{ BInt{R} } ) = UncheckedBInt{R}( rand(rng, R ) )
Base.rand{R,T}(rng::AbstractRNG, ::Type{ BInt{R,T} } ) = UncheckedBInt{R,T}( rand(rng, R ) )

import Base.<
import Base.<=
import Base.>
import Base.>=
import Base.==
import Base.!=

 <{   T1<:Integer,R2,T2<:Integer}( x::T1         , y::BInt{R2,T2} ) = x   < y.n
 <{R1,T1<:Integer,   T2<:Integer}( x::BInt{R1,T1}, y::T2          ) = x.n < y
 <{R1,T1<:Integer,R2,T2<:Integer}( x::BInt{R1,T1}, y::BInt{R2,T2} ) = x.n < y.n

<={   T1<:Integer,R2,T2<:Integer}( x::T1         , y::BInt{R2,T2} ) = x   <= y.n
<={R1,T1<:Integer,   T2<:Integer}( x::BInt{R1,T1}, y::T2          ) = x.n <= y
<={R1,T1<:Integer,R2,T2<:Integer}( x::BInt{R1,T1}, y::BInt{R2,T2} ) = x.n <= y.n

 >{   T1<:Integer,R2,T2<:Integer}( x::T1         , y::BInt{R2,T2} ) = x   > y.n
 >{R1,T1<:Integer,   T2<:Integer}( x::BInt{R1,T1}, y::T2          ) = x.n > y
 >{R1,T1<:Integer,R2,T2<:Integer}( x::BInt{R1,T1}, y::BInt{R2,T2} ) = x.n > y.n

>={   T1<:Integer,R2,T2<:Integer}( x::T1         , y::BInt{R2,T2} ) = x   >= y.n
>={R1,T1<:Integer,   T2<:Integer}( x::BInt{R1,T1}, y::T2          ) = x.n >= y
>={R1,T1<:Integer,R2,T2<:Integer}( x::BInt{R1,T1}, y::BInt{R2,T2} ) = x.n >= y.n

=={   T1<:Integer,R2,T2<:Integer}( x::T1         , y::BInt{R2,T2} ) = x   == y.n
=={R1,T1<:Integer,   T2<:Integer}( x::BInt{R1,T1}, y::T2          ) = x.n == y
=={R1,T1<:Integer,R2,T2<:Integer}( x::BInt{R1,T1}, y::BInt{R2,T2} ) = x.n == y.n

!={   T1<:Integer,R2,T2<:Integer}( x::T1         , y::BInt{R2,T2} ) = x   != y.n
!={R1,T1<:Integer,   T2<:Integer}( x::BInt{R1,T1}, y::T2          ) = x.n != y
!={R1,T1<:Integer,R2,T2<:Integer}( x::BInt{R1,T1}, y::BInt{R2,T2} ) = x.n != y.n

# sign is supposed to return the same type as the arg
# so interpret this as returning a BInt
@generated function Base.sign{R}( x::BInt{R} )
        l,u = get_range_bounds( R )

        l1 = Base.sign(l)
        u1 = Base.sign(u)

        RangeOut = make_range( l1, u1 )
        Tout = default_int_type( range_union( R, RangeOut ) )

        return :( UncheckedBInt{ $RangeOut, $Tout }( $Tout( sign(x.n ) ) ) )
end

Base.signbit( x::BInt ) = Base.signbit( x.n )
Base.ispow2( x::BInt ) = Base.ispow2( x.n )

@generated function Base.nextpow2{R}( x::BInt{R} )
        l,u = get_range_bounds( R )

        l1 = Base.nextpow2(l)
        u1 = Base.nextpow2(u)

        RangeOut = make_range( l1, u1 )
        Tout = default_int_type( range_union( R, RangeOut ) )

        return :( UncheckedBInt{ $RangeOut, $Tout }( $Tout( nextpow2(x.n ) ) ) )
end

@generated function Base.prevpow2{R}( x::BInt{R} )
        l,u = get_range_bounds( R )

        l1 = Base.prevpow2(l)
        u1 = Base.prevpow2(u)

        RangeOut = make_range( l1, u1 )
        Tout = default_int_type( range_union( R, RangeOut ) )

        return :( UncheckedBInt{ $RangeOut, $Tout }( $Tout( prevpow2(x.n ) ) ) )
end

# could do subset of round, trunc etc