text:
# initialization

# CUDA packages require complex initialization (discover CUDA, download artifacts, etc)
# that can't happen at module load time, so defer that to run time upon actual use.

const configured = Threads.Atomic{Int}(-1)   # -1=unconfigured, -2=configuring,
                                             # 0=failed, 1=configured

"""
    functional(show_reason=false)

Check if the package has been configured successfully and is ready to use.

This call is intended for packages that support conditionally using an available GPU. If you
fail to check whether CUDA is functional, actual use of functionality might warn and error.
"""
function functional(show_reason::Bool=false)
    if configured[] < 0
        Bool(_functional(show_reason))
    else
        Bool(configured[])
    end
end

const configure_lock = ReentrantLock()
@noinline function _functional(show_reason::Bool=false)
    Base.@lock configure_lock begin
        if configured[] == -1
            configured[] = -2
            configured[] = try
                __runtime_init__()
                1
            catch ex
                show_reason && @error("Error during initialization of CUDA.jl", exception=(ex,catch_backtrace()))
                0
            end
        elseif configured[] == -2
            # claim that we're initialized _during_ initialization (this can only happen
            # on the same task that's currently initializing, due to the reentrant lock)
            1
        else
            # we got stuck on the lock while another thread was initializing
            configured[]
        end
    end
end

# macro to guard code that only can run after the package has successfully initialized
macro after_init(ex)
    quote
        if !functional(true)
            error("""CUDA.jl did not successfully initialize, and is not usable.
                     If you did not see any other error message, try again in a new session
                     with the JULIA_DEBUG environment variable set to 'CUDA'.""")
        end
        $(esc(ex))
    end
end

## deferred initialization API

const __libcuda = Sys.iswindows() ? "nvcuda" : ( Sys.islinux() ? "libcuda.so.1" : "libcuda" )
libcuda() = @after_init(__libcuda)

# load-time initialization: only perform mininal checks here
function __init__()
    if Base.libllvm_version != LLVM.version()
        error("LLVM $(LLVM.version()) incompatible with Julia's LLVM $(Base.libllvm_version)")
    end

    # enable generation of FMA instructions to mimic behavior of nvcc
    LLVM.clopts("-nvptx-fma-level=1")

    precompiling = ccall(:jl_generating_output, Cint, ()) != 0
    if !precompiling
        eval(overrides)
    end
end

function __runtime_init__()
    if haskey(ENV, "_") && basename(ENV["_"]) == "rr"
        error("Running under rr, which is incompatible with CUDA")
    end

    @debug "Initializing CUDA driver"
    res = ccall((:cuInit, __libcuda), CUresult, (UInt32,), 0)
    if res == 0xffffffff
        error("Cannot use the CUDA stub libraries. You either don't have the NVIDIA driver installed, or it is not properly discoverable.")
    elseif res != SUCCESS
        throw_api_error(res)
    end

    if version() < v"10.1"
        @warn "This version of CUDA.jl only supports NVIDIA drivers for CUDA 10.1 or higher (yours is for CUDA $(version()))"
    end

    __init_dependencies__() || error("Could not find a suitable CUDA installation")

    if toolkit_release() < v"10.1"
        @warn "This version of CUDA.jl only supports CUDA 10.1 or higher (your toolkit provides CUDA $(toolkit_release()))"
    elseif toolkit_release() > release()
        @warn """You are using CUDA toolkit $(toolkit_release()) with a driver that only supports up to $(release()).
                 It is recommended to upgrade your driver, or switch to automatic installation of CUDA."""
    end

    if has_cudnn()
        cudnn_release = VersionNumber(CUDNN.version().major, CUDNN.version().minor)
        if cudnn_release < v"8.0"
            @warn "This version of CUDA.jl only supports CUDNN 8.0 or higher"
        end
    end

    if has_cutensor()
        cutensor_release = VersionNumber(CUTENSOR.version().major, CUTENSOR.version().minor)
        if !(v"1.0" <= cutensor_release <= v"1.2")
            @warn "This version of CUDA.jl only supports CUTENSOR 1.0 to 1.2"
        end
    end

    resize!(__device_contexts, ndevices())
    fill!(__device_contexts, nothing)

    __init_compatibility__()

    __init_pool__()

    CUBLAS.__runtime_init__()
    has_cudnn() && CUDNN.__runtime_init__()
    return
end

## convenience functions

# TODO: update docstrings

export has_cuda, has_cuda_gpu

"""
    has_cuda()::Bool

Check whether the local system provides an installation of the CUDA driver and toolkit.
Use this function if your code loads packages that require CUDA.jl.

Note that CUDA-dependent packages might still fail to load if the installation is broken,
so it's recommended to guard against that and print a warning to inform the user:

```
using CUDA
if has_cuda()
    try
        using CuArrays
    catch ex
        @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex,catch_backtrace())
    end
end
```
"""
has_cuda(show_reason::Bool=false) = functional(show_reason)

"""
    has_cuda_gpu()::Bool

Check whether the local system provides an installation of the CUDA driver and toolkit, and
if it contains a CUDA-capable GPU. See [`has_cuda`](@ref) for more details.

Note that this function initializes the CUDA API in order to check for the number of GPUs.
"""
has_cuda_gpu(show_reason::Bool=false) = has_cuda(show_reason) && length(devices()) > 0


text:
module GridVisualize

using Printf
using LinearAlgebra

using DocStringExtensions
using OrderedCollections
using ElasticArrays
using StaticArrays
using Colors
using ColorSchemes
using GeometryBasics
using PkgVersion

using ExtendableGrids

include("dispatch.jl")
include("common.jl")
include("pyplot.jl")
include("makie.jl")
include("vtkview.jl")
include("meshcat.jl")
include("plots.jl")
include("plutovista.jl")

export scalarplot,scalarplot!
export gridplot,gridplot!
export save,reveal,backend!
export isplots,isvtkview,ispyplot,ismakie,isplutovista
export GridVisualizer, SubVisualizer
export plottertype, available_kwargs
export default_plotter!, default_plotter
export PyPlotType,MakieType,PlotsType,VTKViewType, PlutoVistaType

end

ext: jl


text:
"""
    MetropolisHastings

Sampling of the initial conditions using the Metropolis-Hastings Markov chain Monte Carlo method.

Included within is the ability to sample the canonical distribution for adiabatic classical and
ring polymer systems.

Usage involves creating an instance of an `AbstractSystem{MonteCarlo}` and calling
`run_monte_carlo_sampling`.
"""
module MetropolisHastings

using Random
using StatsBase: AnalyticWeights, sample
using ProgressMeter: @showprogress
using Distributions: Normal

using NQCBase: NQCBase, Atoms
using NQCDynamics:
    NQCDynamics,
    AbstractSimulation,
    Simulation,
    RingPolymerSimulation,
    RingPolymers,
    Estimators,
    DynamicsUtils,
    ndofs

export run_monte_carlo_sampling

abstract type MonteCarloParameters end

"""
Parameters for Monte carlo simulations.
"""
mutable struct MonteCarlo{T} <: MonteCarloParameters
    Eᵢ::T
    Eₚ::T
    Δ::Dict{Symbol, T}
    steps::UInt
    moveable_atoms::AnalyticWeights
    extra_function::Function
    function MonteCarlo{T}(Δ::Dict{Symbol, T},
        n_atoms::Real, passes::Real, fix::Vector{<:Integer}, extra::Function) where {T<:AbstractFloat}

        steps = n_atoms * passes

        moveables = ones(n_atoms)
        moveables[fix] .= 0
        moveable_atoms = AnalyticWeights(moveables)

        new(0.0, 0.0, Δ, steps, moveable_atoms, extra)
    end
end

mutable struct PathIntegralMonteCarlo{T} <: MonteCarloParameters
    Eᵢ::T
    Eₚ::T
    Δ::Dict{Symbol, T}
    steps::UInt
    moveable_atoms::AnalyticWeights
    segment_length::UInt
    pass_length::UInt
    extra_function::Function
    function PathIntegralMonteCarlo{T}(Δ::Dict{Symbol, T},
        n_atoms::Real, passes::Real, fix::Vector{<:Integer},
        segment::Real, n_beads::Integer, extra::Function) where {T<:AbstractFloat}

        steps = n_atoms * passes

        moveables = ones(n_atoms)
        moveables[fix] .= 0
        moveable_atoms = AnalyticWeights(moveables)

        segment_length = max(1, round(Int, segment * n_beads))
        pass_length = n_beads ÷ segment_length

        new(0.0, 0.0, Δ, steps, moveable_atoms, segment_length, pass_length, extra)
    end
end

"""Container for storing simulation quantities
"""
mutable struct MonteCarloOutput{T<:AbstractFloat,S}
    R::Vector{S}
    acceptance::Dict{Symbol, T}
    total_moves::Dict{Symbol, UInt}
    energy::Vector{T}
    function MonteCarloOutput(R0::S, atoms::Atoms{N,T}) where {N,T,S}
        output = typeof(R0)[]
        acceptance = Dict{Symbol, T}()
        total_moves = Dict{Symbol, UInt}()
        for element in atoms.types
            acceptance[element] = 0.0
            total_moves[element] = 0
        end
        energy = T[]
        new{T,S}(output, acceptance, total_moves, energy)
    end
end

"""
    run_monte_carlo_sampling(sim::AbstractSimulation, monte::MonteCarloParameters, R0)

Perform Monte Carlo sampling for the system defined by the `sim` and `monte` parameters.

From the initial positions specified `R0` the system will be explored using the
Metropolis-Hastings algorithm.
"""
function run_monte_carlo_sampling(sim::Simulation, R0::Matrix{T},
    Δ::Dict{Symbol,T}, passes::Real; fix::Vector{<:Integer}=Int[], extra_function::Function=x->true) where {T}

    Rᵢ = copy(R0) # Current positions
    Rₚ = zero(Rᵢ) # Proposed positions
    monte = MonteCarlo{T}(Δ, length(sim.atoms), passes, fix, extra_function)
    monte.Eᵢ = DynamicsUtils.classical_potential_energy(sim, Rᵢ)
    output = MonteCarloOutput(Rᵢ, sim.atoms)

    run_main_loop!(sim, monte, Rᵢ, Rₚ, output)

    for key in keys(output.acceptance)
        output.acceptance[key] /= output.total_moves[key]
    end
    output
end

function run_monte_carlo_sampling(sim::RingPolymerSimulation, R0::AbstractArray{T,3},
    Δ::Dict{Symbol,T}, passes::Real;
    fix::Vector{<:Integer}=Int[], segment::Real=1, extra_function::Function=x->true) where {T}

    Rᵢ = copy(R0) # Current positions
    Rₚ = zero(Rᵢ) # Proposed positions
    monte = PathIntegralMonteCarlo{T}(Δ, length(sim.atoms), passes, fix, segment, length(sim.beads), extra_function)
    monte.Eᵢ = DynamicsUtils.classical_potential_energy(sim, Rᵢ)
    output = MonteCarloOutput(Rᵢ, sim.atoms)

    run_main_loop!(sim, monte, Rᵢ, Rₚ, output)

    for key in keys(output.acceptance)
        output.acceptance[key] /= output.total_moves[key]
    end
    output
end

"""
    run_main_loop!(system::System{MonteCarlo}, Rᵢ::Matrix, Rₚ::Matrix, output::MonteCarloOutput)

Main loop for classical systems.
"""
function run_main_loop!(sim::Simulation, monte::MonteCarlo, Rᵢ::Matrix, Rₚ::Matrix, output::MonteCarloOutput)
    @showprogress 0.1 "Sampling... " for i=1:convert(Int, monte.steps)
        atom = sample(range(sim.atoms), monte.moveable_atoms)
        propose_move!(sim, monte, Rᵢ, Rₚ, atom)
        @views for i in axes(Rₚ, 2) # atoms
            if monte.moveable_atoms[i] > 0
                NQCBase.apply_cell_boundaries!(sim.cell, Rₚ[:,i])
            end
        end
        assess_proposal!(sim, monte, Rᵢ, Rₚ, output, atom)
    end
end

