text:
# Tiger Problem
# POMDP, Robust POMDP, and Robust ρPOMDP formulations

struct TigerPOMDP <: POMDP{Symbol, Symbol, Symbol}
    discount::Float64
end

struct TigerRPOMDP <: RPOMDP{Symbol, Symbol, Symbol}
    discount::Float64
    usize::Float64
end

struct TigerIPOMDP <: IPOMDP{Symbol, Symbol, Symbol}
    discount::Float64
    inforeward::Vector{Vector{Float64}}
end

struct TigerRIPOMDP <: RIPOMDP{Symbol, Symbol, Symbol}
    discount::Float64
    usize::Float64
    inforeward::Vector{Vector{Float64}}
end

const discount_tiger = 0.95
const states_tiger = [:tigerleft, :tigerright]
const actions_tiger = [:listen, :openleft, :openright]
const observations_tiger = [:tigerleft, :tigerright]
const inforeward_tiger = [[1.0, -1.0], [-1.0, 1.0], [0.2, 0.2]]

TigerPOMDP() = TigerPOMDP(discount_tiger)

TigerRPOMDP(err) = TigerRPOMDP(discount_tiger, err)
TigerRPOMDP() = TigerRPOMDP(discount_tiger, 0.025)

TigerIPOMDP(alphas) = TigerIPOMDP(discount_tiger, alphas)
TigerIPOMDP() = TigerIPOMDP(discount_tiger, inforeward_tiger)

TigerRIPOMDP(err::Float64, alphas::Vector{Vector{Float64}}) = TigerRIPOMDP(discount_tiger, err, alphas)
TigerRIPOMDP() = TigerRIPOMDP(discount_tiger, 0.025, inforeward_tiger)

