text:
################################################################################
# pset.jl
#
# Implementation of a persistent set type that uses the PTree type internally.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2020-2021 Noah C. Benson

import Base.IdSet

# #PLinearSet ==================================================================
# We construct these objects in an unfortunaetly complex way in order to reduce
# code duplication. The two kinds of persistent sets are very similar, just
# using different equality tests, so we can keep most of the code the same.
# See PDict.jl for information about this macro-based construction, used here.
macro _plinset_code(name::Symbol, eq, h)
    return quote
        struct $name{T} <: AbstractPSet{T}
            elements::Union{Nothing, Vector{T}}
            function $name{T}() where {T}
                return new{T}(nothing)
            end
            function $name{T}(u::Vector{T}, ::Val{:safe}) where {T}
                return length(u) == 0 ? new{T}(nothing) : new{T}(u)
            end
            function $name{T}(u::Vector{T}, ::Val{:unique}) where {T}
                return length(u) == 0 ? new{T}(nothing) : new{T}(copy(u))
            end
            function $name{T}(u::Vector{T}) where {T}
                (length(u) == 0) && return new{T}(nothing)
                v = Vector{T}(undef, length(u))
                kk = 0
                for ii in 1:length(u)
                    ui = u[ii]
                    found = false
                    for jj in 1:k
                        if v[ii] == ui
                            found = true
                            break
                        end
                    end
                    found && continue
                    k += 1
                    v[k] = ui
                end
                resize!(v, k)
                return new{T}(v)
            end
        end
        function $name{T}(d::$name{T}) where {T}
            return d
        end
        function $name{T}(t0, ts...) where {T}
            return $name{T}(T[t0, ts...])
        end
        function $name{T}(itr) where {T}
            return $name{T}(T[itr...])
        end
        # Document the equal/hash types.
        equalfn(::Type{T}) where {T <: $name} = $eq
        hashfn(::Type{T}) where {T <: $name} = $h
        equalfn(::$name) = $eq
        hashfn(::$name) = $h
        # Generic construction.
        $name() = $name{Any}()
        $name(::Tuple{}) = $name{Any}()
        $name(itr) = $name(itr, Base.IteratorEltype(itr))
        $name(itr, ::Base.HasEltype) = $name{eltype(itr)}(itr)
        $name(itr, ::Base.EltypeUnknown) = begin
            T = typejoin(map(typeof, itr)...)
            return $name{T}(itr)
        end
        # Base methods.
        Base.empty(s::$name{T}, ::Type{S}=T) where {T,S} = $name{S}()
        Base.length(s::$name) = begin
            ks = getfield(s, :elements)
            return (ks === nothing ? 0 : length(ks))
        end
        Base.IteratorEltype(::Type{$name{T}}) where {T} = Base.HasEltype()
        Base.eltype(::Type{$name{T}}) where {T} = T
        Base.eltype(::$name{T}) where {T} = T
        Base.IteratorSize(::Type{$name{T}}) where {T} = Base.HasLength()
        Base.iterate(u::$name{T}) where {T} = begin
            els = getfield(u, :elements)
            (els === nothing) && return nothing
            return iterate(els)
        end
        Base.iterate(u::$name{T}, state) where {T} = begin
            els = getfield(u, :elements)
            return iterate(els, state)
        end
        Base.in(x::S, u::$name{T}) where {S,T} = begin
            els = getfield(u, :elements)
            (els === nothing) && return false
            els::Vector{T}
            for el in els
                $eq(x,el) && return true
            end
            return false
        end
        # And the Air methods.
        Air.push(u::$name{T}, x::S) where {T,S} = begin
            els = getfield(u, :elements)
            (els === nothing) && return $name{T}(x)
            for el in els
                $eq(el,x) && return u
            end
            return $name{T}(push(els, x), Val{:safe}())
        end
        Air.delete(u::$name{T}, x::S) where {T,S} = begin
            els = getfield(u, :elements)
            (els === nothing) && return u
            els::Vector{T}
            n = length(els)
            if n == 1
                for ii in 1:n
                    if $eq(els[ii], x)
                        return $name{T}()
                    end
                end
            else
                for ii in 1:n
                    if $eq(els[ii], x)
                        return $name{T}(delete(elss, ii), Val{:safe}())
                    end
                end
            end
            return u
        end
        Base.propertynames(u::$name) = (:count, :elements)
        Base.getproperty(u::$name{T}, s::Symbol) where {T} = begin
            if s == :count
                return length(u)
            elseif s == :elements
                els = getfield(u, :elements)
                return els === nothing ? () : tuple(els...)
            else
                throw(ArgumentError("no such property $s of type $(typeof(u))"))
            end
        end
    end |> esc
end
# Macro-write the code fo these two types.
@_plinset_code PLinearSet isequal hash
@_plinset_code PIdLinearSet (===) objectid
# Document the types.
@doc """
    PLinearSet{T}

A persistent set type that uses a simple array representation internally;
accordingly `PLinearSet` does not peform persistent operations such as
`setindex`, `push`, and `delete` efficiently. It is intended to be used by Air
internally as a set that stores items whose hash values are identical.

See also: [`PIdLinearSet`](@ref), [`PSet`](@ref).
""" PLinearSet

@doc """
    PIdLinearSet{T}

A persistent `IdSet` type that uses a simple array representation internally;
accordingly `PIdLinearSet` does not peform persistent operations such as
`setindex`, `push`, and `delete` efficiently. It is intended to be used by Air
internally as a set that stores items whose hash values are identical.

See also: [`PLinearSet`](@ref), [`PSet`](@ref).
""" PIdLinearSet

# #PSet ========================================================================
# The PSet type employs both PTree and PLinearSet to create an unordered set.
macro _pset_code(name::Symbol, eq, h, linset)
    return quote
        struct $name{T} <: AbstractPSet{T}
            count::Int
            root::PTree{$linset{T}}
        end
        # Document the equal/hash types.
        equalfn(::Type{T}) where {T <: $name} = $eq
        hashfn(::Type{T}) where {T <: $name} = $h
        equalfn(::$name) = $eq
        hashfn(::$name) = $h
        # Construction.
        $name{T}() where {T} = $name{T}(0, PTree{$linset{T}}())
        $name{T}(itr) where {T} = reduce(push, itr, init=$name{T}())
        $name() = $name{Any}()
        $name(itr) = $name(itr, Base.IteratorEltype(itr))
        $name(itr, ::Base.HasEltype) = $name{eltype(itr)}(itr)
        $name(itr, ::Base.EltypeUnknown) = begin
            T = typejoin(map(typeof, itr)...)
            return $name{T}(itr)
        end
        # Base methods.
        Base.empty(s::$name{T}) where {T} = $name{T}()
        Base.empty(s::$name{T}, ::Type{S}) where {T,S} = $name{S}()
        Base.length(s::$name) = getfield(s, :count)
        Base.IteratorSize(::Type{$name{T}}) where {T} = Base.HasLength()
        Base.IteratorEltype(::Type{$name{T}}) where {T} = Base.HasEltype()
        Base.eltype(::Type{$name{T}}) where {T,N} = T
        Base.eltype(u::$name{T}) where {T,N} = T
        Base.iterate(u::$name{T}) where {T} = begin
            (getfield(u, :count) == 0) && return nothing
            (lst,titer) = iterate(getfield(u, :root))
            lst = lst[2]
            (el,liter) = iterate(lst)
            return (el,(lst,liter,titer))
        end
        Base.iterate(u::$name{T}, tup) where {T} = begin
            (lst,liter,titer) = tup
            if liter !== nothing
                r = iterate(lst, liter)
                if r !== nothing
                    (el,liter) = r
                    return (el,(lst,liter,titer))
                end
            end
            q = iterate(getfield(u, :root), titer)
            (q === nothing) && return q
            (lst,titer) = q
            lst = lst[2]
            (el,liter) = iterate(lst)
            return (el,(lst,liter,titer))
        end
        Base.in(x::S, u::$name{T}) where {T,S} = begin
            hh = $h(x)
            uu = get(getfield(u, :root), hh, nothing)
            return uu === nothing ? false : in(x, uu::$linset{T})
        end
        # And the Air methods.
        Air.push(u::$name{T}, x::S) where {T,S} = begin
            hh = $h(x)
            root = getfield(u, :root)
            uu = get(root, hh, nothing)
            if uu === nothing
                uu = $linset{T}(T[x], Val{:safe}())
            else
                in(x, uu) && return u
                uu = push(uu, x)
            end
            tr = setindex(root, uu, hh)
            return $name{T}(getfield(u, :count) + 1, tr)
        end
        Air.delete(u::$name{T}, x::S) where {T,S} = begin
            hh = $h(x)
            uu = get(getfield(u, :root), hh, nothing)
            (uu === nothing) && return u
            vv = delete(uu, x)
            (uu === vv) && return u
            if length(vv) == 0
                return $name{T}(getfield(u, :count) - 1, delete(getfield(u, :root), hh))
            else
                return $name{T}(getfield(u, :count) - 1, setindex(getfield(u, :root), vv, hh))
            end
        end
    end |> esc
end
# Declare the types.
@_pset_code PSet isequal hash PLinearSet
@_pset_code PIdSet (===) objectid PIdLinearSet
# The show function is best written in terms of an AbstractPSet.
Base.show(io::IO, s::AbstractPSet{T}) where {T} = begin
    S = typeof(s)
    print(io, get(io, :typeinfo, Any) == s ? String(nameof(S)) : S)
    print(io, "(")
    isempty(s) || Base.show_vector(io, s)
    print(io, ")")
end
# Document them.
@doc """
    PSet{T}

A persistent set type that uses an array hash-mapped trie structure to allow for
efficient updating of the set using the persistent operations defined by Air
such as `setindex`, `push`, and `delete`.

See also: `Set`, [`PIdSet`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> PSet()
PSet{Any}()

julia> :a in PSet([:a, :b, :c])
true

julia> :d in PSet([:a, :b, :c])
false

julia> :d in push(PSet(), :d)
true
```
""" PSet
@doc """
    PIdSet{T}

A persistent identity-based set type that uses an array hash-mapped trie
structure to allow for efficient updating of the set using the persistent
operations defined by Air such as `push`, and `delete`.

See also: `Set`, [`PIdSet`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> PIdSet()
PIdSet{Any}()

julia> :a in PIdSet([:a, :b, :c])
true

julia> :d in PIdSet([:a, :b, :c])
false

julia> :d in push(PIdSet(), :d)
true
```
""" PIdSet

# Export these types.
export PSet, PIdSet


text:
function reshop_solve(mdl::Ptr{reshop_model}, mdl_solver::Ptr{reshop_model}, ctx_dest::Ptr{context}, solver_name::String)

    CONFIG[:solver_log] && hack_solver_log()

    # Set the modeltype
    res = ccall((:myo_analyze_modeltype, libreshop), Cint, (Ptr{reshop_model}, Ptr{filter_ops}), mdl, C_NULL)
    res != 0 && error("return code $res from ReSHOP")

    # TODO(xhub) HIGH why is true here?
    res = ccall((:rhp_process, libreshop), Cint, (Ptr{reshop_model}, Ptr{reshop_model}), mdl, mdl_solver)
    res != 0 && error("return code $res from ReSHOP")

    if CONFIG[:export_gms]
        ccall((:gams_setsolverstr, libreshop), Cint, (Ptr{context}, Cstring), ctx_dest, "CONVERTD")
        ccall((:rhp_solve, libreshop), Cint, (Ptr{reshop_model},), mdl_solver)
    end

    ccall((:gams_setsolverstr, libreshop), Cint, (Ptr{ReSHOP.context}, Cstring), ctx_dest, solver_name)
    res = ccall((:rhp_solve, libreshop), Cint, (Ptr{reshop_model},), mdl_solver)

    return res
end

ext: jl


text:
@testset "1.2.4.2 (d x)^m (a x^q+b x^n+c x^(2 n-q))^p" begin
    (a, b, c, d, e, m, n, q, x, ) = @variables a b c d e m n q x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*problems*of*the*form*x^m*(a*x^q+b*x^n+c*x^(2*n-q))^p=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(a*x^2+b*x^3+c*x^4)^p=#

    #= ::Subsection::Closed:: =#
    #=x^m*(a*x^2+b*x^3+c*x^4)^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^2*(a*x^2 + b*x^3 + c*x^4), x, 2, (a*x^5)/5 + (b*x^6)/6 + (c*x^7)/7]
    @test_int [x*(a*x^2 + b*x^3 + c*x^4), x, 2, (a*x^4)/4 + (b*x^5)/5 + (c*x^6)/6]
    @test_int [a*x^2 + b*x^3 + c*x^4, x, 1, (a*x^3)/3 + (b*x^4)/4 + (c*x^5)/5]
    @test_int [(a*x^2 + b*x^3 + c*x^4)/x, x, 2, (a*x^2)/2 + (b*x^3)/3 + (c*x^4)/4]
    @test_int [(a*x^2 + b*x^3 + c*x^4)/x^2, x, 2, a*x + (b*x^2)/2 + (c*x^3)/3]

    @test_int [x^2*(a*x^2 + b*x^3 + c*x^4)^2, x, 3, (a^2*x^7)/7 + (a*b*x^8)/4 + ((b^2 + 2*a*c)*x^9)/9 + (b*c*x^10)/5 + (c^2*x^11)/11]
    @test_int [x*(a*x^2 + b*x^3 + c*x^4)^2, x, 3, (a^2*x^6)/6 + (2*a*b*x^7)/7 + ((b^2 + 2*a*c)*x^8)/8 + (2*b*c*x^9)/9 + (c^2*x^10)/10]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^2, x, 3, (a^2*x^5)/5 + (a*b*x^6)/3 + ((b^2 + 2*a*c)*x^7)/7 + (b*c*x^8)/4 + (c^2*x^9)/9]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^2/x, x, 3, (a^2*x^4)/4 + (2*a*b*x^5)/5 + ((b^2 + 2*a*c)*x^6)/6 + (2*b*c*x^7)/7 + (c^2*x^8)/8]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^2/x^2, x, 3, (a^2*x^3)/3 + (a*b*x^4)/2 + ((b^2 + 2*a*c)*x^5)/5 + (b*c*x^6)/3 + (c^2*x^7)/7]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^5/(a*x^2 + b*x^3 + c*x^4), x, 7, -((b*x)/c^2) + x^2/(2*c) + (b*(b^2 - 3*a*c)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(c^3*sqrt(b^2 - 4*a*c)) + ((b^2 - a*c)*log(a + b*x + c*x^2))/(2*c^3)]
    @test_int [x^4/(a*x^2 + b*x^3 + c*x^4), x, 6, x/c - ((b^2 - 2*a*c)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(c^2*sqrt(b^2 - 4*a*c)) - (b*log(a + b*x + c*x^2))/(2*c^2)]
    @test_int [x^3/(a*x^2 + b*x^3 + c*x^4), x, 5, (b*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(c*sqrt(b^2 - 4*a*c)) + log(a + b*x + c*x^2)/(2*c)]
    @test_int [x^2/(a*x^2 + b*x^3 + c*x^4), x, 3, (-2*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/sqrt(b^2 - 4*a*c)]
    @test_int [x/(a*x^2 + b*x^3 + c*x^4), x, 7, (b*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a*sqrt(b^2 - 4*a*c)) + log(x)/a - log(a + b*x + c*x^2)/(2*a)]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(-1), x, 8, -(1/(a*x)) - ((b^2 - 2*a*c)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^2*sqrt(b^2 - 4*a*c)) - (b*log(x))/a^2 + (b*log(a + b*x + c*x^2))/(2*a^2)]
    @test_int [1/(x*(a*x^2 + b*x^3 + c*x^4)), x, 8, -1/(2*a*x^2) + b/(a^2*x) + (b*(b^2 - 3*a*c)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^3*sqrt(b^2 - 4*a*c)) + ((b^2 - a*c)*log(x))/a^3 - ((b^2 - a*c)*log(a + b*x + c*x^2))/(2*a^3)]
    @test_int [1/(x^2*(a*x^2 + b*x^3 + c*x^4)), x, 8, -1/(3*a*x^3) + b/(2*a^2*x^2) - (b^2 - a*c)/(a^3*x) - ((b^4 - 4*a*b^2*c + 2*a^2*c^2)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^4*sqrt(b^2 - 4*a*c)) - (b*(b^2 - 2*a*c)*log(x))/a^4 + (b*(b^2 - 2*a*c)*log(a + b*x + c*x^2))/(2*a^4)]

    @test_int [x^8/(a*x^2 + b*x^3 + c*x^4)^2, x, 8, (2*(b^2 - 3*a*c)*x)/(c^2*(b^2 - 4*a*c)) - (b*x^2)/(c*(b^2 - 4*a*c)) + (x^3*(2*a + b*x))/((b^2 - 4*a*c)*(a + b*x + c*x^2)) - (2*(b^4 - 6*a*b^2*c + 6*a^2*c^2)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(c^3*(b^2 - 4*a*c)^(3/2)) - (b*log(a + b*x + c*x^2))/c^3]
    @test_int [x^7/(a*x^2 + b*x^3 + c*x^4)^2, x, 7, -((b*x)/(c*(b^2 - 4*a*c))) + (x^2*(2*a + b*x))/((b^2 - 4*a*c)*(a + b*x + c*x^2)) + (b*(b^2 - 6*a*c)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(c^2*(b^2 - 4*a*c)^(3/2)) + log(a + b*x + c*x^2)/(2*c^2)]
    @test_int [x^6/(a*x^2 + b*x^3 + c*x^4)^2, x, 4, (x*(2*a + b*x))/((b^2 - 4*a*c)*(a + b*x + c*x^2)) + (4*a*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(b^2 - 4*a*c)^(3/2)]
    @test_int [x^5/(a*x^2 + b*x^3 + c*x^4)^2, x, 4, (2*a + b*x)/((b^2 - 4*a*c)*(a + b*x + c*x^2)) - (2*b*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(b^2 - 4*a*c)^(3/2)]
    @test_int [x^4/(a*x^2 + b*x^3 + c*x^4)^2, x, 4, -((b + 2*c*x)/((b^2 - 4*a*c)*(a + b*x + c*x^2))) + (4*c*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(b^2 - 4*a*c)^(3/2)]
    @test_int [x^3/(a*x^2 + b*x^3 + c*x^4)^2, x, 8, (b^2 - 2*a*c + b*c*x)/(a*(b^2 - 4*a*c)*(a + b*x + c*x^2)) + (b*(b^2 - 6*a*c)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^2*(b^2 - 4*a*c)^(3/2)) + log(x)/a^2 - log(a + b*x + c*x^2)/(2*a^2)]
    @test_int [x^2/(a*x^2 + b*x^3 + c*x^4)^2, x, 8, (-2*(b^2 - 3*a*c))/(a^2*(b^2 - 4*a*c)*x) + (b^2 - 2*a*c + b*c*x)/(a*(b^2 - 4*a*c)*x*(a + b*x + c*x^2)) - (2*(b^4 - 6*a*b^2*c + 6*a^2*c^2)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^3*(b^2 - 4*a*c)^(3/2)) - (2*b*log(x))/a^3 + (b*log(a + b*x + c*x^2))/a^3]
    @test_int [x^1/(a*x^2 + b*x^3 + c*x^4)^2, x, 8, -((3*b^2 - 8*a*c)/(2*a^2*(b^2 - 4*a*c)*x^2)) + (b*(3*b^2 - 11*a*c))/(a^3*(b^2 - 4*a*c)*x) + (b^2 - 2*a*c + b*c*x)/(a*(b^2 - 4*a*c)*x^2*(a + b*x + c*x^2)) + (b*(3*b^4 - 20*a*b^2*c + 30*a^2*c^2)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^4*(b^2 - 4*a*c)^(3/2)) + ((3*b^2 - 2*a*c)*log(x))/a^4 - ((3*b^2 - 2*a*c)*log(a + b*x + c*x^2))/(2*a^4)]
    @test_int [x^0/(a*x^2 + b*x^3 + c*x^4)^2, x, 8, -((2*(2*b^2 - 5*a*c))/(3*a^2*(b^2 - 4*a*c)*x^3)) + (b*(2*b^2 - 7*a*c))/(a^3*(b^2 - 4*a*c)*x^2) - (2*(2*b^4 - 9*a*b^2*c + 5*a^2*c^2))/(a^4*(b^2 - 4*a*c)*x) + (b^2 - 2*a*c + b*c*x)/(a*(b^2 - 4*a*c)*x^3*(a + b*x + c*x^2)) - (2*(2*b^6 - 15*a*b^4*c + 30*a^2*b^2*c^2 - 10*a^3*c^3)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^5*(b^2 - 4*a*c)^(3/2)) - (2*b*(2*b^2 - 3*a*c)*log(x))/a^5 + (b*(2*b^2 - 3*a*c)*log(a + b*x + c*x^2))/a^5]
    @test_int [1/(x*(a*x^2 + b*x^3 + c*x^4)^2), x, 8, -((5*b^2 - 12*a*c)/(4*a^2*(b^2 - 4*a*c)*x^4)) + (b*(5*b^2 - 17*a*c))/(3*a^3*(b^2 - 4*a*c)*x^3) - (5*b^4 - 22*a*b^2*c + 12*a^2*c^2)/(2*a^4*(b^2 - 4*a*c)*x^2) + (b*(5*b^4 - 27*a*b^2*c + 29*a^2*c^2))/(a^5*(b^2 - 4*a*c)*x) + (b^2 - 2*a*c + b*c*x)/(a*(b^2 - 4*a*c)*x^4*(a + b*x + c*x^2)) + (b*(5*b^6 - 42*a*b^4*c + 105*a^2*b^2*c^2 - 70*a^3*c^3)*atanh((b + 2*c*x)/sqrt(b^2 - 4*a*c)))/(a^6*(b^2 - 4*a*c)^(3/2)) + ((5*b^4 - 12*a*b^2*c + 3*a^2*c^2)*log(x))/a^6 - ((5*b^4 - 12*a*b^2*c + 3*a^2*c^2)*log(a + b*x + c*x^2))/(2*a^6)]

    #= ::Subsection::Closed:: =#
    #=x^m*(a*x^2+b*x^3+c*x^4)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^2*sqrt(a*x^2 + b*x^3 + c*x^4), x, 8, (b*(35*b^2 - 116*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(960*c^3) - ((105*b^4 - 460*a*b^2*c + 256*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(1920*c^4*x) - ((7*b^2 - 16*a*c)*x*sqrt(a*x^2 + b*x^3 + c*x^4))/(240*c^2) + (x^2*(b + 8*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(40*c) + (b*(7*b^2 - 12*a*c)*(b^2 - 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(256*c^(9/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [x*sqrt(a*x^2 + b*x^3 + c*x^4), x, 7, -(((5*b^2 - 12*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(96*c^2)) + (b*(15*b^2 - 52*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(192*c^3*x) + (x*(b + 6*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(24*c) - ((b^2 - 4*a*c)*(5*b^2 - 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(128*c^(7/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [sqrt(a*x^2 + b*x^3 + c*x^4), x, 5, -((b*(b + 2*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(8*c^2*x)) + ((a + b*x + c*x^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(3*c*x) + (b*(b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(16*c^(5/2)*x*sqrt(a + b*x + c*x^2))]
    @test_int [sqrt(a*x^2 + b*x^3 + c*x^4)/x, x, 4, ((b + 2*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*c*x) - ((b^2 - 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(8*c^(3/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [sqrt(a*x^2 + b*x^3 + c*x^4)/x^2, x, 7, sqrt(a*x^2 + b*x^3 + c*x^4)/x - (sqrt(a)*x*sqrt(a + b*x + c*x^2)*atanh((2*a + b*x)/(2*sqrt(a)*sqrt(a + b*x + c*x^2))))/sqrt(a*x^2 + b*x^3 + c*x^4) + (b*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(2*sqrt(c)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [sqrt(a*x^2 + b*x^3 + c*x^4)/x^3, x, 7, -(sqrt(a*x^2 + b*x^3 + c*x^4)/x^2) - (b*x*sqrt(a + b*x + c*x^2)*atanh((2*a + b*x)/(2*sqrt(a)*sqrt(a + b*x + c*x^2))))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4)) + (sqrt(c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/sqrt(a*x^2 + b*x^3 + c*x^4)]
    @test_int [sqrt(a*x^2 + b*x^3 + c*x^4)/x^4, x, 5, -(sqrt(a*x^2 + b*x^3 + c*x^4)/(2*x^3)) - (b*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*a*x^2) + ((b^2 - 4*a*c)*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(8*a^(3/2))]
    @test_int [sqrt(a*x^2 + b*x^3 + c*x^4)/x^5, x, 6, -(sqrt(a*x^2 + b*x^3 + c*x^4)/(3*x^4)) - (b*sqrt(a*x^2 + b*x^3 + c*x^4))/(12*a*x^3) + ((3*b^2 - 8*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(24*a^2*x^2) - (b*(b^2 - 4*a*c)*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(16*a^(5/2))]
    @test_int [sqrt(a*x^2 + b*x^3 + c*x^4)/x^6, x, 7, -(sqrt(a*x^2 + b*x^3 + c*x^4)/(4*x^5)) - (b*sqrt(a*x^2 + b*x^3 + c*x^4))/(24*a*x^4) + ((5*b^2 - 12*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(96*a^2*x^3) - (b*(15*b^2 - 52*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(192*a^3*x^2) + ((b^2 - 4*a*c)*(5*b^2 - 4*a*c)*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(128*a^(7/2))]

    @test_int [x*(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 10, ((1155*b^6 - 8988*a*b^4*c + 18896*a^2*b^2*c^2 - 6720*a^3*c^3)*sqrt(a*x^2 + b*x^3 + c*x^4))/(286720*c^5) - (b*(3465*b^6 - 30660*a*b^4*c + 81648*a^2*b^2*c^2 - 58816*a^3*c^3)*sqrt(a*x^2 + b*x^3 + c*x^4))/(573440*c^6*x) - (b*(231*b^4 - 1560*a*b^2*c + 2416*a^2*c^2)*x*sqrt(a*x^2 + b*x^3 + c*x^4))/(71680*c^4) + ((99*b^4 - 568*a*b^2*c + 560*a^2*c^2)*x^2*sqrt(a*x^2 + b*x^3 + c*x^4))/(35840*c^3) - (x^3*(b*(11*b^2 + 68*a*c) + 10*c*(11*b^2 - 28*a*c)*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(4480*c^2) + (x*(3*b + 14*c*x)*(a*x^2 + b*x^3 + c*x^4)^(3/2))/(112*c) + (3*(b^2 - 4*a*c)^2*(33*b^4 - 72*a*b^2*c + 16*a^2*c^2)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(32768*c^(13/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 10, -((b*(105*b^4 - 728*a*b^2*c + 1168*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(17920*c^4)) + ((315*b^6 - 2520*a*b^4*c + 5488*a^2*b^2*c^2 - 2048*a^3*c^3)*sqrt(a*x^2 + b*x^3 + c*x^4))/(35840*c^5*x) + ((7*b^2 - 32*a*c)*(3*b^2 - 4*a*c)*x*sqrt(a*x^2 + b*x^3 + c*x^4))/(4480*c^3) - (b*(9*b^2 - 44*a*c)*x^2*sqrt(a*x^2 + b*x^3 + c*x^4))/(2240*c^2) + (x^3*(b^2 + 24*a*c + 10*b*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(280*c) + (1/7)*x*(a*x^2 + b*x^3 + c*x^4)^(3/2) - (3*b*(b^2 - 4*a*c)^2*(3*b^2 - 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(2048*c^(11/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x, x, 8, ((35*b^4 - 216*a*b^2*c + 240*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(3840*c^3) - (b*(105*b^4 - 760*a*b^2*c + 1296*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(7680*c^4*x) - (x*(b*(7*b^2 + 12*a*c) + 6*c*(7*b^2 - 20*a*c)*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(960*c^2) + ((3*b + 10*c*x)*(a*x^2 + b*x^3 + c*x^4)^(3/2))/(60*c*x) + ((b^2 - 4*a*c)^2*(7*b^2 - 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(1024*c^(9/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^2, x, 6, (3*b*(b^2 - 4*a*c)*(b + 2*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(128*c^3*x) - (b*(b + 2*c*x)*(a*x^2 + b*x^3 + c*x^4)^(3/2))/(16*c^2*x^3) + (a*x^2 + b*x^3 + c*x^4)^(5/2)/(5*c*x^5) - (3*b*(b^2 - 4*a*c)^2*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(256*c^(7/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^3, x, 5, -((3*(b^2 - 4*a*c)*(b + 2*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(64*c^2*x)) + ((b + 2*c*x)*(a*x^2 + b*x^3 + c*x^4)^(3/2))/(8*c*x^3) + (3*(b^2 - 4*a*c)^2*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(128*c^(5/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^4, x, 8, ((b^2 + 8*a*c + 2*b*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(8*c*x) + (a*x^2 + b*x^3 + c*x^4)^(3/2)/(3*x^3) - (a^(3/2)*x*sqrt(a + b*x + c*x^2)*atanh((2*a + b*x)/(2*sqrt(a)*sqrt(a + b*x + c*x^2))))/sqrt(a*x^2 + b*x^3 + c*x^4) - (b*(b^2 - 12*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(16*c^(3/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^5, x, 8, (3*(3*b + 2*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*x) - (a*x^2 + b*x^3 + c*x^4)^(3/2)/x^4 - (3*sqrt(a)*b*x*sqrt(a + b*x + c*x^2)*atanh((2*a + b*x)/(2*sqrt(a)*sqrt(a + b*x + c*x^2))))/(2*sqrt(a*x^2 + b*x^3 + c*x^4)) + (3*(b^2 + 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(8*sqrt(c)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^6, x, 8, -((3*(b - 2*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*x^2)) - (a*x^2 + b*x^3 + c*x^4)^(3/2)/(2*x^5) - (3*(b^2 + 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((2*a + b*x)/(2*sqrt(a)*sqrt(a + b*x + c*x^2))))/(8*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4)) + (3*b*sqrt(c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(2*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^7, x, 9, ((b^2 - 8*a*c + 2*b*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(8*a*x^2) - (a*x^2 + b*x^3 + c*x^4)^(3/2)/(3*x^6) - (b*(a*x^2 + b*x^3 + c*x^4)^(3/2))/(4*a*x^5) + (b*(b^2 - 12*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((2*a + b*x)/(2*sqrt(a)*sqrt(a + b*x + c*x^2))))/(16*a^(3/2)*sqrt(a*x^2 + b*x^3 + c*x^4)) + (c^(3/2)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/sqrt(a*x^2 + b*x^3 + c*x^4)]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^8, x, 7, -(((b^2 - 12*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(32*a*x^3)) + (b*(3*b^2 - 20*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(64*a^2*x^2) - ((b + 6*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(8*x^4) - (a*x^2 + b*x^3 + c*x^4)^(3/2)/(4*x^7) - (3*(b^2 - 4*a*c)^2*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(128*a^(5/2))]
    @test_int [(a*x^2 + b*x^3 + c*x^4)^(3/2)/x^9, x, 8, -(((b^2 - 8*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(80*a*x^4)) + (b*(5*b^2 - 28*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(320*a^2*x^3) - ((15*b^4 - 100*a*b^2*c + 128*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(640*a^3*x^2) - (3*(b + 4*c*x)*sqrt(a*x^2 + b*x^3 + c*x^4))/(40*x^5) - (a*x^2 + b*x^3 + c*x^4)^(3/2)/(5*x^8) + (3*b*(b^2 - 4*a*c)^2*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(256*a^(7/2))]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^3/sqrt(a*x^2 + b*x^3 + c*x^4), x, 6, -((3*b*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*c^2*x)) + sqrt(a*x^2 + b*x^3 + c*x^4)/(2*c) + ((3*b^2 - 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(8*c^(5/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [x^2/sqrt(a*x^2 + b*x^3 + c*x^4), x, 4, sqrt(a*x^2 + b*x^3 + c*x^4)/(c*x) - (b*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(2*c^(3/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [x/sqrt(a*x^2 + b*x^3 + c*x^4), x, 3, (x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(sqrt(c)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [1/sqrt(a*x^2 + b*x^3 + c*x^4), x, 2, -(atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4)))/sqrt(a))]
    @test_int [1/(x*sqrt(a*x^2 + b*x^3 + c*x^4)), x, 3, -(sqrt(a*x^2 + b*x^3 + c*x^4)/(a*x^2)) + (b*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(2*a^(3/2))]
    @test_int [1/(x^2*sqrt(a*x^2 + b*x^3 + c*x^4)), x, 5, -(sqrt(a*x^2 + b*x^3 + c*x^4)/(2*a*x^3)) + (3*b*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*a^2*x^2) - ((3*b^2 - 4*a*c)*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(8*a^(5/2))]

    @test_int [x^7/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 8, (2*x^4*(2*a + b*x))/((b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4)) + ((5*b^2 - 12*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(2*c^2*(b^2 - 4*a*c)) - (b*(15*b^2 - 52*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*c^3*(b^2 - 4*a*c)*x) - (2*b*x*sqrt(a*x^2 + b*x^3 + c*x^4))/(c*(b^2 - 4*a*c)) + (3*(5*b^2 - 4*a*c)*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(8*c^(7/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [x^6/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 7, (2*x^3*(2*a + b*x))/((b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4)) - (2*b*sqrt(a*x^2 + b*x^3 + c*x^4))/(c*(b^2 - 4*a*c)) + ((3*b^2 - 8*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(c^2*(b^2 - 4*a*c)*x) - (3*b*x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(2*c^(5/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [x^5/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 6, (2*x^2*(2*a + b*x))/((b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4)) - (2*b*sqrt(a*x^2 + b*x^3 + c*x^4))/(c*(b^2 - 4*a*c)*x) + (x*sqrt(a + b*x + c*x^2)*atanh((b + 2*c*x)/(2*sqrt(c)*sqrt(a + b*x + c*x^2))))/(c^(3/2)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [x^4/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 1, (2*x*(2*a + b*x))/((b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))]
    @test_int [x^3/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 1, -((2*x*(b + 2*c*x))/((b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4)))]
    @test_int [x^2/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 3, (2*x*(b^2 - 2*a*c + b*c*x))/(a*(b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4)) - atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4)))/a^(3/2)]
    @test_int [x/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 5, (2*(b^2 - 2*a*c + b*c*x))/(a*(b^2 - 4*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4)) - ((3*b^2 - 8*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(a^2*(b^2 - 4*a*c)*x^2) + (3*b*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(2*a^(5/2))]
    @test_int [1/(a*x^2 + b*x^3 + c*x^4)^(3/2), x, 6, (2*(b^2 - 2*a*c + b*c*x))/(a*(b^2 - 4*a*c)*x*sqrt(a*x^2 + b*x^3 + c*x^4)) - ((5*b^2 - 12*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(2*a^2*(b^2 - 4*a*c)*x^3) + (b*(15*b^2 - 52*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*a^3*(b^2 - 4*a*c)*x^2) - (3*(5*b^2 - 4*a*c)*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(8*a^(7/2))]
    @test_int [1/(x*(a*x^2 + b*x^3 + c*x^4)^(3/2)), x, 7, (2*(b^2 - 2*a*c + b*c*x))/(a*(b^2 - 4*a*c)*x^2*sqrt(a*x^2 + b*x^3 + c*x^4)) - ((7*b^2 - 16*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(3*a^2*(b^2 - 4*a*c)*x^4) + (b*(35*b^2 - 116*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(12*a^3*(b^2 - 4*a*c)*x^3) - ((105*b^4 - 460*a*b^2*c + 256*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(24*a^4*(b^2 - 4*a*c)*x^2) + (5*b*(7*b^2 - 12*a*c)*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(16*a^(9/2))]
    @test_int [1/(x^2*(a*x^2 + b*x^3 + c*x^4)^(3/2)), x, 8, (2*(b^2 - 2*a*c + b*c*x))/(a*(b^2 - 4*a*c)*x^3*sqrt(a*x^2 + b*x^3 + c*x^4)) - ((9*b^2 - 20*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(4*a^2*(b^2 - 4*a*c)*x^5) + (b*(21*b^2 - 68*a*c)*sqrt(a*x^2 + b*x^3 + c*x^4))/(8*a^3*(b^2 - 4*a*c)*x^4) - ((105*b^4 - 448*a*b^2*c + 240*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(32*a^4*(b^2 - 4*a*c)*x^3) + (b*(315*b^4 - 1680*a*b^2*c + 1808*a^2*c^2)*sqrt(a*x^2 + b*x^3 + c*x^4))/(64*a^5*(b^2 - 4*a*c)*x^2) - (15*(21*b^4 - 56*a*b^2*c + 16*a^2*c^2)*atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4))))/(128*a^(11/2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(a*x+b*x^3+c*x^5)^p=#

    #= ::Subsection::Closed:: =#
    #=x^m*(a*x+b*x^3+c*x^5)^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^m*(a*x + b*x^3 + c*x^5), x, 2, (a*x^(2 + m))/(2 + m) + (b*x^(4 + m))/(4 + m) + (c*x^(6 + m))/(6 + m)]

    @test_int [x^2*(a*x + b*x^3 + c*x^5), x, 2, (a*x^4)/4 + (b*x^6)/6 + (c*x^8)/8]
    @test_int [x^1*(a*x + b*x^3 + c*x^5), x, 2, (a*x^3)/3 + (b*x^5)/5 + (c*x^7)/7]
    @test_int [x^0*(a*x + b*x^3 + c*x^5), x, 1, (a*x^2)/2 + (b*x^4)/4 + (c*x^6)/6]
    @test_int [(a*x + b*x^3 + c*x^5)/x^1, x, 2, a*x + (b*x^3)/3 + (c*x^5)/5]
    @test_int [(a*x + b*x^3 + c*x^5)/x^2, x, 2, (b*x^2)/2 + (c*x^4)/4 + a*log(x)]
    @test_int [(a*x + b*x^3 + c*x^5)/x^3, x, 2, -(a/x) + b*x + (c*x^3)/3]

    @test_int [x^m*(a*x + b*x^3 + c*x^5)^2, x, 3, (a^2*x^(3 + m))/(3 + m) + (2*a*b*x^(5 + m))/(5 + m) + ((b^2 + 2*a*c)*x^(7 + m))/(7 + m) + (2*b*c*x^(9 + m))/(9 + m) + (c^2*x^(11 + m))/(11 + m)]

    @test_int [x^2*(a*x + b*x^3 + c*x^5)^2, x, 3, (a^2*x^5)/5 + (2*a*b*x^7)/7 + ((b^2 + 2*a*c)*x^9)/9 + (2*b*c*x^11)/11 + (c^2*x^13)/13]
    @test_int [x^1*(a*x + b*x^3 + c*x^5)^2, x, 4, (a^2*x^4)/4 + (a*b*x^6)/3 + ((b^2 + 2*a*c)*x^8)/8 + (b*c*x^10)/5 + (c^2*x^12)/12]
    @test_int [x^0*(a*x + b*x^3 + c*x^5)^2, x, 3, (a^2*x^3)/3 + (2*a*b*x^5)/5 + ((b^2 + 2*a*c)*x^7)/7 + (2*b*c*x^9)/9 + (c^2*x^11)/11]
    @test_int [(a*x + b*x^3 + c*x^5)^2/x^1, x, 4, (a^2*x^2)/2 + (a*b*x^4)/2 + ((b^2 + 2*a*c)*x^6)/6 + (b*c*x^8)/4 + (c^2*x^10)/10]
    @test_int [(a*x + b*x^3 + c*x^5)^2/x^2, x, 3, a^2*x + (2*a*b*x^3)/3 + ((b^2 + 2*a*c)*x^5)/5 + (2*b*c*x^7)/7 + (c^2*x^9)/9]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^8/(a*x + b*x^3 + c*x^5), x, 8, -((b*x^2)/(2*c^2)) + x^4/(4*c) + (b*(b^2 - 3*a*c)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*c^3*sqrt(b^2 - 4*a*c)) + ((b^2 - a*c)*log(a + b*x^2 + c*x^4))/(4*c^3)]
    @test_int [x^7/(a*x + b*x^3 + c*x^5), x, 6, -((b*x)/c^2) + x^3/(3*c) + ((b^2 - a*c - (b*(b^2 - 3*a*c))/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(sqrt(2)*c^(5/2)*sqrt(b - sqrt(b^2 - 4*a*c))) + ((b^2 - a*c + (b*(b^2 - 3*a*c))/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(sqrt(2)*c^(5/2)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [x^6/(a*x + b*x^3 + c*x^5), x, 7, x^2/(2*c) - ((b^2 - 2*a*c)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*c^2*sqrt(b^2 - 4*a*c)) - (b*log(a + b*x^2 + c*x^4))/(4*c^2)]
    @test_int [x^5/(a*x + b*x^3 + c*x^5), x, 5, x/c - ((b - (b^2 - 2*a*c)/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(sqrt(2)*c^(3/2)*sqrt(b - sqrt(b^2 - 4*a*c))) - ((b + (b^2 - 2*a*c)/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(sqrt(2)*c^(3/2)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [x^4/(a*x + b*x^3 + c*x^5), x, 6, (b*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*c*sqrt(b^2 - 4*a*c)) + log(a + b*x^2 + c*x^4)/(4*c)]
    @test_int [x^3/(a*x + b*x^3 + c*x^5), x, 4, -((sqrt(b - sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(sqrt(2)*sqrt(c)*sqrt(b^2 - 4*a*c))) + (sqrt(b + sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(sqrt(2)*sqrt(c)*sqrt(b^2 - 4*a*c))]
    @test_int [x^2/(a*x + b*x^3 + c*x^5), x, 4, -(atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c))/sqrt(b^2 - 4*a*c))]
    @test_int [x/(a*x + b*x^3 + c*x^5), x, 4, (sqrt(2)*sqrt(c)*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(sqrt(b^2 - 4*a*c)*sqrt(b - sqrt(b^2 - 4*a*c))) - (sqrt(2)*sqrt(c)*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(sqrt(b^2 - 4*a*c)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [1/(a*x + b*x^3 + c*x^5), x, 8, (b*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*a*sqrt(b^2 - 4*a*c)) + log(x)/a - log(a + b*x^2 + c*x^4)/(4*a)]
    @test_int [1/(x*(a*x + b*x^3 + c*x^5)), x, 5, -(1/(a*x)) - (sqrt(c)*(1 + b/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(sqrt(2)*a*sqrt(b - sqrt(b^2 - 4*a*c))) - (sqrt(c)*(1 - b/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(sqrt(2)*a*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [1/(x^2*(a*x + b*x^3 + c*x^5)), x, 9, -1/(2*a*x^2) - ((b^2 - 2*a*c)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*a^2*sqrt(b^2 - 4*a*c)) - (b*log(x))/a^2 + (b*log(a + b*x^2 + c*x^4))/(4*a^2)]

    @test_int [x^11/(a*x + b*x^3 + c*x^5)^2, x, 9, ((b^2 - 3*a*c)*x^2)/(c^2*(b^2 - 4*a*c)) - (b*x^4)/(2*c*(b^2 - 4*a*c)) + (x^6*(2*a + b*x^2))/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) - ((b^4 - 6*a*b^2*c + 6*a^2*c^2)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(c^3*(b^2 - 4*a*c)^(3/2)) - (b*log(a + b*x^2 + c*x^4))/(2*c^3)]
    @test_int [x^10/(a*x + b*x^3 + c*x^5)^2, x, 7, ((3*b^2 - 10*a*c)*x)/(2*c^2*(b^2 - 4*a*c)) - (b*x^3)/(2*c*(b^2 - 4*a*c)) + (x^5*(2*a + b*x^2))/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) - ((3*b^3 - 13*a*b*c - (3*b^4 - 19*a*b^2*c + 20*a^2*c^2)/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*c^(5/2)*(b^2 - 4*a*c)*sqrt(b - sqrt(b^2 - 4*a*c))) - ((3*b^3 - 13*a*b*c + (3*b^4 - 19*a*b^2*c + 20*a^2*c^2)/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*c^(5/2)*(b^2 - 4*a*c)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [x^9/(a*x + b*x^3 + c*x^5)^2, x, 8, -((b*x^2)/(2*c*(b^2 - 4*a*c))) + (x^4*(2*a + b*x^2))/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) + (b*(b^2 - 6*a*c)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*c^2*(b^2 - 4*a*c)^(3/2)) + log(a + b*x^2 + c*x^4)/(4*c^2)]
    @test_int [x^8/(a*x + b*x^3 + c*x^5)^2, x, 6, -((b*x)/(2*c*(b^2 - 4*a*c))) + (x^3*(2*a + b*x^2))/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) + ((b^2 - 6*a*c - (b*(b^2 - 8*a*c))/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*c^(3/2)*(b^2 - 4*a*c)*sqrt(b - sqrt(b^2 - 4*a*c))) + ((b^2 - 6*a*c + (b*(b^2 - 8*a*c))/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*c^(3/2)*(b^2 - 4*a*c)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [x^7/(a*x + b*x^3 + c*x^5)^2, x, 5, (x^2*(2*a + b*x^2))/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) + (2*a*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(b^2 - 4*a*c)^(3/2)]
    @test_int [x^6/(a*x + b*x^3 + c*x^5)^2, x, 5, (x*(2*a + b*x^2))/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) + ((b - (b^2 + 4*a*c)/sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*sqrt(c)*(b^2 - 4*a*c)*sqrt(b - sqrt(b^2 - 4*a*c))) + ((b^2 + 4*a*c + b*sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*sqrt(c)*(b^2 - 4*a*c)^(3/2)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [x^5/(a*x + b*x^3 + c*x^5)^2, x, 5, (2*a + b*x^2)/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) - (b*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(b^2 - 4*a*c)^(3/2)]
    @test_int [x^4/(a*x + b*x^3 + c*x^5)^2, x, 5, -((x*(b + 2*c*x^2))/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4))) + (sqrt(c)*(2*b - sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(sqrt(2)*(b^2 - 4*a*c)^(3/2)*sqrt(b - sqrt(b^2 - 4*a*c))) - (sqrt(c)*(2*b + sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(sqrt(2)*(b^2 - 4*a*c)^(3/2)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [x^3/(a*x + b*x^3 + c*x^5)^2, x, 5, -((b + 2*c*x^2)/(2*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4))) + (2*c*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(b^2 - 4*a*c)^(3/2)]
    @test_int [x^2/(a*x + b*x^3 + c*x^5)^2, x, 5, (x*(b^2 - 2*a*c + b*c*x^2))/(2*a*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) + (sqrt(c)*(b^2 - 12*a*c + b*sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*a*(b^2 - 4*a*c)^(3/2)*sqrt(b - sqrt(b^2 - 4*a*c))) - (sqrt(c)*(b^2 - 12*a*c - b*sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*a*(b^2 - 4*a*c)^(3/2)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [x^1/(a*x + b*x^3 + c*x^5)^2, x, 9, (b^2 - 2*a*c + b*c*x^2)/(2*a*(b^2 - 4*a*c)*(a + b*x^2 + c*x^4)) + (b*(b^2 - 6*a*c)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*a^2*(b^2 - 4*a*c)^(3/2)) + log(x)/a^2 - log(a + b*x^2 + c*x^4)/(4*a^2)]
    @test_int [x^0/(a*x + b*x^3 + c*x^5)^2, x, 6, -((3*b^2 - 10*a*c)/(2*a^2*(b^2 - 4*a*c)*x)) + (b^2 - 2*a*c + b*c*x^2)/(2*a*(b^2 - 4*a*c)*x*(a + b*x^2 + c*x^4)) - (sqrt(c)*(3*b^3 - 16*a*b*c + (3*b^2 - 10*a*c)*sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*a^2*(b^2 - 4*a*c)^(3/2)*sqrt(b - sqrt(b^2 - 4*a*c))) + (sqrt(c)*(3*b^3 - 16*a*b*c - (3*b^2 - 10*a*c)*sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*a^2*(b^2 - 4*a*c)^(3/2)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [1/(x^1*(a*x + b*x^3 + c*x^5)^2), x, 9, -((b^2 - 3*a*c)/(a^2*(b^2 - 4*a*c)*x^2)) + (b^2 - 2*a*c + b*c*x^2)/(2*a*(b^2 - 4*a*c)*x^2*(a + b*x^2 + c*x^4)) - ((b^4 - 6*a*b^2*c + 6*a^2*c^2)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(a^3*(b^2 - 4*a*c)^(3/2)) - (2*b*log(x))/a^3 + (b*log(a + b*x^2 + c*x^4))/(2*a^3)]
    @test_int [1/(x^2*(a*x + b*x^3 + c*x^5)^2), x, 7, -((5*b^2 - 14*a*c)/(6*a^2*(b^2 - 4*a*c)*x^3)) + (b*(5*b^2 - 19*a*c))/(2*a^3*(b^2 - 4*a*c)*x) + (b^2 - 2*a*c + b*c*x^2)/(2*a*(b^2 - 4*a*c)*x^3*(a + b*x^2 + c*x^4)) + (sqrt(c)*(5*b^4 - 29*a*b^2*c + 28*a^2*c^2 + b*(5*b^2 - 19*a*c)*sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b - sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*a^3*(b^2 - 4*a*c)^(3/2)*sqrt(b - sqrt(b^2 - 4*a*c))) - (sqrt(c)*(5*b^4 - 29*a*b^2*c + 28*a^2*c^2 - b*(5*b^2 - 19*a*c)*sqrt(b^2 - 4*a*c))*atan((sqrt(2)*sqrt(c)*x)/sqrt(b + sqrt(b^2 - 4*a*c))))/(2*sqrt(2)*a^3*(b^2 - 4*a*c)^(3/2)*sqrt(b + sqrt(b^2 - 4*a*c)))]
    @test_int [1/(x^3*(a*x + b*x^3 + c*x^5)^2), x, 9, -((3*b^2 - 8*a*c)/(4*a^2*(b^2 - 4*a*c)*x^4)) + (b*(3*b^2 - 11*a*c))/(2*a^3*(b^2 - 4*a*c)*x^2) + (b^2 - 2*a*c + b*c*x^2)/(2*a*(b^2 - 4*a*c)*x^4*(a + b*x^2 + c*x^4)) + (b*(3*b^4 - 20*a*b^2*c + 30*a^2*c^2)*atanh((b + 2*c*x^2)/sqrt(b^2 - 4*a*c)))/(2*a^4*(b^2 - 4*a*c)^(3/2)) + ((3*b^2 - 2*a*c)*log(x))/a^4 - ((3*b^2 - 2*a*c)*log(a + b*x^2 + c*x^4))/(4*a^4)]

    #= ::Subsection::Closed:: =#
    #=x^m*(a*x+b*x^3+c*x^5)^(p/2)=#

    @test_int [x/sqrt(a*x + b*x^3 + c*x^5), x, 3, (2*x^2*sqrt(1 + (2*c*x^2)/(b - sqrt(b^2 - 4*a*c)))*sqrt(1 + (2*c*x^2)/(b + sqrt(b^2 - 4*a*c)))*AppellF1(3/4, 1/2, 1/2, 7/4, -((2*c*x^2)/(b - sqrt(b^2 - 4*a*c))), -((2*c*x^2)/(b + sqrt(b^2 - 4*a*c)))))/(3*sqrt(a*x + b*x^3 + c*x^5))]

    #= ::Subsection::Closed:: =#
    #=x^(m/2)*(a*x+b*x^3+c*x^5)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [x^(3/2)*sqrt(a*x + b*x^3 + c*x^5), x, 5, -((2*(b^2 - 3*a*c)*x^(3/2)*(a + b*x^2 + c*x^4))/(15*c^(3/2)*(sqrt(a) + sqrt(c)*x^2)*sqrt(a*x + b*x^3 + c*x^5))) + (sqrt(x)*(b + 3*c*x^2)*sqrt(a*x + b*x^3 + c*x^5))/(15*c) + (2*a^(1/4)*(b^2 - 3*a*c)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.E(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(15*c^(7/4)*sqrt(a*x + b*x^3 + c*x^5)) - (a^(1/4)*(2*b^2 + sqrt(a)*b*sqrt(c) - 6*a*c)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(30*c^(7/4)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [sqrt(x)*sqrt(a*x + b*x^3 + c*x^5), x, 5, ((b + 2*c*x^2)*sqrt(a*x + b*x^3 + c*x^5))/(8*c*sqrt(x)) - ((b^2 - 4*a*c)*sqrt(x)*sqrt(a + b*x^2 + c*x^4)*atanh((b + 2*c*x^2)/(2*sqrt(c)*sqrt(a + b*x^2 + c*x^4))))/(16*c^(3/2)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [sqrt(a*x + b*x^3 + c*x^5)/sqrt(x), x, 5, (b*x^(3/2)*(a + b*x^2 + c*x^4))/(3*sqrt(c)*(sqrt(a) + sqrt(c)*x^2)*sqrt(a*x + b*x^3 + c*x^5)) + (1/3)*sqrt(x)*sqrt(a*x + b*x^3 + c*x^5) - (a^(1/4)*b*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.E(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(3*c^(3/4)*sqrt(a*x + b*x^3 + c*x^5)) + (a^(1/4)*(b + 2*sqrt(a)*sqrt(c))*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(6*c^(3/4)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [sqrt(a*x + b*x^3 + c*x^5)/x^(3/2), x, 8, sqrt(a*x + b*x^3 + c*x^5)/(2*sqrt(x)) - (sqrt(a)*sqrt(x)*sqrt(a + b*x^2 + c*x^4)*atanh((2*a + b*x^2)/(2*sqrt(a)*sqrt(a + b*x^2 + c*x^4))))/(2*sqrt(a*x + b*x^3 + c*x^5)) + (b*sqrt(x)*sqrt(a + b*x^2 + c*x^4)*atanh((b + 2*c*x^2)/(2*sqrt(c)*sqrt(a + b*x^2 + c*x^4))))/(4*sqrt(c)*sqrt(a*x + b*x^3 + c*x^5))]

    @test_int [x^(3/2)*(a*x + b*x^3 + c*x^5)^(3/2), x, 8, ((15*b^4 - 100*a*b^2*c + 128*a^2*c^2)*sqrt(a*x + b*x^3 + c*x^5))/(1280*c^3*sqrt(x)) - (x^(3/2)*(b*(5*b^2 - 4*a*c) + 4*c*(5*b^2 - 16*a*c)*x^2)*sqrt(a*x + b*x^3 + c*x^5))/(640*c^2) + (sqrt(x)*(3*b + 8*c*x^2)*(a*x + b*x^3 + c*x^5)^(3/2))/(80*c) - (3*b*(b^2 - 4*a*c)^2*sqrt(x)*sqrt(a + b*x^2 + c*x^4)*atanh((b + 2*c*x^2)/(2*sqrt(c)*sqrt(a + b*x^2 + c*x^4))))/(512*c^(7/2)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [sqrt(x)*(a*x + b*x^3 + c*x^5)^(3/2), x, 6, ((8*b^4 - 57*a*b^2*c + 84*a^2*c^2)*x^(3/2)*(a + b*x^2 + c*x^4))/(315*c^(5/2)*(sqrt(a) + sqrt(c)*x^2)*sqrt(a*x + b*x^3 + c*x^5)) - (sqrt(x)*(b*(4*b^2 - 9*a*c) + 6*c*(2*b^2 - 7*a*c)*x^2)*sqrt(a*x + b*x^3 + c*x^5))/(315*c^2) + ((3*b + 7*c*x^2)*(a*x + b*x^3 + c*x^5)^(3/2))/(63*c*sqrt(x)) - (a^(1/4)*(8*b^4 - 57*a*b^2*c + 84*a^2*c^2)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.E(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(315*c^(11/4)*sqrt(a*x + b*x^3 + c*x^5)) + (a^(1/4)*(8*b^4 - 57*a*b^2*c + 84*a^2*c^2 + 4*sqrt(a)*b*sqrt(c)*(b^2 - 6*a*c))*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(630*c^(11/4)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [(a*x + b*x^3 + c*x^5)^(3/2)/sqrt(x), x, 6, -((3*(b^2 - 4*a*c)*(b + 2*c*x^2)*sqrt(a*x + b*x^3 + c*x^5))/(128*c^2*sqrt(x))) + ((b + 2*c*x^2)*(a*x + b*x^3 + c*x^5)^(3/2))/(16*c*x^(3/2)) + (3*(b^2 - 4*a*c)^2*sqrt(x)*sqrt(a + b*x^2 + c*x^4)*atanh((b + 2*c*x^2)/(2*sqrt(c)*sqrt(a + b*x^2 + c*x^4))))/(256*c^(5/2)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [(a*x + b*x^3 + c*x^5)^(3/2)/x^(3/2), x, 6, -((2*b*(b^2 - 8*a*c)*x^(3/2)*(a + b*x^2 + c*x^4))/(35*c^(3/2)*(sqrt(a) + sqrt(c)*x^2)*sqrt(a*x + b*x^3 + c*x^5))) + (sqrt(x)*(b^2 + 10*a*c + 3*b*c*x^2)*sqrt(a*x + b*x^3 + c*x^5))/(35*c) + (a*x + b*x^3 + c*x^5)^(3/2)/(7*sqrt(x)) + (2*a^(1/4)*b*(b^2 - 8*a*c)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.E(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(35*c^(7/4)*sqrt(a*x + b*x^3 + c*x^5)) - (a^(1/4)*(sqrt(a)*sqrt(c)*(b^2 - 20*a*c) + 2*b*(b^2 - 8*a*c))*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(70*c^(7/4)*sqrt(a*x + b*x^3 + c*x^5))]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [x^(3/2)/sqrt(a*x + b*x^3 + c*x^5), x, 4, (sqrt(x)*sqrt(a + b*x^2 + c*x^4)*atanh((b + 2*c*x^2)/(2*sqrt(c)*sqrt(a + b*x^2 + c*x^4))))/(2*sqrt(c)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [sqrt(x)/sqrt(a*x + b*x^3 + c*x^5), x, 2, (sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(2*a^(1/4)*c^(1/4)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [1/(sqrt(x)*sqrt(a*x + b*x^3 + c*x^5)), x, 2, -(atanh((sqrt(x)*(2*a + b*x^2))/(2*sqrt(a)*sqrt(a*x + b*x^3 + c*x^5)))/(2*sqrt(a)))]
    @test_int [1/(x^(3/2)*sqrt(a*x + b*x^3 + c*x^5)), x, 6, (sqrt(c)*x^(3/2)*(a + b*x^2 + c*x^4))/(a*(sqrt(a) + sqrt(c)*x^2)*sqrt(a*x + b*x^3 + c*x^5)) - sqrt(a*x + b*x^3 + c*x^5)/(a*x^(3/2)) - (c^(1/4)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.E(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(a^(3/4)*sqrt(a*x + b*x^3 + c*x^5)) + (c^(1/4)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(2*a^(3/4)*sqrt(a*x + b*x^3 + c*x^5))]

    @test_int [x^(3/2)/(a*x + b*x^3 + c*x^5)^(3/2), x, 5, (x^(3/2)*(b^2 - 2*a*c + b*c*x^2))/(a*(b^2 - 4*a*c)*sqrt(a*x + b*x^3 + c*x^5)) - (b*sqrt(c)*x^(3/2)*(a + b*x^2 + c*x^4))/(a*(b^2 - 4*a*c)*(sqrt(a) + sqrt(c)*x^2)*sqrt(a*x + b*x^3 + c*x^5)) + (b*c^(1/4)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.E(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(a^(3/4)*(b^2 - 4*a*c)*sqrt(a*x + b*x^3 + c*x^5)) - (c^(1/4)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(2*a^(3/4)*(b - 2*sqrt(a)*sqrt(c))*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [sqrt(x)/(a*x + b*x^3 + c*x^5)^(3/2), x, 3, (sqrt(x)*(b^2 - 2*a*c + b*c*x^2))/(a*(b^2 - 4*a*c)*sqrt(a*x + b*x^3 + c*x^5)) - atanh((sqrt(x)*(2*a + b*x^2))/(2*sqrt(a)*sqrt(a*x + b*x^3 + c*x^5)))/(2*a^(3/2))]
    @test_int [1/(sqrt(x)*(a*x + b*x^3 + c*x^5)^(3/2)), x, 6, (b^2 - 2*a*c + b*c*x^2)/(a*(b^2 - 4*a*c)*sqrt(x)*sqrt(a*x + b*x^3 + c*x^5)) + (2*sqrt(c)*(b^2 - 3*a*c)*x^(3/2)*(a + b*x^2 + c*x^4))/(a^2*(b^2 - 4*a*c)*(sqrt(a) + sqrt(c)*x^2)*sqrt(a*x + b*x^3 + c*x^5)) - (2*(b^2 - 3*a*c)*sqrt(a*x + b*x^3 + c*x^5))/(a^2*(b^2 - 4*a*c)*x^(3/2)) - (2*c^(1/4)*(b^2 - 3*a*c)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.E(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(a^(7/4)*(b^2 - 4*a*c)*sqrt(a*x + b*x^3 + c*x^5)) + (c^(1/4)*(2*b^2 + sqrt(a)*b*sqrt(c) - 6*a*c)*sqrt(x)*(sqrt(a) + sqrt(c)*x^2)*sqrt((a + b*x^2 + c*x^4)/(sqrt(a) + sqrt(c)*x^2)^2)*Elliptic.F(2*atan((c^(1/4)*x)/a^(1/4)), (1/4)*(2 - b/(sqrt(a)*sqrt(c)))))/(2*a^(7/4)*(b^2 - 4*a*c)*sqrt(a*x + b*x^3 + c*x^5))]
    @test_int [1/(x^(3/2)*(a*x + b*x^3 + c*x^5)^(3/2)), x, 5, (b^2 - 2*a*c + b*c*x^2)/(a*(b^2 - 4*a*c)*x^(3/2)*sqrt(a*x + b*x^3 + c*x^5)) - ((3*b^2 - 8*a*c)*sqrt(a*x + b*x^3 + c*x^5))/(2*a^2*(b^2 - 4*a*c)*x^(5/2)) + (3*b*atanh((sqrt(x)*(2*a + b*x^2))/(2*sqrt(a)*sqrt(a*x + b*x^3 + c*x^5))))/(4*a^(5/2))]

    @test_int [x^((3*(n - 1))/2)/(a*x^(n - 1) + b*x^n + c*x^(n + 1))^(3/2), x, 1, -((2*x^((1/2)*(-1 + n))*(b + 2*c*x))/((b^2 - 4*a*c)*sqrt(a*x^(-1 + n) + b*x^n + c*x^(1 + n))))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(d+e*x^2)^n*(a*x+b*x^3+c*x^5)^p=#

    #= ::Subsection:: =#
    #=x^m*(d+e*x^2)^n*(a*x+b*x^3+c*x^5)^p=#

    #= ::Subsection::Closed:: =#
    #=x^m*(d+e*x^2)^n*(a*x+b*x^3+c*x^5)^(p/2)=#

    @test_int [x*(d + e*x^2)/sqrt(a*x + b*x^3 + c*x^5), x, 7, (2*d*x^2*sqrt(1 + (2*c*x^2)/(b - sqrt(b^2 - 4*a*c)))*sqrt(1 + (2*c*x^2)/(b + sqrt(b^2 - 4*a*c)))*AppellF1(3/4, 1/2, 1/2, 7/4, -((2*c*x^2)/(b - sqrt(b^2 - 4*a*c))), -((2*c*x^2)/(b + sqrt(b^2 - 4*a*c)))))/(3*sqrt(a*x + b*x^3 + c*x^5)) + (2*e*x^4*sqrt(1 + (2*c*x^2)/(b - sqrt(b^2 - 4*a*c)))*sqrt(1 + (2*c*x^2)/(b + sqrt(b^2 - 4*a*c)))*AppellF1(7/4, 1/2, 1/2, 11/4, -((2*c*x^2)/(b - sqrt(b^2 - 4*a*c))), -((2*c*x^2)/(b + sqrt(b^2 - 4*a*c)))))/(7*sqrt(a*x + b*x^3 + c*x^5))]

    #= ::Subsection:: =#
    #=x^(m/2)*(d+e*x^2)^n*(a*x+b*x^3+c*x^5)^(p/2)=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*(a*x^q+b*x^n+c*x^(2*n-q))^p*when*m=q/2-1=#

    @test_int [1/sqrt(3*x^2 - 3*x^4 + x^6), x, 2, -(atanh((x*(6 - 3*x^2))/(2*sqrt(3)*sqrt(3*x^2 - 3*x^4 + x^6)))/(2*sqrt(3)))]
    @test_int [1/sqrt(x^2*(3 - 3*x^2 + x^4)), x, 3, -(atanh((x*(6 - 3*x^2))/(2*sqrt(3)*sqrt(3*x^2 - 3*x^4 + x^6)))/(2*sqrt(3)))]
    @test_int [1/sqrt(1 - (1 - x^2)^3), x, 3, -(atanh((x*(6 - 3*x^2))/(2*sqrt(3)*sqrt(3*x^2 - 3*x^4 + x^6)))/(2*sqrt(3)))]

    @test_int [sqrt(3*x^2 - 3*x^4 + x^6), x, 5, -(((3 - 2*x^2)*sqrt(3*x^2 - 3*x^4 + x^6))/(8*x)) - (3*sqrt(3*x^2 - 3*x^4 + x^6)*asinh((3 - 2*x^2)/sqrt(3)))/(16*x*sqrt(3 - 3*x^2 + x^4))]
    @test_int [sqrt(x^2*(3 - 3*x^2 + x^4)), x, 6, -(((3 - 2*x^2)*sqrt(3*x^2 - 3*x^4 + x^6))/(8*x)) - (3*sqrt(3*x^2 - 3*x^4 + x^6)*asinh((3 - 2*x^2)/sqrt(3)))/(16*x*sqrt(3 - 3*x^2 + x^4))]
    @test_int [sqrt(1 - (1 - x^2)^3), x, 6, -(((3 - 2*x^2)*sqrt(3*x^2 - 3*x^4 + x^6))/(8*x)) - (3*sqrt(3*x^2 - 3*x^4 + x^6)*asinh((3 - 2*x^2)/sqrt(3)))/(16*x*sqrt(3 - 3*x^2 + x^4))]

    @test_int [1/(x*sqrt(a + b*x + c*x^2)), x, 2, -(atanh((2*a + b*x)/(2*sqrt(a)*sqrt(a + b*x + c*x^2)))/sqrt(a))]
    @test_int [1/sqrt(x^2*(a + b*x + c*x^2)), x, 3, -(atanh((x*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^2 + b*x^3 + c*x^4)))/sqrt(a))]
    @test_int [1/(sqrt(x)*sqrt(x*(a + b*x + c*x^2))), x, 3, -(atanh((sqrt(x)*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x + b*x^2 + c*x^3)))/sqrt(a))]
    @test_int [sqrt(x)/sqrt(x^3*(a + b*x + c*x^2)), x, 3, -(atanh((x^(3/2)*(2*a + b*x))/(2*sqrt(a)*sqrt(a*x^3 + b*x^4 + c*x^5)))/sqrt(a))]

    @test_int [1/(x*sqrt(a + b*x^2 + c*x^4)), x, 3, -(atanh((2*a + b*x^2)/(2*sqrt(a)*sqrt(a + b*x^2 + c*x^4)))/(2*sqrt(a)))]
    @test_int [1/sqrt(x^2*(a + b*x^2 + c*x^4)), x, 3, -(atanh((x*(2*a + b*x^2))/(2*sqrt(a)*sqrt(a*x^2 + b*x^4 + c*x^6)))/(2*sqrt(a)))]
    @test_int [1/(sqrt(x)*sqrt(x*(a + b*x^2 + c*x^4))), x, 3, -(atanh((sqrt(x)*(2*a + b*x^2))/(2*sqrt(a)*sqrt(a*x + b*x^3 + c*x^5)))/(2*sqrt(a)))]
    @test_int [sqrt(x)/(sqrt(x^3*(a + b*x^2 + c*x^4))), x, 3, -(atanh((x^(3/2)*(2*a + b*x^2))/(2*sqrt(a)*sqrt(a*x^3 + b*x^5 + c*x^7)))/(2*sqrt(a)))]

    @test_int [1/(x*sqrt(3 - 3*x^2 + x^4)), x, 3, -(atanh((sqrt(3)*(2 - x^2))/(2*sqrt(3 - 3*x^2 + x^4)))/(2*sqrt(3)))]
    @test_int [1/sqrt(x^2*(3 - 3*x^2 + x^4)), x, 3, -(atanh((x*(6 - 3*x^2))/(2*sqrt(3)*sqrt(3*x^2 - 3*x^4 + x^6)))/(2*sqrt(3)))]
    @test_int [1/(sqrt(x)*sqrt(x*(3 - 3*x + x^2))), x, 3, -(atanh((sqrt(3)*(2 - x)*sqrt(x))/(2*sqrt(3*x - 3*x^2 + x^3)))/sqrt(3))]

    @test_int [x^(q/2 - 1)/sqrt(a*x^q + b*x^n + c*x^(2*n - q)), x, 2, -(atanh((x^(q/2)*(2*a + b*x^(n - q)))/(2*sqrt(a)*sqrt(b*x^n + c*x^(2*n - q) + a*x^q)))/(sqrt(a)*(n - q)))]
    end

max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.4 Improper/1.2.4.2 (d x)^m (a x^q+b x^n+c x^(2 n-q))^p.jl


text:
# Definitely inspired by jraph: https://github.com/deepmind/jraph/blob/master/jraph/_src/graph.py

abstract type AbstractGraphTuple end

# Keep structs immutable to avoid weird "side-effects"
struct GraphTuple <: AbstractGraphTuple
    nodes
    edges
    globals
    senders::AbstractVector{<:Integer}
    receivers::AbstractVector{<:Integer}
end

function GraphTuple(;nodes=nothing, edges=nothing, globals=nothing,
                    senders, receivers)

    # Sanity Check
    @assert length(senders) == length(receivers)

    return GraphTuple(nodes, edges, globals, senders, receivers)
end

# quick question: do we need accessor methods?
nodes(g::AbstractGraphTuple) = g.nodes
edges(g::AbstractGraphTuple) = g.edges
globals(g::AbstractGraphTuple) = g.globals

senders(g::AbstractGraphTuple) = g.senders
receivers(g::AbstractGraphTuple) = g.receivers

# TODO: does this fully make sense?
nnodes(g::AbstractGraphTuple) = size(nodes(g), 2)
nedges(g::AbstractGraphTuple) = g |> senders |> length

updatenodes(g::AbstractGraphTuple, v) = @set g.nodes = v
updateedges(g::AbstractGraphTuple, v) = @set g.edges = v
updateglobals(g::AbstractGraphTuple, v) = @set g.globals = v

struct BatchGraphTuple <: AbstractGraphTuple
    nodes
    edges
    globals
    senders::AbstractVector{<:Integer}
    receivers::AbstractVector{<:Integer}

    node2graph::AbstractVector{<:Integer}
    edge2graph::AbstractVector{<:Integer}
    numgraphs::Integer
end

# this is one ugly function, please refactor ASAP
function batchgraphs(gs::Vector{GraphTuple}, nodespergraph, edgespergraph, symmetrize=false)
    nodesv = hcat([nodes(g) for g in gs]...)
    edgesv = hcat([edges(g) for g in gs]...)
    globalsv = hcat([globals(g) for g in gs]...)

    numgraphs = length(gs)
    nodes2graph = vcat(fill.(1:numgraphs, nodespergraph)...)
    edges2graph = vcat(fill.(1:numgraphs, edgespergraph)...)

    shift = circshift(cumsum(nodespergraph), 1)
    shift[1] = 0
    edgeshift = vcat(fill.(shift, edgespergraph)...)

    sendersv = vcat([senders(g) for g in gs]...) .+ edgeshift
    receiversv = vcat([receivers(g) for g in gs]...) .+ edgeshift

    if symmetrize
        sendersvsym = vcat([sendersv, receiversv]...)
        receiversvsym = vcat([receiversv, sendersv]...)
        sendersv, receiversv = sendersvsym, receiversvsym

        edges2graph = repeat(edges2graph, 2)
        edgesv = repeat(edgesv, 1, 2)
    end

    return BatchGraphTuple(nodesv, edgesv, globalsv,
                            sendersv, receiversv,
                            nodes2graph, edges2graph, numgraphs)
end

Flux.@functor GraphTuple
Flux.@functor BatchGraphTuple


text:
@testset "example notebooks" begin
    notebookdir = joinpath(@__DIR__, "..", "notebooks")
    excludes = String[]
    printinterval = 60 # seconds
    printcallback = timer -> Compat.@info "Running notebook tests."
    timertask = if VERSION < v"0.7-"
        Timer(printcallback, printinterval, printinterval)
    else
        Timer(printcallback, printinterval, interval=printinterval)
    end
    for file in readdir(notebookdir)
        path = joinpath(notebookdir, file)
        path in excludes && continue
        name, ext = splitext(file)
        lowercase(ext) == ".ipynb" || continue

        @eval module $(gensym()) # Each notebook is run in its own module.
        using Compat
        using Compat.Test
        using NBInclude
        @testset "$($name)" begin
            @nbinclude($path; regex = r"^((?!\#NBSKIP).)*$"s) # Use #NBSKIP in a cell to skip it during tests.
        end
        end # module
    end
    close(timertask)
end


text:
# =========================================================================================
@doc raw"""
    _etol(::Float64)
    _etol(::Type{Float64})
    _etol(::Float32)
    _etol(::Type{Float32})

Common tolerance values internally used by NESSie.

# Return type
`T`
""" _etol
const _etol_f64 = 1.45e-8
const _etol_f32 = 3.45f-4

@inline _etol(::Float64) = _etol_f64
@inline _etol(::Type{Float64}) = _etol_f64
@inline _etol(::Float32) = _etol_f32
@inline _etol(::Type{Float32}) = _etol_f32

# =========================================================================================
"""
Vacuum permittivity

# Unit
``\\frac{F}{m}``
"""
const ε0 = 1 / (4π * 1e-7 * 299792458^2)

# =========================================================================================
"""
Geometric quantity ``σ(ξ)``

```math
σ(ξ) = \\lim_{ε→0} \\frac{1}{4πε²} ∫_{r ∈ Ω:|r-ξ|=ε} dΓᵣ = \\frac{1}{2}
```
for almost all ``ξ ∈ Γ`` (cf. [[Ste03]](@ref Bibliography)).
"""
const σ  = 0.5

# =========================================================================================
"""
``10^{10}`` times the elementary charge (for ``Å → m`` conversion)

# Unit
``C``
"""
const ec = 1.602176e-9

# =========================================================================================
@doc raw"""
    potprefactor(T::Type{Float64} = Float64)
    potprefactor(T::Type{Float32})

Common prefactor for all potentials ``φ_Ω`` and ``φ_Σ``:

```math
\frac{1.602 ⋅ 10^{-19}}{10^{-10} ⋅ 4π  ⋅ ε₀} ≈ 1.145 ⋅ 4π
```

# Return type
`T`
""" potprefactor
for T in [:Float64, :Float32]
    varname = Symbol("potprefactor_", T)
    @eval begin
        const $(varname) = $(T)(ec / 4π / ε0)
        potprefactor(::Type{$(T)}) = $(varname)
    end
end
potprefactor() = potprefactor(Float64)

# =========================================================================================
"""
    struct Option{T <: AbstractFloat}
        εΩ::T       # dielectric constant of the solute
        εΣ::T       # dielectric constant of the solvent
        ε∞::T       # large-scale (bulk) solvent response
        λ ::T       # correlation length scale [λ] = Å
    end

System parameters
"""
mutable struct Option{T <: AbstractFloat}
    "dielectric constant of the solute"
    εΩ::T
    "dielectric constant of the solvent"
    εΣ::T
    "large-scale (bulk) solvent response"
    ε∞::T
    "correlation length scale [λ] = Å"
    λ::T
end

# =========================================================================================
@doc raw"""
    defaultopt(T::Type{Float64} = Float64)
    defaultopt(T::Type{Float32})

Default system parameters

# Return type
[`Option{T}`](@ref NESSie.Option)

# Default values
 * ``ε_Ω = 2``
 * ``ε_Σ = 78``
 * ``ε_∞ = 1.8``
 * ``λ = 20 Å``
""" defaultopt
for T in [:Float64, :Float32]
    varname = Symbol("defaultopt_", T)
    @eval begin
        const $(varname) = Option($(T)[2, 78, 1.8, 20]...)
        defaultopt(::Type{$(T)}) = $(varname)
    end
end
defaultopt() = defaultopt(Float64)

# =========================================================================================
"""
    yukawa{T}(opt::Option{T})

Exponent ``1/Λ`` for the fundamental solution of the yukawa operator

```math
Λ := λ\\sqrt{\\frac{ε_∞}{ε_Σ}}
```

# Return type
`T`
"""
yukawa(opt::Option{T}) where T = √(opt.εΣ/opt.ε∞)/opt.λ


text:
# test codes for
#   * matrix vec conversions

X = rand(3, 3)
Y = vecs(X) |> mat
@test X == Y
Y[1] = -1.0
@test X == Y

# the positions test is also a test for the py-reference thing

# check the positions-as-ref thing
at = Atoms("Al") * (2,3,2)
X = unsafe_positions(at)
X[1] = JVec(rand(3))
@test X == positions(at)

# test set_positions!
Y = copy(X)
Y[3] = JVec(rand(3))
@test Y != positions(at)
set_positions!(at, Y)
@test Y == positions(at)


text:
# This file is a part of JuliaFEM/TetGen.
# License is MIT: see https://github.com/JuliaFEM/TetGen.jl/blob/master/LICENSE.md

using TetGen
using FactCheck
using Logging
@Logging.configure(level=DEBUG)

facts("Testing if somebody used print, println(), @sprint in src directory") do
  # TODO: make better reqular expression. Currently it will match all print words
  lines_with_print = Dict()
  src_dir = joinpath(Pkg.dir("TetGen"),"src")
  src = readdir(src_dir)
  for file_name in src
    fil = open(joinpath(src_dir,file_name),"r")
    for (line_number,line) in enumerate(readlines(fil))
      if ismatch(r"print",line)
        lines_with_print[file_name * ":$line_number"] = "print found in line: $line_number"
      end
    end
    close(fil)
  end
  @fact lines_with_print => isempty "Instead of println() use Logging.jl package"
end

facts("Testing if we have non ascii characters in the src files") do
  # TODO: we should allow Greeck letters in documentation. Thus this test should skip docstrings
  lines_with_non_ascii = []
  src_dir = joinpath(Pkg.dir("TetGen"),"src")
  src = readdir(src_dir)
  for file_name in src
    fil = open(joinpath(src_dir,file_name),"r")
    for (line_number,line) in enumerate(readlines(fil))
      if ismatch(r"[^\x00-\x7F]",line)
        push!(lines_with_non_ascii, file_name * ":line $line_number")
      end
    end
    close(fil)
  end
  @fact lines_with_non_ascii => isempty "non ascii charecters found in src -> test is failing"
end

facts("Looking the [src,test] folders *.jl files header information") do
  files_no_license = []
  pkg_dir = Pkg.dir("TetGen")
  dirs = ["test";"src";"dev"]
  for folder in dirs
    for file_name in readdir(joinpath(pkg_dir,folder))
      if file_name[end-2:end] == ".jl"
        fil = open(joinpath(pkg_dir,folder,file_name),"r")
        head = readall(fil)
      else
        continue
      end
      if ~ismatch(r"This file is a part of JuliaFEM/TetGen",head)
        push!(files_no_license,"$file_name is missing JuliaFEM/TetGen statement")
      end
      if ~ismatch(r"MIT",head)
        push!(files_no_license,"$file_name is missing MIT statement")
      end
      if ~ismatch(r"https://github.com/JuliaFEM/TetGen.jl/blob/master/LICENSE.md",head)
        push!(files_no_license,"$file_name is missing reference to LICENSE.md")
      end
      close(fil)
    end
  end
  out_str = """License information missing or incorrect. Please use these two lines in the beginning of each file:

# This file is a part of JuliaFEM/TetGen.
# License is MIT: see https://github.com/JuliaFEM/TetGen.jl/blob/master/LICENSE.md

  """
  @fact files_no_license => isempty out_str
end


text:
# Momentum-oriented technical indicator functions

"""
```
aroon(hl::Matrix{T}; n::Int64=25)::Array{Float64}
```

Aroon up/down/oscillator

*Output*

- Column 1: Aroon Up
- Column 2: Aroon Down
- Column 3: Aroon Oscillator
"""
function aroon(hl::AbstractMatrix{T}; n::Int64=25)::Matrix{Float64} where {T<:Real}
    @assert size(hl,2) == 2 "Argument `hl` must have exactly 2 columns."
    @assert n < size(hl,1) "Argument `n` must be less than the number of rows in argument `hl`."
    out = zeros(T, (size(hl,1),3))
    @inbounds for i in n:size(hl,1)
        out[i,1] = 100.0 * (n - findmax(hl[i-n+1:i,1])[2]) / 25.0
    end
    @inbounds for i in n:size(hl,1)
        out[i,2] = 100.0 * (n - findmin(hl[i-n+1:i,2])[2]) / 25.0
    end
    out[:,3] = out[:,1]-out[:,2]  # aroon oscillator
    out[1:n-1,:] .= NaN
    return out
end

"""
```
donch(hl::Matrix{T}; n::Int64=10, inclusive::Bool=true)::Array{Float64}
```

Donchian channel (if inclusive is set to true, will include current bar in calculations.)

*Output*

- Column 1: Lowest low of last `n` periods
- Column 2: Average of highest high and lowest low of last `n` periods
- Column 3: Highest high of last `n` periods
"""
function donch(hl::AbstractMatrix{T}; n::Int64=10, inclusive::Bool=true)::Matrix{Float64} where {T<:Real}
    @assert size(hl,2) == 2 "Argument `hl` must have exactly 2 columns."
    local lower::Array{T} = runmin(hl[:,2], n=n, cumulative=false, inclusive=inclusive)
    local upper::Array{T} = runmax(hl[:,1], n=n, cumulative=false, inclusive=inclusive)
    local middle::Array{T} = (lower .+ upper) ./ 2.0
    return [lower middle upper]
end

"""
```
momentum(x::Array{T}; n::Int64=1)::Array{Float64}
```

Momentum indicator (price now vs price `n` periods back)
"""
function momentum(x::AbstractArray{T}; n::Int64=1)::Array{Float64} where {T<:Real}
    @assert n>0 "Argument n must be positive."
    return diffn(x, n=n)
end

"""
```
roc(x::Array{T}; n::Int64=1)::Array{Float64}
```

Rate of change indicator (percent change between i'th observation and (i-n)'th observation)
"""
function roc(x::AbstractArray{T}; n::Int64=1)::Array{Float64} where {T<:Real}
    @assert n<size(x,1) && n>0 "Argument n out of bounds."
    out = zeros(size(x)) .* NaN
    @inbounds for i = (n+1):size(x,1)
        out[i] = x[i]/x[i-n] - 1.0
    end
    return out * 100.0
end

"""
```
macd(x::Array{T}; nfast::Int64=12, nslow::Int64=26, nsig::Int64=9)::Array{Float64}
```

Moving average convergence-divergence

*Output*

- Column 1: MACD
- Column 2: MACD Signal Line
- Column 3: MACD Histogram
"""
function macd(x::AbstractArray{T}; nfast::Int64=12, nslow::Int64=26, nsig::Int64=9,
    fastMA::Function=ema, slowMA::Function=ema, signalMA::Function=sma)::Matrix{Float64} where {T<:Real}
    out = zeros(T, (length(x),3))
    out[:,1] = fastMA(x, n=nfast) - slowMA(x, n=nslow)
    out[:,2] = signalMA(out[:,1], n=nsig)
    out[:,3] = out[:,1] - out[:,2]
    return out
end

"""
```
rsi(x::Array{T}; n::Int64=14, ma::Function=ema, args...)::Array{Float64}
```

Relative strength index
"""
function rsi(x::AbstractArray{T}; n::Int64=14, ma::Function=ema, args...)::Array{Float64} where {T<:Real}
    @assert n<size(x,1) && n>0 "Argument n is out of bounds."
    N = size(x,1)
    ups = zeros(N)
    dns = zeros(N)
    zro = 0.0
    dx = [NaN; ndims(x) > 1 ? diff(x, dims=1) : diff(x)]
    @inbounds for i=2:N
        if dx[i] > zro
            ups[i] = dx[i]
        elseif dx[i] < zro
            dns[i] = -dx[i]
        end
    end
    rs = [NaN; ma(ups[2:end], n=n; args...) ./ ma(dns[2:end], n=n; args...)]
    return 100.0 .- 100.0 ./ (1.0 .+ rs)
end

"""
```
adx(hlc::Array{T}; n::Int64=14, wilder=true)::Array{Float64}
```

Average directional index

*Output*

- Column 1: DI+
- Column 2: DI-
- Column 3: ADX
"""
function adx(hlc::AbstractArray{T}; n::Int64=14, ma::Function=ema, args...)::Matrix{Float64} where {T<:Real}
    @assert n<size(hlc,1) && n>0 "Argument n is out of bounds."
    if size(hlc,2) != 3
        error("HLC array must have three columns")
    end
    N = size(hlc,1)
    updm = zeros(N)
    dndm = zeros(N)
    updm[1] = dndm[1] = NaN
    @inbounds for i = 2:N
        upmove = hlc[i,1] - hlc[i-1,1]
        dnmove = hlc[i-1,2] - hlc[i,2]
        if upmove > dnmove && upmove > 0.0
            updm[i] = upmove
        elseif dnmove > upmove && dnmove > 0.0
            dndm[i] = dnmove
        end
    end
    dip = [NaN; ma(updm[2:N], n=n; args...)] ./ atr(hlc, n=n) * 100.0
    dim = [NaN; ma(dndm[2:N], n=n; args...)] ./ atr(hlc, n=n) * 100.0
    dmx = abs.(dip-dim) ./ (dip+dim)
    adx = [fill(NaN,n); ma(dmx[n+1:N], n=n; args...)] * 100.0
    return [dip dim adx]
end

"""
```
psar(hl::Array{T}; af_min::T=0.02, af_max::T=0.2, af_inc::T=af_min)::Array{Float64}
```

Parabolic stop and reverse (SAR)

*Arguments*
- `hl`: 2D array of high and low prices in first and second columns respectively
- `af_min`: starting/initial value for acceleration factor
- `af_max`: maximum acceleration factor (accel factor capped at this value)
- `af_inc`: increment to the acceleration factor (speed of increase in accel factor)
"""
function psar(hl::AbstractArray{T}; af_min::T=0.02, af_max::T=0.2, af_inc::T=af_min)::Array{Float64} where {T<:Real}
    @assert af_min<1.0 && af_min>0.0 "Argument af_min must be in [0,1]."
    @assert af_max<1.0 && af_max>0.0 "Argument af_max must be in [0,1]."
    @assert af_inc<1.0 && af_inc>0.0 "Argument af_inc must be in [0,1]."
    @assert size(hl,2) == 2 "Argument hl must have 2 columns."
    ls0 = 1
    ls = 0
    af0 = af_min
    af = 0.0
    ep0 = hl[1,1]
    ep = 0.0
    maxi = 0.0
    mini = 0.0
    sar = zeros(T,size(hl,1))
    sar[1] = hl[1,2] - std(hl[:,1]-hl[:,2])
    @inbounds for i = 2:size(hl,1)
        ls = ls0
        ep = ep0
        af = af0
        mini = min(hl[i-1,2], hl[i,2])
        maxi = max(hl[i-1,1], hl[i,1])
        # Long/short signals and local extrema
        if (ls == 1)
            ls0 = hl[i,2] > sar[i-1] ? 1 : -1
            ep0 = max(maxi, ep)
        else
            ls0 = hl[i,1] < sar[i-1] ? -1 : 1
            ep0 = min(mini, ep)
        end
        # Acceleration vector
        if ls0 == ls  # no signal change
            sar[i] = sar[i-1] + af*(ep-sar[i-1])
            af0 = (af == af_max) ? af_max : (af + af_inc)
            if ls0 == 1  # current long signal
                af0 = (ep0 > ep) ? af0 : af
                sar[i] = min(sar[i], mini)
            else  # current short signal
                af0 = (ep0 < ep) ? af0 : af
                sar[i] = max(sar[i], maxi)
            end
        else  # new signal
            af0 = af_min
            sar[i] = ep0
        end
    end
    return sar
end

"""
```
kst(x::Array{T};
    nroc::Array{Int64}=[10,15,20,30], navg::Array{Int64}=[10,10,10,15],
    wgts::Array{Int64}=collect(1:length(nroc)), ma::Function=sma)::Array{Float64}

```

KST (Know Sure Thing) -- smoothed and summed rates of change
"""
function kst(x::AbstractArray{T}; nroc::AbstractArray{Int64}=[10,15,20,30], navg::AbstractArray{Int64}=[10,10,10,15],
    wgts::AbstractArray{Int64}=collect(1:length(nroc)), ma::Function=sma)::Array{Float64} where {T<:Real}
    @assert length(nroc) == length(navg)
    @assert all(nroc.>0) && all(nroc.<size(x,1))
    @assert all(navg.>0) && all(navg.<size(x,1))
    N = length(x)
    k = length(nroc)
    out = zeros(size(x))
    @inbounds for j = 1:k
        out += ma(roc(x, n=nroc[j]), n=navg[j]) * wgts[j]
    end
    return out
end

"""
```
wpr(hlc::Matrix{T}, n::Int64=14)::Array{Float64}
```

Williams %R
"""
function wpr(hlc::AbstractMatrix{T}; n::Int64=14)::Array{Float64} where {T<:Real}
    hihi = runmax(hlc[:,1], n=n, cumulative=false)
    lolo = runmin(hlc[:,2], n=n, cumulative=false)
    return -100 * (hihi - hlc[:,3]) ./ (hihi - lolo)
end

"""
```
cci(hlc::Matrix{T}; n::Int64=20, c::T=0.015, ma::Function=sma)::Array{Float64}
```

Commodity channel index
"""
function cci(hlc::AbstractMatrix{T}; n::Int64=20, c::T=0.015, ma::Function=sma, args...)::Array{Float64} where {T<:Real}
    tp = (hlc[:,1] + hlc[:,2] + hlc[:,3]) ./ 3.0
    dev = runmad(tp, n=n, cumulative=false, fun=mean)
    avg = ma(tp, n=n; args...)
    return (tp - avg) ./ (c * dev)
end

"""
```
stoch(hlc::Matrix{T}; nK::Int64=14, nD::Int64=3, kind::Symbol=:fast, ma::Function=sma, args...)::Matrix{Float64}
```

Stochastic oscillator (fast or slow)
"""
function stoch(hlc::AbstractMatrix{T}; nK::Int64=14, nD::Int64=3,
    kind::Symbol=:fast, ma::Function=sma, args...)::Matrix{Float64} where {T<:Real}
    @assert kind == :fast || kind == :slow "Argument `kind` must be either :fast or :slow"
    @assert nK<size(hlc,1) && nK>0 "Argument `nK` out of bounds."
    @assert nD<size(hlc,1) && nD>0 "Argument `nD` out of bounds."
    hihi = runmax(hlc[:,1], n=nK, cumulative=false)
    lolo = runmin(hlc[:,2], n=nK, cumulative=false)
    out = zeros(T, (size(hlc,1),2))
    out[:,1] = (hlc[:,3]-lolo) ./ (hihi-lolo) * 100.0
    out[:,2] = ma(out[:,1], n=nD; args...)
    if kind == :slow
        out[:,1] = out[:,2]
        out[:,2] = ma(out[:,1], n=nD; args...)
    end
    return out
end

"""
```
smi(hlc::Matrix{T}; n::Int64=13, nFast::Int64=2, nSlow::Int64=25, nSig::Int64=9,
    maFast::Function=ema, maSlow::Function=ema, maSig::Function=sma)::Matrix{Float64}
```

SMI (stochastic momentum oscillator)
"""
function smi(hlc::AbstractMatrix{T}; n::Int64=13, nFast::Int64=2, nSlow::Int64=25, nSig::Int64=9,
    maFast::Function=ema, maSlow::Function=ema, maSig::Function=sma)::Matrix{Float64} where {T<:Real}
    hihi = runmax(hlc[:,1], n=n, cumulative=false)
    lolo = runmin(hlc[:,2], n=n, cumulative=false)
    hldif = hihi-lolo
    delta = hlc[:,3] - (hihi+lolo) / 2.0
    numer = maSlow(maFast(delta, n=nFast), n=nSlow)
    denom = maSlow(maFast(hldif, n=nFast), n=nSlow) / 2.0
    out = zeros(T, (size(hlc,1),2))
    out[:,1] = 100.0*(numer./denom)
    out[:,2] = maSig(out[:,1], n=nSig)
    return out
end


text:
using NeighbourLists
using Test

# ---- FLAGS -----

# whether to run performance tests
performance = true

# check whether on CI
isCI = haskey(ENV, "CI")
notCI = !isCI

# TODO: switch the JuLIP test to an ASE test
# check whether we have JuLIP
# hasjulip = true
# try
#    using JuLIP
# catch
#    hasjulip = false
# end

# ----------------- TESTS -------------------

println("# threads = $(Base.Threads.nthreads())")

@testset "NeighbourLists" begin
   @testset "Aux" begin include("test_aux.jl") end
   @testset "CellList" begin include("test_celllist.jl") end

   # pointless until we switch to comparing against ASE / matscipy
   # if hasjulip
   #    @testset "JuLIP" begin include("test_julip.jl") end
   # end
end


text:
"""
   Implementing the ideas in

      @ARTICLE{1993JCoPh.104..211D,
          author = {{Dang-Vu}, H. and {Delcarte}, C.},
           title = "{An Accurate Solution of the Poisson Equation by the Chebyshev Collocation Method}",
         journal = {Journal of Computational Physics},
            year = 1993,
           month = jan,
          volume = {104},
          number = {1},
           pages = {211-220},
             doi = {10.1006/jcph.1993.1021},
          adsurl = {https://ui.adsabs.harvard.edu/abs/1993JCoPh.104..211D},
         adsnote = {Provided by the SAO/NASA Astrophysics Data System}
      }
"""
module ChebyshevSpec

include("CustomTypes.jl")
using .CustomTypes
using SparseArrays

export mat_invD1, mat_invD2

function mat_invD1(nx::myI)
   X = Vector{myF}()
   Y = Vector{myF}()
   V = Vector{myF}()

   c(n) = (n==0)      ? 2.0 : 1.0
   e(n) = (n<=(nx-1)) ? 1.0 : 0.0

   for i=2:(nx-1)
      n=i-1
      push!(X,i)
      push!(Y,i-1)
      push!(V,tomyF(c(n-1)*0.5/n))

      push!(X,i)
      push!(Y,i+1)
      push!(V,tomyF(-e(n+1)*0.5/n))
   end

   return dropzeros(sparse(X,Y,V))
end

function mat_invD2(nx::myI)
   X = Vector{myF}()
   Y = Vector{myF}()
   V = Vector{myF}()

   c(n) = (n==0)      ? 2.0 : 1.0
   e(n) = (n<=(nx-1)) ? 1.0 : 0.0

   for i=3:(nx-2)
      n=i-1
      push!(X,i)
      push!(Y,i-2)
      push!(V,tomyF(c(n-2)*0.25/(n*(n-1))))

      push!(X,i)
      push!(Y,i)
      push!(V,tomyF(-e(n+2)*0.5/((n+1)*(n-1))))

      push!(X,i)
      push!(Y,i+2)
      push!(V,tomyF(e(n+4)*0.25/(n*(n+1))))
   end
   n=nx-1
   push!(X,nx-1)
   push!(Y,nx-3)
   push!(V,tomyF(c(n-2)*0.25/(n*(n-1))))

   n=nx
   push!(X,nx)
   push!(Y,nx-2)
   push!(V,tomyF(c(n-2)*0.25/(n*(n-1))))

   return dropzeros(sparse(X,Y,V))
end

end


text:
@testset "RationalPoly" begin
    #@test RationalPoly{true, Int, Int}(1) == 1
    #@test typeof(RationalPoly{false, Int, Int}(1)) == RationalPoly{false, Int, Int}
    #@inferred RationalPoly{true, Int, Int}(1)
    Mod.@polyvar x
    @test 1 / x isa RationalPoly
    @test 1.0 / x isa RationalPoly
    @inferred 1 / x
    @inferred 1.0 / x
    @test eltype([2.0x / x^2, (x+x) / (1 + 2x^2)]) <: RationalPoly
    @test 2 * (1/x * (1-x)) + (1/x * x) * (1/x^2 * x^2) - (1-x)/x == (1-x)/x + 1
    @test (1/x + 1/x) / 2 == ((1 / (x^2 - 1) + (x+1)) - (x+1)) * ((x^2 - 1) / x)
    #@test typeof(zero(1/x)) == Term{true, Int}
    @test iszero(zero(1/x))
    @test zero(1/x) == 0
    @test one(1/x) == 1
    #@test typeof(zero(RationalPoly{true, Float64, Int})) == Polynomial{true, Float64}
    #@test iszero(zero(RationalPoly{true, Float64, Int}))
    #@test typeof(x / 2) == Term{true, Float64}
    #@test typeof((x + x^2) / 3.0) == Polynomial{true, Float64}
    @test nothing != x / (x^2 + 1)
    @test (x^2 + 1) / (2x) != nothing
    @test Dict{Int,Int}() != x / (x^2 + 1)
    @test (x^2 + 1) / (2x) != Dict{Int,Int}()
    @test numerator(x / x^2) == x
    @test denominator(x / x^2) == x^2
    @test inv(x / x^2) == x
    @test x / x^2 == inv(x)
    @test isone(((x+1) / (x-1)) / ((x+1) / (x-1)))
    @test ((x+1)^2 / (x-1)) / ((x+1) / (x-1)) == x+1
    poly = (x+1)/(x+2.0)
    RType = typeof(poly)
    @test RType(true) == one(RType)
    @test RType(false) == zero(RType)
    @test one(RType) isa RType
    @test zero(RType) isa RType
    @test one(poly) isa RType
    @test zero(poly) isa RType
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/solving-difference-equation-part-2/67057

# %%
using DynamicalSystems, DelimitedFiles, .Threads,  BenchmarkTools, ProgressMeter
## Components of a test DiscreteDynamicalSystem

function dds_constructor(u0 = [0.5, 0.7];  r=1.0, k=2.0)
    return DiscreteDynamicalSystem(dds_rule, u0, [r, k], dds_jac)
end
## equations of motion:
function dds_rule(x, par, n)
   r, k = par
   a, mu, d = 5.0, 0.5, 0.2
   dx = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
   dy = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
   return @SVector [dx, dy]
end
## Jacobian:
function dds_jac(x, par, n)
    r, k = par;
    a = 5.0; mu = 0.5; d = 0.2;
    J11 = exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)) - x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1))*(r/k - (2*x[1]*x[2])/(x[1]^2 + a)^2)
    J12 = -(x[1]*exp(- x[2]/(x[1]^2 + a) - r*(x[1]/k - 1)))/(x[1]^2 + a)
    J21 = x[2]*exp((mu*x[1])/(x[1]^2 + a) - d)*(mu/(x[1]^2 + a) - (2*mu*x[1]^2)/(x[1]^2 + a)^2)
    J22 = exp((mu*x[1])/(x[1]^2 + a) - d)
    return @SMatrix [J11  J12; J21 J22]
end

# %%
function traj!(tr, f, p, x0, T; Ttr=0)
    x = x0
    for i in 1:Ttr
        x = f(x, p, i-1)
    end
    @inbounds tr[1] = x
    for i in 1:T
        @inbounds tr[i+1] = f(tr[i], p, Ttr+i-1)
    end
    tr
end

# %%
function traj2d!(tr1, tr2, f, p, x0, T; Ttr=0)
    x = x0
    for i in 1:Ttr
        x = f(x, p, i-1)
    end
    @inbounds tr1[1], tr2[1] = x
    for i in 1:T
        @inbounds tr1[i+1], tr2[i+1] = f((tr1[i], tr2[i]), p, Ttr+i-1)
    end
    tr1, tr2
end

# %%
traj2d!(zeros(21), zeros(21), ((a, b), p, t) -> (b, a+b), nothing, (0, 1), 20)

# %%
function seqper_new(x; tol=1e-3)       # function to calculate periodicity
    n = length(x)
    @inbounds for k in 2:(n ÷ 2 + 1)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:n) && return k - 1
        end
    end
    return n
end

# %%
function meshh(x,y)    # create a meshgrid of x and y
    len_x = length(x);
    len_y = length(y);
    xmesh = Float64[]
    ymesh = Float64[]
    for i in 1:len_x
        for j in 1:len_y
            push!(xmesh, x[i])
            push!(ymesh, y[j])
        end
    end
    return  xmesh, ymesh
end

x = 1:3
y = 1:4
@show meshh(x, y)
meshh(x, y) .== vec.(reim(complex.(x', y)))

# %%
function meshh!(xmesh, ymesh, x, y)    # create a meshgrid of x and y
    len_x, len_y = length(x), length(y)
    for i in 1:len_x
        for j in 1:len_y
            xmesh[len_y*(i-1) + j] = x[i]
            ymesh[len_y*(i-1) + j] = y[j]
        end
    end
end

x = 1:3
y = 1:4
xmesh = Vector{Float64}(undef, length(x)*length(y))
ymesh = Vector{Float64}(undef, length(x)*length(y))
meshh!(xmesh, ymesh, x, y)
@show xmesh, ymesh
meshh(x, y) .== (xmesh, ymesh)

# %%
function isoperiodic_test_org(ds, par_area::Matrix{Float64}, NIter::Int64,
                         nxblock::Int64, nyblock::Int64, NTr::Int64, xpts::Int64, ypts::Int64)

    p1_st =  par_area[1];          # beginning of first parameter
    p1_nd =  par_area[2];          # end of first parameter
    p2_st =  par_area[3];          # beginning of second parameter
    p2_nd =  par_area[4];          # end of second parameter

    p1_block = range(p1_st, p1_nd, length=nxblock + 1);    # divide par_area in blocks
    p2_block = range(p2_st, p2_nd, length=nyblock + 1);    # divide par_area in blocks
    l_bipar = xpts * ypts;                                 # total number of pts in each block
    total_blocks = nxblock * nyblock                       # total number of blocks
    sol_last = Array{Float64,2}(undef, NIter + 1, l_bipar);   # stores x values of solution for l_bipar (r,k) paris

    dss = [deepcopy(ds) for _ in 1:nthreads()]
    periods = Vector{Int}(undef, l_bipar * total_blocks)
    par1mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    par2mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    prog = Progress(total_blocks)
    for ii in 1:nxblock
        par1range = range(p1_block[ii], p1_block[ii+1], length=xpts+1)[1:end-1]
        for jj in 1:nyblock
            par2range = range(p2_block[jj], p2_block[jj+1], length=ypts+1)[1:end-1]
            Threads.@threads for i = 1:xpts
                tid = threadid()
                for j = 1:ypts
                     set_parameter!(dss[tid], [par1range[i], par2range[j]]) # change the parameter values
                     tr = trajectory(dss[tid], NIter; Ttr=NTr)              # find the solution
                     sol_last[:, ypts*(i-1)+j] = tr[:,1]                    # store x-values
                end
            end
            a = l_bipar*(nyblock*(ii-1)+jj-1)
            periods[a+1:a+l_bipar] = seqper_new.(eachcol(sol_last), tol=0.001) # seqper_new calculates periodicity
            par1mesh[a+1:a+l_bipar], par2mesh[a+1:a+l_bipar] = meshh(par1range, par2range) # meshgrid of parameter values
            next!(prog)
        end
    end
    par1mesh, par2mesh, periods
end

# %%
function isoperiodic_test_rev1(ds, par_area, NIter, nxblock, nyblock, NTr, xpts, ypts)

    p1_st = par_area[1]          # beginning of first parameter
    p1_nd = par_area[2]          # end of first parameter
    p2_st = par_area[3]          # beginning of second parameter
    p2_nd = par_area[4]          # end of second parameter

    p1_block = range(p1_st, p1_nd, length = nxblock + 1)   # divide par_area in blocks
    p2_block = range(p2_st, p2_nd, length = nyblock + 1)   # divide par_area in blocks
    l_bipar = xpts * ypts                                  # total number of pts in each block
    total_blocks = nxblock * nyblock                       # total number of blocks
    sol_last = Matrix{Float64}(undef, NIter + 1, l_bipar)  # stores x values of solution for l_bipar (r,k) paris

    param = [zeros(2) for _ in 1:nthreads()]
    dss = [deepcopy(ds) for _ in 1:nthreads()]
    periods = Vector{Int}(undef, l_bipar * total_blocks)
    par1mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    par2mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    prog = Progress(total_blocks)
    @inbounds for ii in 1:nxblock
        par1range = range(p1_block[ii], p1_block[ii+1], length=xpts+1)[1:end-1]
        for jj in 1:nyblock
            par2range = range(p2_block[jj], p2_block[jj+1], length=ypts+1)[1:end-1]
            Threads.@threads for i = 1:xpts
                tid = threadid()
                @inbounds for j = 1:ypts
                    param[tid] .= (par1range[i], par2range[j])
                    set_parameter!(dss[tid], param[tid])                    # change the parameter values
                    tr = trajectory(dss[tid], NIter; Ttr=NTr)                # find the solution
                    sol_last[:, ypts*(i-1)+j] .= first.(tr.data) # store x-values
                end
            end
            a = l_bipar*(nyblock*(ii-1)+jj-1)
            periods[a+1:a+l_bipar] .= seqper_new.(eachcol(sol_last), tol=0.001) # seqper_new calculates periodicity
            @views meshh!(par1mesh[a+1:a+l_bipar], par2mesh[a+1:a+l_bipar], par1range, par2range) # meshgrid of parameter values
            next!(prog)
        end
    end
    par1mesh, par2mesh, periods
end

# %%
function isoperiodic_test_rev2_old(init, par_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
    u0 = SVector{2}(init)
    p1_st = par_area[1]          # beginning of first parameter
    p1_nd = par_area[2]          # end of first parameter
    p2_st = par_area[3]          # beginning of second parameter
    p2_nd = par_area[4]          # end of second parameter

    p1_block = range(p1_st, p1_nd, length = nxblock + 1)   # divide par_area in blocks
    p2_block = range(p2_st, p2_nd, length = nyblock + 1)   # divide par_area in blocks
    l_bipar = xpts * ypts                                  # total number of pts in each block
    total_blocks = nxblock * nyblock                       # total number of blocks
    sol_last = Matrix{Float64}(undef, NIter + 1, l_bipar)  # stores x values of solution for l_bipar (r,k) paris

    periods = Vector{Int}(undef, l_bipar * total_blocks)
    par1mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    par2mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    tr = [Vector{SVector{2, Float64}}(undef, NIter + 1) for _ in 1:nthreads()]
    prog = Progress(total_blocks)
    @inbounds for ii in 1:nxblock
        par1range = range(p1_block[ii], p1_block[ii+1], length=xpts+1)[1:end-1]
        for jj in 1:nyblock
            par2range = range(p2_block[jj], p2_block[jj+1], length=ypts+1)[1:end-1]
            Threads.@threads for i = 1:xpts
                tid = threadid()
                @inbounds for j = 1:ypts
                    traj!(tr[tid], dds_rule, (par1range[i], par2range[j]), u0, NIter; Ttr=NTr) # find the solution
                    sol_last[:, ypts*(i-1)+j] .= first.(tr[tid])   # store x-values
                end
            end
            a = l_bipar*(nyblock*(ii-1)+jj-1)
            periods[a+1:a+l_bipar] .= seqper_new.(eachcol(sol_last), tol=0.001) # seqper_new calculates periodicity
            @views meshh!(par1mesh[a+1:a+l_bipar], par2mesh[a+1:a+l_bipar], par1range, par2range) # meshgrid of parameter values
            next!(prog)
        end
    end
    par1mesh, par2mesh, periods
end

# %%
function isoperiodic_test_rev2(init, par_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
    u0 = SVector{2}(init)
    p1_st = par_area[1]          # beginning of first parameter
    p1_nd = par_area[2]          # end of first parameter
    p2_st = par_area[3]          # beginning of second parameter
    p2_nd = par_area[4]          # end of second parameter

    p1_block = range(p1_st, p1_nd, length = nxblock + 1)   # divide par_area in blocks
    p2_block = range(p2_st, p2_nd, length = nyblock + 1)   # divide par_area in blocks
    l_bipar = xpts * ypts                                  # total number of pts in each block
    total_blocks = nxblock * nyblock                       # total number of blocks

    periods = Vector{Int}(undef, l_bipar * total_blocks)
    par1mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    par2mesh = Vector{Float64}(undef, l_bipar * total_blocks)
    tr1 = [Vector{Float64}(undef, NIter + 1) for _ in 1:nthreads()]
    tr2 = [Vector{Float64}(undef, NIter + 1) for _ in 1:nthreads()]
    prog = Progress(total_blocks)
    @inbounds for ii in 1:nxblock
        par1range = range(p1_block[ii], p1_block[ii+1], length=xpts+1)[1:end-1]
        for jj in 1:nyblock
            par2range = range(p2_block[jj], p2_block[jj+1], length=ypts+1)[1:end-1]
            a = l_bipar*(nyblock*(ii-1)+jj-1)
            Threads.@threads for i = 1:xpts
                tid = threadid()
                @inbounds for j = 1:ypts
                    traj2d!(tr1[tid], tr2[tid],
                        dds_rule, (par1range[i], par2range[j]), u0, NIter; Ttr=NTr) # find the solution
                    periods[a + ypts*(i-1)+j] = seqper_new(tr1[tid], tol=0.001)
                end
            end
            @views meshh!(par1mesh[a+1:a+l_bipar], par2mesh[a+1:a+l_bipar], par1range, par2range) # meshgrid of parameter values
            next!(prog)
        end
    end
    par1mesh, par2mesh, periods
end

# %%
##  Test

parameter_area = [1.0 5.0 2.0 5.0]
nxblock = 1
nyblock = 1
xpts = 100
ypts = 100
NIter = 2000
NTr = 50000
init = [0.4, 0.5]
ds = dds_constructor(init)
##
println("********** Minor correction of the original:")
result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Revised 1:")
result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Revised 2:")
result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Confirmation of equivalence:")
@show result_org .== result_rev1;
@show result_rev1 .== result_rev2;

# %% tags=[]
##  Test

parameter_area = [1.0 5.0 2.0 5.0]
nxblock = 2
nyblock = 2
xpts = 100
ypts = 100
NIter = 2000
NTr = 50000
init = [0.4, 0.5]
ds = dds_constructor(init)
##
println("********** Minor correction of the original:")
result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Revised 1:")
result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Revised 2:")
result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Confirmation of equivalence:")
@show result_org .== result_rev1;
@show result_rev1 .== result_rev2;

# %%
##  Test

parameter_area = [1.0 5.0 2.0 5.0]
nxblock = 10
nyblock = 10
xpts = 10
ypts = 10
NIter = 2000
NTr = 50000
init = [0.4, 0.5]
ds = dds_constructor(init)
##
println("********** Minor correction of the original:")
result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
#result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
#result_org = @time isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Revised 1:")
result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
#result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
#result_rev1 = @time isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Revised 2:")
result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
#result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
#result_rev2 = @time isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)
println("********** Confirmation of equivalence:")
@show result_org .== result_rev1;
@show result_rev1 .== result_rev2;

# %%
using Plots
xoxp1(x) = x/(x+1)

r, k, period = @time isoperiodic_test_rev2([0.4, 0.5], [1.0 5.0 2.0 5.0], 2000, 1, 1, 50000, 80, 60)
r = reshape(r, 60, 80)
k = reshape(k, 60, 80)
period = reshape(period, 60, 80)
heatmap(vec(r[1,:]), vec(k[:,1]), xoxp1.(period); xlabel="r", ylabel="k", title="period/(period + 1)")

# %%
@btime isoperiodic_test_rev2($([0.4, 0.5]), $([1.0 5.0 2.0 5.0]), 2000, 1, 1, 50000, 80, 60)

# %%
@code_warntype isoperiodic_test_org(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)

# %%
@code_warntype isoperiodic_test_rev1(ds, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)

# %%
@code_warntype isoperiodic_test_rev2(init, parameter_area, NIter, nxblock, nyblock, NTr, xpts, ypts)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
export invariantmeasure, InvariantDistribution

import LinearAlgebra: norm

# Invariant measures (measures over the partition do not have to sum to 1)
abstract type InvariantMeasure end
abstract type AbstractRectangularInvariantMeasure <: InvariantMeasure end
abstract type AbstractTriangulationInvariantMeasure <: InvariantMeasure end

# Invariant distributions (measures over the partition must sum to 1)
abstract type AbstractInvariantDistribution{TO <: TransferOperatorApproximation, T} end

"""
    invariantmeasure(to::TransferOperatorApproximation) → iv::InvariantDistribution
    invariantmeasure(pts, method::SingleGrid) → iv::InvariantDistribution
    invariantmeasure(pts, method::SimplexPoint) → iv::InvariantDistribution
    invariantmeasure(pts, method::SimplexExact) → iv::InvariantDistribution

Compute an invariant probability distribution `iv` from an approximation to the
transfer operator over a partitioning of the input `pts` using one of the
available `method`s.

Alternatively, compute the distribution from a precomputed transfer operator approximation.

## Examples

### `SingleGrid` estimator

Use the `SingleGrid` estimator for long time series.

```julia
using DelayEmbeddings, PerronFrobenius
n = 20000
x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
τs, js = (0, -2, -3), (1, 1, 1)
pts = genembed(x, τs, js);

# Direct approach
b = RectangularBinning(5)
invariantmeasure(pts, SingleGrid(b))

# Indirect approach
to = transferoperator(pts, SingleGrid(b))
invariantmeasure(to)
```

### Triangulation estimators

The `SimplexApprox` estimator can be useful for very short time series.
The `SimplexExact` estimator can also be useful for very short time series.
It computes transition probabilities from exact simplex intersections, so
is much slower than the `SimplexApprox` estimator.

```julia
using DelayEmbeddings, PerronFrobenius
n = 30
x = [cos(2π*(t + sin(t)/30)) for t in 1:n] .+ 0.2rand(n)
τs, js = (0, -2, -3), (1, 1, 1)
pts = genembed(x, τs, js);

# Use approximate simplex intersections. Sample each simplex using
# at least 300 points
invariantmeasure(pts, SimplexPoint(), n = 300)

# Use exact simplex intersections.
invariantmeasure(pts, SimplexExact())
```
"""
function invariantmeasure end

"""
    InvariantDistribution(to::TransferOperatorApproximation, dist, nonzero_inds)

An invariant distribution `dist` computed from the transfer operator approximation `to`.
`nonzero_inds` contains the indices of `dist` with nonzero entries.
"""
struct InvariantDistribution{TO <: TransferOperatorApproximation, T} <: AbstractInvariantDistribution{TO, T}
    to::TO
    dist::Vector{T} # Distribution over the simplices
    nonzero_inds::Vector{Int} # indices of nonzero entries
end

function Base.summary(iv::InvariantDistribution{TO, T}) where {TO, T}
    method = iv.to.g.method
    data = iv.to.g.pts
    s = size(iv.to.M, 1)
    l = length(iv.to.M)
    percent_nonzero = @sprintf("%.4f", count(iv.to.M .> 0.0)/length(iv.to.M) * 100)
    f = length(iv.nonzero_inds)/length(iv.dist)
    return "InvariantDistribution over $s partition elements ($(f*100)% with positive measure)"
end

function Base.show(io::IO, iv::InvariantDistribution{TO, T}) where {TO, T}
    s = summary(iv)
    print(io, s)
end

function invariantmeasure(to::TransferOperatorApproximation;
        N::Int = 200, tol::Float64 = 1e-8, δ::Float64 = 1e-8)
    M = to.M
    #=
    # Start with a random distribution `Ρ` (big rho). Normalise it so that it
    # sums to 1 and forms a true probability distribution over the simplices.
    =#
    Ρ = rand(Float64, 1, size(M, 1))
    Ρ .= Ρ ./ sum(Ρ, dims = 2)

    #=
    # Start estimating the invariant distribution. We could either do this by
    # finding the left-eigenvector of M, or by repeated application of M on Ρ
    # until the distribution converges. Here, we use the latter approach,
    # meaning that we iterate until Ρ doesn't change substantially between
    # iterations.
    =#
    distribution = Ρ * M

    distance = norm(distribution - Ρ) / norm(Ρ)

    check = floor(Int, 1 / δ)
    check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N))
    check_pts = check_pts[check_pts .> 0]
    num_checkpts = size(check_pts, 1)
    check_pts_counter = 1

    counter = 1
    while counter <= N && distance >= tol
        counter += 1
        Ρ = distribution

        # Apply the Markov matrix to the current state of the distribution
        distribution = Ρ * M

        if (check_pts_counter <= num_checkpts &&
           counter == check_pts[check_pts_counter])

            check_pts_counter += 1
            colsum_distribution = sum(distribution, dims = 2)[1]
            if abs(colsum_distribution - 1) > δ
                distribution = distribution ./ colsum_distribution
            end
        end

        distance = norm(distribution - Ρ) / norm(Ρ)
    end
    distribution = dropdims(distribution, dims = 1)

    # Do the last normalisation and check
    colsum_distribution = sum(distribution)

    if abs(colsum_distribution - 1) > δ
        distribution = distribution ./ colsum_distribution
    end

    # Find partition elements with strictly positive measure.
    Δ = δ/size(M, 1)
    positive_measure_inds = findall(distribution .> Δ)

    return InvariantDistribution(to, distribution, positive_measure_inds)
end


text:
module ApproximationGrids

### Imports
using FFTW

### Exports
export ChebyshevGrid, RegularGrid,
       interpolate, jacobian, nodes, weights

### Implemetation
abstract type ApproximationGrid end

# Objects' string representations
function Base.show(io::IO, g::T) where {T <: ApproximationGrid}
    print(io, "$(nameof(T))(($(g.a),$(g.b)), $(g.n))")
end

include("regular.jl")
include("chebyshev.jl")

end # module


text:
module FriendlyDoLoop

using Reexport

include("Foos.jl")

include("Bars.jl")

end # module