text:
ExpCacheDict = Dict{Pair{ProbΔNode, LogisticΔNode}, Array{Float64, 2}}
MomentCacheDict = Dict{Tuple{ProbΔNode, LogisticΔNode, Int64}, Array{Float64, 2}}

struct ExpectationCache
    f::ExpCacheDict
    fg::ExpCacheDict
end
ExpectationCache() = ExpectationCache(ExpCacheDict(), ExpCacheDict())

struct MomentCache
    f::ExpCacheDict
    fg::MomentCacheDict
end
MomentCache() = MomentCache( ExpCacheDict(),  MomentCacheDict())

# Find a better way to cache n_choose_k values
max_k = 31
choose_cache = [ 1.0 * binomial(i,j) for i=0:max_k+1, j=0:max_k+1 ]
@inline function choose(n::Int, m::Int)
    return choose_cache[n+1, m+1]
end

# On Tractable Computation of Expected Predictions (https://arxiv.org/abs/1910.02182)
"""
Missing values should be denoted by -1
"""
function Expectation(pc::ProbΔ, lc::LogisticΔ, data::XData{Int8})
    # 1. Get probability of each observation
    fc, log_likelihoods = marginal_log_likelihood_per_instance(pc, data)
    p_observed = exp.( log_likelihoods )

    # 2. Expectation w.r.t. P(x_m, x_o)
    cache = ExpectationCache()
    results_unnormalized = exp_g(pc[end], lc[end-1], data, cache) # skipping the bias node of lc

    # 3. Expectation w.r.t P(x_m | x_o)
    results = transpose(results_unnormalized) ./ p_observed

    # 4. Add Bias terms
    biases = lc[end].thetas
    results .+= biases

    results, cache
end

function Moment(pc::ProbΔ, lc::LogisticΔ, data::XData{Int8}, moment::Int)
    # 1. Get probability of each observation
    fc, log_likelihoods = marginal_log_likelihood_per_instance(pc, data)
    p_observed = exp.( log_likelihoods )

    # 2. Moment w.r.t. P(x_m, x_o)
    cache = MomentCache()
    biases = lc[end].thetas
    results_unnormalized = zeros(num_examples(data), classes(lc[end]))

    for z = 0:moment-1
        results_unnormalized .+= choose(moment, z) .* (biases .^ (z)) .* transpose(moment_g(pc[end], lc[end-1], data, moment - z, cache))
    end

    # 3. Moment w.r.t P(x_m | x_o)
    results = results_unnormalized ./ p_observed

    # 4. Add Bias^moment terms
    results .+= biases .^ (moment)

    results, cache
end

function ExpectationUpward(pc::ProbΔ, lc::LogisticΔ, data::XData{Int8})
        # 1. Get probability of each observation
        fc, log_likelihoods = marginal_log_likelihood_per_instance(pc, data)
        p_observed = exp.( log_likelihoods )

        # 2. Expectation w.r.t. P(x_m, x_o)
        exps_flow = exp_pass_up(pc, lc, data)
        results_unnormalized = exps_flow[end].fg

        # 3. Expectation w.r.t P(x_m | x_o)
        results = transpose(results_unnormalized) ./ p_observed

        # 4. Add Bias terms
        biases = lc[end].thetas
        results .+= biases

        results, exps_flow
end

# exp_f (pr-constraint) is originally from:
#   Arthur Choi, Guy Van den Broeck, and Adnan Darwiche. Tractable learning for structured probability spaces: A case study in learning preference distributions. In Proceedings of IJCAI, 2015.

