text:
using Distributions, Bijectors
using ForwardDiff
using Tracker
using Tracker: TrackedReal, TrackedArray, track, @grad, data

import Base: inv, ∘

import Random: AbstractRNG
import Distributions: logpdf, rand, rand!, _rand!, _logpdf, params

#######################################
# AD stuff "extracted" from Turing.jl #
#######################################

abstract type ADBackend end
struct ForwardDiffAD <: ADBackend end
struct TrackerAD <: ADBackend end

const ADBACKEND = Ref(:forward_diff)
function setadbackend(backend_sym)
    @assert backend_sym == :forward_diff || backend_sym == :reverse_diff
    backend_sym == :forward_diff
    ADBACKEND[] = backend_sym
end

ADBackend() = ADBackend(ADBACKEND[])
ADBackend(T::Symbol) = ADBackend(Val(T))
function ADBackend(::Val{T}) where {T}
    if T === :forward_diff
        return ForwardDiffAD
    else
        return TrackerAD
    end
end

######################
# Bijector interface #
######################

"Abstract type for a bijector."
abstract type AbstractBijector end

"Abstract type of bijectors with fixed dimensionality."
abstract type Bijector{N} <:AbstractBijector end

dimension(b::Bijector{N}) where {N} = N

Broadcast.broadcastable(b::Bijector) = Ref(b)

"""
    isclosedform(b::Bijector)::bool
    isclosedform(b⁻¹::Inversed{<:Bijector})::bool

Returns `true` or `false` depending on whether or not evaluation of `b`
has a closed-form implementation.

Most bijectors have closed-form evaluations, but there are cases where
this is not the case. For example the *inverse* evaluation of `PlanarLayer`
requires an iterative procedure to evaluate and thus is not differentiable.
"""
isclosedform(b::Bijector) = true

"""
Abstract type for a `Bijector{N}` making use of auto-differentation (AD) to
implement `jacobian` and, by impliciation, `logabsdetjac`.
"""
abstract type ADBijector{AD, N} <: Bijector{N} end

"""
    inv(b::Bijector)
    Inversed(b::Bijector)

A `Bijector` representing the inverse transform of `b`.
"""
# TODO: can we do something like `Bijector{N}` instead?
struct Inversed{B <: Bijector, N} <: Bijector{N}
    orig::B

    Inversed(b::B) where {N, B<:Bijector{N}} = new{B, N}(b)
end

inv(b::Bijector) = Inversed(b)
inv(ib::Inversed{<:Bijector}) = ib.orig

"""
    logabsdetjac(b::Bijector, x)
    logabsdetjac(ib::Inversed{<:Bijector}, y)

Computes the log(abs(det(J(b(x))))) where J is the jacobian of the transform.
Similarily for the inverse-transform.

Default implementation for `Inversed{<:Bijector}` is implemented as
`- logabsdetjac` of original `Bijector`.
"""
logabsdetjac(ib::Inversed{<:Bijector}, y) = - logabsdetjac(ib.orig, ib(y))

"""
    forward(b::Bijector, x)
    forward(ib::Inversed{<:Bijector}, y)

Computes both `transform` and `logabsdetjac` in one forward pass, and
returns a named tuple `(rv=b(x), logabsdetjac=logabsdetjac(b, x))`.

This defaults to the call above, but often one can re-use computation
in the computation of the forward pass and the computation of the
`logabsdetjac`. `forward` allows the user to take advantange of such
efficiencies, if they exist.
"""
forward(b::Bijector, x) = (rv=b(x), logabsdetjac=logabsdetjac(b, x))
forward(ib::Inversed{<:Bijector}, y) = (
    rv=ib(y),
    logabsdetjac=logabsdetjac(ib, y)
)

# AD implementations
function jacobian(b::ADBijector{<:ForwardDiffAD}, x::Real)
    return ForwardDiff.derivative(b, x)
end
function jacobian(b::Inversed{<:ADBijector{<:ForwardDiffAD}}, y::Real)
    return ForwardDiff.derivative(b, y)
end
function jacobian(b::ADBijector{<:ForwardDiffAD}, x::AbstractVector{<:Real})
    return ForwardDiff.jacobian(b, x)
end
function jacobian(b::Inversed{<:ADBijector{<:ForwardDiffAD}}, y::AbstractVector{<:Real})
    return ForwardDiff.jacobian(b, y)
end

function jacobian(b::ADBijector{<:TrackerAD}, x::Real)
    return Tracker.data(Tracker.gradient(b, x)[1])
end
function jacobian(b::Inversed{<:ADBijector{<:TrackerAD}}, y::Real)
    return Tracker.data(Tracker.gradient(b, y)[1])
end
function jacobian(b::ADBijector{<:TrackerAD}, x::AbstractVector{<:Real})
    # We extract `data` so that we don't returne a `Tracked` type
    return Tracker.data(Tracker.jacobian(b, x))
end
function jacobian(b::Inversed{<:ADBijector{<:TrackerAD}}, y::AbstractVector{<:Real})
    # We extract `data` so that we don't returne a `Tracked` type
    return Tracker.data(Tracker.jacobian(b, y))
end

struct SingularJacobianException{B} <: Exception where {B<:Bijector}
    b::B
end
function Base.showerror(io::IO, e::SingularJacobianException)
    print(io, "jacobian of $(e.b) is singular")
end

# TODO: allow batch-computation, especially for univariate case?
"Computes the absolute determinant of the Jacobian of the inverse-transformation."
function logabsdetjac(b::ADBijector, x::Real)
    res = log(abs(jacobian(b, x)))
    return isfinite(res) ? res : throw(SingularJacobianException(b))
end

function logabsdetjac(b::ADBijector, x::AbstractVector{<:Real})
    fact = lu(jacobian(b, x), check=false)
    return issuccess(fact) ? log(abs(det(fact))) : throw(SingularJacobianException(b))
end

"""
    logabsdetjacinv(b::Bijector, y)

Just an alias for `logabsdetjac(inv(b), y)`.
"""
logabsdetjacinv(b::Bijector, y) = logabsdetjac(inv(b), y)

###############
# Composition #
###############

"""
    Composed(ts::A)

    ∘(b1::Bijector{N}, b2::Bijector{N})::Composed{<:Tuple}
    composel(ts::Bijector{N}...)::Composed{<:Tuple}
    composer(ts::Bijector{N}...)::Composed{<:Tuple}

where `A` refers to either
- `Tuple{Vararg{<:Bijector{N}}}`: a tuple of bijectors of dimensionality `N`
- `AbstractArray{<:Bijector{N}}`: an array of bijectors of dimensionality `N`

A `Bijector` representing composition of bijectors. `composel` and `composer` results in a
`Composed` for which application occurs from left-to-right and right-to-left, respectively.

Note that all the alternative ways of constructing a `Composed` returns a `Tuple` of bijectors.
This ensures type-stability of implementations of all relating methdos, e.g. `inv`.

If you want to use an `Array` as the container instead you can do

    Composed([b1, b2, ...])

In general this is not advised since you lose type-stability, but there might be cases
where this is desired, e.g. if you have a insanely large number of bijectors to compose.

# Examples
It's important to note that `∘` does what is expected mathematically, which means that the
bijectors are applied to the input right-to-left, e.g. first applying `b2` and then `b1`:
```
(b1 ∘ b2)(x) == b1(b2(x))     # => true
```
But in the `Composed` struct itself, we store the bijectors left-to-right, so that
```
cb1 = b1 ∘ b2                  # => Composed.ts == (b2, b1)
cb2 = composel(b2, b1)         # => Composed.ts == (b2, b1)
cb1(x) == cb2(x) == b1(b2(x))  # => true
```
"""
struct Composed{A, N} <: Bijector{N}
    ts::A

    Composed(bs::C) where {N, C<:Tuple{Vararg{<:Bijector{N}}}} = new{C, N}(bs)
    Composed(bs::A) where {N, A<:AbstractArray{<:Bijector{N}}} = new{A, N}(bs)
end

isclosedform(b::Composed) = all(isclosedform.(b.ts))

"""
    composel(ts::Bijector...)::Composed{<:Tuple}

Constructs `Composed` such that `ts` are applied left-to-right.
"""
composel(ts::Bijector{N}...) where {N} = Composed(ts)

"""
    composer(ts::Bijector...)::Composed{<:Tuple}

Constructs `Composed` such that `ts` are applied right-to-left.
"""
composer(ts::Bijector{N}...) where {N} = Composed(reverse(ts))

# The transformation of `Composed` applies functions left-to-right
# but in mathematics we usually go from right-to-left; this reversal ensures that
# when we use the mathematical composition ∘ we get the expected behavior.
# TODO: change behavior of `transform` of `Composed`?
@generated function ∘(b1::Bijector{N1}, b2::Bijector{N2}) where {N1, N2}
    if N1 == N2
        return :(composel(b2, b1))
    else
        return :(throw(DimensionMismatch("$(typeof(b1)) expects $(N1)-dim but $(typeof(b2)) expects $(N2)-dim")))
    end
end

inv(ct::Composed) = composer(map(inv, ct.ts)...)

# # TODO: should arrays also be using recursive implementation instead?
function (cb::Composed{<:AbstractArray{<:Bijector}})(x)
    res = x
    for b ∈ cb.ts
        res = b(res)
    end

    return res
end

# recursive implementation like this allows type-inference
_transform(x, b1::Bijector, b2::Bijector) = b2(b1(x))
_transform(x, b::Bijector, bs::Bijector...) = _transform(b(x), bs...)
(cb::Composed{<:Tuple})(x) = _transform(x, cb.ts...)

function _logabsdetjac(x, b1::Bijector, b2::Bijector)
    res = forward(b1, x)
    return logabsdetjac(b2, res.rv) + res.logabsdetjac
end
function _logabsdetjac(x, b1::Bijector, bs::Bijector...)
    res = forward(b1, x)
    return _logabsdetjac(res.rv, bs...) + res.logabsdetjac
end
logabsdetjac(cb::Composed, x) = _logabsdetjac(x, cb.ts...)

# Recursive implementation of `forward`
# NOTE: we need this one in the case where `length(cb.ts) == 2`
# in which case forward(...) immediately calls `_forward(::NamedTuple, b::Bijector)`
function _forward(f::NamedTuple, b::Bijector)
    y, logjac = forward(b, f.rv)
    return (rv=y, logabsdetjac=logjac + f.logabsdetjac)
end
function _forward(f::NamedTuple, b1::Bijector, b2::Bijector)
    f1 = forward(b1, f.rv)
    f2 = forward(b2, f1.rv)
    return (rv=f2.rv, logabsdetjac=f2.logabsdetjac + f1.logabsdetjac + f.logabsdetjac)
end
function _forward(f::NamedTuple, b::Bijector, bs::Bijector...)
    f1 = forward(b, f.rv)
    f_ = (rv=f1.rv, logabsdetjac=f1.logabsdetjac + f.logabsdetjac)
    return _forward(f_, bs...)
end
_forward(x, b::Bijector, bs::Bijector...) = _forward(forward(b, x), bs...)
forward(cb::Composed{<:Tuple}, x) = _forward(x, cb.ts...)

function forward(cb::Composed, x)
    rv, logjac = forward(cb.ts[1], x)

    for t in cb.ts[2:end]
        res = forward(t, rv)
        rv = res.rv
        logjac = res.logabsdetjac + logjac
    end
    return (rv=rv, logabsdetjac=logjac)
end

-###########
-# Stacked #
-###########
const ZeroOrOneDimBijector = Union{Bijector{0}, Bijector{1}}

