text:
module JungleHelperNodedCrumblePlatform

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/NodedCrumblePlatform" NodedCrumblePlatform(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, texture::String="default",
    nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const placements = Ahorn.PlacementDict(
    "Crumble Block (Noded, $(uppercasefirst(texture))) (Jungle Helper)" => Ahorn.EntityPlacement(
        NodedCrumblePlatform,
        "rectangle",
        Dict{String, Any}(
            "texture" => texture
        )
    ) for texture in Maple.crumble_block_textures
)

Ahorn.editingOptions(entity::NodedCrumblePlatform) = Dict{String, Any}(
    "texture" => Maple.crumble_block_textures
)

Ahorn.nodeLimits(entity::NodedCrumblePlatform) = 0, -1

Ahorn.minimumSize(entity::NodedCrumblePlatform) = 8, 0
Ahorn.resizable(entity::NodedCrumblePlatform) = true, false

function Ahorn.selection(entity::NodedCrumblePlatform)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    res = Ahorn.Rectangle[Ahorn.Rectangle(x, y, width, 8)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.Rectangle(nx, ny, width, 8))
    end

    return res
end

function renderCrumbleBlock(ctx::Ahorn.Cairo.CairoContext, x::Int, y::Int, width::Int, texture::String, alpha::Number=1.0)
    tilesWidth = div(width, 8)

    for i in 0:floor(Int, tilesWidth / 4)
        Ahorn.drawImage(ctx, texture, x + 32 * i, y, 0, 0, min(32, width - 32 * i), 8; alpha=alpha)
    end
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::NodedCrumblePlatform)
    texture = get(entity.data, "texture", "default")
    texture = "objects/crumbleBlock/$texture"

    x, y = Ahorn.position(entity)
    px, py = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)
        Ahorn.drawArrow(ctx, px + width / 2, py + 4, nx + width / 2, ny + 4, Ahorn.colors.selection_selected_fc, headLength=6)
        px, py = nx, ny
    end

    Ahorn.drawArrow(ctx, px + width / 2, py + 4, x + width / 2, y + 4, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::NodedCrumblePlatform, room::Maple.Room)
    texture = get(entity.data, "texture", "default")
    texture = "objects/crumbleBlock/$texture"

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))

    renderCrumbleBlock(ctx, x, y, width, texture)

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)
        renderCrumbleBlock(ctx, nx, ny, width, texture, 0.5)
    end
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
@testset "777.swap-adjacent-in-lr-string.jl" begin
    @test can_transform("RXXLRXRXL", "XRLXXRRLX") == true
    @test can_transform("XXXL", "LXXL") == false
end

ext: jl


text:
using ExprParsers

"""
    @insert_into_runhandlers outer_handler @syntax_eff begin
      # ...
    end

Next to simple Effects which can be directly composed down to plain values and reconstructed again from plain values,
there are also a couple of more elaborate Effects, which values cannot be extracted without providing further context.

Callables are a good example. It is impossible to extract the values of a callable without calling it, or without
constructing another Callable around it.

With all these example the typical flow would be like
```julia

Callable(function (args...; kwargs...)
  callablehandler = CallableHandler(args...; kwargs...)

  SomeOtherNeededContext() do info
    otherhandler = SomeOtherHandler(info)

    # ... possibly further nestings

    @runhandlers (callablehandler, otherhandler #= possible further handlers =#) @syntax_eff begin
      # ...
    end
  end
end)
```

`@inser_into_runhandlers` can be used to simplify and even separate these outer handlers from one another, so that
they can be used as composable interchangeable macros.

For example, here the implementation of `@runcallable` (ignoring macro hygiene)
```julia
macro runcallable(expr)
  :(Callable(function(args...; kwargs...)
    @insert_into_runhandlers(CallableHandler(args...; kwargs...), \$expr)
  end))
end
```
This will search for an existing call to `runhandlers` within the given `expr`, and if found, inserts the
CallableHandler similar to the motivating example above. If no `runhandlers` is found, it will create a new one.

This way, you can compose the nested code-structures very easily.
You only have to be careful, that you always run all outer effects at once, in one single statement,
so that `@insert_into_runhandlers` can actually find the right `runhandlers`.
"""
macro insert_into_runhandlers(handler, expr)
  expr = macroexpand(__module__, expr)
  esc(ExtensibleEffects.insert_into_runhandlers!(handler, expr))
end

insert_into_runhandlers!(handler, symbol::Symbol) = :(ExtensibleEffects.runhandlers($handler, $symbol))
function insert_into_runhandlers!(handler, expr::Expr)
  found_runhandlers = _insert_into_runhandlers!(handler, expr)
  found_runhandlers ? expr : :(ExtensibleEffects.runhandlers($handler, $expr))
end

"""
inserts handler into runhandlers expression and returns true if at least one runhandlers was found
"""
_insert_into_runhandlers!(handler, value) = false
function _insert_into_runhandlers!(handler, expr::Expr)
  if isrunhandlerscall(expr)
    handlerexpr = expr.args[2]  # first argument to macrocall/call
    if handlerexpr isa Expr && handlerexpr.head == :tuple
      insert!(handlerexpr.args, 1, handler)
    else
      handlerexpr_new = :(($handler, $handlerexpr))
      expr.args[2] = handlerexpr_new
    end
    true
  else
    any(expr.args) do subexpr
      _insert_into_runhandlers!(handler, subexpr)
    end
  end
end

const nesteddotparser = EP.NestedDot()
splitdots(value) = [value]
function splitdots(expr::Expr)
  parsed = parse_expr(nesteddotparser, expr)
  insert!(parsed.properties, 1, parsed.base)
end
isrunhandlerscall(value) = false
isrunhandlerscall(expr::Expr) = (
  (expr.head == :macrocall && splitdots(expr.args[1])[end] == Symbol("@runhandlers"))
  || (expr.head == :call && splitdots(expr.args[1])[end] == :runhandlers)
)

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
export write_mol, write_ff, write_mol_and_ff

"""
    write_mol(f, mol::Molecule)

Write the molecular topology of `mol` into `f`. `f` can be an I/O stream or a filename.
"""
function write_mol(io, mol::Molecule)
    println(io, mol.comment)
    println(io, length(mol.coords), " atoms")
    println(io, length(mol.bonds), " bonds")
    println(io, length(mol.angles), " angles")
    println(io, length(mol.diheds), " dihedrals")
    println(io, length(mol.improps), " impropers")

    print_coords(io, mol.coords)
    print_types(io, mol.types)
    print_masses(io, mol.masses)
    print_charges(io, mol.charges)

    !isempty(mol.bonds) && print_bonds(io, mol.bonds)
    !isempty(mol.angles) && print_angles(io, mol.angles)
    !isempty(mol.diheds) && print_dihedrals(io, mol.diheds)
    !isempty(mol.improps) && print_impropers(io, mol.improps)
end

function write_mol(fname::AbstractString, mol::Molecule)
    open(fname, "w") do io
        write_mol(io, mol)
    end
end

"""
    write_mol(f, ligpargenfile::AbstractString)

Read the molecular topology from file named `ligpargenfile` and write it into `f`. `f` can
    be an I/O stream or a filename.
"""
write_mol(f, mol::AbstractString) = write_mol(f, read_lpg_data(mol))

"""
    write_ff(f, mol::Molecule)

Write the forcefield parameters of `mol` into `f`. `f` can be an I/O stream or a filename.
"""
function write_ff(io, mol::Molecule)
    println(io,
    """
    pair_style        lj/cut/coul/long 12.0
    bond_style        harmonic
    angle_style       harmonic
    dihedral_style    opls
    improper_style    cvff
    special_bonds     lj/coul 0.0 0.0 0.5
    pair_modify       mix geometric
    """)

    for type in 1:length(mol.pair_coeffs)
        itype = findlast(==(type), mol.types)
        m = mol.masses[itype]
        join(io, ("mass          ", type, m, '\n'), ' ')
    end

    !isempty(mol.pair_coeffs) && println(io)

    for type in 1:length(mol.pair_coeffs)
        ε, σ = mol.pair_coeffs[type]
        join(io, ("pair_coeff    ", type, type, ε, σ, '\n'), ' ')
    end

    !isempty(mol.bond_coeffs) && println(io)

    for btype in 1:length(mol.bond_coeffs)
        k, r0 = mol.bond_coeffs[btype]
        join(io, ("bond_coeff    ", btype, k, r0, '\n'), ' ')
    end

    !isempty(mol.angle_coeffs) && println(io)

    for atype in 1:length(mol.angle_coeffs)
        k, θ0 = mol.angle_coeffs[atype]
        join(io, ("angle_coeff    ", atype, k, θ0, '\n'), ' ')
    end

    !isempty(mol.dihed_coeffs) && println(io)

    for dtype in 1:length(mol.dihed_coeffs)
        c1, c2, c3, c4 = mol.dihed_coeffs[dtype]
        join(io, ("dihedral_coeff    ", dtype, c1, c2, c3, c4, '\n'), ' ')
    end

    !isempty(mol.improper_coeffs) && println(io)

    for itype in 1:length(mol.improper_coeffs)
        k, d, n = mol.improper_coeffs[itype]
        join(io, ("improper_coeff    ", itype, k, d, n, '\n'), ' ')
    end
end

function write_ff(fname::AbstractString, mol::Molecule)
    open(fname, "w") do io
        write_ff(io, mol)
    end
end

"""
    write_ff(f, ligpargenfile::AbstractString)

Read the molecular topology from file named `ligpargenfile` and write its forcefield
    into `f`. `f` can be an I/O stream or a filename.
"""
write_ff(f, mol::AbstractString) = write_ff(f, read_lpg_data(mol))

"""
    write_mol_and_ff(fmask::AbstractString, mol::Molecule)

Write files "fmask.txt" and "fmask.ff" with molecular topology
and forcefield parameters of `mol`, respectively.
"""
function write_mol_and_ff(fmask::AbstractString, mol::Molecule)
    @sync begin
        @async write_ff("$fmask.ff", mol)
        @async write_mol("$fmask.txt", mol)
    end
    return
end

"""
    write_ff(f, ligpargenfile::AbstractString)

Read the molecular topology and forcefield parameters from `ligpargenfile` and write them
    into files "fmask.txt" and "fmask.ff", respectively.
"""
write_mol_and_ff(fmask, mol::AbstractString) = write_mol_and_ff(fmask, read_lpg_data(mol))

function print_coords(io, coords)
    println(io, "\nCoords\n")
    for i in 1:length(coords)
        x, y, z = coords[i]
        join(io, (i, x, y, z, '\n'), ' ')
    end
end

function print_types(io, types)
    println(io, "\nTypes\n")
    for i in 1:length(types)
        join(io, (i, types[i], '\n'), ' ')
    end
end

function print_charges(io, charges)
    println(io, "\nCharges\n")
    for i in 1:length(charges)
        join(io, (i, charges[i], '\n'), ' ')
    end
end

function print_masses(io, masses)
    println(io, "\nMasses\n")
    for i in 1:length(masses)
        join(io, (i, masses[i], '\n'), ' ')
    end
end

function print_bonds(io, bonds)
    println(io, "\nBonds\n")
    for i in 1:length(bonds)
        (a1, a2), type = bonds[i]
        join(io, (i, type, a1, a2, '\n'), ' ')
    end
end

function print_angles(io, angles)
    println(io, "\nAngles\n")
    for i in 1:length(angles)
        (a1, a2, a3), type = angles[i]
        join(io, (i, type, a1, a2, a3, '\n'), ' ')
    end
end

function print_dihedrals(io, dihed)
    println(io, "\nDihedrals\n")
    for i in 1:length(dihed)
        (a1, a2, a3, a4), type = dihed[i]
        join(io, (i, type, a1, a2, a3, a4, '\n'), ' ')
    end
end

function print_impropers(io, impropers)
    println(io, "\nImpropers\n")
    for i in 1:length(impropers)
        (a1, a2, a3, a4), type = impropers[i]
        join(io, (i, type, a1, a2, a3, a4, '\n'), ' ')
    end
end


text:
const STATS = Dict()

# Verifications

