text:
@with_kw immutable MIQPTrajOptParams
    nsteps::Int = 10
    Δt::Float64 = 0.1
    bilinearmethod::Symbol = :Logarithmic1D
    disc_level::Int = 9
    verbose::Bool = true
end

immutable ContactPointDescription{N}
    max_vel::Float64
    kinematic_region::SimpleVRepresentation{N,Float64}
end

immutable BoxRobotWithRotation2D
    θmax::Float64
    I::Float64
    m::Float64
    g::SVector{2, Float64}
    contact_point_descriptions::Dict{Symbol, ContactPointDescription{2}}
end

# max_angle(robot::BoxRobotWithRotation2D) = robot.θmax
# moment_of_inertia(robot::BoxRobotWithRotation2D) = robot.I
# mass(robot::BoxRobotWithRotation2D) = robot.m
# gravitational_acceleration(robot::BoxRobotWithRotation2D) = robot.g
# limb_descriptions(robot::BoxRobotWithRotation2D) = robot.contact_point_descriptions

immutable ContactPointState{N, T}
    pos::SVector{N, T}
    # TODO: contact info
end

immutable ContactPointInput{N, T}
    vel::SVector{N, T}
    force::SVector{N, T}
end

immutable BoxRobotWithRotation2DState{T}
    θ::T
    ω::T
    r::SVector{2, T}
    ṙ::SVector{2, T}
    contact_point_states::Dict{Symbol, ContactPointState{2, T}}
end

immutable MIQPTrajOptDiagnostics{A<:AbstractArray}
    solvetime::Float64
    total_torque_constraint_violation::A
    num_continuous_variables::Int
    num_binary_variables::Int
end

function rms_constraint_violation(diagnostics::MIQPTrajOptDiagnostics)
    n = length(diagnostics.total_torque_constraint_violation)
    sqrt(1 / n * sum(x^2 for x in diagnostics.total_torque_constraint_violation))
end

function miqp_trajopt{E<:EnvironmentRegion}(robot::BoxRobotWithRotation2D, environment::AbstractVector{E},
        initialstate::BoxRobotWithRotation2DState, params::MIQPTrajOptParams)
    nsteps = params.nsteps
    h = params.Δt
    bilinearmethod = params.bilinearmethod
    disc_level = params.disc_level

    coords = Axis{:coord}([:x, :z])
    steps = Axis{:step}(1 : nsteps)
    contacts = Axis{:contact}(collect(keys(robot.contact_point_descriptions)))
    regions = Axis{:region}(1 : length(environment))
    ncontacts = length(contacts)
    nregions = length(regions)

    r0 = initialstate.r
    ṙ0 = initialstate.ṙ
    θ0 = initialstate.θ
    ω0 = initialstate.ω

    θmax = robot.θmax
    I = robot.I
    m = robot.m
    g = robot.g

    model = if bilinearmethod == :minlp
        Model(solver = OsilSolver(solver = "couenne"))
    else
        Model(solver = GurobiSolver(Presolve = 0, OutputFlag = Int(params.verbose)))
    end

    rs = @axis_variable(model, r[coords, steps])
    θs = @axis_variable(model, θ[steps])
    ṙs = @axis_variable(model, ṙ[coords, steps])
    ωs = @axis_variable(model, ω[steps])
    rcs = @axis_variable(model, rc[coords, contacts, steps])
    ṙcs = @axis_variable(model, ṙc[coords, contacts, steps])
    fs = @axis_variable(model, f[coords, contacts, steps])
    Fs = @axis_variable(model, F[coords, steps])
    Ts = @axis_variable(model, T[steps])
    ss = @axis_variable(model, s[coords, contacts, steps])
    wzxs = @axis_variable(model, wzx[contacts, steps])
    wxzs = @axis_variable(model, wxz[contacts, steps])
    γs = @axis_variable(model, γ[contacts, regions, steps(1 : length(steps) - 1)])
    Γs = @axis_variable(model, Γ[contacts, regions])

    # Determine contact force bounds (needed for relaxation of bilinear constraints)
    fmin, fmax = force_bounds(environment)
    for i = 1 : ncontacts, n = 1 : nsteps
        f = fs[contacts(i), steps(n)]
        setlowerbound.(f, fmin)
        setupperbound.(f, fmax)
    end

    # region constraints (4.6)
    zs = AxisArray(Array{Variable, 3}(ncontacts, nregions, nsteps), contacts, regions, steps)
    product_regions = collect(region.position * region.force for region in environment)
    for i = 1 : ncontacts, n = 1 : nsteps
        rc = rcs[contacts(i), steps(n)]
        f = fs[contacts(i), steps(n)]
        zs[contacts(i), steps(n)] = hull_reformulation(model, product_regions, [rc; f])
    end

    # contact point velocity limits
    for i = 1 : length(contacts), n = 1 : length(steps)
        ṙc = ṙcs[contacts(i), steps(n)]

        zfree = JuMP.AffExpr(0)
        for j = 1 : length(regions)
            isfree(environment[j]) && (zfree += zs[contacts(i), steps(n), regions(j)]) # TODO: no need to do this in the inner loop
        end

        max_contact_point_vel = robot.contact_point_descriptions[contacts[i]].max_vel
        @constraints(model, begin
            -max_contact_point_vel * zfree .<= ṙc
            ṙc .<= max_contact_point_vel * zfree
        end)
    end

    # define s variables (contact point locations relative to CoM) (4.13)
    for i = 1 : ncontacts, n = 1 : nsteps
        s = ss[contacts(i), steps(n)]
        rc = rcs[contacts(i), steps(n)]
        r = rs[steps(n)]
        @constraint(model, s .== rc - r)
    end

    # kinematic constraints
    # setlowerbound.(rs[:z], 0.55) # TODO
    setlowerbound.(θs, -θmax)
    setupperbound.(θs, θmax)

    # contact point positions relative to CoM
    # NOTE: adding these increased solve time a lot
    for i = 1 : ncontacts, n = 1 : nsteps
        s = ss[contacts(i), steps(n)]
        polyhedron_constraints(model, robot.contact_point_descriptions[contacts[i]].kinematic_region, s)
    end

    # centroidal dynamics (4.18)
    for n = 1 : nsteps - 1
        @constraints(model, begin
            rs[steps(n + 1)] .== rs[steps(n)] + h * ṙs[steps(n + 1)]
            θs[steps(n + 1)] == θs[steps(n)] + h * ωs[steps(n + 1)]
            ṙs[steps(n + 1)] .== ṙs[steps(n)] + h / m * Fs[steps(n + 1)]
            ωs[steps(n + 1)] == ωs[steps(n)] + h / I * Ts[steps(n + 1)]
        end)
    end

    for n = 1 : nsteps
        @constraints(model, begin
            Fs[steps(n)] .== m * g + sum(fs[steps(n), contacts(i)] for i = 1 : length(contacts))
            Ts[steps(n)] == sum(wzxs[steps(n), contacts(i)] - wxzs[steps(n), contacts(i)] for i = 1 : length(contacts))
        end)
    end

    # contact point dynamics
    for i = 1 : ncontacts, n = 1 : nsteps - 1
        @constraints(model, begin
            rcs[contacts(i), steps(n + 1)] .== rcs[contacts(i), steps(n)] + h * ṙcs[contacts(i), steps(n + 1)]
        end)
    end

    # convex bilinear term approximations (4.12, 4.13)
    for i = 1 : ncontacts
        kinematic_region = robot.contact_point_descriptions[contacts[i]].kinematic_region
        sbounds = axis_aligned_bounding_box(kinematic_region) # needed for relaxation of bilinear terms
        smin = AxisArray(sbounds[1], coords)
        smax = AxisArray(sbounds[2], coords)

        for n = 1 : nsteps
            wxz = wxzs[contacts(i), steps(n)]
            wzx = wzxs[contacts(i), steps(n)]
            s = ss[contacts(i), steps(n)]
            f = fs[contacts(i), steps(n)]
            setlowerbound.(s, smin)
            setupperbound.(s, smax)
            if bilinearmethod == :HRepConvexHull
                srange = AxisArray(linspace.(smin, smax, disc_level), coords)
                fmin, fmax = AxisArray(getlowerbound.(f), coords), AxisArray(getupperbound.(f), coords)
                piecewise_mccormick_envelope_constraints(model, s[:x], f[:z], wxz, srange[:x], fmin[:z]..fmax[:z])
                piecewise_mccormick_envelope_constraints(model, s[:z], f[:x], wzx, srange[:z], fmin[:x]..fmax[:x])
            else
                @constraint(model, wxz == s[:x] * f[:z])
                @constraint(model, wzx == s[:z] * f[:x])
            end
        end
    end

    # initial conditions
    @constraints(model, begin
        rs[steps(1)] .== r0
        θs[steps(1)] == θ0
        ṙs[steps(1)] .== ṙ0
        ωs[steps(1)] == ω0
    end)

    for i = 1 : ncontacts
        rc0 = initialstate.contact_point_states[contacts[i]].pos
        @constraint(model, rcs[contacts(i), steps(1)] .== rc0)
        # TODO: initial region assignments (maybe).
    end

    # final conditions
    @constraints(model, begin
        # rs[steps(length(steps))] .== [0.81; 0.9]
        θs[steps(length(steps))] == 0
        ṙs[steps(length(steps))] .== 0
        ωs[steps(length(steps))] == 0
    end)

    # section 4.4
    # objective function + slack parameterization
    @constraint(model, Γs .== sum(γs[steps(n)] for n = 1 : length(steps) - 1));
    for n = 1 : length(steps) - 1
        @constraints(model, begin
            -γs[steps(n)] .<= zs[steps(n + 1)] - zs[steps(n)]
            zs[steps(n + 1)] - zs[steps(n)] .<= γs[steps(n)]
        end)
    end
    setlowerbound.(γs, 0.)
    setupperbound.(γs, 1.)

    δs = rs .- r0
    @objective(model, Min, 1e3 * sum(Γ) + 1e-3 * vecdot(Fs, Fs) + 1e-5 * vecdot(fs, fs) + 1e-2 * vecdot(δs, δs) + vecdot(θs, θs))

    # relax bilinear constraints
    if bilinearmethod ∉ [:HRepConvexHull; :minlp]
        relaxbilinear!(model, method = bilinearmethod, disc_level = disc_level)
    end

    # solve
    status, solvetime = @timed solve(model)
    @assert status == :Optimal

    # compile results
    # states
    states = AxisArray(Vector{BoxRobotWithRotation2DState{Float64}}(length(steps)), steps)
    for n = 1 : nsteps
        θ = getvalue(θs[steps(n)])
        ω = getvalue(ωs[steps(n)])
        r = SVector{2}(getvalue.(rs[steps(n)]))
        ṙ = SVector{2}(getvalue.(ṙs[steps(n)]))

        contact_point_states = Dict{Symbol, ContactPointState{2, Float64}}()
        for i = 1 : ncontacts
            rc = SVector{2}(getvalue.(rcs[contacts(i), steps(n)]))
            contact_point_states[contacts[i]] = ContactPointState(rc)
        end

        states[steps(n)] = BoxRobotWithRotation2DState(θ, ω, r, ṙ, contact_point_states)
    end

    # inputs
    createinput = (contact, step) -> begin
        ṙc = SVector{2}(getvalue(ṙcs[contact, step]))
        f = SVector{2}(getvalue(fs[contact, step]))
        input = ContactPointInput(ṙc, f)
    end
    # inputs = Dict(contacts[i] => AxisArray(map(n -> createinput(contacts(i), steps(n)), 1 : nsteps), steps) for i = 1 : ncontacts)
    inputs = AxisArray(map(n -> Dict(contacts[i] => createinput(contacts(i), steps(n)) for i = 1 : ncontacts), 1 : nsteps), steps)

    # diagnostics
    zx_constraint_violations = getvalue(AxisArray(wzxs .- (ss[:z] .* fs[:x]), wzxs.axes))
    xz_constraint_violations = getvalue(AxisArray(wxzs .- (ss[:x] .* fs[:z]), wxzs.axes))
    torque_constraint_violations = AxisArray(vec(sum((zx_constraint_violations - xz_constraint_violations), 1)), steps) # TODO: use Axes

    num_continuous_variables = sum(model.colCat .== :Cont)
    num_binary_variables = sum(model.colCat .== :Bin)
    diagnostics = MIQPTrajOptDiagnostics(solvetime, torque_constraint_violations, num_continuous_variables, num_binary_variables)

    states, inputs, diagnostics
end


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module Gates

using QuantumCircuits.QCircuits.QBase
using QuantumCircuits.QCircuits.Math
using QuantumCircuits.QCircuits.Registers

import Base: show, length, inv
import QuantumCircuits.QCircuits.QBase: tomatrix, setparameters!, simplify,
       standardGateError, decompose, bindparameters!

export X, Y, Z, S, Sd, T, Td, H, CX, U3, Rx, Ry, Rz, U, Sx, Sxd,
      getqubits, getqubitsids, toU3,
      ChromosomeGate, DNAGate, Parameters, getArgs, decompose,
      ParameterT, getvalue

"Parameter type"
const ParamT = Number

mutable struct ParameterT
    value::Union{ParamT, Nothing}
end
ParameterT() = ParameterT(nothing)
ParameterT(param::ParameterT) = ParameterT(param.value)
Base.show(io::IO, p::ParameterT) = print(io, "$(p.value)")
"Comparator"
Base.:(==)(p1::ParameterT, p2::ParameterT) = p1.value == p2.value

const Parameter = Union{ParamT, ParameterT}

getvalue(p::ParamT) = p
function getvalue(p::ParameterT)
    @assert !isnothing(p.value) "Unable get unset parameter."
    return p.value
end

#const ParamIdx = Int
# "The gates parameters structure, one plece for all parameters in cirquit."
# struct Parameters
#     params::Vector{ParamT}
#
#     Parameters() = new(ParamT[])
# end
# "Function return the next free parameter."
# function getnext(p::Parameters)
#     push!(p.params, zero(ParamT))
#
#     return length(p.params)
# end
# "Get/Set index"
# Base.getindex(p::Parameters, idx::ParamIdx) = p.params[idx]
# Base.setindex!(p::Parameters, v::ParamT, idx::ParamIdx) = p.params[idx] = v

"Gate which can be use in chromosome"
#abstract type ChromosomeGate <: QuantumGate end
#abstract type DNAGate <: ChromosomeGate end

#"Single Qubit Gate"
#abstract type SingleQubitGate <: QuantumGate end

"Quantum Universal Gate"
abstract type UniversalGate <: QuantumGate end

"Comparator"
Base.:(==)(c1::QuantumGate, c2::QuantumGate) = false

"Single qubit create macro"
macro singleQubitGate(name, rawdoc)
    eval(quote
        @doc $rawdoc
        struct $name <: QuantumGate qubit::Qubit end

        Base.:(==)(g1::$name, g2::$name) = g1.qubit == g2.qubit
        Base.hash(g::$name, h::UInt) = hash(g.qubit, h)

        getArgs(g::$name) = (getid(g.qubit), )
    end)
end

# Single qubit quantum gates
@singleQubitGate(X,
raw"""
    X(qubit::Qubit)
Single-qubit Pauli-X gate (``\sigma_x``), equivalent to [U3](@ref)(``\pi,0,\pi``)
**Matrix Representation**
```math
X = \begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.x(0)
```
""")
@singleQubitGate(Y,
raw"""
    Y(qubit::Qubit)
Single-qubit Pauli-Y gate (``\sigma_y``), equivalent to [U3](@ref)(``\pi,\frac{\pi}{2},\frac{\pi}{2}``)
**Matrix Representation**
```math
Y = \begin{pmatrix}
0 & -i \\
i & 0
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.y(0)
```
""")
@singleQubitGate(Z,
raw"""
    Z(qubit::Qubit)
Single-qubit Pauli-Z gate (``\sigma_z``), equivalent to [U3](@ref)(``0,0,\pi``)
**Matrix Representation**
```math
Z = \begin{pmatrix}
1 & 0 \\
0 & -1
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.z(0)
```
""")
@singleQubitGate(S,
raw"""
    S(qubit::Qubit)
Single-qubit S gate, equivalent to [U3](@ref)(``0,0,\frac{\pi}{2}``). This
gate is also referred to a square-root of Pauli-[Z](@ref).
**Matrix Representation**
```math
S = \begin{pmatrix}
1 & 0 \\
0 & i
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.s(0)
```
""")
@singleQubitGate(Sd,
raw"""
Sd(qubit::Qubit)
Single-qubit, hermitian conjugate of the [S](@ref). This is also an alternative square root of
the [Z](@ref).
**Matrix Representation**
```math
S^{\dagger} = \begin{pmatrix}
1 & 0 \\
0 & -i
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.sd(0)
```
""")
@singleQubitGate(H,
raw"""
    H(qubit::Qubit)
Single-qubit Hadamard gate, which is a ``\pi`` rotation about the X+Z axis, thus equivalent to [U3](@ref)(``\frac{\pi}{2},0,\pi``)
**Matrix Representation**
```math
H = \frac{1}{\sqrt{2}}
        \begin{pmatrix}
            1 & 1 \\
            1 & -1
        \end{pmatrix}
```
""")
@singleQubitGate(T,
raw"""
    T(qubit::Qubit)
Single-qubit T gate, equivalent to [U3](@ref)(``0,0,\frac{\pi}{4}``). This
gate is also referred to as a ``\frac{\pi}{8}`` gate or as a fourth-root of Pauli-[Z](@ref).
**Matrix Representation**
```math
T = \begin{pmatrix}
1 & 0 \\
0 & e^{i\pi/4}
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.t(0)
```
""")
@singleQubitGate(Td,
raw"""
    Td(qubit::Qubit)
Single-qubit, hermitian conjugate of the [T](@ref). This gate is equivalent to [U3](@ref)(``0,0,-\frac{\pi}{4}``). This
gate is also referred to as the fourth-root of Pauli-[Z](@ref).
**Matrix Representation**
```math
T^{\dagger} = \begin{pmatrix}
1 & 0 \\
0 & e^{-i\pi/4}
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.td(0)
```
""")
@singleQubitGate(Sx,
raw"""
    Sx(qubit::Qubit)
Single-qubit square root of pauli-[X](@ref).
**Matrix Representation**
```math
\sqrt{X} = \frac{1}{2} \begin{pmatrix}
1 + i & 1 - i \\
1 - i & 1 + i
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.sx(0)
```
""")
@singleQubitGate(Sxd,
raw"""
    Sxd(qubit::Qubit)
Single-qubit hermitian conjugate of the square root of pauli-[X](@ref), or the [Sx](@ref).
**Matrix Representation**
```math
\sqrt{X}^{\dagger} = \frac{1}{2} \begin{pmatrix}
1 - i & 1 + i \\
1 + i & 1 - i
\end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.sxd(0)
```
""")

@doc raw"""
    CX(control::Qubit, target::Qubit)
Two-qubit controlled NOT gate with control and target on first and second qubits, respectively. This is also
called the controlled X gate.
**Circuit Representation**
```
q_0: ──■──
     ┌─┴─┐
q_1: ┤ X ├
     └───┘
```
**Matrix Representation**
```math
CX = \begin{pmatrix}
    1 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 0 & 1 & 0
    \end{pmatrix}
```

```julia
qc = QCircuit(2)
qc.cx(0, 1)
```
"""
struct CX <: QuantumGate
    control::Qubit
    target::Qubit
end
Base.:(==)(g1::CX, g2::CX) = g1.control == g2.control && g1.target == g2.target
Base.hash(g::CX, h::UInt) = hash((g.control, g.target), h)
function Base.:(>)(x::CX, y::CX)
    if x.control == y.control
        return x.target > y.target
    else
        return x.control > y.control
    end
end
getArgs(g::CX) = (getid(g.control), getid(g.target))

"Rotation gates"
abstract type RotationGate <: QuantumGate end

"Rotation macro"
macro rotationGate(name, rawdoc)
    eval(quote
        @doc $rawdoc
        mutable struct $name <: RotationGate
            qubit::Qubit
            θ::Parameter
        end
        $name(qubit::Integer) = $name(qubit, ParameterT(rand() * 2π))

        Base.:(==)(g1::$name, g2::$name) = g1.qubit == g2.qubit && g1.θ == g2.θ
        Base.hash(g::$name, h::UInt) = hash((g.qubit, g.θ), h)

        $name(qubits::Vector{<:Integer}, θ::ParamT) = [$name(q, θ) for q in qubits]
        getArgs(g::$name) = (getid(g.qubit), g.θ)
    end)
end

@rotationGate(Rx,
raw"""
    Rx(qubit::Qubit, θ::Parameter)
A single-qubit Pauli gate which represents rotation about the X axis.
**Matrix Representation**
```math
\newcommand{\th}{\frac{\theta}{2}}
Rx(\theta) = exp(-i \th X) =
    \begin{pmatrix}
        \cos{\th}   & -i\sin{\th} \\
        -i\sin{\th} & \cos{\th}
    \end{pmatrix}
```
```julia
qc = QCircuit(1)
qc.rx(0, π/2)
```

You can also create a gate without adding parameter, in that case, it will be initialized by random values and can be used in Quantum Machine Learning as the loss function parameters.
```julia
using QuantumCircuits
using QuantumCircuits.QML
using QuantumCircuits.Execute

# Expected circuit
expqc = QCircuit(1)
expqc.x(0)
expmat = tomatrix(expqc)

# The ansact
qc = QCircuit(1)
qc.rx(0)

# Find the parameters which fit the expected unitary matrix in the best way.
findparam(expmat, qc)

```
""")
@rotationGate(Ry,
raw"""
    Ry(qubit::Qubit, θ::Parameter)
A single-qubit Pauli gate which represents rotation about the Y axis.
**Matrix Representation**
```math
\newcommand{\th}{\frac{\theta}{2}}
Ry(\theta) = exp(-i \th Y) =
    \begin{pmatrix}
        \cos{\th} & -\sin{\th} \\
        \sin{\th} & \cos{\th}
    \end{pmatrix}
```
```julia
qc = QCircuit(1)
qc.ry(0, π/2)
```

You can also create a gate without adding parameter, in that case, it will be initialized by random values and can be used in Quantum Machine Learning as the loss function parameters.
```julia
using QuantumCircuits
using QuantumCircuits.QML
using QuantumCircuits.Execute

# Expected circuit
expqc = QCircuit(1)
expqc.y(0)
expmat = tomatrix(expqc)

# The ansact
qc = QCircuit(1)
qc.ry(0)

# Find the parameters which fit the expected unitary matrix in the best way.
findparam(expmat, qc)

```
""")
@rotationGate(Rz,
raw"""
    Rz(qubit::Qubit, θ::Parameter)
A single-qubit Pauli gate which represents rotation about the Z axis.

**Matrix Representation**
```math
\newcommand{\th}{\frac{\theta}{2}}
Rz(\theta) = exp(-i\th Z) =
\begin{pmatrix}
    e^{-i\th} & 0 \\
    0 & e^{i\th}
\end{pmatrix}
```julia
qc = QCircuit(1)
qc.rz(0, π/2)
```

You can also create a gate without adding parameter, in that case, it will be initialized by random values and can be used in Quantum Machine Learning as the loss function parameters.
```julia
using QuantumCircuits
using QuantumCircuits.QML
using QuantumCircuits.Execute

# Expected circuit
expqc = QCircuit(1)
expqc.z(0)
expmat = tomatrix(expqc)

# The ansact
qc = QCircuit(1)
qc.rz(0)

# Find the parameters which fit the expected unitary matrix in the best way.
findparam(expmat, qc)

```
""")

appendparams!(param, gate::RotationGate) = appendparams!(param, gate.θ)
appendRandParams!(param, gate::RotationGate) = appendRandParams!(param, gate.θ)
Base.length(gate::RotationGate) = getlength(gate.θ)
function setparameters!(gate::RotationGate, params)
    @assert length(params) == length(gate) "RotationGate has $(length(gate)) parameters but $(length(param)) was provided."
    setparameter!(gate.θ, params, 1)
end

@doc raw"""
    U3(qubit::Qubit, θ::Parameter, ϕ::Parameter, λ::Parameter)
Universal single-qubit rotation gate with three Euler angles, ``\theta``, ``\phi`` and ``\lambda``.
**Matrix Representation**
```math
\newcommand{\th}{\frac{\theta}{2}}
U3(\theta, \phi, \lambda) =
    \begin{pmatrix}
        \cos(\th)          & -e^{i\lambda}\sin(\th) \\
        e^{i\phi}\sin(\th) & e^{i(\phi+\lambda)}\cos(\th)
    \end{pmatrix}
```

```julia
qc = QCircuit(1)
qc.u3(0, 1, 2, 3)
```

You can also create a gate without adding parameters, in that case, they will be initialized by random values and can be used in Quantum Machine Learning as the loss function parameters.
```julia
using QuantumCircuits
using QuantumCircuits.QML
using QuantumCircuits.Execute

# Expected circuit
expqc = QCircuit(1)
expqc.x(0)
expmat = tomatrix(expqc)

# The ansact
qc = QCircuit(1)
qc.u3(0)

# Find the parameters which fit the expected unitary matrix in the best way.
findparam(expmat, qc)

```
"""
mutable struct U3 <: UniversalGate
    qubit::Qubit
    θ::Parameter
    ϕ::Parameter
    λ::Parameter
end
U3(qubit::Qubit) = U3(qubit, ([ParameterT(p) for p in rand(3) * 2π])...)
getArgs(g::U3) = (getid(g.qubit), g.θ, g.ϕ, g.λ)

Base.show(io::IO, g::U3) = print(io, "U3($(g.qubit), $(g.θ), $(g.ϕ), $(g.λ))")

Base.:(==)(g1::U3, g2::U3) = g1.qubit == g2.qubit && g1.θ == g2.θ && g1.ϕ == g2.ϕ && g1.λ == g2.λ
Base.hash(g::U3, h::UInt) = hash((g.qubit, g.θ, g.ϕ, g.λ), h)

bindparameters!(::QuantumGate) = nothing
bindparameters!(gate::RotationGate) = gate.θ = getvalue(gate.θ)
function bindparameters!(gate::U3)
    gate.θ = getvalue(gate.θ)
    gate.ϕ = getvalue(gate.ϕ)
    gate.λ = getvalue(gate.λ)
end

isparam(::Any) = false
isparam(::ParameterT) = true
onlyParameters(::QuantumGate) = false
onlyParameters(gate::U3) = isparam(gate.θ) && isparam(gate.ϕ) && isparam(gate.λ)

# "The universar single qubit unitary gate macro"
# macro universalGate(name)
#     eval(quote
#         mutable struct $name <: UniversalGate
#             qubit::Qubit
#             θ::ParamT
#             ϕ::ParamT
#             λ::ParamT
#         end
#         $name(qubit::Integer) = $name(qubit, (rand(3) * 2π)...)
#     end)
# end

"The universar single qubit unitary gate"
mutable struct U <: UniversalGate
    qubit::Qubit
    β::Union{ParamT, Nothing}
    γ::Union{ParamT, Nothing}
    δ::Union{ParamT, Nothing}
end
U(qubit::Integer) = U(qubit, nothing, nothing, nothing)
getArgs(g::U) = (getid(g.qubit), g.β, g.γ, g.δ)
Base.:(==)(g1::U, g2::U) = g1.qubit == g2.qubit && g1.β == g2.β && g1.γ == g2.γ && g1.δ == g2.δ
Base.hash(g::U, h::UInt) = hash((g.qubit, g.β, g.γ, g.δ), h)

function setparameter!(param::ParameterT, params, i)
    param.value = params[i]
    return i + 1
end
setparameter!(::ParamT, ::Any, i) = i

function setparameters!(gate::U3, params)
    @assert length(params) == length(gate) "U3 has $(length(gate)) parameters but $(length(params)) was provided."
        i = setparameter!(gate.θ, params, 1)
        i = setparameter!(gate.ϕ, params, i)
        i = setparameter!(gate.λ, params, i)
end

function setparameters!(gate::U, params)
    @assert length(params) == 3 "U has 3 parameters but $(length(params)) was provided."
    gate.β = params[1]
    gate.γ = params[2]
    gate.δ = params[3]
end

"Return the qubits on which operate the gate"
getqubits(gate::QuantumGate) = (gate.qubit,) # All Single Qubit Gates has qubit field
getqubits(gate::CX) = (gate.control, gate.target)

"Return the qubits ids on which operate the gate"
getqubitsids(gate::QuantumGate) = (getid(gate.qubit),) # All Single Qubit Gates has qubit field
getqubitsids(gate::CX) = (getid(gate.control), getid(gate.target))

"Show method"
Base.show(io::IO, gate::QuantumGate) = print(io, "$(typeof(gate))($(gate.qubit))")
Base.show(io::IO, gate::CX) = print(io, "CX($(gate.control), $(gate.target))")

"Convert gate to U3 gate"
toU3(gate::X) = U3(gate.qubit, π, 0, π)

"Return the gate matrix"
u3matrix(θ, ϕ, λ) = [cos(θ/2)+0im              -cis(λ) * sin(θ/2);
                     cis(ϕ) * sin(θ/2)   cis(λ + ϕ) * cos(θ/2)]
u2matrix(ϕ, λ) = u3matrix(π/2, ϕ, λ)
u1matrix(λ) = u3matrix(0, 0, λ)

# u3matrix(θ, -π/2, π/2)
Rxmatrix(θ) = [cos(θ/2)+0im -1im * sin(θ/2); -1im * sin(θ/2) cos(θ/2)]
# u3matrix(θ, 0, 0)
Rymatrix(θ) = [cos(θ/2)+0im -sin(θ/2); sin(θ/2) cos(θ/2)]
# exp(-1im * θ/2) * u1matrix(θ)
Rzmatrix(θ) = [cis(-θ/2) 0; 0 cis(θ/2)]
# Universary matrix
umatrix(β, γ, δ) = Rzmatrix(β) * Rymatrix(γ) * Rzmatrix(δ)

# Initial stata
const _0 = [1+0im; 0]
const _1 = [0im; 1]

const I0X0I = _0 * adjoint(_0)
const I1X1I = _1 * adjoint(_1)

# Gate matrixs
const Xmatrix = round.(u3matrix(π, 0, π), digits=2)
const Ymatrix = round.(u3matrix(π, π/2, π/2), digits=2)
const Zmatrix = round.(u1matrix(π), digits=2)
const Imatrix = [1.0 + 0im 0; 0 1]
const Hmatrix = u2matrix(0, π)
const Smatrix = round.(u1matrix(π/2), digits=2)
const Sdmatrix = round.(u1matrix(-π/2), digits=2)
const Tmatrix = u1matrix(π/4)
const Tdmatrix = u1matrix(-π/4)
const CX0matrix = [1.0 + 0im 0 0 0; 0 0 0 1; 0 0 1 0; 0 1 0 0]
const CX1matrix = [1.0 + 0im 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]

tomatrix(::X) = Xmatrix
tomatrix(::Y) = Ymatrix
tomatrix(::Z) = Zmatrix
tomatrix(::S) = Smatrix
tomatrix(::Sd) = Sdmatrix
tomatrix(::T) = Tmatrix
tomatrix(::Td) = Tdmatrix
tomatrix(::H) = Hmatrix
tomatrix(gate::U3) = u3matrix(getvalue(gate.θ), getvalue(gate.ϕ), getvalue(gate.λ))
function tomatrix(::U3, param)
    @assert length(param) == 3 "U3 has 3 parameters but $(length(param)) was provided."
    u3matrix(param[1], param[2], param[3])
end

tomatrix(gate::Rx) = Rxmatrix(getvalue(gate.θ))
tomatrix(gate::Ry) = Rymatrix(getvalue(gate.θ))
tomatrix(gate::Rz) = Rzmatrix(getvalue(gate.θ))
function tomatrix(::Rx, param)
    @assert length(param) == 1 "Rx has 1 parameters but $(length(param)) was provided."
    Rxmatrix(param[1])
end
function tomatrix(::Ry, param)
    @assert length(param) == 1 "Ry has 1 parameters but $(length(param)) was provided."
    Rymatrix(param[1])
end
function tomatrix(::Rz, param)
    @assert length(param) == 1 "Rz has 1 parameters but $(length(param)) was provided."
    Rzmatrix(param[1])
end

tomatrix(gate::U) = umatrix(gate.β, gate.γ, gate.δ)
function tomatrix(::U, param)
    @assert length(param) == 3 "U has 3 parameters but $(length(param)) was provided."
    umatrix(param[1], param[2], param[3])
end

function tomatrix(gate::CX)
    control = getid(gate.control)
    target = getid(gate.target)
    if abs(control - target) == 1
        if control < target
            return CX0matrix
        else
            return CX1matrix
        end
    else
        # https://quantumcomputing.stackexchange.com/questions/4252/how-to-derive-the-cnot-matrix-for-a-3-qubit-system-where-the-control-target-qu
        if control < target
            s0 = kron(Imatrix, I0X0I)
            s1 = kron(Imatrix, I1X1I)
            e0 = Imatrix
            e1 = Xmatrix
        else
            s0 = kron(Imatrix, Imatrix)
            s1 = kron(Imatrix, Xmatrix)
            e0 = I0X0I
            e1 = I1X1I
        end

        additional_update = abs(target - control) - 2
        for i in 1:additional_update
            s0 = kron(Imatrix, s0)
            s1 = kron(Imatrix, s1)
        end

        return kron(e0, s0) + kron(e1, s1)
    end
end

function tomatrix(qubits::Integer, gate::QuantumGate, params=nothing)
    # get gate unitary
    if !isnothing(params)
        gate_unitary = tomatrix(gate, params)
    else
        gate_unitary = tomatrix(gate)
    end

    # create kroneker product
    if qubits == 1
        return gate_unitary
    end

    # Get rest
    idxs = getqubitsids(gate)
    minidx = minimum(idxs)
    maxidx = qubits - maximum(idxs) - 1
    #println("idxs: $idxs, minidx: $minidx, maxidx: $maxidx")

    return kron(eye(2^maxidx), gate_unitary, eye(2^minidx))
end

"Method return the error of standard gate"
standardGateError(gate::QuantumGate, params=nothing) = 0.0

function standardGateError(gate::UniversalGate, params=nothing)
    # get gate unitary
    if !isnothing(params)
        m = tomatrix(gate, params)
    else
        m = tomatrix(gate)
    end

    # err = 1.0
    # err *= unitary_error(m, Xmatrix)
    # err *= unitary_error(m, Ymatrix)
    # err *= unitary_error(m, Zmatrix)
    # err *= unitary_error(m, Imatrix)
    # err *= unitary_error(m, Hmatrix)
    # err *= unitary_error(m, Smatrix)
    # err *= unitary_error(m, Sdmatrix)
    # err *= unitary_error(m, Tmatrix)
    # err *= unitary_error(m, Tdmatrix)
    #
    # return log(err + 1.0)

    return minimum([unitary_error(m, Xmatrix),
                    unitary_error(m, Ymatrix),
                    unitary_error(m, Zmatrix),
                    unitary_error(m, Imatrix),
                    unitary_error(m, Hmatrix),
                    unitary_error(m, Smatrix),
                    unitary_error(m, Sdmatrix),
                    unitary_error(m, Tmatrix),
                    unitary_error(m, Tdmatrix)])
end

appendparams!(::Any, ::Any) = nothing
appendparams!(param, value::ParameterT) = push!(param, value.value)
appendparams!(::Any, ::ParamT) = nothing
function appendparams!(param, gate::U3)
    appendparams!(param, gate.θ)
    appendparams!(param, gate.ϕ)
    appendparams!(param, gate.λ)
    nothing
end
function appendparams!(param, gate::U)
    push!(param, gate.β)
    push!(param, gate.γ)
    push!(param, gate.δ)
    nothing
end

appendRandParams!(param, ::ParameterT) = append!(param, rand(1) * 2π)
appendRandParams!(::Any, value::ParamT) = nothing
appendRandParams!(::Any, gate) = nothing
function appendRandParams!(param, gate::U3)
    appendRandParams!(param, gate.θ)
    appendRandParams!(param, gate.ϕ)
    appendRandParams!(param, gate.λ)
    nothing
end
appendRandParams!(param, ::U) = append!(param, rand(3) * 2π)

getlength(value::ParameterT) = 1
getlength(value::ParamT) = 0
Base.length(gate) = 0
Base.length(gate::U3) = getlength(gate.θ) + getlength(gate.ϕ) + getlength(gate.λ)
Base.length(gate::U) = 3

################################################################################
const maxerr = 5e-5

"The library of single qubit gates"
const singlequbitsgates = Dict(Xmatrix => X,
                               Ymatrix => Y,
                               Zmatrix => Z,
                               Smatrix => S,
                               Sdmatrix => Sd,
                               Tmatrix => T,
                               Tdmatrix => Td,
                               Hmatrix => H)
# TODO  + I !!!!!

decompose(gate::QuantumGate) = [gate]

"The simplify the gates."
simplify(gate::QuantumGate) = [gate]

function simplify(gate::UniversalGate)
    exp = tomatrix(gate)

    # check for identify
    if unitary_error(Imatrix, exp) < maxerr
        return []
    end

    # single qubits gate
    for (matrix, mgate) in singlequbitsgates
        if unitary_error(matrix, exp) < maxerr
            # println("found: $(mgate(gate.qubit))")
            return [mgate(gate.qubit)]
        end
    end

    # two qubits gate
    for (matrix, mgate) in singlequbitsgates
        for (matrix2, mgate2) in singlequbitsgates
            # println("$mgate -> err: $(unitary_error(matrix, exp)) == $mgate")
            # println(matrix)
            # println(exp)
            if unitary_error(matrix2 * matrix, exp) < maxerr
                # println("found: $(mgate(gate.qubit))")
                return [mgate(gate.qubit), mgate2(gate.qubit)]
            end
        end
    end

    return [gate]
end

inv(gate::X) = gate
inv(gate::H) = gate
inv(gate::Z) = gate
inv(gate::Y) = gate
inv(gate::Rx) = Rx(gate.qubit, neg(gate.θ))
inv(gate::Ry) = Ry(gate.qubit, neg(gate.θ))
inv(gate::Rz) = Rz(gate.qubit, neg(gate.θ))
inv(gate::CX) = gate
inv(gate::T) = Td(gate.qubit)
inv(gate::S) = Sd(gate.qubit)
inv(gate::Td) = T(gate.qubit)
inv(gate::Sd) = S(gate.qubit)
inv(gate::U3) = U3(gate.qubit, neg(gate.θ), neg(gate.λ), neg(gate.ϕ))

function neg(p::ParameterT)
    @assert !isnothing(p.value) "Unable inv unset parameter."
    return ParameterT(-p.value)
end
neg(value::ParamT) = -value

end  # module Gates

ext: jl


text:
using GeometricAlgebra:
	bits_to_indices, indices_to_bits,
	bits_to_linear_index, linear_index_to_bits,
	sign_from_swaps, factor_from_squares

@testset "bits <-> indices" begin
	for bits ∈ [0b000, 0b101, 0b110, 0b1101001, rand(UInt16), rand(UInt64)]
		# expect to fail for UInt128
		@test indices_to_bits(bits_to_indices(bits)) == bits
	end
end

@testset "bits <-> linear index" begin
	cases = [
		(1, 0) => 0b000,
		(1, 1) => 0b001,
		(2, 1) => 0b010,
		(1, 2) => 0b011,
		(3, 1) => 0b100,
		(2, 2) => 0b101,
		(3, 2) => 0b110,
		(1, 3) => 0b111,
		(1, 8) => 0b01111_1111,
		(2, 8) => 0b10111_1111,
	]
	@testset "no. $(i) of grade $k: $bits" for ((i, k), bits) ∈ cases
		@test linear_index_to_bits(i, k) == bits
		@test bits_to_linear_index(bits) == i
	end

	for _ ∈ 1:10, bits ∈ (rand(UInt8), rand(UInt16))
		# warning: bits_to_linear_index() is typically slow for UInt32 or higher
		@test linear_index_to_bits(bits_to_linear_index(bits), grade(bits)) == bits
	end
end

@testset "sign_from_swaps()" begin
	# note that typographically bits are written 'backwards'; v124 == 0b1011
	@test sign_from_swaps(0b01, 0b10) == +1
	@test sign_from_swaps(0b10, 0b01) == -1
	@test sign_from_swaps(0b11, 0b01) == -1
	@test sign_from_swaps(0b111, 0b01) == +1
end

@testset "factor_from_squares()" begin
	@test factor_from_squares((1, 1, 1), 0b000) == 1
	@test factor_from_squares((1, 1, 1), 0b111) == 1
	@test factor_from_squares((-1, 1, 1), 0b001) == -1
	@test factor_from_squares((-1, 1, 1), 0b010) == +1
	@test factor_from_squares((-1, 1, 1), 0b111) == -1
	@test factor_from_squares((-1, -1, -1), 0b101) == +1
	@test factor_from_squares((-1, -1, -1), 0b111) == -1
end


text:
# constructors
export mirror_symmetric_qubit_3povm, asymmetric_qubit_3povm
export trine_qubit_povm, sic_qubit_povm
export planar_symmetric_qubit_povm

export sqrt_povm
export naimark_dilation

"""
    mirror_symmetric_qubit_3povm( θ :: Real ) :: POVM{Float64}

Constructs a `POVM` aligned with the three mirror symmetric qubit states.
The first measurement is aligned with the ``|0\\rangle`` state and the remaining
two are symmetric across the z-axis.

A `DomainError` is thrown if argument `θ` ∉ [π/4,π/2].
"""
function mirror_symmetric_qubit_3povm(θ::Real) :: POVM{Float64}
    if !( π/4 <= θ <= π/2)
        throw(DomainError(θ, "angle must be in range [π/4,π/2]"))
    end

    # POVM weights
    γ1 = 1 - cot(θ)^2
    γ2 = 0.5*(1 + cot(θ)^2)

    mirror_symmetric_qubits = mirror_symmetric_qubit_states(θ)

    π1 = γ1*mirror_symmetric_qubits[1]
    π2 = γ2*mirror_symmetric_qubits[2]
    π3 = γ2*mirror_symmetric_qubits[3]

    POVM([π1, π2, π3])
end

"""
    asymmetric_qubit_3povm( θ1::Real, θ2::Real ) :: POVM

Constructs a general non-orthogonal 3-element `POVM`.

Inputs:
* `θ1 ∈ (0,π/2] or [-π/2,0)`, the angle element 2 makes with ``|0\\rangle`` state.
* `θ2 ∈ [-π/2,0) or (0,π/2]`, the angle element 3 makes with ``|0\\rangle`` state.

A `DomainError` is thrown if `θ1` and `θ2` are not in the valid ranges. Note that
one angle must be positive and the other negative.
"""
function asymmetric_qubit_3povm(θ1::Real,θ2::Real) :: POVM{Float64}

    if !(((0 < θ1 <= π/2) && (-π/2 <= θ2 < 0)) || ((0 < θ2 <= π/2) && (-π/2 <= θ1 < 0)))
        throw(DomainError((θ1,θ2), "angles (θ1, θ2) are not in the valid range"))
    end

    # povm scaling terms derived analytically
    γ2 = 1/(sin(θ1)^2 - sin(θ1)*cos(θ1)*tan(θ2))
    γ3 = -γ2*(sin(θ1)*cos(θ1))/(sin(θ2)*cos(θ2))
    γ1 = 1 - γ2*cos(θ1)^2 - γ3*cos(θ2)^2

    π1 = γ1*[1 0;0 0]
    π2 = γ2*[cos(θ1)^2 cos(θ1)*sin(θ1); cos(θ1)*sin(θ1) sin(θ1)^2]
    π3 = γ3*[cos(θ2)^2 cos(θ2)*sin(θ2); cos(θ2)*sin(θ2) sin(θ2)^2]

    POVM([π1, π2, π3])
end

"""
    planar_symmetric_qubit_povm( n :: Int64 ) :: POVM{Float64}

Constructs an `n`-element `POVM{Float64}` from the [`planar_symmetric_qubit_states`](@ref).
Each state is multipled by a factor of `2/n` to satisfy the completeness relation.

A `DomainError` is thrown if `n ≥ 2` is not satisfied.
"""
function planar_symmetric_qubit_povm(n :: Int64) :: POVM{Float64}
    POVM(2/n*planar_symmetric_qubit_states(n))
end

"""
    trine_qubit_povm() :: POVM{Float64}

The `POVM` with elements parallel to the trine qubit states (see [`trine_qubit_kets`](@ref)).
"""
trine_qubit_povm() :: POVM{Float64} = POVM(2/3 * trine_qubit_states())

"""
    sic_qubit_povm() :: POVM{Complex{Float64}}

The `POVM` with elements parallel to the symmetric informationally complete (SIC) qubit states
(see [`sic_qubit_kets`](@ref)).
"""
sic_qubit_povm() :: POVM{Complex{Float64}} = POVM(0.5 * sic_qubit_states())

"""
    sqrt_povm(priors :: Probabilities, states :: Vector{<:State}) :: POVM

Returns the "pretty good" square-root povm for the given density operator
states and priors.
"""
function sqrt_povm(priors :: Probabilities, states :: Vector{<:State}) :: POVM
    ρ_mix = sum(priors .* states)
    ρ_sqrt = sqrt(inv(ρ_mix))

    Π = map(
        (ρ) -> ρ_sqrt * ρ * ρ_sqrt,
        (priors .* states)
    )

    POVM(Π)
end

"""
    _naimark_kraus_operators(Π::POVM) :: Array{Array{Complex{Float64},2},1}

Returns the Kraus operators for the provided POVM. In general, the Kraus operators
form a continuum and are non-unique. This method applies the construction

```math
k_i = \\sqrt{\\Pi_i}\\otimes |i\\rangle
```
"""
function _naimark_kraus_operators(Π::POVM) :: Array{Array{Complex{Float64},2},1}
    num_el = length(Π)

    map(i -> kron(sqrt(Π[i]), Matrix(I,num_el,num_el)[:,i]), 1:num_el)
end

"""
    naimark_dilation( Π::POVM )

Returns the dilated projectors which are equivalent to the provided POVM. During
measurement, the measured state must be tensored with the ancilla.
"""
function naimark_dilation(Π::POVM)
    n = length(Π)
    d = size(Π[1])[1]

    ancilla = zeros(n,n)
    ancilla[1,1] = 1

    k = _naimark_kraus_operators(Π)
    V = sum(k)

    orth_comp = nullspace(V')

    U_cols = []
    null_id = 1
    for i in 1:n*d
        if floor((i-1)/n) == (i-1)/n
            push!(U_cols,V[:,convert(Int,(i-1)/n + 1)])
        else
            push!(U_cols, orth_comp[:,null_id])
            null_id = null_id + 1
        end

    end

    U = hcat(U_cols...)

    p_ancilla = map( i -> begin
        p = zeros(n,n)
        p[i,i] = 1

        return p
    end, 1:n )

    projectors = map(i -> U'*kron(Matrix(1I,d,d),p_ancilla[i])*U, 1:n)

    Dict(
        "projectors" => projectors,
        "ancilla" => ancilla
    )
end


text:
module Reproduce

export
    create_info!,
    create_info,
    create_custom_info!,
    create_custom_info,
    get_save_dir,
    get_hash,
    get_git_info
include("parse.jl")

export ItemCollection, search, details
include("search.jl")

abstract type AbstractArgIter end

export ArgIterator, ArgLooper
include("args_iter.jl")
include("args_looper.jl")

export Experiment, create_experiment_dir, add_experiment, post_experiment
include("exp_utils.jl")

export job
include("parallel_experiment.jl")
# greeting() = println("Hello World!")

end # module


text:
@testset "4.7.1 (c trig)^m (d trig)^n" begin
    (a, b, c, d, m, n, x, ) = @variables a b c d m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*Trig(a+b*x)^m*Trig(2*a+2*b*x)^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*sin(a+b*x)^m*sin(2*a+2*b*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(a+b*x)^m*sin(2*a+2*b*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^7, x, 4, (128*sin(a + b*x)^9)/(9*b) - (384*sin(a + b*x)^11)/(11*b) + (384*sin(a + b*x)^13)/(13*b) - (128*sin(a + b*x)^15)/(15*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^6, x, 4, -((64*cos(a + b*x)^7)/(7*b)) + (64*cos(a + b*x)^9)/(3*b) - (192*cos(a + b*x)^11)/(11*b) + (64*cos(a + b*x)^13)/(13*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^5, x, 4, (32*sin(a + b*x)^7)/(7*b) - (64*sin(a + b*x)^9)/(9*b) + (32*sin(a + b*x)^11)/(11*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^4, x, 4, -((16*cos(a + b*x)^5)/(5*b)) + (32*cos(a + b*x)^7)/(7*b) - (16*cos(a + b*x)^9)/(9*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^3, x, 4, (8*sin(a + b*x)^5)/(5*b) - (8*sin(a + b*x)^7)/(7*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^2, x, 4, -((4*cos(a + b*x)^3)/(3*b)) + (4*cos(a + b*x)^5)/(5*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^1, x, 1, sin(a + b*x)/(2*b) - sin(3*a + 3*b*x)/(6*b)]
    @test_int [sin(a + b*x)*csc(2*a + 2*b*x)^1, x, 2, atanh(sin(a + b*x))/(2*b)]
    @test_int [sin(a + b*x)*csc(2*a + 2*b*x)^2, x, 4, -(atanh(cos(a + b*x))/(4*b)) + sec(a + b*x)/(4*b)]
    @test_int [sin(a + b*x)*csc(2*a + 2*b*x)^3, x, 5, (3*atanh(sin(a + b*x)))/(16*b) - (3*csc(a + b*x))/(16*b) + (csc(a + b*x)*sec(a + b*x)^2)/(16*b)]
    @test_int [sin(a + b*x)*csc(2*a + 2*b*x)^4, x, 6, -((5*atanh(cos(a + b*x)))/(32*b)) + (5*sec(a + b*x))/(32*b) + (5*sec(a + b*x)^3)/(96*b) - (csc(a + b*x)^2*sec(a + b*x)^3)/(32*b)]
    @test_int [sin(a + b*x)*csc(2*a + 2*b*x)^5, x, 7, (35*atanh(sin(a + b*x)))/(256*b) - (35*csc(a + b*x))/(256*b) - (35*csc(a + b*x)^3)/(768*b) + (7*csc(a + b*x)^3*sec(a + b*x)^2)/(256*b) + (csc(a + b*x)^3*sec(a + b*x)^4)/(128*b)]

    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^5, x, 5, (4*sin(a + b*x)^8)/b - (32*sin(a + b*x)^10)/(5*b) + (8*sin(a + b*x)^12)/(3*b)]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^4, x, 6, (3*x)/16 - (3*cos(2*a + 2*b*x)*sin(2*a + 2*b*x))/(32*b) - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^3)/(16*b) - sin(2*a + 2*b*x)^5/(20*b)]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^3, x, 4, (4*sin(a + b*x)^6)/(3*b) - sin(a + b*x)^8/b]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^2, x, 5, x/4 - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x))/(8*b) - sin(2*a + 2*b*x)^3/(12*b)]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^1, x, 3, sin(a + b*x)^4/(2*b)]
    @test_int [sin(a + b*x)^2*csc(2*a + 2*b*x)^1, x, 2, -(log(cos(a + b*x))/(2*b))]
    @test_int [sin(a + b*x)^2*csc(2*a + 2*b*x)^2, x, 3, tan(a + b*x)/(4*b)]
    @test_int [sin(a + b*x)^2*csc(2*a + 2*b*x)^3, x, 4, log(tan(a + b*x))/(8*b) + tan(a + b*x)^2/(16*b)]
    @test_int [sin(a + b*x)^2*csc(2*a + 2*b*x)^4, x, 4, -(cot(a + b*x)/(16*b)) + tan(a + b*x)/(8*b) + tan(a + b*x)^3/(48*b)]
    @test_int [sin(a + b*x)^2*csc(2*a + 2*b*x)^5, x, 5, -(cot(a + b*x)^2/(64*b)) + (3*log(tan(a + b*x)))/(32*b) + (3*tan(a + b*x)^2)/(64*b) + tan(a + b*x)^4/(128*b)]

    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^5, x, 4, (32*sin(a + b*x)^9)/(9*b) - (64*sin(a + b*x)^11)/(11*b) + (32*sin(a + b*x)^13)/(13*b)]
    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^4, x, 4, -((16*cos(a + b*x)^5)/(5*b)) + (48*cos(a + b*x)^7)/(7*b) - (16*cos(a + b*x)^9)/(3*b) + (16*cos(a + b*x)^11)/(11*b)]
    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^3, x, 4, (8*sin(a + b*x)^7)/(7*b) - (8*sin(a + b*x)^9)/(9*b)]
    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^2, x, 4, -((4*cos(a + b*x)^3)/(3*b)) + (8*cos(a + b*x)^5)/(5*b) - (4*cos(a + b*x)^7)/(7*b)]
    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^1, x, 3, (2*sin(a + b*x)^5)/(5*b)]
    @test_int [sin(a + b*x)^3*csc(2*a + 2*b*x)^1, x, 4, atanh(sin(a + b*x))/(2*b) - sin(a + b*x)/(2*b)]
    @test_int [sin(a + b*x)^3*csc(2*a + 2*b*x)^2, x, 3, sec(a + b*x)/(4*b)]
    @test_int [sin(a + b*x)^3*csc(2*a + 2*b*x)^3, x, 3, atanh(sin(a + b*x))/(16*b) + (sec(a + b*x)*tan(a + b*x))/(16*b)]
    @test_int [sin(a + b*x)^3*csc(2*a + 2*b*x)^4, x, 5, -(atanh(cos(a + b*x))/(16*b)) + sec(a + b*x)/(16*b) + sec(a + b*x)^3/(48*b)]
    @test_int [sin(a + b*x)^3*csc(2*a + 2*b*x)^5, x, 6, (15*atanh(sin(a + b*x)))/(256*b) - (15*csc(a + b*x))/(256*b) + (5*csc(a + b*x)*sec(a + b*x)^2)/(256*b) + (csc(a + b*x)*sec(a + b*x)^4)/(128*b)]

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^8, x, 4, -((256*cos(a + b*x)^9)/(9*b)) + (768*cos(a + b*x)^11)/(11*b) - (768*cos(a + b*x)^13)/(13*b) + (256*cos(a + b*x)^15)/(15*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^7, x, 4, (128*sin(a + b*x)^7)/(7*b) - (128*sin(a + b*x)^9)/(3*b) + (384*sin(a + b*x)^11)/(11*b) - (128*sin(a + b*x)^13)/(13*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^6, x, 4, -((64*cos(a + b*x)^7)/(7*b)) + (128*cos(a + b*x)^9)/(9*b) - (64*cos(a + b*x)^11)/(11*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^5, x, 4, (32*sin(a + b*x)^5)/(5*b) - (64*sin(a + b*x)^7)/(7*b) + (32*sin(a + b*x)^9)/(9*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^4, x, 4, -((16*cos(a + b*x)^5)/(5*b)) + (16*cos(a + b*x)^7)/(7*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^3, x, 4, (8*sin(a + b*x)^3)/(3*b) - (8*sin(a + b*x)^5)/(5*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^2, x, 3, -((4*cos(a + b*x)^3)/(3*b))]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^1, x, 2, (2*sin(a + b*x))/b]
    @test_int [csc(a + b*x)*csc(2*a + 2*b*x)^1, x, 4, atanh(sin(a + b*x))/(2*b) - csc(a + b*x)/(2*b)]
    @test_int [csc(a + b*x)*csc(2*a + 2*b*x)^2, x, 5, -((3*atanh(cos(a + b*x)))/(8*b)) + (3*sec(a + b*x))/(8*b) - (csc(a + b*x)^2*sec(a + b*x))/(8*b)]
    @test_int [csc(a + b*x)*csc(2*a + 2*b*x)^3, x, 6, (5*atanh(sin(a + b*x)))/(16*b) - (5*csc(a + b*x))/(16*b) - (5*csc(a + b*x)^3)/(48*b) + (csc(a + b*x)^3*sec(a + b*x)^2)/(16*b)]
    @test_int [csc(a + b*x)*csc(2*a + 2*b*x)^4, x, 7, -((35*atanh(cos(a + b*x)))/(128*b)) + (35*sec(a + b*x))/(128*b) + (35*sec(a + b*x)^3)/(384*b) - (7*csc(a + b*x)^2*sec(a + b*x)^3)/(128*b) - (csc(a + b*x)^4*sec(a + b*x)^3)/(64*b)]

    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^8, x, 9, (5*x)/8 + (5*cos(a + b*x)*sin(a + b*x))/(8*b) + (5*cos(a + b*x)^3*sin(a + b*x))/(12*b) + (cos(a + b*x)^5*sin(a + b*x))/(3*b) + (2*cos(a + b*x)^7*sin(a + b*x))/(7*b) - (16*cos(a + b*x)^9*sin(a + b*x))/(7*b) - (160*cos(a + b*x)^9*sin(a + b*x)^3)/(21*b) - (128*cos(a + b*x)^9*sin(a + b*x)^5)/(7*b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^7, x, 5, -((16*cos(a + b*x)^8)/b) + (128*cos(a + b*x)^10)/(5*b) - (32*cos(a + b*x)^12)/(3*b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^6, x, 7, (3*x)/4 + (3*cos(a + b*x)*sin(a + b*x))/(4*b) + (cos(a + b*x)^3*sin(a + b*x))/(2*b) + (2*cos(a + b*x)^5*sin(a + b*x))/(5*b) - (12*cos(a + b*x)^7*sin(a + b*x))/(5*b) - (32*cos(a + b*x)^7*sin(a + b*x)^3)/(5*b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^5, x, 4, -((16*cos(a + b*x)^6)/(3*b)) + (4*cos(a + b*x)^8)/b]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^4, x, 5, x + (cos(a + b*x)*sin(a + b*x))/b + (2*cos(a + b*x)^3*sin(a + b*x))/(3*b) - (8*cos(a + b*x)^5*sin(a + b*x))/(3*b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^3, x, 3, -((2*cos(a + b*x)^4)/b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^2, x, 3, 2*x + (2*cos(a + b*x)*sin(a + b*x))/b]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^1, x, 2, (2*log(sin(a + b*x)))/b]
    @test_int [csc(a + b*x)^2*csc(2*a + 2*b*x)^1, x, 4, -(cot(a + b*x)^2/(4*b)) + log(tan(a + b*x))/(2*b)]
    @test_int [csc(a + b*x)^2*csc(2*a + 2*b*x)^2, x, 4, -(cot(a + b*x)/(2*b)) - cot(a + b*x)^3/(12*b) + tan(a + b*x)/(4*b)]
    @test_int [csc(a + b*x)^2*csc(2*a + 2*b*x)^3, x, 5, -((3*cot(a + b*x)^2)/(16*b)) - cot(a + b*x)^4/(32*b) + (3*log(tan(a + b*x)))/(8*b) + tan(a + b*x)^2/(16*b)]
    @test_int [csc(a + b*x)^2*csc(2*a + 2*b*x)^4, x, 4, -((3*cot(a + b*x))/(8*b)) - cot(a + b*x)^3/(12*b) - cot(a + b*x)^5/(80*b) + tan(a + b*x)/(4*b) + tan(a + b*x)^3/(48*b)]
    @test_int [csc(a + b*x)^2*csc(2*a + 2*b*x)^5, x, 5, -((5*cot(a + b*x)^2)/(32*b)) - (5*cot(a + b*x)^4)/(128*b) - cot(a + b*x)^6/(192*b) + (5*log(tan(a + b*x)))/(16*b) + (5*tan(a + b*x)^2)/(64*b) + tan(a + b*x)^4/(128*b)]
    @test_int [csc(a + b*x)^2*csc(2*a + 2*b*x)^6, x, 4, -((5*cot(a + b*x))/(16*b)) - (5*cot(a + b*x)^3)/(64*b) - (3*cot(a + b*x)^5)/(160*b) - cot(a + b*x)^7/(448*b) + (15*tan(a + b*x))/(64*b) + tan(a + b*x)^3/(32*b) + tan(a + b*x)^5/(320*b)]

    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^10, x, 4, -((1024*cos(a + b*x)^11)/(11*b)) + (3072*cos(a + b*x)^13)/(13*b) - (1024*cos(a + b*x)^15)/(5*b) + (1024*cos(a + b*x)^17)/(17*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^9, x, 4, (512*sin(a + b*x)^7)/(7*b) - (2048*sin(a + b*x)^9)/(9*b) + (3072*sin(a + b*x)^11)/(11*b) - (2048*sin(a + b*x)^13)/(13*b) + (512*sin(a + b*x)^15)/(15*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^8, x, 4, -((256*cos(a + b*x)^9)/(9*b)) + (512*cos(a + b*x)^11)/(11*b) - (256*cos(a + b*x)^13)/(13*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^7, x, 4, (128*sin(a + b*x)^5)/(5*b) - (384*sin(a + b*x)^7)/(7*b) + (128*sin(a + b*x)^9)/(3*b) - (128*sin(a + b*x)^11)/(11*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^6, x, 4, -((64*cos(a + b*x)^7)/(7*b)) + (64*cos(a + b*x)^9)/(9*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^5, x, 4, (32*sin(a + b*x)^3)/(3*b) - (64*sin(a + b*x)^5)/(5*b) + (32*sin(a + b*x)^7)/(7*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^4, x, 3, -((16*cos(a + b*x)^5)/(5*b))]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^3, x, 3, (8*sin(a + b*x))/b - (8*sin(a + b*x)^3)/(3*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^2, x, 4, -((4*atanh(cos(a + b*x)))/b) + (4*cos(a + b*x))/b]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^1, x, 3, -((2*csc(a + b*x))/b)]
    @test_int [csc(a + b*x)^3*csc(2*a + 2*b*x)^1, x, 5, atanh(sin(a + b*x))/(2*b) - csc(a + b*x)/(2*b) - csc(a + b*x)^3/(6*b)]
    @test_int [csc(a + b*x)^3*csc(2*a + 2*b*x)^2, x, 6, -((15*atanh(cos(a + b*x)))/(32*b)) + (15*sec(a + b*x))/(32*b) - (5*csc(a + b*x)^2*sec(a + b*x))/(32*b) - (csc(a + b*x)^4*sec(a + b*x))/(16*b)]
    @test_int [csc(a + b*x)^3*csc(2*a + 2*b*x)^3, x, 6, (7*atanh(sin(a + b*x)))/(16*b) - (7*csc(a + b*x))/(16*b) - (7*csc(a + b*x)^3)/(48*b) - (7*csc(a + b*x)^5)/(80*b) + (csc(a + b*x)^5*sec(a + b*x)^2)/(16*b)]
    @test_int [csc(a + b*x)^3*csc(2*a + 2*b*x)^4, x, 8, -((105*atanh(cos(a + b*x)))/(256*b)) + (105*sec(a + b*x))/(256*b) + (35*sec(a + b*x)^3)/(256*b) - (21*csc(a + b*x)^2*sec(a + b*x)^3)/(256*b) - (3*csc(a + b*x)^4*sec(a + b*x)^3)/(128*b) - (csc(a + b*x)^6*sec(a + b*x)^3)/(96*b)]

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

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^(5/2), x, 4, -((5*asin(cos(a + b*x) - sin(a + b*x)))/(32*b)) + (5*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(32*b) - (5*cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(16*b) + (5*sin(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(24*b) - (cos(a + b*x)*sin(2*a + 2*b*x)^(5/2))/(6*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^(3/2), x, 3, -((3*asin(cos(a + b*x) - sin(a + b*x)))/(16*b)) - (3*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(16*b) + (3*sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(8*b) - (cos(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(4*b)]
    @test_int [sin(a + b*x)*sin(2*a + 2*b*x)^(1/2), x, 2, -(asin(cos(a + b*x) - sin(a + b*x))/(4*b)) + log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/(4*b) - (cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(2*b)]
    @test_int [sin(a + b*x)/sin(2*a + 2*b*x)^(1/2), x, 1, -(asin(cos(a + b*x) - sin(a + b*x))/(2*b)) - log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/(2*b)]
    @test_int [sin(a + b*x)/sin(2*a + 2*b*x)^(3/2), x, 1, sin(a + b*x)/(b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [sin(a + b*x)/sin(2*a + 2*b*x)^(5/2), x, 2, sin(a + b*x)/(3*b*sin(2*a + 2*b*x)^(3/2)) - (2*cos(a + b*x))/(3*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [sin(a + b*x)/sin(2*a + 2*b*x)^(7/2), x, 3, sin(a + b*x)/(5*b*sin(2*a + 2*b*x)^(5/2)) - (4*cos(a + b*x))/(15*b*sin(2*a + 2*b*x)^(3/2)) + (8*sin(a + b*x))/(15*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [sin(a + b*x)/sin(2*a + 2*b*x)^(9/2), x, 4, sin(a + b*x)/(7*b*sin(2*a + 2*b*x)^(7/2)) - (6*cos(a + b*x))/(35*b*sin(2*a + 2*b*x)^(5/2)) + (8*sin(a + b*x))/(35*b*sin(2*a + 2*b*x)^(3/2)) - (16*cos(a + b*x))/(35*b*sqrt(sin(2*a + 2*b*x)))]

    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^(7/2), x, 4, (5*Elliptic.F(a - pi/4 + b*x, 2))/(42*b) - (5*cos(2*a + 2*b*x)*sqrt(sin(2*a + 2*b*x)))/(42*b) - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(5/2))/(14*b) - sin(2*a + 2*b*x)^(9/2)/(18*b)]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^(5/2), x, 3, (3*Elliptic.E(a - pi/4 + b*x, 2))/(10*b) - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(3/2))/(10*b) - sin(2*a + 2*b*x)^(7/2)/(14*b)]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^(3/2), x, 3, Elliptic.F(a - pi/4 + b*x, 2)/(6*b) - (cos(2*a + 2*b*x)*sqrt(sin(2*a + 2*b*x)))/(6*b) - sin(2*a + 2*b*x)^(5/2)/(10*b)]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^(1/2), x, 2, Elliptic.E(a - pi/4 + b*x, 2)/(2*b) - sin(2*a + 2*b*x)^(3/2)/(6*b)]
    @test_int [sin(a + b*x)^2/sin(2*a + 2*b*x)^(1/2), x, 2, Elliptic.F(a - pi/4 + b*x, 2)/(2*b) - sqrt(sin(2*a + 2*b*x))/(2*b)]
    @test_int [sin(a + b*x)^2/sin(2*a + 2*b*x)^(3/2), x, 2, -(Elliptic.E(a - pi/4 + b*x, 2)/(2*b)) + sin(a + b*x)^2/(b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [sin(a + b*x)^2/sin(2*a + 2*b*x)^(5/2), x, 2, Elliptic.F(a - pi/4 + b*x, 2)/(6*b) + sin(a + b*x)^2/(3*b*sin(2*a + 2*b*x)^(3/2))]
    @test_int [sin(a + b*x)^2/sin(2*a + 2*b*x)^(7/2), x, 3, -((3*Elliptic.E(a - pi/4 + b*x, 2))/(10*b)) + sin(a + b*x)^2/(5*b*sin(2*a + 2*b*x)^(5/2)) - (3*cos(2*a + 2*b*x))/(10*b*sqrt(sin(2*a + 2*b*x)))]

    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^(3/2), x, 4, -((7*asin(cos(a + b*x) - sin(a + b*x)))/(64*b)) - (7*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(64*b) + (7*sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(32*b) - (7*cos(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(48*b) - (sin(a + b*x)*sin(2*a + 2*b*x)^(5/2))/(12*b)]
    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^(1/2), x, 3, -((5*asin(cos(a + b*x) - sin(a + b*x)))/(32*b)) + (5*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(32*b) - (5*cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(16*b) - (sin(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(8*b)]
    @test_int [sin(a + b*x)^3/sin(2*a + 2*b*x)^(1/2), x, 2, -((3*asin(cos(a + b*x) - sin(a + b*x)))/(8*b)) - (3*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(8*b) - (sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(4*b)]
    @test_int [sin(a + b*x)^3/sin(2*a + 2*b*x)^(3/2), x, 3, asin(cos(a + b*x) - sin(a + b*x))/(4*b) - log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/(4*b) + sin(a + b*x)/(b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [sin(a + b*x)^3/sin(2*a + 2*b*x)^(5/2), x, 1, sin(a + b*x)^3/(3*b*sin(2*a + 2*b*x)^(3/2))]
    @test_int [sin(a + b*x)^3/sin(2*a + 2*b*x)^(7/2), x, 2, sin(a + b*x)^3/(5*b*sin(2*a + 2*b*x)^(5/2)) + sin(a + b*x)/(5*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [sin(a + b*x)^3/sin(2*a + 2*b*x)^(9/2), x, 3, sin(a + b*x)^3/(7*b*sin(2*a + 2*b*x)^(7/2)) + (2*sin(a + b*x))/(21*b*sin(2*a + 2*b*x)^(3/2)) - (4*cos(a + b*x))/(21*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [sin(a + b*x)^3/sin(2*a + 2*b*x)^(11/2), x, 4, sin(a + b*x)^3/(9*b*sin(2*a + 2*b*x)^(9/2)) + sin(a + b*x)/(15*b*sin(2*a + 2*b*x)^(5/2)) - (4*cos(a + b*x))/(45*b*sin(2*a + 2*b*x)^(3/2)) + (8*sin(a + b*x))/(45*b*sqrt(sin(2*a + 2*b*x)))]

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^(7/2), x, 5, -((5*asin(cos(a + b*x) - sin(a + b*x)))/(16*b)) - (5*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(16*b) + (5*sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(8*b) - (5*cos(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(12*b) + (sin(a + b*x)*sin(2*a + 2*b*x)^(5/2))/(3*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^(5/2), x, 4, -((3*asin(cos(a + b*x) - sin(a + b*x)))/(8*b)) + (3*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(8*b) - (3*cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(4*b) + (sin(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(2*b)]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^(3/2), x, 3, -(asin(cos(a + b*x) - sin(a + b*x))/(2*b)) - log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/(2*b) + (sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/b]
    @test_int [csc(a + b*x)*sin(2*a + 2*b*x)^(1/2), x, 2, -(asin(cos(a + b*x) - sin(a + b*x))/b) + log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/b]
    @test_int [csc(a + b*x)/sin(2*a + 2*b*x)^(1/2), x, 1, -((csc(a + b*x)*sqrt(sin(2*a + 2*b*x)))/b)]
    @test_int [csc(a + b*x)/sin(2*a + 2*b*x)^(3/2), x, 3, -((2*cos(a + b*x))/(3*b*sin(2*a + 2*b*x)^(3/2))) + (4*sin(a + b*x))/(3*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [csc(a + b*x)/sin(2*a + 2*b*x)^(5/2), x, 4, -((2*cos(a + b*x))/(5*b*sin(2*a + 2*b*x)^(5/2))) + (8*sin(a + b*x))/(15*b*sin(2*a + 2*b*x)^(3/2)) - (16*cos(a + b*x))/(15*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [csc(a + b*x)/sin(2*a + 2*b*x)^(7/2), x, 5, -((2*cos(a + b*x))/(7*b*sin(2*a + 2*b*x)^(7/2))) + (12*sin(a + b*x))/(35*b*sin(2*a + 2*b*x)^(5/2)) - (16*cos(a + b*x))/(35*b*sin(2*a + 2*b*x)^(3/2)) + (32*sin(a + b*x))/(35*b*sqrt(sin(2*a + 2*b*x)))]

    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^(9/2), x, 4, (6*Elliptic.E(a - pi/4 + b*x, 2))/(5*b) - (2*cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(3/2))/(5*b) - (2*cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(7/2))/(7*b) + (csc(a + b*x)^2*sin(2*a + 2*b*x)^(11/2))/(7*b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^(7/2), x, 4, (2*Elliptic.F(a - pi/4 + b*x, 2))/(3*b) - (2*cos(2*a + 2*b*x)*sqrt(sin(2*a + 2*b*x)))/(3*b) - (2*cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(5/2))/(5*b) + (csc(a + b*x)^2*sin(2*a + 2*b*x)^(9/2))/(5*b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^(5/2), x, 3, (2*Elliptic.E(a - pi/4 + b*x, 2))/b - (2*cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(3/2))/(3*b) + (csc(a + b*x)^2*sin(2*a + 2*b*x)^(7/2))/(3*b)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^(3/2), x, 3, (2*Elliptic.F(a - pi/4 + b*x, 2))/b - (2*cos(2*a + 2*b*x)*sqrt(sin(2*a + 2*b*x)))/b + (csc(a + b*x)^2*sin(2*a + 2*b*x)^(5/2))/b]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^(1/2), x, 2, -((2*Elliptic.E(a - pi/4 + b*x, 2))/b) - (csc(a + b*x)^2*sin(2*a + 2*b*x)^(3/2))/b]
    @test_int [csc(a + b*x)^2/sin(2*a + 2*b*x)^(1/2), x, 2, (2*Elliptic.F(a - pi/4 + b*x, 2))/(3*b) - (csc(a + b*x)^2*sqrt(sin(2*a + 2*b*x)))/(3*b)]
    @test_int [csc(a + b*x)^2/sin(2*a + 2*b*x)^(3/2), x, 3, -((6*Elliptic.E(a - pi/4 + b*x, 2))/(5*b)) - (6*cos(2*a + 2*b*x))/(5*b*sqrt(sin(2*a + 2*b*x))) - csc(a + b*x)^2/(5*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [csc(a + b*x)^2/sin(2*a + 2*b*x)^(5/2), x, 3, (10*Elliptic.F(a - pi/4 + b*x, 2))/(21*b) - (10*cos(2*a + 2*b*x))/(21*b*sin(2*a + 2*b*x)^(3/2)) - csc(a + b*x)^2/(7*b*sin(2*a + 2*b*x)^(3/2))]
    @test_int [csc(a + b*x)^2/sin(2*a + 2*b*x)^(7/2), x, 4, -((14*Elliptic.E(a - pi/4 + b*x, 2))/(15*b)) - (14*cos(2*a + 2*b*x))/(45*b*sin(2*a + 2*b*x)^(5/2)) - csc(a + b*x)^2/(9*b*sin(2*a + 2*b*x)^(5/2)) - (14*cos(2*a + 2*b*x))/(15*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [csc(a + b*x)^2/sin(2*a + 2*b*x)^(9/2), x, 4, (30*Elliptic.F(a - pi/4 + b*x, 2))/(77*b) - (18*cos(2*a + 2*b*x))/(77*b*sin(2*a + 2*b*x)^(7/2)) - csc(a + b*x)^2/(11*b*sin(2*a + 2*b*x)^(7/2)) - (30*cos(2*a + 2*b*x))/(77*b*sin(2*a + 2*b*x)^(3/2))]

    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^(9/2), x, 7, -((7*asin(cos(a + b*x) - sin(a + b*x)))/(8*b)) + (7*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(8*b) - (7*cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(4*b) + (7*sin(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(6*b) - (14*cos(a + b*x)*sin(2*a + 2*b*x)^(5/2))/(15*b) + (4*sin(a + b*x)*sin(2*a + 2*b*x)^(7/2))/(5*b) + (csc(a + b*x)^3*sin(2*a + 2*b*x)^(11/2))/(5*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^(7/2), x, 6, -((5*asin(cos(a + b*x) - sin(a + b*x)))/(4*b)) - (5*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(4*b) + (5*sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(2*b) - (5*cos(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(3*b) + (4*sin(a + b*x)*sin(2*a + 2*b*x)^(5/2))/(3*b) + (csc(a + b*x)^3*sin(2*a + 2*b*x)^(9/2))/(3*b)]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^(5/2), x, 5, -((3*asin(cos(a + b*x) - sin(a + b*x)))/b) + (3*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/b - (6*cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/b + (4*sin(a + b*x)*sin(2*a + 2*b*x)^(3/2))/b + (csc(a + b*x)^3*sin(2*a + 2*b*x)^(7/2))/b]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^(3/2), x, 4, (2*asin(cos(a + b*x) - sin(a + b*x)))/b + (2*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/b - (4*sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/b - (csc(a + b*x)^3*sin(2*a + 2*b*x)^(5/2))/b]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^(1/2), x, 1, -((csc(a + b*x)^3*sin(2*a + 2*b*x)^(3/2))/(3*b))]
    @test_int [csc(a + b*x)^3/sin(2*a + 2*b*x)^(1/2), x, 2, -((4*csc(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(5*b)) - (csc(a + b*x)^3*sqrt(sin(2*a + 2*b*x)))/(5*b)]
    @test_int [csc(a + b*x)^3/sin(2*a + 2*b*x)^(3/2), x, 4, -((16*cos(a + b*x))/(21*b*sin(2*a + 2*b*x)^(3/2))) - csc(a + b*x)^3/(7*b*sqrt(sin(2*a + 2*b*x))) + (32*sin(a + b*x))/(21*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [csc(a + b*x)^3/sin(2*a + 2*b*x)^(5/2), x, 5, -((8*cos(a + b*x))/(15*b*sin(2*a + 2*b*x)^(5/2))) - csc(a + b*x)^3/(9*b*sin(2*a + 2*b*x)^(3/2)) + (32*sin(a + b*x))/(45*b*sin(2*a + 2*b*x)^(3/2)) - (64*cos(a + b*x))/(45*b*sqrt(sin(2*a + 2*b*x)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(a+b*x)^m*sin(2*a+2*b*x)^n*with*m*symbolic=#

    @test_int [sin(a + b*x)^3*sin(2*a + 2*b*x)^m, x, 2, ((cos(a + b*x)^2)^((1 - m)/2)*HypergeometricFunctions._₂F₁((1 - m)/2, (4 + m)/2, (6 + m)/2, sin(a + b*x)^2)*sin(a + b*x)^3*sin(2*a + 2*b*x)^m*tan(a + b*x))/(b*(4 + m))]
    @test_int [sin(a + b*x)^2*sin(2*a + 2*b*x)^m, x, 2, ((cos(a + b*x)^2)^((1 - m)/2)*HypergeometricFunctions._₂F₁((1 - m)/2, (3 + m)/2, (5 + m)/2, sin(a + b*x)^2)*sin(a + b*x)^2*sin(2*a + 2*b*x)^m*tan(a + b*x))/(b*(3 + m))]
    @test_int [sin(a + b*x)^1*sin(2*a + 2*b*x)^m, x, 2, ((cos(a + b*x)^2)^((1 - m)/2)*HypergeometricFunctions._₂F₁((1 - m)/2, (2 + m)/2, (4 + m)/2, sin(a + b*x)^2)*sin(a + b*x)*sin(2*a + 2*b*x)^m*tan(a + b*x))/(b*(2 + m))]
    @test_int [csc(a + b*x)^1*sin(2*a + 2*b*x)^m, x, 2, ((cos(a + b*x)^2)^((1 - m)/2)*HypergeometricFunctions._₂F₁((1 - m)/2, m/2, (2 + m)/2, sin(a + b*x)^2)*sec(a + b*x)*sin(2*a + 2*b*x)^m)/(b*m)]
    @test_int [csc(a + b*x)^2*sin(2*a + 2*b*x)^m, x, 2, -(((cos(a + b*x)^2)^((1 - m)/2)*csc(a + b*x)*HypergeometricFunctions._₂F₁((1 - m)/2, (1/2)*(-1 + m), (1 + m)/2, sin(a + b*x)^2)*sec(a + b*x)*sin(2*a + 2*b*x)^m)/(b*(1 - m)))]
    @test_int [csc(a + b*x)^3*sin(2*a + 2*b*x)^m, x, 2, -(((cos(a + b*x)^2)^((1 - m)/2)*csc(a + b*x)^2*HypergeometricFunctions._₂F₁((1 - m)/2, (1/2)*(-2 + m), m/2, sin(a + b*x)^2)*sec(a + b*x)*sin(2*a + 2*b*x)^m)/(b*(2 - m)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*sin(2*a+2*b*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*sin(2*a+2*b*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^7, x, 4, -((128*cos(a + b*x)^9)/(9*b)) + (384*cos(a + b*x)^11)/(11*b) - (384*cos(a + b*x)^13)/(13*b) + (128*cos(a + b*x)^15)/(15*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^6, x, 4, (64*sin(a + b*x)^7)/(7*b) - (64*sin(a + b*x)^9)/(3*b) + (192*sin(a + b*x)^11)/(11*b) - (64*sin(a + b*x)^13)/(13*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^5, x, 4, -((32*cos(a + b*x)^7)/(7*b)) + (64*cos(a + b*x)^9)/(9*b) - (32*cos(a + b*x)^11)/(11*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^4, x, 4, (16*sin(a + b*x)^5)/(5*b) - (32*sin(a + b*x)^7)/(7*b) + (16*sin(a + b*x)^9)/(9*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^3, x, 4, -((8*cos(a + b*x)^5)/(5*b)) + (8*cos(a + b*x)^7)/(7*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^2, x, 4, (4*sin(a + b*x)^3)/(3*b) - (4*sin(a + b*x)^5)/(5*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^1, x, 1, -(cos(a + b*x)/(2*b)) - cos(3*a + 3*b*x)/(6*b)]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^1, x, 2, -(atanh(cos(a + b*x))/(2*b))]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^2, x, 4, atanh(sin(a + b*x))/(4*b) - csc(a + b*x)/(4*b)]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^3, x, 5, -((3*atanh(cos(a + b*x)))/(16*b)) + (3*sec(a + b*x))/(16*b) - (csc(a + b*x)^2*sec(a + b*x))/(16*b)]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^4, x, 6, (5*atanh(sin(a + b*x)))/(32*b) - (5*csc(a + b*x))/(32*b) - (5*csc(a + b*x)^3)/(96*b) + (csc(a + b*x)^3*sec(a + b*x)^2)/(32*b)]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^5, x, 7, -((35*atanh(cos(a + b*x)))/(256*b)) + (35*sec(a + b*x))/(256*b) + (35*sec(a + b*x)^3)/(768*b) - (7*csc(a + b*x)^2*sec(a + b*x)^3)/(256*b) - (csc(a + b*x)^4*sec(a + b*x)^3)/(128*b)]

    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^5, x, 5, -((4*cos(a + b*x)^8)/b) + (32*cos(a + b*x)^10)/(5*b) - (8*cos(a + b*x)^12)/(3*b)]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^4, x, 6, (3*x)/16 - (3*cos(2*a + 2*b*x)*sin(2*a + 2*b*x))/(32*b) - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^3)/(16*b) + sin(2*a + 2*b*x)^5/(20*b)]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^3, x, 4, -((4*cos(a + b*x)^6)/(3*b)) + cos(a + b*x)^8/b]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^2, x, 5, x/4 - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x))/(8*b) + sin(2*a + 2*b*x)^3/(12*b)]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^1, x, 3, -(cos(a + b*x)^4/(2*b))]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^1, x, 2, log(sin(a + b*x))/(2*b)]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^2, x, 3, -(cot(a + b*x)/(4*b))]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^3, x, 4, -(cot(a + b*x)^2/(16*b)) + log(tan(a + b*x))/(8*b)]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^4, x, 4, -(cot(a + b*x)/(8*b)) - cot(a + b*x)^3/(48*b) + tan(a + b*x)/(16*b)]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^5, x, 5, -((3*cot(a + b*x)^2)/(64*b)) - cot(a + b*x)^4/(128*b) + (3*log(tan(a + b*x)))/(32*b) + tan(a + b*x)^2/(64*b)]

    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^5, x, 4, -((32*cos(a + b*x)^9)/(9*b)) + (64*cos(a + b*x)^11)/(11*b) - (32*cos(a + b*x)^13)/(13*b)]
    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^4, x, 4, (16*sin(a + b*x)^5)/(5*b) - (48*sin(a + b*x)^7)/(7*b) + (16*sin(a + b*x)^9)/(3*b) - (16*sin(a + b*x)^11)/(11*b)]
    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^3, x, 4, -((8*cos(a + b*x)^7)/(7*b)) + (8*cos(a + b*x)^9)/(9*b)]
    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^2, x, 4, (4*sin(a + b*x)^3)/(3*b) - (8*sin(a + b*x)^5)/(5*b) + (4*sin(a + b*x)^7)/(7*b)]
    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^1, x, 3, -((2*cos(a + b*x)^5)/(5*b))]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^1, x, 4, -(atanh(cos(a + b*x))/(2*b)) + cos(a + b*x)/(2*b)]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^2, x, 3, -(csc(a + b*x)/(4*b))]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^3, x, 3, -(atanh(cos(a + b*x))/(16*b)) - (cot(a + b*x)*csc(a + b*x))/(16*b)]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^4, x, 5, atanh(sin(a + b*x))/(16*b) - csc(a + b*x)/(16*b) - csc(a + b*x)^3/(48*b)]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^5, x, 6, -((15*atanh(cos(a + b*x)))/(256*b)) + (15*sec(a + b*x))/(256*b) - (5*csc(a + b*x)^2*sec(a + b*x))/(256*b) - (csc(a + b*x)^4*sec(a + b*x))/(128*b)]

    #= ::Subsubsection:: =#
    #=m<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*sin(2*a+2*b*x)^(n/2)=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^(5/2), x, 4, -((5*asin(cos(a + b*x) - sin(a + b*x)))/(32*b)) - (5*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(32*b) + (5*sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(16*b) - (5*cos(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(24*b) + (sin(a + b*x)*sin(2*a + 2*b*x)^(5/2))/(6*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^(3/2), x, 3, -((3*asin(cos(a + b*x) - sin(a + b*x)))/(16*b)) + (3*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(16*b) - (3*cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(8*b) + (sin(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(4*b)]
    @test_int [cos(a + b*x)*sin(2*a + 2*b*x)^(1/2), x, 2, -(asin(cos(a + b*x) - sin(a + b*x))/(4*b)) - log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/(4*b) + (sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(2*b)]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^(1/2), x, 1, -(asin(cos(a + b*x) - sin(a + b*x))/(2*b)) + log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/(2*b)]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^(3/2), x, 1, -(cos(a + b*x)/(b*sqrt(sin(2*a + 2*b*x))))]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^(5/2), x, 2, -(cos(a + b*x)/(3*b*sin(2*a + 2*b*x)^(3/2))) + (2*sin(a + b*x))/(3*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^(7/2), x, 3, -(cos(a + b*x)/(5*b*sin(2*a + 2*b*x)^(5/2))) + (4*sin(a + b*x))/(15*b*sin(2*a + 2*b*x)^(3/2)) - (8*cos(a + b*x))/(15*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [cos(a + b*x)/sin(2*a + 2*b*x)^(9/2), x, 4, -(cos(a + b*x)/(7*b*sin(2*a + 2*b*x)^(7/2))) + (6*sin(a + b*x))/(35*b*sin(2*a + 2*b*x)^(5/2)) - (8*cos(a + b*x))/(35*b*sin(2*a + 2*b*x)^(3/2)) + (16*sin(a + b*x))/(35*b*sqrt(sin(2*a + 2*b*x)))]

    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^(7/2), x, 4, (5*Elliptic.F(a - pi/4 + b*x, 2))/(42*b) - (5*cos(2*a + 2*b*x)*sqrt(sin(2*a + 2*b*x)))/(42*b) - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(5/2))/(14*b) + sin(2*a + 2*b*x)^(9/2)/(18*b)]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^(5/2), x, 3, (3*Elliptic.E(a - pi/4 + b*x, 2))/(10*b) - (cos(2*a + 2*b*x)*sin(2*a + 2*b*x)^(3/2))/(10*b) + sin(2*a + 2*b*x)^(7/2)/(14*b)]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^(3/2), x, 3, Elliptic.F(a - pi/4 + b*x, 2)/(6*b) - (cos(2*a + 2*b*x)*sqrt(sin(2*a + 2*b*x)))/(6*b) + sin(2*a + 2*b*x)^(5/2)/(10*b)]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^(1/2), x, 2, Elliptic.E(a - pi/4 + b*x, 2)/(2*b) + sin(2*a + 2*b*x)^(3/2)/(6*b)]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^(1/2), x, 2, Elliptic.F(a - pi/4 + b*x, 2)/(2*b) + sqrt(sin(2*a + 2*b*x))/(2*b)]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^(3/2), x, 2, -(Elliptic.E(a - pi/4 + b*x, 2)/(2*b)) - cos(a + b*x)^2/(b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^(5/2), x, 2, Elliptic.F(a - pi/4 + b*x, 2)/(6*b) - cos(a + b*x)^2/(3*b*sin(2*a + 2*b*x)^(3/2))]
    @test_int [cos(a + b*x)^2/sin(2*a + 2*b*x)^(7/2), x, 3, -((3*Elliptic.E(a - pi/4 + b*x, 2))/(10*b)) - cos(a + b*x)^2/(5*b*sin(2*a + 2*b*x)^(5/2)) - (3*cos(2*a + 2*b*x))/(10*b*sqrt(sin(2*a + 2*b*x)))]

    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^(3/2), x, 4, -((7*asin(cos(a + b*x) - sin(a + b*x)))/(64*b)) + (7*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(64*b) - (7*cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(32*b) + (7*sin(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(48*b) + (cos(a + b*x)*sin(2*a + 2*b*x)^(5/2))/(12*b)]
    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^(1/2), x, 3, -((5*asin(cos(a + b*x) - sin(a + b*x)))/(32*b)) - (5*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(32*b) + (5*sin(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(16*b) + (cos(a + b*x)*sin(2*a + 2*b*x)^(3/2))/(8*b)]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^(1/2), x, 2, -((3*asin(cos(a + b*x) - sin(a + b*x)))/(8*b)) + (3*log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x))))/(8*b) + (cos(a + b*x)*sqrt(sin(2*a + 2*b*x)))/(4*b)]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^(3/2), x, 3, asin(cos(a + b*x) - sin(a + b*x))/(4*b) + log(cos(a + b*x) + sin(a + b*x) + sqrt(sin(2*a + 2*b*x)))/(4*b) - cos(a + b*x)/(b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^(5/2), x, 1, -(cos(a + b*x)^3/(3*b*sin(2*a + 2*b*x)^(3/2)))]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^(7/2), x, 2, -(cos(a + b*x)^3/(5*b*sin(2*a + 2*b*x)^(5/2))) - cos(a + b*x)/(5*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^(9/2), x, 3, -(cos(a + b*x)^3/(7*b*sin(2*a + 2*b*x)^(7/2))) - (2*cos(a + b*x))/(21*b*sin(2*a + 2*b*x)^(3/2)) + (4*sin(a + b*x))/(21*b*sqrt(sin(2*a + 2*b*x)))]
    @test_int [cos(a + b*x)^3/sin(2*a + 2*b*x)^(11/2), x, 4, -(cos(a + b*x)^3/(9*b*sin(2*a + 2*b*x)^(9/2))) - cos(a + b*x)/(15*b*sin(2*a + 2*b*x)^(5/2)) + (4*sin(a + b*x))/(45*b*sin(2*a + 2*b*x)^(3/2)) - (8*cos(a + b*x))/(45*b*sqrt(sin(2*a + 2*b*x)))]

    #= 2*cos(x)/sqrt(sin(2*x)) == csc(x)*sqrt(sin(2*x)) =#
    @test_int [cos(x)/sqrt(sin(2*x)), x, 1, (-(1/2))*asin(cos(x) - sin(x)) + (1/2)*log(cos(x) + sin(x) + sqrt(sin(2*x)))]
    @test_int [csc(x)*sqrt(sin(2*x)), x, 2, -asin(cos(x) - sin(x)) + log(cos(x) + sin(x) + sqrt(sin(2*x)))]

    #= ::Subsubsection:: =#
    #=m<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*sin(2*a+2*b*x)^n*with*m*symbolic=#

    @test_int [cos(a + b*x)^3*sin(2*a + 2*b*x)^m, x, 2, -((1/(b*(4 + m)))*(cos(a + b*x)^3*cot(a + b*x)*HypergeometricFunctions._₂F₁((1 - m)/2, (4 + m)/2, (6 + m)/2, cos(a + b*x)^2)*(sin(a + b*x)^2)^((1 - m)/2)*sin(2*a + 2*b*x)^m))]
    @test_int [cos(a + b*x)^2*sin(2*a + 2*b*x)^m, x, 2, -((1/(b*(3 + m)))*(cos(a + b*x)^2*cot(a + b*x)*HypergeometricFunctions._₂F₁((1 - m)/2, (3 + m)/2, (5 + m)/2, cos(a + b*x)^2)*(sin(a + b*x)^2)^((1 - m)/2)*sin(2*a + 2*b*x)^m))]
    @test_int [cos(a + b*x)^1*sin(2*a + 2*b*x)^m, x, 2, -((1/(b*(2 + m)))*(cos(a + b*x)*cot(a + b*x)*HypergeometricFunctions._₂F₁((1 - m)/2, (2 + m)/2, (4 + m)/2, cos(a + b*x)^2)*(sin(a + b*x)^2)^((1 - m)/2)*sin(2*a + 2*b*x)^m))]

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

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

    @test_int [cos(a + b*x)^2*sin(a + b*x)^3*sin(2*a + 2*b*x)^2, x, 4, -((4*cos(a + b*x)^5)/(5*b)) + (8*cos(a + b*x)^7)/(7*b) - (4*cos(a + b*x)^9)/(9*b)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*sin(a+b*x)^n*sin(2*a+2*b*x)^(p/2)=#

    #= ::Title:: =#
    #=Integrands*of*the*form*Trig(a+b*x)^m*Trig(c+d*x)^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*sin(a+b*x)^m*Trig(c+d*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(a+b*x)^m*sin(c+d*x)^n=#

    @test_int [sin(a + b*x)*sin(c + d*x)^n, x, 10, -((2^(-1 - n)*ℯ^(im*(a - c*n) + im*(b - d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, (b - d*n)/(2*d), (1/2)*(2 + b/d - n), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(b - d*n))) - (2^(-1 - n)*ℯ^((-im)*(a + c*n) - im*(b + d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -((b + d*n)/(2*d)), 1 - (b + d*n)/(2*d), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(b + d*n))]

    @test_int [sin(a + b*x)*sin(c + d*x)^3, x, 6, -(sin(a - 3*c + (b - 3*d)*x)/(8*(b - 3*d))) + (3*sin(a - c + (b - d)*x))/(8*(b - d)) - (3*sin(a + c + (b + d)*x))/(8*(b + d)) + sin(a + 3*c + (b + 3*d)*x)/(8*(b + 3*d))]
    @test_int [sin(a + b*x)*sin(c + d*x)^2, x, 5, -(cos(a + b*x)/(2*b)) + cos(a - 2*c + (b - 2*d)*x)/(4*(b - 2*d)) + cos(a + 2*c + (b + 2*d)*x)/(4*(b + 2*d))]
    @test_int [sin(a + b*x)*sin(c + d*x)^1, x, 4, sin(a - c + (b - d)*x)/(2*(b - d)) - sin(a + c + (b + d)*x)/(2*(b + d))]
    @test_int [sin(a + b*x)*csc(c + b*x)^1, x, 3, x*cos(a - c) + (log(sin(c + b*x))*sin(a - c))/b]
    @test_int [sin(a + b*x)*csc(c + b*x)^2, x, 4, -((atanh(cos(c + b*x))*cos(a - c))/b) - (csc(c + b*x)*sin(a - c))/b]
    @test_int [sin(a + b*x)*csc(c + b*x)^3, x, 5, -((cos(a - c)*cot(c + b*x))/b) - (csc(c + b*x)^2*sin(a - c))/(2*b)]
    @test_int [sin(a + b*x)*csc(c + b*x)^4, x, 5, -((atanh(cos(c + b*x))*cos(a - c))/(2*b)) - (cos(a - c)*cot(c + b*x)*csc(c + b*x))/(2*b) - (csc(c + b*x)^3*sin(a - c))/(3*b)]
    @test_int [sin(a + b*x)*csc(c + b*x)^5, x, 5, -((cos(a - c)*cot(c + b*x))/b) - (cos(a - c)*cot(c + b*x)^3)/(3*b) - (csc(c + b*x)^4*sin(a - c))/(4*b)]
    @test_int [sin(a + b*x)*csc(c + b*x)^6, x, 6, -((3*atanh(cos(c + b*x))*cos(a - c))/(8*b)) - (3*cos(a - c)*cot(c + b*x)*csc(c + b*x))/(8*b) - (cos(a - c)*cot(c + b*x)*csc(c + b*x)^3)/(4*b) - (csc(c + b*x)^5*sin(a - c))/(5*b)]

    @test_int [sin(a + b*x)^2*sin(c + d*x)^n, x, 15, -((im*2^(-2 - n)*ℯ^((-im)*(2*a + c*n) - im*(2*b + d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁((1/2)*(-((2*b)/d) - n), -n, (1/2)*(2 - (2*b)/d - n), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(2*b + d*n))) + (im*2^(-2 - n)*ℯ^(im*(2*a - c*n) + im*(2*b - d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁((1/2)*((2*b)/d - n), -n, (1/2)*(2 + (2*b)/d - n), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(2*b - d*n)) + (im*2^(-1 - n)*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -(n/2), 1 - n/2, ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*(c + d*x)))^n*(d*n))]

    @test_int [sin(a + b*x)^2*sin(c + d*x)^1, x, 5, -(cos(2*a - c + (2*b - d)*x)/(4*(2*b - d))) - cos(c + d*x)/(2*d) + cos(2*a + c + (2*b + d)*x)/(4*(2*b + d))]
    @test_int [sin(a + b*x)^2*sin(c + d*x)^2, x, 6, x/4 - sin(2*a + 2*b*x)/(8*b) + sin(2*(a - c) + 2*(b - d)*x)/(16*(b - d)) - sin(2*c + 2*d*x)/(8*d) + sin(2*(a + c) + 2*(b + d)*x)/(16*(b + d))]
    @test_int [sin(a + b*x)^2*sin(c + d*x)^3, x, 8, cos(2*a - 3*c + (2*b - 3*d)*x)/(16*(2*b - 3*d)) - (3*cos(2*a - c + (2*b - d)*x))/(16*(2*b - d)) - (3*cos(c + d*x))/(8*d) + cos(3*c + 3*d*x)/(24*d) + (3*cos(2*a + c + (2*b + d)*x))/(16*(2*b + d)) - cos(2*a + 3*c + (2*b + 3*d)*x)/(16*(2*b + 3*d))]

    @test_int [sin(a + b*x)^3*sin(c + d*x)^n, x, 18, (2^(-3 - n)*ℯ^(im*(3*a - c*n) + im*(3*b - d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁((1/2)*((3*b)/d - n), -n, (1/2)*(2 + (3*b)/d - n), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(3*b - d*n)) - (3*2^(-3 - n)*ℯ^(im*(a - c*n) + im*(b - d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, (b - d*n)/(2*d), (1/2)*(2 + b/d - n), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(b - d*n)) - (3*2^(-3 - n)*ℯ^((-im)*(a + c*n) - im*(b + d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -((b + d*n)/(2*d)), 1 - (b + d*n)/(2*d), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(b + d*n)) + (2^(-3 - n)*ℯ^((-im)*(3*a + c*n) - im*(3*b + d*n)*x + im*n*(c + d*x))*(im/ℯ^(im*(c + d*x)) - im*ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -((3*b + d*n)/(2*d)), (1/2)*(2 - (3*b)/d - n), ℯ^(2*im*(c + d*x))))/((1 - ℯ^(2*im*c + 2*im*d*x))^n*(3*b + d*n))]

    @test_int [sin(a + b*x)^3*sin(c + d*x)^1, x, 6, (3*sin(a - c + (b - d)*x))/(8*(b - d)) - sin(3*a - c + (3*b - d)*x)/(8*(3*b - d)) - (3*sin(a + c + (b + d)*x))/(8*(b + d)) + sin(3*a + c + (3*b + d)*x)/(8*(3*b + d))]
    @test_int [sin(a + b*x)^3*sin(c + d*x)^2, x, 8, -((3*cos(a + b*x))/(8*b)) + cos(3*a + 3*b*x)/(24*b) + (3*cos(a - 2*c + (b - 2*d)*x))/(16*(b - 2*d)) - cos(3*a - 2*c + (3*b - 2*d)*x)/(16*(3*b - 2*d)) + (3*cos(a + 2*c + (b + 2*d)*x))/(16*(b + 2*d)) - cos(3*a + 2*c + (3*b + 2*d)*x)/(16*(3*b + 2*d))]
    @test_int [sin(a + b*x)^3*sin(c + d*x)^3, x, 10, -((3*sin(a - 3*c + (b - 3*d)*x))/(32*(b - 3*d))) + (9*sin(a - c + (b - d)*x))/(32*(b - d)) + sin(3*(a - c) + 3*(b - d)*x)/(96*(b - d)) - (3*sin(3*a - c + (3*b - d)*x))/(32*(3*b - d)) - (9*sin(a + c + (b + d)*x))/(32*(b + d)) - sin(3*(a + c) + 3*(b + d)*x)/(96*(b + d)) + (3*sin(3*a + c + (3*b + d)*x))/(32*(3*b + d)) + (3*sin(a + 3*c + (b + 3*d)*x))/(32*(b + 3*d))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(a+b*x)^m*cos(c+d*x)^n=#

    @test_int [sin(a + b*x)*cos(c + d*x)^n, x, 8, -((2^(-1 - n)*ℯ^(im*(a - c*n) + im*(b - d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, (b - d*n)/(2*d), (1/2)*(2 + b/d - n), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(b - d*n))) - (2^(-1 - n)*ℯ^((-im)*(a + c*n) - im*(b + d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -((b + d*n)/(2*d)), 1 - (b + d*n)/(2*d), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(b + d*n))]

    @test_int [sin(a + b*x)*cos(c + d*x)^3, x, 6, -(cos(a - 3*c + (b - 3*d)*x)/(8*(b - 3*d))) - (3*cos(a - c + (b - d)*x))/(8*(b - d)) - (3*cos(a + c + (b + d)*x))/(8*(b + d)) - cos(a + 3*c + (b + 3*d)*x)/(8*(b + 3*d))]
    @test_int [sin(a + b*x)*cos(c + d*x)^2, x, 5, -(cos(a + b*x)/(2*b)) - cos(a - 2*c + (b - 2*d)*x)/(4*(b - 2*d)) - cos(a + 2*c + (b + 2*d)*x)/(4*(b + 2*d))]
    @test_int [sin(a + b*x)*cos(c + d*x)^1, x, 4, -(cos(a - c + (b - d)*x)/(2*(b - d))) - cos(a + c + (b + d)*x)/(2*(b + d))]
    @test_int [sin(a + b*x)*sec(c + b*x)^1, x, 3, -((cos(a - c)*log(cos(c + b*x)))/b) + x*sin(a - c)]
    @test_int [sin(a + b*x)*sec(c + b*x)^2, x, 4, (cos(a - c)*sec(c + b*x))/b + (atanh(sin(c + b*x))*sin(a - c))/b]
    @test_int [sin(a + b*x)*sec(c + b*x)^3, x, 5, (cos(a - c)*sec(c + b*x)^2)/(2*b) + (sin(a - c)*tan(c + b*x))/b]
    @test_int [sin(a + b*x)*sec(c + b*x)^4, x, 5, (cos(a - c)*sec(c + b*x)^3)/(3*b) + (atanh(sin(c + b*x))*sin(a - c))/(2*b) + (sec(c + b*x)*sin(a - c)*tan(c + b*x))/(2*b)]
    @test_int [sin(a + b*x)*sec(c + b*x)^5, x, 5, (cos(a - c)*sec(c + b*x)^4)/(4*b) + (sin(a - c)*tan(c + b*x))/b + (sin(a - c)*tan(c + b*x)^3)/(3*b)]
    @test_int [sin(a + b*x)*sec(c + b*x)^6, x, 6, (cos(a - c)*sec(c + b*x)^5)/(5*b) + (3*atanh(sin(c + b*x))*sin(a - c))/(8*b) + (3*sec(c + b*x)*sin(a - c)*tan(c + b*x))/(8*b) + (sec(c + b*x)^3*sin(a - c)*tan(c + b*x))/(4*b)]

    @test_int [sin(a + b*x)^2*cos(c + d*x)^n, x, 11, -((im*2^(-2 - n)*ℯ^((-im)*(2*a + c*n) - im*(2*b + d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁((1/2)*(-((2*b)/d) - n), -n, (1/2)*(2 - (2*b)/d - n), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(2*b + d*n))) + (im*2^(-2 - n)*ℯ^(im*(2*a - c*n) + im*(2*b - d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁((1/2)*((2*b)/d - n), -n, (1/2)*(2 + (2*b)/d - n), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(2*b - d*n)) + (im*2^(-1 - n)*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -(n/2), 1 - n/2, -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*(c + d*x)))^n*(d*n))]

    @test_int [sin(a + b*x)^2*cos(c + d*x)^1, x, 5, -(sin(2*a - c + (2*b - d)*x)/(4*(2*b - d))) + sin(c + d*x)/(2*d) - sin(2*a + c + (2*b + d)*x)/(4*(2*b + d))]
    @test_int [sin(a + b*x)^2*cos(c + d*x)^2, x, 6, x/4 - sin(2*a + 2*b*x)/(8*b) - sin(2*(a - c) + 2*(b - d)*x)/(16*(b - d)) + sin(2*c + 2*d*x)/(8*d) - sin(2*(a + c) + 2*(b + d)*x)/(16*(b + d))]
    @test_int [sin(a + b*x)^2*cos(c + d*x)^3, x, 8, -(sin(2*a - 3*c + (2*b - 3*d)*x)/(16*(2*b - 3*d))) - (3*sin(2*a - c + (2*b - d)*x))/(16*(2*b - d)) + (3*sin(c + d*x))/(8*d) + sin(3*c + 3*d*x)/(24*d) - (3*sin(2*a + c + (2*b + d)*x))/(16*(2*b + d)) - sin(2*a + 3*c + (2*b + 3*d)*x)/(16*(2*b + 3*d))]

    @test_int [sin(a + b*x)^3*cos(c + d*x)^n, x, 14, (2^(-3 - n)*ℯ^(im*(3*a - c*n) + im*(3*b - d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁((1/2)*((3*b)/d - n), -n, (1/2)*(2 + (3*b)/d - n), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(3*b - d*n)) - (3*2^(-3 - n)*ℯ^(im*(a - c*n) + im*(b - d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, (b - d*n)/(2*d), (1/2)*(2 + b/d - n), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(b - d*n)) - (3*2^(-3 - n)*ℯ^((-im)*(a + c*n) - im*(b + d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -((b + d*n)/(2*d)), 1 - (b + d*n)/(2*d), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(b + d*n)) + (2^(-3 - n)*ℯ^((-im)*(3*a + c*n) - im*(3*b + d*n)*x + im*n*(c + d*x))*(ℯ^((-im)*(c + d*x)) + ℯ^(im*(c + d*x)))^n*HypergeometricFunctions._₂F₁(-n, -((3*b + d*n)/(2*d)), (1/2)*(2 - (3*b)/d - n), -ℯ^(2*im*(c + d*x))))/((1 + ℯ^(2*im*c + 2*im*d*x))^n*(3*b + d*n))]

    @test_int [sin(a + b*x)^3*cos(c + d*x)^1, x, 6, -((3*cos(a - c + (b - d)*x))/(8*(b - d))) + cos(3*a - c + (3*b - d)*x)/(8*(3*b - d)) - (3*cos(a + c + (b + d)*x))/(8*(b + d)) + cos(3*a + c + (3*b + d)*x)/(8*(3*b + d))]
    @test_int [sin(a + b*x)^3*cos(c + d*x)^2, x, 8, -((3*cos(a + b*x))/(8*b)) + cos(3*a + 3*b*x)/(24*b) - (3*cos(a - 2*c + (b - 2*d)*x))/(16*(b - 2*d)) + cos(3*a - 2*c + (3*b - 2*d)*x)/(16*(3*b - 2*d)) - (3*cos(a + 2*c + (b + 2*d)*x))/(16*(b + 2*d)) + cos(3*a + 2*c + (3*b + 2*d)*x)/(16*(3*b + 2*d))]
    @test_int [sin(a + b*x)^3*cos(c + d*x)^3, x, 10, -((3*cos(a - 3*c + (b - 3*d)*x))/(32*(b - 3*d))) - (9*cos(a - c + (b - d)*x))/(32*(b - d)) + cos(3*(a - c) + 3*(b - d)*x)/(96*(b - d)) + (3*cos(3*a - c + (3*b - d)*x))/(32*(3*b - d)) - (9*cos(a + c + (b + d)*x))/(32*(b + d)) + cos(3*(a + c) + 3*(b + d)*x)/(96*(b + d)) + (3*cos(3*a + c + (3*b + d)*x))/(32*(3*b + d)) - (3*cos(a + 3*c + (b + 3*d)*x))/(32*(b + 3*d))]

    @test_int [cos(a + b*x)/sin(c + b*x)^1, x, 3, (cos(a - c)*log(sin(c + b*x)))/b - x*sin(a - c)]
    @test_int [cos(a + b*x)/sin(c + b*x)^2, x, 4, -((cos(a - c)*csc(c + b*x))/b) + (atanh(cos(c + b*x))*sin(a - c))/b]
    @test_int [cos(a + b*x)/sin(c + b*x)^3, x, 5, -((cos(a - c)*csc(c + b*x)^2)/(2*b)) + (cot(c + b*x)*sin(a - c))/b]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*sin(a+b*x)^m*tan(c+d*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [sin(a + b*x)*tan(c + b*x)^3, x, 9, -((3*atanh(sin(c + b*x))*cos(a - c))/(2*b)) + (sec(c + b*x)*sin(a - c))/b + sin(a + b*x)/b + (cos(a - c)*sec(c + b*x)*tan(c + b*x))/(2*b)]
    @test_int [sin(a + b*x)*tan(c + b*x)^2, x, 6, cos(a + b*x)/b + (cos(a - c)*sec(c + b*x))/b + (atanh(sin(c + b*x))*sin(a - c))/b]
    @test_int [sin(a + b*x)*tan(c + b*x)^1, x, 3, (atanh(sin(c + b*x))*cos(a - c))/b - sin(a + b*x)/b]
    @test_int [sin(a + b*x)*cot(c + b*x)^1, x, 3, -((atanh(cos(c + b*x))*sin(a - c))/b) + sin(a + b*x)/b]
    @test_int [sin(a + b*x)*cot(c + b*x)^2, x, 6, -((atanh(cos(c + b*x))*cos(a - c))/b) + cos(a + b*x)/b - (csc(c + b*x)*sin(a - c))/b]
    @test_int [sin(a + b*x)*cot(c + b*x)^3, x, 9, -((cos(a - c)*csc(c + b*x))/b) + (3*atanh(cos(c + b*x))*sin(a - c))/(2*b) - (cot(c + b*x)*csc(c + b*x)*sin(a - c))/(2*b) - sin(a + b*x)/b]

    @test_int [sin(a + b*x)*tan(c + d*x), x, 6, im/(ℯ^(im*(a + b*x))*(2*b)) + (im*ℯ^(im*(a + b*x)))/(2*b) - (im*HypergeometricFunctions._₂F₁(1, -(b/(2*d)), 1 - b/(2*d), -ℯ^(2*im*(c + d*x))))/(ℯ^(im*(a + b*x))*b) - (im*ℯ^(im*(a + b*x))*HypergeometricFunctions._₂F₁(1, b/(2*d), 1 + b/(2*d), -ℯ^(2*im*(c + d*x))))/b]
    @test_int [sin(a + b*x)*cot(c + d*x), x, 6, -(im/(ℯ^(im*(a + b*x))*(2*b))) - (im*ℯ^(im*(a + b*x)))/(2*b) + (im*HypergeometricFunctions._₂F₁(1, -(b/(2*d)), 1 - b/(2*d), ℯ^(2*im*(c + d*x))))/(ℯ^(im*(a + b*x))*b) + (im*ℯ^(im*(a + b*x))*HypergeometricFunctions._₂F₁(1, b/(2*d), 1 + b/(2*d), ℯ^(2*im*(c + d*x))))/b]

    #= ::Subsubsection:: =#
    #=m<0=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*Trig(c+d*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*cos(c+d*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [cos(a + b*x)*cos(c + d*x)^3, x, 6, sin(a - 3*c + (b - 3*d)*x)/(8*(b - 3*d)) + (3*sin(a - c + (b - d)*x))/(8*(b - d)) + (3*sin(a + c + (b + d)*x))/(8*(b + d)) + sin(a + 3*c + (b + 3*d)*x)/(8*(b + 3*d))]
    @test_int [cos(a + b*x)*cos(c + d*x)^2, x, 5, sin(a + b*x)/(2*b) + sin(a - 2*c + (b - 2*d)*x)/(4*(b - 2*d)) + sin(a + 2*c + (b + 2*d)*x)/(4*(b + 2*d))]
    @test_int [cos(a + b*x)*cos(c + d*x)^1, x, 4, sin(a - c + (b - d)*x)/(2*(b - d)) + sin(a + c + (b + d)*x)/(2*(b + d))]
    @test_int [cos(a + b*x)*sec(c + b*x)^1, x, 3, x*cos(a - c) + (log(cos(c + b*x))*sin(a - c))/b]
    @test_int [cos(a + b*x)*sec(c + b*x)^2, x, 4, (atanh(sin(c + b*x))*cos(a - c))/b - (sec(c + b*x)*sin(a - c))/b]
    @test_int [cos(a + b*x)*sec(c + b*x)^3, x, 5, -((sec(c + b*x)^2*sin(a - c))/(2*b)) + (cos(a - c)*tan(c + b*x))/b]

    @test_int [cos(a + b*x)^2*cos(c + d*x)^3, x, 8, sin(2*a - 3*c + (2*b - 3*d)*x)/(16*(2*b - 3*d)) + (3*sin(2*a - c + (2*b - d)*x))/(16*(2*b - d)) + (3*sin(c + d*x))/(8*d) + sin(3*c + 3*d*x)/(24*d) + (3*sin(2*a + c + (2*b + d)*x))/(16*(2*b + d)) + sin(2*a + 3*c + (2*b + 3*d)*x)/(16*(2*b + 3*d))]
    @test_int [cos(a + b*x)^2*cos(c + d*x)^2, x, 6, x/4 + sin(2*a + 2*b*x)/(8*b) + sin(2*(a - c) + 2*(b - d)*x)/(16*(b - d)) + sin(2*c + 2*d*x)/(8*d) + sin(2*(a + c) + 2*(b + d)*x)/(16*(b + d))]

    @test_int [cos(a + b*x)^3*cos(c + d*x)^3, x, 10, (3*sin(a - 3*c + (b - 3*d)*x))/(32*(b - 3*d)) + (9*sin(a - c + (b - d)*x))/(32*(b - d)) + sin(3*(a - c) + 3*(b - d)*x)/(96*(b - d)) + (3*sin(3*a - c + (3*b - d)*x))/(32*(3*b - d)) + (9*sin(a + c + (b + d)*x))/(32*(b + d)) + sin(3*(a + c) + 3*(b + d)*x)/(96*(b + d)) + (3*sin(3*a + c + (3*b + d)*x))/(32*(3*b + d)) + (3*sin(a + 3*c + (b + 3*d)*x))/(32*(b + 3*d))]

    #= ::Subsubsection:: =#
    #=m<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*tan(c+d*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [cos(a + b*x)*tan(c + b*x)^3, x, 9, cos(a + b*x)/b + (cos(a - c)*sec(c + b*x))/b + (3*atanh(sin(c + b*x))*sin(a - c))/(2*b) - (sec(c + b*x)*sin(a - c)*tan(c + b*x))/(2*b)]
    @test_int [cos(a + b*x)*tan(c + b*x)^2, x, 6, (atanh(sin(c + b*x))*cos(a - c))/b - (sec(c + b*x)*sin(a - c))/b - sin(a + b*x)/b]
    @test_int [cos(a + b*x)*tan(c + b*x)^1, x, 3, -(cos(a + b*x)/b) - (atanh(sin(c + b*x))*sin(a - c))/b]
    @test_int [cos(a + b*x)*cot(c + b*x)^1, x, 3, -((atanh(cos(c + b*x))*cos(a - c))/b) + cos(a + b*x)/b]
    @test_int [cos(a + b*x)*cot(c + b*x)^2, x, 6, -((cos(a - c)*csc(c + b*x))/b) + (atanh(cos(c + b*x))*sin(a - c))/b - sin(a + b*x)/b]
    @test_int [cos(a + b*x)*cot(c + b*x)^3, x, 9, (3*atanh(cos(c + b*x))*cos(a - c))/(2*b) - cos(a + b*x)/b - (cos(a - c)*cot(c + b*x)*csc(c + b*x))/(2*b) + (csc(c + b*x)*sin(a - c))/b]

    @test_int [cos(a + b*x)*tan(c + d*x), x, 6, 1/(ℯ^(im*(a + b*x))*(2*b)) - ℯ^(im*(a + b*x))/(2*b) - HypergeometricFunctions._₂F₁(1, -(b/(2*d)), 1 - b/(2*d), -ℯ^(2*im*(c + d*x)))/(ℯ^(im*(a + b*x))*b) + (ℯ^(im*(a + b*x))*HypergeometricFunctions._₂F₁(1, b/(2*d), 1 + b/(2*d), -ℯ^(2*im*(c + d*x))))/b]
    @test_int [cos(a + b*x)*cot(c + d*x), x, 6, -(1/(ℯ^(im*(a + b*x))*(2*b))) + ℯ^(im*(a + b*x))/(2*b) + HypergeometricFunctions._₂F₁(1, -(b/(2*d)), 1 - b/(2*d), ℯ^(2*im*(c + d*x)))/(ℯ^(im*(a + b*x))*b) - (ℯ^(im*(a + b*x))*HypergeometricFunctions._₂F₁(1, b/(2*d), 1 + b/(2*d), ℯ^(2*im*(c + d*x))))/b]

    #= ::Subsubsection:: =#
    #=m<0=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*cos(a+b*x)^m*cot(c+d*x)^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tan(a+b*x)^m*Trig(c+d*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tan(a+b*x)^m*tan(c+d*x)^n=#

    #= ::Subsubsection:: =#
    #=m>0=#

    #= ::Subsubsection:: =#
    #=m<0=#
    end

max_stars_repo_path: test/4 Trig functions/4.7 Miscellaneous/4.7.1 (c trig)^m (d trig)^n.jl