"""
    Stacked(bs)
    Stacked(bs, ranges)
    stack(bs::Bijector{0}...) # where `0` means 0-dim `Bijector`

A `Bijector` which stacks bijectors together which can then be applied to a vector
where `bs[i]::Bijector` is applied to `x[ranges[i]]::UnitRange{Int}`.

# Arguments
- `bs` can be either a `Tuple` or an `AbstractArray` of 0- and/or 1-dimensional bijectors
  - If `bs` is a `Tuple`, implementations are type-stable using generated functions
  - If `bs` is an `AbstractArray`, implementations are _not_ type-stable and use iterative methods
- `ranges` needs to be an iterable consisting of `UnitRange{Int}`
  - `length(bs) == length(ranges)` needs to be true.

# Examples
```
b1 = Logit(0.0, 1.0)
b2 = Identity{0}()
b = stack(b1, b2)
b([0.0, 1.0]) == [b1(0.0), 1.0]  # => true
```
"""
struct Stacked{Bs, N} <: Bijector{1} where N
    bs::Bs
    ranges::NTuple{N, UnitRange{Int}}

    function Stacked(
        bs::C,
        ranges::NTuple{N, UnitRange{Int}}
    ) where {N, C<:Tuple{Vararg{<:ZeroOrOneDimBijector, N}}}
        return new{C, N}(bs, ranges)
    end

    function Stacked(
        bs::A,
        ranges::NTuple{N, UnitRange{Int}}
    ) where {N, A<:AbstractArray{<:Bijector}}
        @assert length(bs) == N "number of bijectors is not same as number of ranges"
        @assert all(isa.(bs, ZeroOrOneDimBijector))
        return new{A, N}(bs, ranges)
    end
end
Stacked(bs, ranges::AbstractArray) = Stacked(bs, tuple(ranges...))
Stacked(bs) = Stacked(bs, tuple([i:i for i = 1:length(bs)]...))

isclosedform(b::Stacked) = all(isclosedform.(b.bs))

stack(bs::Bijector{0}...) = Stacked(bs)

inv(sb::Stacked) = Stacked(inv.(sb.bs), sb.ranges)

# TODO: Is there a better approach to this?
@generated function _transform(x, rs::NTuple{N, UnitRange{Int}}, bs::Bijector...) where N
    exprs = []
    for i = 1:N
        push!(exprs, :(bs[$i](x[rs[$i]])))
    end

    return :(vcat($(exprs...)))
end
_transform(x, rs::NTuple{1, UnitRange{Int}}, b::Bijector) = b(x)

function (sb::Stacked{<:Tuple})(x::AbstractVector{<:Real})
    y = _transform(x, sb.ranges, sb.bs...)
    @assert size(y) == size(x) "x is size $(size(x)) but y is $(size(y))"
    return y
end
function (sb::Stacked{<:AbstractArray, N})(x::AbstractVector{<:Real}) where {N}
    y = vcat([sb.bs[i](x[sb.ranges[i]]) for i = 1:N]...)
    @assert size(y) == size(x) "x is size $(size(x)) but y is $(size(y))"
    return y
end

(sb::Stacked)(x::AbstractMatrix{<:Real}) = mapslices(z -> sb(z), x; dims = 1)

# TODO: implement custom adjoint since we can exploit block-diagonal nature of `Stacked`
function (sb::Stacked)(x::TrackedArray{A, 2}) where {A}
    return Tracker.collect(hcat([sb(x[:, i]) for i = 1:size(x, 2)]...))
end

@generated function logabsdetjac(
    b::Stacked{<:Tuple, N},
    x::AbstractVector{<:Real}
) where {N}
    exprs = []
    for i = 1:N
        push!(exprs, :(sum(logabsdetjac(b.bs[$i], x[b.ranges[$i]]))))
    end

    return :(sum([$(exprs...), ]))
end
function logabsdetjac(
    b::Stacked{<:AbstractArray, N},
    x::AbstractVector{<:Real}
) where {N}
    # TODO: drop the `sum` when we have dimensionality
    return sum([sum(logabsdetjac(b.bs[i], x[b.ranges[i]])) for i = 1:N])
end
function logabsdetjac(b::Stacked, x::AbstractMatrix{<: Real})
    return vec(mapslices(z -> logabsdetjac(b, z), x; dims = 1))
end
function logabsdetjac(b::Stacked, x::TrackedArray{A, 2}) where {A}
    return Tracker.collect(vec(mapslices(z -> logabsdetjac(b, z), x; dims = 1)))
end

# Generates something similar to:
#
# quote
#     (y_1, _logjac) = forward(b.bs[1], x[b.ranges[1]])
#     logjac = sum(_logjac)
#     (y_2, _logjac) = forward(b.bs[2], x[b.ranges[2]])
#     logjac += sum(_logjac)
#     return (rv = vcat(y_1, y_2), logabsdetjac = logjac)
# end
@generated function forward(b::Stacked{T, N}, x::AbstractVector) where {N, T<:Tuple}
    expr = Expr(:block)
    y_names = []

    push!(expr.args, :((y_1, _logjac) = forward(b.bs[1], x[b.ranges[1]])))
    # TODO: drop the `sum` when we have dimensionality
    push!(expr.args, :(logjac = sum(_logjac)))
    push!(y_names, :y_1)
    for i = 2:length(T.parameters)
        y_name = Symbol("y_$i")
        push!(expr.args, :(($y_name, _logjac) = forward(b.bs[$i], x[b.ranges[$i]])))

        # TODO: drop the `sum` when we have dimensionality
        push!(expr.args, :(logjac += sum(_logjac)))

        push!(y_names, y_name)
    end

    push!(expr.args, :(return (rv = vcat($(y_names...)), logabsdetjac = logjac)))
    return expr
end

function forward(sb::Stacked{<:AbstractArray, N}, x::AbstractVector) where {N}
    ys = []
    logjacs = []
    for i = 1:N
        y, logjac = forward(sb.bs[i], x[sb.ranges[i]])
        push!(ys, y)
        # TODO: drop the `sum` when we have dimensionality
        push!(logjacs, sum(logjac))
    end

    return (rv = vcat(ys...), logabsdetjac = sum(logjacs))
end

##############################
# Example bijector: Identity #
##############################

struct Identity{N} <: Bijector{N} end
(::Identity)(x) = x
inv(b::Identity) = b

logabsdetjac(::Identity, x::Real) = zero(eltype(x))
@generated function logabsdetjac(
    b::Identity{N1},
    x::AbstractArray{T2, N2}
) where {N1, T2, N2}
    if N1 == N2
        return :(zero(eltype(x)))
    elseif N1 + 1 == N2
        return :(zeros(eltype(x), size(x, $N2)))
    else
        return :(throw(MethodError(logabsdetjac, (b, x))))
    end
end

###############################
# Example: Logit and Logistic #
###############################
using StatsFuns: logit, logistic

struct Logit{T<:Real} <: Bijector{0}
    a::T
    b::T
end

(b::Logit)(x) = @. logit((x - b.a) / (b.b - b.a))
(ib::Inversed{<:Logit{<:Real}})(y) = @. (ib.orig.b - ib.orig.a) * logistic(y) + ib.orig.a

logabsdetjac(b::Logit{<:Real}, x) = @. - log((x - b.a) * (b.b - x) / (b.b - b.a))

#############
# Exp & Log #
#############

struct Exp{N} <: Bijector{N} end
struct Log{N} <: Bijector{N} end

Exp() = Exp{0}()
Log() = Log{0}()

(b::Exp)(y) = @. exp(y)
(b::Log)(x) = @. log(x)

inv(b::Exp{N}) where {N} = Log{N}()
inv(b::Log{N}) where {N} = Exp{N}()

logabsdetjac(b::Exp{0}, x::Real) = x
logabsdetjac(b::Exp{0}, x::AbstractVector) = x
logabsdetjac(b::Exp{1}, x::AbstractVector) = sum(x)
logabsdetjac(b::Exp{1}, x::AbstractMatrix) = vec(sum(x; dims = 1))

logabsdetjac(b::Log{0}, x::Real) = -log(x)
logabsdetjac(b::Log{0}, x::AbstractVector) = -log.(x)
logabsdetjac(b::Log{1}, x::AbstractVector) = - sum(log.(x))
logabsdetjac(b::Log{1}, x::AbstractMatrix) = - vec(sum(log.(x); dims = 1))

#################
# Shift & Scale #
#################
struct Shift{T, N} <: Bijector{N}
    a::T
end

Shift(a::T; dim::Type{Val{D}} = Val{0}) where {T<:Real, D} = Shift{T, D}(a)
Shift(a::A; dim::Type{Val{D}} = Val{N}) where {T, D, N, A<:AbstractArray{T, N}} = Shift{A, N}(a)

(b::Shift)(x) = b.a + x
(b::Shift{<:Real})(x::AbstractArray) = b.a .+ x
(b::Shift{<:AbstractVector})(x::AbstractMatrix) = b.a .+ x

inv(b::Shift{T, N}) where {T, N} = Shift{T, N}(-b.a)

# FIXME: implement custom adjoint to ensure we don't get tracking
logabsdetjac(b::Shift{T, N}, x) where {T, N} = _logabsdetjac_shift(b.a, x, Val{N})

_logabsdetjac_shift(a::Real, x::Real, ::Type{Val{0}}) = zero(eltype(x))
_logabsdetjac_shift(a::Real, x::AbstractVector{T}, ::Type{Val{0}}) where {T<:Real} = zeros(T, length(x))
_logabsdetjac_shift(a::T1, x::AbstractVector{T2}, ::Type{Val{1}}) where {T1<:Union{Real, AbstractVector}, T2<:Real} = zero(T2)
_logabsdetjac_shift(a::T1, x::AbstractMatrix{T2}, ::Type{Val{1}}) where {T1<:Union{Real, AbstractVector}, T2<:Real} = zeros(T2, size(x, 2))

function _logabsdetjac_shift(a::TrackedReal, x::Real, ::Type{Val{0}})
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val{0}))
end
function _logabsdetjac_shift(a::TrackedReal, x::AbstractVector{T}, ::Type{Val{0}}) where {T<:Real}
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val{0}))
end
function _logabsdetjac_shift(a::T1, x::AbstractVector{T2}, ::Type{Val{1}}) where {T1<:Union{TrackedReal, TrackedVector}, T2<:Real}
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val{1}))
end
function _logabsdetjac_shift(a::T1, x::AbstractMatrix{T2}, ::Type{Val{1}}) where {T1<:Union{TrackedReal, TrackedVector}, T2<:Real}
    return Tracker.param(_logabsdetjac_shift(data(a), data(x), Val{1}))
end

struct Scale{T, N} <: Bijector{N}
    a::T
end

Scale(a::T; dim::Type{Val{D}} = Val{0}) where {T<:Real, D} = Scale{T, D}(a)
Scale(a::A; dim::Type{Val{D}} = Val{N}) where {T, D, N, A<:AbstractArray{T, N}} = Scale{A, D}(a)

(b::Scale)(x) = b.a .* x
(b::Scale{<:Real})(x::AbstractArray) = b.a .* x
(b::Scale{<:AbstractMatrix})(x::AbstractArray) = b.a * x
(b::Scale{<:AbstractVector{<:Real}, 2})(x::AbstractMatrix{<:Real}) = b.a .* x

inv(b::Scale{T, D}) where {T, D} = Scale(inv(b.a); dim = Val{D})
inv(b::Scale{<:AbstractVector, D}) where {D} = Scale(inv.(b.a); dim = Val{D})

# We're going to implement custom adjoint for this
logabsdetjac(b::Scale{T, N}, x) where {T, N} = _logabsdetjac_scale(b.a, x, Val{N})

_logabsdetjac_scale(a::Real, x::Real, ::Type{Val{0}}) = log(abs(a))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Type{Val{0}}) = fill(log(abs(a)), length(x))
_logabsdetjac_scale(a::Real, x::AbstractVector, ::Type{Val{1}}) = log(abs(a)) * length(x)
_logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Type{Val{1}}) = fill(log(abs(a)) * size(x, 1), size(x, 2))
_logabsdetjac_scale(a::AbstractVector, x::AbstractVector, ::Type{Val{1}}) = sum(log.(abs.(a)))
_logabsdetjac_scale(a::AbstractVector, x::AbstractMatrix, ::Type{Val{1}}) = fill(sum(log.(abs.(a))), size(x, 2))
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractVector, ::Type{Val{1}}) = log(abs(det(a)))
_logabsdetjac_scale(a::AbstractMatrix, x::AbstractMatrix{T}, ::Type{Val{1}}) where {T} = log(abs(det(a))) * ones(T, size(x, 2))

# Adjoints for 0-dim and 1-dim `Scale` using `Real`
function _logabsdetjac_scale(a::TrackedReal, x::Real, ::Type{Val{0}})
    return track(_logabsdetjac_scale, a, data(x), Val{0})
end
@grad function _logabsdetjac_scale(a::Real, x::Real, ::Type{Val{0}})
    return _logabsdetjac_scale(data(a), data(x), Val{0}), Δ -> (inv(data(a)) .* Δ, nothing, nothing)
end

# Need to treat `AbstractVector` and `AbstractMatrix` separately due to ambiguity errors
function _logabsdetjac_scale(a::TrackedReal, x::AbstractVector, ::Type{Val{0}})
    return track(_logabsdetjac_scale, a, data(x), Val{0})
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractVector, ::Type{Val{0}})
    da = data(a)
    J = fill(inv.(da), length(x))
    return _logabsdetjac_scale(da, data(x), Val{0}), Δ -> (transpose(J) * Δ, nothing, nothing)
end

function _logabsdetjac_scale(a::TrackedReal, x::AbstractMatrix, ::Type{Val{0}})
    return track(_logabsdetjac_scale, a, data(x), Val{0})
end
@grad function _logabsdetjac_scale(a::Real, x::AbstractMatrix, ::Type{Val{0}})
    da = data(a)
    J = fill(size(x, 1) / da, size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val{0}), Δ -> (transpose(J) * Δ, nothing, nothing)
end

# adjoints for 1-dim and 2-dim `Scale` using `AbstractVector`
function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Type{Val{1}})
    return track(_logabsdetjac_scale, a, data(x), Val{1})
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractVector, ::Type{Val{1}})
    da = data(a)
    J = sum(inv.(da))
    return _logabsdetjac_scale(da, data(x), Val{1}), Δ -> (transpose(J) * Δ, nothing, nothing)
end

function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Type{Val{1}})
    return track(_logabsdetjac_scale, a, data(x), Val{1})
end
@grad function _logabsdetjac_scale(a::TrackedVector, x::AbstractMatrix, ::Type{Val{1}})
    da = data(a)
    J = sum(inv.(da)) .* ones(size(x, 2))
    return _logabsdetjac_scale(da, data(x), Val{1}), Δ -> (transpose(J) * Δ, nothing, nothing)
end

# TODO: implement analytical gradient for scaling a vector using a matrix
# function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Type{Val{1}})
#     track(_logabsdetjac_scale, a, data(x), Val{1})
# end
# @grad function _logabsdetjac_scale(a::TrackedMatrix, x::AbstractVector, ::Type{Val{1}})
#     throw
# end

####################
# Simplex bijector #
####################
struct SimplexBijector{T} <: Bijector{1} where {T} end
SimplexBijector(proj::Bool) = SimplexBijector{Val{proj}}()
SimplexBijector() = SimplexBijector(true)

# The following implementations are basically just copy-paste from `invlink` and
# `link` for `SimplexDistributions` but dropping the dependence on the `Distribution`.
function _clamp(x::T, b::Union{SimplexBijector, Inversed{<:SimplexBijector}}) where {T}
    bounds = (zero(T), one(T))
    clamped_x = clamp(x, bounds...)
    DEBUG && @debug "x = $x, bounds = $bounds, clamped_x = $clamped_x"
    return clamped_x
end

function (b::SimplexBijector{Val{proj}})(x::AbstractVector{T}) where {T, proj}
    y, K = similar(x), length(x)
    @assert K > 1 "x needs to be of length greater than 1"

    ϵ = _eps(T)
    sum_tmp = zero(T)
    @inbounds z = x[1] * (one(T) - 2ϵ) + ϵ # z ∈ [ϵ, 1-ϵ]
    @inbounds y[1] = StatsFuns.logit(z) + log(T(K - 1))
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k - 1]
        # z ∈ [ϵ, 1-ϵ]
        # x[k] = 0 && sum_tmp = 1 -> z ≈ 1
        z = (x[k] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
        y[k] = StatsFuns.logit(z) + log(T(K - k))
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if proj
        y[K] = zero(T)
    else
        y[K] = one(T) - sum_tmp - x[K]
    end

    return y
end

# Vectorised implementation of the above.
function (b::SimplexBijector{Val{proj}})(X::AbstractMatrix{T}) where {T<:Real, proj}
    Y, K, N = similar(X), size(X, 1), size(X, 2)
    @assert K > 1 "x needs to be of length greater than 1"

    ϵ = _eps(T)
    @inbounds @simd for n in 1:size(X, 2)
        sum_tmp = zero(T)
        z = X[1, n] * (one(T) - 2ϵ) + ϵ
        Y[1, n] = StatsFuns.logit(z) + log(T(K - 1))
        for k in 2:(K - 1)
            sum_tmp += X[k - 1, n]
            z = (X[k, n] + ϵ)*(one(T) - 2ϵ)/((one(T) + ϵ) - sum_tmp)
            Y[k, n] = StatsFuns.logit(z) + log(T(K - k))
        end
        sum_tmp += X[K-1, n]
        if proj
            Y[K, n] = zero(T)
        else
            Y[K, n] = one(T) - sum_tmp - X[K, n]
        end
    end

    return Y
end

function (ib::Inversed{<:SimplexBijector{Val{proj}}})(y::AbstractVector{T}) where {T, proj}
    x, K = similar(y), length(y)
    @assert K > 1 "x needs to be of length greater than 1"

    ϵ = _eps(T)
    @inbounds z = StatsFuns.logistic(y[1] - log(T(K - 1)))
    @inbounds x[1] = _clamp((z - ϵ) / (one(T) - 2ϵ), ib.orig)
    sum_tmp = zero(T)
    @inbounds @simd for k = 2:(K - 1)
        z = StatsFuns.logistic(y[k] - log(T(K - k)))
        sum_tmp += x[k-1]
        x[k] = _clamp(((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ, ib.orig)
    end
    @inbounds sum_tmp += x[K - 1]
    @inbounds if proj
        x[K] = _clamp(one(T) - sum_tmp, ib.orig)
    else
        x[K] = _clamp(one(T) - sum_tmp - y[K], ib.orig)
    end

    return x
end

# Vectorised implementation of the above.
function (ib::Inversed{<:SimplexBijector{Val{proj}}})(
    Y::AbstractMatrix{T}
) where {T<:Real, proj}
    X, K, N = similar(Y), size(Y, 1), size(Y, 2)
    @assert K > 1 "x needs to be of length greater than 1"

    ϵ = _eps(T)
    @inbounds @simd for n in 1:size(X, 2)
        sum_tmp, z = zero(T), StatsFuns.logistic(Y[1, n] - log(T(K - 1)))
        X[1, n] = _clamp((z - ϵ) / (one(T) - 2ϵ), ib.orig)
        for k in 2:(K - 1)
            z = StatsFuns.logistic(Y[k, n] - log(T(K - k)))
            sum_tmp += X[k - 1]
            X[k, n] = _clamp(((one(T) + ϵ) - sum_tmp) / (one(T) - 2ϵ) * z - ϵ, ib.orig)
        end
        sum_tmp += X[K - 1, n]
        if proj
            X[K, n] = _clamp(one(T) - sum_tmp, ib.orig)
        else
            X[K, n] = _clamp(one(T) - sum_tmp - Y[K, n], ib.orig)
        end
    end

    return X
end

function logabsdetjac(b::SimplexBijector, x::AbstractVector{T}) where T
    ϵ = _eps(T)
    lp = zero(T)

    K = length(x)

    sum_tmp = zero(eltype(x))
    @inbounds z = x[1]
    lp += log(z + ϵ) + log((one(T) + ϵ) - z)
    @inbounds @simd for k in 2:(K - 1)
        sum_tmp += x[k-1]
        z = x[k] / ((one(T) + ϵ) - sum_tmp)
        lp += log(z + ϵ) + log((one(T) + ϵ) - z) + log((one(T) + ϵ) - sum_tmp)
    end

    return - lp
end

function logabsdetjac(b::SimplexBijector, x::AbstractMatrix{<:Real})
    return vec(mapslices(z -> logabsdetjac(b, z), x; dims = 1))
end

#######################################################
# Constrained to unconstrained distribution bijectors #
#######################################################
"""
    DistributionBijector(d::Distribution)
    DistributionBijector{<:ADBackend, D}(d::Distribution)

This is the default `Bijector` for a distribution.

It uses `link` and `invlink` to compute the transformations, and `AD` to compute
the `jacobian` and `logabsdetjac`.
"""
struct DistributionBijector{AD, D, N} <: ADBijector{AD, N} where {D<:Distribution}
    dist::D
end
function DistributionBijector(dist::D) where {D<:UnivariateDistribution}
    DistributionBijector{ADBackend(), D, 0}(dist)
end
function DistributionBijector(dist::D) where {D<:MultivariateDistribution}
    DistributionBijector{ADBackend(), D, 1}(dist)
end
function DistributionBijector(dist::D) where {D<:MatrixDistribution}
    DistributionBijector{ADBackend(), D, 2}(dist)
end

# Simply uses `link` and `invlink` as transforms with AD to get jacobian
(b::DistributionBijector)(x) = link(b.dist, x)
(ib::Inversed{<:DistributionBijector})(y) = invlink(ib.orig.dist, y)

# Transformed distributions
struct TransformedDistribution{D, B, V} <: Distribution{V, Continuous} where {D<:Distribution{V, Continuous}, B<:Bijector}
    dist::D
    transform::B

    TransformedDistribution(d::UnivariateDistribution, b::Bijector{0}) = new{typeof(d), typeof(b), Univariate}(d, b)
    TransformedDistribution(d::MultivariateDistribution, b::Bijector{1}) = new{typeof(d), typeof(b), Multivariate}(d, b)
    TransformedDistribution(d::MatrixDistribution, b::Bijector{2}) = new{typeof(d), typeof(b), Matrixvariate}(d, b)
end

const UnivariateTransformed = TransformedDistribution{<:Distribution, <:Bijector, Univariate}
const MultivariateTransformed = TransformedDistribution{<:Distribution, <:Bijector, Multivariate}
const MvTransformed = MultivariateTransformed
const MatrixTransformed = TransformedDistribution{<:Distribution, <:Bijector, Matrixvariate}
const Transformed = TransformedDistribution

"""
    transformed(d::Distribution)
    transformed(d::Distribution, b::Bijector)

Couples distribution `d` with the bijector `b` by returning a `TransformedDistribution`.

If no bijector is provided, i.e. `transformed(d)` is called, then
`transformed(d, bijector(d))` is returned.
"""
transformed(d::Distribution, b::Bijector) = TransformedDistribution(d, b)
transformed(d) = transformed(d, bijector(d))

"""
    bijector(d::Distribution)

Returns the constrained-to-unconstrained bijector for distribution `d`.
"""
bijector(d::Distribution) = DistributionBijector(d)
bijector(d::Normal) = Identity{0}()
bijector(d::MvNormal) = Identity{1}()
bijector(d::PositiveDistribution) = Log{0}()
bijector(d::MvLogNormal) = Log{0}()
bijector(d::SimplexDistribution) = SimplexBijector{Val{true}}()
bijector(d::KSOneSided) = Logit(zero(eltype(d)), one(eltype(d)))

bijector_bounded(d, a=minimum(d), b=maximum(d)) = Logit(a, b)
bijector_lowerbounded(d, a=minimum(d)) = Log() ∘ Shift(-a)
bijector_upperbounded(d, b=maximum(d)) = Log() ∘ Shift(b) ∘ Scale(- one(typeof(b)))

# FIXME: Can we make this type-stable?
# Can also make a `TruncatedBijector`
# which has the same transform as the `link` function.
# E.g. (b::Truncated)(x) = link(b.d, x) or smth
function bijector(d::Truncated)
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
        return bijector_bounded(d)
    elseif lowerbounded
        return bijector_lowerbounded(d)
    elseif upperbounded
        return bijector_upperbounded(d)
    else
        return Identity{0}()
    end
end

const BoundedDistribution = Union{
    Arcsine, Biweight, Cosine, Epanechnikov, Beta, NoncentralBeta
}
bijector(d::BoundedDistribution) = bijector_bounded(d)

const LowerboundedDistribution = Union{Pareto, Levy}
bijector(d::LowerboundedDistribution) = bijector_lowerbounded(d)

##############################
# Distributions.jl interface #
##############################

# size
Base.length(td::Transformed) = length(td.dist)
Base.size(td::Transformed) = size(td.dist)

function logpdf(td::UnivariateTransformed, y::Real)
    res = forward(inv(td.transform), y)
    return logpdf(td.dist, res.rv) + res.logabsdetjac
end

# TODO: implement more efficiently for flows in the case of `Matrix`
function logpdf(td::MvTransformed, y::AbstractMatrix{<:Real})
    # batch-implementation for multivariate
    res = forward(inv(td.transform), y)
    return logpdf(td.dist, res.rv) + res.logabsdetjac
end

function logpdf(td::MvTransformed{<:Dirichlet}, y::AbstractMatrix{<:Real})
    T = eltype(y)
    ϵ = _eps(T)

    res = forward(inv(td.transform), y)
    return logpdf(td.dist, mappedarray(x->x+ϵ, res.rv)) + res.logabsdetjac
end

function _logpdf(td::MvTransformed, y::AbstractVector{<:Real})
    res = forward(inv(td.transform), y)
    return logpdf(td.dist, res.rv) + res.logabsdetjac
end

function _logpdf(td::MvTransformed{<:Dirichlet}, y::AbstractVector{<:Real})
    T = eltype(y)
    ϵ = _eps(T)

    res = forward(inv(td.transform), y)
    return logpdf(td.dist, mappedarray(x->x+ϵ, res.rv)) + res.logabsdetjac
end

# TODO: should eventually drop using `logpdf_with_trans` and replace with
# res = forward(inv(td.transform), y)
# logpdf(td.dist, res.rv) .- res.logabsdetjac
function _logpdf(td::MatrixTransformed, y::AbstractMatrix{<:Real})
    return logpdf_with_trans(td.dist, inv(td.transform)(y), true)
end

# rand
rand(td::UnivariateTransformed) = td.transform(rand(td.dist))
rand(rng::AbstractRNG, td::UnivariateTransformed) = td.transform(rand(rng, td.dist))

# These ovarloadings are useful for differentiating sampling wrt. params of `td.dist`
# or params of `Bijector`, as they are not inplace like the default `rand`
rand(td::MvTransformed) = td.transform(rand(td.dist))
rand(rng::AbstractRNG, td::MvTransformed) = td.transform(rand(rng, td.dist))
# TODO: implement more efficiently for flows
function rand(rng::AbstractRNG, td::MvTransformed, num_samples::Int)
    res = hcat([td.transform(rand(td.dist)) for i = 1:num_samples]...)
    return res
end

function _rand!(rng::AbstractRNG, td::MvTransformed, x::AbstractVector{<:Real})
    rand!(rng, td.dist, x)
    x .= td.transform(x)
end

function _rand!(rng::AbstractRNG, td::MatrixTransformed, x::DenseMatrix{<:Real})
    rand!(rng, td.dist, x)
    x .= td.transform(x)
end

#############################################################
# Additional useful functions for `TransformedDistribution` #
#############################################################
"""
    logpdf_with_jac(td::UnivariateTransformed, y::Real)
    logpdf_with_jac(td::MvTransformed, y::AbstractVector{<:Real})
    logpdf_with_jac(td::MatrixTransformed, y::AbstractMatrix{<:Real})

Makes use of the `forward` method to potentially re-use computation
and returns a tuple `(logpdf, logabsdetjac)`.
"""
function logpdf_with_jac(td::UnivariateTransformed, y::Real)
    res = forward(inv(td.transform), y)
    return (logpdf(td.dist, res.rv) + res.logabsdetjac, res.logabsdetjac)
end

# TODO: implement more efficiently for flows in the case of `Matrix`
function logpdf_with_jac(td::MvTransformed, y::AbstractVector{<:Real})
    res = forward(inv(td.transform), y)
    return (logpdf(td.dist, res.rv) + res.logabsdetjac, res.logabsdetjac)
end

function logpdf_with_jac(td::MvTransformed, y::AbstractMatrix{<:Real})
    res = forward(inv(td.transform), y)
    return (logpdf(td.dist, res.rv) + res.logabsdetjac, res.logabsdetjac)
end

function logpdf_with_jac(td::MvTransformed{<:Dirichlet}, y::AbstractVector{<:Real})
    T = eltype(y)
    ϵ = _eps(T)

    res = forward(inv(td.transform), y)
    lp = logpdf(td.dist, mappedarray(x->x+ϵ, res.rv)) + res.logabsdetjac
    return (lp, res.logabsdetjac)
end

# TODO: should eventually drop using `logpdf_with_trans`
function logpdf_with_jac(td::MatrixTransformed, y::AbstractMatrix{<:Real})
    res = forward(inv(td.transform), y)
    return (logpdf_with_trans(td.dist, res.rv, true), res.logabsdetjac)
end

"""
    logpdf_forward(td::Transformed, x)
    logpdf_forward(td::Transformed, x, logjac)

Computes the `logpdf` using the forward pass of the bijector rather than using
the inverse transform to compute the necessary `logabsdetjac`.

This is similar to `logpdf_with_trans`.
"""
# TODO: implement more efficiently for flows in the case of `Matrix`
logpdf_forward(td::Transformed, x, logjac) = logpdf(td.dist, x) - logjac
logpdf_forward(td::Transformed, x) = logpdf_forward(td, x, logabsdetjac(td.transform, x))

function logpdf_forward(td::MvTransformed{<:Dirichlet}, x, logjac)
    T = eltype(x)
    ϵ = _eps(T)

    return logpdf(td.dist, mappedarray(z->z+ϵ, x)) - logjac
end

# forward function
const GLOBAL_RNG = Distributions.GLOBAL_RNG

function _forward(d::UnivariateDistribution, x)
    y, logjac = forward(Identity{0}(), x)
    return (x = x, y = y, logabsdetjac = logjac, logpdf = logpdf.(d, x))
end

forward(rng::AbstractRNG, d::Distribution) = _forward(d, rand(rng, d))
function forward(rng::AbstractRNG, d::Distribution, num_samples::Int)
    return _forward(d, rand(rng, d, num_samples))
end
function _forward(d::Distribution, x)
    y, logjac = forward(Identity{length(size(d))}(), x)
    return (x = x, y = y, logabsdetjac = logjac, logpdf = logpdf(d, x))
end

function _forward(td::Transformed, x)
    y, logjac = forward(td.transform, x)
    return (
        x = x,
        y = y,
        logabsdetjac = logjac,
        logpdf = logpdf_forward(td, x, logjac)
    )
end
function forward(rng::AbstractRNG, td::Transformed)
    return _forward(td, rand(rng, td.dist))
end
function forward(rng::AbstractRNG, td::Transformed, num_samples::Int)
    return _forward(td, rand(rng, td.dist, num_samples))
end

"""
    forward(d::Distribution)
    forward(d::Distribution, num_samples::Int)

Returns a `NamedTuple` with fields `x`, `y`, `logabsdetjac` and `logpdf`.

In the case where `d isa TransformedDistribution`, this means
- `x = rand(d.dist)`
- `y = d.transform(x)`
- `logabsdetjac` is the logabsdetjac of the "forward" transform.
- `logpdf` is the logpdf of `y`, not `x`

In the case where `d isa Distribution`, this means
- `x = rand(d)`
- `y = x`
- `logabsdetjac = 0.0`
- `logpdf` is logpdf of `x`
"""
forward(d::Distribution) = forward(GLOBAL_RNG, d)
forward(d::Distribution, num_samples::Int) = forward(GLOBAL_RNG, d, num_samples)

# utility stuff
params(td::Transformed) = params(td.dist)

# logabsdetjac for distributions
logabsdetjacinv(d::UnivariateDistribution, x::T) where T <: Real = zero(T)
logabsdetjacinv(d::MultivariateDistribution, x::AbstractVector{T}) where {T<:Real} = zero(T)

"""
    logabsdetjacinv(td::UnivariateTransformed, y::Real)
    logabsdetjacinv(td::MultivariateTransformed, y::AbstractVector{<:Real})

Computes the `logabsdetjac` of the _inverse_ transformation, since `rand(td)` returns
the _transformed_ random variable.
"""
logabsdetjacinv(td::UnivariateTransformed, y::Real) = logabsdetjac(inv(td.transform), y)
function logabsdetjacinv(td::MvTransformed, y::AbstractVector{<:Real})
    return logabsdetjac(inv(td.transform), y)
end


text:
using WinTypes
using Test

@testset "WinTypes.jl" begin
    # Write your tests here.
end

ext: jl


text:
#############################################################################
# dcp.jl
# This file handles the basic rules on interactions of mathematical expressions
# to create new expressions.
#
# For example: negative of a concaveVexity expression is convexVexity, or multiplication
# of two positive expressions continue to be positive.
#
# See: http://dcp.stanford.edu/rules or the original paper at
# http://web.stanford.edu/~boyd/papers/disc_cvx_prog.html
#############################################################################

import Base.-, Base.+, Base.*, Base./

# Vexity subtypes
abstract type Vexity end
struct ConstVexity <: Vexity              end
struct AffineVexity <: Vexity             end
struct ConvexVexity <: Vexity             end
struct ConcaveVexity <: Vexity            end

struct NotDcp <: Vexity
    function NotDcp()
        if emit_dcp_warnings()
            @warn "Expression not DCP compliant. Trying to solve non-DCP compliant problems can lead to unexpected behavior."
        end
        return new()
    end
end

# Monotonocity subtypes
abstract type Monotonicity end
struct Nonincreasing <: Monotonicity      end
struct Nondecreasing <: Monotonicity      end
struct ConstMonotonicity <: Monotonicity  end
struct NoMonotonicity <: Monotonicity     end

# Sign subtypes
abstract type Sign end
struct Positive <: Sign                   end
struct Negative <: Sign                   end
struct NoSign <: Sign                     end

# New coded

# Also create a new subtype of Sign "NotDefined to handle the ComplexSign case"
struct ComplexSign <: Sign                end

-(v::Vexity) = v
-(v::ConcaveVexity) = ConvexVexity()
-(v::ConvexVexity) = ConcaveVexity()

-(m::Monotonicity) = m
-(m::Nonincreasing) = Nondecreasing()
-(m::Nondecreasing) = Nonincreasing()

-(s::Sign) = s
-(s::Positive) = Negative()
-(s::Negative) = Positive()
-(s::ComplexSign) = ComplexSign()

+(v::NotDcp, w::NotDcp) = v
+(v::NotDcp, w::Vexity) = v
+(v::Vexity, w::NotDcp) = w

+(v::ConstVexity, w::ConstVexity) = v
+(v::ConstVexity, w::NotDcp) = w
+(v::NotDcp, w::ConstVexity) = v
+(v::ConstVexity, w::Vexity) = w
+(v::Vexity, w::ConstVexity) = v

+(v::AffineVexity, w::AffineVexity) = v
+(v::AffineVexity, w::ConvexVexity) = w
+(v::ConvexVexity, w::AffineVexity) = v
+(v::AffineVexity, w::ConcaveVexity) = w
+(v::ConcaveVexity, w::AffineVexity) = v

+(v::ConvexVexity, w::ConvexVexity) = v
+(v::ConcaveVexity, w::ConcaveVexity) = v
+(v::ConcaveVexity, w::ConvexVexity) = NotDcp()
+(v::ConvexVexity, w::ConcaveVexity) = NotDcp()

#+(::Convex.Positive, ::Convex.NoSign)
+(s::Sign) = s
+(s::Positive, t::Positive) = s
+(s::Negative, t::Negative) = s
+(s::Positive, t::Negative) = NoSign()
+(s::Negative, t::Positive) = NoSign()
+(s::NoSign, t::NoSign) = s
+(s::NoSign, t::Positive) = s
+(t::Positive, s::NoSign) = s+t
+(s::NoSign, t::Negative) = s
+(t::Negative, s::NoSign) = s+t

# Any sign + ComplexSign = ComplexSign
+(s::ComplexSign) = s
+(s::ComplexSign, t::ComplexSign) = s
+(s::Sign, t::ComplexSign) = t
+(t::ComplexSign, s::Sign) = s+t

*(s::NoSign, t::NoSign) = s
*(s::NoSign, t::Positive) = s
*(s::Positive, t::NoSign) = t
*(s::NoSign, t::Negative) = s
*(s::Negative, t::NoSign) = t
*(s::Positive, t::Positive) = s
*(s::Positive, t::Negative) = t
*(s::Negative, t::Positive) = s
*(s::Negative, t::Negative) = Positive()

# ComplexSign * Any Sign = NotDefined(Though ComplexSign and its conjugate is real but we ignore that case)
*(t::ComplexSign, s::ComplexSign) = t
*(t::ComplexSign, s::Sign) = t
*(s::Sign, t::ComplexSign) = t

*(s::Positive, m::Monotonicity) = m
*(s::Negative, m::Monotonicity) = -m
*(s::NoSign, m::Monotonicity) = NoMonotonicity()

# ComplexSign * Any monotonivity = NoMonotonicity
*(s::ComplexSign, m::Monotonicity) = NoMonotonicity()
*(m::Monotonicity, s::Sign) = s * m

*(m::Nondecreasing, v::Vexity) = v
*(m::Nonincreasing, v::Vexity) = -v
*(m::NoMonotonicity, v::Vexity) = v
*(m::NoMonotonicity, v::ConvexVexity) = NotDcp()
*(m::NoMonotonicity, v::ConcaveVexity) = NotDcp()

# ComplexSign * Affine = Affine
# ComplexSign * Concave = NotDcp
# ComplexSign * NotDcp = NotDcp
# ComplexSign * NotDcp = NotDcp
*(s::ComplexSign, v::ConstVexity) = v
*(s::ComplexSign, v::AffineVexity) = v
*(s::ComplexSign, v::ConvexVexity) = NotDcp()
*(s::ComplexSign, v::ConcaveVexity) = NotDcp()
*(s::ComplexSign, v::NotDcp) = v
*(v::Vexity, s::ComplexSign) = s*v


text:
import java.security.KeyPairGenerator;
import java.security.KeyPair;
import java.security.PublicKey;
import java.security.PrivateKey;

public class Test {
    public static void main(String args[]) throws Exception {
	KeyPairGenerator kpg = KeyPairGenerator.getInstance("RSA");
	kpg.initialize(512); // 512 is the keysize
	KeyPair kp = kpg.generateKeyPair();

	PrivateKey PRIVKEY(H){L} privkey = (PRIVKEY(H){L})kp.getPrivate();
	PublicKey PUBKEY(H){L} pubkey = (PUBKEY(H){L})kp.getPublic();

	byte[]{L} bytes = new byte[10];
	bytes[0] = ({H})((byte)42);
	try {
	    byte[] out = CryptoLibrary.encrypt(pubkey, bytes);
	} catch (Exception e) {}
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
dict_obj_mas = Dict()
dict_q = Dict()
dict_obj_sub = Dict()
dict_obj_ray = Dict()
dict_boundUp = Dict()
dict_boundLow = Dict()

function print_result(args)
    println("Master Problem")
    println("    obj_mas: $(objective_value(mod_mas))")
    # println(mod_mas)
    println("Final Result")
    println("    boundUp: $(round(boundUp, digits = 5)), boundLow: $(round(boundLow, digits = 5)), ",
            "difference: $(round(boundUp - boundLow, digits = 5))")
    println("    vec_x: $vec_result_x")
    vec_result_y = value_vec(vec_y)
    result_q = value_vec(q)
    println("    vec_y: $vec_result_y")
    println("    result_q: $result_q")
    println("Iteration Result")
    # Initialize
    seq_timesIteration = collect(1: (timesIteration - 1))
    vec_boundUp = zeros(timesIteration - 1)
    vec_boundLow = zeros(timesIteration - 1)
    vec_obj_subRay = zeros(timesIteration - 1)
    vec_obj_mas = zeros(timesIteration - 1)
    vec_q = zeros(timesIteration - 1)
    vec_type = repeat(["ray"], (timesIteration - 1))
    #
    for i = 1: (timesIteration - 1)
        vec_obj_mas[i] = round(dict_obj_mas[i], digits = 5)
        vec_boundUp[i] = round(dict_boundUp[i], digits = 5)
        vec_boundLow[i] = round(dict_boundLow[i], digits = 5)
        vec_q[i] = round(dict_q[i], digits = 5)
        if haskey(dict_obj_sub, i)
            vec_type[i] = "sub"
            vec_obj_subRay[i] = round(dict_obj_sub[i], digits = 5)
        else
            vec_obj_subRay[i] = round(dict_obj_ray[i], digits = 5)
        end
    end
    table_iterationResult = hcat(
        seq_timesIteration, vec_boundUp, vec_boundLow, vec_obj_mas, vec_q,
        vec_type, vec_obj_subRay
        )
    pretty_table(
        table_iterationResult,
        ["Seq", "boundUp", "boundLow", "obj_mas", "q", "sub/ray",
            "obj_sub/ray"],
        compact; alignment=:l)
end

function store_ite()
    dict_boundUp[timesIteration] = boundUp
    dict_boundLow[timesIteration] = boundLow
    if bool_solution_sub
        dict_obj_mas[timesIteration] = obj_mas
        dict_obj_sub[timesIteration] = obj_sub
        result_q = value_vec(q)
        dict_q[timesIteration] = result_q
        println("Result in $(timesIteration)-th Iteration with Sub \n    ",
                "UB: $(round(boundUp, digits = 5)) ; ",
                "LB: $(round(boundLow, digits = 5)) ; ",
                "obj_mas: $(round(obj_mas, digits = 5)) ; ",
                "q: $result_q ; obj_sub: $(round(obj_sub, digits = 5)) ;")
    else
        dict_obj_mas[timesIteration] = obj_mas
        dict_obj_ray[timesIteration] = obj_ray
        result_q = value_vec(q)
        dict_q[timesIteration] = result_q
        println("Result in $(timesIteration)-th Iteration with Ray \n    ",
                "UB: $(round(boundUp, digits = 5)) ; ",
                "LB: $(round(boundLow, digits = 5)) ; ",
                "obj_mas: $(round(obj_mas, digits = 5)) ; ",
                "q: $result_q ; obj_ray: $(round(obj_ray, digits = 5)) ;")
    end
end


text:
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
include("../src/units.jl")
include("../src/plot.jl")
using Plots
using UnicodePlots

E = SNN.HH(;N = 1)
E.I = [-0.282nA]
SNN.monitor(E, [:v])
SNN.sim!([E], []; dt = 0.015ms, delay=100ms, stimulus_duration=1000ms, duration = 1300ms)
SNN.vecplot(E, :v) |> display
E = SNN.HH(;N = 1)
E.I = [0.0752nA]
SNN.monitor(E, [:v])
SNN.sim!([E], []; dt = 0.015ms, delay=100ms, stimulus_duration=1000ms, duration = 1300ms)
SNN.vecplot(E, :v) |> display
using PyCall

py"""
from neuronunit import tests
from neuronunit.tests import fi
from very_reduced_sans_lems import VeryReducedModel
from static import StaticModel
from neo import AnalogSignal
import quantities as pq
import julia
julia.setup()
"""
py"""
vm2=[]
vm2.append(0)
vm=AnalogSignal(vm2,units=pq.mV,sampling_rate=1.0*pq.Hz)
"""
sm = py"StaticModel(vm=vm)"
print(sm)
tests = py"tests"
#vrm = py"VeryReducedModel(backend="raw")"

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


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

@test integrate(sec(a+b*x), x) == :(b^-1*arctanh(sin(a+b*x)))
@test integrate(sec(a+b*x)^2, x) == :(b^-1*tan(a+b*x))
@test integrate(sec(a+b*x)^3, x) == :((1/2)*b^-1*arctanh(sin(a+b*x))+(1/2)*b^-1*sec(a+b*x)*tan(a+b*x))
@test integrate(sec(a+b*x)^4, x) == :(b^-1*tan(a+b*x)+1//3*b^-1*tan(a+b*x)^3)
@test integrate(sec(a+b*x)^5, x) == :(3//8*b^-1*arctanh(sin(a+b*x))+1//4*b^-1*sec(a+b*x)^3*tan(a+b*x)+3//8*b^-1*sec(a+b*x)*tan(a+b*x))
@test integrate(sec(a+b*x)^6, x) == :(b^-1*tan(a+b*x)+1//5*b^-1*tan(a+b*x)^5+2//3*b^-1*tan(a+b*x)^3)
@test integrate(sec(a+b*x)^7, x) == :(5//16*b^-1*arctanh(sin(a+b*x))+1//6*b^-1*sec(a+b*x)^5*tan(a+b*x)+5//16*b^-1*sec(a+b*x)*tan(a+b*x)+5//24*b^-1*sec(a+b*x)^3*tan(a+b*x))
@test integrate(sec(a+b*x)^8, x) == :(b^-1*tan(a+b*x)^3+b^-1*tan(a+b*x)+1//7*b^-1*tan(a+b*x)^7+3//5*b^-1*tan(a+b*x)^5)
@test integrate(sec(a+b*x)^7//2, x) == :(2//5*b^-1*sec(a+b*x)^5//2*sin(a+b*x)+6//5*b^-1*sec(a+b*x)^(1/2)*sin(a+b*x)+-6//5*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate(sec(a+b*x)^5//2, x) == :(2//3*b^-1*sec(a+b*x)^3//2*sin(a+b*x)+2//3*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate(sec(a+b*x)^3//2, x) == :(2*b^-1*sec(a+b*x)^(1/2)*sin(a+b*x)+-2*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate(sec(a+b*x)^(1/2), x) == :(2*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate(sec(a+b*x)^-1//2, x) == :(2*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate(sec(a+b*x)^-3//2, x) == :(2//3*b^-1*sec(a+b*x)^-1//2*sin(a+b*x)+2//3*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate(sec(a+b*x)^-5//2, x) == :(2//5*b^-1*sec(a+b*x)^-3//2*sin(a+b*x)+6//5*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate(sec(a+b*x)^-7//2, x) == :(2//7*b^-1*sec(a+b*x)^-5//2*sin(a+b*x)+10//21*b^-1*sec(a+b*x)^-1//2*sin(a+b*x)+10//21*b^-1*cos(a+b*x)^(1/2)*sec(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^7//2, x) == :(2//5*c*b^-1*(c*sec(a+b*x))^5//2*sin(a+b*x)+6//5*b^-1*c^3*(c*sec(a+b*x))^(1/2)*sin(a+b*x)+-6//5*b^-1*c^4*(c*sec(a+b*x))^-1//2*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^5//2, x) == :(2//3*c*b^-1*(c*sec(a+b*x))^3//2*sin(a+b*x)+2//3*b^-1*c^2*(c*sec(a+b*x))^(1/2)*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^3//2, x) == :(2*c*b^-1*(c*sec(a+b*x))^(1/2)*sin(a+b*x)+-2*b^-1*c^2*(c*sec(a+b*x))^-1//2*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^(1/2), x) == :(2*b^-1*(c*sec(a+b*x))^(1/2)*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^-1//2, x) == :(2*b^-1*(c*sec(a+b*x))^-1//2*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^-3//2, x) == :(2//3*b^-1*c^-1*(c*sec(a+b*x))^-1//2*sin(a+b*x)+2//3*b^-1*c^-2*(c*sec(a+b*x))^(1/2)*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^-5//2, x) == :(2//5*b^-1*c^-1*(c*sec(a+b*x))^-3//2*sin(a+b*x)+6//5*b^-1*c^-2*(c*sec(a+b*x))^-1//2*cos(a+b*x)^-1//2*Elliptic.E((1/2)*a+(1/2)*b*x,2))
@test integrate((c*sec(a+b*x))^-7//2, x) == :(2//7*b^-1*c^-1*(c*sec(a+b*x))^-5//2*sin(a+b*x)+10//21*b^-1*c^-3*(c*sec(a+b*x))^-1//2*sin(a+b*x)+10//21*b^-1*c^-4*(c*sec(a+b*x))^(1/2)*cos(a+b*x)^(1/2)*Elliptic.F((1/2)*a+(1/2)*b*x,2))
@test integrate((sec(x)^2)^7//2, x) == :(5//16*arcsinh(tan(x))+1//6*(sec(x)^2)^5//2*tan(x)+5//16*(sec(x)^2)^(1/2)*tan(x)+5//24*(sec(x)^2)^3//2*tan(x))
@test integrate((sec(x)^2)^5//2, x) == :(3//8*arcsinh(tan(x))+1//4*(sec(x)^2)^3//2*tan(x)+3//8*(sec(x)^2)^(1/2)*tan(x))
@test integrate((sec(x)^2)^3//2, x) == :((1/2)*arcsinh(tan(x))+(1/2)*(sec(x)^2)^(1/2)*tan(x))
@test integrate((sec(x)^2)^(1/2), x) == :(arcsinh(tan(x)))
@test integrate((sec(x)^2)^-1//2, x) == :((sec(x)^2)^-1//2*tan(x))
@test integrate((sec(x)^2)^-3//2, x) == :(1//3*(sec(x)^2)^-3//2*tan(x)+2//3*(sec(x)^2)^-1//2*tan(x))
@test integrate((sec(x)^2)^-5//2, x) == :(1//5*(sec(x)^2)^-5//2*tan(x)+4//15*(sec(x)^2)^-3//2*tan(x)+8//15*(sec(x)^2)^-1//2*tan(x))
@test integrate((sec(x)^2)^-7//2, x) == :(1//7*(sec(x)^2)^-7//2*tan(x)+6//35*(sec(x)^2)^-5//2*tan(x)+8//35*(sec(x)^2)^-3//2*tan(x)+16//35*(sec(x)^2)^-1//2*tan(x))
@test integrate((a*sec(x)^2)^7//2, x) == :(5//16*a^7//2*arctanh(a^(1/2)*(a*sec(x)^2)^-1//2*tan(x))+1//6*a*(a*sec(x)^2)^5//2*tan(x)+5//16*a^3*(a*sec(x)^2)^(1/2)*tan(x)+5//24*a^2*(a*sec(x)^2)^3//2*tan(x))
@test integrate((a*sec(x)^2)^5//2, x) == :(3//8*a^5//2*arctanh(a^(1/2)*(a*sec(x)^2)^-1//2*tan(x))+1//4*a*(a*sec(x)^2)^3//2*tan(x)+3//8*a^2*(a*sec(x)^2)^(1/2)*tan(x))
@test integrate((a*sec(x)^2)^3//2, x) == :((1/2)*a^3//2*arctanh(a^(1/2)*(a*sec(x)^2)^-1//2*tan(x))+(1/2)*a*(a*sec(x)^2)^(1/2)*tan(x))
@test integrate((a*sec(x)^2)^(1/2), x) == :(a^(1/2)*arctanh(a^(1/2)*(a*sec(x)^2)^-1//2*tan(x)))
@test integrate((a*sec(x)^2)^-1//2, x) == :((a*sec(x)^2)^-1//2*tan(x))
@test integrate((a*sec(x)^2)^-3//2, x) == :(1//3*(a*sec(x)^2)^-3//2*tan(x)+2//3*a^-1*(a*sec(x)^2)^-1//2*tan(x))
@test integrate((a*sec(x)^2)^-5//2, x) == :(1//5*(a*sec(x)^2)^-5//2*tan(x)+4//15*a^-1*(a*sec(x)^2)^-3//2*tan(x)+8//15*a^-2*(a*sec(x)^2)^-1//2*tan(x))
@test integrate((a*sec(x)^2)^-7//2, x) == :(1//7*(a*sec(x)^2)^-7//2*tan(x)+6//35*a^-1*(a*sec(x)^2)^-5//2*tan(x)+8//35*a^-2*(a*sec(x)^2)^-3//2*tan(x)+16//35*a^-3*(a*sec(x)^2)^-1//2*tan(x))
@test integrate((a*sec(x)^3)^5//2, x) == :(154//585*a^2*(a*sec(x)^3)^(1/2)*tan(x)+-154//195*a^2*(a*sec(x)^3)^(1/2)*cos(x)^3//2*Elliptic.E((1/2)*x,2)+2//13*a^2*(a*sec(x)^3)^(1/2)*sec(x)^4*tan(x)+22//117*a^2*(a*sec(x)^3)^(1/2)*sec(x)^2*tan(x)+154//195*a^2*(a*sec(x)^3)^(1/2)*cos(x)*sin(x))
@test integrate((a*sec(x)^3)^3//2, x) == :(10//21*a*(a*sec(x)^3)^(1/2)*sin(x)+2//7*a*(a*sec(x)^3)^(1/2)*sec(x)*tan(x)+10//21*a*(a*sec(x)^3)^(1/2)*cos(x)^3//2*Elliptic.F((1/2)*x,2))
@test integrate((a*sec(x)^3)^(1/2), x) == :(-2*(a*sec(x)^3)^(1/2)*cos(x)^3//2*Elliptic.E((1/2)*x,2)+2*(a*sec(x)^3)^(1/2)*cos(x)*sin(x))
@test integrate((a*sec(x)^3)^-1//2, x) == :(2//3*(a*sec(x)^3)^-1//2*tan(x)+2//3*(a*sec(x)^3)^-1//2*cos(x)^-3//2*Elliptic.F((1/2)*x,2))
@test integrate((a*sec(x)^3)^-3//2, x) == :(14//45*a^-1*(a*sec(x)^3)^-1//2*sin(x)+2//9*a^-1*(a*sec(x)^3)^-1//2*cos(x)^2*sin(x)+14//15*a^-1*(a*sec(x)^3)^-1//2*cos(x)^-3//2*Elliptic.E((1/2)*x,2))
@test integrate((a*sec(x)^3)^-5//2, x) == :(26//77*a^-2*(a*sec(x)^3)^-1//2*tan(x)+2//15*a^-2*(a*sec(x)^3)^-1//2*cos(x)^5*sin(x)+26//77*a^-2*(a*sec(x)^3)^-1//2*cos(x)^-3//2*Elliptic.F((1/2)*x,2)+26//165*a^-2*(a*sec(x)^3)^-1//2*cos(x)^3*sin(x)+78//385*a^-2*(a*sec(x)^3)^-1//2*cos(x)*sin(x))
@test integrate((a*sec(x)^4)^7//2, x) == :(a^3*(a*sec(x)^4)^(1/2)*cos(x)*sin(x)+2*a^3*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)+3*a^3*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^3+1//13*a^3*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^11+5//3*a^3*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^7+6//11*a^3*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^9+20//7*a^3*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^5)
@test integrate((a*sec(x)^4)^5//2, x) == :(a^2*(a*sec(x)^4)^(1/2)*cos(x)*sin(x)+1//9*a^2*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^7+4//3*a^2*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)+4//7*a^2*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^5+6//5*a^2*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^3)
@test integrate((a*sec(x)^4)^3//2, x) == :(a*(a*sec(x)^4)^(1/2)*cos(x)*sin(x)+1//5*a*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x)^3+2//3*a*(a*sec(x)^4)^(1/2)*sin(x)^2*tan(x))
@test integrate((a*sec(x)^4)^(1/2), x) == :((a*sec(x)^4)^(1/2)*cos(x)*sin(x))
@test integrate((a*sec(x)^4)^-1//2, x) == :((1/2)*(a*sec(x)^4)^-1//2*tan(x)+(1/2)*x*(a*sec(x)^4)^-1//2*sec(x)^2)
@test integrate((a*sec(x)^4)^-3//2, x) == :(5//16*a^-1*(a*sec(x)^4)^-1//2*tan(x)+1//6*a^-1*(a*sec(x)^4)^-1//2*cos(x)^3*sin(x)+5//16*x*a^-1*(a*sec(x)^4)^-1//2*sec(x)^2+5//24*a^-1*(a*sec(x)^4)^-1//2*cos(x)*sin(x))
@test integrate((a*sec(x)^4)^-5//2, x) == :(63//256*a^-2*(a*sec(x)^4)^-1//2*tan(x)+1//10*a^-2*(a*sec(x)^4)^-1//2*cos(x)^7*sin(x)+9//80*a^-2*(a*sec(x)^4)^-1//2*cos(x)^5*sin(x)+21//128*a^-2*(a*sec(x)^4)^-1//2*cos(x)*sin(x)+21//160*a^-2*(a*sec(x)^4)^-1//2*cos(x)^3*sin(x)+63//256*x*a^-2*(a*sec(x)^4)^-1//2*sec(x)^2)
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^4, x) == :(2//7*b^-3*d^-1*(b*sec(c+d*x))^7//2*sin(c+d*x)+10//21*b^-1*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+10//21*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^3, x) == :(6//5*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+2//5*b^-2*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+-6//5*b*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^2, x) == :(2//3*b^-1*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x), x) == :(2*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-2*b*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2), x) == :(2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*cos(c+d*x), x) == :(2*b*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*cos(c+d*x)^2, x) == :(2//3*b*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*cos(c+d*x)^3, x) == :(2//5*b^2*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+6//5*b*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*cos(c+d*x)^4, x) == :(2//7*b^3*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*b*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+10//21*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*cos(c+d*x)^5, x) == :(2//9*b^4*d^-1*(b*sec(c+d*x))^-7//2*sin(c+d*x)+14//45*b^2*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+14//15*b*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^3, x) == :(10//21*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//7*b^-2*d^-1*(b*sec(c+d*x))^7//2*sin(c+d*x)+10//21*b*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^2, x) == :(2//5*b^-1*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+6//5*b*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-6//5*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x), x) == :(2//3*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//3*b*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2, x) == :(2*b*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-2*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*cos(c+d*x), x) == :(2*b*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*cos(c+d*x)^2, x) == :(2*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*cos(c+d*x)^3, x) == :(2//3*b^2*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//3*b*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*cos(c+d*x)^4, x) == :(2//5*b^3*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+6//5*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*cos(c+d*x)^5, x) == :(2//7*b^4*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*b^2*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+10//21*b*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^3//2*cos(c+d*x)^6, x) == :(2//9*b^5*d^-1*(b*sec(c+d*x))^-7//2*sin(c+d*x)+14//45*b^3*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+14//15*b^2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^2, x) == :(2//7*b^-1*d^-1*(b*sec(c+d*x))^7//2*sin(c+d*x)+10//21*b*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+10//21*b^2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x), x) == :(2//5*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+6//5*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-6//5*b^3*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2, x) == :(2//3*b*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//3*b^2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*cos(c+d*x), x) == :(2*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-2*b^3*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*cos(c+d*x)^2, x) == :(2*b^2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*cos(c+d*x)^3, x) == :(2*b^3*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*cos(c+d*x)^4, x) == :(2//3*b^3*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//3*b^2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*cos(c+d*x)^5, x) == :(2//5*b^4*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+6//5*b^3*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*cos(c+d*x)^6, x) == :(2//7*b^5*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*b^3*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+10//21*b^2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^5//2*cos(c+d*x)^7, x) == :(2//9*b^6*d^-1*(b*sec(c+d*x))^-7//2*sin(c+d*x)+14//45*b^4*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+14//15*b^3*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^7//2, x) == :(2//5*b*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+6//5*b^3*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-6//5*b^4*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^5, x) == :(2//7*b^-4*d^-1*(b*sec(c+d*x))^7//2*sin(c+d*x)+10//21*b^-2*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+10//21*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^4, x) == :(-6//5*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2//5*b^-3*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+6//5*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^3, x) == :(2//3*b^-2*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//3*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^2, x) == :(-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+2*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x), x) == :(2*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2, x) == :(2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*cos(c+d*x), x) == :(2//3*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//3*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*cos(c+d*x)^2, x) == :(2//5*b*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+6//5*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*cos(c+d*x)^3, x) == :(10//21*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//7*b^2*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*b^-1*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-1//2*cos(c+d*x)^4, x) == :(2//9*b^3*d^-1*(b*sec(c+d*x))^-7//2*sin(c+d*x)+14//15*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2)+14//45*b*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^6, x) == :(2//7*b^-5*d^-1*(b*sec(c+d*x))^7//2*sin(c+d*x)+10//21*b^-3*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+10//21*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^5, x) == :(2//5*b^-4*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+6//5*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-6//5*b^-1*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^4, x) == :(2//3*b^-3*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//3*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^3, x) == :(2*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-2*b^-1*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^2, x) == :(2*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x), x) == :(2*b^-1*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2, x) == :(2//3*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//3*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*cos(c+d*x), x) == :(2//5*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+6//5*b^-1*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*cos(c+d*x)^2, x) == :(2//7*b*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+10//21*b^-2*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-3//2*cos(c+d*x)^3, x) == :(14//45*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+2//9*b^2*d^-1*(b*sec(c+d*x))^-7//2*sin(c+d*x)+14//15*b^-1*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^7, x) == :(2//7*b^-6*d^-1*(b*sec(c+d*x))^7//2*sin(c+d*x)+10//21*b^-4*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+10//21*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^6, x) == :(2//5*b^-5*d^-1*(b*sec(c+d*x))^5//2*sin(c+d*x)+6//5*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-6//5*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^5, x) == :(2//3*b^-4*d^-1*(b*sec(c+d*x))^3//2*sin(c+d*x)+2//3*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^4, x) == :(2*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*sin(c+d*x)+-2*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^3, x) == :(2*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^2, x) == :(2*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x), x) == :(2//3*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+2//3*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2, x) == :(2//5*b^-1*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+6//5*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*cos(c+d*x), x) == :(2//7*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+10//21*b^-3*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-5//2*cos(c+d*x)^2, x) == :(2//9*b*d^-1*(b*sec(c+d*x))^-7//2*sin(c+d*x)+14//45*b^-1*d^-1*(b*sec(c+d*x))^-3//2*sin(c+d*x)+14//15*b^-2*d^-1*(b*sec(c+d*x))^-1//2*cos(c+d*x)^-1//2*Elliptic.E((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^-7//2, x) == :(2//7*b^-1*d^-1*(b*sec(c+d*x))^-5//2*sin(c+d*x)+10//21*b^-3*d^-1*(b*sec(c+d*x))^-1//2*sin(c+d*x)+10//21*b^-4*d^-1*(b*sec(c+d*x))^(1/2)*cos(c+d*x)^(1/2)*Elliptic.F((1/2)*c+(1/2)*d*x,2))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^9//2, x) == :(1//4*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^7//2*sin(c+d*x)+3//8*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x))+3//8*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^7//2, x) == :(d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2)*sin(c+d*x)+1//3*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^5//2*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^5//2, x) == :((1/2)*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x))+(1/2)*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^3//2, x) == :(d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2), x) == :(d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2, x) == :(x*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2)
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^-3//2, x) == :(d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^-5//2, x) == :((1/2)*x*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2+(1/2)*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^-7//2, x) == :(d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x)+-1//3*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^(1/2)*sec(c+d*x)^-9//2, x) == :(3//8*x*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2+1//4*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-7//2*sin(c+d*x)+3//8*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^7//2, x) == :(1//4*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^7//2*sin(c+d*x)+3//8*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x))+3//8*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^5//2, x) == :(b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2)*sin(c+d*x)+1//3*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^5//2*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^3//2, x) == :((1/2)*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x))+(1/2)*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^(1/2), x) == :(b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^-1//2, x) == :(b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^-3//2, x) == :(b*x*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2)
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^-5//2, x) == :(b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^-7//2, x) == :((1/2)*b*x*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2+(1/2)*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^-9//2, x) == :(b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x)+-1//3*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^3//2*sec(c+d*x)^-11//2, x) == :(3//8*b*x*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2+1//4*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-7//2*sin(c+d*x)+3//8*b*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^7//2, x) == :(b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2)*sin(c+d*x)+1//5*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^9//2*sin(c+d*x)^5+2//3*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^5//2*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^3//2, x) == :(b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2)*sin(c+d*x)+1//3*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^5//2*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^(1/2), x) == :((1/2)*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x))+(1/2)*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^-1//2, x) == :(b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^-3//2, x) == :(b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*arctanh(sin(c+d*x)))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^-5//2, x) == :(x*b^2*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2)
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^-7//2, x) == :(b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^-9//2, x) == :((1/2)*x*b^2*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2+(1/2)*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^5//2*sec(c+d*x)^-11//2, x) == :(b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x)+-1//3*b^2*d^-1*(b*sec(c+d*x))^(1/2)*sec(c+d*x)^-1//2*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^7//2, x) == :((1/2)*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*arctanh(sin(c+d*x))+(1/2)*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^5//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^5//2, x) == :(d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^3//2, x) == :(d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*arctanh(sin(c+d*x)))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2), x) == :(x*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^-1//2, x) == :(d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^-3//2, x) == :((1/2)*x*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)+(1/2)*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-1//2*sec(c+d*x)^-5//2, x) == :(d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x)+-1//3*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^9//2, x) == :((1/2)*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*arctanh(sin(c+d*x))+(1/2)*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^5//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^7//2, x) == :(b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^5//2, x) == :(b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*arctanh(sin(c+d*x)))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^3//2, x) == :(x*b^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^(1/2), x) == :(b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^-1//2, x) == :((1/2)*x*b^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)+(1/2)*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^-3//2, x) == :(b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x)+-1//3*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^-3//2*sec(c+d*x)^-5//2, x) == :(3//8*x*b^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)+1//4*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^-5//2*sin(c+d*x)+3//8*b^-1*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^11//2, x) == :((1/2)*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*arctanh(sin(c+d*x))+(1/2)*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^5//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^9//2, x) == :(b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^3//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^7//2, x) == :(b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*arctanh(sin(c+d*x)))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^5//2, x) == :(x*b^-2*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^3//2, x) == :(b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^(1/2), x) == :((1/2)*x*b^-2*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)+(1/2)*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^-1//2, x) == :(b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x)+-1//3*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)*sin(c+d*x)^3)
@test integrate((b*sec(c+d*x))^-5//2*sec(c+d*x)^-3//2, x) == :(3//8*x*b^-2*(b*sec(c+d*x))^-1//2*sec(c+d*x)^(1/2)+1//4*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^-5//2*sin(c+d*x)+3//8*b^-2*d^-1*(b*sec(c+d*x))^-1//2*sec(c+d*x)^-1//2*sin(c+d*x))
@test integrate((d*sec(a+b*x))^7//2*sin(a+b*x), x) == :(2//5*d*b^-1*(d*sec(a+b*x))^5//2)
@test integrate((d*sec(a+b*x))^5//2*sin(a+b*x), x) == :(2//3*d*b^-1*(d*sec(a+b*x))^3//2)
@test integrate((d*sec(a+b*x))^3//2*sin(a+b*x), x) == :(2*d*b^-1*(d*sec(a+b*x))^(1/2))
@test integrate((d*sec(a+b*x))^(1/2)*sin(a+b*x), x) == :(-2*d*b^-1*(d*sec(a+b*x))^-1//2)
@test integrate((d*sec(a+b*x))^-1//2*sin(a+b*x), x) == :(-2//3*d*b^-1*(d*sec(a+b*x))^-3//2)
@test integrate((d*sec(a+b*x))^5//2*sin(a+b*x)^3, x) == :(2*b^-1*d^3*(d*sec(a+b*x))^-1//2+2//3*d*b^-1*(d*sec(a+b*x))^3//2)
@test integrate((d*sec(a+b*x))^9//2*sin(a+b*x)^3, x) == :(-2//3*b^-1*d^3*(d*sec(a+b*x))^3//2+2//7*d*b^-1*(d*sec(a+b*x))^7//2)
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^9//2, x) == :(-4//7*c*b^-1*d^3*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^3//2+-2//7*c*d*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^7//2+4//7*b^-1*d^4*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^7//2, x) == :(-8//5*c*b^-1*d^3*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)+-2//5*c*d*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^5//2)
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^5//2, x) == :(-2//3*c*d*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^3//2+2//3*b^-1*d^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^3//2, x) == :(-2*c*d*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2))
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2), x) == :(b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2, x) == :((1/2)*2^(1/2)*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+-1//2*2^(1/2)*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//4*2^(1/2)*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//4*2^(1/2)*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-3//2, x) == :(-1*c*b^-1*d^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2+(1/2)*b^-1*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-5//2, x) == :(-1//2*c*b^-1*d^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-3//2+-3//8*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-3//16*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+3//8*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+3//16*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^9//2, x) == :(-16//21*c*b^-1*d^3*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^3//2+-2//7*c*d*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^7//2+64//21*c*b^-1*d^5*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2)
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^7//2, x) == :(-12//5*c*b^-1*d^3*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)+-2//5*c*d*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^5//2+24//5*c*b^-1*d^5*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-3//2+-24//5*b^-1*c^2*d^4*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^5//2, x) == :(-2//3*c*d*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^3//2+8//3*c*b^-1*d^3*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2)
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^3//2, x) == :(-2*c*d*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)+4*c*b^-1*d^3*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-3//2+-4*b^-1*c^2*d^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^(1/2), x) == :(2*c*d*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2)
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-1//2, x) == :(2*c*d*b^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-3//2+-2*b^-1*c^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-3//2, x) == :(2*c*b^-1*d^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2+(1/2)*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//2*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+-1//4*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//4*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-5//2, x) == :(2*c*b^-1*d^-1*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-3//2+-3*b^-1*c^2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^9//2, x) == :(-20//21*c*b^-1*d^3*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^3//2+-2//7*c*d*b^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^7//2+40//21*c*b^-1*d^5*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-1//2+40//21*b^-1*c^2*d^4*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^7//2, x) == :(-64//15*b^-1*c^3*d^3*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)+-2//5*c*d*b^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^5//2+16//15*c*b^-1*d^3*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^(1/2))
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^5//2, x) == :(-2//3*c*d*b^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^3//2+4//3*c*b^-1*d^3*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-1//2+4//3*b^-1*c^2*d^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^3//2, x) == :(-8//3*d*b^-1*c^3*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)+2//3*c*d*b^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^(1/2))
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^(1/2), x) == :(2//3*c*d*b^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-1//2+2//3*b^-1*c^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^-1//2, x) == :(2//3*c*d*b^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-3//2)
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^-3//2, x) == :(2//3*c*b^-1*d^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-1//2+-1//3*b^-1*c^2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^5//2*(d*csc(a+b*x))^-5//2, x) == :(2//3*c*b^-1*d^-1*(c*sec(a+b*x))^3//2*(d*csc(a+b*x))^-3//2+(1/2)*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//2*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+-1//4*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//4*2^(1/2)*b^-1*c^2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^9//2, x) == :(-8//21*c*b^-1*d^3*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^3//2+-2//7*c*d*b^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^7//2)
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^7//2, x) == :(-4//5*c*b^-1*d^3*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^(1/2)+-2//5*c*d*b^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^5//2+-4//5*b^-1*d^4*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^5//2, x) == :(-2//3*c*d*b^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^3//2)
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^3//2, x) == :(-2*c*d*b^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^(1/2)+-2*b^-1*d^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2), x) == :((1/2)*2^(1/2)*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+-1//2*2^(1/2)*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//4*2^(1/2)*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//4*2^(1/2)*b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2, x) == :(b^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-3//2, x) == :(-1//2*c*b^-1*d^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-1//2+-1//8*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//16*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//8*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+1//16*2^(1/2)*b^-1*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-5//2, x) == :(-1//3*c*b^-1*d^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-3//2+(1/2)*b^-1*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^11//2, x) == :(-2//9*d*b^-1*c^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^9//2+2//45*b^-1*c^-1*d^3*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^5//2+8//45*b^-1*c^-1*d^5*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^9//2, x) == :(-2//7*d*b^-1*c^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^7//2+2//21*b^-1*c^-1*d^3*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^3//2+-2//21*b^-1*c^-2*d^4*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^7//2, x) == :(-2//5*c*d*b^-1*(c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^5//2)
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^5//2, x) == :(-2//3*d*b^-1*c^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^3//2+-1//3*b^-1*c^-2*d^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^3//2, x) == :(-2*d*b^-1*c^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)+(1/2)*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//2*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+-1//4*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//4*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^(1/2), x) == :(d*b^-1*c^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2+(1/2)*b^-1*c^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-1//2, x) == :((1/2)*d*b^-1*c^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-3//2+-1//8*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//16*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//8*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+1//16*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-3//2, x) == :(-1//3*c*b^-1*d^-1*(c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^-1//2+1//6*b^-1*c^-1*d^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2+1//12*b^-1*c^-2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^(1/2)*sin(2a+2*b*x)^(1/2)*Elliptic.F(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-5//2, x) == :(-1//4*c*b^-1*d^-1*(c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^-3//2+3//16*b^-1*c^-1*d^-1*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-3//2+-3//64*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-3//128*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+3//64*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+3//128*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^(1/2)*(d*csc(a+b*x))^-1//2*tan(a+b*x)^-1//2*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^9//2, x) == :(-2//7*c*d*b^-1*(c*sec(a+b*x))^-7//2*(d*csc(a+b*x))^7//2)
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^7//2, x) == :(-2//5*d*b^-1*c^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^5//2+6//5*b^-1*c^-1*d^3*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^(1/2)+6//5*b^-1*c^-2*d^4*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^5//2, x) == :(-2//3*d*b^-1*c^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^3//2+(1/2)*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-1//2*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+-1//4*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+1//4*2^(1/2)*b^-1*c^-2*d^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^3//2, x) == :(-2*d*b^-1*c^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^(1/2)+-3*b^-1*c^-2*d^2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^(1/2), x) == :((1/2)*d*b^-1*c^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-1//2+-3//8*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-3//16*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+3//8*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+3//16*2^(1/2)*b^-1*c^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^-1//2, x) == :(1//3*d*b^-1*c^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-3//2+(1/2)*b^-1*c^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^-3//2, x) == :(-1//4*c*b^-1*d^-1*(c*sec(a+b*x))^-7//2*(d*csc(a+b*x))^-1//2+1//16*b^-1*c^-1*d^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-1//2+-3//64*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-3//128*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+3//64*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+3//128*2^(1/2)*b^-1*c^-2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^-5//2, x) == :(-1//5*c*b^-1*d^-1*(c*sec(a+b*x))^-7//2*(d*csc(a+b*x))^-3//2+1//10*b^-1*c^-1*d^-1*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-3//2+3//20*b^-1*c^-2*d^-2*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^-1//2*sin(2a+2*b*x)^-1//2*Elliptic.E(a+-1//4*pi+b*x,2))
@test integrate((c*sec(a+b*x))^-5//2*(d*csc(a+b*x))^-7//2, x) == :(-5//48*c*b^-1*d^-3*(c*sec(a+b*x))^-7//2*(d*csc(a+b*x))^-1//2+-1//6*c*b^-1*d^-1*(c*sec(a+b*x))^-7//2*(d*csc(a+b*x))^-5//2+5//192*b^-1*c^-1*d^-3*(c*sec(a+b*x))^-3//2*(d*csc(a+b*x))^-1//2+-5//256*2^(1/2)*b^-1*c^-2*d^-4*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+-1*2^(1/2)*tan(a+b*x)^(1/2))+-5//512*2^(1/2)*b^-1*c^-2*d^-4*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+-1*2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x))+5//256*2^(1/2)*b^-1*c^-2*d^-4*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*arctan(1+2^(1/2)*tan(a+b*x)^(1/2))+5//512*2^(1/2)*b^-1*c^-2*d^-4*(c*sec(a+b*x))^-1//2*(d*csc(a+b*x))^(1/2)*tan(a+b*x)^(1/2)*log(1+2^(1/2)*tan(a+b*x)^(1/2)+tan(a+b*x)))
@test integrate((b*csc(e+f*x))^n*cos(e+f*x), x) == :(b*f^-1*(b*csc(e+f*x))^(-1+n)*(1+-1n)^-1)
@test integrate((b*csc(e+f*x))^n*cos(e+f*x)^3, x) == :(b*f^-1*(b*csc(e+f*x))^(-1+n)*(1+-1n)^-1+-1*b^3*f^-1*(b*csc(e+f*x))^(-3+n)*(3+-1n)^-1)
@test integrate((b*csc(e+f*x))^n*cos(e+f*x)^5, x) == :(b*f^-1*(b*csc(e+f*x))^(-1+n)*(1+-1n)^-1+b^5*f^-1*(b*csc(e+f*x))^(-5+n)*(5+-1n)^-1+-2*b^3*f^-1*(b*csc(e+f*x))^(-3+n)*(3+-1n)^-1)

max_stars_repo_path: test/4 Trig functions/4.5 Secant/4.5.0 (a sec)^m (b trg)^n.jl


text:
mutable struct ExtendedJumpArray{T,T2} <: AbstractArray{Float64,1}
  u::T
  jump_u::T2
end

Base.length(A::ExtendedJumpArray) = length(A.u) + length(A.jump_u)
Base.size(A::ExtendedJumpArray) = (length(A),)
@inline function Base.getindex(A::ExtendedJumpArray,i::Int)
  i <= length(A.u) ? A.u[i] : A.jump_u[i-length(A.u)]
end
@inline function Base.getindex(A::ExtendedJumpArray,I...)
  A[CartesianIndices(A.u,I...)]
end
@inline function Base.getindex(A::ExtendedJumpArray,I::CartesianIndex{1})
  A[I[1]]
end
@inline Base.setindex!(A::ExtendedJumpArray,v,I...) = (A[CartesianIndices(A.u,I...)] = v)
@inline Base.setindex!(A::ExtendedJumpArray,v,I::CartesianIndex{1}) = (A[I[1]] = v)
@inline function Base.setindex!(A::ExtendedJumpArray,v,i::Int)
  i <= length(A.u) ? (A.u[i] = v) : (A.jump_u[i-length(A.u)] = v)
end

Base.IndexStyle(::Type{<:ExtendedJumpArray}) = IndexLinear()
Base.similar(A::ExtendedJumpArray) = ExtendedJumpArray(similar(A.u),similar(A.jump_u))
Base.similar(A::ExtendedJumpArray,::Type{S}) where {S} = ExtendedJumpArray(similar(A.u,S),similar(A.jump_u,S))
Base.zero(A::ExtendedJumpArray) = fill!(similar(A),0)

# Required for non-diagonal noise
LinearAlgebra.mul!(c::ExtendedJumpArray,A::AbstractVecOrMat,u::AbstractVector) = mul!(c.u,A,u)

# Ignore axes
Base.similar(A::ExtendedJumpArray,::Type{S},axes::Tuple{Base.OneTo{Int}}) where {S} = ExtendedJumpArray(similar(A.u,S),similar(A.jump_u,S))

function recursivecopy!(dest::T, src::T) where T<:ExtendedJumpArray
  recursivecopy!(dest.u,src.u)
  recursivecopy!(dest.jump_u,src.jump_u)
end
Base.show(io::IO,A::ExtendedJumpArray) = show(io,A.u)
TreeViews.hastreeview(x::ExtendedJumpArray) = true
plot_indices(A::ExtendedJumpArray) = eachindex(A.u)

###### Broadcast overloading

const ExtendedJumpArrayStyle = Broadcast.ArrayStyle{ExtendedJumpArray}
Base.BroadcastStyle(::Type{<:ExtendedJumpArray}) = Broadcast.ArrayStyle{ExtendedJumpArray}()
Base.BroadcastStyle(::Broadcast.ArrayStyle{ExtendedJumpArray},::Broadcast.ArrayStyle) = Broadcast.ArrayStyle{ExtendedJumpArray}()
Base.BroadcastStyle(::Broadcast.ArrayStyle,::Broadcast.ArrayStyle{ExtendedJumpArray}) = Broadcast.ArrayStyle{ExtendedJumpArray}()
Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{ExtendedJumpArray}},::Type{ElType}) where ElType = similar(bc)

add_idxs1(A,x,expr) = expr
add_idxs1(A,::Type{T},expr) where {T<:ExtendedJumpArray} = :($(expr).u)
add_idxs1(A,::Type{T},expr) where {T<:AbstractArray} = :(@view($(expr)[1:length(A.u)]))

add_idxs2(A,x,expr) = expr
add_idxs2(A,::Type{T},expr) where {T<:ExtendedJumpArray} = :($(expr).jump_u)
add_idxs2(A,::Type{T},expr) where {T<:AbstractArray} = :(@view($(expr)[1:length(A.jump_u)]))

function Base.copy(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle})
    ret = Broadcast.flatten(bc)
    __broadcast(ret.f,ret.args...)
end

function Base.copyto!(dest::AbstractArray, bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle})
    ret = Broadcast.flatten(bc)
    __broadcast!(ret.f,dest,ret.args...)
end

@generated function __broadcast(f,A::ExtendedJumpArray,B...)
  exs1 = ((add_idxs1(A,B[i],:(B[$i])) for i in eachindex(B))...,)
  exs2 = ((add_idxs2(A,B[i],:(B[$i])) for i in eachindex(B))...,)
  res = quote
      ExtendedJumpArray(broadcast(f,A.u,$(exs1...)),broadcast(f,A.jump_u,$(exs2...)))
  end
  res
end

@generated function __broadcast!(f,A::ExtendedJumpArray,B...)
  exs1 = ((add_idxs1(A,B[i],:(B[$i])) for i in eachindex(B))...,)
  exs2 = ((add_idxs2(A,B[i],:(B[$i])) for i in eachindex(B))...,)
  res = quote
      broadcast!(f,A.u,$(exs1...));broadcast!(f,A.jump_u,$(exs2...))
      A
  end
  res
end