"""
    spy(f::Function)

Run the function without any patches. Using `spy` is more efficient than [`apply`](@ref)
due to the zero patch overhead.
"""
function spy(f::Function)
    reset_statistics()
    f()
end

"""
    called_exactly_once(f::Function, args...; kwargs...)

Return `true` if the function `f` with the provided `args` and `kwargs` was called
exactly once. This can be used after any executing a test script with mockable functions.
"""
called_exactly_once(f::Function, args...; kwargs...) = count_calls(f, args; kwargs...) == 1

"""
    called_at_least_once(f::Function, args...; kwargs...)

Return `true` if the function `f` with the provided `args` and `kwargs` was called
at least once. This can be used after any executing a test script with mockable functions.
"""
called_at_least_once(f::Function, args...; kwargs...) = count_calls(f, args; kwargs...) >= 1

"""
    was_not_called(f::Function, args...; kwargs...)

Return `true` if the function `f` with the provided `args` and `kwargs` was not called
at all. This can be used after any executing a test script with mockable functions.
"""
was_not_called(f::Function, args...; kwargs...)        = count_calls(f, args; kwargs...) == 0

# Spy recordings

materialize_kwargs(kwargs) = Any[p for p in kwargs]

function record_call(f::Function, args; kwargs...)
    kwargs = materialize_kwargs(kwargs)
    @debug "record_call" f args kwargs
    key = (f, args, kwargs)
    if haskey(STATS, key)
        STATS[key] += 1
    else
        STATS[key] = 1
    end
end

function count_calls(f::Function, args; kwargs...)
    kwargs = materialize_kwargs(kwargs)
    @debug "count_calls" f args kwargs
    key = (f, args, kwargs)
    return get(STATS, key, 0)
end

function reset_statistics()
    empty!(STATS)
end


text:
"""
Return the minimum distance between
two array of arrays.

Arguments
-----------
:Type{MinDistance}
    must be data type MinDistance

array_1:Array
    input 1

array_2:Array
    input 2
"""
function gage(::Type{MinDistance},array_1::Array, array_2::Array)
    if length(array_1) == 0
        @debug print_with_color(:red, "WARNING! input array 1 is empty for gage(MinDistance,...)\n\n")
        return 0.
    end

    if length(array_2) == 0
        @debug print_with_color(:red, "WARNING! input array 2 is empty for gage(MinDistance,...)\n\n")
        return 0.
    end

    if issubtype(typeof(array_1[1]),Number)
        array_1 = Array[array_1]
    end

    if issubtype(typeof(array_2[1]),Number)
        array_2 = Array[array_2]
    end

    return minimum(Distances.pairwise(Distances.Euclidean(), hcat(array_1...), hcat(array_2...)))
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/21/2014

#ACASX implementation: based on interfacing to the the Julia ADD

module ACASX_ADD_Impl

export
    addObserver,

    ACASX_ADD,
    ACASXInput,
    ACASXOutput,
    initialize,
    update

import Compat.ASCIIString

using AbstractCollisionAvoidanceSystemImpl
using AbstractCollisionAvoidanceSystemInterfaces
using CommonInterfaces

using AbstractCASCoordImpl

using ACASXCommonImpl
using CASCoordination
using CASInterface
using RLESUtils, Observers

type ACASX_ADD <: AbstractCollisionAvoidanceSystem
  my_id::Int64 #aircraft number
  max_intruders::Int64
  version::ASCIIString
  equipage::EQUIPAGE
  input::Input
  output::Output
  quant::Int64
  coord::AbstractCASCoord

  function ACASX_ADD(aircraft_id::Int64,
                 num_aircraft::Int, coord::AbstractCASCoord, equipage::EQUIPAGE=EQUIPAGE_TCAS)
    cas = new()
    cas.my_id = aircraft_id
    cas.max_intruders = num_aircraft - 1
    cas.version = "n/a"
    cas.equipage = equipage
    cas.coord = coord
    cas.input = Input(cas.max_intruders)
    cas.output = Output(cas.max_intruders)
    cas.quant = equipage == EQUIPAGE_TCAS ? 25 : 100
    setRecord(cas.coord, cas.my_id,
              ACASXCoordRecord(cas.my_id, equipage, cas.quant, cas.max_intruders))
    #reset(cas.casShared)
    return cas
  end
end

function update(cas::ACASX_ADD, input::ACASXInput)
  ACASXCommonImpl.update_from_coord!(input, cas.coord, cas.my_id)

  if cas.equipage == EQUIPAGE_TCAS
    #call custom update here...
    #my_custom_update()

    ACASXCommonImpl.update_to_coord!(cas.coord, cas.my_id, cas.output)
  end

  cas.input = input #keep a record
  return cas.output
end

function initialize(cas::ACASX_ADD)
  ACASXCommonImpl.initalize(cas)
end

end #module


text:
using Hecke
using  Hecke.GenericRound2
import Hecke.GenericRound2: radical_basis_power, radical_basis_trace, radical_basis_power_non_perfect
export FfOrdIdl

################################################################################
#
#  Constructors
#
################################################################################

mutable struct FfOrdIdl
  order::GenericRound2.Order
  basis::Vector{GenericRound2.OrderElem}
  basis_matrix::MatElem
  basis_mat_inv::FakeFracFldMat
  norm::RingElem
  minimum::RingElem
  is_prime::Int            # 0: don't know
                           # 1 known to be prime
                           # 2 known to be not prime
  iszero::Int              # as above
  is_principal::Int        # as above

  gen_one::RingElem
  gen_two::GenericRound2.OrderElem

  princ_gen::GenericRound2.OrderElem

  splitting_type::Tuple{Int, Int}
                         #ordered as ramification index, inertia degree

  function FfOrdIdl(O::GenericRound2.Order)
    r = new()
    r.order = O
    r.is_prime = 0
    r.iszero = 0
    r.is_principal = 0
    return r
  end

  function FfOrdIdl(O::GenericRound2.Order, M::MatElem)
    # create ideal of O with basis_matrix M
    r = FfOrdIdl(O)
    r.basis_matrix = M
    r.basis_mat_inv = FakeFracFldMat(pseudo_inv(M))
    return r
  end

  function FfOrdIdl(O::GenericRound2.Order, x::GenericRound2.OrderElem)
    # create ideal of O generated by x
    r = FfOrdIdl(O)
    r.princ_gen = x
    r.is_principal = 1

    if iszero(x)
       r.iszero = 1
    end

    r.norm = O.R(norm(O.F(x)))
    r.gen_one = r.norm
    r.gen_two = x
    return r
  end

  function FfOrdIdl(O::GenericRound2.Order, x::RingElem)
    return FfOrdIdl(O,O(x))
  end

  function FfOrdIdl(O::GenericRound2.Order, T:: Vector)
  V = hnf(vcat([representation_matrix(O(O.F(x))) for x in T]),:lowerleft)
  d = ncols(V)
  n = length(T)
  return FfOrdIdl(O, V[(n-1)*d+1:n*d,1:d])
  end

  function FfOrdIdl(O::GenericRound2.Order, p::RingElem, a::GenericRound2.OrderElem)
    #  create ideal of O generated by x
    r = FfOrdIdl(O)
    r.gen_one = p
    r.gen_two = a
    return r
  end

end

function AbstractAlgebra.zero(a::FfOrdIdl)
  O = a.order
  return FfOrdIdl(O,O(0))
end

AbstractAlgebra.iszero(I::FfOrdIdl) = (I.iszero == 1)

has_2_elem(A::FfOrdIdl) = isdefined(A, :gen_two)

################################################################################
#
#  Basic field access
#
################################################################################

@doc Markdown.doc"""
    order(x::FfOrdIdl) -> Order

Return the order, of which $x$ is an ideal.
"""
Hecke.order(a::FfOrdIdl) = a.order

###########################################################################################
#
#   Basis
#
###########################################################################################

@doc Markdown.doc"""
    has_basis(A::FfOrdIdl) -> Bool

Return whether $A$ has a basis already computed.
"""
has_basis(A::FfOrdIdl) = isdefined(A, :basis)

function assure_has_basis(A::FfOrdIdl)
  if isdefined(A, :basis)
    return nothing
  else
    assure_has_basis_matrix(A)
    O = order(A)
    M = A.basis_matrix
    Ob = basis(O)
    B = Vector{elem_type(O)}(undef, degree(O))
    y = O()
    for i in 1:degree(O)
      z = O()
      for k in 1:degree(O)
        mul!(y,O(M[i, k]), Ob[k])
        add!(z, z, y)
      end
      B[i] = z
    end
    A.basis = B
    return nothing
  end
end

@doc Markdown.doc"""
    basis(A::FfOrdIdl) -> Vector{OrderElem}

Return the basis of $A$.
"""
function Hecke.basis(A::FfOrdIdl; copy::Bool = true)
  assure_has_basis(A)
  if copy
    return deepcopy(A.basis)
  else
    return A.basis
  end
end

###########################################################################################
#
#   Basis Matrix
#
###########################################################################################

@doc Markdown.doc"""
    basis_matrix(A::FfOrdIdl) -> Mat

Return the basis matrix of $A$.
"""
function Hecke.basis_matrix(A::FfOrdIdl; copy::Bool = true)
  assure_has_basis_matrix(A)
  if copy
    return deepcopy(A.basis_matrix)
  else
    return A.basis_matrix
  end
end

function assure_has_basis_matrix(A::FfOrdIdl)
  if isdefined(A, :basis_matrix)
    return nothing
  end
  O = order(A)
  n = degree(O.F)

  if iszero(A)
    A.basis_matrix = zero_matrix(base_ring(O), n, n)
    return nothing
  end

  if has_princ_gen(A)
    A.basis_matrix = representation_matrix(A.princ_gen)
    return nothing
  end

  @hassert :NfOrd 1 has_2_elem(A)

  V = hnf(vcat([representation_matrix(x) for x in [O(A.gen_one),A.gen_two]]),:lowerleft)
  d = ncols(V)
  A.basis_matrix = V[d+1:2*d,1:d]
  return nothing
end

################################################################################
#
#  Basis matrix inverse
#
################################################################################

@doc Markdown.doc"""
    has_basis_mat_inv(A::FfOrdIdl) -> Bool

Return whether $A$ knows its inverse basis matrix.
"""
has_basis_mat_inv(A::FfOrdIdl) = isdefined(A, :basis_mat_inv)

@doc Markdown.doc"""
    basis_mat_inv(A::FfOrdIdl) -> FakeFracFldMat

Return the inverse of the basis matrix of $A$.
"""
function Hecke.basis_mat_inv(A::FfOrdIdl; copy::Bool = true) where T
  assure_has_basis_mat_inv(A)
  if copy
    return deepcopy(A.basis_mat_inv)
  else
    return A.basis_mat_inv
  end
end

function assure_has_basis_mat_inv(A::FfOrdIdl)
  if isdefined(A, :basis_mat_inv)
    return nothing
  else
    A.basis_mat_inv = FakeFracFldMat(pseudo_inv(basis_matrix(A)))
    return nothing
  end
end

###########################################################################################
#
#   Misc
#
###########################################################################################

(O::GenericRound2.Order)(p::PolyElem) = O(O.F(p))
Hecke.iscommutative(O::GenericRound2.Order) = true

Nemo.elem_type(::Type{GenericRound2.Order}) = GenericRound2.OrderElem

function Hecke.hnf(x::T, shape::Symbol =:upperright) where {T <: MatElem}
  if shape == :lowerleft
    h = hnf(reverse_cols(x))
    reverse_cols!(h)
    reverse_rows!(h)
    return h::T
  end
  return hnf(x)::T
end

################################################################################
#
#  Binary Operations
#
################################################################################

function Base.:(+)(a::FfOrdIdl, b::FfOrdIdl)
  check_parent(a, b)
  O = order(a)

  if iszero(a)
    return b
  end
  if iszero(b)
    return a
  end

  V = hnf(vcat(basis_matrix(a),basis_matrix(b)),:lowerleft)
  d = ncols(V)
  return FfOrdIdl(a.order, V[d+1:2*d,1:d])
end

Base.:(==)(a::FfOrdIdl, b::FfOrdIdl) = hnf(basis_matrix(a),:lowerleft) == hnf(basis_matrix(b),:lowerleft)
Base.isequal(a::FfOrdIdl, b::FfOrdIdl) = a == b

function Base.:(*)(a::FfOrdIdl, b::FfOrdIdl)
  O = order(a)
  Ma = basis_matrix(a)
  Mb = basis_matrix(b)
  V = hnf(vcat([Mb*representation_matrix(O([Ma[i,o] for o in 1:ncols(Ma)])) for i in 1:ncols(Ma)]),:lowerleft)
  d = ncols(V)
  return FfOrdIdl(O, V[d*(d-1)+1:d^2,1:d])
end

@doc Markdown.doc"""
    intersect(x::FfOrdIdl, y::FfOrdIdl) -> FfOrdIdl

Returns $x \cap y$.
"""
#TODO: Check for new hnf
function Base.intersect(a::FfOrdIdl, b::FfOrdIdl)
  M1 = hcat(basis_matrix(a), basis_matrix(a))
  d = nrows(M1)
  M2 = hcat(basis_matrix(b), zero_matrix(M1.base_ring,d,d))
  M = vcat(M1, M2)
  H = sub(hnf(M), d+1:2*d, d+1:2*d)
  return FfOrdIdl(a.order, H)
end

################################################################################
#
#  Powering
#
################################################################################

function Base.:(^)(A::FfOrdIdl, e::Int)
  O = order(A)
  if e == 0
    return FfOrdIdl(O, one(O))
  elseif e == 1
    return A
  end
    return Base.power_by_squaring(A, e)
end

################################################################################
#
#  Ad hoc multiplication
#
################################################################################

function Base.:*(x::GenericRound2.OrderElem, y::FfOrdIdl)
  parent(x) !== order(y) && error("Orders of element and ideal must be equal")
  return FfOrdIdl(parent(x), x) * y
end

Base.:*(x::FfOrdIdl, y::GenericRound2.OrderElem) = y * x

function Hecke.colon(a::FfOrdIdl, b::FfOrdIdl)

  O = order(a)
  n = degree(O)
  if isdefined(b, :gens)
    B = b.gens
  else
    B = basis(b)
  end

  bmatinv = basis_mat_inv(a, copy = false)

  n = representation_matrix(B[1])*bmatinv
  m = numerator(n)
  d = denominator(n)
  for i in 2:length(B)
    n = representation_matrix(B[i])*bmatinv
    mm = n.num
    dd = n.den
    l = lcm(dd, d)
    if l == d && l == dd
      m = hcat(m, mm)
    elseif l == d
      m = hcat(m, div(d, dd)*mm)
    elseif l == dd
      m = hcat(div(dd, d)*m, mm)
      d = dd
    else
      m = hcat(m*div(l, d), mm*div(l, dd))
      d = l
    end
  end
  m = transpose(m)
  m = hnf(m)
  # m is upper right HNF
  m = transpose(sub(m, 1:degree(O), 1:degree(O)))
  b = inv(FakeFracFldMat(m, d))
  return FfOrdFracIdl(O, b)
end

################################################################################
#
#  Exact Division
#
################################################################################

function Hecke.divexact(A::FfOrdIdl, b::RingElem)
  if iszero(A)
    return A
  end
  O = order(A)
  b = Hecke.AbstractAlgebra.MPolyFactor.make_monic(b)

  B = FfOrdIdl(O, divexact(basis_matrix(A), b))
  if false && has_basis_mat_inv(A)
    error("not defined at all")
    B.basis_mat_inv = b*A.basis_mat_inv
  end
  if has_princ_gen(A)
    B.princ_gen = O(divexact(O.F(A.princ_gen), b))
  end

  return B
end

function has_minimum(A::FfOrdIdl)
  return isdefined(A, :minimum)
end

function Hecke.minimum(A::FfOrdIdl; copy::Bool = true)
  assure_has_minimum(A)
  if copy
    return deepcopy(A.minimum)
  else
    return A.minimum
  end
end

function assure_has_minimum(A::FfOrdIdl)
  if has_minimum(A)
    return nothing
  end

  O = order(A)
  M = basis_matrix(A, copy = false)
  d = prod([M[i, i] for i = 1:nrows(M)])
  v = change_base_ring(O.R, coordinates(O(d)))
  fl, s = can_solve_with_solution(M, v, side = :left)
  @assert fl
  den = denominator(s[1]//d)
  for i = 2:ncols(s)
    den = lcm(den, denominator(s[i]//d))
  end
  A.minimum = Hecke.AbstractAlgebra.MPolyFactor.make_monic(den)
  return nothing
end

################################################################################
#
#  Norm
#
################################################################################

@doc Markdown.doc"""
    has_norm(A::FfOrdIdl) -> Bool

Return whether $A$ knows its norm.
"""
function has_norm(A::FfOrdIdl)
  return isdefined(A, :norm)
end

function assure_has_norm(A::FfOrdIdl)
  if has_norm(A)
    return nothing
  end

  if isdefined(A, :basis_matrix)
    A.norm = Hecke.AbstractAlgebra.MPolyFactor.make_monic(det(basis_matrix(A)))
    return nothing
  end

  if has_princ_gen(A)
    A.norm = Hecke.AbstractAlgebra.MPolyFactor.make_monic(norm(A.princ_gen))
    return nothing
  end

  assure_has_basis_matrix(A)
  A.norm = Hecke.AbstractAlgebra.MPolyFactor.make_monic(det(basis_matrix(A)))
  return nothing
end

@doc Markdown.doc"""
    norm(A::FfOrdIdl) -> RingElem

Return the norm of $A$, that is, the cardinality of $\mathcal O/A$, where
$\mathcal O$ is the order of $A$.
"""
function Hecke.norm(A::FfOrdIdl; copy::Bool = true)
  assure_has_norm(A)
  if copy
    return deepcopy(A.norm)
  else
    return A.norm
  end
end

################################################################################
#
#  Numerator and Denominator
#
################################################################################

function Hecke.numerator(a::Generic.FunctionFieldElem, O::GenericRound2.Order)
  return integral_split(a,O)[1]
end

function Hecke.denominator(a::Generic.FunctionFieldElem, O::GenericRound2.Order)
  return integral_split(a,O)[2]
end

################################################################################
#
#  Principal Generator
#
################################################################################

@doc Markdown.doc"""
    has_princ_gen(A::FfOrdIdl) -> Bool

Returns whether $A$ knows if it is generated by one element.
"""
function has_princ_gen(A::FfOrdIdl)
  return isdefined(A, :princ_gen)
end

################################################################################
#
#  Reduction of element modulo ideal
#
################################################################################

function Hecke.mod(x::GenericRound2.OrderElem, y::FfOrdIdl)
  parent(x) !== order(y) && error("Orders of element and ideal must be equal")
  # this function assumes that HNF is lower left
  # !!! This must be changed as soon as HNF has a different shape

  O = order(y)
  d = degree(O)
  a = change_base_ring(O.R,coordinates(x)).entries

  c = basis_matrix(y)
  t = O.R(0)
  for i in 1:d
    t = div(a[i], c[i,i])
    for j in 1:i
      a[j] = a[j] - t*c[i,j]
    end
  end
  z = O([a[i] for i in 1:lenth(a)])
  return z
end

################################################################################
#
#  Residue field
#
################################################################################

function Hecke.ResidueField(R::GFPPolyRing, p::gfp_poly)
  K, _ = FiniteField(p,"o")
  return K, MapFromFunc(x->K(x), y->R(y), R, K)
end

################################################################################
#
#  Factorization
#
################################################################################

function Hecke.index(O::GenericRound2.Order)
  index = O.R(1)
  if isdefined(O, :itrans)
    index = O.R(det(O.itrans))
  end
  return index
end

function prime_dec_nonindex(O::GenericRound2.Order, p::PolyElem, degree_limit::Int = 0, lower_limit::Int = 0)
  K = ResidueField(parent(p),p)[1]
  fact = factor(poly_to_residue(K,O.F.pol))
  result = []
  for (fac,e) in fact
    facnew = change_coefficient_ring(O.F.base_ring, fac)
    I = FfOrdIdl(O,p,O(facnew))
    I.is_prime = 1
    f = degree(fac)
    I.splitting_type = e, f
    I.norm = p^f
    I.minimum = p
    push!(result,(I,e))
  end
  return result
end

function poly_to_residue(K::AbstractAlgebra.Field, poly:: AbstractAlgebra.Generic.Poly{AbstractAlgebra.Generic.Rat{T}}) where T
  if poly == 0
    return K(0)
  else
    P, y = PolynomialRing(K,"y")
    coeffs = coefficients(poly)
    return sum([K(numerator(coeffs[i]))//K(denominator(coeffs[i]))*y^i for i in (0:length(poly)-1)])
  end
end

function Hecke.valuation(p::FfOrdIdl,A::FfOrdIdl)
  O = order(A)
  e = 0
  if has_2_elem(p)
    beta = Hecke.numerator(inv(O.F(p.gen_two)),O)
    newA = FfOrdFracIdl(beta*A,p.gen_one)
    while isintegral(newA)
      e += 1
      newA = FfOrdFracIdl(numerator(beta*newA),p.gen_one)
    end
  else
    newA = Hecke.colon(A,p)
    while isintegral(newA)
      e+=1
      newA = Hecke.colon(newA,FfOrdFracIdl(p))
    end
  end
  return e
end

function Hecke.factor(A::FfOrdIdl)
  O = A.order
  N = norm(A)
  factors = factor(N)
  primes = Dict{FfOrdIdl,Int}()
  for (f,e) in factors
    for (p,r) in prime_decomposition(O,f)
      primes[p] = valuation(p,A)
    end
  end
  return primes
end

function Hecke.prime_decomposition(O::GenericRound2.Order, p::RingElem, degree_limit::Int = degree(O), lower_limit::Int = 0; cached::Bool = true)
  if !(divides(index(O), p)[1])
    return prime_dec_nonindex(O, p, degree_limit, lower_limit)
  else
    return prime_dec_gen(O, p, degree_limit, lower_limit)
  end
end

function prime_dec_gen(O::GenericRound2.Order, p::RingElem, degree_limit::Int = degree(O), lower_limit::Int = 0)
  Ip = pradical(O, p)
  lp = _decomposition(O, FfOrdIdl(O, p), Ip, FfOrdIdl(O, one(O)), p)
  #=z = Tuple{ideal_type(O), Int}[]
  for (Q, e) in lp
    if degree(Q) <= degree_limit && degree(Q) >= lower_limit
      push!(z, (Q, e))
    end
  end
  return z
  =#
  return lp
end

function Hecke.pradical(O::GenericRound2.Order, p::RingElem)

  t = ResidueField(parent(p), p)

  if isa(t, Tuple)
    R, mR = t
  else
    R = t
    mR = MapFromFunc(x->R(x), y->lift(y), parent(p), R)
  end
#  @assert characteristic(F) == 0 || (isfinite(F) && characteristic(F) > degree(O))
  if characteristic(R) == 0 || characteristic(R) > degree(O)
    @vprint :NfOrd 1 "using trace-radical for $p\n"
    rad = radical_basis_trace
  elseif isa(R, Generic.RationalFunctionField)
    @vprint :NfOrd 1 "non-perfect case for radical for $p\n"
    rad = radical_basis_power_non_perfect
  else
    @vprint :NfOrd 1 "using radical-by-power for $p\n"
    rad = radical_basis_power
  end
  return FfOrdIdl(O,rad(O,p))
end

function _decomposition(O::GenericRound2.Order, I::FfOrdIdl, Ip::FfOrdIdl, T::FfOrdIdl, p::RingElem)
  #I is an ideal lying over p
  #T is contained in the product of all the prime ideals lying over p that do not appear in the factorization of I
  #Ip is the p-radical
  Ip1 = Ip + I
  A, OtoA = AlgAss(O, Ip1, p)
  AtoO = pseudo_inv(OtoA)
  ideals , AA = _from_algs_to_ideals(A, OtoA, AtoO, Ip1, p)
  for j in 1:length(ideals)
    P = ideals[j][1]
    f = P.splitting_type[2]
    e = valuation(P,FfOrdIdl(O,p))
    P.splitting_type = e, f
    ideals[j] = (P,e)
  end
  return ideals
end

function Hecke.AlgAss(O::GenericRound2.Order, I::FfOrdIdl, p::RingElem)
  @assert order(I) === O

  n = degree(O)
  bmatI = basis_matrix(I)

  basis_elts = Vector{Int}()
  for i = 1:n
    if iscoprime(bmatI[i, i], p)
      continue
    end

    push!(basis_elts, i)
  end

  r = length(basis_elts)
  FQ, phi = ResidueField(O.R,p)
  phi_inv = inv(phi)

  if r == 0
    A = _zero_algebra(FQ)

    local _image_zero

    let A = A
      function _image_zero(a::GenericRound2.OrderElem)
        return A()
      end
    end

    local _preimage_zero

    let O = O
      function _preimage_zero(a::AlgAssElem)
        return O()
      end
    end

    OtoA = Hecke.AbsOrdToAlgAssMor{typeof(O), elem_type(FQ)}(O, A, _image_zero, _preimage_zero)
    return A, OtoA
  end

  BO = basis(O)
  mult_table = Array{elem_type(FQ), 3}(undef, r, r, r)
  for i = 1:r
    M = representation_matrix(BO[basis_elts[i]])
    #TODO: write version of reduce rows mod hnf
    if r != degree(O)
      M = reduce_rows_mod_hnf!(M, bmatI, basis_elts)
    end
    for j = 1:r
      for k = 1:r
        mult_table[i, j, k] = FQ(M[basis_elts[j], basis_elts[k]])
      end
    end
  end

  if isone(BO[1])
    one = zeros(FQ, r)
    one[1] = FQ(1)
    A = AlgAss(FQ, mult_table, one)
  else
    A = AlgAss(FQ, mult_table)
  end
  if iscommutative(O)
    A.iscommutative = 1
  end

  local _image

  let I = I, A = A, basis_elts = basis_elts, FQ = FQ
    function _image(a::GenericRound2.OrderElem)
      c = coordinates(mod(a, I))
      return A([ FQ(numerator(c[i]))//FQ(denominator(c[i])) for i in basis_elts ])
    end
  end

  local _preimage

  let BO = BO, basis_elts = basis_elts, r = r
    function _preimage(a::AlgAssElem)
      return O(coefficients(a))
    end
  end

  OtoA = OrderToAlgAssMor{typeof(O), elem_type(FQ)}(O, A, _image, _preimage)

  return A, OtoA
end

# Reduces the rows of M in `rows` modulo N in place.
# Assumes that N is in lowerleft HNF.
function reduce_rows_mod_hnf!(M::MatElem, N::MatElem, rows::Vector{Int})
  for i in rows
    for j = ncols(M):-1:1
      if iszero(M[i, j])
        continue
      end

      t = div(M[i, j], N[j, j])
      for k = 1:j
        M[i, k] = M[i, k] - t*N[j, k]
      end
    end
  end
  return M
end

###############################################################################
#
#  Decomposition type using polygons
#
###############################################################################

function _from_algs_to_ideals(A::AlgAss{T}, OtoA::Map, AtoO::Map, Ip1, p::RingElem) where {T}

  O = order(Ip1)
  n = degree(O)
  R = O.R
  @vprint :NfOrd 1 "Splitting the algebra\n"
  AA = Hecke.decompose(A)
  @vprint :NfOrd 1 "Done \n"
  ideals = Vector{Tuple{typeof(Ip1), Int}}(undef, length(AA))
  N = basis_matrix(Ip1, copy = false)
  list_bases = Vector{Vector{Vector{elem_type(R)}}}(undef, length(AA))
  for i = 1:length(AA)
    l = Vector{Vector{elem_type(R)}}(undef, dim(AA[i][1]))
    for j = 1:length(l)
      v = change_base_ring(O.R,coordinates(AtoO(AA[i][2](AA[i][1][j]))))
      l[j] = [v[o] for o in 1:length(v)]
    end
    list_bases[i] = l
  end
  for i = 1:length(AA)
    B = AA[i][1]
    BtoA = AA[i][2]
    #we need the kernel of the projection map from A to B.
    #This is given by the basis of all the other components.
    f = dim(B)
    N1 = vcat(N, zero_matrix(R, dim(A) - f, n))
    t = 1
    for j = 1:length(AA)
      if j == i
        continue
      end
      for s = 1:length(list_bases[j])
        b = list_bases[j][s]
        for j = 1:degree(O)
          N1[n + t, j] = b[j]
        end
        t += 1
      end
    end
    N1 = view(hnf(N1, :lowerleft), nrows(N1) - degree(O) + 1:nrows(N1), 1:degree(O))
    P = FfOrdIdl(O, N1)
    P.minimum = p
    P.norm = p^f
    P.splitting_type = (0, f)
    P.is_prime = 1
    fromOtosimplealgebra = compose(OtoA, pseudo_inv(BtoA))
    #compute_residue_field_data!(P, fromOtosimplealgebra)
    ideals[i] = (P, 0)
  end
  return ideals, AA
end

################################################################################
#
#  OrderToAlgAssMor type
#
################################################################################

mutable struct OrderToAlgAssMor{S, T} <: Map{S, AlgAss{T}, Hecke.HeckeMap, OrderToAlgAssMor}
  header::Hecke.MapHeader

  function OrderToAlgAssMor{S, T}(O::S, A::AlgAss{T}, _image::Function, _preimage::Function) where {S <: GenericRound2.Order, T}
    z = new{S, T}()
    z.header = Hecke.MapHeader(O, A, _image, _preimage)
    return z
  end
end

function OrderToAlgAssMor(O::GenericRound2.Order, A::AlgAss{T}, _image, _preimage) where {T}
  return AbsOrdToAlgAssMor{typeof(O), T}(O, A, _image, _preimage)
end


text:
######################################################## model
F(z::BorderedArray,θ::AbstractVector) = F(z.u,(θ=θ,p=z.p))
function F(u::AbstractVector,parameters::NamedTuple)
	@unpack θ,p = parameters
	θ = Dict(

		# morphogens
		"c₆" => 10^p,
		"c₁₂" => 10^θ[1],

		# dissociation constants
		"R⁻P76" => 10^θ[2],
		"S⁻P76" => 10^θ[3],
		"R⁻P81" => 10^θ[4],
		"S⁻P81" => 10^θ[5],

		# crosstalk
		"S⁻C6" => 10^θ[6],
		"R⁻C12" => 10^θ[7],

		# leaky expression
		"ε₇₆" => 10^θ[8],
		"ε₈₁" => 10^θ[9],

		# hill coefficients
		"nᴿ" => 10^θ[10],
		"nˢ" => 10^θ[11],
		"nᴸ" => 10^θ[12],
		"nᵀ" => 10^θ[13],

		# synthesis rates
		"aᴿ" => 10^θ[14],
		"aˢ" => 10^θ[15],
		"aᴸ" => 10^θ[16],
		"aᵀ" => 10^θ[17],

		# degredations
		"dᴿ" => 10^θ[18],
		"dˢ" => 10^θ[19],
		"dᴸ" => 10^θ[20],
		"dᵀ" => 10^θ[21],
	)

	return F( sinh.(u), θ )
end

##############################  receiver binding
bound(x, cₓ, cₑ, ε, nₓ)  = x * x * ( abs(cₓ) ^ nₓ + abs(ε*cₑ) ^ nₓ )
hill(x, n)  = 1 / (1 + abs(x)^n)

##############################  promoter activations
function P76( u::AbstractVector, θ::Dict )
	@unpack R⁻C12,S⁻C6, nˢ,nᴿ, R⁻P76,S⁻P76,ε₇₆, c₆,c₁₂ = θ
	R,S,L,T = u

	activation = R⁻P76 * bound( R, c₆, c₁₂, R⁻C12, nᴿ ) + S⁻P76 * bound( S, c₁₂, c₆, S⁻C6, nˢ )
	return  (ε₇₆ + activation) / (1 + activation)
end

function P81( u::AbstractVector, θ::Dict )
	@unpack R⁻C12,S⁻C6, nˢ,nᴿ, R⁻P81,S⁻P81,ε₈₁, c₆,c₁₂ = θ
	R,S,L,T = u

	activation = R⁻P81 * bound( R, c₆, c₁₂, R⁻C12, nᴿ ) + S⁻P81 * bound( S, c₁₂, c₆, S⁻C6, nˢ )
	return (ε₈₁ + activation) / (1 + activation)
end

function F( u::AbstractVector, θ::Dict )
	@unpack aᴿ,aˢ,aᴸ,aᵀ, nᵀ,nᴸ, dᴿ,dˢ,dᴸ,dᵀ = θ

	f = first(u)*first(values(θ))
	F = similar(u,typeof(f))

	# calculate reaction rates
	F[1] = aᴿ * hill(u[4],nᵀ) - (1+dᴿ) * u[1]
	F[2] = aˢ * hill(u[3],nᴸ) - (1+dˢ) * u[2]

	F[3] = aᴸ * P76(u,θ) - (1+dᴸ) * u[3]
	F[4] = aᵀ * P81(u,θ) - (1+dᵀ) * u[4]

	return F
end

######################################################### targets and initial guess
X = StateSpace( 4, 0:0.01:3, [1,2] )
θ = SizedVector{21}(0.0,

	# dissociation constants
	-1.465,
	-3.782,
	-3.129,
	1.017,

	# crosstalk
	-8.0,
	-7.0,

	# leaky expression
	-1.943,
	-2.573,

	# hill coefficients
	-0.209,
	-0.489,
	-0.236,
	0.217,

	# synthesis rates
	1.222,
	0.456,
	3.394,
	1.828,

	# degredations
	1.390,
	0.933,
	0.0,
	0.0,
)


text:
# the result arrays are indexed by type Ti and have type T
struct WignerF{T<:Real, Ti} <: AbstractWignerF{T, Ti}
    j₂::T
    j₃::T
    m₂::T
    m₃::T
    nₘᵢₙ::Ti
    nₘₐₓ::Ti
end
function WignerF(::Type{T}, j₂::Ti, j₃::Ti, m₂::Ti, m₃::Ti) where {T<:Real, Ti}
    WignerF{T, Ti}(j₂, j₃, m₂, m₃, max(abs(j₂ - j₃), abs(m₂ + m₃)), j₂ + j₃)
end
function WignerF(::Type{T}, j₂, j₃, m₂, m₃) where {T}
    nₘᵢₙ = max(abs(j₂ - j₃), abs(m₂ + m₃))
    nₘₐₓ = j₂ + j₃
    # if isinteger(nₘᵢₙ) && isinteger(nₘₐₓ)
    #     return WignerF{T, Int}(j₂, j₃, m₂, m₃, nₘᵢₙ, nₘₐₓ)
    # else
    return WignerF{T, HalfInt}(j₂, j₃, m₂, m₃, nₘᵢₙ, nₘₐₓ)  # better to just be type stable
    # end
end
WignerF(j₂, j₃, m₂, m₃) = WignerF(Float64, j₂, j₃, m₂, m₃)

"""
    get_wigner_array(w::AbstractWigner{T}) where {T}

Utility function for getting an OffsetArray with indices from jₘᵢₙ to jₘₐₓ.

# Arguments
- `w::AbstractWigner{T}`: contains the quantum numbers and dispatches on the kind of symbol

# Returns
- `OffsetArray{T}`: an array for wigner symbols
"""
get_wigner_array(w::AbstractWigner{T,Ti}) where {T,Ti} = (
    WignerSymbolVector(zeros(T, length(w.nₘᵢₙ:w.nₘₐₓ)), w.nₘᵢₙ:w.nₘₐₓ))

"""
    rψ!(w::AbstractWigner{T}, n::Integer, iterates::AbstractVector{T}) where T

Backward recurrence scheme defined by equation 2 in Luscombe and Luban 1998. Iteratively
generates the ratio rψ(n) = ψ(n) / ψ(n+1) in the `iterates` vector.

# Arguments
- `w::AbstractWigner{T}`: contains the quantum numbers and dispatches on the kind of symbol
- `nmid::Integer`: current index of the recurrence
- `ψ::AbstractVector{T}`: store the values of rψ here during recursion.

# Returns
- `stop::Int`: the index the iteration stopped
"""
function rψ!(w::AbstractWigner{T,Ti}, nmid::Ti, ψ::AbstractVector{T}) where {T,Ti}
    for n in w.nₘₐₓ:-1:nmid
        if n == w.nₘₐₓ
            ψ[n] = -Zψ(w, n) / Yψ(w, n)
        else
            ψ[n] = -Zψ(w, n) / (Yψ(w, n) + Xψ(w, n) * ψ[n+1])
        end
        if !isfinite(ψ[n])
            ψ[n] = zero(T)
            return n
        end
        if abs(ψ[n]) ≥ 1.0
            return n
        end
    end
    return nmid
end

function sψ!(w::AbstractWigner{T,Ti}, nmid::Ti, ψ::AbstractVector{T}) where {T,Ti}
    for n in w.nₘᵢₙ:nmid
        if n == w.nₘᵢₙ
            ψ[n] = -Xψ(w, n) / Yψ(w, n)
        else
            ψ[n] = -Xψ(w, n) / (Yψ(w, n) + Zψ(w, n) * ψ[n-1])
        end
        if !isfinite(ψ[n])
            ψ[n] = zero(T)
            return n
        end
        if abs(ψ[n]) ≥ one(T)
            return n
        end
    end
    return nmid
end

"""
Three-term recurrence relation for the classical region.
"""
function ψauxplus!(w::AbstractWignerF{T, Ti},
                   n₋::Ti, nc::Ti, ψ::AbstractVector{T}) where {T, Ti}
    start_index = n₋
    if n₋ == w.nₘᵢₙ
        ψ[w.nₘᵢₙ] = one(T)  # set an arbitrary value for this
        if iszero(w.nₘᵢₙ)
            ψ[w.nₘᵢₙ+1] = -(w.m₃-w.m₂+2B(w,w.nₘᵢₙ)) / A(w,1)
        else
            ψ[w.nₘᵢₙ+1] = -(B(w,w.nₘᵢₙ)) / (w.nₘᵢₙ* A(w, w.nₘᵢₙ+1))
        end
        start_index = w.nₘᵢₙ+1
    end

    for n in start_index:(nc-1)
        Xn = Xψ(w, n)
        if iszero(Xn)
            ψ[n+1] = zero(T)
        else
            ψ[n+1] = -(Yψ(w, n) * ψ[n] + Zψ(w, n) * ψ[n - 1]) / Xn
        end
    end
end

function selection_rules(w::AbstractWignerF{T}) where T
    small = 10 * eps(T)
    return (
        (abs(w.m₂) < w.j₂ + small) &&
        (abs(w.m₃) < w.j₃ + small) &&
        isinteger(w.m₂ - w.j₂) &&
        isinteger(w.m₃ - w.j₃)
    )
end

A(w::AbstractWignerF, j) = sqrt((j^2 - (w.j₂ - w.j₃)^2) *
    ((w.j₂ + w.j₃ + 1)^2 - j^2) * (j^2 - (w.m₂ + w.m₃)^2))
B(w::AbstractWignerF, j) = (2 * j + 1) * ((w.m₂ + w.m₃) * (w.j₂ * (w.j₂ + 1) -
    w.j₃ * (w.j₃ + 1)) - (w.m₂ - w.m₃) * j * (j + 1))
Xψ(w::AbstractWignerF, j) = j * A(w, j+1)
Yψ(w::AbstractWignerF, j) = B(w, j)
Zψ(w::AbstractWignerF, j) = (j + 1) * A(w, j)
function normalization(w::AbstractWignerF{T}, ψ₀::AbstractVector{T}) where T
    norm = zero(T)
    for j in eachindex(ψ₀)
        norm += (2j + 1) * ψ₀[j]^2
    end
    return sqrt(abs(norm))
end
f_jmax_sgn(w::AbstractWignerF) = iseven(Int(w.j₂ - w.j₃ + w.m₂ + w.m₃)) ? 1 : -1

"""
    wigner3j_f(::Type{T}, j₂, j₃, m₂, m₃) where {T<:Real}

Computes all allowed j₁ given fixed j₂, j₃, m₂, m₃, m₁=-m₂-m₃.

# Arguments
- `T::Type{<:Real}`: output array type
- `j₂`: quantum number
- `j₃`: quantum number
- `m₂`: quantum number
- `m₃`: quantum number

# Returns
- `Tuple{Vector{Int}, Vector{T}}`: j₁ values and wigner symbols
"""
function wigner3j_f(::Type{T}, j₂, j₃, m₂, m₃) where {T<:Real}
    w = WignerF(T, j₂, j₃, m₂, m₃)  # establish the problem variables and method
    if !selection_rules(w)
        # throw(DomainError("Bad combination of inputs."))
        return OffsetArray(T[], 1:0)
    end
    w3j = get_wigner_array(w)  # construct an OffsetArray with the right indices and type
    wigner3j_f!(w, w3j)  # call the mutating function
    return w3j
end
wigner3j_f(j₂, j₃, m₂, m₃) = wigner3j_f(Float64, j₂, j₃, m₂, m₃)

function wigner3j_f!(w::AbstractWignerF{T,Ti}, w3j::AbstractVector{T}) where {T,Ti}

    # special case, if it's just one nontrivial result
    if length(w3j) == 1
        w3j[w.nₘᵢₙ] = f_jmax_sgn(w) / sqrt(w.nₘᵢₙ + w.j₂ + w.j₃+1)
        return
    end

    nmid = Ti( (w.nₘᵢₙ + w.nₘₐₓ) / 2 )
    fill!(parent(w3j), zero(T))

    # special case that performs an outwards classical solution if m₁ = m₂ = m₃ = 0
    # and skips the symbols with odd ∑jᵢ since those are zero.
    if iszero(w.m₂) && iszero(w.m₃) && isinteger(nmid)
        nmid = Int(floor((w.nₘᵢₙ + w.nₘₐₓ) / 2 ))
        nmid = iseven(nmid - w.nₘᵢₙ) ? nmid : nmid + 1  # ensure start index is even
        if w.nₘᵢₙ < nmid < w.nₘₐₓ  # check to make sure it's not the ends
            return classical_wigner3j_m0!(w, w3j)
        end
    end

    # otherwise attempt the general solution
    nmid = Ti(ceil((w.nₘᵢₙ + w.nₘₐₓ)/2))
    # attempt the non-classical two term nonlinear recurrences
    n₊ = rψ!(w, nmid, w3j)
    n₋ = sψ!(w, nmid, w3j)
    for k in (n₊+1):w.nₘₐₓ  # product the ratios upwards
        w3j[k] *= w3j[k-1]
    end
    for k in (n₋-1):-1:w.nₘᵢₙ  # product the ratios downwards
        w3j[k] *= w3j[k+1]
    end
    # now use the three-term classical recurrence in the middle
    ψn₋ = w3j[n₋]
    ψn₊ = w3j[n₊]
    ψauxplus!(w, n₋, n₊, w3j)
    # rescale the three parts
    if n₋ > w.nₘᵢₙ
        scale_middle = ψn₋ / w3j[n₋]
        for j in w.nₘᵢₙ:(n₋+1)
            w3j[j] *= scale_middle
        end
    end
    if n₊ < w.nₘₐₓ
        scale_end = w3j[n₊] / ψn₊
        for j in (n₊+1):w.nₘₐₓ
            w3j[j] *= scale_end
        end
    end
    # normalize the results
    norm = one(T) / normalization(w, w3j)
    if sign(w3j[w.nₘₐₓ]) != f_jmax_sgn(w)
        norm *= -1
    end
    symbols = parent(w3j)
    symbols .*= norm
end

"""
Special case iteration for mᵢ=0. In this case, j₁ must be an integer too, so the index
type is restricted to Int.
"""
function f_to_min_m0!(w::AbstractWignerF{T, Int},
                      nmid::Int, ψ::AbstractVector{T}) where {T}
    ψ[nmid] = one(T)
    ψ[nmid+1] = zero(T)
    for n in (nmid+1):2:(w.nₘₐₓ-1)
        Xn = Xψ(w, n)
        @assert abs(Xn) > zero(T)
        ψ[n+1] = -(Yψ(w, n) * ψ[n] + Zψ(w, n) * ψ[n - 1]) / Xn
    end
end

function f_to_max_m0!(w::AbstractWignerF{T, Int},
                      nmid::Int, ψ::AbstractVector{T}) where {T}
    ψ[nmid] = one(T)
    ψ[nmid-1] = zero(T)
    for n in (nmid-1):-2:(w.nₘᵢₙ+1)
        Zn = Zψ(w, n)
        @assert abs(Zn) > zero(T)
        ψ[n-1] = -(Yψ(w, n) * ψ[n] + Xψ(w, n) * ψ[n + 1]) / Zn
    end
end

"""
    classical_wigner3j_m0!(::Type{T}, j₂, j₃, m₂, m₃) where {T<:Real}

Computes all allowed j₁ given fixed j₂, j₃, m₁, m₂, m₃, subject to m₁ = m₂ = m₃ = 0. This
applies the classical three-term recurrence relation and iterates two at a time, since in
this special case all symbols with odd ∑jᵢ are zero. Unlike other Wigner symbols, this
special case requires iterating outwards, as one must recur towards increasing |fⱼ| for
stability.

# Arguments
- `T::Type{<:Real}`: output array type
- `j₂::Tn`: quantum number
- `j₃::Tn`: quantum number
- `m₂::Tn`: quantum number
"""
function classical_wigner3j_m0!(w::AbstractWignerF{T,Int}, w3j::AbstractVector{T}) where T
    nmid = Int(floor((w.nₘᵢₙ + w.nₘₐₓ) / 2 ))
    nmid = iseven(nmid - w.nₘᵢₙ) ? nmid : nmid + 1  # ensure start index is even
    f_to_min_m0!(w, nmid, w3j)
    f_to_max_m0!(w, nmid, w3j)

    norm = one(T) / normalization(w, w3j)
    if sign(w3j[w.nₘₐₓ]) != f_jmax_sgn(w)
        norm *= -1
    end
    symbols = parent(w3j)
    symbols .*= norm
end

## ----- todo: 6j

# struct WignerH{T} <: AbstractWignerH{T}
#     j₂::T
#     j₃::T
#     l₁::T
#     l₂::T
#     l₃::T
#     nₘᵢₙ::T
#     nₘₐₓ::T

#     WignerH(::Type{T}, j₂, j₃, l₁, l₂, l₃) where {T} = new{T}(
#         j₂, j₃, l₁, l₂, l₃, max(abs(j₂ - j₃), abs(l₂ - l₃)), min(j₂ + j₃, l₂ + l₃))
# end

# E(w::AbstractWignerH, j) = sqrt((j^2 - (w.j₂ - w.j₃)^2) * ((w.j₂ + w.j₃ + 1)^2 - j^2) *
#     (j^2 - (w.l₂ - w.l₃)^2) * ((w.l₂ + w.l₃ + 1)^2 - j^2))
# F(w::AbstractWignerH, j) = (2j+1) * (
#     j * (j + 1)  * (-j * (j+1) + w.j₂ * (w.j₂ + 1) + w.j₃ * (w.j₃ + 1) - 2 * w.l₁ * (w.l₁ + 1)) +
#     w.l₂ * (w.l₂ + 1) * (j * (j+1) + w.j₂ * (w.j₂ + 1) - w.j₃ * (w.j₃ + 1)) +
#     w.l₃ * (w.l₃ + 1) * (j * (j+1) - w.j₂ * (w.j₂ + 1) + w.j₃ * (w.j₃ + 1))
# )
# Xψ(w::AbstractWignerH, j) = j * E(w, j+1)
# Yψ(w::AbstractWignerH, j) = F(w, j)
# Zψ(w::AbstractWignerH, j)= (j+1) * E(w, j)


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

@test integrate(x^5*(a+b*sec(c+d*x^2)), x) == :(1//6*a*x^6+b*d^-3*Polylogarithms.polylog(3,im*exp(im*(c+d*x^2)))+-1*b*d^-3*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^2)))+im*b*d^-2*x^2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^2)))+-1*im*b*d^-1*x^4*arctan(exp(im*(c+d*x^2)))+-1*im*b*d^-2*x^2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^2))))
@test integrate(x^3*(a+b*sec(c+d*x^2)), x) == :(1//4*a*x^4+(1/2)*im*b*d^-2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^2)))+-1//2*im*b*d^-2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^2)))+-1*im*b*d^-1*x^2*arctan(exp(im*(c+d*x^2))))
@test integrate(x*(a+b*sec(c+d*x^2)), x) == :((1/2)*a*x^2+(1/2)*b*d^-1*arctanh(sin(c+d*x^2)))
@test integrate(x^5*(a+b*sec(c+d*x^2))^2, x) == :(1//6*a^2*x^6+b^2*d^-2*x^2*log(1+exp(2*im*(c+d*x^2)))+(1/2)*b^2*d^-1*x^4*tan(c+d*x^2)+-2*a*b*d^-3*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^2)))+2*a*b*d^-3*Polylogarithms.polylog(3,im*exp(im*(c+d*x^2)))+-1//2*im*b^2*d^-1*x^4+-1//2*im*b^2*d^-3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^2)))+-2*im*a*b*d^-1*x^4*arctan(exp(im*(c+d*x^2)))+-2*im*a*b*d^-2*x^2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^2)))+2*im*a*b*d^-2*x^2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^2))))
@test integrate(x^3*(a+b*sec(c+d*x^2))^2, x) == :(1//4*a^2*x^4+(1/2)*b^2*d^-2*log(cos(c+d*x^2))+(1/2)*b^2*d^-1*x^2*tan(c+d*x^2)+im*a*b*d^-2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^2)))+-1*im*a*b*d^-2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^2)))+-2*im*a*b*d^-1*x^2*arctan(exp(im*(c+d*x^2))))
@test integrate(x*(a+b*sec(c+d*x^2))^2, x) == :((1/2)*a^2*x^2+(1/2)*b^2*d^-1*tan(c+d*x^2)+a*b*d^-1*arctanh(sin(c+d*x^2)))
@test integrate(x*sec(a+b*x^2)^7, x) == :(5//32*b^-1*arctanh(sin(a+b*x^2))+1//12*b^-1*sec(a+b*x^2)^5*tan(a+b*x^2)+5//32*b^-1*sec(a+b*x^2)*tan(a+b*x^2)+5//48*b^-1*sec(a+b*x^2)^3*tan(a+b*x^2))
@test integrate(x^5*(a+b*sec(c+d*x^2))^-1, x) == :(1//6*a^-1*x^6+im*b*a^-1*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+b*a^-1*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1*im*b*a^-1*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1*b*a^-1*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+(1/2)*im*b*a^-1*d^-1*x^4*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1//2*im*b*a^-1*d^-1*x^4*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2))))
@test integrate(x^3*(a+b*sec(c+d*x^2))^-1, x) == :(1//4*a^-1*x^4+(1/2)*b*a^-1*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1//2*b*a^-1*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+(1/2)*im*b*a^-1*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1//2*im*b*a^-1*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2))))
@test integrate(x*(a+b*sec(c+d*x^2))^-1, x) == :((1/2)*a^-1*x^2+-1*b*a^-1*d^-1*(a+b)^-1//2*(a+-1b)^-1//2*arctanh((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*c+(1/2)*d*x^2)))
@test integrate(x^5*(a+b*sec(c+d*x^2))^-2, x) == :(1//6*a^-2*x^6+im*a^-2*b^3*d^-3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+a^-2*b^2*d^-2*x^2*(a^2+-1*b^2)^-1*log(1+a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^2)))+a^-2*b^2*d^-2*x^2*(a^2+-1*b^2)^-1*log(1+a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^2)))+a^-2*b^3*d^-2*x^2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1*im*a^-2*b^2*d^-3*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1*im*a^-2*b^2*d^-3*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1*im*a^-2*b^3*d^-3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1*a^-2*b^3*d^-2*x^2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-2*im*b*a^-2*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-2*b*a^-2*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+2*im*b*a^-2*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+2*b*a^-2*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1//2*im*a^-2*b^2*d^-1*x^4*(a^2+-1*b^2)^-1+im*b*a^-2*d^-1*x^4*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+(1/2)*im*a^-2*b^3*d^-1*x^4*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+(1/2)*a^-1*b^2*d^-1*x^4*(b+a*cos(c+d*x^2))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^2)+-1*im*b*a^-2*d^-1*x^4*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1//2*im*a^-2*b^3*d^-1*x^4*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2))))
@test integrate(x^3*(a+b*sec(c+d*x^2))^-2, x) == :(1//4*a^-2*x^4+b*a^-2*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+(1/2)*a^-2*b^2*d^-2*(a^2+-1*b^2)^-1*log(b+a*cos(c+d*x^2))+(1/2)*a^-2*b^3*d^-2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1*b*a^-2*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1//2*a^-2*b^3*d^-2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+im*b*a^-2*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+(1/2)*im*a^-2*b^3*d^-1*x^2*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+(1/2)*a^-1*b^2*d^-1*x^2*(b+a*cos(c+d*x^2))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^2)+-1*im*b*a^-2*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2)))+-1//2*im*a^-2*b^3*d^-1*x^2*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^2))))
@test integrate(x*(a+b*sec(c+d*x^2))^-2, x) == :((1/2)*a^-2*x^2+(1/2)*a^-1*b^2*d^-1*(a+b*sec(c+d*x^2))^-1*(a^2+-1*b^2)^-1*tan(c+d*x^2)+-1*b*a^-2*d^-1*(a+b)^-3//2*(a+-1b)^-3//2*(-1*b^2+2*a^2)*arctanh((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*c+(1/2)*d*x^2)))
@test integrate(x^3*(a+b*sec(c+d*x^(1/2))), x) == :(1//4*a*x^4+-10080*im*b*d^-8*Polylogarithms.polylog(8,-1*im*exp(im*(c+d*x^(1/2))))+-10080*b*d^-7*x^(1/2)*Polylogarithms.polylog(7,-1*im*exp(im*(c+d*x^(1/2))))+-1680*b*d^-5*x^3//2*Polylogarithms.polylog(5,im*exp(im*(c+d*x^(1/2))))+-84*b*d^-3*x^5//2*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+84*b*d^-3*x^5//2*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+1680*b*d^-5*x^3//2*Polylogarithms.polylog(5,-1*im*exp(im*(c+d*x^(1/2))))+10080*im*b*d^-8*Polylogarithms.polylog(8,im*exp(im*(c+d*x^(1/2))))+10080*b*d^-7*x^(1/2)*Polylogarithms.polylog(7,im*exp(im*(c+d*x^(1/2))))+-5040*im*b*x*d^-6*Polylogarithms.polylog(6,im*exp(im*(c+d*x^(1/2))))+-420*im*b*d^-4*x^2*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-14*im*b*d^-2*x^3*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-4*im*b*d^-1*x^7//2*arctan(exp(im*(c+d*x^(1/2))))+14*im*b*d^-2*x^3*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2))))+420*im*b*d^-4*x^2*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2))))+5040*im*b*x*d^-6*Polylogarithms.polylog(6,-1*im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*sec(c+d*x^(1/2))), x) == :(1//3*a*x^3+-240*im*b*d^-6*Polylogarithms.polylog(6,im*exp(im*(c+d*x^(1/2))))+-240*b*d^-5*x^(1/2)*Polylogarithms.polylog(5,im*exp(im*(c+d*x^(1/2))))+-40*b*d^-3*x^3//2*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+40*b*d^-3*x^3//2*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+240*im*b*d^-6*Polylogarithms.polylog(6,-1*im*exp(im*(c+d*x^(1/2))))+240*b*d^-5*x^(1/2)*Polylogarithms.polylog(5,-1*im*exp(im*(c+d*x^(1/2))))+-120*im*b*x*d^-4*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-10*im*b*d^-2*x^2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-4*im*b*d^-1*x^5//2*arctan(exp(im*(c+d*x^(1/2))))+10*im*b*d^-2*x^2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2))))+120*im*b*x*d^-4*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*sec(c+d*x^(1/2))), x) == :((1/2)*a*x^2+-12*im*b*d^-4*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-12*b*d^-3*x^(1/2)*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+12*im*b*d^-4*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2))))+12*b*d^-3*x^(1/2)*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+-6*im*b*x*d^-2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-4*im*b*d^-1*x^3//2*arctan(exp(im*(c+d*x^(1/2))))+6*im*b*x*d^-2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^3*(a+b*sec(c+d*x^(1/2)))^2, x) == :(1//4*a^2*x^4+315//2*b^2*d^-8*Polylogarithms.polylog(7,-1*exp(2*im*(c+d*x^(1/2))))+-315*x*b^2*d^-6*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^(1/2))))+-2*im*b^2*d^-1*x^7//2+2*b^2*d^-1*x^7//2*tan(c+d*x^(1/2))+14*b^2*d^-2*x^3*log(1+exp(2*im*(c+d*x^(1/2))))+105*b^2*d^-4*x^2*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-20160*im*a*b*d^-8*Polylogarithms.polylog(8,-1*im*exp(im*(c+d*x^(1/2))))+-20160*a*b*d^-7*x^(1/2)*Polylogarithms.polylog(7,-1*im*exp(im*(c+d*x^(1/2))))+-3360*a*b*d^-5*x^3//2*Polylogarithms.polylog(5,im*exp(im*(c+d*x^(1/2))))+-315*im*b^2*d^-7*x^(1/2)*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^(1/2))))+-168*a*b*d^-3*x^5//2*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+-42*im*b^2*d^-3*x^5//2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+168*a*b*d^-3*x^5//2*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+210*im*b^2*d^-5*x^3//2*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+3360*a*b*d^-5*x^3//2*Polylogarithms.polylog(5,-1*im*exp(im*(c+d*x^(1/2))))+20160*im*a*b*d^-8*Polylogarithms.polylog(8,im*exp(im*(c+d*x^(1/2))))+20160*a*b*d^-7*x^(1/2)*Polylogarithms.polylog(7,im*exp(im*(c+d*x^(1/2))))+-10080*im*a*b*x*d^-6*Polylogarithms.polylog(6,im*exp(im*(c+d*x^(1/2))))+-840*im*a*b*d^-4*x^2*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-28*im*a*b*d^-2*x^3*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-8*im*a*b*d^-1*x^7//2*arctan(exp(im*(c+d*x^(1/2))))+28*im*a*b*d^-2*x^3*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2))))+840*im*a*b*d^-4*x^2*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2))))+10080*im*a*b*x*d^-6*Polylogarithms.polylog(6,-1*im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*sec(c+d*x^(1/2)))^2, x) == :(1//3*a^2*x^3+-15*b^2*d^-6*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^(1/2))))+-2*im*b^2*d^-1*x^5//2+2*b^2*d^-1*x^5//2*tan(c+d*x^(1/2))+10*b^2*d^-2*x^2*log(1+exp(2*im*(c+d*x^(1/2))))+30*x*b^2*d^-4*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-480*im*a*b*d^-6*Polylogarithms.polylog(6,im*exp(im*(c+d*x^(1/2))))+-480*a*b*d^-5*x^(1/2)*Polylogarithms.polylog(5,im*exp(im*(c+d*x^(1/2))))+-80*a*b*d^-3*x^3//2*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+-20*im*b^2*d^-3*x^3//2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+30*im*b^2*d^-5*x^(1/2)*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+80*a*b*d^-3*x^3//2*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+480*im*a*b*d^-6*Polylogarithms.polylog(6,-1*im*exp(im*(c+d*x^(1/2))))+480*a*b*d^-5*x^(1/2)*Polylogarithms.polylog(5,-1*im*exp(im*(c+d*x^(1/2))))+-240*im*a*b*x*d^-4*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-20*im*a*b*d^-2*x^2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-8*im*a*b*d^-1*x^5//2*arctan(exp(im*(c+d*x^(1/2))))+20*im*a*b*d^-2*x^2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2))))+240*im*a*b*x*d^-4*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*sec(c+d*x^(1/2)))^2, x) == :((1/2)*a^2*x^2+3*b^2*d^-4*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-2*im*b^2*d^-1*x^3//2+2*b^2*d^-1*x^3//2*tan(c+d*x^(1/2))+6*x*b^2*d^-2*log(1+exp(2*im*(c+d*x^(1/2))))+-24*im*a*b*d^-4*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-24*a*b*d^-3*x^(1/2)*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+-6*im*b^2*d^-3*x^(1/2)*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+24*im*a*b*d^-4*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2))))+24*a*b*d^-3*x^(1/2)*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+-12*im*a*b*x*d^-2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-8*im*a*b*d^-1*x^3//2*arctan(exp(im*(c+d*x^(1/2))))+12*im*a*b*x*d^-2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^3*(a+b*sec(c+d*x^(1/2)))^-1, x) == :(1//4*a^-1*x^4+-10080*b*a^-1*d^-8*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(8,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10080*b*a^-1*d^-8*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(8,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-5040*b*x*a^-1*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-420*b*a^-1*d^-4*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-14*b*a^-1*d^-2*x^3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+14*b*a^-1*d^-2*x^3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+420*b*a^-1*d^-4*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+5040*b*x*a^-1*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10080*im*b*a^-1*d^-7*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(7,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-1680*im*b*a^-1*d^-5*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-84*im*b*a^-1*d^-3*x^5//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*b*a^-1*d^-1*x^7//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*b*a^-1*d^-1*x^7//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+84*im*b*a^-1*d^-3*x^5//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+1680*im*b*a^-1*d^-5*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10080*im*b*a^-1*d^-7*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(7,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*sec(c+d*x^(1/2)))^-1, x) == :(1//3*a^-1*x^3+-240*b*a^-1*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+240*b*a^-1*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-120*b*x*a^-1*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10*b*a^-1*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10*b*a^-1*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+120*b*x*a^-1*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-240*im*b*a^-1*d^-5*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-40*im*b*a^-1*d^-3*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*b*a^-1*d^-1*x^5//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*b*a^-1*d^-1*x^5//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+40*im*b*a^-1*d^-3*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+240*im*b*a^-1*d^-5*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*sec(c+d*x^(1/2)))^-1, x) == :((1/2)*a^-1*x^2+-12*b*a^-1*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+12*b*a^-1*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-6*b*x*a^-1*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+6*b*x*a^-1*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-12*im*b*a^-1*d^-3*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*b*a^-1*d^-1*x^3//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*b*a^-1*d^-1*x^3//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+12*im*b*a^-1*d^-3*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^3*(a+b*sec(c+d*x^(1/2)))^-2, x) == :(1//4*a^-2*x^4+-20160*b*a^-2*d^-8*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(8,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10080*a^-2*b^3*d^-8*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(8,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10080*a^-2*b^2*d^-8*(a^2+-1*b^2)^-1*Polylogarithms.polylog(7,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10080*a^-2*b^2*d^-8*(a^2+-1*b^2)^-1*Polylogarithms.polylog(7,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10080*a^-2*b^3*d^-8*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(8,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+20160*b*a^-2*d^-8*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(8,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10080*b*x*a^-2*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-5040*x*a^-2*b^2*d^-6*(a^2+-1*b^2)^-1*Polylogarithms.polylog(5,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-5040*x*a^-2*b^2*d^-6*(a^2+-1*b^2)^-1*Polylogarithms.polylog(5,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-5040*x*a^-2*b^3*d^-6*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(6,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-840*b*a^-2*d^-4*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-420*a^-2*b^3*d^-4*x^2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-28*b*a^-2*d^-2*x^3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-14*a^-2*b^3*d^-2*x^3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^2*d^-1*x^7//2*(a^2+-1*b^2)^-1+14*a^-2*b^2*d^-2*x^3*(a^2+-1*b^2)^-1*log(1+a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+14*a^-2*b^2*d^-2*x^3*(a^2+-1*b^2)^-1*log(1+a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+14*a^-2*b^3*d^-2*x^3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+28*b*a^-2*d^-2*x^3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+420*a^-2*b^2*d^-4*x^2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+420*a^-2*b^2*d^-4*x^2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+420*a^-2*b^3*d^-4*x^2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+840*b*a^-2*d^-4*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+5040*x*a^-2*b^3*d^-6*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(6,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10080*b*x*a^-2*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-20160*im*b*a^-2*d^-7*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(7,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10080*im*a^-2*b^2*d^-7*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(6,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10080*im*a^-2*b^2*d^-7*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(6,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10080*im*a^-2*b^3*d^-7*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(7,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-3360*im*b*a^-2*d^-5*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-1680*im*a^-2*b^3*d^-5*x^3//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-168*im*b*a^-2*d^-3*x^5//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-84*im*a^-2*b^2*d^-3*x^5//2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-84*im*a^-2*b^2*d^-3*x^5//2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-84*im*a^-2*b^3*d^-3*x^5//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*b*a^-2*d^-1*x^7//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^3*d^-1*x^7//2*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*a^-2*b^3*d^-1*x^7//2*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*a^-1*b^2*d^-1*x^7//2*(b+a*cos(c+d*x^(1/2)))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^(1/2))+4*im*b*a^-2*d^-1*x^7//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+84*im*a^-2*b^3*d^-3*x^5//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+168*im*b*a^-2*d^-3*x^5//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+1680*im*a^-2*b^2*d^-5*x^3//2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(4,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+1680*im*a^-2*b^2*d^-5*x^3//2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(4,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+1680*im*a^-2*b^3*d^-5*x^3//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+3360*im*b*a^-2*d^-5*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10080*im*a^-2*b^3*d^-7*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(7,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+20160*im*b*a^-2*d^-7*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(7,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*sec(c+d*x^(1/2)))^-2, x) == :(1//3*a^-2*x^3+-480*b*a^-2*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-240*a^-2*b^2*d^-6*(a^2+-1*b^2)^-1*Polylogarithms.polylog(5,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-240*a^-2*b^2*d^-6*(a^2+-1*b^2)^-1*Polylogarithms.polylog(5,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-240*a^-2*b^3*d^-6*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(6,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+240*a^-2*b^3*d^-6*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(6,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+480*b*a^-2*d^-6*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(6,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-240*b*x*a^-2*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-120*x*a^-2*b^3*d^-4*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-20*b*a^-2*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-10*a^-2*b^3*d^-2*x^2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^2*d^-1*x^5//2*(a^2+-1*b^2)^-1+10*a^-2*b^2*d^-2*x^2*(a^2+-1*b^2)^-1*log(1+a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10*a^-2*b^2*d^-2*x^2*(a^2+-1*b^2)^-1*log(1+a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+10*a^-2*b^3*d^-2*x^2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+20*b*a^-2*d^-2*x^2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+120*x*a^-2*b^2*d^-4*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+120*x*a^-2*b^2*d^-4*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+120*x*a^-2*b^3*d^-4*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+240*b*x*a^-2*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-480*im*b*a^-2*d^-5*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-240*im*a^-2*b^3*d^-5*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-80*im*b*a^-2*d^-3*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-40*im*a^-2*b^2*d^-3*x^3//2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-40*im*a^-2*b^2*d^-3*x^3//2*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-40*im*a^-2*b^3*d^-3*x^3//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*b*a^-2*d^-1*x^5//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^3*d^-1*x^5//2*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*a^-2*b^3*d^-1*x^5//2*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*a^-1*b^2*d^-1*x^5//2*(b+a*cos(c+d*x^(1/2)))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^(1/2))+4*im*b*a^-2*d^-1*x^5//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+40*im*a^-2*b^3*d^-3*x^3//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+80*im*b*a^-2*d^-3*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+240*im*a^-2*b^2*d^-5*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(4,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+240*im*a^-2*b^2*d^-5*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(4,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+240*im*a^-2*b^3*d^-5*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+480*im*b*a^-2*d^-5*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*sec(c+d*x^(1/2)))^-2, x) == :((1/2)*a^-2*x^2+-24*b*a^-2*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-12*a^-2*b^3*d^-4*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+12*a^-2*b^2*d^-4*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+12*a^-2*b^2*d^-4*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+12*a^-2*b^3*d^-4*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+24*b*a^-2*d^-4*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-12*b*x*a^-2*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-6*x*a^-2*b^3*d^-2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^2*d^-1*x^3//2*(a^2+-1*b^2)^-1+6*x*a^-2*b^2*d^-2*(a^2+-1*b^2)^-1*log(1+a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+6*x*a^-2*b^2*d^-2*(a^2+-1*b^2)^-1*log(1+a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+6*x*a^-2*b^3*d^-2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+12*b*x*a^-2*d^-2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-24*im*b*a^-2*d^-3*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-12*im*a^-2*b^2*d^-3*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-12*im*a^-2*b^2*d^-3*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-12*im*a^-2*b^3*d^-3*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*b*a^-2*d^-1*x^3//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^3*d^-1*x^3//2*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*a^-2*b^3*d^-1*x^3//2*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*a^-1*b^2*d^-1*x^3//2*(b+a*cos(c+d*x^(1/2)))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^(1/2))+4*im*b*a^-2*d^-1*x^3//2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+12*im*a^-2*b^3*d^-3*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+24*im*b*a^-2*d^-3*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^3//2*(a+b*sec(c+d*x^(1/2))), x) == :(2//5*a*x^5//2+-48*b*d^-5*Polylogarithms.polylog(5,im*exp(im*(c+d*x^(1/2))))+48*b*d^-5*Polylogarithms.polylog(5,-1*im*exp(im*(c+d*x^(1/2))))+-24*b*x*d^-3*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+24*b*x*d^-3*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+-48*im*b*d^-4*x^(1/2)*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-8*im*b*d^-2*x^3//2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-4*im*b*d^-1*x^2*arctan(exp(im*(c+d*x^(1/2))))+8*im*b*d^-2*x^3//2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2))))+48*im*b*d^-4*x^(1/2)*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^(1/2)*(a+b*sec(c+d*x^(1/2))), x) == :(2//3*a*x^3//2+-4*b*d^-3*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+4*b*d^-3*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+-4*im*b*x*d^-1*arctan(exp(im*(c+d*x^(1/2))))+-4*im*b*d^-2*x^(1/2)*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+4*im*b*d^-2*x^(1/2)*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^-1//2*(a+b*sec(c+d*x^(1/2))), x) == :(2*a*x^(1/2)+2*b*d^-1*arctanh(sin(c+d*x^(1/2))))
@test integrate(x^3//2*(a+b*sec(c+d*x^(1/2)))^2, x) == :(2//5*a^2*x^5//2+-96*a*b*d^-5*Polylogarithms.polylog(5,im*exp(im*(c+d*x^(1/2))))+-2*im*b^2*d^-1*x^2+2*b^2*d^-1*x^2*tan(c+d*x^(1/2))+6*im*b^2*d^-5*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+8*b^2*d^-2*x^3//2*log(1+exp(2*im*(c+d*x^(1/2))))+12*b^2*d^-4*x^(1/2)*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+96*a*b*d^-5*Polylogarithms.polylog(5,-1*im*exp(im*(c+d*x^(1/2))))+-48*a*b*x*d^-3*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+-12*im*x*b^2*d^-3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+48*a*b*x*d^-3*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+-96*im*a*b*d^-4*x^(1/2)*Polylogarithms.polylog(4,-1*im*exp(im*(c+d*x^(1/2))))+-16*im*a*b*d^-2*x^3//2*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+-8*im*a*b*d^-1*x^2*arctan(exp(im*(c+d*x^(1/2))))+16*im*a*b*d^-2*x^3//2*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2))))+96*im*a*b*d^-4*x^(1/2)*Polylogarithms.polylog(4,im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^(1/2)*(a+b*sec(c+d*x^(1/2)))^2, x) == :(2//3*a^2*x^3//2+-8*a*b*d^-3*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^(1/2))))+-2*im*x*b^2*d^-1+-2*im*b^2*d^-3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+2*x*b^2*d^-1*tan(c+d*x^(1/2))+4*b^2*d^-2*x^(1/2)*log(1+exp(2*im*(c+d*x^(1/2))))+8*a*b*d^-3*Polylogarithms.polylog(3,im*exp(im*(c+d*x^(1/2))))+-8*im*a*b*x*d^-1*arctan(exp(im*(c+d*x^(1/2))))+-8*im*a*b*d^-2*x^(1/2)*Polylogarithms.polylog(2,im*exp(im*(c+d*x^(1/2))))+8*im*a*b*d^-2*x^(1/2)*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^(1/2)))))
@test integrate(x^-1//2*(a+b*sec(c+d*x^(1/2)))^2, x) == :(2*a^2*x^(1/2)+2*b^2*d^-1*tan(c+d*x^(1/2))+4*a*b*d^-1*arctanh(sin(c+d*x^(1/2))))
@test integrate(x^3//2*(a+b*sec(c+d*x^(1/2)))^-1, x) == :(2//5*a^-1*x^5//2+-48*im*b*a^-1*d^-5*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-48*b*a^-1*d^-4*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-8*b*a^-1*d^-2*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+8*b*a^-1*d^-2*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*im*b*a^-1*d^-5*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*b*a^-1*d^-4*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-24*im*b*x*a^-1*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*b*a^-1*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*b*a^-1*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+24*im*b*x*a^-1*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^(1/2)*(a+b*sec(c+d*x^(1/2)))^-1, x) == :(2//3*a^-1*x^3//2+-4*im*b*a^-1*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*b*a^-1*d^-2*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+4*im*b*a^-1*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+4*b*a^-1*d^-2*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*b*x*a^-1*d^-1*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*b*x*a^-1*d^-1*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^-1//2*(a+b*sec(c+d*x^(1/2)))^-1, x) == :(2*a^-1*x^(1/2)+-4*b*a^-1*d^-1*(a+b)^-1//2*(a+-1b)^-1//2*arctanh((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*c+(1/2)*d*x^(1/2))))
@test integrate(x^3//2*(a+b*sec(c+d*x^(1/2)))^-2, x) == :(2//5*a^-2*x^5//2+-96*im*b*a^-2*d^-5*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-96*b*a^-2*d^-4*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-48*im*a^-2*b^3*d^-5*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-48*a^-2*b^3*d^-4*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-16*b*a^-2*d^-2*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-8*a^-2*b^3*d^-2*x^3//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^2*d^-1*x^2*(a^2+-1*b^2)^-1+8*a^-2*b^2*d^-2*x^3//2*(a^2+-1*b^2)^-1*log(1+a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+8*a^-2*b^2*d^-2*x^3//2*(a^2+-1*b^2)^-1*log(1+a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+8*a^-2*b^3*d^-2*x^3//2*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+16*b*a^-2*d^-2*x^3//2*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*im*a^-2*b^2*d^-5*(a^2+-1*b^2)^-1*Polylogarithms.polylog(4,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*im*a^-2*b^2*d^-5*(a^2+-1*b^2)^-1*Polylogarithms.polylog(4,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*im*a^-2*b^3*d^-5*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(5,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*a^-2*b^2*d^-4*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*a^-2*b^2*d^-4*x^(1/2)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(3,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*a^-2*b^3*d^-4*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(4,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+96*im*b*a^-2*d^-5*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(5,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+96*b*a^-2*d^-4*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(4,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-48*im*b*x*a^-2*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-24*im*x*a^-2*b^2*d^-3*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-24*im*x*a^-2*b^2*d^-3*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-24*im*x*a^-2*b^3*d^-3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*b*a^-2*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*a^-2*b^3*d^-1*x^2*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*a^-2*b^3*d^-1*x^2*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*a^-1*b^2*d^-1*x^2*(b+a*cos(c+d*x^(1/2)))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^(1/2))+4*im*b*a^-2*d^-1*x^2*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+24*im*x*a^-2*b^3*d^-3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+48*im*b*x*a^-2*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^(1/2)*(a+b*sec(c+d*x^(1/2)))^-2, x) == :(2//3*a^-2*x^3//2+-8*im*b*a^-2*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-8*b*a^-2*d^-2*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*a^-2*b^2*d^-3*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*a^-2*b^2*d^-3*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*a^-2*b^3*d^-3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*a^-2*b^3*d^-2*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*x*a^-2*b^2*d^-1*(a^2+-1*b^2)^-1+4*im*a^-2*b^3*d^-3*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+4*a^-2*b^2*d^-2*x^(1/2)*(a^2+-1*b^2)^-1*log(1+a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+4*a^-2*b^2*d^-2*x^(1/2)*(a^2+-1*b^2)^-1*log(1+a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+4*a^-2*b^3*d^-2*x^(1/2)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+8*im*b*a^-2*d^-3*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+8*b*a^-2*d^-2*x^(1/2)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-4*im*b*x*a^-2*d^-1*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+-2*im*x*a^-2*b^3*d^-1*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*im*x*a^-2*b^3*d^-1*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2))))+2*x*a^-1*b^2*d^-1*(b+a*cos(c+d*x^(1/2)))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^(1/2))+4*im*b*x*a^-2*d^-1*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^(1/2)))))
@test integrate(x^-1//2*(a+b*sec(c+d*x^(1/2)))^-2, x) == :(2*a^-2*x^(1/2)+2*a^-1*b^2*d^-1*(a+b*sec(c+d*x^(1/2)))^-1*(a^2+-1*b^2)^-1*tan(c+d*x^(1/2))+-4*b*a^-2*d^-1*(a+b)^-3//2*(a+-1b)^-3//2*(-1*b^2+2*a^2)*arctanh((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*c+(1/2)*d*x^(1/2))))
@test integrate((e*x)^(-1+n)*(a+b*sec(c+d*x^n)), x) == :(a*e^-1*n^-1*(e*x)^n+b*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^n*arctanh(sin(c+d*x^n)))
@test integrate((e*x)^(-1+2n)*(a+b*sec(c+d*x^n)), x) == :((1/2)*a*e^-1*n^-1*(e*x)^(2n)+im*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^n)))+-1*im*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*Polylogarithms.polylog(2,im*exp(im*(c+d*x^n)))+-2*im*b*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*arctan(exp(im*(c+d*x^n))))
@test integrate((e*x)^(-1+3n)*(a+b*sec(c+d*x^n)), x) == :(1//3*a*e^-1*n^-1*(e*x)^(3n)+-2*b*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^n)))+2*b*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*Polylogarithms.polylog(3,im*exp(im*(c+d*x^n)))+-2*im*b*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*arctan(exp(im*(c+d*x^n)))+-2*im*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*Polylogarithms.polylog(2,im*exp(im*(c+d*x^n)))+2*im*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^n))))
@test integrate((e*x)^(-1+n)*(a+b*sec(c+d*x^n))^2, x) == :(a^2*e^-1*n^-1*(e*x)^n+b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^n*tan(c+d*x^n)+2*a*b*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^n*arctanh(sin(c+d*x^n)))
@test integrate((e*x)^(-1+2n)*(a+b*sec(c+d*x^n))^2, x) == :((1/2)*a^2*e^-1*n^-1*(e*x)^(2n)+b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*tan(c+d*x^n)+b^2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*log(cos(c+d*x^n))+-4*im*a*b*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*arctan(exp(im*(c+d*x^n)))+-2*im*a*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*Polylogarithms.polylog(2,im*exp(im*(c+d*x^n)))+2*im*a*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^n))))
@test integrate((e*x)^(-1+3n)*(a+b*sec(c+d*x^n))^2, x) == :(1//3*a^2*e^-1*n^-1*(e*x)^(3n)+b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*tan(c+d*x^n)+-1*im*b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)+2*b^2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*log(1+exp(2*im*(c+d*x^n)))+-1*im*b^2*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^n)))+-4*a*b*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*Polylogarithms.polylog(3,-1*im*exp(im*(c+d*x^n)))+4*a*b*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*Polylogarithms.polylog(3,im*exp(im*(c+d*x^n)))+-4*im*a*b*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*arctan(exp(im*(c+d*x^n)))+-4*im*a*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*Polylogarithms.polylog(2,im*exp(im*(c+d*x^n)))+4*im*a*b*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*Polylogarithms.polylog(2,-1*im*exp(im*(c+d*x^n))))
@test integrate((e*x)^(-1+n)*(a+b*sec(c+d*x^n))^-1, x) == :(a^-1*e^-1*n^-1*(e*x)^n+-2*b*a^-1*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^n*(a+b)^-1//2*(a+-1b)^-1//2*arctanh((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*c+(1/2)*d*x^n)))
@test integrate((e*x)^(-1+2n)*(a+b*sec(c+d*x^n))^-1, x) == :((1/2)*a^-1*e^-1*n^-1*(e*x)^(2n)+b*a^-1*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-1*b*a^-1*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+im*b*a^-1*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-1*im*b*a^-1*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n))))
@test integrate((e*x)^(-1+3n)*(a+b*sec(c+d*x^n))^-1, x) == :(1//3*a^-1*e^-1*n^-1*(e*x)^(3n)+-2*b*a^-1*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*b*a^-1*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+im*b*a^-1*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-1*im*b*a^-1*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*im*b*a^-1*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*im*b*a^-1*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n))))
@test integrate((e*x)^(-1+n)*(a+b*sec(c+d*x^n))^-2, x) == :(a^-2*e^-1*n^-1*(e*x)^n+a^-1*b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^n*(a+b*sec(c+d*x^n))^-1*(a^2+-1*b^2)^-1*tan(c+d*x^n)+-2*b*a^-2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^n*(a+b)^-3//2*(a+-1b)^-3//2*(-1*b^2+2*a^2)*arctanh((a+b)^-1//2*(a+-1b)^(1/2)*tan((1/2)*c+(1/2)*d*x^n)))
@test integrate((e*x)^(-1+2n)*(a+b*sec(c+d*x^n))^-2, x) == :((1/2)*a^-2*e^-1*n^-1*(e*x)^(2n)+a^-2*b^2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*(a^2+-1*b^2)^-1*log(b+a*cos(c+d*x^n))+a^-2*b^3*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-1*a^-2*b^3*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*b*a^-2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*b*a^-2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+im*a^-2*b^3*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+a^-1*b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*(b+a*cos(c+d*x^n))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^n)+-1*im*a^-2*b^3*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*im*b*a^-2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*im*b*a^-2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(2n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n))))
@test integrate((e*x)^(-1+3n)*(a+b*sec(c+d*x^n))^-2, x) == :(1//3*a^-2*e^-1*n^-1*(e*x)^(3n)+-1*im*a^-2*b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(a^2+-1*b^2)^-1+-4*b*a^-2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*a^-2*b^3*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*a^-2*b^2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(a^2+-1*b^2)^-1*log(1+a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*a^-2*b^2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(a^2+-1*b^2)^-1*log(1+a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*a^-2*b^3*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(2,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+4*b*a^-2*d^-2*e^-1*n^-1*x^(-2n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(2,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+im*a^-2*b^3*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(b^2+-1*a^2)^-3//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+a^-1*b^2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(b+a*cos(c+d*x^n))^-1*(a^2+-1*b^2)^-1*sin(c+d*x^n)+-1*im*a^-2*b^3*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(b^2+-1*a^2)^-3//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-4*im*b*a^-2*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*im*b*a^-2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*im*a^-2*b^2*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*im*a^-2*b^2*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(a^2+-1*b^2)^-1*Polylogarithms.polylog(2,-1*a*(b+-1*im*(a^2+-1*b^2)^(1/2))^-1*exp(im*(c+d*x^n)))+-2*im*a^-2*b^3*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*im*b*a^-2*d^-1*e^-1*n^-1*x^(-1n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*log(1+a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+2*im*a^-2*b^3*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(b^2+-1*a^2)^-3//2*Polylogarithms.polylog(3,-1*a*(b+(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n)))+4*im*b*a^-2*d^-3*e^-1*n^-1*x^(-3n)*(e*x)^(3n)*(b^2+-1*a^2)^-1//2*Polylogarithms.polylog(3,-1*a*(b+-1*(b^2+-1*a^2)^(1/2))^-1*exp(im*(c+d*x^n))))

max_stars_repo_path: test/4 Trig functions/4.5 Secant/4.5.11 (e x)^m (a+b sec(c+d x^n))^p.jl