"""Main loop for ring polymer systems.
"""
function run_main_loop!(sim::RingPolymerSimulation, monte::PathIntegralMonteCarlo, Rᵢ::Array{T,3}, Rₚ::Array{T,3},
    output::MonteCarloOutput) where {T}
    @showprogress 0.1 "Sampling... " for i=1:convert(Int, monte.steps)
        if !all(monte.moveable_atoms .== 0.0)
            atom = sample(range(sim.atoms), monte.moveable_atoms)
            propose_centroid_move!(sim, monte, Rᵢ, Rₚ, atom)
            NQCBase.apply_cell_boundaries!(sim.cell, Rₚ, sim.beads)
            assess_proposal!(sim, monte, Rᵢ, Rₚ, output, atom)
            if length(sim.beads) > 1
                if atom ∈ sim.beads.quantum_atoms
                    for j=1:monte.pass_length
                        propose_normal_mode_move!(sim, monte, Rᵢ, Rₚ, atom)
                        assess_proposal!(sim, monte, Rᵢ, Rₚ, output, atom)
                    end
                end
            end
        else
            if length(sim.beads) > 1
                atom = sample(sim.beads.quantum_atoms)
                for j=1:monte.pass_length
                    propose_normal_mode_move!(sim, monte, Rᵢ, Rₚ, atom)
                    assess_proposal!(sim, monte, Rᵢ, Rₚ, output, atom)
                end
            end
        end
        assess_proposal!(sim, monte, Rᵢ, Rₚ, output, atom)
    end
end

"""
    propose_move!(system::System{MonteCarlo}, Rᵢ::Matrix{T}, Rₚ::Matrix{T}) where {T<:AbstractFloat}

Propose simple cartesian move for a single atom.
"""
function propose_move!(sim::Simulation, monte::MonteCarlo, Rᵢ::Matrix, Rₚ::Matrix, atom::Integer)
    Rₚ .= Rᵢ
    apply_random_perturbation!(sim.atoms, monte, Rₚ, atom, ndofs(sim))
end

"""
    propose_centroid_move!(system::RingPolymerSystem{MonteCarlo}, Rᵢ::Array{T, 3}, Rₚ::Array{T, 3}) where {T<:AbstractFloat}

Propose a move for the ring polymer centroid for one atom.
"""
function propose_centroid_move!(sim::RingPolymerSimulation, monte::PathIntegralMonteCarlo, Rᵢ::Array{T, 3}, Rₚ::Array{T, 3}, atom::Integer) where {T<:AbstractFloat}
    Rₚ .= Rᵢ
    RingPolymers.transform_to_normal_modes!(sim.beads, Rₚ)
    @views apply_random_perturbation!(sim.atoms, monte, Rₚ[:,:,1], atom, ndofs(sim)) # Perturb the centroid mode
    @views if atom ∉ sim.beads.quantum_atoms # Ensure replicas are identical if not quantum
        for i=2:length(sim.beads)
            Rₚ[:, atom, i] .= Rₚ[:, atom, 1]
        end
    end
    RingPolymers.transform_from_normal_modes!(sim.beads, Rₚ)
end

"""
    apply_random_perturbation!(system::AbstractSystem{MonteCarlo}, R::AbstractMatrix, atom::Integer)

Randomly perturb the xyz coordinates of a single atom.
"""
function apply_random_perturbation!(atoms::Atoms, monte::MonteCarloParameters, R::AbstractMatrix, atom::Integer, DoFs::Integer)
    R[:, atom] .+= (rand(DoFs) .- 0.5) .* monte.Δ[atoms.types[atom]] / sqrt(DoFs)
end

"""
    propose_normal_mode_move!(system::RingPolymerSystem{MonteCarlo}, Rᵢ::Array{T, 3}, Rₚ::Array{T, 3}) where {T}

Propose a move for a single normal mode for a single atom.
"""
function propose_normal_mode_move!(sim::RingPolymerSimulation, monte::PathIntegralMonteCarlo, Rᵢ::Array{T,3},  Rₚ::Array{T, 3}, i::Integer) where {T}
    Rₚ .= Rᵢ
    RingPolymers.transform_to_normal_modes!(sim.beads, Rₚ)
    for _=1:monte.segment_length
        mode = sample(2:length(sim.beads))
        @views sample_mode!(sim.beads, sim.atoms.masses[i], mode, Rₚ[:, i, mode])
    end
    RingPolymers.transform_from_normal_modes!(sim.beads, Rₚ)
end

function sample_mode!(beads::RingPolymers.RingPolymerParameters, mass::AbstractFloat, mode::Integer, R::AbstractArray)
    σ = sqrt(beads.ω_n / mass) / (beads.normal_mode_springs[mode] * mass)
    rand!(Normal(0, σ), R)
end

"""
    assess_proposal!(system::AbstractSystem{MonteCarlo}, Rᵢ, Rₚ, output, i)

Update the energy, check for acceptance, and update the output.
"""
function assess_proposal!(sim::AbstractSimulation, monte::MonteCarloParameters, Rᵢ, Rₚ, output, atom::Integer)
    monte.Eₚ = DynamicsUtils.classical_potential_energy(sim, Rₚ)
    output.total_moves[sim.atoms.types[atom]] += 1
    if monte.extra_function(Rₚ)
        if acceptance_probability(sim, monte) > rand()
            monte.Eᵢ = monte.Eₚ
            Rᵢ .= Rₚ
            output.acceptance[sim.atoms.types[atom]] += 1
        end
    end
    write_output!(output, Rᵢ, monte.Eᵢ)
end

"""
    acceptance_probability(system::AbstractSystem{MonteCarlo})

Return the Metropolis-Hastings acceptance probability.
"""
function acceptance_probability(sim::Simulation, monte::MonteCarlo)
    acceptance_probability(monte.Eₚ, monte.Eᵢ, NQCDynamics.get_temperature(sim))
end
function acceptance_probability(sim::RingPolymerSimulation, monte::PathIntegralMonteCarlo)
    acceptance_probability(monte.Eₚ, monte.Eᵢ, sim.beads.ω_n)
end
acceptance_probability(Eₚ::T, Eᵢ::T, kT::T) where {T<:AbstractFloat} = min(1, exp(-(Eₚ - Eᵢ)/kT))

"""
    write_output!(output::MonteCarloOutput, Rᵢ::AbstractArray, energy::AbstractFloat, i::Integer)

Store the current configuration and associated energy.
"""
function write_output!(output::MonteCarloOutput, Rᵢ::AbstractArray, energy::AbstractFloat)
    push!(output.R, copy(Rᵢ))
    push!(output.energy, copy(energy))
end

end # module


text:
# --- PRIVATE METHODS --------------------------------------------------------------------------------------- #

# array based -
function _build_ternary_lattice_intrinsic_value_array(contractSet::Set{PSAbstractAsset}, ternaryPriceArray::Array{Float64,1})::PSResult

    # initialize -
    intrinsic_value_array = similar(ternaryPriceArray)

    # ok, so lets go through each price, compute the intrinsic value - and then store
    for (index,price) in enumerate(ternaryPriceArray)

        # ok, get the underlying price -
        underlying_price_value = ternaryPriceArray[index]

        # compute the intrinsic value -
        result = intrinsic_value(contractSet,underlying_price_value)
        if (isa(result.value,Exception) == true)
            return result
        end
        iv_value = result.value

        # add the iv to the array - index should work out
        intrinsic_value_array[index] = iv_value
    end

    # return -
    return PSResult{Array{Float64,1}}(intrinsic_value_array)
end

function _build_ternary_lattice_underlying_price_array(basePrice::Float64, volatility::Float64, timeToExercise::Int)::PSResult

    # compute up and down perturbations -
    numberOfLevels = timeToExercise + 1     # assumption: our time unit is 1 day *always* = is this legit?
    Δt = (1.0)
    U = exp(volatility * sqrt(2*Δt))
    D = 1.0 / U
    C = 1.0

    # compute price array -
    number_of_elements = Int64(((3^numberOfLevels) - 1)/2)
    priceArray = zeros(number_of_elements)
    priceArray[1] = basePrice

    # populate the ternary price tree -
    for index = 1:number_of_elements

        # get the basePrice -
        basePrice = priceArray[index]

        # compute the prices -
        down_price = basePrice*D
        up_price = basePrice*U
        center_price = basePrice*C

        # add the prices to the array -
        left_index = 3*(index - 1) + 2
        center_index = 3*(index - 1) + 3
        right_index = 3*(index - 1) + 4

        # note - we need to check that we don't write into the array past the end
        if (left_index<=number_of_elements)
            priceArray[left_index] = up_price
        end

        if (center_index<=number_of_elements)
            priceArray[center_index] = center_price
        end

        if (right_index<=number_of_elements)
            priceArray[right_index] = down_price
        end
    end

    # return the price array -
    return PSResult{Array{Float64,1}}(priceArray)
end

function _build_ternary_lattice_option_value_array(intrinsicValueArray::Array{Float64,1}, latticeModel::PSTernaryLatticeModel; earlyExcercise::Bool = false)::PSResult

    # initialize -
    contract_price_array = copy(intrinsicValueArray)

    # get stuff from the lattice model -
    σ = latticeModel.volatility
    timeToExercise = latticeModel.timeToExercise
    riskFreeRate = latticeModel.riskFreeRate
    dividendRate = latticeModel.dividendRate

    # compute the probability -
    numberOfLevels = timeToExercise + 1     # assumption: our time unit is 1 day *always* = is this legit?
    Δt = (1.0)
    T1 = exp((riskFreeRate - dividendRate)*(Δt/2))
    T2 = exp(-σ*sqrt(Δt/2))
    T3 = exp(σ*sqrt(Δt/2))
    pup = ((T1 - T2)/(T3 - T2))^2
    pdown = ((T3 - T1)/(T3 - T2))^2
    DF = exp(-riskFreeRate*Δt)
    U = exp(σ*sqrt(2*Δt))
    D = 1.0 / U
    C = 1.0

    # create a index table -
    # number_of_elements = Int64(((3^numberOfLevels) - 1)/2)
    number_of_elements = numberOfLevels
    index_table = zeros(numberOfLevels, 5)
    backwards_index_array = range(numberOfLevels, step=-1,stop=1) |> collect
    for (forward_index, backward_index) in enumerate(backwards_index_array)

        # add the prices to the array -
        left_index = 3*(backward_index - 1) + 2
        center_index = 3*(backward_index - 1) + 3
        right_index = 3*(backward_index - 1) + 4

        # populate the index table -
        index_table[forward_index,1] = backward_index
        index_table[forward_index,2] = left_index
        index_table[forward_index,3] = center_index
        index_table[forward_index,4] = right_index
    end

    # @show index_table

    # ok, so now lets compute the value for the nodes -
    for compute_index = 1:number_of_elements

        # get the indexs -
        parent_node_index = Int(index_table[compute_index,1])
        child_left_index = Int(index_table[compute_index,2])
        child_center_index = Int(index_table[compute_index,3])
        child_right_index = Int(index_table[compute_index,4])

        #@show (parent_node_index,child_left_index,child_center_index,child_right_index)

        # compute -
        contract_price = DF*(pup*contract_price_array[child_left_index]+(pdown)*contract_price_array[child_right_index]+(1-pup-pdown)*contract_price_array[child_center_index])
        if (earlyExcercise == false)
            contract_price_array[parent_node_index] = contract_price
        else
            excercise_value = contract_price_array[parent_node_index]
            contract_price_array[parent_node_index] = max(excercise_value,contract_price)
        end
    end

   # setup the results tuple -
   results_tuple = (option_index_table=index_table, option_contract_price_array=contract_price_array, U=U, D=D, PUP=pup, PDOWN=pdown, PCENTER=(1-pup-pdown), DF=DF)

   # return -
   return PSResult(results_tuple)
end
# ----------------------------------------------------------------------------------------------------------- #

# --- PUBLIC METHODS ---------------------------------------------------------------------------------------- #

"""
    option_contract_price(contractSet::Set{PSAbstractAsset}, latticeModel::PSTernaryLatticeModel, baseUnderlyingPrice::Float64;
        earlyExercise::Bool = false)::PooksoftBase.PSResult

Estimate the price of a contract using a ternary lattice pricing model.
"""
function option_contract_price(contractSet::Set{PSAbstractAsset}, latticeModel::PSTernaryLatticeModel, baseUnderlyingPrice::Float64;
    earlyExercise::Bool = false)::PSResult

    # initialize -
    option_contract_price = 0.0

    # we need to get some stuff from the lattice model -
    volatility = latticeModel.volatility
    timeToExercise = latticeModel.timeToExercise

    # compute the price array -
    result = _build_ternary_lattice_underlying_price_array(baseUnderlyingPrice, volatility, timeToExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    lattice_price_array = result.value

    # compute the intrinsic value array -
    result = _build_ternary_lattice_intrinsic_value_array(contractSet, lattice_price_array)
    if (isa(result.value,Exception) == true)
        return result
    end
    iv_array = result.value

    # ok, let's build the option value array -
    result = _build_ternary_lattice_option_value_array(iv_array, latticeModel; earlyExcercise = earlyExercise)
    if (isa(result.value,Exception) == true)
        return result
    end
    cost_calc_tuple = result.value

    # create a named tuple and return the results -
    results_tuple = (lattice_price_array=lattice_price_array, intrinsic_value_array = iv_array, cost_calculation_result=cost_calc_tuple)

    # return -
    return PSResult(results_tuple)
end
# ----------------------------------------------------------------------------------------------------------- #


text:
#
# include("RigidBody.jl")
# include("Joint.jl")
# include("OrientationConversion.jl")

function extF(t::T, j::Vector{Joint})::Vector{extForces} where T<:Real
    # Function to generate external forces
    # that may or may not be functions of time
    extFList = zeroExtForceVec(length(j)+1)

    # For reference, calculate total mass of robot
    TotalMass = (0.0);
    for i=1:length(j)
        TotalMass += j[i].RB2.m
    end

    # Pendulum
    # extFList[2] = extForces([0.0 0 0.1],[0.0 0 0],[0.0 0.0 0.0]);

    # Double Pendulum Test
    # extFList[2] = zeroExtForce()
    # extFList[3] = extForces([0.0 0 0],[0.0 0 0],[0.0 0.0 0.0])

    # Cube Prop Test
    # extFList[2] = extForces([0.0 0.0 0.0], zeros(1,3), [0.0 0.0 0.0])
    # extFList[3] = extForces([0.0 0.0 0.0], zeros(1,3), [0.0 0.0 0.1])

    # QuadTest
    # extFList[3] = extForces([0.0 0.0 0.0], zeros(1,3), [0.0 0.0 0.1])#*sin(t)])
    # extFList[4] = extForces([0.0 0.0 0.0], zeros(1,3), [0.0 0.0 0.1])#*sin(t)])
    # extFList[5] = extForces([0.0 0.0 0.0], zeros(1,3), [0.0 0.0 -0.1*sin(t)])
    # extFList[6] = extForces([0.0 0.0 0.0], zeros(1,3), [0.0 0.0 -0.1*sin(t)])

    # nRotors = length(j) - 1;
    # C_T = 1.27e-6;
    # Ω_h = sqrt(TotalMass*(9.81)/nRotors/C_T);
    # flt = ones(nRotors);  # fault vector for the n Rotors
    # flt = [1.0;0.5;1.0;0.5];
    # for i=3:j[end].RB2.bodyID
    #     extFList[i] = genRotorF(Ω_h,flt[i-2]);
    # end

    # grav = [0.0,0.0,-9.806]

    return extFList
end

function genJointF(t::T,j::Joint)::Tuple{Vector{T}, Vector{T}} where T<:Real
    # For actuated joints.
    # j.jointF gives us the actuation force and torque applied ...
    # ... on the child link in the body frame of the parent link
    β1 = j.RB1.x[4:7]
    β2 = j.RB2.x[4:7]

    Fb1 = -j.jointF # Joint Force acting on parent link
    Fb2 = Vector{T}(undef,6) # Force acting on child link
    Fb2[1:3] = quat2dcm(β2)*transpose(quat2dcm(β1))*j.jointF[1:3]
    Fb2[4:6] = quat2dcm(β2)*transpose(quat2dcm(β1))*j.jointF[4:6]

    if j==2
        print("Fb = \t")
        println(Fb2)
    end
    return (Fb1,Fb2)

end

function genRotorF(Ω,f)
    # f: fault signal ∈ [0,1]. 0 indicates total fault and 1 indicates nominal performance
    C_T = 1.27e-6; # thrust coefficient
    T = C_T*Ω^2; # thrust
    F = f*[0;0;T];
    prop_Fτ = extForces(permutedims(F), zeros(1,3), zeros(1,3));
    return prop_Fτ;
end

## old stuff
# function genMotorF(t::Float64,j::Joint)::Tuple{Vector{Float64}, Vector{Float64}}
#     # For prop joints.
#     # input PWM value from the user (to be supplied from a different function in the future)
#     # lookupTable fn gives us the actuation force and torque applied
#     # ... on the child link in its body frame
#     β1 = j.RB1.x[4:7]
#     β2 = j.RB2.x[4:7]
#
#     FArr =     [1000      0      0;
#                 1100      0      0;
#                 1200 0.0261 0.1564;
#                 1300 0.0719 0.4443;
#                 1400 0.1147 0.7171;
#                 1500 0.1585 0.9877;
#                 1600 0.2138 1.3191]
#
#     pwm = 1400
#     # println("Enter PWM.")
#     # inp = parse(Int64,readline())
#     # if (inp==0) pwm = 1000 else pwm = inp end
#
#     Fb1 = zeros(6) # Joint Force acting on parent link
#     Fb2 = zeros(6) # Joint Force acting on child link
#
#     # Get Force and torque generated on the props in the prop frame
#     Fb2[3], Fb2[6] = lookupTable(FArr,pwm)
#     Fb1[1:3] = quat2dcm(β1)*transpose(quat2dcm(β2))*(-Fb2)[1:3]
#     Fb1[4:6] = quat2dcm(β1)*transpose(quat2dcm(β2))*(-Fb2)[4:6]
#
#     return (Fb1,Fb2)
# end
# ##
# # Testing lookup table implementation for PWM to F,τ
#
# # const τArr = [1000      0;
# #         1100      0;
# #         1200 0.1564;
# #         1300 0.4443;
# #         1400 0.7171;
# #         1500 0.9877;
# #         1600 1.3191]
#
# function lookupTable(FArr::Matrix{Float64},pwm::Real):: Tuple{Float64,Float64}
#
#     a = findfirst(isequal(pwm),FArr)
#     F = FArr[a[1],2]
#     τ = FArr[a[1],3]
#
#     return (F,τ)
# end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
@testset "1.1.2.8 P(x) (c x)^m (a+b x^2)^p" begin
    (A, B, C, D, F, a, b, c, d, e, f, m, x, ) = @variables A B C D F a b c d e f m x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(c*x)^m*Pq(x)*(a+b*x^2)^p=#

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

    #= ::Subsection:: =#
    #=Integrands*of*the*form*x^m*(A+B*x)*(a+b*x^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(A+B*x)*(a+b*x^2)^(p/2)=#

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

    @test_int [x^3*(A + B*x)*sqrt(a + b*x^2), x, 6, (a^2*B*x*sqrt(a + b*x^2))/(16*b^2) + (A*x^2*(a + b*x^2)^(3/2))/(5*b) + (B*x^3*(a + b*x^2)^(3/2))/(6*b) - (a*(16*A + 15*B*x)*(a + b*x^2)^(3/2))/(120*b^2) + (a^3*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(16*b^(5/2))]
    @test_int [x^2*(A + B*x)*sqrt(a + b*x^2), x, 5, -(a*A*x*sqrt(a + b*x^2))/(8*b) + (B*x^2*(a + b*x^2)^(3/2))/(5*b) - ((8*a*B - 15*A*b*x)*(a + b*x^2)^(3/2))/(60*b^2) - (a^2*A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(8*b^(3/2))]
    @test_int [x^1*(A + B*x)*sqrt(a + b*x^2), x, 4, -(a*B*x*sqrt(a + b*x^2))/(8*b) + ((4*A + 3*B*x)*(a + b*x^2)^(3/2))/(12*b) - (a^2*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(8*b^(3/2))]
    @test_int [x^0*(A + B*x)*sqrt(a + b*x^2), x, 4, (A*x*sqrt(a + b*x^2))/2 + (B*(a + b*x^2)^(3/2))/(3*b) + (a*A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(2*sqrt(b))]
    @test_int [(A + B*x)*sqrt(a + b*x^2)/x^1, x, 7, ((2*A + B*x)*sqrt(a + b*x^2))/2 + (a*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(2*sqrt(b)) - sqrt(a)*A*atanh(sqrt(a + b*x^2)/sqrt(a))]
    @test_int [(A + B*x)*sqrt(a + b*x^2)/x^2, x, 7, -(((A - B*x)*sqrt(a + b*x^2))/x) + A*sqrt(b)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)) - sqrt(a)*B*atanh(sqrt(a + b*x^2)/sqrt(a))]
    @test_int [(A + B*x)*sqrt(a + b*x^2)/x^3, x, 7, -((A + 2*B*x)*sqrt(a + b*x^2))/(2*x^2) + sqrt(b)*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)) - (A*b*atanh(sqrt(a + b*x^2)/sqrt(a)))/(2*sqrt(a))]

    @test_int [x^3*(A + B*x)*(a + b*x^2)^(3/2), x, 7, (3*a^3*B*x*sqrt(a + b*x^2))/(128*b^2) + (a^2*B*x*(a + b*x^2)^(3/2))/(64*b^2) + (A*x^2*(a + b*x^2)^(5/2))/(7*b) + (B*x^3*(a + b*x^2)^(5/2))/(8*b) - (a*(32*A + 35*B*x)*(a + b*x^2)^(5/2))/(560*b^2) + (3*a^4*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(128*b^(5/2))]
    @test_int [x^2*(A + B*x)*(a + b*x^2)^(3/2), x, 6, -(a^2*A*x*sqrt(a + b*x^2))/(16*b) - (a*A*x*(a + b*x^2)^(3/2))/(24*b) + (B*x^2*(a + b*x^2)^(5/2))/(7*b) - ((12*a*B - 35*A*b*x)*(a + b*x^2)^(5/2))/(210*b^2) - (a^3*A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(16*b^(3/2))]
    @test_int [x^1*(A + B*x)*(a + b*x^2)^(3/2), x, 5, -(a^2*B*x*sqrt(a + b*x^2))/(16*b) - (a*B*x*(a + b*x^2)^(3/2))/(24*b) + ((6*A + 5*B*x)*(a + b*x^2)^(5/2))/(30*b) - (a^3*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(16*b^(3/2))]
    @test_int [x^0*(A + B*x)*(a + b*x^2)^(3/2), x, 5, (3*a*A*x*sqrt(a + b*x^2))/8 + (A*x*(a + b*x^2)^(3/2))/4 + (B*(a + b*x^2)^(5/2))/(5*b) + (3*a^2*A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(8*sqrt(b))]
    @test_int [(A + B*x)*(a + b*x^2)^(3/2)/x^1, x, 8, (a*(8*A + 3*B*x)*sqrt(a + b*x^2))/8 + ((4*A + 3*B*x)*(a + b*x^2)^(3/2))/12 + (3*a^2*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(8*sqrt(b)) - a^(3/2)*A*atanh(sqrt(a + b*x^2)/sqrt(a))]
    @test_int [(A + B*x)*(a + b*x^2)^(3/2)/x^2, x, 8, ((2*a*B + 3*A*b*x)*sqrt(a + b*x^2))/2 - ((3*A - B*x)*(a + b*x^2)^(3/2))/(3*x) + (3*a*A*sqrt(b)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/2 - a^(3/2)*B*atanh(sqrt(a + b*x^2)/sqrt(a))]
    @test_int [(A + B*x)*(a + b*x^2)^(3/2)/x^3, x, 8, (-3*(a*B - A*b*x)*sqrt(a + b*x^2))/(2*x) - ((A - B*x)*(a + b*x^2)^(3/2))/(2*x^2) + (3*a*sqrt(b)*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/2 - (3*sqrt(a)*A*b*atanh(sqrt(a + b*x^2)/sqrt(a)))/2]

    @test_int [x^3*(A + B*x)*(a + b*x^2)^(5/2), x, 8, (3*a^4*B*x*sqrt(a + b*x^2))/(256*b^2) + (a^3*B*x*(a + b*x^2)^(3/2))/(128*b^2) + (a^2*B*x*(a + b*x^2)^(5/2))/(160*b^2) + (A*x^2*(a + b*x^2)^(7/2))/(9*b) + (B*x^3*(a + b*x^2)^(7/2))/(10*b) - (a*(160*A + 189*B*x)*(a + b*x^2)^(7/2))/(5040*b^2) + (3*a^5*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(256*b^(5/2))]
    @test_int [x^2*(A + B*x)*(a + b*x^2)^(5/2), x, 7, (-5*a^3*A*x*sqrt(a + b*x^2))/(128*b) - (5*a^2*A*x*(a + b*x^2)^(3/2))/(192*b) - (a*A*x*(a + b*x^2)^(5/2))/(48*b) + (B*x^2*(a + b*x^2)^(7/2))/(9*b) - ((16*a*B - 63*A*b*x)*(a + b*x^2)^(7/2))/(504*b^2) - (5*a^4*A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(128*b^(3/2))]
    @test_int [x^1*(A + B*x)*(a + b*x^2)^(5/2), x, 6, (-5*a^3*B*x*sqrt(a + b*x^2))/(128*b) - (5*a^2*B*x*(a + b*x^2)^(3/2))/(192*b) - (a*B*x*(a + b*x^2)^(5/2))/(48*b) + ((8*A + 7*B*x)*(a + b*x^2)^(7/2))/(56*b) - (5*a^4*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(128*b^(3/2))]
    @test_int [x^0*(A + B*x)*(a + b*x^2)^(5/2), x, 6, (5*a^2*A*x*sqrt(a + b*x^2))/16 + (5*a*A*x*(a + b*x^2)^(3/2))/24 + (A*x*(a + b*x^2)^(5/2))/6 + (B*(a + b*x^2)^(7/2))/(7*b) + (5*a^3*A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(16*sqrt(b))]
    @test_int [(A + B*x)*(a + b*x^2)^(5/2)/x^1, x, 9, (a^2*(16*A + 5*B*x)*sqrt(a + b*x^2))/16 + (a*(8*A + 5*B*x)*(a + b*x^2)^(3/2))/24 + ((6*A + 5*B*x)*(a + b*x^2)^(5/2))/30 + (5*a^3*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(16*sqrt(b)) - a^(5/2)*A*atanh(sqrt(a + b*x^2)/sqrt(a))]
    @test_int [(A + B*x)*(a + b*x^2)^(5/2)/x^2, x, 9, (a*(8*a*B + 15*A*b*x)*sqrt(a + b*x^2))/8 + ((4*a*B + 15*A*b*x)*(a + b*x^2)^(3/2))/12 - ((5*A - B*x)*(a + b*x^2)^(5/2))/(5*x) + (15*a^2*A*sqrt(b)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/8 - a^(5/2)*B*atanh(sqrt(a + b*x^2)/sqrt(a))]
    @test_int [(A + B*x)*(a + b*x^2)^(5/2)/x^3, x, 9, (5*a*b*(4*A + 3*B*x)*sqrt(a + b*x^2))/8 - (5*(3*a*B - 2*A*b*x)*(a + b*x^2)^(3/2))/(12*x) - ((2*A - B*x)*(a + b*x^2)^(5/2))/(4*x^2) + (15*a^2*sqrt(b)*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/8 - (5*a^(3/2)*A*b*atanh(sqrt(a + b*x^2)/sqrt(a)))/2]

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

    @test_int [x^3*(A + B*x)/sqrt(a + b*x^2), x, 5, (A*x^2*sqrt(a + b*x^2))/(3*b) + (B*x^3*sqrt(a + b*x^2))/(4*b) - (a*(16*A + 9*B*x)*sqrt(a + b*x^2))/(24*b^2) + (3*a^2*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(8*b^(5/2))]
    @test_int [x^2*(A + B*x)/sqrt(a + b*x^2), x, 4, (B*x^2*sqrt(a + b*x^2))/(3*b) - ((4*a*B - 3*A*b*x)*sqrt(a + b*x^2))/(6*b^2) - (a*A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(2*b^(3/2))]
    @test_int [x^1*(A + B*x)/sqrt(a + b*x^2), x, 3, ((2*A + B*x)*sqrt(a + b*x^2))/(2*b) - (a*B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(2*b^(3/2))]
    @test_int [x^0*(A + B*x)/sqrt(a + b*x^2), x, 3, (B*sqrt(a + b*x^2))/b + (A*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/sqrt(b)]
    @test_int [(A + B*x)/(x^1*sqrt(a + b*x^2)), x, 6, (B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/sqrt(b) - (A*atanh(sqrt(a + b*x^2)/sqrt(a)))/sqrt(a)]
    @test_int [(A + B*x)/(x^2*sqrt(a + b*x^2)), x, 4, -((A*sqrt(a + b*x^2))/(a*x)) - (B*atanh(sqrt(a + b*x^2)/sqrt(a)))/sqrt(a)]
    @test_int [(A + B*x)/(x^3*sqrt(a + b*x^2)), x, 5, -(A*sqrt(a + b*x^2))/(2*a*x^2) - (B*sqrt(a + b*x^2))/(a*x) + (A*b*atanh(sqrt(a + b*x^2)/sqrt(a)))/(2*a^(3/2))]

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

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(A+B*x)*(A^2-B^2*x^2)^(p/2)=#

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

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

    @test_int [x*(1 - x)/sqrt(1 - x^2), x, 2, (-(1/2))*(2 - x)*sqrt(1 - x^2) - asin(x)/2]
    @test_int [(x - x^2)/sqrt(1 - x^2), x, 3, (-(1/2))*(2 - x)*sqrt(1 - x^2) - asin(x)/2]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(A+B*x+C*x^2)*(a+b*x^2)^p=#

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(A+B*x+C*x^2)*(a+b*x^2)^(p/2)=#

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

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

    @test_int [x^7*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 7, -((x^7*(a*B - (A*b - a*C)*x))/(7*a*b*(a + b*x^2)^(7/2))) - (x^5*(7*a*B - (A*b - 8*a*C)*x))/(35*a*b^2*(a + b*x^2)^(5/2)) - (x^3*(35*a*B - 6*(A*b - 8*a*C)*x))/(105*a*b^3*(a + b*x^2)^(3/2)) - (x*(35*a*B - 8*(A*b - 8*a*C)*x))/(35*a*b^4*sqrt(a + b*x^2)) - (16*(A*b - 8*a*C)*sqrt(a + b*x^2))/(35*a*b^5) + (B*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/b^(9/2)]
    @test_int [x^6*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 6, -((x^6*(a*B - (A*b - a*C)*x))/(7*a*b*(a + b*x^2)^(7/2))) - (x^4*(6*B + 7*C*x))/(35*b^2*(a + b*x^2)^(5/2)) - (x^2*(24*B + 35*C*x))/(105*b^3*(a + b*x^2)^(3/2)) - (16*B + 35*C*x)/(35*b^4*sqrt(a + b*x^2)) + (C*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/b^(9/2)]
    @test_int [x^5*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 5, -((x^5*(a*B - (A*b - a*C)*x))/(7*a*b*(a + b*x^2)^(7/2))) - (x^4*(A*b + 6*a*C - 5*b*B*x))/(35*a*b^2*(a + b*x^2)^(5/2)) + (4*(A*b + 6*a*C))/(105*b^4*(a + b*x^2)^(3/2)) - (4*(A*b + 6*a*C))/(35*a*b^4*sqrt(a + b*x^2))]
    @test_int [x^4*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 4, -((x^4*(a*B - (A*b - a*C)*x))/(7*a*b*(a + b*x^2)^(7/2))) - (x^2*(4*a*B + (2*A*b + 5*a*C)*x))/(35*a*b^2*(a + b*x^2)^(5/2)) - (8*a*B + 3*(2*A*b + 5*a*C)*x)/(105*a*b^3*(a + b*x^2)^(3/2)) + ((2*A*b + 5*a*C)*x)/(35*a^2*b^3*sqrt(a + b*x^2))]
    @test_int [x^3*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 4, -((x^3*(a*B - (A*b - a*C)*x))/(7*a*b*(a + b*x^2)^(7/2))) - (x*(3*a*B + (3*A*b + 4*a*C)*x))/(35*a*b^2*(a + b*x^2)^(5/2)) - (2*(3*A*b + 4*a*C) - 3*b*B*x)/(105*a*b^3*(a + b*x^2)^(3/2)) + (2*B*x)/(35*a^2*b^2*sqrt(a + b*x^2))]
    @test_int [x^2*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 4, -((x^2*(a*B - (A*b - a*C)*x))/(7*a*b*(a + b*x^2)^(7/2))) - (2*a*B + (4*A*b + 3*a*C)*x)/(35*a*b^2*(a + b*x^2)^(5/2)) + ((4*A*b + 3*a*C)*x)/(105*a^2*b^2*(a + b*x^2)^(3/2)) + (2*(4*A*b + 3*a*C)*x)/(105*a^3*b^2*sqrt(a + b*x^2))]
    @test_int [x^1*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 4, -((x*(a*B - (A*b - a*C)*x))/(7*a*b*(a + b*x^2)^(7/2))) - (5*A*b + 2*a*C - b*B*x)/(35*a*b^2*(a + b*x^2)^(5/2)) + (4*B*x)/(105*a^2*b*(a + b*x^2)^(3/2)) + (8*B*x)/(105*a^3*b*sqrt(a + b*x^2))]
    @test_int [x^0*(A + B*x + C*x^2)/(a + b*x^2)^(9/2), x, 5, -((a*B - (A*b - a*C)*x)/(7*a*b*(a + b*x^2)^(7/2))) + ((6*A*b + a*C)*x)/(35*a^2*b*(a + b*x^2)^(5/2)) + (4*(6*A*b + a*C)*x)/(105*a^3*b*(a + b*x^2)^(3/2)) + (8*(6*A*b + a*C)*x)/(105*a^4*b*sqrt(a + b*x^2))]
    @test_int [(A + B*x + C*x^2)/(x^1*(a + b*x^2)^(9/2)), x, 8, (A*b - a*C + b*B*x)/(7*a*b*(a + b*x^2)^(7/2)) + (7*A + 6*B*x)/(35*a^2*(a + b*x^2)^(5/2)) + (35*A + 24*B*x)/(105*a^3*(a + b*x^2)^(3/2)) + (35*A + 16*B*x)/(35*a^4*sqrt(a + b*x^2)) - (A*atanh(sqrt(a + b*x^2)/sqrt(a)))/a^(9/2)]
    @test_int [(A + B*x + C*x^2)/(x^2*(a + b*x^2)^(9/2)), x, 8, (B - ((A*b)/a - C)*x)/(7*a*(a + b*x^2)^(7/2)) + (7*B - ((13*A*b)/a - 6*C)*x)/(35*a^2*(a + b*x^2)^(5/2)) + (35*B - 3*((29*A*b)/a - 8*C)*x)/(105*a^3*(a + b*x^2)^(3/2)) + (35*B - ((93*A*b)/a - 16*C)*x)/(35*a^4*sqrt(a + b*x^2)) - (A*sqrt(a + b*x^2))/(a^5*x) - (B*atanh(sqrt(a + b*x^2)/sqrt(a)))/a^(9/2)]
    @test_int [(A + B*x + C*x^2)/(x^3*(a + b*x^2)^(9/2)), x, 9, -((a*((A*b)/a - C) + b*B*x)/(7*a^2*(a + b*x^2)^(7/2))) - (7*(2*A*b - a*C) + 13*b*B*x)/(35*a^3*(a + b*x^2)^(5/2)) - (35*(3*A*b - a*C) + 87*b*B*x)/(105*a^4*(a + b*x^2)^(3/2)) - (35*(4*A*b - a*C) + 93*b*B*x)/(35*a^5*sqrt(a + b*x^2)) - (A*sqrt(a + b*x^2))/(2*a^5*x^2) - (B*sqrt(a + b*x^2))/(a^5*x) + ((9*A*b - 2*a*C)*atanh(sqrt(a + b*x^2)/sqrt(a)))/(2*a^(11/2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c*x)^m*(A+B*x+C*x^2)*(a+b*x^2)^p*when*m*symbolic=#

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

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

    @test_int [(c*x)^m*(A + 0*x + 0*x^2)/(a + b*x^2), x, 2, (A*(c*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*c*(1 + m))]
    @test_int [(c*x)^m*(A + B*x + 0*x^2)/(a + b*x^2), x, 3, (A*(c*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*c*(1 + m)) + (B*(c*x)^(2 + m)*HypergeometricFunctions._₂F₁(1, (2 + m)/2, (4 + m)/2, -((b*x^2)/a)))/(a*c^2*(2 + m))]
    @test_int [(c*x)^m*(A + 0*x + C*x^2)/(a + b*x^2), x, 2, (C*(c*x)^(1 + m))/(b*c*(1 + m)) + ((A*b - a*C)*(c*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*b*c*(1 + m))]
    @test_int [(c*x)^m*(A + B*x + C*x^2)/(a + b*x^2), x, 5, (C*(c*x)^(1 + m))/(b*c*(1 + m)) + ((A*b - a*C)*(c*x)^(1 + m)*HypergeometricFunctions._₂F₁(1, (1 + m)/2, (3 + m)/2, -((b*x^2)/a)))/(a*b*c*(1 + m)) + (B*(c*x)^(2 + m)*HypergeometricFunctions._₂F₁(1, (2 + m)/2, (4 + m)/2, -((b*x^2)/a)))/(a*c^2*(2 + m))]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(A+B*x+C*x^2+D*x^3)*(a+b*x^2)^p=#

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

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

    @test_int [x^3*(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3), x, 2, (a^2*A*x^4)/4 + (a^2*B*x^5)/5 + (a*(2*A*b + a*C)*x^6)/6 + (a*(2*b*B + a*D)*x^7)/7 + (b*(A*b + 2*a*C)*x^8)/8 + (b*(b*B + 2*a*D)*x^9)/9 + (b^2*C*x^10)/10 + (b^2*D*x^11)/11]
    @test_int [x^2*(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3), x, 2, (a^2*A*x^3)/3 + (a^2*B*x^4)/4 + (a*(2*A*b + a*C)*x^5)/5 + (a*(2*b*B + a*D)*x^6)/6 + (b*(A*b + 2*a*C)*x^7)/7 + (b*(b*B + 2*a*D)*x^8)/8 + (b^2*C*x^9)/9 + (b^2*D*x^10)/10]
    @test_int [x^1*(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3), x, 3, (a^2*B*x^3)/3 + (a^2*C*x^4)/4 + (a*(2*b*B + a*D)*x^5)/5 + (a*b*C*x^6)/3 + (b*(b*B + 2*a*D)*x^7)/7 + (b^2*C*x^8)/8 + (b^2*D*x^9)/9 + (A*(a + b*x^2)^3)/(6*b)]
    @test_int [x^0*(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3), x, 3, a^2*A*x + (a*(2*A*b + a*C)*x^3)/3 + (a^2*D*x^4)/4 + (b*(A*b + 2*a*C)*x^5)/5 + (a*b*D*x^6)/3 + (b^2*C*x^7)/7 + (b^2*D*x^8)/8 + (B*(a + b*x^2)^3)/(6*b)]
    @test_int [(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3)/x^1, x, 3, a^2*B*x + a*A*b*x^2 + (a*(2*b*B + a*D)*x^3)/3 + (A*b^2*x^4)/4 + (b*(b*B + 2*a*D)*x^5)/5 + (b^2*D*x^7)/7 + (C*(a + b*x^2)^3)/(6*b) + a^2*A*log(x)]
    @test_int [(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3)/x^2, x, 3, -((a^2*A)/x) + a*(2*A*b + a*C)*x + a*b*B*x^2 + (b*(A*b + 2*a*C)*x^3)/3 + (b^2*B*x^4)/4 + (b^2*C*x^5)/5 + (D*(a + b*x^2)^3)/(6*b) + a^2*B*log(x)]
    @test_int [(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3)/x^3, x, 2, -(a^2*A)/(2*x^2) - (a^2*B)/x + a*(2*b*B + a*D)*x + (b*(A*b + 2*a*C)*x^2)/2 + (b*(b*B + 2*a*D)*x^3)/3 + (b^2*C*x^4)/4 + (b^2*D*x^5)/5 + a*(2*A*b + a*C)*log(x)]
    @test_int [(a + b*x^2)^2*(A + B*x + C*x^2 + D*x^3)/x^4, x, 2, -(a^2*A)/(3*x^3) - (a^2*B)/(2*x^2) - (a*(2*A*b + a*C))/x + b*(A*b + 2*a*C)*x + (b*(b*B + 2*a*D)*x^2)/2 + (b^2*C*x^3)/3 + (b^2*D*x^4)/4 + a*(2*b*B + a*D)*log(x)]

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

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

    @test_int [x^4*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2), x, 5, -((a*(A*b - a*C)*x)/b^3) - (a*(b*B - a*D)*x^2)/(2*b^3) + ((A*b - a*C)*x^3)/(3*b^2) + ((b*B - a*D)*x^4)/(4*b^2) + (C*x^5)/(5*b) + (D*x^6)/(6*b) + (a^(3/2)*(A*b - a*C)*atan((sqrt(b)*x)/sqrt(a)))/b^(7/2) + (a^2*(b*B - a*D)*log(a + b*x^2))/(2*b^4)]
    @test_int [x^3*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2), x, 5, -((a*(b*B - a*D)*x)/b^3) + ((A*b - a*C)*x^2)/(2*b^2) + ((b*B - a*D)*x^3)/(3*b^2) + (C*x^4)/(4*b) + (D*x^5)/(5*b) + (a^(3/2)*(b*B - a*D)*atan((sqrt(b)*x)/sqrt(a)))/b^(7/2) - (a*(A*b - a*C)*log(a + b*x^2))/(2*b^3)]
    @test_int [x^2*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2), x, 5, ((A*b - a*C)*x)/b^2 + ((b*B - a*D)*x^2)/(2*b^2) + (C*x^3)/(3*b) + (D*x^4)/(4*b) - (sqrt(a)*(A*b - a*C)*atan((sqrt(b)*x)/sqrt(a)))/b^(5/2) - (a*(b*B - a*D)*log(a + b*x^2))/(2*b^3)]
    @test_int [x^1*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2), x, 5, ((b*B - a*D)*x)/b^2 + (C*x^2)/(2*b) + (D*x^3)/(3*b) - (sqrt(a)*(b*B - a*D)*atan((sqrt(b)*x)/sqrt(a)))/b^(5/2) + ((A*b - a*C)*log(a + b*x^2))/(2*b^2)]
    @test_int [x^0*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2), x, 5, (C*x)/b + (D*x^2)/(2*b) + ((A*b - a*C)*atan((sqrt(b)*x)/sqrt(a)))/(sqrt(a)*b^(3/2)) + ((b*B - a*D)*log(a + b*x^2))/(2*b^2)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^1*(a + b*x^2)), x, 5, (D*x)/b + ((b*B - a*D)*atan((sqrt(b)*x)/sqrt(a)))/(sqrt(a)*b^(3/2)) + (A*log(x))/a - ((A*b - a*C)*log(a + b*x^2))/(2*a*b)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^2*(a + b*x^2)), x, 5, -(A/(a*x)) - ((A*b - a*C)*atan((sqrt(b)*x)/sqrt(a)))/(a^(3/2)*sqrt(b)) + (B*log(x))/a - ((b*B - a*D)*log(a + b*x^2))/(2*a*b)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^3*(a + b*x^2)), x, 5, -A/(2*a*x^2) - B/(a*x) - ((b*B - a*D)*atan((sqrt(b)*x)/sqrt(a)))/(a^(3/2)*sqrt(b)) - ((A*b - a*C)*log(x))/a^2 + ((A*b - a*C)*log(a + b*x^2))/(2*a^2)]

    @test_int [x^4*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^2, x, 6, ((3*A*b - 5*a*C)*x)/(2*b^3) + ((2*b*B - 3*a*D)*x^2)/(2*b^3) - ((3*A*b - 5*a*C)*x^3)/(6*a*b^2) + (D*x^4)/(4*b^2) - (x^4*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(2*a*b*(a + b*x^2)) - (sqrt(a)*(3*A*b - 5*a*C)*atan((sqrt(b)*x)/sqrt(a)))/(2*b^(7/2)) - (a*(2*b*B - 3*a*D)*log(a + b*x^2))/(2*b^4)]
    @test_int [x^3*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^2, x, 6, ((3*b*B - 5*a*D)*x)/(2*b^3) - ((A*b - 2*a*C)*x^2)/(2*a*b^2) + (D*x^3)/(3*b^2) - (x^3*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(2*a*b*(a + b*x^2)) - (sqrt(a)*(3*b*B - 5*a*D)*atan((sqrt(b)*x)/sqrt(a)))/(2*b^(7/2)) + ((A*b - 2*a*C)*log(a + b*x^2))/(2*b^3)]
    @test_int [x^2*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^2, x, 6, -((A*b - 3*a*C)*x)/(2*a*b^2) + (D*x^2)/(2*b^2) - (x^2*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(2*a*b*(a + b*x^2)) + ((A*b - 3*a*C)*atan((sqrt(b)*x)/sqrt(a)))/(2*sqrt(a)*b^(5/2)) + ((b*B - 2*a*D)*log(a + b*x^2))/(2*b^3)]
    @test_int [x^1*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^2, x, 6, (D*x)/b^2 - (x*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(2*a*b*(a + b*x^2)) + ((b*B - 3*a*D)*atan((sqrt(b)*x)/sqrt(a)))/(2*sqrt(a)*b^(5/2)) + (C*log(a + b*x^2))/(2*b^2)]
    @test_int [x^0*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^2, x, 4, -(a*(B - (a*D)/b) - (A*b - a*C)*x)/(2*a*b*(a + b*x^2)) + ((A*b + a*C)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(3/2)*b^(3/2)) + (D*log(a + b*x^2))/(2*b^2)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^1*(a + b*x^2)^2), x, 6, (A*b - a*C + (b*B - a*D)*x)/(2*a*b*(a + b*x^2)) + ((b*B + a*D)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(3/2)*b^(3/2)) + (A*log(x))/a^2 - (A*log(a + b*x^2))/(2*a^2)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^2*(a + b*x^2)^2), x, 6, -(A/(a^2*x)) + (b*B - a*D - b*((A*b)/a - C)*x)/(2*a*b*(a + b*x^2)) - ((3*A*b - a*C)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(5/2)*sqrt(b)) + (B*log(x))/a^2 - (B*log(a + b*x^2))/(2*a^2)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^3*(a + b*x^2)^2), x, 6, -A/(2*a^2*x^2) - B/(a^2*x) - ((A*b)/a - C + ((b*B)/a - D)*x)/(2*a*(a + b*x^2)) - ((3*b*B - a*D)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(5/2)*sqrt(b)) - ((2*A*b - a*C)*log(x))/a^3 + ((2*A*b - a*C)*log(a + b*x^2))/(2*a^3)]

    @test_int [x^4*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^3, x, 7, (-3*(A*b - 5*a*C)*x)/(8*a*b^3) - ((b*B - 3*a*D)*x^2)/(2*a*b^3) - (x^4*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(4*a*b*(a + b*x^2)^2) + (x^3*(A*b - 5*a*C + 4*(b*B - 2*a*D)*x))/(8*a*b^2*(a + b*x^2)) + (3*(A*b - 5*a*C)*atan((sqrt(b)*x)/sqrt(a)))/(8*sqrt(a)*b^(7/2)) + ((b*B - 3*a*D)*log(a + b*x^2))/(2*b^4)]
    @test_int [x^3*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^3, x, 6, (-3*(b*B - 5*a*D)*x)/(8*a*b^3) - (x^3*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(4*a*b*(a + b*x^2)^2) - (x^2*(4*a*C - (3*b*B - 7*a*D)*x))/(8*a*b^2*(a + b*x^2)) + (3*(b*B - 5*a*D)*atan((sqrt(b)*x)/sqrt(a)))/(8*sqrt(a)*b^(7/2)) + (C*log(a + b*x^2))/(2*b^3)]
    @test_int [x^2*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^3, x, 5, -(x^2*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(4*a*b*(a + b*x^2)^2) - (x*(A*b + 3*a*C - 2*(b*B - 3*a*D)*x))/(8*a*b^2*(a + b*x^2)) + ((A*b + 3*a*C)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(3/2)*b^(5/2)) + (D*log(a + b*x^2))/(2*b^3)]
    @test_int [x^1*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^3, x, 4, -(x*(a*(B - (a*D)/b) - (A*b - a*C)*x))/(4*a*b*(a + b*x^2)^2) - (2*(A*b + a*C) - (b*B - 5*a*D)*x)/(8*a*b^2*(a + b*x^2)) + ((b*B + 3*a*D)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(3/2)*b^(5/2))]
    @test_int [x^0*(A + B*x + C*x^2 + D*x^3)/(a + b*x^2)^3, x, 3, -(a*(B - (a*D)/b) - (A*b - a*C)*x)/(4*a*b*(a + b*x^2)^2) - (4*a^2*D - b*(3*A*b + a*C)*x)/(8*a^2*b^2*(a + b*x^2)) + ((3*A*b + a*C)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(5/2)*b^(3/2))]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^1*(a + b*x^2)^3), x, 7, (A*b - a*C + (b*B - a*D)*x)/(4*a*b*(a + b*x^2)^2) + (4*A*b + (3*b*B + a*D)*x)/(8*a^2*b*(a + b*x^2)) + ((3*b*B + a*D)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(5/2)*b^(3/2)) + (A*log(x))/a^3 - (A*log(a + b*x^2))/(2*a^3)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^2*(a + b*x^2)^3), x, 7, -(A/(a^3*x)) + (b*B - a*D - b*((A*b)/a - C)*x)/(4*a*b*(a + b*x^2)^2) + (4*B - ((7*A*b)/a - 3*C)*x)/(8*a^2*(a + b*x^2)) - (3*(5*A*b - a*C)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(7/2)*sqrt(b)) + (B*log(x))/a^3 - (B*log(a + b*x^2))/(2*a^3)]
    @test_int [(A + B*x + C*x^2 + D*x^3)/(x^3*(a + b*x^2)^3), x, 7, -A/(2*a^3*x^2) - B/(a^3*x) - ((A*b)/a - C + ((b*B)/a - D)*x)/(4*a*(a + b*x^2)^2) - (4*(2*A*b - a*C) + (7*b*B - 3*a*D)*x)/(8*a^3*(a + b*x^2)) - (3*(5*b*B - a*D)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(7/2)*sqrt(b)) - ((3*A*b - a*C)*log(x))/a^4 + ((3*A*b - a*C)*log(a + b*x^2))/(2*a^4)]

    @test_int [(-x + 4*x^3)/(5 + x^2)^2, x, 4, 21/(2*(5 + x^2)) + 2*log(5 + x^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c*x)^m*P3(x^2)*(a+b*x^2)^(p/2)=#

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

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

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

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

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

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

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

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

    #= ::Title:: =#
    #=Integrands*of*the*form*(c*x)^m*Pq(x^2)*(a+b*x^2)^p=#

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

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

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

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

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

    @test_int [x^6*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2, x, 5, -((a*(5*b^3*c - 7*a*b^2*d + 9*a^2*b*e - 11*a^3*f)*x)/(2*b^6)) + ((5*b^3*c - 7*a*b^2*d + 9*a^2*b*e - 11*a^3*f)*x^3)/(6*b^5) - ((5*b^3*c - 7*a*b^2*d + 9*a^2*b*e - 11*a^3*f)*x^5)/(10*a*b^4) + ((b*e - 2*a*f)*x^7)/(7*b^3) + (f*x^9)/(9*b^2) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x^7)/(2*a*(a + b*x^2)) + (a^(3/2)*(5*b^3*c - 7*a*b^2*d + 9*a^2*b*e - 11*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*b^(13/2))]
    @test_int [x^4*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2, x, 5, ((3*b^3*c - 5*a*b^2*d + 7*a^2*b*e - 9*a^3*f)*x)/(2*b^5) - ((3*b^3*c - 5*a*b^2*d + 7*a^2*b*e - 9*a^3*f)*x^3)/(6*a*b^4) + ((b*e - 2*a*f)*x^5)/(5*b^3) + (f*x^7)/(7*b^2) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x^5)/(2*a*(a + b*x^2)) - (sqrt(a)*(3*b^3*c - 5*a*b^2*d + 7*a^2*b*e - 9*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*b^(11/2))]
    @test_int [x^2*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2, x, 5, -(((b^3*c - 3*a*b^2*d + 5*a^2*b*e - 7*a^3*f)*x)/(2*a*b^4)) + ((b*e - 2*a*f)*x^3)/(3*b^3) + (f*x^5)/(5*b^2) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x^3)/(2*a*(a + b*x^2)) + ((b^3*c - 3*a*b^2*d + 5*a^2*b*e - 7*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*sqrt(a)*b^(9/2))]
    @test_int [x^0*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2, x, 4, ((b*e - 2*a*f)*x)/b^3 + (f*x^3)/(3*b^2) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x)/(2*a*(a + b*x^2)) + ((b^3*c + a*b^2*d - 3*a^2*b*e + 5*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(3/2)*b^(7/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2/x^2, x, 4, -(c/(a^2*x)) + (f*x)/b^2 - (((b*c)/a - d + (a*e)/b - (a^2*f)/b^2)*x)/(2*a*(a + b*x^2)) - ((3*b^3*c - a*b^2*d - a^2*b*e + 3*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(5/2)*b^(5/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2/x^4, x, 4, -(c/(3*a^2*x^3)) + (2*b*c - a*d)/(a^3*x) + (((b^2*c)/a^2 - (b*d)/a + e - (a*f)/b)*x)/(2*a*(a + b*x^2)) + ((5*b^3*c - 3*a*b^2*d + a^2*b*e + a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(7/2)*b^(3/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2/x^6, x, 4, -(c/(5*a^2*x^5)) + (2*b*c - a*d)/(3*a^3*x^3) - (3*b^2*c - 2*a*b*d + a^2*e)/(a^4*x) - ((b^3*c - a*b^2*d + a^2*b*e - a^3*f)*x)/(2*a^4*(a + b*x^2)) - ((7*b^3*c - 5*a*b^2*d + 3*a^2*b*e - a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(9/2)*sqrt(b))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2/x^8, x, 4, -(c/(7*a^2*x^7)) + (2*b*c - a*d)/(5*a^3*x^5) - (3*b^2*c - 2*a*b*d + a^2*e)/(3*a^4*x^3) + (4*b^3*c - 3*a*b^2*d + 2*a^2*b*e - a^3*f)/(a^5*x) + (b*(b^3*c - a*b^2*d + a^2*b*e - a^3*f)*x)/(2*a^5*(a + b*x^2)) + (sqrt(b)*(9*b^3*c - 7*a*b^2*d + 5*a^2*b*e - 3*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(11/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^2/x^10, x, 4, -(c/(9*a^2*x^9)) + (2*b*c - a*d)/(7*a^3*x^7) - (3*b^2*c - 2*a*b*d + a^2*e)/(5*a^4*x^5) + (4*b^3*c - 3*a*b^2*d + 2*a^2*b*e - a^3*f)/(3*a^5*x^3) - (b*(5*b^3*c - 4*a*b^2*d + 3*a^2*b*e - 2*a^3*f))/(a^6*x) - (b^2*(b^3*c - a*b^2*d + a^2*b*e - a^3*f)*x)/(2*a^6*(a + b*x^2)) - (b^(3/2)*(11*b^3*c - 9*a*b^2*d + 7*a^2*b*e - 5*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(2*a^(13/2))]

    @test_int [x^8*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3, x, 6, -((a*(15*b^3*c - 27*a*b^2*d + 43*a^2*b*e - 63*a^3*f)*x)/(4*b^7)) + ((5*b^3*c - 9*a*b^2*d + 15*a^2*b*e - 23*a^3*f)*x^3)/(6*b^6) - ((5*b^3*c - 9*a*b^2*d + 17*a^2*b*e - 29*a^3*f)*x^5)/(20*a*b^5) + ((b*e - 3*a*f)*x^7)/(7*b^4) + (f*x^9)/(9*b^3) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x^9)/(4*a*(a + b*x^2)^2) - (a^2*(5*b^3*c - 9*a*b^2*d + 13*a^2*b*e - 17*a^3*f)*x)/(8*b^7*(a + b*x^2)) + (a^(3/2)*(35*b^3*c - 63*a*b^2*d + 99*a^2*b*e - 143*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*b^(15/2))]
    @test_int [x^6*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3, x, 6, ((3*b^3*c - 7*a*b^2*d + 13*a^2*b*e - 21*a^3*f)*x)/(2*b^6) - ((3*b^3*c - 7*a*b^2*d + 15*a^2*b*e - 27*a^3*f)*x^3)/(12*a*b^5) + ((b*e - 3*a*f)*x^5)/(5*b^4) + (f*x^7)/(7*b^3) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x^7)/(4*a*(a + b*x^2)^2) + (a*(3*b^3*c - 7*a*b^2*d + 11*a^2*b*e - 15*a^3*f)*x)/(8*b^6*(a + b*x^2)) - (sqrt(a)*(15*b^3*c - 35*a*b^2*d + 63*a^2*b*e - 99*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*b^(13/2))]
    @test_int [x^4*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3, x, 6, -(((b^3*c - 5*a*b^2*d + 13*a^2*b*e - 25*a^3*f)*x)/(4*a*b^5)) + ((b*e - 3*a*f)*x^3)/(3*b^4) + (f*x^5)/(5*b^3) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x^5)/(4*a*(a + b*x^2)^2) - ((b^3*c - 5*a*b^2*d + 9*a^2*b*e - 13*a^3*f)*x)/(8*b^5*(a + b*x^2)) + ((3*b^3*c - 15*a*b^2*d + 35*a^2*b*e - 63*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*sqrt(a)*b^(11/2))]
    @test_int [x^2*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3, x, 6, ((b*e - 3*a*f)*x)/b^4 + (f*x^3)/(3*b^3) + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x^3)/(4*a*(a + b*x^2)^2) - ((b^3*c + 3*a*b^2*d - 7*a^2*b*e + 11*a^3*f)*x)/(8*a*b^4*(a + b*x^2)) + ((b^3*c + 3*a*b^2*d - 15*a^2*b*e + 35*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(3/2)*b^(9/2))]
    @test_int [x^0*(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3, x, 4, (f*x)/b^3 + ((c - (a*(b^2*d - a*b*e + a^2*f))/b^3)*x)/(4*a*(a + b*x^2)^2) + ((3*b^3*c + a*b^2*d - 5*a^2*b*e + 9*a^3*f)*x)/(8*a^2*b^3*(a + b*x^2)) + ((3*b^3*c + a*b^2*d + 3*a^2*b*e - 15*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(5/2)*b^(7/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3/x^2, x, 4, -(c/(a^3*x)) - (((b*c)/a - d + (a*e)/b - (a^2*f)/b^2)*x)/(4*a*(a + b*x^2)^2) - ((7*b^3*c - 3*a*b^2*d - a^2*b*e + 5*a^3*f)*x)/(8*a^3*b^2*(a + b*x^2)) - ((15*b^3*c - 3*a*b^2*d - a^2*b*e - 3*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(7/2)*b^(5/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3/x^4, x, 5, -(c/(3*a^3*x^3)) + (3*b*c - a*d)/(a^4*x) + (((b^2*c)/a^2 - (b*d)/a + e - (a*f)/b)*x)/(4*a*(a + b*x^2)^2) + ((11*b^3*c - 7*a*b^2*d + 3*a^2*b*e + a^3*f)*x)/(8*a^4*b*(a + b*x^2)) + ((35*b^3*c - 15*a*b^2*d + 3*a^2*b*e + a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(9/2)*b^(3/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3/x^6, x, 5, -(c/(5*a^3*x^5)) + (3*b*c - a*d)/(3*a^4*x^3) - (6*b^2*c - 3*a*b*d + a^2*e)/(a^5*x) - ((b^3*c - a*b^2*d + a^2*b*e - a^3*f)*x)/(4*a^4*(a + b*x^2)^2) - ((15*b^3*c - 11*a*b^2*d + 7*a^2*b*e - 3*a^3*f)*x)/(8*a^5*(a + b*x^2)) - ((63*b^3*c - 35*a*b^2*d + 15*a^2*b*e - 3*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(11/2)*sqrt(b))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3/x^8, x, 5, -(c/(7*a^3*x^7)) + (3*b*c - a*d)/(5*a^4*x^5) - (6*b^2*c - 3*a*b*d + a^2*e)/(3*a^5*x^3) + (10*b^3*c - 6*a*b^2*d + 3*a^2*b*e - a^3*f)/(a^6*x) + (b*(b^3*c - a*b^2*d + a^2*b*e - a^3*f)*x)/(4*a^5*(a + b*x^2)^2) + (b*(19*b^3*c - 15*a*b^2*d + 11*a^2*b*e - 7*a^3*f)*x)/(8*a^6*(a + b*x^2)) + (sqrt(b)*(99*b^3*c - 63*a*b^2*d + 35*a^2*b*e - 15*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(13/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(a + b*x^2)^3/x^10, x, 5, -(c/(9*a^3*x^9)) + (3*b*c - a*d)/(7*a^4*x^7) - (6*b^2*c - 3*a*b*d + a^2*e)/(5*a^5*x^5) + (10*b^3*c - 6*a*b^2*d + 3*a^2*b*e - a^3*f)/(3*a^6*x^3) - (b*(15*b^3*c - 10*a*b^2*d + 6*a^2*b*e - 3*a^3*f))/(a^7*x) - (b^2*(b^3*c - a*b^2*d + a^2*b*e - a^3*f)*x)/(4*a^6*(a + b*x^2)^2) - (b^2*(23*b^3*c - 19*a*b^2*d + 15*a^2*b*e - 11*a^3*f)*x)/(8*a^7*(a + b*x^2)) - (b^(3/2)*(143*b^3*c - 99*a*b^2*d + 63*a^2*b*e - 35*a^3*f)*atan((sqrt(b)*x)/sqrt(a)))/(8*a^(15/2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c*x)^m*P3(x^2)*(a+b*x^2)^(p/2)=#

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

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

    @test_int [x^5*(c + d*x^2 + e*x^4 + f*x^6)/sqrt(a + b*x^2), x, 3, (a^2*(b^3*c - a*b^2*d + a^2*b*e - a^3*f)*sqrt(a + b*x^2))/b^6 - (a*(2*b^3*c - 3*a*b^2*d + 4*a^2*b*e - 5*a^3*f)*(a + b*x^2)^(3/2))/(3*b^6) + ((b^3*c - 3*a*b^2*d + 6*a^2*b*e - 10*a^3*f)*(a + b*x^2)^(5/2))/(5*b^6) + ((b^2*d - 4*a*b*e + 10*a^2*f)*(a + b*x^2)^(7/2))/(7*b^6) + ((b*e - 5*a*f)*(a + b*x^2)^(9/2))/(9*b^6) + (f*(a + b*x^2)^(11/2))/(11*b^6)]
    @test_int [x^3*(c + d*x^2 + e*x^4 + f*x^6)/sqrt(a + b*x^2), x, 3, -((a*(b^3*c - a*b^2*d + a^2*b*e - a^3*f)*sqrt(a + b*x^2))/b^5) + ((b^3*c - 2*a*b^2*d + 3*a^2*b*e - 4*a^3*f)*(a + b*x^2)^(3/2))/(3*b^5) + ((b^2*d - 3*a*b*e + 6*a^2*f)*(a + b*x^2)^(5/2))/(5*b^5) + ((b*e - 4*a*f)*(a + b*x^2)^(7/2))/(7*b^5) + (f*(a + b*x^2)^(9/2))/(9*b^5)]
    @test_int [x^1*(c + d*x^2 + e*x^4 + f*x^6)/sqrt(a + b*x^2), x, 3, ((b^3*c - a*b^2*d + a^2*b*e - a^3*f)*sqrt(a + b*x^2))/b^4 + ((b^2*d - 2*a*b*e + 3*a^2*f)*(a + b*x^2)^(3/2))/(3*b^4) + ((b*e - 3*a*f)*(a + b*x^2)^(5/2))/(5*b^4) + (f*(a + b*x^2)^(7/2))/(7*b^4)]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^1*sqrt(a + b*x^2)), x, 5, ((b^2*d - a*b*e + a^2*f)*sqrt(a + b*x^2))/b^3 + ((b*e - 2*a*f)*(a + b*x^2)^(3/2))/(3*b^3) + (f*(a + b*x^2)^(5/2))/(5*b^3) - (c*atanh(sqrt(a + b*x^2)/sqrt(a)))/sqrt(a)]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^3*sqrt(a + b*x^2)), x, 6, ((b*e - a*f)*sqrt(a + b*x^2))/b^2 - (c*sqrt(a + b*x^2))/(2*a*x^2) + (f*(a + b*x^2)^(3/2))/(3*b^2) + ((b*c - 2*a*d)*atanh(sqrt(a + b*x^2)/sqrt(a)))/(2*a^(3/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^5*sqrt(a + b*x^2)), x, 6, (f*sqrt(a + b*x^2))/b - (c*sqrt(a + b*x^2))/(4*a*x^4) + ((3*b*c - 4*a*d)*sqrt(a + b*x^2))/(8*a^2*x^2) - ((3*b^2*c - 4*a*b*d + 8*a^2*e)*atanh(sqrt(a + b*x^2)/sqrt(a)))/(8*a^(5/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^7*sqrt(a + b*x^2)), x, 6, -((c*sqrt(a + b*x^2))/(6*a*x^6)) + ((5*b*c - 6*a*d)*sqrt(a + b*x^2))/(24*a^2*x^4) - ((5*b^2*c - 6*a*b*d + 8*a^2*e)*sqrt(a + b*x^2))/(16*a^3*x^2) + ((5*b^3*c - 6*a*b^2*d + 8*a^2*b*e - 16*a^3*f)*atanh(sqrt(a + b*x^2)/sqrt(a)))/(16*a^(7/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^9*sqrt(a + b*x^2)), x, 7, -((c*sqrt(a + b*x^2))/(8*a*x^8)) + ((7*b*c - 8*a*d)*sqrt(a + b*x^2))/(48*a^2*x^6) - ((35*b^2*c - 40*a*b*d + 48*a^2*e)*sqrt(a + b*x^2))/(192*a^3*x^4) + ((35*b^3*c - 40*a*b^2*d + 48*a^2*b*e - 64*a^3*f)*sqrt(a + b*x^2))/(128*a^4*x^2) - (b*(35*b^3*c - 40*a*b^2*d + 48*a^2*b*e - 64*a^3*f)*atanh(sqrt(a + b*x^2)/sqrt(a)))/(128*a^(9/2))]

    @test_int [x^4*(c + d*x^2 + e*x^4 + f*x^6)/sqrt(a + b*x^2), x, 7, -((a*(96*b^3*c - 80*a*b^2*d + 70*a^2*b*e - 63*a^3*f)*x*sqrt(a + b*x^2))/(256*b^5)) + ((96*b^3*c - 80*a*b^2*d + 70*a^2*b*e - 63*a^3*f)*x^3*sqrt(a + b*x^2))/(384*b^4) + ((80*b^2*d - 70*a*b*e + 63*a^2*f)*x^5*sqrt(a + b*x^2))/(480*b^3) + ((10*b*e - 9*a*f)*x^7*sqrt(a + b*x^2))/(80*b^2) + (f*x^9*sqrt(a + b*x^2))/(10*b) + (a^2*(96*b^3*c - 80*a*b^2*d + 70*a^2*b*e - 63*a^3*f)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(256*b^(11/2))]
    @test_int [x^2*(c + d*x^2 + e*x^4 + f*x^6)/sqrt(a + b*x^2), x, 6, ((64*b^3*c - 48*a*b^2*d + 40*a^2*b*e - 35*a^3*f)*x*sqrt(a + b*x^2))/(128*b^4) + ((48*b^2*d - 40*a*b*e + 35*a^2*f)*x^3*sqrt(a + b*x^2))/(192*b^3) + ((8*b*e - 7*a*f)*x^5*sqrt(a + b*x^2))/(48*b^2) + (f*x^7*sqrt(a + b*x^2))/(8*b) - (a*(64*b^3*c - 48*a*b^2*d + 40*a^2*b*e - 35*a^3*f)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(128*b^(9/2))]
    @test_int [x^0*(c + d*x^2 + e*x^4 + f*x^6)/sqrt(a + b*x^2), x, 5, ((8*b^2*d - 6*a*b*e + 5*a^2*f)*x*sqrt(a + b*x^2))/(16*b^3) + ((6*b*e - 5*a*f)*x^3*sqrt(a + b*x^2))/(24*b^2) + (f*x^5*sqrt(a + b*x^2))/(6*b) + ((16*b^3*c - 8*a*b^2*d + 6*a^2*b*e - 5*a^3*f)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(16*b^(7/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^2*sqrt(a + b*x^2)), x, 6, -((c*sqrt(a + b*x^2))/(a*x)) + ((4*b*e - 3*a*f)*x*sqrt(a + b*x^2))/(8*b^2) + (f*x^3*sqrt(a + b*x^2))/(4*b) + ((8*b^2*d - 4*a*b*e + 3*a^2*f)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(8*b^(5/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^4*sqrt(a + b*x^2)), x, 6, -((c*sqrt(a + b*x^2))/(3*a*x^3)) + ((2*b*c - 3*a*d)*sqrt(a + b*x^2))/(3*a^2*x) + (f*x*sqrt(a + b*x^2))/(2*b) + ((2*b*e - a*f)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(2*b^(3/2))]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^6*sqrt(a + b*x^2)), x, 6, -((c*sqrt(a + b*x^2))/(5*a*x^5)) + ((4*b*c - 5*a*d)*sqrt(a + b*x^2))/(15*a^2*x^3) - ((8*b^2*c - 10*a*b*d + 15*a^2*e)*sqrt(a + b*x^2))/(15*a^3*x) + (f*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/sqrt(b)]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^8*sqrt(a + b*x^2)), x, 5, -((c*sqrt(a + b*x^2))/(7*a*x^7)) + ((6*b*c - 7*a*d)*sqrt(a + b*x^2))/(35*a^2*x^5) - ((24*b^2*c - 28*a*b*d + 35*a^2*e)*sqrt(a + b*x^2))/(105*a^3*x^3) + ((48*b^3*c - 56*a*b^2*d + 70*a^2*b*e - 105*a^3*f)*sqrt(a + b*x^2))/(105*a^4*x)]
    @test_int [(c + d*x^2 + e*x^4 + f*x^6)/(x^10*sqrt(a + b*x^2)), x, 6, -((c*sqrt(a + b*x^2))/(9*a*x^9)) + ((8*b*c - 9*a*d)*sqrt(a + b*x^2))/(63*a^2*x^7) - ((16*b^2*c - 18*a*b*d + 21*a^2*e)*sqrt(a + b*x^2))/(105*a^3*x^5) + ((64*b^3*c - 72*a*b^2*d + 84*a^2*b*e - 105*a^3*f)*sqrt(a + b*x^2))/(315*a^4*x^3) - (2*b*(64*b^3*c - 72*a*b^2*d + 84*a^2*b*e - 105*a^3*f)*sqrt(a + b*x^2))/(315*a^5*x)]

    @test_int [x^8*(A + B*x^2 + C*x^4 + D*x^6)/(a + b*x^2)^(9/2), x, 11, ((A - (a*(b^2*B - a*b*C + a^2*D))/b^3)*x^9)/(7*a*(a + b*x^2)^(7/2)) - ((2*A*b^3 - a*(9*b^2*B - 16*a*b*C + 23*a^2*D))*x^9)/(35*a^2*b^3*(a + b*x^2)^(5/2)) - ((16*A*b^3 - 3*a*(24*b^2*B - 66*a*b*C + 143*a^2*D))*x^7)/(210*a^2*b^4*(a + b*x^2)^(3/2)) + (D*x^9)/(6*b^3*(a + b*x^2)^(3/2)) - ((16*A*b^3 - 3*a*(24*b^2*B - 66*a*b*C + 143*a^2*D))*x^5)/(30*a^2*b^5*sqrt(a + b*x^2)) - ((16*A*b^3 - 3*a*(24*b^2*B - 66*a*b*C + 143*a^2*D))*x*sqrt(a + b*x^2))/(16*a*b^7) + ((16*A*b^3 - 3*a*(24*b^2*B - 66*a*b*C + 143*a^2*D))*x^3*sqrt(a + b*x^2))/(24*a^2*b^6) + ((16*A*b^3 - 72*a*b^2*B + 198*a^2*b*C - 429*a^3*D)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(16*b^(15/2))]
    @test_int [x^6*(A + B*x^2 + C*x^4 + D*x^6)/(a + b*x^2)^(9/2), x, 10, ((A - (a*(b^2*B - a*b*C + a^2*D))/b^3)*x^7)/(7*a*(a + b*x^2)^(7/2)) + ((b^2*B - 2*a*b*C + 3*a^2*D)*x^7)/(5*a*b^3*(a + b*x^2)^(5/2)) + ((8*b^2*B - 36*a*b*C + 99*a^2*D)*x^5)/(60*a*b^4*(a + b*x^2)^(3/2)) + (D*x^7)/(4*b^3*(a + b*x^2)^(3/2)) + ((8*b^2*B - 36*a*b*C + 99*a^2*D)*x^3)/(12*a*b^5*sqrt(a + b*x^2)) - ((8*b^2*B - 36*a*b*C + 99*a^2*D)*x*sqrt(a + b*x^2))/(8*a*b^6) + ((8*b^2*B - 36*a*b*C + 99*a^2*D)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(8*b^(13/2))]
    @test_int [x^4*(A + B*x^2 + C*x^4 + D*x^6)/(a + b*x^2)^(9/2), x, 9, ((A - (a*(b^2*B - a*b*C + a^2*D))/b^3)*x^5)/(7*a*(a + b*x^2)^(7/2)) + ((2*A*b^3 + a*(5*b^2*B - 12*a*b*C + 19*a^2*D))*x^5)/(35*a^2*b^3*(a + b*x^2)^(5/2)) + (a*(b*C - 3*a*D)*x)/(3*b^5*(a + b*x^2)^(3/2)) - ((4*b*C - 15*a*D)*x)/(3*b^5*sqrt(a + b*x^2)) + (D*x*sqrt(a + b*x^2))/(2*b^5) + ((2*b*C - 9*a*D)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(2*b^(11/2))]
    @test_int [x^2*(A + B*x^2 + C*x^4 + D*x^6)/(a + b*x^2)^(9/2), x, 8, -((a^3*D*x)/(b^4*(a + b*x^2)^(7/2))) + ((A*b^3 - 10*a^3*D)*x^3)/(3*a*b^3*(a + b*x^2)^(7/2)) + ((4*A*b^3 + 3*a*b^2*B - 58*a^3*D)*x^5)/(15*a^2*b^2*(a + b*x^2)^(7/2)) + ((8*A*b^3 + 6*a*b^2*B + 15*a^2*b*C - 176*a^3*D)*x^7)/(105*a^3*b*(a + b*x^2)^(7/2)) + (D*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/b^(9/2), ((A - (a*(b^2*B - a*b*C + a^2*D))/b^3)*x^3)/(7*a*(a + b*x^2)^(7/2)) + ((4*A*b^3 + a*(3*b^2*B - 10*a*b*C + 17*a^2*D))*x^3)/(35*a^2*b^3*(a + b*x^2)^(5/2)) + ((8*A*b^3 + a*(6*b^2*B + 15*a*b*C - 71*a^2*D))*x^3)/(105*a^3*b^3*(a + b*x^2)^(3/2)) - (D*x)/(b^4*sqrt(a + b*x^2)) + (D*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/b^(9/2)]
    @test_int [x^0*(A + B*x^2 + C*x^4 + D*x^6)/(a + b*x^2)^(9/2), x, 5, (A*x)/(a*(a + b*x^2)^(7/2)) + ((6*A*b + a*B)*x^3)/(3*a^2*(a + b*x^2)^(7/2)) + ((24*A*b^2 + a*(4*b*B + 3*a*C))*x^5)/(15*a^3*(a + b*x^2)^(7/2)) + ((48*A*b^3 + a*(8*b^2*B + 6*a*b*C + 15*a^2*D))*x^7)/(105*a^4*(a + b*x^2)^(7/2))]
    @test_int [(A + B*x^2 + C*x^4 + D*x^6)/(x^2*(a + b*x^2)^(9/2)), x, 6, -(A/(a*x*(a + b*x^2)^(7/2))) - ((8*A*b - a*B)*x)/(a^2*(a + b*x^2)^(7/2)) - ((48*A*b^2 - a*(6*b*B + a*C))*x^3)/(3*a^3*(a + b*x^2)^(7/2)) - ((4*b*(48*A*b^2 - a*(6*b*B + a*C)) - 3*a^3*D)*x^5)/(15*a^4*(a + b*x^2)^(7/2)) - (2*b*(4*b*(48*A*b^2 - a*(6*b*B + a*C)) - 3*a^3*D)*x^7)/(105*a^5*(a + b*x^2)^(7/2))]
    @test_int [(A + B*x^2 + C*x^4 + D*x^6)/(x^4*(a + b*x^2)^(9/2)), x, 7, -(A/(3*a*x^3*(a + b*x^2)^(7/2))) + (10*A*b - 3*a*B)/(3*a^2*x*(a + b*x^2)^(7/2)) + ((80*A*b^2 - 3*a*(8*b*B - a*C))*x)/(3*a^3*(a + b*x^2)^(7/2)) + ((160*A*b^3 - a*(48*b^2*B - 6*a*b*C - a^2*D))*x^3)/(3*a^4*(a + b*x^2)^(7/2)) + (4*b*(160*A*b^3 - a*(48*b^2*B - 6*a*b*C - a^2*D))*x^5)/(15*a^5*(a + b*x^2)^(7/2)) + (8*b^2*(160*A*b^3 - a*(48*b^2*B - 6*a*b*C - a^2*D))*x^7)/(105*a^6*(a + b*x^2)^(7/2))]
    @test_int [(A + B*x^2 + C*x^4 + D*x^6)/(x^6*(a + b*x^2)^(9/2)), x, 8, -(A/(5*a*x^5*(a + b*x^2)^(7/2))) + (12*A*b - 5*a*B)/(15*a^2*x^3*(a + b*x^2)^(7/2)) - (24*A*b^2 - a*(10*b*B - 3*a*C))/(3*a^3*x*(a + b*x^2)^(7/2)) - ((192*A*b^3 - a*(80*b^2*B - 24*a*b*C + 3*a^2*D))*x)/(21*a^4*(a + b*x^2)^(7/2)) - (2*(192*A*b^3 - a*(80*b^2*B - 24*a*b*C + 3*a^2*D))*x)/(35*a^5*(a + b*x^2)^(5/2)) - (8*(192*A*b^3 - a*(80*b^2*B - 24*a*b*C + 3*a^2*D))*x)/(105*a^6*(a + b*x^2)^(3/2)) - (16*(192*A*b^3 - a*(80*b^2*B - 24*a*b*C + 3*a^2*D))*x)/(105*a^7*sqrt(a + b*x^2))]
    @test_int [(A + B*x^2 + C*x^4 + D*x^6)/(x^8*(a + b*x^2)^(9/2)), x, 9, -(A/(7*a*x^7*(a + b*x^2)^(7/2))) + (2*A*b - a*B)/(5*a^2*x^5*(a + b*x^2)^(7/2)) - (24*A*b^2 - a*(12*b*B - 5*a*C))/(15*a^3*x^3*(a + b*x^2)^(7/2)) + (48*A*b^3 - a*(24*b^2*B - 10*a*b*C + 3*a^2*D))/(3*a^4*x*(a + b*x^2)^(7/2)) + (8*b*(48*A*b^3 - a*(24*b^2*B - 10*a*b*C + 3*a^2*D))*x)/(21*a^5*(a + b*x^2)^(7/2)) + (16*b*(48*A*b^3 - a*(24*b^2*B - 10*a*b*C + 3*a^2*D))*x)/(35*a^6*(a + b*x^2)^(5/2)) + (64*b*(48*A*b^3 - a*(24*b^2*B - 10*a*b*C + 3*a^2*D))*x)/(105*a^7*(a + b*x^2)^(3/2)) + (128*b*(48*A*b^3 - a*(24*b^2*B - 10*a*b*C + 3*a^2*D))*x)/(105*a^8*sqrt(a + b*x^2))]
    @test_int [(A + B*x^2 + C*x^4 + D*x^6)/(x^10*(a + b*x^2)^(9/2)), x, 10, -(A/(9*a*x^9*(a + b*x^2)^(7/2))) + (16*A*b - 9*a*B)/(63*a^2*x^7*(a + b*x^2)^(7/2)) - (32*A*b^2 - 9*a*(2*b*B - a*C))/(45*a^3*x^5*(a + b*x^2)^(7/2)) + (128*A*b^3 - 3*a*(24*b^2*B - 12*a*b*C + 5*a^2*D))/(45*a^4*x^3*(a + b*x^2)^(7/2)) - (2*b*(128*A*b^3 - 3*a*(24*b^2*B - 12*a*b*C + 5*a^2*D)))/(9*a^5*x*(a + b*x^2)^(7/2)) - (16*b^2*(128*A*b^3 - 3*a*(24*b^2*B - 12*a*b*C + 5*a^2*D))*x)/(63*a^6*(a + b*x^2)^(7/2)) - (32*b^2*(128*A*b^3 - 3*a*(24*b^2*B - 12*a*b*C + 5*a^2*D))*x)/(105*a^7*(a + b*x^2)^(5/2)) - (128*b^2*(128*A*b^3 - 3*a*(24*b^2*B - 12*a*b*C + 5*a^2*D))*x)/(315*a^8*(a + b*x^2)^(3/2)) - (256*b^2*(128*A*b^3 - 3*a*(24*b^2*B - 12*a*b*C + 5*a^2*D))*x)/(315*a^9*sqrt(a + b*x^2))]

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

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

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(c*x)^m*P4(x^2)*(a+b*x^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c*x)^m*P4(x^2)*(a+b*x^2)^(p/2)=#

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

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

    @test_int [x^2*(A + B*x^2 + C*x^4 + D*x^6 + F*x^8)/(a + b*x^2)^(9/2), x, 10, ((A*b^4 - a*(b^3*B - a*b^2*C + a^2*b*D - a^3*F))*x^3)/(7*a*b^4*(a + b*x^2)^(7/2)) + ((4*A*b^4 + a*(3*b^3*B - 10*a*b^2*C + 17*a^2*b*D - 24*a^3*F))*x^3)/(35*a^2*b^4*(a + b*x^2)^(5/2)) + ((8*A*b^4 + a*(6*b^3*B + 15*a*b^2*C - 71*a^2*b*D + 162*a^3*F))*x^3)/(105*a^3*b^4*(a + b*x^2)^(3/2)) - ((b*D - 4*a*F)*x)/(b^5*sqrt(a + b*x^2)) + (F*x*sqrt(a + b*x^2))/(2*b^5) + ((2*b*D - 9*a*F)*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/(2*b^(11/2))]
    @test_int [x^0*(A + B*x^2 + C*x^4 + D*x^6 + F*x^8)/(a + b*x^2)^(9/2), x, 6, ((A*b^4 - a^4*F)*x)/(a*b^4*(a + b*x^2)^(7/2)) + ((6*A*b^4 + a*b^3*B - 10*a^4*F)*x^3)/(3*a^2*b^3*(a + b*x^2)^(7/2)) + ((24*A*b^4 + a*(4*b^3*B + 3*a*b^2*C - 58*a^3*F))*x^5)/(15*a^3*b^2*(a + b*x^2)^(7/2)) + ((48*A*b^4 + a*(8*b^3*B + 6*a*b^2*C + 15*a^2*b*D - 176*a^3*F))*x^7)/(105*a^4*b*(a + b*x^2)^(7/2)) + (F*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/b^(9/2), ((A*b^4 - a*(b^3*B - a*b^2*C + a^2*b*D - a^3*F))*x)/(7*a*b^4*(a + b*x^2)^(7/2)) + ((6*A*b^4 + a*(b^3*B - 8*a*b^2*C + 15*a^2*b*D - 22*a^3*F))*x)/(35*a^2*b^4*(a + b*x^2)^(5/2)) + ((24*A*b^4 + a*(4*b^3*B + 3*a*b^2*C - 45*a^2*b*D + 122*a^3*F))*x)/(105*a^3*b^4*(a + b*x^2)^(3/2)) + ((48*A*b^4 + a*(8*b^3*B + 6*a*b^2*C + 15*a^2*b*D - 176*a^3*F))*x)/(105*a^4*b^4*sqrt(a + b*x^2)) + (F*atanh((sqrt(b)*x)/sqrt(a + b*x^2)))/b^(9/2)]
    @test_int [(A + B*x^2 + C*x^4 + D*x^6 + F*x^8)/(x^2*(a + b*x^2)^(9/2)), x, 6, -(A/(a*x*(a + b*x^2)^(7/2))) - ((8*A*b - a*B)*x)/(a^2*(a + b*x^2)^(7/2)) - ((48*A*b^2 - a*(6*b*B + a*C))*x^3)/(3*a^3*(a + b*x^2)^(7/2)) - ((192*A*b^3 - a*(24*b^2*B + 4*a*b*C + 3*a^2*D))*x^5)/(15*a^4*(a + b*x^2)^(7/2)) - ((384*A*b^4 - a*(48*b^3*B + 8*a*b^2*C + 6*a^2*b*D + 15*a^3*F))*x^7)/(105*a^5*(a + b*x^2)^(7/2))]
    end

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.8 P(x) (c x)^m (a+b x^2)^p.jl