function exp_f(n::Prob⋁, m::Logistic⋁, data::XData{Int8}, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        value = zeros(1 , num_examples(data) )
        pthetas = [exp(n.log_thetas[i]) for i in 1:length(n.children)]
        @fastmath @simd for i in 1:length(n.children)
            @simd for j in 1:length(m.children)
                value .+= (pthetas[i] .* exp_f(n.children[i], m.children[j], data, cache))
            end
        end
        return value
    end
end

function exp_f(n::Prob⋀, m::Logistic⋀, data::XData{Int8}, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        value = ones(1 , num_examples(data) )
        @fastmath for (i,j) in zip(n.children, m.children)
            value .*= exp_f(i, j, data, cache)
        end
        return value
        # exp_f(n.children[1], m.children[1], data, cache) .* exp_f(n.children[2], m.children[2], data, cache)
    end
end

@inline function exp_f(n::ProbLiteral, m::LogisticLiteral, data::XData{Int8}, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        value = zeros(1 , num_examples(data) )
        var = lit2var(literal(m))
        X = feature_matrix(data)
        if positive(n) && positive(m)
            # value[1, X[:, var] .== -1 ] .= 1.0  # missing observation always agrees
            # value[1, X[:, var] .== 1 ] .= 1.0 # positive observations
            value[1, X[:, var] .!= 0 ] .= 1.0 # positive or missing observations
        elseif negative(n) && negative(m)
            # value[1, X[:, var] .== -1 ] .= 1.0  # missing observation always agrees
            # value[1, X[:, var] .== 0 ] .= 1.0 # negative observations
            value[1, X[:, var] .!= 1 ] .= 1.0 # negative or missing observations
        end
        return value
    end
end

"""
Has to be a Logistic⋁ with only one child, which is a leaf node
"""
@inline function exp_f(n::ProbLiteral, m::Logistic⋁, data::XData{Int8}, cache::Union{ExpectationCache, MomentCache})
    @inbounds get!(cache.f, Pair(n, m)) do
        exp_f(n, m.children[1], data, cache)
    end
end

#######################################################################
######## exp_g, exp_fg
########################################################################

@inline function exp_g(n::Prob⋁, m::Logistic⋁, data::XData{Int8}, cache::ExpectationCache)
    exp_fg(n, m, data, cache) # exp_fg and exp_g are the same for OR nodes
end

# function exp_g(n::Prob⋀, m::Logistic⋀, data::XData{Int8}, cache::ExpectationCache)
#     value = zeros(classes(m) , num_examples(data))
#     @fastmath for (i,j) in zip(n.children, m.children)
#         value .+= exp_fg(i, j, data, cache)
#     end
#     return value
#     # exp_fg(n.children[1], m.children[1], data, cache) .+ exp_fg(n.children[2], m.children[2], data, cache)
# end

function exp_fg(n::Prob⋁, m::Logistic⋁, data::XData{Int8}, cache::ExpectationCache)
    @inbounds get!(cache.fg, Pair(n, m)) do
        value = zeros(classes(m) , num_examples(data) )
        pthetas = [exp(n.log_thetas[i]) for i in 1:length(n.children)]
        @fastmath @simd for i in 1:length(n.children)
            for j in 1:length(m.children)
                value .+= (pthetas[i] .* m.thetas[j,:]) .* exp_f(n.children[i], m.children[j], data, cache)
                value .+= pthetas[i] .* exp_fg(n.children[i], m.children[j], data, cache)
            end
        end
        return value
    end
end

function exp_fg(n::Prob⋀, m::Logistic⋀, data::XData{Int8}, cache::ExpectationCache)
    @inbounds get!(cache.fg, Pair(n, m)) do
        # Assuming 2 children
        value = exp_f(n.children[1], m.children[1], data, cache) .* exp_fg(n.children[2], m.children[2], data, cache)
        value .+= exp_f(n.children[2], m.children[2], data, cache) .* exp_fg(n.children[1], m.children[1], data, cache)
        return value
    end
end

"""
Has to be a Logistic⋁ with only one child, which is a leaf node
"""
@inline function exp_fg(n::ProbLiteral, m::Logistic⋁, data::XData{Int8}, cache::ExpectationCache)
    @inbounds get!(cache.fg, Pair(n, m)) do
        m.thetas[1,:] .* exp_f(n, m, data, cache)
    end
end

@inline function exp_fg(n::ProbLiteral, m::LogisticLiteral, data::XData{Int8}, cache::ExpectationCache)
    #dont know how many classes, boradcasting does the job
    zeros(1 , num_examples(data))
end

#######################################################################
######## moment_g, moment_fg
########################################################################

@inline function moment_g(n::Prob⋁, m::Logistic⋁, data::XData{Int8}, moment::Int, cache::MomentCache)
    get!(cache.fg, (n, m, moment)) do
        moment_fg(n, m, data, moment, cache)
    end
end

"""
Calculating  E[g^k * f]
"""
function moment_fg(n::Prob⋁, m::Logistic⋁, data::XData{Int8}, moment::Int, cache::MomentCache)
    if moment == 0
        return exp_f(n, m, data, cache)
    end

    get!(cache.fg, (n, m, moment)) do
        value = zeros(classes(m) , num_examples(data) )
        pthetas = [exp(n.log_thetas[i]) for i in 1:length(n.children)]
        @fastmath @simd for i in 1:length(n.children)
            for j in 1:length(m.children)
                for z in 0:moment
                    value .+= pthetas[i] .* choose(moment, z) .* m.thetas[j,:].^(moment - z) .* moment_fg(n.children[i], m.children[j], data, z, cache)
                end
            end
        end
        return value
    end
end

@inline function moment_fg(n::ProbLiteral, m::Logistic⋁, data::XData{Int8}, moment::Int, cache::MomentCache)
    get!(cache.fg, (n, m, moment)) do
        m.thetas[1,:].^(moment) .* exp_f(n, m, data, cache)
    end
end

@inline function moment_fg(n::ProbLiteral, m::LogisticLiteral, data::XData{Int8}, moment::Int, cache::MomentCache)
    #dont know how many classes, boradcasting does the job
    if moment == 0
        exp_f(n, m, data, cache)
    else
        zeros(1, num_examples(data))
    end
end

function moment_fg(n::Prob⋀, m::Logistic⋀, data::XData{Int8}, moment::Int, cache::MomentCache)
    if moment == 0
        return exp_f(n, m, data, cache)
    end
    get!(cache.fg, (n, m, moment)) do
        value = moment_fg(n.children[1], m.children[1], data, 0, cache) .* moment_fg(n.children[2], m.children[2], data, moment, cache)

        for z in 1:moment
            value .+= choose(moment, z) .* moment_fg(n.children[1], m.children[1], data, z, cache) .* moment_fg(n.children[2], m.children[2], data, moment - z, cache)
        end
        return value
    end
end


text:
ENV["JuLearnBackend"]="native"
using JuLearn#,Plots
@Eval lr=0.01 δ=0. momentum=0. batchsize=50 ▽max=0.01 λ=0.0
include("dataProvider.jl")
include("predict.jl")

inputBase,outputBase=dataProvider()
net=chain([:d,:h,:h,:l],[1 1;1 1000;1000 1;1 1],:regression);
@time err=backprop(net,inputBase,outputBase,nIter=1000)

# plot(err)
# plot(inputBase',predict(net,inputBase)')
# plot!(inputBase',outputBase')

ext: jl


text:
export PathResult, PathTrackerStatus, PathTracker,
       pathtracker, pathtracker_startsolutions, solution,
       accuracy, residual, start_solution, isfailed, isatinfinity,
       issingular, isnonsingular, isprojective, isaffine, set_parameters!

const pathtracker_supported_keywords = [
    :at_infinity_check, :min_step_size_endgame_start,
    :min_val_accuracy, :samples_per_loop,
    :max_winding_number, :max_affine_norm,
    :overdetermined_min_accuracy, :overdetermined_min_residual,
    :cond_eg_start, :min_cond_at_infinity,
    :t_eg_start, :tol_val_inf_accurate, :tol_val_finite_accurate, :accuracy, :accuracy_eg]

###############
## VALUATION ##
###############

mutable struct Valuation
    v::Vector{Float64}
    v̇::Vector{Float64}
    v̈::Vector{Float64}
    v_data::NTuple{3,Vector{Float64}}
    s_data::NTuple{3,Float64}
end

@enum ValuationVerdict begin
    VALUATION_INDECISIVE
    VALUATION_FINITE
    VALUATION_AT_INFINIY
end

function Valuation(x::ProjectiveVectors.PVector; at_infinity_check::Bool=true)
    if at_infinity_check
        Valuation(length(x) - length(ProjectiveVectors.dims(x)))
    else
        Valuation(length(x))
    end
end
Valuation(x::Vector; kwargs...) = Valuation(length(x))
function Valuation(n::Integer)
    v = zeros(n)
    v̇ = zeros(n)
    v̈ = zeros(n)
    v_data = (zeros(n), zeros(n), zeros(n))
    s_data = (NaN, NaN, NaN)
    Valuation(v, v̇, v̈, v_data, s_data)
end

function Base.show(io::IO, val::Valuation)
    println(io, typeof(val), ":")
    for name in [:v, :v̇, :v̈]
        println(io, " • ", name, " → ", getfield(val, name))
    end
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", v::Valuation) = v

function reset!(val::Valuation)
    val.v .= NaN
    val.v̇ .= NaN
    val.v̈ .= NaN
    val.s_data = (NaN, NaN, NaN)
    val
end

function update!(val::Valuation, z::Vector, ż, s, Δs::Float64, ::Val)
    _update!(val, z, ż, s, Δs)
end
function update!(val::Valuation, z::PVector, ż, s, Δs::Float64, ::Val{false})
    _update!(val, z, ż, s, Δs)
end
function update!(val::Valuation, z::PVector, ż, s, Δs::Float64, ::Val{true})
    _update_affine!(val, z, ż, s, Δs)
end
function _update!(val::Valuation, z::AbstractVector, ż, s, Δs::Float64)
    @unpack v, v_data, s_data = val
    (v₂, v₁, v₀), (s₂, s₁, _) = v_data, s_data
    v .= v₀ .= ν.(z, ż)
    # cycle data around
    val.v_data = (v₀, v₂, v₁)
    val.s_data = (s, s₂, s₁)

    !isnan(s₁) && finite_differences!(val)

    val
end
function _update_affine!(val::Valuation, z::PVector, ż, s, Δs::Float64)
    @unpack v, v_data, s_data = val
    (v₂, v₁, v₀), (s₂, s₁, _) = v_data, s_data
    k = 1
    for (rⱼ, j) in ProjectiveVectors.dimension_indices_homvars(z)
        vⱼ = ν(z[j], ż[j])
        for i in rⱼ
            v[k] = v₀[k] =  ν(z[i], ż[i]) - vⱼ
            k += 1
        end
    end
    # cycle data around
    val.v_data = (v₀, v₂, v₁)
    val.s_data = (s, s₂, s₁)

    !isnan(s₁) && finite_differences!(val)

    val
end

function ν(z::Complex, ż::Complex)
    x, y = reim(z); ẋ, ẏ = reim(ż);
    -(x * ẋ + y * ẏ) / abs2(z)
end

"""
    finite_differences!(val::Valuation)

Use a finite difference scheme to approximate the first and second derivative
of the valuation map.

Since we have a non-uniform grid, we need a more elaborate difference scheme.
The implementation follows the formulas derived in [^BS05]

[^BS05]: Bowen, M. K., and Ronald Smith. "Derivative formulae and errors for non-uniformly spaced points." Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 461.2059 (2005): 1975-1997.
"""
function finite_differences!(val::Valuation)
    @unpack v_data, s_data, v̇, v̈ = val
    v₃, v₂, v₁ = v_data
    s₃, s₂, s₁ = s_data
    s₃₁ = s₃ - s₁
    a₃, a₁ = s₃ - s₂, s₁ - s₂
    for i in eachindex(v̇)
        v̇[i] = (a₃*v₁[i])/(a₁*s₃₁) - ((a₁+a₃)*v₂[i])/(a₁*a₃) - (a₁*v₃[i])/(s₃₁*a₃)
        v̈[i] = -2v₁[i]/((a₁)*s₃₁) + 2v₂[i]/(a₁*a₃) + 2v₃[i]/(s₃₁*a₃)
    end
    nothing
end

@inline function update!(val::Valuation, core_tracker::CoreTracker; at_infinity_check::Bool=true)
    z = core_tracker.state.x
    ż = core_tracker.state.ẋ
    Δs = core_tracker.state.Δs_prev
    s = core_tracker.state.s
    if at_infinity_check
        update!(val, z, ż, s, Δs, Val(true))
    else
        update!(val, z, ż, s, Δs, Val(false))
    end
end

#################
## PATHTRACKER ##
#################

module PathTrackerStatus

    import ..CoreTrackerStatus

    @doc """
        PathTrackerStatus.states

    The possible return codes the path tracker can return are

    * `PathTrackerStatus.success`
    * `PathTrackerStatus.at_infinity`
    * `PathTrackerStatus.terminated_maximal_iterations`
    * `PathTrackerStatus.terminated_invalid_startvalue`
    * `PathTrackerStatus.terminated_step_size_too_small`
    * `PathTrackerStatus.terminated_singularity`
    * `PathTrackerStatus.terminated_ill_conditioned`
    * `PathTrackerStatus.terminated`
    * `PathTrackerStatus.post_check_failed`
    * `PathTrackerStatus.excess_solution`
    """
    @enum states begin
        tracking
        success
        at_infinity
        terminated_invalid_startvalue
        terminated_maximal_iterations
        terminated_step_size_too_small
        terminated_singularity
        terminated_ill_conditioned
        tracker_failed
        post_check_failed
        excess_solution
    end

    function status(code::CoreTrackerStatus.states)
        if code == CoreTrackerStatus.success
            return success
        elseif code == CoreTrackerStatus.terminated_invalid_startvalue
            return terminated_invalid_startvalue
        elseif code == CoreTrackerStatus.terminated_maximal_iterations
            return terminated_maximal_iterations
        elseif code == CoreTrackerStatus.terminated_step_size_too_small
            return terminated_step_size_too_small
        elseif code == CoreTrackerStatus.terminated_singularity
            return terminated_singularity
        elseif code == CoreTrackerStatus.terminated_ill_conditioned
            return terminated_ill_conditioned
        else
            # this shouldn't happen
            return tracker_failed
        end
    end
end

mutable struct PathTrackerOptions
    at_infinity_check::Bool
    samples_per_loop::Int
    max_winding_number::Int
    # The minimal residual a solution needs to have to be considered
    # a solution of the original system (only applied for singular solutions)
    overdetermined_min_residual::Float64
    overdetermined_min_accuracy::Float64
    # minimial condition number where the endgame starts
    cond_eg_start::Float64
    min_cond_at_infinity::Float64
    # maximal t where the endgame starts
    t_eg_start::Float64
    tol_val_inf_accurate::Float64
    tol_val_finite_accurate::Float64
    accuracy::Float64
    accuracy_eg::Float64
end

function PathTrackerOptions(;
            at_infinity_check=true,
            samples_per_loop::Int=12,
            max_winding_number::Int=12,
            overdetermined_min_residual::Float64=1e-3,
            overdetermined_min_accuracy::Float64=1e-4,
            cond_eg_start::Float64=1e4,
            min_cond_at_infinity::Float64=1e7,
            t_eg_start::Float64=0.1,
            tol_val_inf_accurate::Float64=1e-4,
            tol_val_finite_accurate::Float64=1e-3,
            accuracy::Float64=error("You have to set `accuracy`"),
            accuracy_eg::Float64=min(accuracy, 1e-5))

    PathTrackerOptions(at_infinity_check, samples_per_loop, max_winding_number,
                       overdetermined_min_residual,
                       overdetermined_min_accuracy, cond_eg_start,
                       min_cond_at_infinity,
                       t_eg_start,
                       tol_val_inf_accurate, tol_val_finite_accurate, accuracy, accuracy_eg)
end

Base.show(io::IO, opts::PathTrackerOptions) = print_fieldnames(io, opts)
Base.show(io::IO, ::MIME"application/prs.juno.inline", opts::PathTrackerOptions) = opts

mutable struct PathTrackerState{V<:AbstractVector}
    status::PathTrackerStatus.states
    s::Float64
    prediction::V
    solution::V
    solution_accuracy::Float64
    solution_cond::Float64
    winding_number::Int
    val::Valuation
end

function PathTrackerState(x; at_infinity_check::Bool=true)
    status = PathTrackerStatus.tracking
    s = 0.0
    prediction = copy(x)
    solution = copy(x)
    solution_accuracy = solution_cond = NaN
    winding_number = 0
    val = Valuation(x; at_infinity_check=at_infinity_check)

    PathTrackerState(status, s, prediction, solution,
                    solution_accuracy, solution_cond,
                    winding_number, val)
end

function reset!(state::PathTrackerState)
    state.status = PathTrackerStatus.tracking
    state.s = 0.0
    state.prediction .= zero(eltype(state.prediction))
    state.solution .= zero(eltype(state.prediction))
    state.solution_accuracy = state.solution_cond = NaN
    state.winding_number = 0
    reset!(state.val)

    state
end

struct PathTrackerCache{T, V<:AbstractVector{Complex{T}}, S<:AbstractSystem, NC<:AbstractNewtonCache}
    unit_roots::Vector{ComplexF64}
    base_point::V
    target_system::S
    target_residual::Vector{Complex{T}}
    target_jacobian::Jacobian{Complex{T}}
    target_newton_cache::NC
    weighted_ip::WeightedIP{T}
end

function PathTrackerCache(prob::Problem, core_tracker::CoreTracker)
    unit_roots = ComplexF64[]
    base_point = copy(core_tracker.state.x)

    x = currx(core_tracker)
    if is_squared_up_system(core_tracker.homotopy)
        F = squared_up_system(core_tracker.homotopy).F
    else
        F = FixedHomotopy(core_tracker.homotopy, Inf)
    end
    if affine_tracking(core_tracker)
        target_system = F
    elseif pull_back_is_to_affine(prob)
        target_system = PatchedSystem(F, state(EmbeddingPatch(), x))
    else # result is projective
        target_system = PatchedSystem(F, state(OrthogonalPatch(), x))
    end
    target_newton_cache = newton_cache(target_system, x)

    res = evaluate(target_system, currx(core_tracker), target_newton_cache.system_cache)
    jac = similar(res, size(target_system))
    weighted_ip = WeightedIP(x)
    PathTrackerCache(unit_roots, base_point, target_system, res, Jacobian(Random.rand!(jac)), target_newton_cache, weighted_ip)
end

is_squared_up_system(::StraightLineHomotopy{<:Any,<:SquaredUpSystem}) = true
is_squared_up_system(H::LogHomotopy) = is_squared_up_system(H.homotopy)
is_squared_up_system(::AbstractHomotopy) = false

squared_up_system(H::StraightLineHomotopy{<:Any,<:SquaredUpSystem}) = H.target
squared_up_system(H::LogHomotopy) = squared_up_system(H.homotopy)

"""
     PathTracker{Prob<:AbstractProblem, T, V<:AbstractVector{T}, CT<:CoreTracker}

`PathTracker` the way to track single paths. It combines the core path tracking routine
with an endgame, i.e., it can also deal with singular solutions as well as diverging paths.
We call a diverged path a path going to infinity.
By convention a path is always tracked from t₁ > 0 towards 0.
During the path tracking an approximation of the valuation of a Puiseux series expansion of the solution is computed.
This is used to decide whether a path is diverging.
To compute singular solutions Cauchy's integral formula is used. There you have to trace out loops around the solution.
The number of loops necessary to arrive back at the start point is called the *winding number*.
In order to construct a `PathTracker` it is recommended to use the [`pathtracker`](@ref) and
[`pathtracker_startsolutions`](@ref) helper functions.
With a `PathTracker` constructed you can track a single path using the [`track`](@ref) method.
The result of this will be a [`PathResult`](@ref).

## Keyword arguments
`PathTracker` is a wrapper around [`CoreTracker`](@ref)
and thus it is possible to set all options which are available for [`CoreTracker`](@ref).
There are the following `PathTracker` specific options:

### General endgame parameters
* `accuracy_eg::Float64=min(accuracy, 1e-5))`: It is possible to change the accuracy during the path tracking. Usually you want lower the accuracy.
* `cond_eg_start::Float64=1e4`: The endgame is only started if the condition of the Jacobian is larger than this threshold.
* `max_winding_number::Int=12`: This limits the maximal number of loops taken in applying Cauchy's formula.
* `min_cond_at_infinity::Float64=1e7`: A path is declared as going to infinity only if it's Jacobian is also larger than this threshold.
* `samples_per_loop::Int=12`: To compute singular solutions Cauchy's integral formula is used. The accuracy of the solutions increases with the number of samples per loop.
* `t_eg_start::Float64=0.1`: The endgame starts only if `t` is smaller than this threshold.
* `tol_val_inf_accurate::Float64=1e-4`: A valuation which would result in a path declared as going to infinity is only accepted if the estimated accuracy of the valuation is less than this threshold.
* `tol_val_finite_accurate::Float64=1e-3`: A valuation which would result in a proper solution is only accepted if the estimated accuracy of the valuation is less than this threshold. This is only affects solutions where the path has at some point near 0 a condition number larger than `cond_eg_start`.

### Overdetermined system specific
* `overdetermined_min_accuracy=1e-5`: The minimal accuracy a non-singular solution needs to have to be considered a solution of the original system.
* `overdetermined_min_residual=1e-3`: The minimal residual a singular solution needs to have to be considered a solution of the original system.
"""
struct PathTracker{V<:AbstractVector, Prob<:AbstractProblem, PTC<:PathTrackerCache, CT<:CoreTracker}
    problem::Prob
    core_tracker::CT
    state::PathTrackerState{V}
    options::PathTrackerOptions
    cache::PTC
end

function PathTracker(prob::AbstractProblem, x::AbstractVector{<:Number};
                at_infinity_check=default_at_infinity_check(prob),
                accuracy=default_accuracy(prob),
                min_step_size=1e-30, kwargs...)

    core_tracker_supported, optionskwargs = splitkwargs(kwargs, coretracker_supported_keywords)
    core_tracker = CoreTracker(prob, x;
                        log_transform=true, predictor=Pade21(),
                        min_step_size=min_step_size,
                        accuracy=accuracy,
                        core_tracker_supported...)
    state = PathTrackerState(core_tracker.state.x; at_infinity_check=at_infinity_check)
    options = PathTrackerOptions(;at_infinity_check=at_infinity_check,
                                  accuracy=accuracy, optionskwargs...)
    cache = PathTrackerCache(prob, core_tracker)
    PathTracker(prob, core_tracker, state, options, cache)
end

Base.show(io::IO, tracker::PathTracker) = print(io, "PathTracker")
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::PathTracker) = x

seed(PT::PathTracker) = PT.problem.seed
default_at_infinity_check(prob::Problem{AffineTracking}) = true
default_at_infinity_check(prob::Problem{ProjectiveTracking}) = homvars(prob) !== nothing
default_accuracy(prob::Problem) = pull_back_is_to_affine(prob) ? 1e-6 : 1e-8

"""
    track!(tracker::PathTracker, x₁, t₁::Float64=1.0; options...)::PathTrackerStatus.states

Track the path with start solution `x₁` from `t₁` towards `t=0`.

Possible values for the options are
* `accuracy::Float64`
* `max_corrector_iters::Int`
* `max_steps::Int`
* `start_parameters::AbstractVector`
* `target_parameters::AbstractVector`
"""
function track!(tracker::PathTracker, x₁, s₁=0.0, s₀=-log(tracker.core_tracker.options.min_step_size);
        start_parameters::Union{Nothing,<:AbstractVector}=nothing,
        target_parameters::Union{Nothing,<:AbstractVector}=nothing,
        kwargs...)

    @unpack core_tracker, state, options, cache = tracker
    prev_options = set_options!(core_tracker; kwargs...)
    set_parameters!(tracker; start_parameters=start_parameters, target_parameters=target_parameters)

    _track!(tracker, x₁, s₁, s₀)

    set_options!(core_tracker; prev_options...)
    state.status
end

function _track!(tracker::PathTracker, x₁, s₁::Real, s₀::Real)
    @unpack core_tracker, state, options, cache = tracker
    # For performance reasons we single thread blas
    n_blas_threads = single_thread_blas()

    embed!(core_tracker.state.x, tracker.problem, x₁)
    setup!(core_tracker, core_tracker.state.x, s₁, s₀)
    set_accuracy!(core_tracker, options.accuracy)
    reset!(state)
    # Handle the case that the start value is already a solution
    if residual(tracker, core_tracker.state.x) < 1e-14
        state.status = PathTrackerStatus.success
        state.prediction .= core_tracker.state.x
        state.s = Inf
        # Set winding_number to 1 to get into the possible singular solution case in
        # `check_and_refine_solution!`
        state.winding_number = 1
    end

    s_eg_start = -log(options.t_eg_start)
    while state.status == PathTrackerStatus.tracking
        step!(core_tracker)
        state.s = real(currt(core_tracker))
        check_terminated!(core_tracker)

        if core_tracker.state.status == CoreTrackerStatus.success
            state.status = PathTrackerStatus.success
            break
        elseif core_tracker.state.status != CoreTrackerStatus.tracking
            state.status = PathTrackerStatus.status(core_tracker.state.status)
            break
        end

        # We only care if we moved forward
        core_tracker.state.last_step_failed && continue
        # update valuation and associated data
        update!(state.val, core_tracker; at_infinity_check=tracker.options.at_infinity_check)
        verdict = judge(state.val, options)
        # If we are too early, we don't work on the endgame
        state.s > s_eg_start || verdict == VALUATION_FINITE || continue
        # If the condition number is too low we also do not care
        start_eg(core_tracker, options) || continue
        # possibly reduce desired accuracy
        if tracker.options.at_infinity_check &&
           verdict == VALUATION_AT_INFINIY &&
           start_infinity_cutoff(core_tracker, options)

            state.status = PathTrackerStatus.at_infinity
            break
        end

        if verdict == VALUATION_FINITE
            set_accuracy!(core_tracker, options.accuracy_eg)
            retcode = predict_with_cauchy_integral_method!(state, core_tracker, options, cache)
            if retcode == :success
                state.status = PathTrackerStatus.success
                break
            elseif retcode == :max_winding_number
                continue
            else # path tracker failed during the loops -> break
                state.status = PathTrackerStatus.tracker_failed
                break
            end
        end
    end

    check_and_refine_solution!(tracker)
    # We have to set the number of blas threads to the previous value
    n_blas_threads > 1 && set_num_BLAS_threads(n_blas_threads)
    state.status
end

function start_eg(core_tracker, options)
    C = max(cond(core_tracker), exp10(digits_lost(core_tracker)))
    C > options.cond_eg_start
end

function start_infinity_cutoff(core_tracker, options)
    C = max(cond(core_tracker), exp10(digits_lost(core_tracker)))
    C > options.min_cond_at_infinity
end

"""
    judge(val::Valuation, J::Jacobian)::ValuationVerdict

Judge the current valuation.
"""
function judge(val::Valuation, options::PathTrackerOptions)
    @unpack v, v̇, v̈ = val
    finite = true
    for i in eachindex(v)
        # A coordinate goes to infinity if its valuation is negative.
        # We argue that a valuation is "stable" if it's first and second
        # derivative are "small".
        if v[i] < -0.1 &&
            abs(v̇[i]) < options.tol_val_inf_accurate &&
            abs(v̈[i]) < options.tol_val_inf_accurate

            return VALUATION_AT_INFINIY
        end

        if v[i] < -options.tol_val_finite_accurate ||
            !(abs(v̇[i]) < options.tol_val_finite_accurate &&
              abs(v̈[i]) < options.tol_val_finite_accurate)

            finite = false
        end
    end
    finite && return VALUATION_FINITE

    VALUATION_INDECISIVE
end

function at_infinity_post_check(val::Valuation, options::PathTrackerOptions)
    for (vᵢ, v̇ᵢ) in zip(val.v, val.v̇)
        if vᵢ < -0.1 && abs(v̇ᵢ) < sqrt(options.tol_val_inf_accurate)
            return true
        end
    end
    false
end

"""
    track(tracker::PathTracker, x₁, t₁::Float64=1.0; path_number::Int=1, details::Symbol=:default, options...)::PathResult

Track the path with start solution `x₁` from `t₁` towards `t=0`. The `details` options controls
the level of details of the informations available in the [`PathResult`](@ref).

Possible values for the options are
* `accuracy::Float64`
* `max_corrector_iters::Int`
* `max_steps::Int`
* `start_parameters::AbstractVector`
* `target_parameters::AbstractVector`
"""
function track(tracker::PathTracker, x₁, t₁=0.0, t₀=-log(tracker.core_tracker.options.min_step_size); path_number::Int=1, details::Symbol=:default, kwargs...)
    track!(tracker, x₁, t₁, t₀; kwargs...)
    PathResult(tracker, x₁, path_number; details=details)
end

function set_options!(core_tracker::CoreTracker;
                      accuracy::Float64=core_tracker.options.accuracy,
                      max_corrector_iters::Int=core_tracker.options.max_corrector_iters,
                      max_steps::Int=core_tracker.options.max_steps)
    core_tracker.options.accuracy = accuracy
    core_tracker.options.max_corrector_iters = max_corrector_iters
    core_tracker.options.max_steps = max_steps
    (accuracy=accuracy, max_corrector_iters=max_corrector_iters, max_steps=max_steps)
end

"""
    set_parameters!(tracker::PathTracker; start_parameters=nothing, target_parameters=nothing)

Set the parameters of a parameter homotopy.
"""
@inline function set_parameters!(tracker::PathTracker; start_parameters=nothing, target_parameters=nothing)
    if start_parameters !== nothing
        set_start_parameters!(basehomotopy(tracker.core_tracker.homotopy), start_parameters)
    end
    if target_parameters !== nothing
        set_target_parameters!(basehomotopy(tracker.core_tracker.homotopy), target_parameters)
    end
    nothing
end

vector_at_infinity(x::AbstractVector, tol) = false
vector_at_infinity(x::PVector, tol) = ProjectiveVectors.norm_affine_chart(x) > tol

##########
# Cauchy #
##########

function set_eg_accuracy!(core_tracker::CoreTracker, options::PathTrackerOptions)
    set_accuracy!(core_tracker, options.accuracy_eg)
end

"""
    predict_with_cauchy_integral_method!(state, core_tracker, options, cache)

Try to predict the value of `x(0)` using [Cauchy's integral formula](https://en.wikipedia.org/wiki/Cauchy%27s_integral_formula)
At each iteration we are at some point ``(x, t)``. We then track the polygon defined
by ``te^{i2πk/n}`` until we end again at ``x``. Here ``n`` is the number of samples we take
per loop.
Returns a symbol indicating whether the prediction was successfull.
"""
function predict_with_cauchy_integral_method!(state, core_tracker, options, cache)
    @unpack prediction = state
    @unpack unit_roots, base_point = cache
    @unpack samples_per_loop, max_winding_number = options

    initial_segment = core_tracker.state.segment
    initial_s = core_tracker.state.s
    initial_Δs = core_tracker.state.Δs
    initial_Δs_prev = core_tracker.state.Δs_prev
    initial_step_size = core_tracker.options.initial_step_size
    s = real(currt(core_tracker))

    base_point .= currx(core_tracker)
    prediction .= zero(eltype(state.prediction))

    # during the loop we fix the affine patch
    fix_patch!(core_tracker)

    m = k = 1
    ∂θ = 2π / samples_per_loop
    core_tracker.options.initial_step_size = ∂θ #0.5∂θ
    while m ≤ max_winding_number
        θⱼ = 0.0
        for j=1:samples_per_loop
            θⱼ₋₁ = θⱼ
            θⱼ += ∂θ

            retcode = track!(core_tracker, currx(core_tracker), s + im*θⱼ₋₁, s + im*θⱼ; loop=true)
            if retcode != CoreTrackerStatus.success
                # during the loop we fixed the affine patch
                unfix_patch!(core_tracker)
                core_tracker.state.segment = initial_segment
                core_tracker.options.initial_step_size = initial_step_size
                core_tracker.state.s = initial_s
                core_tracker.state.Δs = initial_Δs
                core_tracker.state.Δs_prev = initial_Δs_prev
                core_tracker.state.status = CoreTrackerStatus.tracking
                return Symbol(retcode)
            end

            prediction .+= currx(core_tracker)
        end

        if distance(base_point, currx(core_tracker), inner(core_tracker)) < 4core_tracker.options.accuracy
            break
        end

        m += 1
    end
    core_tracker.state.segment = initial_segment
    core_tracker.options.initial_step_size = initial_step_size
    core_tracker.state.s = initial_s
    core_tracker.state.Δs = initial_Δs
    core_tracker.state.Δs_prev = initial_Δs_prev
    core_tracker.state.status = CoreTrackerStatus.tracking
    # we have to undo the fixing of the patch
    unfix_patch!(core_tracker)

    if m > max_winding_number
        return :max_winding_number
    end
    state.winding_number = m
    prediction ./= m * samples_per_loop
    :success
end

fix_patch!(tracker::CoreTracker) = tracker.options.update_patch = false
unfix_patch!(tracker::CoreTracker) = tracker.options.update_patch = true

"""
    type_of_x(pathtracker)

Returns the type of `x`.
"""
type_of_x(tracker::PathTracker) = typeof(tracker.core_tracker.state.x)

"""
    check_and_refine_solution!(pathtracker)

In the case of success, we store the solution in `state.solution`
and try to refine the solution to the desired accuracy.
This also makes sure that our accuracy is correct after a possible pull back.
"""
function check_and_refine_solution!(tracker::PathTracker)
    @unpack core_tracker, state, options, cache = tracker

    # The tracking failed. Let's look at the last valuation and see whether we can
    # resonably classify this as going to infinity
    if (state.status == PathTrackerStatus.terminated_ill_conditioned ||
        state.status == PathTrackerStatus.terminated_step_size_too_small ||
        state.status == PathTrackerStatus.terminated_singularity) &&
        at_infinity_post_check(state.val, options)

       state.status = PathTrackerStatus.at_infinity
   end

    if state.status ≠ PathTrackerStatus.success
        state.solution .= currx(core_tracker)
        return nothing
    end

    # We have to differentiate

    # 1.a) Non-singular solution of regular system
    # 1.b) Singular solution of regular system (Cauchy engame used)
    # 2.a) Non-singular solution of squared up system
    # 2.b) Singular solution of squared up system (Cauchy engame used)

    is_singular = state.winding_number > 0
    if is_singular
        state.solution .= state.prediction
        # check that solution is indeed singular
        if !pull_back_is_to_affine(tracker.problem)
            LinearAlgebra.normalize!(state.solution)
        end

        state.solution_cond = condition_jacobian(tracker)
        # If cond is not so high, maybe we don't have a singular solution in the end?
        if state.winding_number == 1 && state.solution_cond < 1e8
            result = correct!(core_tracker.state.x̄, core_tracker, state.solution, Inf;
                use_qr=true, max_iters=3,
                accuracy=core_tracker.options.refinement_accuracy)
            if isconverged(result)
                state.winding_number = 0
                @goto non_singular_case
            end
        elseif state.winding_number > 1 && state.solution_cond < 1e10 && residual(tracker) > 100
            state.status = PathTrackerStatus.post_check_failed
        end

        # In the case of a squared up system we now have to get rid of the
        # excess solutions, but since we don't have Newton's method at hand
        # we simply rely non the residual.
        if is_squared_up_system(core_tracker.homotopy) &&
           residual(tracker) > options.overdetermined_min_residual
             state.status = PathTrackerStatus.excess_solution
        end
    # 1.a) + 2.a)
    else
        # First we refine the obtained solution if possible
        result = correct!(core_tracker.state.x̄, core_tracker, currx(core_tracker), Inf;
            use_qr=true, max_iters=3,
            accuracy=core_tracker.options.refinement_accuracy)
        @label non_singular_case

        state.solution_cond = core_tracker.state.jacobian.cond

        if isconverged(result) && core_tracker.state.jacobian.corank_proposal == 0
            state.solution .= core_tracker.state.x̄
        elseif (!isconverged(result) || core_tracker.state.jacobian.corank_proposal > 0) &&
               at_infinity_post_check(state.val, options)

            state.solution .= currx(core_tracker)
            state.status = PathTrackerStatus.at_infinity
            return nothing
        else
            state.solution .= currx(core_tracker)
            state.status = PathTrackerStatus.post_check_failed
            return nothing
        end

        # We want to check that if we present an affine solution to a user
        # that this is still in a convergent region of the affine target system
        # This covers that we have a squared up system and when we tracked in projective space
        if pull_back_is_to_affine(tracker.problem) &&
            (is_squared_up_system(core_tracker.homotopy) ||
             !affine_tracking(core_tracker))

            # We start with bringing it on the affine patch
            if !affine_tracking(core_tracker)
                ProjectiveVectors.affine_chart!(state.solution)
            end
            init_auto_scaling!(cache.weighted_ip, state.solution, AutoScalingOptions())
            target_result = newton!(cache.base_point, cache.target_system, state.solution,
                            cache.weighted_ip, cache.target_newton_cache;
                            use_qr=true,
                            tol=core_tracker.options.refinement_accuracy, miniters=1, maxiters=2)
            if isconverged(target_result)
                state.solution .= cache.base_point
            elseif is_squared_up_system(core_tracker.homotopy)
                state.status = PathTrackerStatus.excess_solution
            else
                state.status = PathTrackerStatus.post_check_failed
            end
            state.solution_accuracy = target_result.accuracy

            if !pull_back_is_to_affine(tracker.problem)
                LinearAlgebra.normalize!(state.solution)
                changepatch!(cache.target_system.patch, state.solution)
                state.solution_accuracy = result.accuracy
            else
                state.solution_accuracy = result.accuracy
            end
        # We have a purely projective solution. Here we only have to handle the case
        # of overdetermined systems.
        elseif !pull_back_is_to_affine(tracker.problem) &&
                is_squared_up_system(core_tracker.homotopy)

            LinearAlgebra.normalize!(state.solution)
            changepatch!(cache.target_system.patch, state.solution)
            target_proj_result = newton!(cache.base_point, cache.target_system, state.solution,
                            euclidean_norm, cache.target_newton_cache;
                            use_qr=true,
                            tol=core_tracker.options.refinement_accuracy, miniters=1, maxiters=2)
            if isconverged(target_proj_result)
                state.solution .= cache.base_point
                state.solution_accuracy = result.accuracy
                LinearAlgebra.normalize!(state.solution)
            else
                state.status = PathTrackerStatus.excess_solution
            end

        end
    end
end

#############################
# Convencience constructors #
#############################

const pathtracker_startsolutions_supported_keywords = [
    problem_startsolutions_supported_keywords;
    coretracker_supported_keywords;
    pathtracker_supported_keywords]

"""
    pathtracker_startsolutions(args...; kwargs...)

Construct a [`PathTracker`](@ref) and start solutions in the same way [`solve`](@ref) does it.
This also takes the same input arguments as `solve`. This is convenient if you want
to investigate single paths.
"""
function pathtracker_startsolutions(args...; kwargs...)
    invalid = invalid_kwargs(kwargs, pathtracker_startsolutions_supported_keywords)
    check_kwargs_empty(invalid, pathtracker_startsolutions_supported_keywords)
    supported, rest = splitkwargs(kwargs, problem_startsolutions_supported_keywords)
    prob, startsolutions = problem_startsolutions(args...; supported...)
    tracker_startsolutions(prob, startsolutions; rest...)
end

function tracker_startsolutions(prob::Problem, startsolutions; kwargs...)
    tracker = PathTracker(prob, start_solution_sample(startsolutions); kwargs...)
    (tracker=tracker, startsolutions=startsolutions)
end

"""
    pathtracker(args...; kwargs...)

Construct a [`PathTracker`](@ref) in the same way [`solve`](@ref) does it.
This also takes the same input arguments as `solve` with the exception that you do not need to specify startsolutions.

## Examples

### Obtain single solution
We want to construct a path tracker to track a parameterized system `f` with parameters `p`
from the parameters `a` to `b`.
```julia
tracker = pathtracker(f, parameters=p, p₁=a, p₀=b)
```
You then can obtain a single solution at `b` by using
```julia
x_b = solution(track(tracker, x_a))
```
"""
pathtracker(args...; kwargs...) = first(pathtracker_startsolutions(args...; kwargs...))

############
## Result ##
############

solution(tracker::PathTracker) = pull_back(tracker.problem, tracker.state.solution)

function winding_number(tracker::PathTracker)
    tracker.state.winding_number == 0 ? nothing : tracker.state.winding_number
end

"""
    PathResult{V<:AbstractVector}

A `PathResult` is the result of tracking of a path using [`PathTracker`](@ref).
Its fields are

* `return_code`: One of `:success`, `:at_infinity` or any error code in [`PathTrackerStatus.states`](@ref) converted to a `Symbol`.
* `solution::V`: The solution vector.
* `t::Float64`: The value of `t` at which `solution` was computed. Note that if `return_code` is `:at_infinity`, then `t` is the value when this was decided.
* `accuracy::Union{Nothing, Float64}`: An approximation of ``||x-x^*||₂`` where ``x`` is the computed solution and ``x^*`` is the true solution.
* `residual::Union{Nothing, Float64}`: The value of the 2-norm of `H(solution, 0)`.
* `condition_jacobian::Union{Nothing, Float64}`: This is the condition number of the row-equilibrated Jacobian at the solution. A high condition number indicates a singularity.
* `winding_number:Union{Nothing, Int}`: The estimated winding number. This is a lower bound on the multiplicity of the solution.
* `start_solution::Union{Nothing, Int}`: The start solution of the path.
* `accepted_steps::Int`: The number of accepted steps during the path tracking.
* `rejected_steps::Int`: The number of rejected steps during the path tracking.
* `valuation::Union{Nothing, Vector{Float64}}`: An approximation of the valuation of the Puiseux series expansion of `x(t)`.
* `valuation_accuracy::Union{Nothing, Vector{Float64}}`: An estimate of the accuracy of the valuation of the Puiseux series expansion of `x(t)`.

     PathResult(tracker::PathTracker, start_solution=nothing, path_number::Union{Nothing,Int}=nothing; details=:default)

Possible `details` values are `:minimal` (minimal details), `:default` (default) and `:extensive` (all information possible).
"""
struct PathResult{V<:AbstractVector}
    return_code::Symbol
    solution::V
    t::Float64
    accuracy::Union{Nothing, Float64}
    residual::Union{Nothing, Float64} # level 1+
    condition_jacobian::Union{Nothing, Float64}
    winding_number::Union{Nothing, Int}
    path_number::Union{Nothing, Int}
    start_solution::Union{Nothing, V} # level 1+
    # performance stats
    accepted_steps::Int
    rejected_steps::Int
    valuation::Union{Nothing, Vector{Float64}} # level 2+
    valuation_accuracy::Union{Nothing, Vector{Float64}} # level 2+
end

function PathResult(tracker::PathTracker, start_solution, path_number::Union{Nothing,Int}=nothing; details::Symbol=:default)
    @unpack state, core_tracker, cache = tracker
    details_level = detailslevel(details)
    return_code = Symbol(state.status)
    windingnumber = winding_number(tracker)
    x = solution(tracker)
    # accuracy
    if !isnan(state.solution_accuracy)
        accuracy = state.solution_accuracy
    else
        accuracy = nothing
    end

    if return_code == :success
        t = 0.0
    else
        t = exp(-state.s)
    end
    # condition
    if isnan(state.solution_cond)
        condition_jac = cond(core_tracker)
    else
        condition_jac = state.solution_cond
    end
    # residual
    if return_code == :success && details_level ≥ 1
        res = residual(tracker)
    else
        res = nothing
    end

    if details_level ≥ 1
        # mimic the behaviour in track! to get a start solution of the same type as x
        embed!(cache.base_point, tracker.problem, start_solution)
        startsolution = pull_back(tracker.problem, cache.base_point; regauge=false)
    else
        startsolution = nothing
    end

    accepted_steps = core_tracker.state.accepted_steps
    rejected_steps = core_tracker.state.rejected_steps

    if details_level == 2
        valuation = copy(tracker.state.val.v)
        valuation_accuracy = abs.(tracker.state.val.v̇)
    else
        valuation = nothing
        valuation_accuracy = nothing
    end

    PathResult(return_code, x, t, accuracy, res, condition_jac,
               windingnumber, path_number,
               startsolution, accepted_steps, rejected_steps,
               valuation, valuation_accuracy)
end

function detailslevel(details::Symbol)
    if details == :minimal
        return 0
    elseif details == :extensive
        return 2
    else
        return 1
    end
end

function residual(tracker::PathTracker, x=tracker.state.solution)
    evaluate!(tracker.cache.target_residual, tracker.cache.target_system, x,
              tracker.cache.target_newton_cache.system_cache)
    euclidean_norm(tracker.cache.target_residual)
end

function condition_jacobian(tracker::PathTracker, x=tracker.state.solution)
    jac = tracker.cache.target_jacobian
    jacobian!(jac.J, tracker.cache.target_system, x,
              tracker.cache.target_newton_cache.system_cache)
    updated_jacobian!(jac; update_infos=true)
    jac.cond
end

"""
    result_type(tracker::PathTracker)

Returns the type of result `track` will return.
"""
result_type(tracker::PathTracker) = PathResult{typeof(solution(tracker))}

function Base.show(io::IO, r::PathResult)
    iscompact = get(io, :compact, false)
    if iscompact || haskey(io, :typeinfo)
        println(io, "• return_code: $(r.return_code)")
        if r.return_code != PathTrackerStatus.success
            println(io, " • t: $(r.t)")
        end
        println(io, " • solution: ", r.solution)
        r.accuracy !== nothing &&
            println(io, " • accuracy: $(Printf.@sprintf "%.3e" r.accuracy)")
        r.winding_number !== nothing &&
            println(io, " • winding_number: $(r.winding_number)")
        r.path_number !== nothing &&
            println(io, " • path_number: ", r.path_number)
    else
        println(io, "PathResult")
        println(io, "=================")
        println(io, " • return_code: $(r.return_code)")
        if r.return_code != PathTrackerStatus.success
            println(io, " • t: $(r.t)")
        end
        println(io, " • solution: ", r.solution)
        r.accuracy !== nothing &&
            println(io, " • accuracy: $(Printf.@sprintf "%.3e" r.accuracy)")
        r.residual !== nothing &&
            println(io, " • residual: $(Printf.@sprintf "%.3e" r.residual)")
        r.winding_number !== nothing &&
            println(io, " • winding_number: $(r.winding_number)")
        r.condition_jacobian !== nothing &&
            println(io, " • condition_jacobian: $(Printf.@sprintf "%.3e" r.condition_jacobian)")
        r.path_number !== nothing &&
            println(io, " • path_number: ", r.path_number)
    end
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::PathResult) = x

"""
    solution(pathresult)

Get the solution of the path.
"""
solution(r::PathResult) = r.solution

"""
    accuracy(pathresult)

Get the accuracy of the solution ``x`` of the path, i.e., ``||H(x, 0)||₂``.
"""
accuracy(r::PathResult) = r.accuracy

"""
    residual(pathresult)

Get the residual of the solution ``x`` of the path, i.e., ``||H(x, 0)||₂``.
"""
residual(r::PathResult) = r.residual

"""
    start_solution(pathresult)

Get the start solution of the solution ``x`` of the path.
"""
start_solution(r::PathResult) = r.start_solution

"""
    issuccess(pathresult)

Checks whether the path is successfull.
"""
issuccess(r::PathResult) = r.return_code == :success

"""
    isfailed(pathresult)

Checks whether the path failed.
"""
isfailed(r::PathResult) =!(r.return_code == :at_infinity || r.return_code == :success)

"""
    isatinfinity(pathresult)

Checks whether the path goes to infinity.
"""
isatinfinity(r::PathResult) = r.return_code == :at_infinity

"""
    isfinite(pathresult)

Checks whether the path result is finite.
"""
Base.isfinite(r::PathResult) = r.return_code == :success # we don't check isaffine to make other code easier

"""
    issingular(pathresult; tol=1e10)

Checks whether the path result is singular. This is true if
the winding number is larger than  1 or if the condition number of the Jacobian
is larger than `tol`.
"""
issingular(r::PathResult; tol=1e10) = issingular(r, tol)
function issingular(r::PathResult, tol::Real)
    (unpack(r.winding_number, 0) ≥ 1 || unpack(r.condition_jacobian, 1.0) > tol) && LinearAlgebra.issuccess(r)
end

"""
    isnonsingular(pathresult; tol=1e10)

Checks whether the path result is non-singular. This is true if
it is not singular.
"""
isnonsingular(r::PathResult; kwargs...) = !issingular(r; kwargs...) && LinearAlgebra.issuccess(r)
isnonsingular(r::PathResult, tol::Real) = !issingular(r, tol) && LinearAlgebra.issuccess(r)

"""
    isreal(pathresult; tol=1e-6)

We consider a result as `real` if the 2-norm of the imaginary part of the solution is at most `tol`.
"""
Base.isreal(r::PathResult; tol=1e-6) = isreal(r, tol)
Base.isreal(r::PathResult, tol::Real) = isrealvector(r.solution, tol)

"""
    isprojective(pathresult)

Return`s true if the solution is a projective vector.
"""
isprojective(r::PathResult{<:PVector}) = true
isprojective(r::PathResult) = false

"""
    isaffine(pathresult)

Return`s true if the solution is an affine vector.
"""
isaffine(r::PathResult) = !isprojective(r)


text:
import ACE.SphericalHarmonics: SHBasis, index_y

@doc raw"""
`struct Ylm1pBasis <: OneParticleBasis`

One-particle basis of the form
```math
\phi_{lm}({\bm r}) = Y_l^m(\hat{\br r})
```
Fundamental building block of ACE basis sets of the form
```math
   R_{nl}^{\mu_i, \mu_j}(r_{ij}) Y_l^m(\hat{\bm r})
```
This type basically just translates the `SHBasis` into a valid one-particle
basis.
"""
mutable struct Ylm1pBasis{T, VSYM, LSYM, MSYM} <: OneParticleBasis{Complex{T}}
   SH::SHBasis{T}  # SH = Ylm
end

# ---------------------- Implementation of Ylm1pBasis

Ylm1pBasis(maxL::Integer, T = Float64; kwargs...) =
      Ylm1pBasis((SHBasis(maxL, T)); kwargs...)

Ylm1pBasis(SH::SHBasis{T}; varsym = :rr, lsym = :l, msym = :m)  where {T} =
   Ylm1pBasis{T, varsym, lsym, msym}(SH)

Base.length(basis::Ylm1pBasis) = length(basis.SH)

_varsym(::Ylm1pBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = VSYM
_lsym(::Ylm1pBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = LSYM
_msym(::Ylm1pBasis{T, VSYM, LSYM, MSYM}) where {T, VSYM, LSYM, MSYM} = MSYM
_l(b, basis::Ylm1pBasis) = getproperty(b, _lsym(basis))
_m(b, basis::Ylm1pBasis) = getproperty(b, _msym(basis))
_lm(b, basis::Ylm1pBasis) = _l(b, basis), _m(b, basis)

# -> TODO : figure out how to do this well!!!
# Base.rand(basis::Ylm1pBasis) =
#       AtomState(rand(basis.zlist.list), ACE.Random.rand_vec(basis.J))

function get_spec(basis::Ylm1pBasis{T, VS, L, M}) where {T, VS, L, M}
   @assert length(basis) == (basis.SH.maxL + 1)^2
   spec = Vector{NamedTuple{(L, M), Tuple{Int, Int}}}(undef, length(basis))
   for l = 0:basis.SH.maxL, m = -l:l
      spec[index_y(l, m)] = NamedTuple{(L, M), Tuple{Int, Int}}( (l, m) )
   end
   return spec
end

==(P1::Ylm1pBasis, P2::Ylm1pBasis) =  ACE._allfieldsequal(P1, P2)

write_dict(basis::Ylm1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_Ylm1pBasis",
          "SH" => write_dict(basis.SH),
          "varsym" => _varsym(basis),
          "lsym" => _lsym(basis),
          "msym" => _msym(basis) )

read_dict(::Val{:ACE_Ylm1pBasis}, D::Dict) =
      Ylm1pBasis(read_dict(D["SH"]),
                 varsym = D["varsym"],
                 lsym = D["lsym"],
                 msym = D["msym"] )

fltype(basis::Ylm1pBasis{T}) where T = Complex{T}
rfltype(basis::Ylm1pBasis{T}) where T = T

gradtype(::Ylm1pBasis{T}) where {T} = SVector{3, Complex{T}}

symbols(basis::Ylm1pBasis) = [_lsym(basis), _msym(basis)]

function indexrange(basis::Ylm1pBasis)
   maxl = basis.SH.maxL
   # note we create a stupid tensor product domain and then make up for it
   # by using an index filter during the basis generation process
   return NamedTuple{(_lsym(basis), _msym(basis))}( (0:maxl, -maxl:maxl) )
end

function isadmissible(b, basis::Ylm1pBasis)
   l, m = _lm(b, basis)
   return (0 <= l <= basis.SH.maxL) && (-l <= m <= l)
end

# ---------------------------  Evaluation code
#

alloc_B(basis::Ylm1pBasis, args...) = alloc_B(basis.SH)

alloc_dB(basis::Ylm1pBasis, args...) = alloc_dB(basis.SH)

alloc_temp(basis::Ylm1pBasis, args...) = alloc_temp(basis.SH)

alloc_temp_d(basis::Ylm1pBasis, args...) = alloc_temp_d(basis.SH, args...)

_rr(X, basis::Ylm1pBasis) = getproperty(X, _varsym(basis))

evaluate!(B, tmp, basis::Ylm1pBasis, X::AbstractState) =
      evaluate!(B, tmp, basis.SH, _rr(X, basis))

evaluate_d!(dB, tmpd, basis::Ylm1pBasis, X::AbstractState) =
      evaluate_d!(dB, tmpd, basis.SH, _rr(X, basis))

evaluate_ed!(B, dB, tmpd, basis::Ylm1pBasis, X::AbstractState) =
      evaluate_ed!(B, dB, tmpd, basis.SH, _rr(X, basis))

degree(b, Ylm::Ylm1pBasis) = _l(b, Ylm)

get_index(Ylm::Ylm1pBasis, b) = index_y(_l(b, Ylm), _m(b, Ylm))

#
# function add_into_A_dA!(A, dA, tmpd, basis::RnYlm1pBasis, R, iz::Integer, iz0::Integer)
#    r = norm(R)
#    R̂ = R / r
#    # evaluate the r-basis and the R̂-basis for the current neighbour at R
#    evaluate_d!(tmpd.BJ, tmpd.dBJ, tmpd.tmpdJ, basis.J, r)
#    evaluate_d!(tmpd.BY, tmpd.dBY, tmpd.tmpdY, basis.SH, R)
#    # add the contributions to the A_zklm, ∇A
#    @inbounds for (i, nlm) in enumerate(basis.spec)
#       iY = index_y(nlm.l, nlm.m)
#       A[i] += tmpd.BJ[nlm.n] * tmpd.BY[iY]
#       dA[i] = (tmpd.dBJ[nlm.n] * tmpd.BY[iY]) * R̂ + tmpd.BJ[nlm.n] * tmpd.dBY[iY]
#    end
#    return nothing
# end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((e+f*x^2)^4*(a+b*x^2)*(c+d*x^2), x) == :(1//3*e^3*x^3*(a*d*e+b*c*e+4*a*c*f)+1//5*e^2*x^5*(b*e*(d*e+4*c*f)+2*a*f*(2*d*e+3*c*f))+1//9*f^2*x^9*(a*f*(c*f+4*d*e)+2*b*e*(2*c*f+3*d*e))+1//11*f^3*x^11*(a*d*f+b*c*f+4*b*d*e)+a*c*x*e^4+1//13*b*d*f^4*x^13+2//7*e*f*x^7*(a*f*(2*c*f+3*d*e)+b*e*(2*d*e+3*c*f)))
@test integrate((e+f*x^2)^3*(a+b*x^2)*(c+d*x^2), x) == :(1//3*e^2*x^3*(a*d*e+b*c*e+3*a*c*f)+1//5*e*x^5*(b*e*(d*e+3*c*f)+3*a*f*(c*f+d*e))+1//7*f*x^7*(a*f*(c*f+3*d*e)+3*b*e*(c*f+d*e))+1//9*f^2*x^9*(a*d*f+b*c*f+3*b*d*e)+a*c*x*e^3+1//11*b*d*f^3*x^11)
@test integrate((e+f*x^2)^2*(a+b*x^2)*(c+d*x^2), x) == :(x^5*(1//5*a*f*(c*f+2*d*e)+1//5*b*e*(d*e+2*c*f))+1//3*e*x^3*(a*d*e+b*c*e+2*a*c*f)+1//7*f*x^7*(a*d*f+b*c*f+2*b*d*e)+a*c*x*e^2+1//9*b*d*f^2*x^9)
@test integrate((a+b*x^2)*(c+d*x^2)*(e+f*x^2), x) == :(x^3*(1//3*a*c*f+1//3*a*d*e+1//3*b*c*e)+x^5*(1//5*a*d*f+1//5*b*c*f+1//5*b*d*e)+a*c*e*x+1//7*b*d*f*x^7)
@test integrate((e+f*x^2)^-1*(a+b*x^2)*(c+d*x^2), x) == :(-1//3*x*f^-2*(-3*b*c*f+-2*a*d*f+3*b*d*e)+1//3*d*x*f^-1*(a+b*x^2)+e^-1//2*f^-5//2*(b*e+-1*a*f)*(d*e+-1*c*f)*arctan(x*e^-1//2*f^(1/2)))
@test integrate((e+f*x^2)^-2*(a+b*x^2)*(c+d*x^2), x) == :(-1//2*e^-3//2*f^-5//2*(b*e*(-1*c*f+3*d*e)+-1*a*f*(c*f+d*e))*arctan(x*e^-1//2*f^(1/2))+(1/2)*b*x*e^-1*f^-2*(-1*c*f+3*d*e)+-1//2*x*e^-1*f^-1*(e+f*x^2)^-1*(a+b*x^2)*(d*e+-1*c*f))
@test integrate((e+f*x^2)^-3*(a+b*x^2)*(c+d*x^2), x) == :(1//8*e^-5//2*f^-5//2*(a*f*(d*e+3*c*f)+b*e*(c*f+3*d*e))*arctan(x*e^-1//2*f^(1/2))+-1//8*x*e^-2*f^-2*(e+f*x^2)^-1*(b*e*(c*f+3*d*e)+-1*a*f*(d*e+3*c*f))+-1//4*x*e^-1*f^-1*(e+f*x^2)^-2*(a+b*x^2)*(d*e+-1*c*f))
@test integrate((e+f*x^2)^-4*(a+b*x^2)*(c+d*x^2), x) == :(1//16*e^-7//2*f^-5//2*(a*f*(d*e+5*c*f)+b*e*(c*f+d*e))*arctan(x*e^-1//2*f^(1/2))+-1//24*x*e^-2*f^-2*(e+f*x^2)^-2*(-1*a*f*(d*e+5*c*f)+3*b*e*(c*f+d*e))+1//16*x*e^-3*f^-2*(e+f*x^2)^-1*(a*f*(d*e+5*c*f)+b*e*(c*f+d*e))+-1//6*x*e^-1*f^-1*(e+f*x^2)^-3*(a+b*x^2)*(d*e+-1*c*f))
@test integrate((c+d*x^2)^2*(e+f*x^2)^3*(a+b*x^2), x) == :(x^7*(1//7*a*f*(c^2*f^2+3*d^2*e^2+6*c*d*e*f)+1//7*b*e*(d^2*e^2+3*c^2*f^2+6*c*d*e*f))+1//5*e*x^5*(a*(d^2*e^2+3*c^2*f^2+6*c*d*e*f)+b*c*e*(2*d*e+3*c*f))+1//9*f*x^9*(b*(c^2*f^2+3*d^2*e^2+6*c*d*e*f)+a*d*f*(2*c*f+3*d*e))+a*x*c^2*e^3+1//3*c*e^2*x^3*(b*c*e+2*a*d*e+3*a*c*f)+1//11*d*f^2*x^11*(a*d*f+2*b*c*f+3*b*d*e)+1//13*b*d^2*f^3*x^13)
@test integrate((c+d*x^2)^2*(e+f*x^2)^2*(a+b*x^2), x) == :(x^5*(1//5*a*(c^2*f^2+d^2*e^2+4*c*d*e*f)+2//5*b*c*e*(c*f+d*e))+x^7*(1//7*b*(c^2*f^2+d^2*e^2+4*c*d*e*f)+2//7*a*d*f*(c*f+d*e))+a*x*c^2*e^2+1//3*c*e*x^3*(2*a*(c*f+d*e)+b*c*e)+1//9*d*f*x^9*(2*b*(c*f+d*e)+a*d*f)+1//11*b*d^2*f^2*x^11)
@test integrate((c+d*x^2)^2*(a+b*x^2)*(e+f*x^2), x) == :(x^5*(1//5*a*d*(d*e+2*c*f)+1//5*b*c*(c*f+2*d*e))+1//3*c*x^3*(a*c*f+b*c*e+2*a*d*e)+1//7*d*x^7*(a*d*f+b*d*e+2*b*c*f)+a*e*x*c^2+1//9*b*f*d^2*x^9)
@test integrate((c+d*x^2)^2*(e+f*x^2)^-1*(a+b*x^2), x) == :(-1//15*x*f^-3*(-1*b*(8*c^2*f^2+15*d^2*e^2+-25*c*d*e*f)+5*a*d*f*(-5*c*f+3*d*e))+-1//15*x*f^-2*(c+d*x^2)*(-5*a*d*f+-4*b*c*f+5*b*d*e)+1//5*b*x*f^-1*(c+d*x^2)^2+-1*e^-1//2*f^-7//2*(d*e+-1*c*f)^2*(b*e+-1*a*f)*arctan(x*e^-1//2*f^(1/2)))
@test integrate((c+d*x^2)^2*(e+f*x^2)^-2*(a+b*x^2), x) == :((1/2)*e^-3//2*f^-7//2*(d*e+-1*c*f)*(b*e*(-1*c*f+5*d*e)+-1*a*f*(c*f+3*d*e))*arctan(x*e^-1//2*f^(1/2))+-1//6*d*x*e^-1*f^-3*(b*e*(-13*c*f+15*d*e)+-3*a*f*(-1*c*f+3*d*e))+-1//2*x*e^-1*f^-1*(c+d*x^2)^2*(e+f*x^2)^-1*(b*e+-1*a*f)+1//6*d*x*e^-1*f^-2*(c+d*x^2)*(-3*a*f+5*b*e))
@test integrate((c+d*x^2)^2*(e+f*x^2)^-3*(a+b*x^2), x) == :(-1//8*e^-5//2*f^-7//2*(b*e*(-1*c^2*f^2+15*d^2*e^2+-6*c*d*e*f)+-1*a*f*(3*c^2*f^2+3*d^2*e^2+2*c*d*e*f))*arctan(x*e^-1//2*f^(1/2))+1//8*d*x*e^-2*f^-3*(b*e*(-1*c*f+15*d*e)+-3*a*f*(c*f+d*e))+-1//4*x*e^-1*f^-1*(c+d*x^2)^2*(e+f*x^2)^-2*(b*e+-1*a*f)+-1//8*x*e^-2*f^-2*(e+f*x^2)^-1*(c+d*x^2)*(b*e*(-1*c*f+5*d*e)+-1*a*f*(d*e+3*c*f)))
@test integrate((c+d*x^2)^2*(e+f*x^2)^-4*(a+b*x^2), x) == :(1//16*e^-7//2*f^-7//2*(a*f*(d^2*e^2+5*c^2*f^2+2*c*d*e*f)+b*e*(c^2*f^2+5*d^2*e^2+2*c*d*e*f))*arctan(x*e^-1//2*f^(1/2))+-1//48*x*e^-3*f^-3*(e+f*x^2)^-1*(a*f*(-15*c^2*f^2+3*d^2*e^2+4*c*d*e*f)+b*e*(-3*c^2*f^2+15*d^2*e^2+-4*c*d*e*f))+-1//6*x*e^-1*f^-1*(c+d*x^2)^2*(e+f*x^2)^-3*(b*e+-1*a*f)+-1//24*x*e^-2*f^-2*(e+f*x^2)^-2*(c+d*x^2)*(d*e*(a*f+5*b*e)+-1*c*f*(b*e+5*a*f)))
@test integrate((c+d*x^2)^3*(e+f*x^2)^3*(a+b*x^2), x) == :(x^7*(1//7*a*(c^3*f^3+d^3*e^3+9*c*f*d^2*e^2+9*d*e*c^2*f^2)+3//7*b*c*e*(c^2*f^2+d^2*e^2+3*c*d*e*f))+x^9*(1//9*b*(c^3*f^3+d^3*e^3+9*c*f*d^2*e^2+9*d*e*c^2*f^2)+1//3*a*d*f*(c^2*f^2+d^2*e^2+3*c*d*e*f))+a*x*c^3*e^3+1//3*c^2*e^2*x^3*(3*a*(c*f+d*e)+b*c*e)+1//13*d^2*f^2*x^13*(3*b*(c*f+d*e)+a*d*f)+1//15*b*d^3*f^3*x^15+3//5*c*e*x^5*(a*(c^2*f^2+d^2*e^2+3*c*d*e*f)+b*c*e*(c*f+d*e))+3//11*d*f*x^11*(b*(c^2*f^2+d^2*e^2+3*c*d*e*f)+a*d*f*(c*f+d*e)))
@test integrate((c+d*x^2)^3*(e+f*x^2)^2*(a+b*x^2), x) == :(x^7*(1//7*a*d*(d^2*e^2+3*c^2*f^2+6*c*d*e*f)+1//7*b*c*(c^2*f^2+3*d^2*e^2+6*c*d*e*f))+1//5*c*x^5*(a*(c^2*f^2+3*d^2*e^2+6*c*d*e*f)+b*c*e*(2*c*f+3*d*e))+1//9*d*x^9*(b*(d^2*e^2+3*c^2*f^2+6*c*d*e*f)+a*d*f*(2*d*e+3*c*f))+a*x*c^3*e^2+1//3*e*c^2*x^3*(b*c*e+2*a*c*f+3*a*d*e)+1//11*f*d^2*x^11*(a*d*f+2*b*d*e+3*b*c*f)+1//13*b*d^3*f^2*x^13)
@test integrate((c+d*x^2)^3*(a+b*x^2)*(e+f*x^2), x) == :(1//3*c^2*x^3*(a*c*f+b*c*e+3*a*d*e)+1//5*c*x^5*(b*c*(c*f+3*d*e)+3*a*d*(c*f+d*e))+1//7*d*x^7*(a*d*(d*e+3*c*f)+3*b*c*(c*f+d*e))+1//9*d^2*x^9*(a*d*f+b*d*e+3*b*c*f)+a*e*x*c^3+1//11*b*f*d^3*x^11)
@test integrate((c+d*x^2)^3*(e+f*x^2)^-1*(a+b*x^2), x) == :(1//105*x*f^-4*(-1*b*(-48*c^3*f^3+105*d^3*e^3+-280*c*f*d^2*e^2+231*d*e*c^2*f^2)+7*a*d*f*(15*d^2*e^2+33*c^2*f^2+-40*c*d*e*f))+-1//35*x*f^-2*(c+d*x^2)^2*(-7*a*d*f+-6*b*c*f+7*b*d*e)+-1//105*x*f^-3*(c+d*x^2)*(-1*b*(24*c^2*f^2+35*d^2*e^2+-63*c*d*e*f)+7*a*d*f*(-9*c*f+5*d*e))+1//7*b*x*f^-1*(c+d*x^2)^3+e^-1//2*f^-9//2*(d*e+-1*c*f)^3*(b*e+-1*a*f)*arctan(x*e^-1//2*f^(1/2)))
@test integrate((c+d*x^2)^3*(e+f*x^2)^-2*(a+b*x^2), x) == :(-1//2*e^-3//2*f^-9//2*(d*e+-1*c*f)^2*(b*e*(-1*c*f+7*d*e)+-1*a*f*(c*f+5*d*e))*arctan(x*e^-1//2*f^(1/2))+-1//30*d*x*e^-1*f^-4*(-1*b*e*(81*c^2*f^2+105*d^2*e^2+-190*c*d*e*f)+5*a*f*(3*c^2*f^2+15*d^2*e^2+-22*c*d*e*f))+-1//2*x*e^-1*f^-1*(c+d*x^2)^3*(e+f*x^2)^-1*(b*e+-1*a*f)+-1//30*d*x*e^-1*f^-3*(c+d*x^2)*(b*e*(-33*c*f+35*d*e)+-5*a*f*(-3*c*f+5*d*e))+1//10*d*x*e^-1*f^-2*(c+d*x^2)^2*(-5*a*f+7*b*e))
@test integrate((c+d*x^2)^3*(e+f*x^2)^-3*(a+b*x^2), x) == :(1//8*e^-5//2*f^-9//2*(d*e+-1*c*f)*(b*e*(-1*c^2*f^2+35*d^2*e^2+-10*c*d*e*f)+-3*a*f*(c^2*f^2+5*d^2*e^2+2*c*d*e*f))*arctan(x*e^-1//2*f^(1/2))+1//24*d*x*e^-2*f^-4*(-1*b*e*(3*c^2*f^2+105*d^2*e^2+-100*c*d*e*f)+3*a*f*(-3*c^2*f^2+15*d^2*e^2+-4*c*d*e*f))+-1//4*x*e^-1*f^-1*(c+d*x^2)^3*(e+f*x^2)^-2*(b*e+-1*a*f)+-1//8*x*e^-2*f^-2*(c+d*x^2)^2*(e+f*x^2)^-1*(b*e*(-1*c*f+7*d*e)+-3*a*f*(c*f+d*e))+1//24*d*x*e^-2*f^-3*(c+d*x^2)*(b*e*(-3*c*f+35*d*e)+-3*a*f*(3*c*f+5*d*e)))
@test integrate((c+d*x^2)^3*(e+f*x^2)^-4*(a+b*x^2), x) == :(-1//16*e^-7//2*f^-9//2*(b*e*(-1*c^3*f^3+35*d^3*e^3+-15*c*f*d^2*e^2+-3*d*e*c^2*f^2)+-1*a*f*(5*c^3*f^3+5*d^3*e^3+3*c*f*d^2*e^2+3*d*e*c^2*f^2))*arctan(x*e^-1//2*f^(1/2))+1//48*d*x*e^-3*f^-4*(b*e*(-3*c^2*f^2+105*d^2*e^2+-10*c*d*e*f)+-1*a*f*(15*c^2*f^2+15*d^2*e^2+14*c*d*e*f))+-1//6*x*e^-1*f^-1*(c+d*x^2)^3*(e+f*x^2)^-3*(b*e+-1*a*f)+-1//24*x*e^-2*f^-2*(c+d*x^2)^2*(e+f*x^2)^-2*(b*e*(-1*c*f+7*d*e)+-1*a*f*(d*e+5*c*f))+-1//48*x*e^-3*f^-3*(e+f*x^2)^-1*(c+d*x^2)*(b*e*(-3*c^2*f^2+35*d^2*e^2+-8*c*d*e*f)+-1*a*f*(5*d^2*e^2+15*c^2*f^2+4*c*d*e*f)))
@test integrate((c+d*x^2)^3//2*(e+f*x^2)^(1/2)*(a+b*x^2), x) == :(1//7*b*x*d^-1*(c+d*x^2)^5//2*(e+f*x^2)^(1/2)+-1//105*x*d^-2*f^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-6*c^3*f^3+8*d^3*e^3+-19*c*f*d^2*e^2+9*d*e*c^2*f^2)+7*a*d*f*(-3*c^2*f^2+2*d^2*e^2+-7*c*d*e*f))+1//35*x*d^-1*f^-1*(c+d*x^2)^3//2*(e+f*x^2)^(1/2)*(b*d*e+-2*b*c*f+7*a*d*f)+1//105*x*d^-1*f^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(-1*b*(4*d^2*e^2+6*c^2*f^2+-6*c*d*e*f)+7*a*d*f*(d*e+3*c*f))+-1//105*d^-1*e^3//2*f^-5//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-3*c^2*f^2+4*d^2*e^2+-9*c*d*e*f)+7*a*d*f*(d*e+-9*c*f))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//105*d^-2*e^(1/2)*f^-5//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-6*c^3*f^3+8*d^3*e^3+-19*c*f*d^2*e^2+9*d*e*c^2*f^2)+7*a*d*f*(-3*c^2*f^2+2*d^2*e^2+-7*c*d*e*f))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(a+b*x^2), x) == :(1//5*b*x*d^-1*(c+d*x^2)^3//2*(e+f*x^2)^(1/2)+1//15*x*d^-1*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(b*d*e+-2*b*c*f+5*a*d*f)+1//15*x*d^-2*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-2*b*(c^2*f^2+d^2*e^2+-1*c*d*e*f)+5*a*d*f*(c*f+d*e))+-1//15*d^-1*e^3//2*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*c*f+b*d*e+-10*a*d*f)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1//15*d^-2*e^(1/2)*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-2*b*(c^2*f^2+d^2*e^2+-1*c*d*e*f)+5*a*d*f*(c*f+d*e))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(a+b*x^2), x) == :(1//3*b*x*d^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)+1//3*x*d^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*d*e+-2*b*c*f+3*a*d*f)+-1//3*d^-2*e^(1/2)*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*d*e+-2*b*c*f+3*a*d*f)*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1//3*c^-1*d^-1*e^3//2*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*c+-3*a*d)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-3//2*(e+f*x^2)^(1/2)*(a+b*x^2), x) == :(-1*x*c^-1*d^-1*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(b*c+-1*a*d)+f*x*c^-1*d^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*a*d+2*b*c)+b*c^-1*d^-1*e^3//2*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1*c^-1*d^-2*e^(1/2)*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*a*d+2*b*c)*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-5//2*(e+f*x^2)^(1/2)*(a+b*x^2), x) == :(-1//3*x*c^-1*d^-1*(c+d*x^2)^-3//2*(e+f*x^2)^(1/2)*(b*c+-1*a*d)+1//3*c^-3//2*d^-3//2*(c*e^-1*(c+d*x^2)^-1*(e+f*x^2))^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-1*(d*e*(b*c+2*a*d)+-1*c*f*(a*d+2*b*c))*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*c*f*d^-1*e^-1)+1//3*c^-2*d^-1*e^3//2*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-1*(b*c+-1*a*d)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-7//2*(e+f*x^2)^(1/2)*(a+b*x^2), x) == :(-1//5*x*c^-1*d^-1*(c+d*x^2)^-5//2*(e+f*x^2)^(1/2)*(b*c+-1*a*d)+1//15*x*c^-2*d^-1*(c+d*x^2)^-3//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-1*(a*d*(-3*c*f+4*d*e)+b*c*(d*e+-2*c*f))+1//15*c^-5//2*d^-3//2*(c*e^-1*(c+d*x^2)^-1*(e+f*x^2))^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-2*(a*d*(3*c^2*f^2+8*d^2*e^2+-13*c*d*e*f)+2*b*c*(c^2*f^2+d^2*e^2+-1*c*d*e*f))*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*c*f*d^-1*e^-1)+-1//15*c^-3*d^-1*e^3//2*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-2*(b*c*(c*f+d*e)+2*a*d*(-3*c*f+2*d*e))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^(1/2)*(e+f*x^2)^3//2*(a+b*x^2), x) == :(1//7*b*x*d^-1*(c+d*x^2)^3//2*(e+f*x^2)^3//2+1//35*x*d^-2*(c+d*x^2)^3//2*(e+f*x^2)^(1/2)*(-4*b*c*f+3*b*d*e+7*a*d*f)+1//105*x*d^-3*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-8*c^3*f^3+6*d^3*e^3+-9*c*f*d^2*e^2+19*d*e*c^2*f^2)+7*a*d*f*(-2*c^2*f^2+3*d^2*e^2+7*c*d*e*f))+1//105*x*d^-2*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(b*(3*d^2*e^2+8*c^2*f^2+-15*c*d*e*f)+14*a*d*f*(-1*c*f+3*d*e))+-1//105*d^-3*e^(1/2)*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-8*c^3*f^3+6*d^3*e^3+-9*c*f*d^2*e^2+19*d*e*c^2*f^2)+7*a*d*f*(-2*c^2*f^2+3*d^2*e^2+7*c*d*e*f))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//105*d^-2*e^3//2*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-4*c^2*f^2+3*d^2*e^2+9*c*d*e*f)+7*a*d*f*(-1*c*f+9*d*e))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-1//2*(e+f*x^2)^3//2*(a+b*x^2), x) == :(1//5*b*x*d^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^3//2+1//15*x*d^-3*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*(3*d^2*e^2+8*c^2*f^2+-13*c*d*e*f)+10*a*d*f*(-1*c*f+2*d*e))+1//15*x*d^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(-4*b*c*f+3*b*d*e+5*a*d*f)+-1//15*d^-3*e^(1/2)*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*(3*d^2*e^2+8*c^2*f^2+-13*c*d*e*f)+10*a*d*f*(-1*c*f+2*d*e))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//15*c^-1*d^-2*e^3//2*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-4*f*c^2+6*c*d*e)+5*a*d*(-1*c*f+3*d*e))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-3//2*(e+f*x^2)^3//2*(a+b*x^2), x) == :(-1*x*c^-1*d^-1*(c+d*x^2)^-1//2*(e+f*x^2)^3//2*(b*c+-1*a*d)+1//3*f*x*c^-1*d^-3*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*c*(-8*c*f+7*d*e)+-3*a*d*(d*e+-2*c*f))+1//3*f*x*c^-1*d^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(-3*a*d+4*b*c)+-1//3*c^-1*d^-3*e^(1/2)*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*c*(-8*c*f+7*d*e)+-3*a*d*(d*e+-2*c*f))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//3*c^-1*d^-2*e^3//2*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-4*b*c*f+3*a*d*f+3*b*d*e)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-5//2*(e+f*x^2)^3//2*(a+b*x^2), x) == :(-1//3*x*c^-1*d^-1*(c+d*x^2)^-3//2*(e+f*x^2)^3//2*(b*c+-1*a*d)+1//3*x*c^-2*d^-2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(a*d*(c*f+2*d*e)+b*c*(d*e+-4*c*f))+-1//3*f*x*c^-2*d^-3*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*c*(d*e+-8*c*f)+2*a*d*(c*f+d*e))+1//3*c^-2*d^-3*e^(1/2)*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*c*(d*e+-8*c*f)+2*a*d*(c*f+d*e))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//3*c^-2*d^-2*e^3//2*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*a*d+4*b*c)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-7//2*(e+f*x^2)^3//2*(a+b*x^2), x) == :(-1//5*x*c^-1*d^-1*(c+d*x^2)^-5//2*(e+f*x^2)^3//2*(b*c+-1*a*d)+1//15*x*c^-2*d^-2*(c+d*x^2)^-3//2*(e+f*x^2)^(1/2)*(d*e*(b*c+4*a*d)+-1*c*f*(a*d+4*b*c))+1//15*c^-5//2*d^-5//2*(c*e^-1*(c+d*x^2)^-1*(e+f*x^2))^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-1*(a*d*(-2*c^2*f^2+8*d^2*e^2+-3*c*d*e*f)+b*c*(-8*c^2*f^2+2*d^2*e^2+3*c*d*e*f))*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*c*f*d^-1*e^-1)+-1//15*c^-3*d^-2*e^3//2*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-1*(a*d*(-1*c*f+4*d*e)+b*c*(d*e+-4*c*f))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-9//2*(e+f*x^2)^3//2*(a+b*x^2), x) == :(-1//7*x*c^-1*d^-1*(c+d*x^2)^-7//2*(e+f*x^2)^3//2*(b*c+-1*a*d)+1//35*x*c^-2*d^-2*(c+d*x^2)^-5//2*(e+f*x^2)^(1/2)*(d*e*(b*c+6*a*d)+-1*c*f*(3*a*d+4*b*c))+1//105*x*c^-3*d^-2*(c+d*x^2)^-3//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-1*(b*c*(-8*c^2*f^2+4*d^2*e^2+c*d*e*f)+3*a*d*(-2*c^2*f^2+8*d^2*e^2+-5*c*d*e*f))+1//105*c^-7//2*d^-5//2*(c*e^-1*(c+d*x^2)^-1*(e+f*x^2))^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-2*(b*c*(8*c^3*f^3+8*d^3*e^3+-5*c*f*d^2*e^2+-5*d*e*c^2*f^2)+6*a*d*(c^3*f^3+8*d^3*e^3+-12*c*f*d^2*e^2+2*d*e*c^2*f^2))*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*c*f*d^-1*e^-1)+-1//105*c^-4*d^-2*e^3//2*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-2*(2*b*c*(2*c^2*f^2+2*d^2*e^2+-1*c*d*e*f)+3*a*d*(c^2*f^2+8*d^2*e^2+-11*c*d*e*f))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^5//2*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(-1//35*x*f^-2*(c+d*x^2)^3//2*(e+f*x^2)^(1/2)*(-7*a*d*f+-5*b*c*f+6*b*d*e)+-1//105*x*f^-3*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(-1*b*(15*c^2*f^2+24*d^2*e^2+-43*c*d*e*f)+28*a*d*f*(d*e+-2*c*f))+1//7*b*x*f^-1*(c+d*x^2)^5//2*(e+f*x^2)^(1/2)+1//105*x*d^-1*f^-3*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-15*c^3*f^3+48*d^3*e^3+-128*c*f*d^2*e^2+103*d*e*c^2*f^2)+7*a*d*f*(8*d^2*e^2+23*c^2*f^2+-23*c*d*e*f))+1//105*e^(1/2)*f^-7//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*e*(24*d^2*e^2+45*c^2*f^2+-61*c*d*e*f)+7*a*f*(4*d^2*e^2+15*c^2*f^2+-11*c*d*e*f))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1//105*d^-1*e^(1/2)*f^-7//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(-15*c^3*f^3+48*d^3*e^3+-128*c*f*d^2*e^2+103*d*e*c^2*f^2)+7*a*d*f*(8*d^2*e^2+23*c^2*f^2+-23*c*d*e*f))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^3//2*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(-1//15*x*f^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(-5*a*d*f+-3*b*c*f+4*b*d*e)+1//5*b*x*f^-1*(c+d*x^2)^3//2*(e+f*x^2)^(1/2)+-1//15*x*d^-1*f^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(3*c^2*f^2+8*d^2*e^2+-13*c*d*e*f)+10*a*d*f*(d*e+-2*c*f))+-1//15*e^(1/2)*f^-5//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(4*d*e^2+-6*c*e*f)+5*a*f*(d*e+-3*c*f))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//15*d^-1*e^(1/2)*f^-5//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(3*c^2*f^2+8*d^2*e^2+-13*c*d*e*f)+10*a*d*f*(d*e+-2*c*f))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(1//3*b*x*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)+-1//3*x*d^-1*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*c*f+-3*a*d*f+2*b*d*e)+-1//3*e^(1/2)*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*e+-3*a*f)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//3*d^-1*e^(1/2)*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*c*f+-3*a*d*f+2*b*d*e)*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-1//2*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(b*x*d^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2+a*c^-1*e^(1/2)*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1*b*d^-1*e^(1/2)*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-3//2*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(-1*c^-1//2*d^-1//2*(c*e^-1*(c+d*x^2)^-1*(e+f*x^2))^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-1*(b*c+-1*a*d)*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*c*f*d^-1*e^-1)+c^-1*e^(1/2)*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-1*(b*e+-1*a*f)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-5//2*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(-1//3*x*c^-1*(c+d*x^2)^-3//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-1*(b*c+-1*a*d)+1//3*c^-3//2*d^-1//2*(c*e^-1*(c+d*x^2)^-1*(e+f*x^2))^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-2*(b*c*(c*f+d*e)+2*a*d*(d*e+-2*c*f))*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*c*f*d^-1*e^-1)+-1//3*c^-2*e^(1/2)*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-2*(a*d*e+-3*a*c*f+2*b*c*e)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-7//2*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(-1//5*x*c^-1*(c+d*x^2)^-5//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-1*(b*c+-1*a*d)+1//15*x*c^-2*(c+d*x^2)^-3//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-2*(b*c*(d*e+3*c*f)+4*a*d*(d*e+-2*c*f))+1//15*c^-5//2*d^-1//2*(c*e^-1*(c+d*x^2)^-1*(e+f*x^2))^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^(1/2)*(d*e+-1*c*f)^-3*(a*d*(8*d^2*e^2+23*c^2*f^2+-23*c*d*e*f)+b*c*(-3*c^2*f^2+2*d^2*e^2+-7*c*d*e*f))*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*c*f*d^-1*e^-1)+-1//15*c^-3*e^(1/2)*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-3*(a*(4*d^2*e^2+15*c^2*f^2+-11*c*d*e*f)+b*c*e*(d*e+-9*c*f))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^5//2*(e+f*x^2)^-3//2*(a+b*x^2), x) == :(-1*x*e^-1*f^-1*(c+d*x^2)^5//2*(e+f*x^2)^-1//2*(b*e+-1*a*f)+-1//15*x*e^-1*f^-3*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-2*b*e*(19*c^2*f^2+24*d^2*e^2+-44*c*d*e*f)+5*a*f*(3*c^2*f^2+8*d^2*e^2+-13*c*d*e*f))+-1//15*d*x*e^-1*f^-3*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(b*e*(-23*c*f+24*d*e)+-5*a*f*(-3*c*f+4*d*e))+-1//15*e^(1/2)*f^-7//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*b*(15*c^2*f^2+24*d^2*e^2+-41*c*d*e*f)+10*a*d*f*(-3*c*f+2*d*e))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//5*d*x*e^-1*f^-2*(c+d*x^2)^3//2*(e+f*x^2)^(1/2)*(-5*a*f+6*b*e)+1//15*e^-1//2*f^-7//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-2*b*e*(19*c^2*f^2+24*d^2*e^2+-44*c*d*e*f)+5*a*f*(3*c^2*f^2+8*d^2*e^2+-13*c*d*e*f))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^3//2*(e+f*x^2)^-3//2*(a+b*x^2), x) == :(-1*x*e^-1*f^-1*(c+d*x^2)^3//2*(e+f*x^2)^-1//2*(b*e+-1*a*f)+-1//3*x*e^-1*f^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*e*(-7*c*f+8*d*e)+-3*a*f*(-1*c*f+2*d*e))+-1//3*e^(1/2)*f^-5//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-3*a*d*f+-3*b*c*f+4*b*d*e)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//3*d*x*e^-1*f^-2*(c+d*x^2)^(1/2)*(e+f*x^2)^(1/2)*(-3*a*f+4*b*e)+1//3*e^-1//2*f^-5//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*e*(-7*c*f+8*d*e)+-3*a*f*(-1*c*f+2*d*e))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^(1/2)*(e+f*x^2)^-3//2*(a+b*x^2), x) == :(x*e^-1*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*a*f+2*b*e)+-1*x*e^-1*f^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*e+-1*a*f)+b*e^(1/2)*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1*e^-1//2*f^-3//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(-1*a*f+2*b*e)*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-1//2*(e+f*x^2)^-3//2*(a+b*x^2), x) == :(e^-1//2*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-1*(b*e+-1*a*f)*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1*c^-1*e^(1/2)*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-1*(b*c+-1*a*d)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-3//2*(e+f*x^2)^-3//2*(a+b*x^2), x) == :(-1*x*c^-1*(c+d*x^2)^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-1*(b*c+-1*a*d)+c^-1*e^(1/2)*f^-1//2*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-2*(b*c*f+b*d*e+-2*a*d*f)*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+-1*c^-1*e^-1//2*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-2*(-1*a*c*f+-1*a*d*e+2*b*c*e)*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((c+d*x^2)^-5//2*(e+f*x^2)^-3//2*(a+b*x^2), x) == :(-1//3*x*c^-1*(c+d*x^2)^-3//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-1*(b*c+-1*a*d)+1//3*x*c^-2*(c+d*x^2)^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-2*(b*c*(d*e+3*c*f)+2*a*d*(d*e+-3*c*f))+-1//3*c^-2*e^(1/2)*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-3*(a*d*(d*e+-9*c*f)+b*c*(3*c*f+5*d*e))*Elliptic.F(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1)+1//3*c^-2*e^-1//2*f^(1/2)*(e*c^-1*(e+f*x^2)^-1*(c+d*x^2))^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(d*e+-1*c*f)^-3*(a*(-3*c^2*f^2+2*d^2*e^2+-7*c*d*e*f)+b*c*e*(d*e+7*c*f))*Elliptic.E(arctan(x*e^-1//2*f^(1/2)),1+-1*d*e*c^-1*f^-1))
@test integrate((a+b*x^2)^-1//2*(c+d*x^2)^-3//2*(e+f*x^2), x) == :(-1*c^-1//2*d^-1//2*(c*a^-1*(c+d*x^2)^-1*(a+b*x^2))^-1//2*(a+b*x^2)^(1/2)*(c+d*x^2)^-1//2*(b*c+-1*a*d)^-1*(d*e+-1*c*f)*Elliptic.E(arctan(x*c^-1//2*d^(1/2)),1+-1*b*c*a^-1*d^-1)+a^-1*c^(1/2)*d^-1//2*(c*a^-1*(c+d*x^2)^-1*(a+b*x^2))^-1//2*(a+b*x^2)^(1/2)*(c+d*x^2)^-1//2*(b*c+-1*a*d)^-1*(b*e+-1*a*f)*Elliptic.F(arctan(x*c^-1//2*d^(1/2)),1+-1*b*c*a^-1*d^-1))
@test integrate((a+-1*b*x^2)^-1//2*(c+d*x^2)^-3//2*(e+f*x^2), x) == :(x*c^-1*(a+-1*b*x^2)^(1/2)*(c+d*x^2)^-1//2*(a*d+b*c)^-1*(d*e+-1*c*f)+f*a^(1/2)*b^-1//2*d^-1*(1+d*c^-1*x^2)^(1/2)*(1+-1*b*a^-1*x^2)^(1/2)*(a+-1*b*x^2)^-1//2*(c+d*x^2)^-1//2*Elliptic.F(arcsin(x*a^-1//2*b^(1/2)),-1*a*d*b^-1*c^-1)+a^(1/2)*b^(1/2)*c^-1*d^-1*(1+d*c^-1*x^2)^-1//2*(1+-1*b*a^-1*x^2)^(1/2)*(a+-1*b*x^2)^-1//2*(c+d*x^2)^(1/2)*(a*d+b*c)^-1*(d*e+-1*c*f)*Elliptic.E(arcsin(x*a^-1//2*b^(1/2)),-1*a*d*b^-1*c^-1))
@test integrate((a+b*x^2)^-1//2*(c+-1*d*x^2)^-3//2*(e+f*x^2), x) == :(x*c^-1*(a+b*x^2)^(1/2)*(c+-1*d*x^2)^-1//2*(a*d+b*c)^-1*(c*f+d*e)+e*c^-1//2*d^-1//2*(1+b*a^-1*x^2)^(1/2)*(1+-1*d*c^-1*x^2)^(1/2)*(a+b*x^2)^-1//2*(c+-1*d*x^2)^-1//2*Elliptic.F(arcsin(x*c^-1//2*d^(1/2)),-1*b*c*a^-1*d^-1)+-1*c^-1//2*d^-1//2*(1+b*a^-1*x^2)^-1//2*(1+-1*d*c^-1*x^2)^(1/2)*(a+b*x^2)^(1/2)*(c+-1*d*x^2)^-1//2*(a*d+b*c)^-1*(c*f+d*e)*Elliptic.E(arcsin(x*c^-1//2*d^(1/2)),-1*b*c*a^-1*d^-1))
@test integrate((a+-1*b*x^2)^-1//2*(c+-1*d*x^2)^-3//2*(e+f*x^2), x) == :(-1*x*c^-1*(a+-1*b*x^2)^(1/2)*(c+-1*d*x^2)^-1//2*(b*c+-1*a*d)^-1*(c*f+d*e)+e*c^-1//2*d^-1//2*(1+-1*b*a^-1*x^2)^(1/2)*(1+-1*d*c^-1*x^2)^(1/2)*(a+-1*b*x^2)^-1//2*(c+-1*d*x^2)^-1//2*Elliptic.F(arcsin(x*c^-1//2*d^(1/2)),b*c*a^-1*d^-1)+c^-1//2*d^-1//2*(1+-1*b*a^-1*x^2)^-1//2*(1+-1*d*c^-1*x^2)^(1/2)*(a+-1*b*x^2)^(1/2)*(c+-1*d*x^2)^-1//2*(b*c+-1*a*d)^-1*(c*f+d*e)*Elliptic.E(arcsin(x*c^-1//2*d^(1/2)),b*c*a^-1*d^-1))
@test integrate((2+d*x^2)^-1//2*(3+f*x^2)^-1//2*(a+b*x^2), x) == :(b*x*d^-1*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2+(1/2)*a*2^(1/2)*f^-1//2*((3+f*x^2)^-1*(2+d*x^2))^-1//2*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*Elliptic.F(arctan(1//3*x*3^(1/2)*f^(1/2)),1+-3//2*d*f^-1)+-1*b*2^(1/2)*d^-1*f^-1//2*((3+f*x^2)^-1*(2+d*x^2))^-1//2*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*Elliptic.E(arctan(1//3*x*3^(1/2)*f^(1/2)),1+-3//2*d*f^-1))
@test integrate((2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*(a+b*x^2), x) == :(1//3*b*x*f^-1*(2+d*x^2)^(1/2)*(3+f*x^2)^(1/2)+-1//3*x*d^-1*f^-1*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*(-2*b*f+6*b*d+-3*a*d*f)+-1*2^(1/2)*f^-3//2*((3+f*x^2)^-1*(2+d*x^2))^-1//2*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*(b+-1*a*f)*Elliptic.F(arctan(1//3*x*3^(1/2)*f^(1/2)),1+-3//2*d*f^-1)+1//3*2^(1/2)*d^-1*f^-3//2*((3+f*x^2)^-1*(2+d*x^2))^-1//2*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*(-2*b*f+6*b*d+-3*a*d*f)*Elliptic.E(arctan(1//3*x*3^(1/2)*f^(1/2)),1+-3//2*d*f^-1))
@test integrate((2+d*x^2)^(1/2)*(3+f*x^2)^(1/2)*(a+b*x^2), x) == :(1//5*b*x*d^-1*(2+d*x^2)^3//2*(3+f*x^2)^(1/2)+1//15*x*d^-1*f^-1*(2+d*x^2)^(1/2)*(3+f*x^2)^(1/2)*(-4*b*f+3*b*d+5*a*d*f)+1//15*x*d^-2*f^-1*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*(-2*b*(4*f^2+9*d^2+-6*d*f)+5*a*d*f*(2*f+3d))+-1//5*2^(1/2)*d^-1*f^-3//2*((3+f*x^2)^-1*(2+d*x^2))^-1//2*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*(2*b*f+3*b*d+-10*a*d*f)*Elliptic.F(arctan(1//3*x*3^(1/2)*f^(1/2)),1+-3//2*d*f^-1)+-1//15*2^(1/2)*d^-2*f^-3//2*((3+f*x^2)^-1*(2+d*x^2))^-1//2*(2+d*x^2)^(1/2)*(3+f*x^2)^-1//2*(-2*b*(4*f^2+9*d^2+-6*d*f)+5*a*d*f*(2*f+3d))*Elliptic.E(arctan(1//3*x*3^(1/2)*f^(1/2)),1+-3//2*d*f^-1))
@test integrate((1+2*c*x^2*((b^2+-4*a*c)^(1/2)+-1b)^-1)^-1//2*(1+2*c*x^2*(-1b+-1*(b^2+-4*a*c)^(1/2))^-1)^-1//2*(-1b+-1*(b^2+-4*a*c)^(1/2)+2*c*x^2), x) == :(-1//2*2^(1/2)*c^-1//2*(b+-1*(b^2+-4*a*c)^(1/2))^(1/2)*(b+(b^2+-4*a*c)^(1/2))*Elliptic.E(arcsin(x*2^(1/2)*c^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))^-1//2),(b+(b^2+-4*a*c)^(1/2))^-1*(b+-1*(b^2+-4*a*c)^(1/2))))
@test integrate((1+2*c*x^2*(b+(b^2+-4*a*c)^(1/2))^-1)^-1//2*(1+2*c*x^2*(b+-1*(b^2+-4*a*c)^(1/2))^-1)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2*c*x^2), x) == :(x*(1+2*c*x^2*(b+(b^2+-4*a*c)^(1/2))^-1)^-1//2*(1+2*c*x^2*(b+-1*(b^2+-4*a*c)^(1/2))^-1)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))+(1/2)*2^(1/2)*c^-1//2*((1+2*c*x^2*(b+(b^2+-4*a*c)^(1/2))^-1)^-1*(1+2*c*x^2*(b+-1*(b^2+-4*a*c)^(1/2))^-1))^-1//2*(1+2*c*x^2*(b+(b^2+-4*a*c)^(1/2))^-1)^-1//2*(1+2*c*x^2*(b+-1*(b^2+-4*a*c)^(1/2))^-1)^(1/2)*(b+(b^2+-4*a*c)^(1/2))^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))*Elliptic.F(arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2),-2*(b+-1*(b^2+-4*a*c)^(1/2))^-1*(b^2+-4*a*c)^(1/2))+-1//2*2^(1/2)*c^-1//2*((1+2*c*x^2*(b+(b^2+-4*a*c)^(1/2))^-1)^-1*(1+2*c*x^2*(b+-1*(b^2+-4*a*c)^(1/2))^-1))^-1//2*(1+2*c*x^2*(b+(b^2+-4*a*c)^(1/2))^-1)^-1//2*(1+2*c*x^2*(b+-1*(b^2+-4*a*c)^(1/2))^-1)^(1/2)*(b+(b^2+-4*a*c)^(1/2))^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2))*Elliptic.E(arctan(x*2^(1/2)*c^(1/2)*(b+(b^2+-4*a*c)^(1/2))^-1//2),-2*(b+-1*(b^2+-4*a*c)^(1/2))^-1*(b^2+-4*a*c)^(1/2)))
@test integrate((c+d*x^2)^(1/2)*(e+f*x^2)^-1*(a+b*x^2), x) == :((1/2)*b*x*f^-1*(c+d*x^2)^(1/2)+-1//2*d^-1//2*f^-2*(-1*b*c*f+-2*a*d*f+2*b*d*e)*arctanh(x*d^(1/2)*(c+d*x^2)^-1//2)+e^-1//2*f^-2*(d*e+-1*c*f)^(1/2)*(b*e+-1*a*f)*arctanh(x*e^-1//2*(c+d*x^2)^-1//2*(d*e+-1*c*f)^(1/2)))
@test integrate((a+b*x^2)^3*(c+d*x^2)^-1*(e+f*x^2)^-1//2, x) == :(b*d^-3*f^-1//2*(b*c+-1*a*d)^2*arctanh(x*f^(1/2)*(e+f*x^2)^-1//2)+-1*c^-1//2*d^-3*(b*c+-1*a*d)^3*(d*e+-1*c*f)^-1//2*arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2))+1//8*b*d^-1*f^-5//2*(3*b^2*e^2+8*a^2*f^2+-8*a*b*e*f)*arctanh(x*f^(1/2)*(e+f*x^2)^-1//2)+(1/2)*b*d^-2*f^-3//2*(b*c+-1*a*d)*(b*e+-2*a*f)*arctanh(x*f^(1/2)*(e+f*x^2)^-1//2)+-3//8*x*b^2*d^-1*f^-2*(e+f*x^2)^(1/2)*(b*e+-2*a*f)+-1//2*x*b^2*d^-2*f^-1*(e+f*x^2)^(1/2)*(b*c+-1*a*d)+1//4*x*b^2*d^-1*f^-1*(e+f*x^2)^(1/2)*(a+b*x^2))
@test integrate((a+b*x^2)^2*(c+d*x^2)^-1*(e+f*x^2)^-1//2, x) == :(c^-1//2*d^-2*(b*c+-1*a*d)^2*(d*e+-1*c*f)^-1//2*arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2))+(1/2)*x*b^2*d^-1*f^-1*(e+f*x^2)^(1/2)+-1*b*d^-2*f^-1//2*(b*c+-1*a*d)*arctanh(x*f^(1/2)*(e+f*x^2)^-1//2)+-1//2*b*d^-1*f^-3//2*(b*e+-2*a*f)*arctanh(x*f^(1/2)*(e+f*x^2)^-1//2))
@test integrate((c+d*x^2)^-1*(e+f*x^2)^-1//2*(a+b*x^2), x) == :(b*d^-1*f^-1//2*arctanh(x*f^(1/2)*(e+f*x^2)^-1//2)+-1*c^-1//2*d^-1*(d*e+-1*c*f)^-1//2*(b*c+-1*a*d)*arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2)))
@test integrate((c+d*x^2)^-1*(e+f*x^2)^-1//2, x) == :(c^-1//2*(d*e+-1*c*f)^-1//2*arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2)))
@test integrate((a+b*x^2)^-1*(c+d*x^2)^-1*(e+f*x^2)^-1//2, x) == :(b*a^-1//2*(b*c+-1*a*d)^-1*(b*e+-1*a*f)^-1//2*arctan(x*a^-1//2*(e+f*x^2)^-1//2*(b*e+-1*a*f)^(1/2))+-1*d*c^-1//2*(b*c+-1*a*d)^-1*(d*e+-1*c*f)^-1//2*arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2)))
@test integrate((a+b*x^2)^-2*(c+d*x^2)^-1*(e+f*x^2)^-1//2, x) == :(c^-1//2*d^2*(b*c+-1*a*d)^-2*(d*e+-1*c*f)^-1//2*arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2))+(1/2)*b*a^-3//2*(b*c+-1*a*d)^-2*(b*e+-1*a*f)^-3//2*(c*e*b^2+4*d*f*a^2+-3*a*b*d*e+-2*a*b*c*f)*arctan(x*a^-1//2*(e+f*x^2)^-1//2*(b*e+-1*a*f)^(1/2))+(1/2)*x*a^-1*b^2*(a+b*x^2)^-1*(e+f*x^2)^(1/2)*(b*c+-1*a*d)^-1*(b*e+-1*a*f)^-1)
@test integrate((1+-1*x^2)^(1/2)*(-1+x^2)^-1*(a+b*x^2)^-1//2, x) == :(-1*(1+b*a^-1*x^2)^(1/2)*(a+b*x^2)^-1//2*Elliptic.F(arcsin(x),-1*b*a^-1))
@test integrate((c+d*x^2)^-1//2*(e+f*x^2)^-2*(a+b*x^2), x) == :(-1//2*e^-3//2*(d*e+-1*c*f)^-3//2*(a*c*f+b*c*e+-2*a*d*e)*arctanh(x*e^-1//2*(c+d*x^2)^-1//2*(d*e+-1*c*f)^(1/2))+(1/2)*x*e^-1*(c+d*x^2)^(1/2)*(e+f*x^2)^-1*(d*e+-1*c*f)^-1*(b*e+-1*a*f))
@test integrate((a+b*x^2)^-3//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2, x) == :(a^-1*e^(1/2)*(a*c^-1*(a+b*x^2)^-1*(c+d*x^2))^-1//2*(a*e^-1*(a+b*x^2)^-1*(e+f*x^2))^(1/2)*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*e+-1*a*f)^-1//2*Elliptic.E(arcsin(x*e^-1//2*(a+b*x^2)^-1//2*(b*e+-1*a*f)^(1/2)),e*c^-1*(b*e+-1*a*f)^-1*(b*c+-1*a*d)))
@test integrate((a+b*x^2)^-1//2*(c+d*x^2)^(1/2)*(e+f*x^2)^-3//2, x) == :(x*e^-1*(a+b*x^2)^(1/2)*(c+d*x^2)^-1//2*(e+f*x^2)^-1//2*(b*e+-1*a*f)^-1*(d*e+-1*c*f)+a^-1*c^3//2*e^-1*(c*a^-1*(c+d*x^2)^-1*(a+b*x^2))^-1//2*(a+b*x^2)^(1/2)*(c+d*x^2)^-1//2*(d*e+-1*c*f)^-1//2*Elliptic.F(arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2)),-1*e*a^-1*(d*e+-1*c*f)^-1*(b*c+-1*a*d))+-1*c^(1/2)*e^-1*(c*a^-1*(c+d*x^2)^-1*(a+b*x^2))^-1//2*(a+b*x^2)^(1/2)*(c+d*x^2)^-1//2*(b*e+-1*a*f)^-1*(d*e+-1*c*f)^(1/2)*Elliptic.E(arctan(x*c^-1//2*(e+f*x^2)^-1//2*(d*e+-1*c*f)^(1/2)),-1*e*a^-1*(d*e+-1*c*f)^-1*(b*c+-1*a*d)))
@test integrate((a+b*x^2)^-1//2*(c+d*x^2)^-1//2*(e+f*x^2)^-1//2, x) == :(c^-1*e^(1/2)*(a*c^-1*(a+b*x^2)^-1*(c+d*x^2))^-1//2*(a*e^-1*(a+b*x^2)^-1*(e+f*x^2))^(1/2)*(c+d*x^2)^(1/2)*(e+f*x^2)^-1//2*(b*e+-1*a*f)^-1//2*Elliptic.F(arcsin(x*e^-1//2*(a+b*x^2)^-1//2*(b*e+-1*a*f)^(1/2)),e*c^-1*(b*e+-1*a*f)^-1*(b*c+-1*a*d)))

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.5 (a+b x^2)^p (c+d x^2)^q (e+f x^2)^r.jl