states(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = states_tiger
actions(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = actions_tiger
observations(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = observations_tiger

n_states(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = 2
n_actions(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = 3
n_observations(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = 2

state_index(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}, s::Symbol) = s == :tigerleft ? 1 : 2
action_index(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}, a::Symbol) = a == :listen ? 1 : a == :openleft ? 2 : 3
observation_index(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}, z::Symbol) = z == :tigerleft ? 1 : 2
obs_index(prob::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}, z::Symbol) = observation_index(prob, z)

initial_state_distribution(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = SparseCat([:tigerleft, :tigerright], [0.5, 0.5])
initial_belief(::Union{TigerPOMDP, TigerRPOMDP, TigerIPOMDP, TigerRIPOMDP}) = [0.5, 0.5]

# Nominal transition function distributions
const pϵ_tiger = 1e-6
const tdist_tiger = [SparseCat(states_tiger, [1.0, 0.0]) SparseCat(states_tiger, [0.5, 0.5]) SparseCat(states_tiger, [0.5, 0.5]);
                    SparseCat(states_tiger, [0.0, 1.0]) SparseCat(states_tiger, [0.5, 0.5]) SparseCat(states_tiger, [0.5, 0.5])]

# Nominal transition function array tarray_tiger[s,a,s'] = Pr(s' | s ,a)
const tarray_tiger = cat(3, [1.0 0.5 0.5; 0.0 0.5 0.5], [0.0 0.5 0.5; 1.0 0.5 0.5]) # P(tigerleft | s, a), P(tigerright | s, a)

# Nominal transitions
transition(prob::Union{TigerIPOMDP, TigerPOMDP}, s::Symbol, a::Symbol) = tdist_tiger[state_index(prob, s), action_index(prob, a)]
transition(prob::Union{TigerIPOMDP, TigerPOMDP}, s::Symbol, a::Symbol, sp::Symbol) = tarray_tiger[state_index(prob, s), action_index(prob, a), state_index(prob, sp)]

# Robust transitions
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PtLower = Plower(t|s,a) = t_lower[s,a,sp]
function transition(prob::Union{TigerRPOMDP, TigerRIPOMDP}, s::Symbol, a::Symbol)
    si, ai = state_index(prob, s), action_index(prob, a)
    plower = max.(tarray_tiger[si,ai,:] - prob.usize, 0.0 + pϵ_tiger)
    pupper = min.(tarray_tiger[si,ai,:] + prob.usize, 1.0 - pϵ_tiger)
    plower, pupper
end
function transition(prob::Union{TigerRPOMDP, TigerRIPOMDP}, s::Symbol, a::Symbol, sp::Symbol)
    si, ai, spi = state_index(prob, s), action_index(prob, a), state_index(prob, sp)
    plower = max.(tarray_tiger[si, ai, spi] - prob.usize, 0.0 + pϵ_tiger)
    pupper = min.(tarray_tiger[si, ai, spi] + prob.usize, 1.0 - pϵ_tiger)
    plower, pupper
end

# Nominal observation function distributions
const odist_tiger = [SparseCat(observations_tiger, [0.85 0.15]) SparseCat(observations_tiger, [0.15 0.85]);
                    SparseCat(observations_tiger, [0.5 0.5]) SparseCat(observations_tiger, [0.5 0.5]);
                    SparseCat(observations_tiger, [0.5 0.5]) SparseCat(observations_tiger, [0.5 0.5])]

# Nominal observation function array oarray_tiger[a,sp,z] = Pr(z|a,sp)
const oarray_tiger = cat(3, [0.85 0.15; 0.5 0.5; 0.5 0.5], [0.15 0.85; 0.5 0.5; 0.5 0.5]) # Pr(tigerleft|a,sp), Pr(tigerright|a,sp)

# Nominal observations
observation(prob::Union{TigerIPOMDP, TigerPOMDP}, a::Symbol, sp::Symbol) = odist_tiger[action_index(prob, a), state_index(prob, sp)]
observation(prob::Union{TigerIPOMDP, TigerPOMDP}, a::Symbol, sp::Symbol, z::Symbol) = oarray_tiger[action_index(prob, a), state_index(prob, sp), observation_index(prob, z)]

# Robust observations
# P(z,t|s,a) = P(z|s,a,t)P(t|s,a)
# PztLower = PzLower * PtLower
# PzLower = Plower(z|s,a,t) = o_lower[a,sp,z]
function observation(prob::Union{TigerRPOMDP, TigerRIPOMDP}, a::Symbol, sp::Symbol)
    ai, spi = action_index(prob, a), state_index(prob, sp)
    plower = max.(oarray_tiger[ai,spi,:] - prob.usize, 0.0 + pϵ_tiger)
    pupper = min.(oarray_tiger[ai,spi,:] + prob.usize, 1.0 - pϵ_tiger)
    plower, pupper
end
function observation(prob::Union{TigerRPOMDP, TigerRIPOMDP}, a::Symbol, sp::Symbol, z::Symbol)
    ai, spi, zi = action_index(prob, a), state_index(prob, sp), observation_index(prob, z)
    plower = max.(oarray_tiger[ai, spi, zi] - prob.usize, 0.0 + pϵ_tiger)
    pupper = min.(oarray_tiger[ai, spi, zi] + prob.usize, 1.0 - pϵ_tiger)
    plower, pupper
end

function dynamics(prob::Union{TigerPOMDP,TigerIPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        pt = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pz = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        p[tind, zind, sind, aind] = pz * pt
    end
    p
end

function dynamics(prob::Union{TigerRPOMDP,TigerRIPOMDP})
    ns = n_states(prob)
    nz = n_observations(prob)
    na = n_actions(prob)
    tarr = states(prob)
    sarr = states(prob)
    aarr = actions(prob)
    zarr = observations(prob)
    p = Array{Float64}(ns,nz,ns,na)
    plower = Array{Float64}(ns,nz,ns,na)
    pupper = Array{Float64}(ns,nz,ns,na)
    for tind in 1:ns, zind in 1:nz, sind in 1:ns, aind in 1:na
        ptl, ptu = transition(prob, sarr[sind], aarr[aind], tarr[tind])
        pzl, pzu = observation(prob, aarr[aind], tarr[tind], zarr[zind])
        plower[tind, zind, sind, aind] = pzl * ptl
        pupper[tind, zind, sind, aind] = pzu * ptu
    end
    plower, pupper
end

function reward(prob::Union{TigerPOMDP, TigerRPOMDP}, s::Symbol, a::Symbol)
    r = 0.0
    if a == :listen
        r = -1.0
    elseif a == :openleft
        r = (s == :tigerleft) ? -100.0 : 10.0
    elseif a == :openright
        r = (s == :tigerright) ? -100.0 : 10.0
    end
    return r
end
reward(prob::Union{TigerPOMDP,TigerRPOMDP}, b::Vector{Float64}, s::Symbol, a::Symbol, sp::Symbol) = reward(prob,s,a)

function reward(prob::Union{TigerIPOMDP, TigerRIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    for α in prob.inforeward
        rmax = max(rmax, dot(α, b))
    end
    rmax
end
reward(prob::Union{TigerIPOMDP,TigerRIPOMDP}, b::Vector{Float64}, s::Symbol, a::Symbol, sp::Symbol) = reward(prob,b,a)

function rewardalpha(prob::Union{TigerIPOMDP, TigerRIPOMDP}, b::Vector{Float64}, a::Symbol)
    rmax = -Inf
    ralpha = nothing
    for α in prob.inforeward
        rnext = dot(α, b)
        if rnext > rmax
            rmax = rnext
            ralpha = copy(α)
        end
    end
    ralpha
end

discount(p::Union{TigerPOMDP,TigerRPOMDP,TigerIPOMDP,TigerRIPOMDP}) = p.discount

# Simulation Functions
function initial_state(prob::Union{TigerPOMDP,TigerIPOMDP,TigerRPOMDP,TigerRIPOMDP}, rng::AbstractRNG)
    d = initial_state_distribution(prob)
    return rand(rng, d)
end

function generate_sor(prob::Union{TigerPOMDP,TigerIPOMDP}, b, s, a, rng::AbstractRNG)
    sp = rand(rng, transition(prob, s, a))
    o = rand(rng, observation(prob, a, sp))
    r = reward(prob, b, s, a, sp)
    sp, o, r
end

function generate_sor(prob::Union{TigerRPOMDP,TigerRIPOMDP}, b, s, a, rng::AbstractRNG)
    tdist = SparseCat(states(prob), psample(transition(prob, s, a)...))
    sp = rand(rng, tdist)
    odist = SparseCat(observations(prob), psample(observation(prob, a, sp)...))
    o = rand(rng, odist)
    r = reward(prob, b, s, a, sp)
    sp, o, r
end


text:
module Filtering

using LinearAlgebra,
    Statistics,
    Optim,
    ProgressMeter,
    StatsBase,
    DiffResults,
    Distributions,
    Logging,
    Random,
    ForwardDiff,
    StatsFuns

include("models.jl")
include("filter.jl")
include("smoother.jl")
include("em.jl")
include("resampling.jl")
include("particle.jl")
include("score_filter.jl")
#include("bare_particle.jl")
#include("gradient_particle.jl") # This may be outdated
#include("online_particle.jl")
#include("threads.jl")
include("optimizers.jl")

# TODO: These models should be moved to a FilteringModels package
include("models/ar1.jl")
include("models/smoothingspline.jl")
include("models/ss_tfi.jl")
include("models/tfi.jl")

end # module

ext: jl


text:
"""
    reflection_coefficient(ω::T, m_wave::MatchPlaneWaveMode,
        source::PlaneSource, material::Material{2,Halfspace{T,2}};

Calculate the reflection coefficient from a matched wave. This requires using both the discrete part and effective wavemode of the matched wave.
"""
function reflection_coefficient(ω::T, m_wave::MatchPlaneWaveMode{T},
    source::PlaneSource{T,2,1,Acoustic{T,2}}, material::Material{2,Halfspace{T,2}};
     # medium::PhysicalMedium{T}, specie::Specie{T};
    kws...) where {T<:AbstractFloat}

    R_eff = reflection_coefficient(ω, m_wave.PlaneWaveModes, source, material; x = m_wave.x_match[end], kws...)
    R_discrete = reflection_coefficient(ω, m_wave.discrete_wave, source, material; kws...)

    return  R_eff + R_discrete
end


text:
using OrdinaryDiffEq, DiffEqDevTools
using Test, Random
Random.seed!(100)

dts = 1 .//2 .^(9:-1:5)
testTol = 0.2

f_dae_linear = (res, du, u, p, t) -> (@. res = du - u)
f_dae_linear_analytic = (du0, u0, p, t) -> @. u0*exp(t)
prob_dae_linear_iip = DAEProblem(
				  DAEFunction(f_dae_linear;analytic=f_dae_linear_analytic),
				  [1.0,1.0], [1.0,1.0], (0.0, 1.0)
				  )

@testset "DAE Solver Convergence Tests (in-place)" begin
	prob = prob_dae_linear_iip

	sim11 = test_convergence(dts,prob,DImplicitEuler())
	@test sim11.𝒪est[:final] ≈ 1 atol=testTol

	sim12 = test_convergence(dts,prob,DImplicitEuler(;autodiff=false))
	@test sim12.𝒪est[:final] ≈ 1 atol=testTol

	sim13 = test_convergence(dts,prob,DABDF2())
	@test sim13.𝒪est[:final] ≈ 2 atol=testTol

	sim14 = test_convergence(dts,prob,DABDF2(;autodiff=false))
	@test sim14.𝒪est[:final] ≈ 2 atol=testTol

	@test_nowarn solve(prob, DFBDF())
end

f_dae_linear = (du, u, p, t) -> (@. du - u)
f_dae_linear_analytic = (du0, u0, p, t) -> @. u0*exp(t)
prob_dae_linear_oop = DAEProblem(
				  DAEFunction(f_dae_linear;analytic=f_dae_linear_analytic),
				  1.0, 1.0, (0.0, 1.0)
				  )

@testset "DAE Solver Convergence Tests (out-of-place)" begin
	prob = prob_dae_linear_oop

	sim21 = test_convergence(dts,prob,DImplicitEuler())
	@test sim21.𝒪est[:final] ≈ 1 atol=testTol

	sim22 = test_convergence(dts,prob,DImplicitEuler(;autodiff=false))
	@test sim22.𝒪est[:final] ≈ 1 atol=testTol

	sim23 = test_convergence(dts,prob,DABDF2())
	@test sim23.𝒪est[:final] ≈ 2 atol=testTol

	sim24 = test_convergence(dts,prob,DABDF2(;autodiff=false))
	@test sim24.𝒪est[:final] ≈ 2 atol=testTol

	@test_nowarn solve(prob, DFBDF())
end


text:
function crawl(::Val{:if}, test, consequent, alternate=nothing)
    return _crawl_if(:if, test, consequent, alternate)
end

function crawl(::Val{:elseif}, test, consequent, alternate=nothing)
    # There's some weirdness in the way Julia parses `if ... elseif ...`.
    # It wraps the test expression in a block (probably to include the
    # LineNumberNode information).
    @assert (test isa Expr && test.head == :block) "Unparseable :elseif expression."
    test_args = filter(x -> !isa(x, LineNumberNode), test.args)
    @assert (length(test_args) == 1) "Invalid test in :elseif expression ($(test_args))."
    test = test_args[1]
    return _crawl_if(:elseif, test, consequent, alternate)
end

function _crawl_if(head::Symbol, test, consequent, alternate)
    body = crawl.(
        alternate === nothing
        ? [test, consequent]
        : [test, consequent, alternate]
    )
    return :(JSAST($(QuoteNode(head)), $(body...)))
end

function deparse(
        ::Val{:if},
        test::JSNode,
        consequent::JSNode,
        alternate::Union{JSNode, Nothing} = nothing
)::JSString
    return jsstring(
        "if (", deparse(test), ") ",
        deparse(consequent),
        _deparse_if_alternate(alternate),
    )
end

# The alternate of an if statement be nothing, an "else" clause (in which case,
# the head of the expression is simply a block), or an "elseif" (in which case
# the head of the expression is :elseif).
_deparse_if_alternate(::Nothing) = jsstring("")
function _deparse_if_alternate(expr::JSNode)
    if expr.head == :elseif
        return jsstring(
            # This is a hack that uses the fact that we can treat the args of an
            # elseif as an if and just change the initial "if" to "elseif".
            " else", deparse(Val(:if), expr.args...)
        )
    elseif expr.head == :block
        return jsstring(" else ", deparse(expr))
    end
    error(
        "Invalid alternate expression in :if expression " *
        "(expected :elseif or :block expression, got $(expr))."
    )
end

# While
function crawl(::Val{:while}, test, body)
    return :(JSAST(:while, $(crawl(test)), $(crawl(body))))
end
function deparse(
        ::Val{:while},
        test::JSNode,
        body::JSNode,
)::JSString
    return jsstring(
        "while (", deparse(test), ") ",
        deparse(body),
    )
end

# For
function crawl(::Val{:for}, iterspec, body)
    return :(JSAST(
        :for,
        $(crawl_iterspec(iterspec)),
        $(crawl(body)),
    ))
end
function deparse(::Val{:for}, iterspec, body)
    return jsstring(
        "for (", deparse(iterspec), ") ",
        deparse(body),
    )
end

"""
Transform a Julia iteration specification into JavaScript.

Importantly, this transforms `for x in iterable` in Julia into
`for x of iterable` in JavaScript (since `of` is more often what is desired).
"""
function crawl_iterspec(iterspec)
    # In the future, we might be able to transform compound iteration
    # specifications (such as `for i = 1:10, j = i:10`) into nested for loops,
    # but for now, we don't support that syntax.
    if iterspec.head != :(=)
        error(
            "Compound iteration specifications in for loops " *
            "are not supported by JSExpr."
        )
    end
    lhs, rhs = iterspec.args
    return :(JSAST(
        :of,
        JSAST(
            :let,
            $(crawl(lhs)),
        ),
        $(crawl(rhs)),
    ))
end

for keyword in (:break, :continue)
    sym = QuoteNode(keyword)
    @eval function crawl(::Val{$sym})
        return Expr(
            :call,
            JSTerminal,
            JSString($(string(keyword))),
        )
    end
end


text:
export convertid

"convid(::Type{T}, id) convert `id` into type `T`"
function convertid end

convertid(::Type{T}, id::T) where {T} = id
convertid(::Type{VectorID}, (i1,)::TupleID{1}) = Int[i1]
convertid(::Type{VectorID}, (i1, i2)::TupleID{2}) = Int[i1, i2]
convertid(::Type{VectorID}, (i1, i2, i3)::TupleID{3}) = Int[i1, i2, i3]
convertid(::Type{VectorID}, (i1, i2, i3, i4)::TupleID{4}) = Int[i1, i2, i3, i4]
convertid(::Type{VectorID}, (i1, i2, i3, i4, i5)::TupleID{5}) = Int[i1, i2, i3, i4, i5]
convertid(::Type{VectorID}, id::TupleID) = Int[id...]
convertid(::Type{VectorID}, ci::CartesianIndex{1}) = Int[ci.I[1]]
convertid(::Type{VectorID}, ci::CartesianIndex{2}) = Int[ci.I[1], ci.I[2]]
convertid(::Type{VectorID}, ci::CartesianIndex{3}) = Int[ci.I[1], ci.I[2], ci.I[3]]

append(x::X, y) where X = append(x, convertid(X, y))

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
@testset "Affine series simulation" begin
    y = collect(1.:30)
    model = structural(y, 2)

    @test isa(model, StateSpaceModel)
    @test model.mode == "time-invariant"

    ss = statespace(model)
    @test ss.filter_type <: KalmanFilter

    sim = simulate(ss, 20, 100)
    media_sim = mean(sim, dims = 3)[:, 1]

    @test size(sim) == (20, 1, 100)
    @test media_sim ≈ collect(31.:50) rtol = 1e-3
    @test var(sim[1, end, :]) < 1e-4
    compare_forecast_simulation(ss, 20, 1000, 1e-3)
end


text:
using Images, IndirectArrays, Test, Suppressor
using Images.MappedArrays: mappedarray

@testset "ColorizedArray" begin
    intensity = [0.1 0.3; 0.2 0.4]
    labels = IndirectArray([1 2; 2 1], [RGB(1,0,0), RGB(0,0,1)])
    colorized_A = @suppress_err ColorizedArray(intensity, labels)
    mapped_A = mappedarray(*, intensity, labels)

    target = intensity .* labels

    for A in (colorized_A, mapped_A)
        @test eltype(A) == RGB{Float64}
        @test size(A) == (2,2)
        @test axes(A) == (Base.OneTo(2), Base.OneTo(2))
        for i = 1:4
            @test A[i] === target[i]
        end
        for j = 1:2, i = 1:2
            @test A[i,j] === target[i,j]
        end
        for (a,t) in zip(A, target)
            @test a === t
        end
    end

    intensity1 = [0.1 0.3 0.6; 0.2 0.4 0.1]
   labels1 = IndirectArray([1 2 3; 2 1 3], [RGB(1,0,0), RGB(0,0,1),RGB(0,1,0)])
   A1 = ColorizedArray(intensity1, labels1)
   target1 = intensity1 .* labels1
   @test eltype(A1) == RGB{Float64}
   @test size(A1) == (2,3)
   @test axes(A1) == (1:2, 1:3)
   for i = 1:6
       @test A1[i] === target1[i]
   end
   for j = 1:3, i = 1:2
       @test A1[i,j] === target1[i,j]
   end
   for (a,t) in zip(A1, target1)
       @test a === t
   end
end

nothing


text:
using Gadfly
using ProgressMeter

# Creates an array of length N where each entry is an independent realization
# of a complex Gaussian with Ez = Ezz = 0 and Ezz^* = 1.
function zGaussian(N=1)
   x, y = rand(N), rand(N)
   sqrt(-log(x)) .* exp(2im * pi * y)
end

function getRHS(nwells, T, KdJ)
   # Diagonal entries (0 for the action)
   center = (nwells + 1) / 2
   diags = push!(.5 * KdJ * ([1:nwells] - center).^2, 0)
   prop = -1im * diagm(diags)
   # Coupling between the wells
   coupls = push!(ones(nwells - 1) * T, 0)
   prop += 1im * (diagm(coupls, 1) + diagm(coupls, -1))

   (t, y) -> prop * y
end

function integrate_rk4(rhs, times, y0)
   y = zeros(typeof(y0[0]), (length(y0), length(times)))

   y[:, 1] = y0
   yrk = zeros(typeof(y0[0]), (length(y0), 4))

   for (i, t) in enumerate(times[2:end])
      dt = t - times[i]
      yrk[:, 1] = dt * rhs(t, y[:, i])
      yrk[:, 2] = dt * rhs(t + dt/2, y[:, i] + .5 * yrk[:, 1])
      yrk[:, 3] = dt * rhs(t + dt/2, y[:, i] + .5 * yrk[:, 2])
      yrk[:, 4] = dt * rhs(t + dt, y[:, i] + yrk[:, 3])

      y[:, i+1] = y[:, i] + 1/6 * (yrk[:, 1] + 2*yrk[:, 2] + 2*yrk[:, 3] + yrk[:, 4])
   end

   y
end

function main(timesteps, dt, nwells, realizations; T=2.0, KdJ=2.0)
   println("GO!")
   shift = convert(Int, floor((nwells - 1) // 2 % 2))
   even_w = filter(n -> isodd(n - shift), [1:nwells])
   odd_w = filter(n -> iseven(n - shift), [1:nwells])
   t = [1:timesteps] * dt
   rhs = getRHS(nwells, T, KdJ)
   norm = zeros(Complex{Float64}, timesteps)
   prog = Progress(realizations, 1)

   for n in [1:realizations]
      z1 = integrate_rk4(rhs, t, push!(zGaussian(nwells), 0))
      z2 = integrate_rk4(rhs, t, push!(zGaussian(nwells), 0))

      prefactor = prod(conj(z1[even_w, 1]) .* z2[even_w, 1])
      # Refactor as dot product
      sum_even = reshape(sum(z1[even_w, :] .* conj(z2[even_w, :]), 1), timesteps)
      sum_odd = reshape(sum(z1[odd_w, :] .* conj(z2[odd_w, :]), 1), timesteps)

      norm += prefactor * exp(sum_even + sum_odd) / realizations
      next!(prog)
   end

   println("DOne")
   t, norm
end

t, norm = main(100, 0.01, 5, 80000)
plot(x=t, y=real(norm), Geom.line)

max_stars_repo_name: dseuss/stuff

max_issues_repo_name: dseuss/stuff

max_forks_repo_name: dseuss/stuff


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

@test integrate((c+d*x^3)^4*(a+b*x^3), x) == :(a*x*c^4+1//4*c^3*x^4*(b*c+4*a*d)+1//13*d^3*x^13*(a*d+4*b*c)+1//16*b*d^4*x^16+1//5*c*d^2*x^10*(2*a*d+3*b*c)+2//7*d*c^2*x^7*(2*b*c+3*a*d))
@test integrate((c+d*x^3)^3*(a+b*x^3), x) == :(a*x*c^3+1//4*c^2*x^4*(b*c+3*a*d)+1//10*d^2*x^10*(a*d+3*b*c)+1//13*b*d^3*x^13+3//7*c*d*x^7*(a*d+b*c))
@test integrate((c+d*x^3)^2*(a+b*x^3), x) == :(a*x*c^2+1//4*c*x^4*(b*c+2*a*d)+1//7*d*x^7*(a*d+2*b*c)+1//10*b*d^2*x^10)
@test integrate((a+b*x^3)*(c+d*x^3), x) == :(x^4*(1//4*a*d+1//4*b*c)+a*c*x+1//7*b*d*x^7)
@test integrate((c+d*x^3)^-1*(a+b*x^3), x) == :(b*x*d^-1+-1//3*c^-2//3*d^-4//3*(b*c+-1*a*d)*log(c^1//3+x*d^1//3)+1//6*c^-2//3*d^-4//3*(b*c+-1*a*d)*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//3*3^(1/2)*c^-2//3*d^-4//3*(b*c+-1*a*d)*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3)))
@test integrate((c+d*x^3)^-2*(a+b*x^3), x) == :(-1//18*c^-5//3*d^-4//3*(b*c+2*a*d)*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//9*c^-5//3*d^-4//3*(b*c+2*a*d)*log(c^1//3+x*d^1//3)+-1//3*x*c^-1*d^-1*(c+d*x^3)^-1*(b*c+-1*a*d)+-1//9*3^(1/2)*c^-5//3*d^-4//3*(b*c+2*a*d)*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3)))
@test integrate((c+d*x^3)^-3*(a+b*x^3), x) == :(-1//54*c^-8//3*d^-4//3*(b*c+5*a*d)*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//27*c^-8//3*d^-4//3*(b*c+5*a*d)*log(c^1//3+x*d^1//3)+-1//6*x*c^-1*d^-1*(c+d*x^3)^-2*(b*c+-1*a*d)+-1//27*3^(1/2)*c^-8//3*d^-4//3*(b*c+5*a*d)*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3))+1//18*x*c^-2*d^-1*(c+d*x^3)^-1*(b*c+5*a*d))
@test integrate((a+b*x^3)^2*(c+d*x^3)^3, x) == :(x*a^2*c^3+1//7*c*x^7*(b^2*c^2+3*a^2*d^2+6*a*b*c*d)+1//10*d*x^10*(a^2*d^2+3*b^2*c^2+6*a*b*c*d)+1//16*b^2*d^3*x^16+1//4*a*c^2*x^4*(2*b*c+3*a*d)+1//13*b*d^2*x^13*(2*a*d+3*b*c))
@test integrate((a+b*x^3)^2*(c+d*x^3)^2, x) == :(x^7*(1//7*a^2*d^2+1//7*b^2*c^2+4//7*a*b*c*d)+x*a^2*c^2+1//13*b^2*d^2*x^13+(1/2)*a*c*x^4*(a*d+b*c)+1//5*b*d*x^10*(a*d+b*c))
@test integrate((a+b*x^3)^2*(c+d*x^3), x) == :(c*x*a^2+1//4*a*x^4*(a*d+2*b*c)+1//7*b*x^7*(b*c+2*a*d)+1//10*d*b^2*x^10)
@test integrate((a+b*x^3)^2*(c+d*x^3)^-1, x) == :(1//4*b^2*d^-1*x^4+-1*b*x*d^-2*(b*c+-2*a*d)+-1//6*c^-2//3*d^-7//3*(b*c+-1*a*d)^2*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//3*c^-2//3*d^-7//3*(b*c+-1*a*d)^2*log(c^1//3+x*d^1//3)+-1//3*3^(1/2)*c^-2//3*d^-7//3*(b*c+-1*a*d)^2*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3)))
@test integrate((a+b*x^3)^2*(c+d*x^3)^-2, x) == :(x*b^2*d^-2+-1//9*c^-5//3*d^-7//3*(a*d+2*b*c)*(-2*a*d+2*b*c)*log(c^1//3+x*d^1//3)+1//3*x*c^-1*d^-2*(c+d*x^3)^-1*(b*c+-1*a*d)^2+1//9*c^-5//3*d^-7//3*(a*d+2*b*c)*(b*c+-1*a*d)*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//9*3^(1/2)*c^-5//3*d^-7//3*(a*d+2*b*c)*(-2*a*d+2*b*c)*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3)))
@test integrate((a+b*x^3)^2*(c+d*x^3)^-3, x) == :(-1//54*c^-8//3*d^-7//3*(2*b^2*c^2+5*a^2*d^2+2*a*b*c*d)*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//27*c^-8//3*d^-7//3*(2*b^2*c^2+5*a^2*d^2+2*a*b*c*d)*log(c^1//3+x*d^1//3)+-1//27*3^(1/2)*c^-8//3*d^-7//3*(2*b^2*c^2+5*a^2*d^2+2*a*b*c*d)*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3))+-1//6*x*c^-1*d^-1*(c+d*x^3)^-2*(a+b*x^3)*(b*c+-1*a*d)+-1//18*x*c^-2*d^-2*(c+d*x^3)^-1*(b*c+-1*a*d)*(4*b*c+5*a*d))
@test integrate((a+b*x^3)^-1*(c+d*x^3)^4, x) == :(1//10*b^-1*d^4*x^10+-1//6*a^-2//3*b^-13//3*(b*c+-1*a*d)^4*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*a^-2//3*b^-13//3*(b*c+-1*a*d)^4*log(a^1//3+x*b^1//3)+1//4*b^-3*d^2*x^4*(a^2*d^2+6*b^2*c^2+-4*a*b*c*d)+1//7*b^-2*d^3*x^7*(-1*a*d+4*b*c)+d*x*b^-4*(-1*a*d+2*b*c)*(a^2*d^2+2*b^2*c^2+-2*a*b*c*d)+-1//3*3^(1/2)*a^-2//3*b^-13//3*(b*c+-1*a*d)^4*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-1*(c+d*x^3)^3, x) == :(1//7*b^-1*d^3*x^7+d*x*b^-3*(a^2*d^2+3*b^2*c^2+-3*a*b*c*d)+-1//6*a^-2//3*b^-10//3*(b*c+-1*a*d)^3*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*a^-2//3*b^-10//3*(b*c+-1*a*d)^3*log(a^1//3+x*b^1//3)+1//4*b^-2*d^2*x^4*(-1*a*d+3*b*c)+-1//3*3^(1/2)*a^-2//3*b^-10//3*(b*c+-1*a*d)^3*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-1*(c+d*x^3)^2, x) == :(1//4*b^-1*d^2*x^4+d*x*b^-2*(-1*a*d+2*b*c)+-1//6*a^-2//3*b^-7//3*(b*c+-1*a*d)^2*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*a^-2//3*b^-7//3*(b*c+-1*a*d)^2*log(a^1//3+x*b^1//3)+-1//3*3^(1/2)*a^-2//3*b^-7//3*(b*c+-1*a*d)^2*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-1*(c+d*x^3), x) == :(d*x*b^-1+-1//6*a^-2//3*b^-4//3*(b*c+-1*a*d)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*a^-2//3*b^-4//3*(b*c+-1*a*d)*log(a^1//3+x*b^1//3)+-1//3*3^(1/2)*a^-2//3*b^-4//3*(b*c+-1*a*d)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-1*(c+d*x^3)^-1, x) == :(-1//3*c^-2//3*d^2//3*(b*c+-1*a*d)^-1*log(c^1//3+x*d^1//3)+-1//6*a^-2//3*b^2//3*(b*c+-1*a*d)^-1*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*a^-2//3*b^2//3*(b*c+-1*a*d)^-1*log(a^1//3+x*b^1//3)+1//6*c^-2//3*d^2//3*(b*c+-1*a*d)^-1*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+-1//3*3^(1/2)*a^-2//3*b^2//3*(b*c+-1*a*d)^-1*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3))+1//3*3^(1/2)*c^-2//3*d^2//3*(b*c+-1*a*d)^-1*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3)))
@test integrate((a+b*x^3)^-1*(c+d*x^3)^-2, x) == :(-1//6*a^-2//3*b^5//3*(b*c+-1*a*d)^-2*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*a^-2//3*b^5//3*(b*c+-1*a*d)^-2*log(a^1//3+x*b^1//3)+-1//3*d*x*c^-1*(c+d*x^3)^-1*(b*c+-1*a*d)^-1+-1//3*3^(1/2)*a^-2//3*b^5//3*(b*c+-1*a*d)^-2*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3))+-1//9*c^-5//3*d^2//3*(b*c+-1*a*d)^-2*(-2*a*d+5*b*c)*log(c^1//3+x*d^1//3)+1//18*c^-5//3*d^2//3*(b*c+-1*a*d)^-2*(-2*a*d+5*b*c)*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//9*3^(1/2)*c^-5//3*d^2//3*(b*c+-1*a*d)^-2*(-2*a*d+5*b*c)*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3)))
@test integrate((a+b*x^3)^-2*(c+d*x^3)^5, x) == :(1//10*b^-2*d^5*x^10+x*b^-5*d^2*(-4*a^3*d^3+10*b^3*c^3+-20*a*d*b^2*c^2+15*b*c*a^2*d^2)+1//4*b^-4*d^3*x^4*(3*a^2*d^2+10*b^2*c^2+-10*a*b*c*d)+1//7*b^-3*d^4*x^7*(-2*a*d+5*b*c)+-1//18*a^-5//3*b^-16//3*(b*c+-1*a*d)^4*(2*b*c+13*a*d)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*x*a^-1*b^-5*(a+b*x^3)^-1*(b*c+-1*a*d)^5+1//9*a^-5//3*b^-16//3*(b*c+-1*a*d)^4*(2*b*c+13*a*d)*log(a^1//3+x*b^1//3)+-1//9*3^(1/2)*a^-5//3*b^-16//3*(b*c+-1*a*d)^4*(2*b*c+13*a*d)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-2*(c+d*x^3)^4, x) == :(1//7*b^-2*d^4*x^7+x*b^-4*d^2*(3*a^2*d^2+6*b^2*c^2+-8*a*b*c*d)+(1/2)*b^-3*d^3*x^4*(-1*a*d+2*b*c)+-1//9*a^-5//3*b^-13//3*(b*c+-1*a*d)^3*(b*c+5*a*d)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*x*a^-1*b^-4*(a+b*x^3)^-1*(b*c+-1*a*d)^4+2//9*a^-5//3*b^-13//3*(b*c+-1*a*d)^3*(b*c+5*a*d)*log(a^1//3+x*b^1//3)+-2//9*3^(1/2)*a^-5//3*b^-13//3*(b*c+-1*a*d)^3*(b*c+5*a*d)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-2*(c+d*x^3)^3, x) == :(1//4*b^-2*d^3*x^4+x*b^-3*d^2*(-2*a*d+3*b*c)+-1//18*a^-5//3*b^-10//3*(b*c+-1*a*d)^2*(2*b*c+7*a*d)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*x*a^-1*b^-3*(a+b*x^3)^-1*(b*c+-1*a*d)^3+1//9*a^-5//3*b^-10//3*(b*c+-1*a*d)^2*(2*b*c+7*a*d)*log(a^1//3+x*b^1//3)+-1//9*3^(1/2)*a^-5//3*b^-10//3*(b*c+-1*a*d)^2*(2*b*c+7*a*d)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-2*(c+d*x^3)^2, x) == :(x*b^-2*d^2+-1//9*a^-5//3*b^-7//3*(b*c+-1*a*d)*(b*c+2*a*d)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*x*a^-1*b^-2*(a+b*x^3)^-1*(b*c+-1*a*d)^2+1//9*a^-5//3*b^-7//3*(b*c+2*a*d)*(-2*a*d+2*b*c)*log(a^1//3+x*b^1//3)+-1//9*3^(1/2)*a^-5//3*b^-7//3*(b*c+2*a*d)*(-2*a*d+2*b*c)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-2*(c+d*x^3), x) == :(-1//18*a^-5//3*b^-4//3*(a*d+2*b*c)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//9*a^-5//3*b^-4//3*(a*d+2*b*c)*log(a^1//3+x*b^1//3)+-1//9*3^(1/2)*a^-5//3*b^-4//3*(a*d+2*b*c)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3))+1//3*x*a^-1*b^-1*(a+b*x^3)^-1*(b*c+-1*a*d))
@test integrate((a+b*x^3)^-2*(c+d*x^3)^-1, x) == :(-1//6*c^-2//3*d^5//3*(b*c+-1*a*d)^-2*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+1//3*c^-2//3*d^5//3*(b*c+-1*a*d)^-2*log(c^1//3+x*d^1//3)+-1//3*3^(1/2)*c^-2//3*d^5//3*(b*c+-1*a*d)^-2*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3))+-1//18*a^-5//3*b^2//3*(b*c+-1*a*d)^-2*(-5*a*d+2*b*c)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+1//3*b*x*a^-1*(a+b*x^3)^-1*(b*c+-1*a*d)^-1+1//9*a^-5//3*b^2//3*(b*c+-1*a*d)^-2*(-5*a*d+2*b*c)*log(a^1//3+x*b^1//3)+-1//9*3^(1/2)*a^-5//3*b^2//3*(b*c+-1*a*d)^-2*(-5*a*d+2*b*c)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3)))
@test integrate((a+b*x^3)^-2*(c+d*x^3)^-2, x) == :(-1//9*a^-5//3*b^5//3*(b*c+-1*a*d)^-3*(b*c+-4*a*d)*log(a^2//3+b^2//3*x^2+-1*x*a^1//3*b^1//3)+-1//9*c^-5//3*d^5//3*(b*c+-1*a*d)^-3*(-1*a*d+4*b*c)*log(c^2//3+d^2//3*x^2+-1*x*c^1//3*d^1//3)+2//9*a^-5//3*b^5//3*(b*c+-1*a*d)^-3*(b*c+-4*a*d)*log(a^1//3+x*b^1//3)+2//9*c^-5//3*d^5//3*(b*c+-1*a*d)^-3*(-1*a*d+4*b*c)*log(c^1//3+x*d^1//3)+-2//9*3^(1/2)*a^-5//3*b^5//3*(b*c+-1*a*d)^-3*(b*c+-4*a*d)*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*b^1//3))+-2//9*3^(1/2)*c^-5//3*d^5//3*(b*c+-1*a*d)^-3*(-1*a*d+4*b*c)*arctan(1//3*3^(1/2)*c^-1//3*(c^1//3+-2*x*d^1//3))+1//3*b*x*a^-1*(a+b*x^3)^-1*(c+d*x^3)^-1*(b*c+-1*a*d)^-1+1//3*d*x*a^-1*c^-1*(c+d*x^3)^-1*(b*c+-1*a*d)^-2*(a*d+b*c))
@test integrate((a+b*x^3)^3*(c+d*x^3)^-13//3, x) == :(1//10*x*c^-1*(a+b*x^3)^3*(c+d*x^3)^-10//3+81//140*x*a^3*c^-4*(c+d*x^3)^-1//3+9//70*a*x*c^-2*(a+b*x^3)^2*(c+d*x^3)^-7//3+27//140*x*a^2*c^-3*(c+d*x^3)^-4//3*(a+b*x^3))
@test integrate((a+b*x^3)^2*(c+d*x^3)^-10//3, x) == :(1//7*x*c^-1*(a+b*x^3)^2*(c+d*x^3)^-7//3+9//14*x*a^2*c^-3*(c+d*x^3)^-1//3+3//14*a*x*c^-2*(c+d*x^3)^-4//3*(a+b*x^3))
@test integrate((c+d*x^3)^-7//3*(a+b*x^3), x) == :(1//4*x*c^-1*(c+d*x^3)^-4//3*(a+b*x^3)+3//4*a*x*c^-2*(c+d*x^3)^-1//3)
@test integrate((c+d*x^3)^-4//3, x) == :(x*c^-1*(c+d*x^3)^-1//3)
@test integrate((a+b*x^3)^-1*(c+d*x^3)^-1//3, x) == :(-1//6*a^-2//3*(b*c+-1*a*d)^-1//3*log(a^2//3+x^2*(c+d*x^3)^-2//3*(b*c+-1*a*d)^2//3+-1*x*a^1//3*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)+1//3*a^-2//3*(b*c+-1*a*d)^-1//3*log(a^1//3+x*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)+-1//3*3^(1/2)*a^-2//3*(b*c+-1*a*d)^-1//3*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)))
@test integrate((a+b*x^3)^-2*(c+d*x^3)^2//3, x) == :(-1//9*c*a^-5//3*(b*c+-1*a*d)^-1//3*log(a^2//3+x^2*(c+d*x^3)^-2//3*(b*c+-1*a*d)^2//3+-1*x*a^1//3*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)+1//3*x*a^-1*(a+b*x^3)^-1*(c+d*x^3)^2//3+2//9*c*a^-5//3*(b*c+-1*a*d)^-1//3*log(a^1//3+x*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)+-2//9*c*3^(1/2)*a^-5//3*(b*c+-1*a*d)^-1//3*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)))
@test integrate((a+b*x^3)^-3*(c+d*x^3)^5//3, x) == :(-5//54*a^-8//3*c^2*(b*c+-1*a*d)^-1//3*log(a^2//3+x^2*(c+d*x^3)^-2//3*(b*c+-1*a*d)^2//3+-1*x*a^1//3*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)+1//6*x*a^-1*(a+b*x^3)^-2*(c+d*x^3)^5//3+5//27*a^-8//3*c^2*(b*c+-1*a*d)^-1//3*log(a^1//3+x*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3)+-5//27*3^(1/2)*a^-8//3*c^2*(b*c+-1*a*d)^-1//3*arctan(1//3*3^(1/2)*a^-1//3*(a^1//3+-2*x*(c+d*x^3)^-1//3*(b*c+-1*a*d)^1//3))+5//18*c*x*a^-2*(a+b*x^3)^-1*(c+d*x^3)^2//3)
@test integrate((a+b*x^3)^(-1+-1*b*c*(-3*a*d+3*b*c)^-1)*(c+d*x^3)^(-1+a*d*(-3*a*d+3*b*c)^-1), x) == :(x*a^-1*c^-1*(a+b*x^3)^(-1*b*c*(-3*a*d+3*b*c)^-1)*(c+d*x^3)^(a*d*(-3*a*d+3*b*c)^-1))
@test integrate((c+d*x^4)^4*(a+b*x^4), x) == :(a*x*c^4+1//5*c^3*x^5*(b*c+4*a*d)+1//17*d^3*x^17*(a*d+4*b*c)+1//21*b*d^4*x^21+2//9*d*c^2*x^9*(2*b*c+3*a*d)+2//13*c*d^2*x^13*(2*a*d+3*b*c))
@test integrate((c+d*x^4)^3*(a+b*x^4), x) == :(a*x*c^3+1//5*c^2*x^5*(b*c+3*a*d)+1//13*d^2*x^13*(a*d+3*b*c)+1//17*b*d^3*x^17+1//3*c*d*x^9*(a*d+b*c))
@test integrate((c+d*x^4)^2*(a+b*x^4), x) == :(a*x*c^2+1//5*c*x^5*(b*c+2*a*d)+1//9*d*x^9*(a*d+2*b*c)+1//13*b*d^2*x^13)
@test integrate((a+b*x^4)*(c+d*x^4), x) == :(x^5*(1//5*a*d+1//5*b*c)+a*c*x+1//9*b*d*x^9)
@test integrate((c+d*x^4)^-1*(a+b*x^4), x) == :(b*x*d^-1+-1//4*2^(1/2)*c^-3//4*d^-5//4*(b*c+-1*a*d)*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+-1//8*2^(1/2)*c^-3//4*d^-5//4*(b*c+-1*a*d)*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4)+1//4*2^(1/2)*c^-3//4*d^-5//4*(b*c+-1*a*d)*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+1//8*2^(1/2)*c^-3//4*d^-5//4*(b*c+-1*a*d)*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4))
@test integrate((c+d*x^4)^-2*(a+b*x^4), x) == :(-1//4*x*c^-1*d^-1*(c+d*x^4)^-1*(b*c+-1*a*d)+-1//16*2^(1/2)*c^-7//4*d^-5//4*(b*c+3*a*d)*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+-1//32*2^(1/2)*c^-7//4*d^-5//4*(b*c+3*a*d)*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4)+1//16*2^(1/2)*c^-7//4*d^-5//4*(b*c+3*a*d)*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+1//32*2^(1/2)*c^-7//4*d^-5//4*(b*c+3*a*d)*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4))
@test integrate((c+d*x^4)^-3*(a+b*x^4), x) == :(-1//8*x*c^-1*d^-1*(c+d*x^4)^-2*(b*c+-1*a*d)+-1//128*2^(1/2)*c^-11//4*d^-5//4*(3*b*c+21*a*d)*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+-1//256*2^(1/2)*c^-11//4*d^-5//4*(3*b*c+21*a*d)*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4)+1//32*x*c^-2*d^-1*(c+d*x^4)^-1*(b*c+7*a*d)+1//128*2^(1/2)*c^-11//4*d^-5//4*(3*b*c+21*a*d)*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+1//256*2^(1/2)*c^-11//4*d^-5//4*(3*b*c+21*a*d)*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4))
@test integrate((a+b*x^4)^2*(c+d*x^4)^4, x) == :(x*a^2*c^4+1//9*c^2*x^9*(b^2*c^2+6*a^2*d^2+8*a*b*c*d)+1//17*d^2*x^17*(a^2*d^2+6*b^2*c^2+8*a*b*c*d)+1//25*b^2*d^4*x^25+2//5*a*c^3*x^5*(b*c+2*a*d)+2//21*b*d^3*x^21*(a*d+2*b*c)+4//13*c*d*x^13*(a^2*d^2+b^2*c^2+3*a*b*c*d))
@test integrate((a+b*x^4)^2*(c+d*x^4)^3, x) == :(x*a^2*c^3+1//9*c*x^9*(b^2*c^2+3*a^2*d^2+6*a*b*c*d)+1//13*d*x^13*(a^2*d^2+3*b^2*c^2+6*a*b*c*d)+1//21*b^2*d^3*x^21+1//5*a*c^2*x^5*(2*b*c+3*a*d)+1//17*b*d^2*x^17*(2*a*d+3*b*c))
@test integrate((a+b*x^4)^2*(c+d*x^4)^2, x) == :(x^9*(1//9*a^2*d^2+1//9*b^2*c^2+4//9*a*b*c*d)+x*a^2*c^2+1//17*b^2*d^2*x^17+2//5*a*c*x^5*(a*d+b*c)+2//13*b*d*x^13*(a*d+b*c))
@test integrate((a+b*x^4)^2*(c+d*x^4), x) == :(c*x*a^2+1//5*a*x^5*(a*d+2*b*c)+1//9*b*x^9*(b*c+2*a*d)+1//13*d*b^2*x^13)
@test integrate((a+b*x^4)^2*(c+d*x^4)^-1, x) == :(1//5*b^2*d^-1*x^5+-1*b*x*d^-2*(b*c+-2*a*d)+-1//4*2^(1/2)*c^-3//4*d^-9//4*(b*c+-1*a*d)^2*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+-1//8*2^(1/2)*c^-3//4*d^-9//4*(b*c+-1*a*d)^2*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4)+1//4*2^(1/2)*c^-3//4*d^-9//4*(b*c+-1*a*d)^2*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+1//8*2^(1/2)*c^-3//4*d^-9//4*(b*c+-1*a*d)^2*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4))
@test integrate((a+b*x^4)^2*(c+d*x^4)^-2, x) == :(x*b^2*d^-2+1//4*x*c^-1*d^-2*(c+d*x^4)^-1*(b*c+-1*a*d)^2+-1//16*2^(1/2)*c^-7//4*d^-9//4*(b*c+-1*a*d)*(3*a*d+5*b*c)*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+-1//32*2^(1/2)*c^-7//4*d^-9//4*(b*c+-1*a*d)*(3*a*d+5*b*c)*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4)+1//16*2^(1/2)*c^-7//4*d^-9//4*(b*c+-1*a*d)*(3*a*d+5*b*c)*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+1//32*2^(1/2)*c^-7//4*d^-9//4*(b*c+-1*a*d)*(3*a*d+5*b*c)*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4))
@test integrate((a+b*x^4)^2*(c+d*x^4)^-3, x) == :(-1//128*2^(1/2)*c^-11//4*d^-9//4*(5*b^2*c^2+21*a^2*d^2+6*a*b*c*d)*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+-1//256*2^(1/2)*c^-11//4*d^-9//4*(5*b^2*c^2+21*a^2*d^2+6*a*b*c*d)*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4)+1//128*2^(1/2)*c^-11//4*d^-9//4*(5*b^2*c^2+21*a^2*d^2+6*a*b*c*d)*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+1//256*2^(1/2)*c^-11//4*d^-9//4*(5*b^2*c^2+21*a^2*d^2+6*a*b*c*d)*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4)+-1//8*x*c^-1*d^-1*(c+d*x^4)^-2*(a+b*x^4)*(b*c+-1*a*d)+-1//32*x*c^-2*d^-2*(c+d*x^4)^-1*(b*c+-1*a*d)*(5*b*c+7*a*d))
@test integrate((a+b*x^4)^-1*(c+d*x^4)^4, x) == :(1//13*b^-1*d^4*x^13+1//5*b^-3*d^2*x^5*(a^2*d^2+6*b^2*c^2+-4*a*b*c*d)+1//9*b^-2*d^3*x^9*(-1*a*d+4*b*c)+d*x*b^-4*(-1*a*d+2*b*c)*(a^2*d^2+2*b^2*c^2+-2*a*b*c*d)+-1//4*2^(1/2)*a^-3//4*b^-17//4*(b*c+-1*a*d)^4*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//8*2^(1/2)*a^-3//4*b^-17//4*(b*c+-1*a*d)^4*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//4*2^(1/2)*a^-3//4*b^-17//4*(b*c+-1*a*d)^4*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//8*2^(1/2)*a^-3//4*b^-17//4*(b*c+-1*a*d)^4*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-1*(c+d*x^4)^3, x) == :(1//9*b^-1*d^3*x^9+d*x*b^-3*(a^2*d^2+3*b^2*c^2+-3*a*b*c*d)+1//5*b^-2*d^2*x^5*(-1*a*d+3*b*c)+-1//4*2^(1/2)*a^-3//4*b^-13//4*(b*c+-1*a*d)^3*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//8*2^(1/2)*a^-3//4*b^-13//4*(b*c+-1*a*d)^3*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//4*2^(1/2)*a^-3//4*b^-13//4*(b*c+-1*a*d)^3*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//8*2^(1/2)*a^-3//4*b^-13//4*(b*c+-1*a*d)^3*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-1*(c+d*x^4)^2, x) == :(1//5*b^-1*d^2*x^5+d*x*b^-2*(-1*a*d+2*b*c)+-1//4*2^(1/2)*a^-3//4*b^-9//4*(b*c+-1*a*d)^2*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//8*2^(1/2)*a^-3//4*b^-9//4*(b*c+-1*a*d)^2*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//4*2^(1/2)*a^-3//4*b^-9//4*(b*c+-1*a*d)^2*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//8*2^(1/2)*a^-3//4*b^-9//4*(b*c+-1*a*d)^2*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-1*(c+d*x^4), x) == :(d*x*b^-1+-1//4*2^(1/2)*a^-3//4*b^-5//4*(b*c+-1*a*d)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//8*2^(1/2)*a^-3//4*b^-5//4*(b*c+-1*a*d)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//4*2^(1/2)*a^-3//4*b^-5//4*(b*c+-1*a*d)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//8*2^(1/2)*a^-3//4*b^-5//4*(b*c+-1*a*d)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-1*(c+d*x^4)^-1, x) == :(-1//4*2^(1/2)*a^-3//4*b^3//4*(b*c+-1*a*d)^-1*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//4*2^(1/2)*c^-3//4*d^3//4*(b*c+-1*a*d)^-1*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+-1//8*2^(1/2)*a^-3//4*b^3//4*(b*c+-1*a*d)^-1*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+-1//8*2^(1/2)*c^-3//4*d^3//4*(b*c+-1*a*d)^-1*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4)+1//4*2^(1/2)*a^-3//4*b^3//4*(b*c+-1*a*d)^-1*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//4*2^(1/2)*c^-3//4*d^3//4*(b*c+-1*a*d)^-1*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+1//8*2^(1/2)*a^-3//4*b^3//4*(b*c+-1*a*d)^-1*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4)+1//8*2^(1/2)*c^-3//4*d^3//4*(b*c+-1*a*d)^-1*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4))
@test integrate((a+b*x^4)^-1*(c+d*x^4)^-2, x) == :(-1//4*d*x*c^-1*(c+d*x^4)^-1*(b*c+-1*a*d)^-1+-1//4*2^(1/2)*a^-3//4*b^7//4*(b*c+-1*a*d)^-2*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//8*2^(1/2)*a^-3//4*b^7//4*(b*c+-1*a*d)^-2*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//4*2^(1/2)*a^-3//4*b^7//4*(b*c+-1*a*d)^-2*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//8*2^(1/2)*a^-3//4*b^7//4*(b*c+-1*a*d)^-2*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4)+-1//16*2^(1/2)*c^-7//4*d^3//4*(b*c+-1*a*d)^-2*(-3*a*d+7*b*c)*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+-1//32*2^(1/2)*c^-7//4*d^3//4*(b*c+-1*a*d)^-2*(-3*a*d+7*b*c)*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4)+1//16*2^(1/2)*c^-7//4*d^3//4*(b*c+-1*a*d)^-2*(-3*a*d+7*b*c)*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+1//32*2^(1/2)*c^-7//4*d^3//4*(b*c+-1*a*d)^-2*(-3*a*d+7*b*c)*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4))
@test integrate((a+b*x^4)^-2*(c+d*x^4)^5, x) == :(1//13*b^-2*d^5*x^13+x*b^-5*d^2*(-4*a^3*d^3+10*b^3*c^3+-20*a*d*b^2*c^2+15*b*c*a^2*d^2)+1//5*b^-4*d^3*x^5*(3*a^2*d^2+10*b^2*c^2+-10*a*b*c*d)+1//9*b^-3*d^4*x^9*(-2*a*d+5*b*c)+1//4*x*a^-1*b^-5*(a+b*x^4)^-1*(b*c+-1*a*d)^5+-1//16*2^(1/2)*a^-7//4*b^-21//4*(b*c+-1*a*d)^4*(3*b*c+17*a*d)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//32*2^(1/2)*a^-7//4*b^-21//4*(b*c+-1*a*d)^4*(3*b*c+17*a*d)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//16*2^(1/2)*a^-7//4*b^-21//4*(b*c+-1*a*d)^4*(3*b*c+17*a*d)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//32*2^(1/2)*a^-7//4*b^-21//4*(b*c+-1*a*d)^4*(3*b*c+17*a*d)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-2*(c+d*x^4)^4, x) == :(1//9*b^-2*d^4*x^9+x*b^-4*d^2*(3*a^2*d^2+6*b^2*c^2+-8*a*b*c*d)+2//5*b^-3*d^3*x^5*(-1*a*d+2*b*c)+1//4*x*a^-1*b^-4*(a+b*x^4)^-1*(b*c+-1*a*d)^4+-1//16*2^(1/2)*a^-7//4*b^-17//4*(b*c+-1*a*d)^3*(3*b*c+13*a*d)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//32*2^(1/2)*a^-7//4*b^-17//4*(b*c+-1*a*d)^3*(3*b*c+13*a*d)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//16*2^(1/2)*a^-7//4*b^-17//4*(b*c+-1*a*d)^3*(3*b*c+13*a*d)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//32*2^(1/2)*a^-7//4*b^-17//4*(b*c+-1*a*d)^3*(3*b*c+13*a*d)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-2*(c+d*x^4)^3, x) == :(1//5*b^-2*d^3*x^5+x*b^-3*d^2*(-2*a*d+3*b*c)+1//4*x*a^-1*b^-3*(a+b*x^4)^-1*(b*c+-1*a*d)^3+-3//16*2^(1/2)*a^-7//4*b^-13//4*(b*c+-1*a*d)^2*(b*c+3*a*d)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-3//32*2^(1/2)*a^-7//4*b^-13//4*(b*c+-1*a*d)^2*(b*c+3*a*d)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+3//16*2^(1/2)*a^-7//4*b^-13//4*(b*c+-1*a*d)^2*(b*c+3*a*d)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+3//32*2^(1/2)*a^-7//4*b^-13//4*(b*c+-1*a*d)^2*(b*c+3*a*d)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-2*(c+d*x^4)^2, x) == :(x*b^-2*d^2+1//4*x*a^-1*b^-2*(a+b*x^4)^-1*(b*c+-1*a*d)^2+-1//16*2^(1/2)*a^-7//4*b^-9//4*(b*c+-1*a*d)*(3*b*c+5*a*d)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//32*2^(1/2)*a^-7//4*b^-9//4*(b*c+-1*a*d)*(3*b*c+5*a*d)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//16*2^(1/2)*a^-7//4*b^-9//4*(b*c+-1*a*d)*(3*b*c+5*a*d)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//32*2^(1/2)*a^-7//4*b^-9//4*(b*c+-1*a*d)*(3*b*c+5*a*d)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-2*(c+d*x^4), x) == :(-1//16*2^(1/2)*a^-7//4*b^-5//4*(a*d+3*b*c)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//32*2^(1/2)*a^-7//4*b^-5//4*(a*d+3*b*c)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//4*x*a^-1*b^-1*(a+b*x^4)^-1*(b*c+-1*a*d)+1//16*2^(1/2)*a^-7//4*b^-5//4*(a*d+3*b*c)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//32*2^(1/2)*a^-7//4*b^-5//4*(a*d+3*b*c)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-2*(c+d*x^4)^-1, x) == :(-1//4*2^(1/2)*c^-3//4*d^7//4*(b*c+-1*a*d)^-2*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+-1//8*2^(1/2)*c^-3//4*d^7//4*(b*c+-1*a*d)^-2*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4)+1//4*b*x*a^-1*(a+b*x^4)^-1*(b*c+-1*a*d)^-1+1//4*2^(1/2)*c^-3//4*d^7//4*(b*c+-1*a*d)^-2*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+1//8*2^(1/2)*c^-3//4*d^7//4*(b*c+-1*a*d)^-2*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4)+-1//16*2^(1/2)*a^-7//4*b^3//4*(b*c+-1*a*d)^-2*(-7*a*d+3*b*c)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//32*2^(1/2)*a^-7//4*b^3//4*(b*c+-1*a*d)^-2*(-7*a*d+3*b*c)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+1//16*2^(1/2)*a^-7//4*b^3//4*(b*c+-1*a*d)^-2*(-7*a*d+3*b*c)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//32*2^(1/2)*a^-7//4*b^3//4*(b*c+-1*a*d)^-2*(-7*a*d+3*b*c)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4))
@test integrate((a+b*x^4)^-2*(c+d*x^4)^-2, x) == :(-1//16*2^(1/2)*a^-7//4*b^7//4*(b*c+-1*a*d)^-3*(-11*a*d+3*b*c)*arctan(1+-1*x*2^(1/2)*a^-1//4*b^1//4)+-1//16*2^(1/2)*c^-7//4*d^7//4*(b*c+-1*a*d)^-3*(-3*a*d+11*b*c)*arctan(1+-1*x*2^(1/2)*c^-1//4*d^1//4)+-1//32*2^(1/2)*a^-7//4*b^7//4*(b*c+-1*a*d)^-3*(-11*a*d+3*b*c)*log(a^(1/2)+b^(1/2)*x^2+-1*x*2^(1/2)*a^1//4*b^1//4)+-1//32*2^(1/2)*c^-7//4*d^7//4*(b*c+-1*a*d)^-3*(-3*a*d+11*b*c)*log(c^(1/2)+d^(1/2)*x^2+-1*x*2^(1/2)*c^1//4*d^1//4)+1//4*b*x*a^-1*(a+b*x^4)^-1*(c+d*x^4)^-1*(b*c+-1*a*d)^-1+1//16*2^(1/2)*a^-7//4*b^7//4*(b*c+-1*a*d)^-3*(-11*a*d+3*b*c)*arctan(1+x*2^(1/2)*a^-1//4*b^1//4)+1//16*2^(1/2)*c^-7//4*d^7//4*(b*c+-1*a*d)^-3*(-3*a*d+11*b*c)*arctan(1+x*2^(1/2)*c^-1//4*d^1//4)+1//32*2^(1/2)*a^-7//4*b^7//4*(b*c+-1*a*d)^-3*(-11*a*d+3*b*c)*log(a^(1/2)+b^(1/2)*x^2+x*2^(1/2)*a^1//4*b^1//4)+1//32*2^(1/2)*c^-7//4*d^7//4*(b*c+-1*a*d)^-3*(-3*a*d+11*b*c)*log(c^(1/2)+d^(1/2)*x^2+x*2^(1/2)*c^1//4*d^1//4)+1//4*d*x*a^-1*c^-1*(c+d*x^4)^-1*(b*c+-1*a*d)^-2*(a*d+b*c))
@test integrate((a+b*x^4)^(1/2)*(a*c+-1*b*c*x^4)^-1, x) == :(1//4*2^(1/2)*a^-1//4*b^-1//4*c^-1*arctan(x*2^(1/2)*a^1//4*b^1//4*(a+b*x^4)^-1//2)+1//4*2^(1/2)*a^-1//4*b^-1//4*c^-1*arctanh(x*2^(1/2)*a^1//4*b^1//4*(a+b*x^4)^-1//2))
@test integrate((a+-1*b*x^4)^(1/2)*(a*c+b*c*x^4)^-1, x) == :((1/2)*a^-1//4*b^-1//4*c^-1*arctan(x*a^-1//4*b^1//4*(a+-1*b*x^4)^-1//2*(a^(1/2)+b^(1/2)*x^2))+(1/2)*a^-1//4*b^-1//4*c^-1*arctanh(x*a^-1//4*b^1//4*(a+-1*b*x^4)^-1//2*(a^(1/2)+-1*b^(1/2)*x^2)))
@test integrate((a+b*x^4)^7//4*(c+d*x^4)^-1, x) == :((1/2)*c^-3//4*d^-2*(b*c+-1*a*d)^7//4*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+(1/2)*c^-3//4*d^-2*(b*c+-1*a*d)^7//4*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//8*b^3//4*d^-2*(-7*a*d+4*b*c)*arctan(x*b^1//4*(a+b*x^4)^-1//4)+-1//8*b^3//4*d^-2*(-7*a*d+4*b*c)*arctanh(x*b^1//4*(a+b*x^4)^-1//4)+1//4*b*x*d^-1*(a+b*x^4)^3//4)
@test integrate((a+b*x^4)^3//4*(c+d*x^4)^-1, x) == :((1/2)*b^3//4*d^-1*arctan(x*b^1//4*(a+b*x^4)^-1//4)+(1/2)*b^3//4*d^-1*arctanh(x*b^1//4*(a+b*x^4)^-1//4)+-1//2*c^-3//4*d^-1*(b*c+-1*a*d)^3//4*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//2*c^-3//4*d^-1*(b*c+-1*a*d)^3//4*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4))
@test integrate((a+b*x^4)^-1//4*(c+d*x^4)^-1, x) == :((1/2)*c^-3//4*(b*c+-1*a*d)^-1//4*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+(1/2)*c^-3//4*(b*c+-1*a*d)^-1//4*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4))
@test integrate((a+b*x^4)^-5//4*(c+d*x^4)^-1, x) == :(-1//2*d*c^-3//4*(b*c+-1*a*d)^-5//4*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//2*d*c^-3//4*(b*c+-1*a*d)^-5//4*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+b*x*a^-1*(a+b*x^4)^-1//4*(b*c+-1*a*d)^-1)
@test integrate((a+b*x^4)^-9//4*(c+d*x^4)^-1, x) == :((1/2)*c^-3//4*d^2*(b*c+-1*a*d)^-9//4*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+(1/2)*c^-3//4*d^2*(b*c+-1*a*d)^-9//4*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+1//5*b*x*a^-1*(a+b*x^4)^-5//4*(b*c+-1*a*d)^-1+1//5*b*x*a^-2*(a+b*x^4)^-1//4*(b*c+-1*a*d)^-2*(-9*a*d+4*b*c))
@test integrate((a+b*x^4)^-13//4*(c+d*x^4)^-1, x) == :(-1//2*c^-3//4*d^3*(b*c+-1*a*d)^-13//4*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//2*c^-3//4*d^3*(b*c+-1*a*d)^-13//4*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+1//9*b*x*a^-1*(a+b*x^4)^-9//4*(b*c+-1*a*d)^-1+1//45*b*x*a^-3*(a+b*x^4)^-1//4*(b*c+-1*a*d)^-3*(32*b^2*c^2+113*a^2*d^2+-100*a*b*c*d)+1//45*b*x*a^-2*(a+b*x^4)^-5//4*(b*c+-1*a*d)^-2*(-17*a*d+8*b*c))
@test integrate((a+b*x^4)^11//4*(c+d*x^4)^-2, x) == :(-1//8*b^7//4*d^-3*(-11*a*d+8*b*c)*arctan(x*b^1//4*(a+b*x^4)^-1//4)+-1//8*b^7//4*d^-3*(-11*a*d+8*b*c)*arctanh(x*b^1//4*(a+b*x^4)^-1//4)+1//8*c^-7//4*d^-3*(b*c+-1*a*d)^7//4*(3*a*d+8*b*c)*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+1//8*c^-7//4*d^-3*(b*c+-1*a*d)^7//4*(3*a*d+8*b*c)*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//4*x*c^-1*d^-1*(a+b*x^4)^7//4*(c+d*x^4)^-1*(b*c+-1*a*d)+1//4*b*x*c^-1*d^-2*(a+b*x^4)^3//4*(-1*a*d+2*b*c))
@test integrate((a+b*x^4)^7//4*(c+d*x^4)^-2, x) == :((1/2)*b^7//4*d^-2*arctan(x*b^1//4*(a+b*x^4)^-1//4)+(1/2)*b^7//4*d^-2*arctanh(x*b^1//4*(a+b*x^4)^-1//4)+-1//8*c^-7//4*d^-2*(b*c+-1*a*d)^3//4*(3*a*d+4*b*c)*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//8*c^-7//4*d^-2*(b*c+-1*a*d)^3//4*(3*a*d+4*b*c)*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//4*x*c^-1*d^-1*(a+b*x^4)^3//4*(c+d*x^4)^-1*(b*c+-1*a*d))
@test integrate((a+b*x^4)^3//4*(c+d*x^4)^-2, x) == :(1//4*x*c^-1*(a+b*x^4)^3//4*(c+d*x^4)^-1+3//8*a*c^-7//4*(b*c+-1*a*d)^-1//4*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+3//8*a*c^-7//4*(b*c+-1*a*d)^-1//4*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4))
@test integrate((a+b*x^4)^-1//4*(c+d*x^4)^-2, x) == :(1//8*c^-7//4*(b*c+-1*a*d)^-5//4*(-3*a*d+4*b*c)*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+1//8*c^-7//4*(b*c+-1*a*d)^-5//4*(-3*a*d+4*b*c)*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//4*d*x*c^-1*(a+b*x^4)^3//4*(c+d*x^4)^-1*(b*c+-1*a*d)^-1)
@test integrate((a+b*x^4)^-5//4*(c+d*x^4)^-2, x) == :(-1//8*d*c^-7//4*(b*c+-1*a*d)^-9//4*(-3*a*d+8*b*c)*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//8*d*c^-7//4*(b*c+-1*a*d)^-9//4*(-3*a*d+8*b*c)*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//4*d*x*c^-1*(a+b*x^4)^-1//4*(c+d*x^4)^-1*(b*c+-1*a*d)^-1+1//4*b*x*a^-1*c^-1*(a+b*x^4)^-1//4*(b*c+-1*a*d)^-2*(a*d+4*b*c))
@test integrate((a+b*x^4)^-9//4*(c+d*x^4)^-2, x) == :(3//8*c^-7//4*d^2*(b*c+-1*a*d)^-13//4*(-1*a*d+4*b*c)*arctan(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+3//8*c^-7//4*d^2*(b*c+-1*a*d)^-13//4*(-1*a*d+4*b*c)*arctanh(x*c^-1//4*(a+b*x^4)^-1//4*(b*c+-1*a*d)^1//4)+-1//4*d*x*c^-1*(a+b*x^4)^-5//4*(c+d*x^4)^-1*(b*c+-1*a*d)^-1+1//20*b*x*a^-1*c^-1*(a+b*x^4)^-5//4*(b*c+-1*a*d)^-2*(4*b*c+5*a*d)+1//20*b*x*a^-2*c^-1*(a+b*x^4)^-1//4*(b*c+-1*a*d)^-3*(-5*a^2*d^2+16*b^2*c^2+-56*a*b*c*d))
@test integrate((1+x^4)^-1//4*(2+x^4)^-1, x) == :(1//4*2^1//4*arctan((1/2)*x*2^3//4*(1+x^4)^-1//4)+1//4*2^1//4*arctanh((1/2)*x*2^3//4*(1+x^4)^-1//4))
@test integrate((a+b*x^4)^-1//4*(a+-1*x^4*(a+-1b))^-1, x) == :((1/2)*a^-5//4*arctan(x*a^1//4*(a+b*x^4)^-1//4)+(1/2)*a^-5//4*arctanh(x*a^1//4*(a+b*x^4)^-1//4))
@test integrate((a+b*x^5)^-1//5*(c+d*x^5)^-1, x) == :(-1//5*c^-4//5*(b*c+-1*a*d)^-1//5*log(c^1//5+-1*x*(a+b*x^5)^-1//5*(b*c+-1*a*d)^1//5)+-1//5*c^-4//5*(5//2+(1/2)*5^(1/2))^(1/2)*(b*c+-1*a*d)^-1//5*arctan((1+-2//5*5^(1/2))^(1/2)+-2*x*2^(1/2)*c^-1//5*(5+5^(1/2))^-1//2*(a+b*x^5)^-1//5*(b*c+-1*a*d)^1//5)+1//5*c^-4//5*(5//2+-1//2*5^(1/2))^(1/2)*(b*c+-1*a*d)^-1//5*arctan((1+2//5*5^(1/2))^(1/2)+x*c^-1//5*(2+2//5*5^(1/2))^(1/2)*(a+b*x^5)^-1//5*(b*c+-1*a*d)^1//5)+1//20*c^-4//5*(b*c+-1*a*d)^-1//5*(1+5^(1/2))*log((a+b*x^5)^-2//5*(2*c^2//5*(a+b*x^5)^2//5+2*x^2*(b*c+-1*a*d)^2//5+x*c^1//5*(a+b*x^5)^1//5*(b*c+-1*a*d)^1//5+x*5^(1/2)*c^1//5*(a+b*x^5)^1//5*(b*c+-1*a*d)^1//5))+1//20*c^-4//5*(b*c+-1*a*d)^-1//5*(1+-1*5^(1/2))*log((a+b*x^5)^-2//5*(2*c^2//5*(a+b*x^5)^2//5+2*x^2*(b*c+-1*a*d)^2//5+x*c^1//5*(a+b*x^5)^1//5*(b*c+-1*a*d)^1//5+-1*x*5^(1/2)*c^1//5*(a+b*x^5)^1//5*(b*c+-1*a*d)^1//5)))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^3, x) == :(x*(a+b*x^-1)^(1/2)*(c+d*x^-1)^3+-7//5*d*(a+b*x^-1)^(1/2)*(c+d*x^-1)^2+a^-1//2*c^2*(b*c+6*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1//15*d*b^-2*(a+b*x^-1)^(1/2)*(-4*a^2*d^2+114*b^2*c^2+b*d*x^-1*(2*a*d+33*b*c)+30*a*b*c*d))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^2, x) == :(-2//3*b^-1*d^2*(a+b*x^-1)^3//2+c*a^-1//2*(b*c+4*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^2*(a+b*x^-1)^3//2+-1*c*a^-1*(a+b*x^-1)^(1/2)*(b*c+4*a*d))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1), x) == :(a^-1//2*(b*c+2*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1*a^-1*(a+b*x^-1)^(1/2)*(b*c+2*a*d)+c*x*a^-1*(a+b*x^-1)^3//2)
@test integrate((a+b*x^-1)^(1/2), x) == :(x*(a+b*x^-1)^(1/2)+b*a^-1//2*arctanh(a^-1//2*(a+b*x^-1)^(1/2)))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^-1, x) == :(x*c^-1*(a+b*x^-1)^(1/2)+a^-1//2*c^-2*(b*c+-2*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+2*c^-2*d^(1/2)*(b*c+-1*a*d)^(1/2)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^-2, x) == :(x*c^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1+a^-1//2*c^-3*(b*c+-4*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+2*d*c^-2*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1+c^-3*d^(1/2)*(b*c+-1*a*d)^-1//2*(-4*a*d+3*b*c)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^-3, x) == :(x*c^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-2+a^-1//2*c^-4*(b*c+-6*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+3//2*d*c^-2*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-2+1//4*c^-4*d^(1/2)*(b*c+-1*a*d)^-3//2*(15*b^2*c^2+24*a^2*d^2+-40*a*b*c*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+1//4*d*c^-3*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1*(b*c+-1*a*d)^-1*(-12*a*d+11*b*c))
@test integrate((a+b*x^-1)^3//2*(c+d*x^-1)^3, x) == :(x*(a+b*x^-1)^3//2*(c+d*x^-1)^3+-3*c^2*(a+b*x^-1)^(1/2)*(b*c+2*a*d)+-9//7*d*(a+b*x^-1)^3//2*(c+d*x^-1)^2+3*a^(1/2)*c^2*(b*c+2*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1//35*d*b^-2*(a+b*x^-1)^3//2*((-2*a*d+26*b*c)*(2*a*d+5*b*c)+3*b*d*x^-1*(2*a*d+19*b*c)))
@test integrate((a+b*x^-1)^3//2*(c+d*x^-1)^2, x) == :(-1*c*(a+b*x^-1)^(1/2)*(3*b*c+4*a*d)+-2//5*b^-1*d^2*(a+b*x^-1)^5//2+c*a^(1/2)*(3*b*c+4*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^2*(a+b*x^-1)^5//2+-1//3*c*a^-1*(a+b*x^-1)^3//2*(3*b*c+4*a*d))
@test integrate((a+b*x^-1)^3//2*(c+d*x^-1), x) == :(-1*(a+b*x^-1)^(1/2)*(2*a*d+3*b*c)+a^(1/2)*(2*a*d+3*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1//3*a^-1*(a+b*x^-1)^3//2*(2*a*d+3*b*c)+c*x*a^-1*(a+b*x^-1)^5//2)
@test integrate((a+b*x^-1)^3//2, x) == :(x*(a+b*x^-1)^3//2+-3*b*(a+b*x^-1)^(1/2)+3*b*a^(1/2)*arctanh(a^-1//2*(a+b*x^-1)^(1/2)))
@test integrate((a+b*x^-1)^3//2*(c+d*x^-1)^-1, x) == :(a*x*c^-1*(a+b*x^-1)^(1/2)+a^(1/2)*c^-2*(-2*a*d+3*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-2*c^-2*d^-1//2*(b*c+-1*a*d)^3//2*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate((a+b*x^-1)^3//2*(c+d*x^-1)^-2, x) == :(a^(1/2)*c^-3*(-4*a*d+3*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1*c^-2*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1*(b*c+-2*a*d)+a*x*c^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1+-1*c^-3*d^-1//2*(b*c+-1*a*d)^(1/2)*(b*c+-4*a*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate((a+b*x^-1)^3//2*(c+d*x^-1)^-3, x) == :(3*a^(1/2)*c^-4*(b*c+-2*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1//2*c^-2*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-2*(b*c+-3*a*d)+-1//4*c^-3*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1*(-12*a*d+3*b*c)+a*x*c^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-2+-1//4*c^-4*d^-1//2*(b*c+-1*a*d)^-1//2*(3*b^2*c^2+24*a^2*d^2+-24*a*b*c*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate((a+b*x^-1)^5//2*(c+d*x^-1)^3, x) == :(x*(a+b*x^-1)^5//2*(c+d*x^-1)^3+-11//9*d*(a+b*x^-1)^5//2*(c+d*x^-1)^2+-1//3*c^2*(a+b*x^-1)^3//2*(5*b*c+6*a*d)+a^3//2*c^2*(5*b*c+6*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1*a*c^2*(a+b*x^-1)^(1/2)*(5*b*c+6*a*d)+-1//315*d*b^-2*(a+b*x^-1)^5//2*(-20*a^2*d^2+938*b^2*c^2+5*b*d*x^-1*(10*a*d+89*b*c)+270*a*b*c*d))
@test integrate((a+b*x^-1)^5//2*(c+d*x^-1)^2, x) == :(-2//7*b^-1*d^2*(a+b*x^-1)^7//2+-1//3*c*(a+b*x^-1)^3//2*(4*a*d+5*b*c)+c*a^3//2*(4*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^2*(a+b*x^-1)^7//2+-1*a*c*(a+b*x^-1)^(1/2)*(4*a*d+5*b*c)+-1//5*c*a^-1*(a+b*x^-1)^5//2*(4*a*d+5*b*c))
@test integrate((a+b*x^-1)^5//2*(c+d*x^-1), x) == :(-1*(a+b*x^-1)^3//2*(2//3*a*d+5//3*b*c)+a^3//2*(2*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-1*a*(a+b*x^-1)^(1/2)*(2*a*d+5*b*c)+-1//5*a^-1*(a+b*x^-1)^5//2*(2*a*d+5*b*c)+c*x*a^-1*(a+b*x^-1)^7//2)
@test integrate((a+b*x^-1)^5//2, x) == :(x*(a+b*x^-1)^5//2+-5//3*b*(a+b*x^-1)^3//2+-5*a*b*(a+b*x^-1)^(1/2)+5*b*a^3//2*arctanh(a^-1//2*(a+b*x^-1)^(1/2)))
@test integrate((a+b*x^-1)^5//2*(c+d*x^-1)^-1, x) == :(a*x*c^-1*(a+b*x^-1)^3//2+a^3//2*c^-2*(-2*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+2*c^-2*d^-3//2*(b*c+-1*a*d)^5//2*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+-1*b*c^-1*d^-1*(a+b*x^-1)^(1/2)*(a*d+2*b*c))
@test integrate((a+b*x^-1)^5//2*(c+d*x^-1)^-2, x) == :(a^3//2*c^-3*(-4*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+a*x*c^-1*(a+b*x^-1)^3//2*(c+d*x^-1)^-1+-1*c^-3*d^-3//2*(b*c+-1*a*d)^3//2*(b*c+4*a*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+c^-2*d^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1*(b*c+-1*a*d)*(b*c+-2*a*d))
@test integrate((a+b*x^-1)^5//2*(c+d*x^-1)^-3, x) == :(a^3//2*c^-4*(-6*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+a*x*c^-1*(a+b*x^-1)^3//2*(c+d*x^-1)^-2+-1//4*c^-4*d^-3//2*(b*c+-1*a*d)^(1/2)*(b^2*c^2+-24*a^2*d^2+8*a*b*c*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+-1//4*c^-3*d^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1*(b^2*c^2+-12*a^2*d^2+7*a*b*c*d)+(1/2)*c^-2*d^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-2*(b*c+-1*a*d)*(b*c+-3*a*d))
@test integrate((a+b*x^-1)^-1//2*(c+d*x^-1)^3, x) == :(-1*a^-3//2*c^2*(b*c+-6*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+c*x*a^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^2+-1//3*d*a^-1*b^-2*(a+b*x^-1)^(1/2)*(-4*a^2*d^2+6*b^2*c^2+b*d*x^-1*(2*a*d+3*b*c)+18*a*b*c*d))
@test integrate((a+b*x^-1)^-1//2*(c+d*x^-1)^2, x) == :(-2*b^-1*d^2*(a+b*x^-1)^(1/2)+x*a^-1*c^2*(a+b*x^-1)^(1/2)+-1*c*a^-3//2*(b*c+-4*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2)))
@test integrate((a+b*x^-1)^-1//2*(c+d*x^-1), x) == :(-1*a^-3//2*(b*c+-2*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+c*x*a^-1*(a+b*x^-1)^(1/2))
@test integrate((a+b*x^-1)^-1//2, x) == :(x*a^-1*(a+b*x^-1)^(1/2)+-1*b*a^-3//2*arctanh(a^-1//2*(a+b*x^-1)^(1/2)))
@test integrate((a+b*x^-1)^-1//2*(c+d*x^-1)^-1, x) == :(x*a^-1*c^-1*(a+b*x^-1)^(1/2)+-1*a^-3//2*c^-2*(b*c+2*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-2*c^-2*d^3//2*(b*c+-1*a*d)^-1//2*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2))
@test integrate((a+b*x^-1)^-1//2*(c+d*x^-1)^-2, x) == :(-1*a^-3//2*c^-3*(b*c+4*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1+-1*c^-3*d^3//2*(b*c+-1*a*d)^-3//2*(-4*a*d+5*b*c)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+d*a^-1*c^-2*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1*(b*c+-1*a*d)^-1*(b*c+-2*a*d))
@test integrate((a+b*x^-1)^-1//2*(c+d*x^-1)^-3, x) == :(-1*a^-3//2*c^-4*(b*c+6*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-2+-1//4*c^-4*d^3//2*(b*c+-1*a*d)^-5//2*(24*a^2*d^2+35*b^2*c^2+-56*a*b*c*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+(1/2)*d*a^-1*c^-2*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-2*(b*c+-1*a*d)^-1*(-3*a*d+2*b*c)+1//4*d*a^-1*c^-3*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1*(b*c+-1*a*d)^-2*(b*c+-4*a*d)*(-3*a*d+4*b*c))
@test integrate((a+b*x^-1)^-3//2*(c+d*x^-1)^3, x) == :(a^-2*b^-2*(a+b*x^-1)^-1//2*((b*c+-2*a*d)*(2*a^2*d^2+3*b^2*c^2+-2*a*b*c*d)+-1*a*b*d^2*x^-1*(b*c+2*a*d))+-3*a^-5//2*c^2*(b*c+-2*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+c*x*a^-1*(a+b*x^-1)^-1//2*(c+d*x^-1)^2)
@test integrate((a+b*x^-1)^-3//2*(c+d*x^-1)^2, x) == :(x*a^-1*c^2*(a+b*x^-1)^-1//2+a^-2*b^-1*(a+b*x^-1)^-1//2*(2*a^2*d^2+b*c*(-4*a*d+3*b*c))+-1*c*a^-5//2*(-4*a*d+3*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2)))
@test integrate((a+b*x^-1)^-3//2*(c+d*x^-1), x) == :(a^-2*(a+b*x^-1)^-1//2*(-2*a*d+3*b*c)+-1*a^-5//2*(-2*a*d+3*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+c*x*a^-1*(a+b*x^-1)^-1//2)
@test integrate((a+b*x^-1)^-3//2, x) == :(x*a^-1*(a+b*x^-1)^-1//2+-3*b*a^-5//2*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+3*b*a^-2*(a+b*x^-1)^-1//2)
@test integrate((a+b*x^-1)^-3//2*(c+d*x^-1)^-1, x) == :(x*a^-1*c^-1*(a+b*x^-1)^-1//2+-1*a^-5//2*c^-2*(2*a*d+3*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+2*c^-2*d^5//2*(b*c+-1*a*d)^-3//2*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+b*a^-2*c^-1*(a+b*x^-1)^-1//2*(b*c+-1*a*d)^-1*(-1*a*d+3*b*c))
@test integrate((a+b*x^-1)^-3//2*(c+d*x^-1)^-2, x) == :(-1*a^-5//2*c^-3*(3*b*c+4*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^-1*(a+b*x^-1)^-1//2*(c+d*x^-1)^-1+c^-3*d^5//2*(b*c+-1*a*d)^-5//2*(-4*a*d+7*b*c)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+b*a^-2*c^-2*(a+b*x^-1)^-1//2*(b*c+-1*a*d)^-2*(2*a^2*d^2+3*b^2*c^2+-2*a*b*c*d)+d*a^-1*c^-2*(a+b*x^-1)^-1//2*(c+d*x^-1)^-1*(b*c+-1*a*d)^-1*(b*c+-2*a*d))
@test integrate((a+b*x^-1)^-3//2*(c+d*x^-1)^-3, x) == :(-1*a^-5//2*c^-4*(3*b*c+6*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^-1*(a+b*x^-1)^-1//2*(c+d*x^-1)^-2+3//4*c^-4*d^5//2*(b*c+-1*a*d)^-7//2*(8*a^2*d^2+21*b^2*c^2+-24*a*b*c*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+(1/2)*d*a^-1*c^-2*(a+b*x^-1)^-1//2*(c+d*x^-1)^-2*(b*c+-1*a*d)^-1*(-3*a*d+2*b*c)+1//4*d*a^-1*c^-3*(a+b*x^-1)^-1//2*(c+d*x^-1)^-1*(b*c+-1*a*d)^-2*(4*b^2*c^2+12*a^2*d^2+-21*a*b*c*d)+3//4*b*a^-2*c^-3*(a+b*x^-1)^-1//2*(b*c+-1*a*d)^-3*(-1*a*d+2*b*c)*(2*b^2*c^2+4*a^2*d^2+-1*a*b*c*d))
@test integrate((a+b*x^-1)^-5//2*(c+d*x^-1)^3, x) == :(-1*a^-7//2*c^2*(-6*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+c*x*a^-1*(a+b*x^-1)^-3//2*(c+d*x^-1)^2+1//3*a^-3*b^-2*x^-1*(a+b*x^-1)^-3//2*(b*c+-1*a*d)*(15*b^3*c^2+-4*x*a^3*d^2+a*c*b^2*(-3d+20*c*x)+-2*b*d*a^2*(3d+5*c*x)))
@test integrate((a+b*x^-1)^-5//2*(c+d*x^-1)^2, x) == :(c*a^-3*(a+b*x^-1)^-1//2*(-4*a*d+5*b*c)+x*a^-1*c^2*(a+b*x^-1)^-3//2+-1*c*a^-7//2*(-4*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+1//3*a^-2*b^-1*(a+b*x^-1)^-3//2*(2*a^2*d^2+b*c*(-4*a*d+5*b*c)))
@test integrate((a+b*x^-1)^-5//2*(c+d*x^-1), x) == :(a^-3*(a+b*x^-1)^-1//2*(-2*a*d+5*b*c)+-1*a^-7//2*(-2*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+1//3*a^-2*(a+b*x^-1)^-3//2*(-2*a*d+5*b*c)+c*x*a^-1*(a+b*x^-1)^-3//2)
@test integrate((a+b*x^-1)^-5//2, x) == :(x*a^-1*(a+b*x^-1)^-3//2+-5*b*a^-7//2*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+5*b*a^-3*(a+b*x^-1)^-1//2+5//3*b*a^-2*(a+b*x^-1)^-3//2)
@test integrate((a+b*x^-1)^-5//2*(c+d*x^-1)^-1, x) == :(x*a^-1*c^-1*(a+b*x^-1)^-3//2+-1*a^-7//2*c^-2*(2*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+-2*c^-2*d^7//2*(b*c+-1*a*d)^-5//2*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+b*a^-3*c^-1*(a+b*x^-1)^-1//2*(b*c+-1*a*d)^-2*(a^2*d^2+5*b^2*c^2+-8*a*b*c*d)+1//3*b*a^-2*c^-1*(a+b*x^-1)^-3//2*(b*c+-1*a*d)^-1*(-3*a*d+5*b*c))
@test integrate((a+b*x^-1)^-5//2*(c+d*x^-1)^-2, x) == :(-1*a^-7//2*c^-3*(4*a*d+5*b*c)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^-1*(a+b*x^-1)^-3//2*(c+d*x^-1)^-1+-1*c^-3*d^7//2*(b*c+-1*a*d)^-7//2*(-4*a*d+9*b*c)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+1//3*b*a^-2*c^-2*(a+b*x^-1)^-3//2*(b*c+-1*a*d)^-2*(5*b^2*c^2+6*a^2*d^2+-6*a*b*c*d)+b*a^-3*c^-2*(a+b*x^-1)^-1//2*(b*c+-1*a*d)^-3*(b*c+-2*a*d)*(a^2*d^2+5*b^2*c^2+-1*a*b*c*d)+d*a^-1*c^-2*(a+b*x^-1)^-3//2*(c+d*x^-1)^-1*(b*c+-1*a*d)^-1*(b*c+-2*a*d))
@test integrate((a+b*x^-1)^-5//2*(c+d*x^-1)^-3, x) == :(-1*a^-7//2*c^-4*(5*b*c+6*a*d)*arctanh(a^-1//2*(a+b*x^-1)^(1/2))+x*a^-1*c^-1*(a+b*x^-1)^-3//2*(c+d*x^-1)^-2+-1//4*c^-4*d^7//2*(b*c+-1*a*d)^-9//2*(24*a^2*d^2+99*b^2*c^2+-88*a*b*c*d)*arctan(d^(1/2)*(a+b*x^-1)^(1/2)*(b*c+-1*a*d)^-1//2)+1//4*b*a^-3*c^-3*(a+b*x^-1)^-1//2*(b*c+-1*a*d)^-4*(12*a^4*d^4+20*b^4*c^4+-56*a*d*b^3*c^3+-35*b*c*a^3*d^3+24*a^2*b^2*c^2*d^2)+1//12*b*a^-2*c^-3*(a+b*x^-1)^-3//2*(b*c+-1*a*d)^-3*(-36*a^3*d^3+20*b^3*c^3+-36*a*d*b^2*c^2+87*b*c*a^2*d^2)+(1/2)*d*a^-1*c^-2*(a+b*x^-1)^-3//2*(c+d*x^-1)^-2*(b*c+-1*a*d)^-1*(-3*a*d+2*b*c)+1//4*d*a^-1*c^-3*(a+b*x^-1)^-3//2*(c+d*x^-1)^-1*(b*c+-1*a*d)^-2*(4*b^2*c^2+12*a^2*d^2+-23*a*b*c*d))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^(1/2), x) == :(x*(a+b*x^-1)^(1/2)*(c+d*x^-1)^(1/2)+-2*b^(1/2)*d^(1/2)*arctanh(b^-1//2*d^(1/2)*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1//2)+a^-1//2*c^-1//2*(a*d+b*c)*arctanh(a^-1//2*c^(1/2)*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1//2))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^-1//2, x) == :(x*c^-1*(a+b*x^-1)^(1/2)*(c+d*x^-1)^(1/2)+a^-1//2*c^-3//2*(b*c+-1*a*d)*arctanh(a^-1//2*c^(1/2)*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1//2))
@test integrate((a+b*x^-1)^(1/2)*(c+d*x^-1)^-3//2, x) == :(a^-1//2*c^-5//2*(b*c+-3*a*d)*arctanh(a^-1//2*c^(1/2)*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1//2)+x*a^-1*c^-1*(a+b*x^-1)^3//2*(c+d*x^-1)^-1//2+-1*a^-1*c^-2*(a+b*x^-1)^(1/2)*(c+d*x^-1)^-1//2*(b*c+-3*a*d))
@test integrate((c+d*x^-2)^-1*(a+b*x^-2), x) == :(a*x*c^-1+c^-3//2*d^-1//2*(b*c+-1*a*d)*arctan(x*c^(1/2)*d^-1//2))
@test integrate((a+b*x^-2)^(1/2)*(c+d*x^-2)^(1/2), x) == :(x*(a+b*x^-2)^(1/2)*(c+d*x^-2)^(1/2)+-2*d*x^-1*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2+2*c^(1/2)*d^(1/2)*(c*a^-1*(c+d*x^-2)^-1*(a+b*x^-2))^-1//2*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2*Elliptic.E(arccot(x*c^(1/2)*d^-1//2),1+-1*b*c*a^-1*d^-1)+-1*a^-1*c^(1/2)*d^-1//2*(c*a^-1*(c+d*x^-2)^-1*(a+b*x^-2))^-1//2*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2*(a*d+b*c)*Elliptic.F(arccot(x*c^(1/2)*d^-1//2),1+-1*b*c*a^-1*d^-1))
@test integrate((a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2, x) == :(x*c^-1*(a+b*x^-2)^(1/2)*(c+d*x^-2)^(1/2)+-1*d*c^-1*x^-1*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2+c^-1//2*d^(1/2)*(c*a^-1*(c+d*x^-2)^-1*(a+b*x^-2))^-1//2*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2*Elliptic.E(arccot(x*c^(1/2)*d^-1//2),1+-1*b*c*a^-1*d^-1)+-1*b*a^-1*c^(1/2)*d^-1//2*(c*a^-1*(c+d*x^-2)^-1*(a+b*x^-2))^-1//2*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2*Elliptic.F(arccot(x*c^(1/2)*d^-1//2),1+-1*b*c*a^-1*d^-1))
@test integrate((a+b*x^-2)^(1/2)*(c+d*x^-2)^-3//2, x) == :(-1*x*c^-1*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2+2*x*c^-2*(a+b*x^-2)^(1/2)*(c+d*x^-2)^(1/2)+-2*d*c^-2*x^-1*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2+2*c^-3//2*d^(1/2)*(c*a^-1*(c+d*x^-2)^-1*(a+b*x^-2))^-1//2*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2*Elliptic.E(arccot(x*c^(1/2)*d^-1//2),1+-1*b*c*a^-1*d^-1)+-1*b*a^-1*c^-1//2*d^-1//2*(c*a^-1*(c+d*x^-2)^-1*(a+b*x^-2))^-1//2*(a+b*x^-2)^(1/2)*(c+d*x^-2)^-1//2*Elliptic.F(arccot(x*c^(1/2)*d^-1//2),1+-1*b*c*a^-1*d^-1))
@test integrate((c+d*x^-3)^-1*(a+b*x^-3), x) == :(a*x*c^-1+-1//6*c^-4//3*d^-2//3*(b*c+-1*a*d)*log(d^2//3+c^2//3*x^2+-1*x*c^1//3*d^1//3)+1//3*c^-4//3*d^-2//3*(b*c+-1*a*d)*log(d^1//3+x*c^1//3)+-1//3*3^(1/2)*c^-4//3*d^-2//3*(b*c+-1*a*d)*arctan(1//3*3^(1/2)*d^-1//3*(d^1//3+-2*x*c^1//3)))
@test integrate((c+d*x^(1/2))^-1*(a+b*x^(1/2)), x) == :(b*x*d^-1+-1*d^-2*x^(1/2)*(-2*a*d+2*b*c)+2*c*d^-3*(b*c+-1*a*d)*log(c+d*x^(1/2)))
@test integrate((1+x^1//3)^-1*(-1+x^1//3), x) == :(x+-6*log(1+x^1//3)+-3*x^2//3+6*x^1//3)
@test integrate((-1+x^2//3)^-1*(1+x^2//3), x) == :(x+-6*arctanh(x^1//3)+6*x^1//3)
@test integrate((16+x^3//4)^-1*(-16+x^3//4), x) == :(x+-128*x^1//4+-128//3*2^1//3*log(x^(1/2)+4*2^2//3+-2*2^1//3*x^1//4)+256//3*2^1//3*log(x^1//4+2*2^1//3)+-256//3*2^1//3*3^(1/2)*arctan(1//6*2^2//3*3^(1/2)*(2^1//3+-1*x^1//4)))
@test integrate((-1+x^-1//3)^-1*(1+x^-1//3), x) == :(-1x+-6*x^1//3+-6*log(1+-1*x^1//3)+-3*x^2//3)
@test integrate((c+d*x^n)^4*(a+b*x^n), x) == :(a*x*c^4+b*d^4*x^(1+5n)*(1+5n)^-1+c^3*x^(1+n)*(1+n)^-1*(b*c+4*a*d)+d^3*x^(1+4n)*(1+4n)^-1*(a*d+4*b*c)+2*c*d^2*x^(1+3n)*(1+3n)^-1*(2*a*d+3*b*c)+2*d*c^2*x^(1+2n)*(1+2n)^-1*(2*b*c+3*a*d))
@test integrate((c+d*x^n)^3*(a+b*x^n), x) == :(a*x*c^3+b*d^3*x^(1+4n)*(1+4n)^-1+c^2*x^(1+n)*(1+n)^-1*(b*c+3*a*d)+d^2*x^(1+3n)*(1+3n)^-1*(a*d+3*b*c)+3*c*d*x^(1+2n)*(1+2n)^-1*(a*d+b*c))
@test integrate((c+d*x^n)^2*(a+b*x^n), x) == :(a*x*c^2+b*d^2*x^(1+3n)*(1+3n)^-1+c*x^(1+n)*(1+n)^-1*(b*c+2*a*d)+d*x^(1+2n)*(1+2n)^-1*(a*d+2*b*c))
@test integrate((a+b*x^n)*(c+d*x^n), x) == :(a*c*x+x^(1+n)*(1+n)^-1*(a*d+b*c)+b*d*x^(1+2n)*(1+2n)^-1)
@test integrate((a+b*x^n)^2*(d+e*x^n)^3, x) == :(x*a^2*d^3+d*x^(1+2n)*(1+2n)^-1*(b^2*d^2+3*a^2*e^2+6*a*b*d*e)+e*x^(1+3n)*(1+3n)^-1*(a^2*e^2+3*b^2*d^2+6*a*b*d*e)+b^2*e^3*x^(1+5n)*(1+5n)^-1+a*d^2*x^(1+n)*(1+n)^-1*(2*b*d+3*a*e)+b*e^2*x^(1+4n)*(1+4n)^-1*(2*a*e+3*b*d))
@test integrate((a+b*x^n)^2*(d+e*x^n)^2, x) == :(x*a^2*d^2+x^(1+2n)*(1+2n)^-1*(a^2*e^2+b^2*d^2+4*a*b*d*e)+b^2*e^2*x^(1+4n)*(1+4n)^-1+2*a*d*x^(1+n)*(1+n)^-1*(a*e+b*d)+2*b*e*x^(1+3n)*(1+3n)^-1*(a*e+b*d))
@test integrate((a+b*x^n)^2*(c+d*x^n), x) == :(c*x*a^2+a*x^(1+n)*(1+n)^-1*(a*d+2*b*c)+b*x^(1+2n)*(1+2n)^-1*(b*c+2*a*d)+d*b^2*x^(1+3n)*(1+3n)^-1)
@test integrate((c+d*x^n)^(-2+-1*n^-1)*(a+b*x^n), x) == :(x*c^-1*(1+n)^-1*(c+d*x^n)^(-1+-1*n^-1)*(a+b*x^n)+a*n*x*c^-2*(1+n)^-1*(c+d*x^n)^(-1*n^-1))
@test integrate((c+d*x^n)^(-1+-1*n^-1), x) == :(x*c^-1*(c+d*x^n)^(-1*n^-1))
@test integrate((a+b*x^n)^(n^-1*(b*c+-1*a*d)^-1*(a*d*n+-1*b*c*(1+n)))*(c+d*x^n)^((b*c*n+-1*a*d*n)^-1*(a*d+a*d*n+-1*b*c*n)), x) == :(x*a^-1*c^-1*(a+b*x^n)^(-1*b*c*n^-1*(b*c+-1*a*d)^-1)*(c+d*x^n)^(a*d*n^-1*(b*c+-1*a*d)^-1))
@test integrate((c+d*x^n)^(-3+-1*n^-1)*(a+b*x^n), x) == :(-1*x*c^-1*d^-1*(1+2n)^-1*(c+d*x^n)^(-2+-1*n^-1)*(b*c+-1*a*d)+x*c^-2*d^-1*(1+n)^-1*(1+2n)^-1*(c+d*x^n)^(-1+-1*n^-1)*(b*c+2*a*d*n)+n*x*c^-3*d^-1*(1+n)^-1*(1+2n)^-1*(c+d*x^n)^(-1*n^-1)*(b*c+2*a*d*n))
@test integrate((c+d*x^n)^(-2+-1*n^-1), x) == :(x*c^-1*(1+n)^-1*(c+d*x^n)^(-1+-1*n^-1)+n*x*c^-2*(1+n)^-1*(c+d*x^n)^(-1*n^-1))
@test integrate(x^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(a+b*x^2), x) == :(a*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)+-1*a*c*arctan(c^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))+1//3*b*d^-2*(c+d*x)^3//2*(-1c+d*x)^3//2)
@test integrate(x^4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(a+b*x^2), x) == :(-1//64*c^6*d^-7*(5*b*c^2+8*a*d^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+1//8*b*d^-2*x^5*(c+d*x)^3//2*(-1c+d*x)^3//2+1//48*d^-4*x^3*(c+d*x)^3//2*(-1c+d*x)^3//2*(5*b*c^2+8*a*d^2)+1//64*x*c^2*d^-6*(c+d*x)^3//2*(-1c+d*x)^3//2*(5*b*c^2+8*a*d^2)+1//128*x*c^4*d^-6*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(5*b*c^2+8*a*d^2))
@test integrate(x^2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(a+b*x^2), x) == :(-1//8*c^4*d^-5*(b*c^2+2*a*d^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+1//6*b*d^-2*x^3*(c+d*x)^3//2*(-1c+d*x)^3//2+1//8*x*d^-4*(c+d*x)^3//2*(-1c+d*x)^3//2*(b*c^2+2*a*d^2)+1//16*x*c^2*d^-4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(b*c^2+2*a*d^2))
@test integrate((c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(a+b*x^2), x) == :(-1//4*c^2*d^-3*(b*c^2+4*a*d^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+1//4*b*x*d^-2*(c+d*x)^3//2*(-1c+d*x)^3//2+1//8*x*d^-2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(b*c^2+4*a*d^2))
@test integrate(x^-2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(a+b*x^2), x) == :(-1*d^-1*(b*c^2+-2*a*d^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+x*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*((1/2)*b+-1*a*c^-2*d^2)+a*c^-2*x^-1*(c+d*x)^3//2*(-1c+d*x)^3//2)
@test integrate(x^-4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(a+b*x^2), x) == :(2*b*d*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+-1*b*x^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)+1//3*a*c^-2*x^-3*(c+d*x)^3//2*(-1c+d*x)^3//2)
@test integrate(x^4*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(1//16*c^-7*(5b+6*a*c^2)*arccosh(c*x)+1//6*b*c^-2*x^5*(1+c*x)^(1/2)*(-1+c*x)^(1/2)+1//16*x*c^-6*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(5b+6*a*c^2)+1//24*c^-4*x^3*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(5b+6*a*c^2))
@test integrate(x^3*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(1//15*c^-6*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(8b+10*a*c^2)+1//5*b*c^-2*x^4*(1+c*x)^(1/2)*(-1+c*x)^(1/2)+1//15*c^-4*x^2*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(4b+5*a*c^2))
@test integrate(x^2*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(1//8*c^-5*(3b+4*a*c^2)*arccosh(c*x)+1//4*b*c^-2*x^3*(1+c*x)^(1/2)*(-1+c*x)^(1/2)+1//8*x*c^-4*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(3b+4*a*c^2))
@test integrate(x*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(1//3*c^-4*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(2b+3*a*c^2)+1//3*b*c^-2*x^2*(1+c*x)^(1/2)*(-1+c*x)^(1/2))
@test integrate((1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :((1/2)*c^-3*(b+2*a*c^2)*arccosh(c*x)+(1/2)*b*x*c^-2*(1+c*x)^(1/2)*(-1+c*x)^(1/2))
@test integrate(x^-1*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(a*arctan((1+c*x)^(1/2)*(-1+c*x)^(1/2))+b*c^-2*(1+c*x)^(1/2)*(-1+c*x)^(1/2))
@test integrate(x^-2*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(b*c^-1*arccosh(c*x)+a*x^-1*(1+c*x)^(1/2)*(-1+c*x)^(1/2))
@test integrate(x^-3*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :((b+(1/2)*a*c^2)*arctan((1+c*x)^(1/2)*(-1+c*x)^(1/2))+(1/2)*a*x^-2*(1+c*x)^(1/2)*(-1+c*x)^(1/2))
@test integrate(x^-4*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(1//3*a*x^-3*(1+c*x)^(1/2)*(-1+c*x)^(1/2)+1//3*x^-1*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(3b+2*a*c^2))
@test integrate(x^-5*(1+c*x)^-1//2*(-1+c*x)^-1//2*(a+b*x^2), x) == :(1//8*c^2*(4b+3*a*c^2)*arctan((1+c*x)^(1/2)*(-1+c*x)^(1/2))+1//4*a*x^-4*(1+c*x)^(1/2)*(-1+c*x)^(1/2)+1//8*x^-2*(1+c*x)^(1/2)*(-1+c*x)^(1/2)*(4b+3*a*c^2))
@test integrate(x^4*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(1//8*c^4*d^-7*(5*b*c^2+6*a*d^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+1//6*b*d^-2*x^5*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)+1//24*d^-4*x^3*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(5*b*c^2+6*a*d^2)+1//16*x*c^2*d^-6*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(5*b*c^2+6*a*d^2))
@test integrate(x^3*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(1//5*b*d^-2*x^4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)+1//15*d^-4*x^2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(4*b*c^2+5*a*d^2)+2//15*c^2*d^-6*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(4*b*c^2+5*a*d^2))
@test integrate(x^2*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(1//4*c^2*d^-5*(3*b*c^2+4*a*d^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+1//4*b*d^-2*x^3*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)+1//8*x*d^-4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(3*b*c^2+4*a*d^2))
@test integrate(x*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(1//3*d^-4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(2*b*c^2+3*a*d^2)+1//3*b*d^-2*x^2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))
@test integrate((c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(d^-3*(b*c^2+2*a*d^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+(1/2)*b*x*d^-2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))
@test integrate(x^-1*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(a*c^-1*arctan(c^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))+b*d^-2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))
@test integrate(x^-2*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(2*b*d^-1*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+a*c^-2*x^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))
@test integrate(x^-3*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :((1/2)*c^-3*(a*d^2+2*b*c^2)*arctan(c^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))+(1/2)*a*c^-2*x^-2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))
@test integrate(x^-4*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(1//3*a*c^-2*x^-3*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)+1//3*c^-4*x^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(2*a*d^2+3*b*c^2))
@test integrate(x^-5*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a+b*x^2), x) == :(1//8*c^-5*d^2*(3*a*d^2+4*b*c^2)*arctan(c^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))+1//4*a*c^-2*x^-4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)+1//8*c^-4*x^-2*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(3*a*d^2+4*b*c^2))
@test integrate(x^4*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(3//4*c^2*d^-7*(4*a*d^2+5*b*c^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+-1//4*d^-4*x^3*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(4*a*d^2+5*b*c^2)+1//4*b*d^-2*x^5*(c+d*x)^-1//2*(-1c+d*x)^-1//2+1//8*x*d^-6*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(12*a*d^2+15*b*c^2))
@test integrate(x^3*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(1//3*d^-6*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(6*a*d^2+8*b*c^2)+-1//3*d^-4*x^2*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(3*a*d^2+4*b*c^2)+1//3*b*d^-2*x^4*(c+d*x)^-1//2*(-1c+d*x)^-1//2)
@test integrate(x^2*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(d^-5*(2*a*d^2+3*b*c^2)*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+-1//2*d^-5*(c+d*x)^-1//2*(-1c+d*x)^(1/2)*(2*a*d^2+3*b*c^2)+(1/2)*b*d^-2*x^3*(c+d*x)^-1//2*(-1c+d*x)^-1//2+-1//2*c*d^-5*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(2*a*d^2+3*b*c^2))
@test integrate(x*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(-1*x^2*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a*c^-2+b*d^-2)+c^-2*d^-4*(c+d*x)^(1/2)*(-1c+d*x)^(1/2)*(a*d^2+2*b*c^2))
@test integrate((c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(2*b*d^-3*arctanh((c+d*x)^-1//2*(-1c+d*x)^(1/2))+-1*x*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a*c^-2+b*d^-2))
@test integrate(x^-1*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(-1*a*c^-3*arctan(c^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))+-1*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(a*c^-2+b*d^-2))
@test integrate(x^-2*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(a*c^-2*x^-1*(c+d*x)^-1//2*(-1c+d*x)^-1//2+-1*x*c^-4*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(b*c^2+2*a*d^2))
@test integrate(x^-3*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(-1//2*c^-5*(2*b*c^2+3*a*d^2)*arctan(c^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))+-1//2*c^-4*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(2*b*c^2+3*a*d^2)+(1/2)*a*c^-2*x^-2*(c+d*x)^-1//2*(-1c+d*x)^-1//2)
@test integrate(x^-4*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(1//3*a*c^-2*x^-3*(c+d*x)^-1//2*(-1c+d*x)^-1//2+1//3*c^-4*x^-1*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(3*b*c^2+4*a*d^2)+-2//3*x*c^-6*d^2*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(3*b*c^2+4*a*d^2))
@test integrate(x^-5*(c+d*x)^-3//2*(-1c+d*x)^-3//2*(a+b*x^2), x) == :(-3//8*c^-7*d^2*(4*b*c^2+5*a*d^2)*arctan(c^-1*(c+d*x)^(1/2)*(-1c+d*x)^(1/2))+-3//8*c^-6*d^2*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(4*b*c^2+5*a*d^2)+1//4*a*c^-2*x^-4*(c+d*x)^-1//2*(-1c+d*x)^-1//2+1//8*c^-4*x^-2*(c+d*x)^-1//2*(-1c+d*x)^-1//2*(4*b*c^2+5*a*d^2))
@test integrate(x^-1*(1+c*x)^-1//2*(-1+c*x)^-1//2*(1+c^2*x^2), x) == :((1+c*x)^(1/2)*(-1+c*x)^(1/2)+arctan((1+c*x)^(1/2)*(-1+c*x)^(1/2)))
@test integrate(x^(-1*(c*b^2+d*a^2)^-1*(d*a^2+2*c*b^2))*(a+b*x)^-1//2*(-1a+b*x)^-1//2*(c+d*x^2), x) == :(x^(-1*c*b^2*(c*b^2+d*a^2)^-1)*(a+b*x)^(1/2)*(-1a+b*x)^(1/2)*(c*a^-2+d*b^-2))
@test integrate((1+x)^-1//2*(-1+x^(1/2))^-1//2*(-1+-1*x^(1/2))^-1//2, x) == :((1+-1x)^(1/2)*(-1+x^(1/2))^-1//2*(-1+-1*x^(1/2))^-1//2*arcsin(x))
@test integrate((a+b*x^(1/2))^-1//2*(a+-1*b*x^(1/2))^-1//2*(a^2+x*b^2)^-1//2, x) == :(-2*b^-2*(a+b*x^(1/2))^-1//2*(a+-1*b*x^(1/2))^-1//2*(a^2+-1*x*b^2)^(1/2)*arctan((a^2+x*b^2)^-1//2*(a^2+-1*x*b^2)^(1/2)))

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.3 (a+b x^n)^p (c+d x^n)^q.jl


text:
function set_boundary(Qbase, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, bdcon)
    """bdcon[i][j]
        i : 境界番号(x-,x+ ,y-,y+)
        j=1-6 : "bd1_con":"2",
                "bd1_rho":"1.0",
                "bd1_u":"300.0",
                "bd1_v":"0.0",
                "bd1_p":"1.0",
                "bd1_T":"300.0",
    """

    #println(bdcon)

    # ---------------------------------
    # bd1 = x-
    # ---------------------------------
    if Int(bdcon[1][1]) == 0
        for j in 1:cellymax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[1,j,k,l] = bdcon[1][l+1]
                end
            end
        end
    elseif Int(bdcon[1][1]) == 1
        for j in 1:cellymax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[1,j,k,l] = Qbase[2,j,k,l]
                end
            end
        end
    elseif Int(bdcon[1][1]) == 2
        println("\n We are preparing now! \n")
        throw(UndefVarError(:x))
        for j in 1:cellymax
            for k in 1:cellzmax
                xvec = vecAx[2,j,k,1]
                yvec = vecAx[2,j,k,2]
                zvec = vecAx[2,j,k,3]
                u = Qbase[2,j,k,2]
                v = Qbase[2,j,k,3]
                w = Qbase[2,j,k,4]

                Qbase[1,j,k,1] = Qbase[2,j,k,1]
                Qbase[1,j,k,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
                Qbase[1,j,k,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
                Qbase[1,j,k,4] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
                Qbase[1,j,k,5] = Qbase[2,j,k,5]
            end
        end
    elseif Int(bdcon[1][1]) == 3
        for j in 1:cellymax
            for k in 1:cellzmax
                u = Qbase[2,j,k,2]
                v = Qbase[2,j,k,3]
                w = Qbase[2,j,k,4]

                Qbase[1,j,k,1] = Qbase[2,j,k,1]
                Qbase[1,j,k,2] = -u
                Qbase[1,j,k,3] = -v
                Qbase[1,j,k,4] = -w
                Qbase[1,j,k,5] = Qbase[2,j,k,5]
            end
        end
    elseif Int(bdcon[1][1]) == 4
        for j in 1:cellymax
            for k in 1:cellzmax
                Qbase[1,j,k,1] = Qbase[cellxmax-1,j,k,1]
                Qbase[1,j,k,2] = Qbase[cellxmax-1,j,k,2]
                Qbase[1,j,k,3] = Qbase[cellxmax-1,j,k,3]
                Qbase[1,j,k,4] = Qbase[cellxmax-1,j,k,4]
                Qbase[1,j,k,5] = Qbase[cellxmax-1,j,k,5]
            end
        end
    elseif Int(bdcon[1][1]) == 4 && Int(bdcon[2][1]) != 4
        println("\n check boundary condition error 4! \n")
        throw(UndefVarError(:x))
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # ---------------------------------
    # bd2 = x+
    # ---------------------------------
    if Int(bdcon[2][1]) == 0
        for j in 1:cellymax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[cellxmax,j,k,l] = bdcon[2][l+1]
                end
            end
        end
    elseif Int(bdcon[2][1]) == 1
        for j in 1:cellymax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[cellxmax,j,k,l] = Qbase[cellxmax-1,j,k,l]
                end
            end
        end
    elseif Int(bdcon[2][1]) == 2
        println("\n We are preparing now! \n")
        throw(UndefVarError(:x))
        for j in 1:cellymax
            xvec = vecAx[cellymax,j,k,1]
            yvec = vecAx[cellymax,j,k,2]
            u = Qbase[cellxmax-1,j,k,2]
            v = Qbase[cellxmax-1,j,k,3]

            Qbase[cellxmax,j,k,1] = Qbase[cellxmax-1,j,k,1]
            Qbase[cellxmax,j,k,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[cellxmax,j,k,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
            Qbase[cellxmax,j,k,4] = Qbase[cellxmax-1,j,k,4]
        end
    elseif Int(bdcon[2][1]) == 3
        for j in 1:cellymax
            for k in 1:cellzmax
                u = Qbase[cellxmax-1,j,k,2]
                v = Qbase[cellxmax-1,j,k,3]
                w = Qbase[cellxmax-1,j,k,4]

                Qbase[cellxmax,j,k,1] = Qbase[cellxmax-1,j,k,1]
                Qbase[cellxmax,j,k,2] = -u
                Qbase[cellxmax,j,k,3] = -v
                Qbase[cellxmax,j,k,4] = -w
                Qbase[cellxmax,j,k,5] = Qbase[cellxmax-1,j,k,5]
            end
        end
    elseif Int(bdcon[2][1]) == 4
        for j in 1:cellymax
            for k in 1:cellzmax
                Qbase[cellxmax,j,k,1] = Qbase[2,j,k,1]
                Qbase[cellxmax,j,k,2] = Qbase[2,j,k,2]
                Qbase[cellxmax,j,k,3] = Qbase[2,j,k,3]
                Qbase[cellxmax,j,k,4] = Qbase[2,j,k,4]
                Qbase[cellxmax,j,k,5] = Qbase[2,j,k,5]
            end
        end
    elseif Int(bdcon[1][1]) != 4 && Int(bdcon[2][1]) == 4
        println("\n check boundary condition error 4! \n")
        throw(UndefVarError(:x))
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # ---------------------------------
    # bd3 = y-
    # ---------------------------------
    if Int(bdcon[3][1]) == 0
        for i in 1:cellxmax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[i,1,k,l] = bdcon[3][l+1]
                end
            end
        end
    elseif Int(bdcon[3][1]) == 1
        for i in 1:cellxmax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[i,1,k,l] = Qbase[i,2,k,l]
                end
            end
        end
    elseif Int(bdcon[3][1]) == 2
        println("\n We are preparing now! \n")
        throw(UndefVarError(:x))
        for i in 1:cellxmax
            for k in 1:cellzmax
                xvec = vecAy[i,2,1]
                yvec = vecAy[i,2,2]
                u = Qbase[i,2,2]
                v = Qbase[i,2,3]

                Qbase[i,1,k,1] = Qbase[i,2,1]
                Qbase[i,1,k,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
                Qbase[i,1,k,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
                Qbase[i,1,k,4] = Qbase[i,2,4]
            end
        end
    elseif Int(bdcon[3][1]) == 3
        for i in 1:cellxmax
            for k in 1:cellzmax
                u = Qbase[i,2,k,2]
                v = Qbase[i,2,k,3]
                w = Qbase[i,2,k,4]

                Qbase[i,1,k,1] = Qbase[i,2,k,1]
                Qbase[i,1,k,2] = -u
                Qbase[i,1,k,3] = -v
                Qbase[i,1,k,4] = -w
                Qbase[i,1,k,5] = Qbase[i,2,k,5]
            end
        end
    elseif Int(bdcon[3][1]) == 4
        for i in 1:cellxmax
            for k in 1:cellzmax
                Qbase[i,1,k,1] = Qbase[i,cellymax-1,k,1]
                Qbase[i,1,k,2] = Qbase[i,cellymax-1,k,2]
                Qbase[i,1,k,3] = Qbase[i,cellymax-1,k,3]
                Qbase[i,1,k,4] = Qbase[i,cellymax-1,k,4]
            end
        end
    elseif Int(bdcon[3][1]) == 4 && Int(bdcon[4][1]) != 4
        println("\n check boundary condition error 4! \n")
        throw(UndefVarError(:x))
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # ---------------------------------
    # bd4 = y+
    # ---------------------------------
    if Int(bdcon[4][1]) == 0
        for i in 1:cellxmax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[i,cellymax,k,l] = bdcon[4][l+1]
                end
            end
        end
    elseif Int(bdcon[4][1]) == 1
        for i in 1:cellxmax
            for k in 1:cellzmax
                for l in 1:5
                    Qbase[i,cellymax,k,l] = Qbase[i,cellxmax-1,k,l]
                end
            end
        end
    elseif Int(bdcon[4][1]) == 2
        println("\n We are preparing now! \n")
        throw(UndefVarError(:x))
        for i in 1:cellxmax
            xvec = vecAy[i,cellymax,1]
            yvec = vecAy[i,cellymax,2]
            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,1] = Qbase[i,cellymax-1,1]
            Qbase[i,cellymax,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,4] = Qbase[i,cellymax-1,4]
        end
    elseif Int(bdcon[4][1]) == 3
        for i in 1:cellxmax
            for k in 1:cellzmax
                u = Qbase[i,cellymax-1,k,2]
                v = Qbase[i,cellymax-1,k,3]
                w = Qbase[i,cellymax-1,k,4]

                Qbase[i,cellymax,k,1] = Qbase[i,cellymax-1,k,1]
                Qbase[i,cellymax,k,2] = -u
                Qbase[i,cellymax,k,3] = -v
                Qbase[i,cellymax,k,4] = -w
                Qbase[i,cellymax,k,5] = Qbase[i,cellymax-1,k,5]
            end
        end
    elseif Int(bdcon[4][1]) == 4
        for i in 1:cellxmax
            for k in 1:cellzmax
                Qbase[i,cellymax,k,1] = Qbase[i,2,k,1]
                Qbase[i,cellymax,k,2] = Qbase[i,2,k,2]
                Qbase[i,cellymax,k,3] = Qbase[i,2,k,3]
                Qbase[i,cellymax,k,4] = Qbase[i,2,k,4]
                Qbase[i,cellymax,k,5] = Qbase[i,2,k,5]
            end
        end
    elseif Int(bdcon[3][1]) != 4 && Int(bdcon[4][1]) == 4
        println("\n check boundary condition error 4! \n")
        throw(UndefVarError(:x))
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # ---------------------------------
    # bd5 = z-
    # ---------------------------------
    if Int(bdcon[5][1]) == 0
        for i in 1:cellxmax
            for j in 1:cellymax
                for l in 1:5
                    Qbase[i,j,1,l] = bdcon[3][l+1]
                end
            end
        end
    elseif Int(bdcon[5][1]) == 1
        for i in 1:cellxmax
            for j in 1:cellymax
                for l in 1:5
                    Qbase[i,j,1,l] = Qbase[i,j,2,l]
                end
            end
        end
    elseif Int(bdcon[5][1]) == 2
        println("\n We are preparing now! \n")
        throw(UndefVarError(:x))
        for i in 1:cellxmax
            for k in 1:cellzmax

            end
        end
    elseif Int(bdcon[5][1]) == 3
        for i in 1:cellxmax
            for j in 1:cellymax
                u = Qbase[i,j,2,2]
                v = Qbase[i,j,2,3]
                w = Qbase[i,j,2,4]

                Qbase[i,j,1,1] = Qbase[i,j,2,1]
                Qbase[i,j,1,2] = -u
                Qbase[i,j,1,3] = -v
                Qbase[i,j,1,4] = -w
                Qbase[i,j,1,5] = Qbase[i,j,2,5]
            end
        end
    elseif Int(bdcon[5][1]) == 4
        for i in 1:cellxmax
            for j in 1:cellzmax
                Qbase[i,j,1,1] = Qbase[i,j,cellzmax-1,1]
                Qbase[i,j,1,2] = Qbase[i,j,cellzmax-1,2]
                Qbase[i,j,1,3] = Qbase[i,j,cellzmax-1,3]
                Qbase[i,j,1,4] = Qbase[i,j,cellzmax-1,4]
            end
        end
    elseif Int(bdcon[5][1]) == 4 && Int(bdcon[6][1]) != 4
        println("\n check boundary condition error 4! \n")
        throw(UndefVarError(:x))
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    # ---------------------------------
    # bd6 = z+
    # ---------------------------------
    if Int(bdcon[6][1]) == 0
        for i in 1:cellxmax
            for j in 1:cellymax
                for l in 1:5
                    Qbase[i,j,cellzmax,l] = bdcon[4][l+1]
                end
            end
        end
    elseif Int(bdcon[6][1]) == 1
        for i in 1:cellxmax
            for j in 1:cellymax
                for l in 1:5
                    Qbase[i,j,cellzmax,l] = Qbase[i,cellxmax-1,k,l]
                end
            end
        end
    elseif Int(bdcon[6][1]) == 2
        println("\n We are preparing now! \n")
        throw(UndefVarError(:x))
        for i in 1:cellxmax
            xvec = vecAy[i,cellymax,1]
            yvec = vecAy[i,cellymax,2]
            u = Qbase[i,cellymax-1,2]
            v = Qbase[i,cellymax-1,3]

            Qbase[i,cellymax,1] = Qbase[i,cellymax-1,1]
            Qbase[i,cellymax,2] = ((-xvec^2+yvec^2)*u-2*xvec*yvec*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,3] = (-2*xvec*yvec*u+(xvec^2-yvec^2)*v)/(xvec^2+yvec^2)
            Qbase[i,cellymax,4] = Qbase[i,cellymax-1,4]
        end
    elseif Int(bdcon[6][1]) == 3
        for i in 1:cellxmax
            for j in 1:cellymax
                u = Qbase[i,j,cellzmax-1,2]
                v = Qbase[i,j,cellzmax-1,3]
                w = Qbase[i,j,cellzmax-1,4]

                Qbase[i,j,cellzmax,1] = Qbase[i,j,cellzmax-1,1]
                Qbase[i,j,cellzmax,2] = -u
                Qbase[i,j,cellzmax,3] = -v
                Qbase[i,j,cellzmax,4] = -w
                Qbase[i,j,cellzmax,5] = Qbase[i,j,cellzmax-1,5]
            end
        end
    elseif Int(bdcon[6][1]) == 4
        for i in 1:cellxmax
            for j in 1:cellymax
                Qbase[i,j,cellzmax,1] = Qbase[i,j,2,1]
                Qbase[i,j,cellzmax,2] = Qbase[i,j,2,2]
                Qbase[i,j,cellzmax,3] = Qbase[i,j,2,3]
                Qbase[i,j,cellzmax,4] = Qbase[i,j,2,4]
                Qbase[i,j,cellzmax,5] = Qbase[i,j,2,5]
            end
        end
    elseif Int(bdcon[5][1]) != 4 && Int(bdcon[6][1]) == 4
        println("\n check boundary condition error 4! \n")
        throw(UndefVarError(:x))
    else
        println("\n check boundary condition ! \n")
        throw(UndefVarError(:x))
    end

    return Qbase
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
using Unicode

global all_paths = []

function read_map(input_file)
    # return a dictionary with nodes as keys, and list of connected nodes as values
    map_dict = Dict("start"=>[],
                    "end"=>[])
    for line in readlines(input_file)
        first,second = split(line,"-")
        if (!).(in(first, keys(map_dict)))
            map_dict[first] = []
        end
        push!(map_dict[first], second)
        # opposite way around
        if (!).(in(second, keys(map_dict)))
            map_dict[second] = []
        end
        push!(map_dict[second], first)
    end
    # by hand don't allow to go back from 'end' to anywhere else
    map_dict["end"] = []
    return map_dict
end

function is_cave_small(cave_name)
    return all(x->Unicode.islowercase(x),cave_name)
end

function make_path_branches(current_path, map_dict)
    current_cave = current_path[length(current_path)]
    if current_cave == "end"
        push!(all_paths,current_path)
        return
    end
    for possibility in map_dict[current_cave]
        path = deepcopy(current_path)
        if (is_cave_small(possibility) == false) || (!).(in(possibility, path))
            push!(path, possibility)
            make_path_branches(path, map_dict)
        end
    end
    return
end

test_map = read_map("input.txt")

make_path_branches(["start"], test_map)
println("length of all paths is ",length(all_paths))


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_reliability(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for reporting dual variable of maximum non-served energy constraint (shadow price of reliability constraint) for each model zone and time step.
"""
function write_reliability(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	# reliability: Dual variable of maximum NSE constraint = shadow value of reliability constraint
	dfReliability = DataFrame(Zone = 1:Z)
	# Dividing dual variable for each hour with corresponding hourly weight to retrieve marginal cost of generation
	if setup["ParameterScale"] == 1
		dfReliability = hcat(dfReliability, DataFrame(transpose(dual.(EP[:cMaxNSE])./inputs["omega"]*ModelScalingFactor), :auto))
	else
		dfReliability = hcat(dfReliability, DataFrame(transpose(dual.(EP[:cMaxNSE])./inputs["omega"]), :auto))
	end

	auxNew_Names=[Symbol("Zone");[Symbol("t$t") for t in 1:T]]
	rename!(dfReliability,auxNew_Names)

	CSV.write(joinpath(path, "reliabilty.csv"), dftranspose(dfReliability, false), writeheader=false)

end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using GNSSTools
using PyPlot
pygui(true)
PyPlot.matplotlib[:rc]("text", usetex=false) # allow tex rendering

σ_n = 1
SNR = 9  # dB
v = Array(0:0.01:10)
pfa = false_alarm_pdf.(v, σ_n)
pd = detection_pdf.(v, σ_n; SNR=SNR)
V_T = 1.9
V_T_idx = findmin(abs.(v .- V_T))[2]
pfa_color = "grey"
pd_color = "lightgrey"

ymax = max(maximum(pfa), maximum(pd))
fig = figure(figsize=(6.5, 1.5))
ax = fig.add_subplot(1,1,1)
ax.plot(v, pfa, pfa_color, label=L"p_{fa}(v)")
ax.plot(v, pd, pd_color, label=L"p_{d}(v)")
ax.vlines(x=1.9, ymin=0, ymax=ymax, color="k", linestyle="--")#, label=L"V_T")
fill_between(v[V_T_idx:end], pd[V_T_idx:end], color=pd_color)
fill_between(v[V_T_idx:end], pfa[V_T_idx:end], color=pfa_color)
ylim(0, ymax+0.05)
xlabel("v")
ylabel("pᵥ(v)")
legend()
subplots_adjust(bottom=0.3, left=0.09, right=0.91)
savefig("figures/ch3_probability_detection_false_alarm.svg", dpi=300)


text:
module example_BER_awgn

# ----------------------------------------------------
# --- Modules
# ----------------------------------------------------
using DigitalComm
using PGFPlotsX

function doPlot(snrVect,ber,qamVect)
	a = 0;
	@pgf a = Axis({
				   ymode	  = "log",
				   height      ="3in",
				   width       ="4in",
				   grid,
				   xlabel      = "SNR [dB]",
				   ylabel      = "Bit Error Rate ",
				   ymax 	   = 1,
				   ymin 	   = 10.0^(-5),
				   title       = "AWGN BER for QAM",
				   legend_style="{at={(0,0)},anchor=south west,legend cell align=left,align=left,draw=white!15!black}"
				   },
				  Plot({color="red",mark="square*"},Table([snrVect,ber[1,:]])),
				  LegendEntry("QPSK"),
				  Plot({color="green",mark="*"},Table([snrVect,ber[2,:]])),
				  LegendEntry("16-QAM"),

				  Plot({color="purple",mark="triangle*"},Table([snrVect,ber[3,:]])),
				  LegendEntry("64-QAM"),
				  Plot({color="blue",mark="diamond*"},Table([snrVect,ber[4,:]])),
				  LegendEntry("256-QAM"),
				  );
	# ---  Adding theoretical curve
	snrLin  = (10.0).^(snrVect/10)
	for qamScheme = qamVect
		ebNo 	= snrLin / log2(qamScheme);
		# This approximation is only valid for high SNR (one symbol error is converted to one bit error with Gray coding).
		berTheo	  = 4 * ( 1 - 1 / sqrt(qamScheme)) / log2(qamScheme) * qFunc.(sqrt.( 2*ebNo * 3 * log2(qamScheme) / (2*(qamScheme-1)  )));
		@pgf push!(a,Plot({color="black"},Table([snrVect,berTheo])));
	end
	display(a);
end

function main()
	# ----------------------------------------------------
	# --- Overall parameters
	# ----------------------------------------------------
	# --- Defining key parameters
	qamVect	        = [4,16,64,256];		# --- Constellation size
	nbSymb			= 1000;					# --- Number of symbols per MC run
	nbIt	        = 10; 					# --- MC runs
	snrVect	        = (-10:30);				# --- SNR range
	# --- Init vector
	nbSNR			= length(snrVect);
	ber				= zeros(Float64,length(qamVect),nbSNR);
	qamSeq			= zeros(Complex{Float64},nbSymb);
	qamNoise		= zeros(Complex{Float64},nbSymb);
	# --- MC run
	for iN = 1 : 1 : length(qamVect)
		# ---  Setting MCS
		mcs = qamVect[iN];
		n	= Int(log2(mcs));
		# ----------------------------------------------------
		# --- Generating data
		# ----------------------------------------------------
		# --- Calculate number of bits
		nbBits		= nbSymb *n;
		# --- Init MC buffers
		bitSeq	    = zeros(UInt8,nbBits);
		bitDec	    = zeros(UInt8,nbBits);
		# ----------------------------------------------------
		# --- Iterative BER measure
		# ----------------------------------------------------
		for k = 1 : 1 : nbSNR
			# --- Update counters
			nbC		= 0;
			nbE		= 0;
			for iN = 1 : 1 : nbIt
				# ----------------------------------------------------
				# --- Tx stage
				# ----------------------------------------------------
				# --- Create random sequence
				# Forcing seed
				genBitSequence!(bitSeq,nbBits);
				# --- QPSK mapping
				bitMappingQAM!(qamSeq,mcs,bitSeq);
				# ----------------------------------------------------
				# --- Channel
				# ----------------------------------------------------
				#  --- AWGN
				# Theoretical power is 1 (normalized constellation)
				addNoise!(qamNoise,qamSeq,snrVect[k],1);
				# ----------------------------------------------------
				# --- Rx Stage: SRRC
				# ----------------------------------------------------
				# --- Binary demapper
				bitDemappingQAM!(bitDec,mcs,qamNoise);
				# --- BER measure
				nbE	 += sum(xor.(bitDec,bitSeq));
				nbC	 += length(bitSeq);
			end
			# --- BER measure
			ber[iN,k]		= nbE / nbC;
		end
	end
	# --- Plotting routine
	doPlot(snrVect,ber,qamVect);
end

end


text:
module VectorSphericalWaves

# Vector spherical harmonics calculated using equations from:
# Mishchenko, M.I., Travis, L.D., and Lacis, A.A. (2002). Scattering, absorption, and emission of light by small particles (Cambridge University Press).

# returned vector spherical wave functions (VSWF) are matrices have the following indices:
# - index for component direction (e.g., 0,1,2 represent r,θ,ϕ)
# - index for order and rank (mn)
# - index for point in space

# inputs
# -- r,θ,ϕ : each is a 1D array, representing spactial spehrical coordinates of points in space, where VSWF will be evaluated
# -- n_max : integer, maximum rank of VSWF
# --

include("utils.jl")
include("vectorSphericalWaves_complex.jl")
include("vectorSphericalWaves_complex_all_indices.jl")
include("vectorSphericalWaves_separateRealImag.jl")

end # module


text:
import accrue.cryptoerase.runtime.Condition;

class C {
  int x = 7;
  void m() {
      final Condition c = new Condition();
      this.x = ({(L/c H)/c H})8; // Sets the level of x to H, since x can't contain an erasure policy
  }
    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase