text:
using SymbolicUtils.Rewriters

_iszero(x) = x isa Number && iszero(x)
_isone(x) = x isa Number && isone(x)
_isinteger(x) = (x isa Number && isinteger(x)) || (x isa Symbolic && symtype(x) <: Integer)
_isreal(x) = (x isa Number && isreal(x)) || (x isa Symbolic && symtype(x) <: Real)

begin
    PLUS_RULES = [
        @rule(~x::isnotflat(+) => flatten_term(+, ~x))
        @rule(~x::needs_sum_ordering => order_sum(~x))
        @ordered_acrule(~a::is_literal_number + ~b::is_literal_number => ~a + ~b)
        @rule(+(~~x::hasrepeats_qterms) => +(merge_sum_common_factors(~~x)...))
        @acrule((~z::_iszero + ~x) => ~x)
        @rule(+(~x) => ~x)
        @acrule( ~x - +(~~a, ~x, ~~b) => -(0, (~~a)..., (~~b)...) )
        @acrule( +(~~a, ~x, ~~b) + -1*~x => +(0, (~~a)..., (~~b)...) )
        @acrule( ~x + -1 * ~x => zero(symtype(~x)) )
        @acrule( *(~~x) + *(-1, ~~x) => 0 )
        @rule(+(~~a, *(~~x), ~~b, *(-1, ~~x), ~~c) => +((~~a)..., (~~b)..., (~~c)...))
    ]

    TIMES_RULES = [
        @rule(~x::isnotflat(*) => flatten_term(*, ~x))
        @rule(~x::needs_hilbert_order_sorting(*) => sort_hilbert_args(*, ~x)) ## beacuse they may not commute
        @rule(*(~~x::hasrepeats) => *(merge_repeats(^, ~~x)...))
        @rule( *(~~a, ~y, (~y)^(~n::is_literal_number), ~~b) => *((~~a)..., (~y)^(~n+1), (~~b)...) )
        @rule( *(~~a, (~y)^(~n::is_literal_number), ~y, ~~b) => *((~~a)..., (~y)^(~n+1), (~~b)...) )
        @rule( *(~~a, (~y)^(~n::is_literal_number), (~y)^(~m::is_literal_number), ~~b) => *((~~a)..., (~y)^(~n+~m), (~~b)...) )
        @ordered_acrule((~z::_isone  * ~x) => ~x)
        @ordered_acrule((~z::_iszero *  ~x) => ~z)
        @rule(*(~x) => ~x)
        @rule( *(~~x, +(~α, ~β), ~~y) => *((~~x)..., ~α, (~~y)...) + *((~~x)..., ~β, (~~y)...) )

        @rule( *(~α, +(~~y)) => +([~α * ζ for ζ in (~~y)]...) )
        @rule( *(+(~~x), ~α, +(~~y)) => +([ξ * ~α * ζ for ξ in (~~x) for ζ in (~~y)]...) )

        @rule( *(~~a, ~x::sym_isa(Number), ~~b) / ~y::sym_isa(Number) => *((~~a)..., ~x / ~y, (~~b)...) )
        @rule( *(~~a, ^(~x::istree, ~n::is_literal_number), ~~b) => *((~~a)..., [~x for i in 1:~n]..., (~~b)...) )

        @rule( *(~~a, ~x::is_boson_destroy, ~y::is_boson_create, ~~b) => *((~~a)..., swap_boson(~x, ~y), (~~b)...) )
        @rule( *(~~a, ~x::is_boson_destroy, (~y::is_boson_create)^(~n::is_literal_number), ~~b) => *((~~a)..., swap_boson(~x, ~y), (~y)^(~n - 1), (~~b)...) )
        @rule( *(~~a, (~x::is_boson_destroy)^(~n::is_literal_number), ~y::is_boson_create, ~~b) => *((~~a)..., (~x)^(~n - 1), swap_boson(~x, ~y), (~~b)...) )
        @rule( *(~~a, (~x::is_boson_destroy)^(~n::is_literal_number), (~y::is_boson_create)^(~m::is_literal_number), ~~b) => *((~~a)..., (~x)^(~n - 1), swap_boson(~x, ~y), (~y)^(~m - 1), (~~b)...) )
    ]

    POW_RULES = [
        @rule( ^(~x::istree, ~n::is_literal_number) => *([~x for i in 1:~n]...) )
        @rule( (((~x)^(~p::_isinteger))^(~q::_isinteger)) => (~x)^((~p)*(~q)) )
        @rule(^(~x, ~z::_iszero) => 1)
        @rule(^(~x, ~z::_isone) => ~x)
        @rule(inv(~x) => 1/(~x))
    ]

    ASSORTED_RULES = [
        @rule(identity(~x) => ~x)
        @rule(-(~x) => -1*~x)
        @rule(-(~x, ~y) => ~x + -1(~y))
        @rule(~x::_isone \ ~y => ~y)
        @rule(~x::_iszero \ ~y => zero(symtype(~y)))
        @rule(~x \ ~y => ~y / (~x))
        @rule(~x / ~x => one(symtype(~x)))
        @rule(~x / 1 => ~x)
        @rule(one(~x) => one(symtype(~x)))
        @rule(zero(~x) => zero(symtype(~x)))
        @rule(conj(~x::_isreal) => ~x)
        @rule(real(~x::_isreal) => ~x)
        @rule(imag(~x::_isreal) => zero(symtype(~x)))
        @rule(ifelse(~x::is_literal_number, ~y, ~z) => ~x ? ~y : ~z)
    ]

    TRIG_EXP_RULES = [
        # @acrule(~r*~x::has_trig_exp + ~r*~y => ~r*(~x + ~y))
        # @acrule(~r*~x::has_trig_exp + -1*~r*~y => ~r*(~x - ~y))
        @acrule(sin(~x)^2 + cos(~x)^2 => one(~x))
        @acrule(sin(~x)^2 + -1        => -1*cos(~x)^2)
        @acrule(cos(~x)^2 + -1        => -1*sin(~x)^2)

        @acrule(cos(~x)^2 + -1*sin(~x)^2 => cos(2 * ~x))
        @acrule(sin(~x)^2 + -1*cos(~x)^2 => -cos(2 * ~x))
        @acrule(cos(~x) * sin(~x) => sin(2 * ~x)/2)

        @acrule(tan(~x)^2 + -1*sec(~x)^2 => one(~x))
        @acrule(-1*tan(~x)^2 + sec(~x)^2 => one(~x))
        @acrule(tan(~x)^2 +  1 => sec(~x)^2)
        @acrule(sec(~x)^2 + -1 => tan(~x)^2)

        @acrule(cot(~x)^2 + -1*csc(~x)^2 => one(~x))
        @acrule(cot(~x)^2 +  1 => csc(~x)^2)
        @acrule(csc(~x)^2 + -1 => cot(~x)^2)

        # @acrule(exp(~x) * exp(~y) => _iszero(~x + ~y) ? 1 : exp(~x + ~y))
        # @rule(exp(~x)^(~y) => exp(~x * ~y))
    ]

    BOOLEAN_RULES = [
        @rule((true | (~x)) => true)
        @rule(((~x) | true) => true)
        @rule((false | (~x)) => ~x)
        @rule(((~x) | false) => ~x)
        @rule((true & (~x)) => ~x)
        @rule(((~x) & true) => ~x)
        @rule((false & (~x)) => false)
        @rule(((~x) & false) => false)

        @rule(!(~x) & ~x => false)
        @rule(~x & !(~x) => false)
        @rule(!(~x) | ~x => true)
        @rule(~x | !(~x) => true)
        @rule(xor(~x, !(~x)) => true)
        @rule(xor(~x, ~x) => false)

        @rule(~x == ~x => true)
        @rule(~x != ~x => false)
        @rule(~x < ~x => false)
        @rule(~x > ~x => false)

        # simplify terms with no symbolic arguments
        # e.g. this simplifies term(isodd, 3, type=Bool)
        # or term(!, false)
        @rule((~f)(~x::is_literal_number) => (~f)(~x))
        # and this simplifies any binary comparison operator
        @rule((~f)(~x::is_literal_number, ~y::is_literal_number) => (~f)(~x, ~y))
    ]

    function quantum_normal_order_simplifier()
        rule_tree = [If(istree, Chain(ASSORTED_RULES)),
                     If(is_operation(+),
                        Chain(PLUS_RULES)),
                     If(is_operation(*),
                        Chain(TIMES_RULES)),
                     If(is_operation(^),
                        Chain(POW_RULES))] |> RestartedChain

        rule_tree
    end

    trig_exp_simplifier(;kw...) = Chain(TRIG_EXP_RULES);

    bool_simplifier() = Chain(BOOLEAN_RULES);

    global normal_order_simplifier
    global serial_normal_order_simplifier
    global threaded_normal_order_simplifier

    function normal_order_simplifier(; kw...)
        IfElse(has_trig_exp,
               Postwalk(IfElse(x->symtype(x) <: QOperator,
                               Chain((quantum_normal_order_simplifier(),
                                      trig_exp_simplifier())),
                               If(x->symtype(x) <: Bool,
                                  bool_simplifier()))
                        ; kw...),
               Postwalk(Chain((If(x->symtype(x) <: QOperator,
                                  quantum_normal_order_simplifier()),
                               If(x->symtype(x) <: Bool,
                                  bool_simplifier())))
                        ; kw...))
    end

    serial_normal_order_simplifier = If(istree, Fixpoint(normal_order_simplifier()))

    threaded_normal_order_simplifier(cutoff) = Fixpoint(normal_order_simplifier(threaded=true,
                                                              thread_cutoff=cutoff));
end


text:
#! /usr/bin/julia

# Rosetta Code, Four bit adder

xor{T<:Bool}(a::T, b::T) = (a&~b)|(~a&b)

halfadder{T<:Bool}(a::T, b::T) = (xor(a,b), a&b)

function fulladder{T<:Bool}(a::T, b::T, c::T=false)
    (s, ca) = halfadder(c, a)
    (s, cb) = halfadder(s, b)
    (s, ca|cb)
end

function adder(a::BitArray{1}, b::BitArray{1}, c0::Bool=false)
    len = length(a)
    length(b) == len || error("Addend width mismatch.")
    c = c0
    s = falses(len)
    for i in 1:len
        (s[i], c) = fulladder(a[i], b[i], c)
    end
    (s, c)
end

function adder{T<:Integer}(m::T, n::T, wid::T=4, c0::Bool=false)
    a = bitpack(digits(m, 2, wid))[1:wid]
    b = bitpack(digits(n, 2, wid))[1:wid]
    adder(a, b, c0)
end

Base.bits(n::BitArray{1}) = join(reverse(int(n)), "")

xavail = trues(15,15)
xcnt = 0
xgoal = 10
println("Testing adder with 4-bit words:")
while xcnt < xgoal
    m = rand(1:15)
    n = rand(1:15)
    xavail[m,n] || continue
    xavail[m,n] = xavail[n,m] = false
    xcnt += 1
    (s, c) = adder(m, n)
    oflow = c ? "*" : ""
    print(@sprintf "    %2d + %2d = %2d => " m n m+n)
    println(@sprintf("%s + %s = %s%s",
                     bits(m)[end-3:end],
                     bits(n)[end-3:end],
                     bits(s), oflow))
end

ext: jl

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using LinearAlgebra: eigvecs, diag, diagind, dot
using StatsBase: sample, Weights
using RingPolymerArrays: get_centroid

using NQCDynamics.Calculators: RingPolymerDiabaticCalculator
using NQCDynamics: RingPolymerSimulation, RingPolymers

function RingPolymerSimulation{FSSH}(atoms::Atoms{T}, model::Model, n_beads::Integer; rescaling=:standard, kwargs...) where {T}
    RingPolymerSimulation(atoms, model, FSSH{T}(NQCModels.nstates(model), rescaling), n_beads; kwargs...)
end

function DynamicsMethods.motion!(du, u, sim::RingPolymerSimulation{<:SurfaceHopping}, t)
    dr = DynamicsUtils.get_positions(du)
    dv = DynamicsUtils.get_velocities(du)
    dσ = DynamicsUtils.get_quantum_subsystem(du)

    r = DynamicsUtils.get_positions(u)
    v = DynamicsUtils.get_velocities(u)
    σ = DynamicsUtils.get_quantum_subsystem(u)

    set_state!(u, sim.method.state) # Make sure the state variables match,

    DynamicsUtils.velocity!(dr, v, r, sim, t)
    Calculators.update_electronics!(sim.calculator, r)
    acceleration!(dv, v, r, sim, t, sim.method.state)
    DynamicsUtils.apply_interbead_coupling!(dv, r, sim)
    DynamicsUtils.set_quantum_derivative!(dσ, v, σ, sim)
end

function get_hopping_nonadiabatic_coupling(sim::RingPolymerSimulation{<:FSSH})
    sim.calculator.centroid_nonadiabatic_coupling
end

function get_hopping_velocity(::RingPolymerSimulation{<:FSSH}, u)
    get_centroid(DynamicsUtils.get_velocities(u))
end

function get_hopping_eigenvalues(sim::RingPolymerSimulation{<:FSSH})
    sim.calculator.centroid_eigen.values
end

function perform_rescaling!(
    sim::RingPolymerSimulation{<:SurfaceHopping}, velocity, γ, d
)
    for I in CartesianIndices(d)
        velocity[I,:] .-= γ * d[I] / masses(sim, I)
    end
    return nothing
end

function DynamicsUtils.classical_hamiltonian(sim::RingPolymerSimulation{<:FSSH}, u)
    kinetic = DynamicsUtils.classical_kinetic_energy(sim, DynamicsUtils.get_velocities(u))
    spring = RingPolymers.get_spring_energy(sim.beads, sim.atoms.masses, DynamicsUtils.get_positions(u))

    all_eigs = Calculators.get_eigen(sim.calculator, DynamicsUtils.get_positions(u))
    potential = sum(eigs.values[u.state] for eigs in all_eigs)
    return kinetic + spring + potential
end


text:
type Source
     isz  :: Int64
     isx  :: Int64
     nz   :: Int64
     nx   :: Int64
     ext  :: Int64
     iflag:: Int64
     ot   :: Float64
     dt   :: Float64
     nt   :: Int64
     Vx_flag  :: Bool
     Vz_flag  :: Bool
     Txx_flag :: Bool
     Tzz_flag :: Bool
     Txz_flag :: Bool
     Vx  :: Array{Float64,1}
     Vz  :: Array{Float64,1}
     Txx :: Array{Float64,1}
     Tzz :: Array{Float64,1}
     Txz :: Array{Float64,1}
end

function Ricker(f0::Float64, dt::Float64)
	  nw = 2.2/f0/dt
	  nw = 2*floor(Int,nw/2)+1
	  nc = floor(Int,nw/2)
	  w  = zeros(nw)
	  k  = collect(1:nw)
    k  = vec(k)
	  alpha = (nc-k+1)*f0*dt*pi
	  beta = alpha.^2
	  w = ((1.-beta.*2).*exp(-beta))
	  return w
end

function InitSource(isz::Int64, isx::Int64, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, f0::Float64, ot::Float64, dt::Float64, flags::Array{Bool,1})
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    if flags[1]
       w_Vx = Ricker(f0, dt)
       nt = length(w_Vx)
    else
       w_Vx = zeros(1)
    end
    if flags[2]
       w_Vz = Ricker(f0, dt)
       nt = length(w_Vz)
    else
       w_Vz = zeros(1)
    end
    if flags[3]
       w_Txx = Ricker(f0, dt)
       nt = length(w_Txx)
    else
       w_Txx = zeros(1)
    end
    if flags[4]
       w_Tzz = Ricker(f0, dt)
       nt = length(w_Tzz)
    else
       w_Tzz = zeros(1)
    end
    if flags[5]
       w_Txz = Ricker(f0, dt)
       nt = length(w_Txz)
    else
       w_Txz = zeros(1)
    end
    src = Source(isz, isx, nz, nx, ext, iflag, ot, dt, nt,
                 flags[1], flags[2], flags[3], flags[4], flags[5],
                 w_Vx, w_Vz, w_Txx, w_Tzz, w_Txz)
    return src
end

function InitMultiSources(isz::Array{Int64,1}, isx::Array{Int64,1}, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, f0::Float64, ot::Array{Float64,1}, dt::Float64, flags::Array{Bool,1})
    ns = length(isz)
    srcs = Array(Source, ns)
    for is = 1 : ns
        srcs[is] = InitSource(isz[is], isx[is], nz, nx, ext, iflag, f0, ot[is], dt, flags)
    end
    return srcs
end

function SrcRange(src::Source)
    tl = src.ot
    tu = src.ot + (src.nt-1)*src.dt
    return tl, tu
end

function SourcesTimeRange(srcs::Array{Source,1})
    tmax = 0.0
    tmin = 0.0
    ns = length(srcs)
    dt = srcs[1].dt
    for is = 1 : ns
        tmp = srcs[is].ot + (srcs[is].nt-1)*dt
        tmp1= srcs[is].ot
        if tmp > tmax
           tmax = tmp
        end
        if tmp1 < tmin
           tmin = tmp1
        end
    end
    return tmin, tmax
end

function AddSource!(spt::SnapShot, src::Source)
    isz= src.isz
    isx= src.isx
    nz = src.nz
    nx = src.nx
    ext= src.ext
    iflag = src.iflag
    dt = src.dt
    if iflag == 1
       Nz = nz + 2*ext
       upper = ext
    elseif iflag ==2
       Nz = nz +   ext
       upper = 0
    end
    Nx = nx + 2*ext
    it  = spt.it
    if src.ot <= (it-1)*dt <= src.ot+(src.nt-1)*dt
       indt = it - round(Int64, src.ot/dt)
       indz = src.isz + upper
       indx = src.isx + ext
       pos = (indx-1) * Nz + indz
       if src.Vx_flag
          spt.Vxx[pos] = spt.Vxx[pos] + src.Vx[indt]/2
          spt.Vxz[pos] = spt.Vxz[pos] + src.Vx[indt]/2
       end
       if src.Vz_flag
          spt.Vzx[pos] = spt.Vzx[pos] + src.Vz[indt]/2
          spt.Vzz[pos] = spt.Vzz[pos] + src.Vz[indt]/2
       end
       if src.Txx_flag
          spt.Txxx[pos] = spt.Txxx[pos] + src.Txx[indt]/2
          spt.Txxz[pos] = spt.Txxz[pos] + src.Txx[indt]/2
       end
       if src.Tzz_flag
          spt.Tzzx[pos] = spt.Tzzx[pos] + src.Tzz[indt]/2
          spt.Tzzz[pos] = spt.Tzzz[pos] + src.Tzz[indt]/2
       end
       if src.Txz_flag
          spt.Txzx[pos] = spt.Txzx[pos] + src.Txz[indt]/2
          spt.Txzz[pos] = spt.Txzz[pos] + src.Txz[indt]/2
       end
    end
    return nothing
end

function AddMultiSources!(spt::SnapShot, srcs::Array{Source,1})
    nz    = srcs[1].nz
    nx    = srcs[1].nx
    ext   = srcs[1].ext
    iflag = srcs[1].iflag
    dt = srcs[1].dt
    if iflag == 1
       Nz = nz + 2*ext
       upper = ext
    elseif iflag ==2
       Nz = nz +   ext
       upper = 0
    end
    Nx = nx + 2*ext
    it = spt.it
    ns = length(srcs)
    for is = 1: ns
        if srcs[is].ot <= (it-1)*dt <= srcs[is].ot+(srcs[is].nt-1)*dt
           indt = it - round(Int64, srcs[is].ot/dt)
           indz = srcs[is].isz + upper
           indx = srcs[is].isx + ext
           pos  = (indx-1) * Nz + indz
           if srcs[is].Vx_flag
              spt.Vxx[pos] = spt.Vxx[pos] + srcs[is].Vx[indt]/2
              spt.Vxz[pos] = spt.Vxz[pos] + srcs[is].Vx[indt]/2
           end
           if srcs[is].Vz_flag
              spt.Vzx[pos] = spt.Vzx[pos] + srcs[is].Vz[indt]/2
              spt.Vzz[pos] = spt.Vzz[pos] + srcs[is].Vz[indt]/2
           end
           if srcs[is].Txx_flag
              spt.Txxx[pos] = spt.Txxx[pos] + srcs[is].Txx[indt]/2
              spt.Txxz[pos] = spt.Txxz[pos] + srcs[is].Txx[indt]/2
           end
           if srcs[is].Tzz_flag
              spt.Tzzx[pos] = spt.Tzzx[pos] + srcs[is].Tzz[indt]/2
              spt.Tzzz[pos] = spt.Tzzz[pos] + srcs[is].Tzz[indt]/2
           end
           if srcs[is].Txz_flag
              spt.Txzx[pos] = spt.Txzx[pos] + srcs[is].Txz[indt]/2
              spt.Txzz[pos] = spt.Txzz[pos] + srcs[is].Txz[indt]/2
           end
        end
    end
    return nothing
end

function srcs2spt(path::String, srcs::Array{Source,1})
    (tl, tu) = SourcesTimeRange(srcs)
    nz = srcs[1].nz; nx = srcs[1].nx; ext = srcs[1].ext; iflag = srcs[1].iflag;
    dt = srcs[1].dt; nt = floor(Int64, tu/dt)+1;
    spt = InitSnapShot(nz, nx, ext, iflag, dt, 1)
    AddMultiSources!(spt, srcs)
    fid = WriteSnapShot(path, spt)
    for it = 2 : nt
        spt = InitSnapShot(nz, nx, ext, iflag, dt, it)
        AddMultiSources!(spt, srcs)
        WriteSnapShot(fid, spt)
    end
    close(fid)
    return nothing
end

function srcs2wfd(path::String, srcs::Array{Source,1})
    (tl, tu) = SourcesTimeRange(srcs)
    nz = srcs[1].nz; nx = srcs[1].nx; ext = srcs[1].ext; iflag = srcs[1].iflag;
    dt = srcs[1].dt; nt = floor(Int64, tu/dt)+1;
    wfd = InitWfd(nz, nx, ext, iflag, dt, 1)
    AddSrcs2Wfd!(wfd, srcs)
    fid = WriteSnapShot(path, spt)
    for it = 2 : nt
        spt = InitSnapShot(nz, nx, ext, iflag, dt, it)
        AddMultiSources!(spt, srcs)
        WriteSnapShot(fid, spt)
    end
    close(fid)
    return nothing
end

function AddSourceBorn!(spt::SnapShot, path::String, dm::Array{Float64,1}, du::Array{Float64,1})
    it = spt.it;
    nz = spt.nz; nx = spt.nx; ext = spt.ext;
    if spt.iflag == 1
       Nz = nz + 2*ext
       upper = ext
     elseif spt.iflag == 2
       Nz = nz +   ext
       upper = 0
    end
    Nx = nx + 2*ext
    (vxx, vxz, vzx, vzz) = ReadPv(path, it)
    pc = vxx+vzz
    sc = vxz+vzx
    # add source to Txx
    tmp = 1/2*(pc.*dm+2*vxx.*du);
    spt.Txxx = spt.Txxx + tmp
    spt.Txxz = spt.Txxz + tmp
    # add source to Tzz
    tmp = 1/2*(pc.*dm+2*vzz.*du)
    spt.Tzzx = spt.Tzzx + tmp
    spt.Tzzz = spt.Tzzz + tmp
    # add source to Txz
    tmp = 1/2*(sc.*du)
    spt.Txzx = spt.Txzx + tmp
    spt.Txzz = spt.Txzz + tmp
    return nothing
end
# function AddSourceBorn!(spt::SnapShot, path::String, dm::Array{Float64,2}, du::Array{Float64,2})
#     it = spt.it;
#     nz = spt.nz; nx = spt.nx; ext = spt.ext;
#     if spt.iflag == 1
#        Nz = nz + 2*ext
#        upper = ext
#      elseif spt.iflag == 2
#        Nz = nz +   ext
#        upper = 0
#     end
#     Nx = nx + 2*ext
#     # add source to Txx
#     tmp = zeros(Nz, Nx)
#     (vxx, vxz, vzx, vzz) = ReadPv(path, it)
#     tmp[upper+1:nz+upper, ext+1:nx+ext] = 1/2*((vxx+vzz).*dm+2*vxx.*du)
#     tmp = vec(tmp)
#     spt.Txxx = spt.Txxx + tmp
#     spt.Txxz = spt.Txxz + tmp
#     # add source to Tzz
#     tmp = zeros(Nz, Nx)
#     tmp[upper+1:nz+upper, ext+1:nx+ext] = 1/2*((vxx+vzz).*dm+2*vzz.*du)
#     tmp = vec(tmp)
#     spt.Tzzx = spt.Tzzx + tmp
#     spt.Tzzz = spt.Tzzz + tmp
#     # add source to Txz
#     tmp = zeros(Nz, Nx)
#     tmp[upper+1:nz+upper, ext+1:nx+ext] = 1/2*((vxz+vzx).*du)
#     tmp = vec(tmp)
#     spt.Txzx = spt.Txzx + tmp
#     spt.Txzz = spt.Txzz + tmp
#     return nothing
# end
function ConvertBornSource2Spts(path_spt::String, path_pv::String, dm::Array{Float64,1}, du::Array{Float64,1})
    (nz, nx, ext, iflag, dt, ns) = InfoPv(path_pv)
    if iflag == 1
       Nz = nz + 2*ext
       upper = ext
     elseif iflag == 2
       Nz = nz +   ext
       upper = 0
    end
    Nx = nx + 2*ext
    # add source to Txx
    fid = open(path_spt, "w");
    write(fid, nz, nx, ext, iflag, dt)
    for it = 1 : ns
        (vxx, vxz, vzx, vzz) = ReadPv(path_pv, it)
        pc = vxx+vzz
        Txx = 1/2*(pc.*dm + 2*vxx.*du)
        Tzz = 1/2*(pc.*dm + 2*vzz.*du)
        Txz = 1/2*(sc.*du)
        write(fid, zeros(Nz*Nx))
        write(fid, zeros(Nz*Nx))
        write(fid, zeros(Nz*Nx))
        write(fid, zeros(Nz*Nx))
        write(fid, Txx         )
        write(fid, Txx         )
        write(fid, Tzz         )
        write(fid, Tzz         )
        write(fid, Txz         )
        write(fid, Txz         )
    end
    close(fid)
    return nothing
end
# function ConvertBornSource2Spts(path_spt::String, path_pv::String, dm::Array{Float64,2}, du::Array{Float64,2})
#     (nz, nx, ext, iflag, dt, ns) = InfoPv(path_pv)
#     if iflag == 1
#        Nz = nz + 2*ext
#        upper = ext
#      elseif iflag == 2
#        Nz = nz +   ext
#        upper = 0
#     end
#     Nx = nx + 2*ext
#     # add source to Txx
#     fid = open(path_spt, "w");
#     write(fid, nz, nx, ext, iflag, dt)
#     for it = 1 : ns
#         (vxx, vxz, vzx, vzz) = ReadPv(path_pv, it)
#         Txx = zeros(Nz, Nx)
#         Txx[upper+1:nz+upper, ext+1:nx+ext] = 1/2*((vxx+vzz).*dm + 2*vxx.*du)
#         Txx = vec(Txx)
#         Tzz = zeros(Nz, Nx)
#         Tzz[upper+1:nz+upper, ext+1:nx+ext] = 1/2*((vxx+vzz).*dm + 2*vzz.*du)
#         Tzz = vec(Tzz)
#         Txz = zeros(Nz, Nx)
#         Txz[upper+1:nz+upper, ext+1:nx+ext] = 1/2*((vxz+vzx).*du)
#         Txz = vec(Txz)
#         write(fid, zeros(Nz*Nx))
#         write(fid, zeros(Nz*Nx))
#         write(fid, zeros(Nz*Nx))
#         write(fid, zeros(Nz*Nx))
#         write(fid, Txx         )
#         write(fid, Txx         )
#         write(fid, Tzz         )
#         write(fid, Tzz         )
#         write(fid, Txz         )
#         write(fid, Txz         )
#     end
#     close(fid)
#     return nothing
# end

max_stars_repo_name: WneleiGao/ElasticWave


text:
plot(exp, 0, 1)
plot!(sqrt, 0, 1, label="sqrt")
legend()


text:
using ArgParse,CairoMakie

CairoMakie.activate!(type="pdf")

function parse_argument()
end

function hypergeo_prob_desity(n,M,N,x)
    binomial(BigInt(M),BigInt(x))*binomial(BigInt(N),BigInt(n-x))/binomial(BigInt(M+N),BigInt(n))
end

function main()
    graph_name = "hypergeo_cumulative_dist"

    xs_max = 50
    xs = 0:xs_max

    graph_params = [Dict("params"=>Dict("n"=>10,"N"=>100,"M"=>100), "color"=>:gray0,  "marker"=>:circle, "label"=>"n=10,M=100,N=100"),
                    Dict("params"=>Dict("n"=>30,"N"=>100,"M"=>100), "color"=>:gray30, "marker"=>:diamond, "label"=>"n=30,M=100,N=100"),
                    Dict("params"=>Dict("n"=>70,"N"=>100,"M"=>100), "color"=>:gray50, "marker"=>:rect,    "label"=>"n=70,M=100,N=100")]

    fig = Figure(resolution = (800, 600))

    axis = Axis(fig[1,1], title = "hypergeometric cumulative distribution function")

    for params in graph_params
        c_s::Array{Float64} = []
        accum::Float64 = 0
        n = params["params"]["n"]
        M = params["params"]["M"]
        N = params["params"]["N"]
        p_s = [hypergeo_prob_desity(n,M,N,x) for x in xs ]
        for p in p_s
            accum += p
            push!(c_s, accum)
        end
        stem!(fig[1,1], xs, c_s, color = params["color"], marker = params["marker"], label=params["label"])
    end

    axis.xlabelsize=12
    axislegend(axis, position = :rb)

    save("$graph_name.pdf", fig)

end

main()

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
@testset "Hearn Problems" begin
    (a, alpha, b, c, d, e, eps, epsilon, h, k, m, mc, p, q, r, x, y, z, ) = @variables a alpha b c d e eps epsilon h k m mc p q r x y z

    #= ::Package:: =#

    #= ::Title:: =#
    #=Anthony*Hearn - Reduce*Integration*Test*Package=#

    #= ::Section::Closed:: =#
    #=Polynomial*and*rational*function*examples=#

    @test_int [1 + x + x^2, x, 1, x + x^2/2 + x^3/3]
    @test_int [x^2*(2*x^2 + x)^2, x, 3, x^5/5 + (2*x^6)/3 + (4*x^7)/7]
    @test_int [x*(x^2 + 2*x + 1), x, 2, x^2/2 + (2*x^3)/3 + x^4/4]
    @test_int [1/x, x, 1, log(x)]
    @test_int [(x + 1)^3/(x - 1)^4, x, 2, 8/(3*(1 - x)^3) - 6/(1 - x)^2 + 6/(1 - x) + log(1 - x)]
    @test_int [1/(x*(x - 1)*(x + 1)^2), x, 2, -(1/(2*(1 + x))) + (1/4)*log(1 - x) - log(x) + (3/4)*log(1 + x)]
    @test_int [(a*x + b)/((x - p)*(x - q)), x, 2, ((b + a*p)*log(p - x))/(p - q) - ((b + a*q)*log(q - x))/(p - q)]
    @test_int [1/(a*x^2 + b*x + c), x, 2, -((2*atanh((b + 2*a*x)/sqrt(b^2 - 4*a*c)))/sqrt(b^2 - 4*a*c))]
    @test_int [(a*x + b)/(1 + x^2), x, 3, b*atan(x) + (1/2)*a*log(1 + x^2)]
    @test_int [1/(x^2 - 2*x + 3), x, 2, -(atan((1 - x)/sqrt(2))/sqrt(2))]

    #= ::Section::Closed:: =#
    #=Rational*function*examples*from*Hardy, Pure*Mathematics, p*253*et*seq.=#

    @test_int [1/((x-1)*(x^2+1))^2, x, 6, 1/(4*(1 - x)) - 1/(4*(1 + x^2)) + atan(x)/4 - (1/2)*log(1 - x) + (1/4)*log(1 + x^2)]
    @test_int [x/((x-a)*(x-b)*(x-c)), x, 2, (a*log(a - x))/((a - b)*(a - c)) - (b*log(b - x))/((a - b)*(b - c)) + (c*log(c - x))/((a - c)*(b - c))]
    @test_int [x/((x^2+a^2)*(x^2+b^2)), x, 4, -(log(a^2 + x^2)/(2*(a^2 - b^2))) + log(b^2 + x^2)/(2*(a^2 - b^2))]
    @test_int [x^2/((x^2+a^2)*(x^2+b^2)), x, 3, (a*atan(x/a))/(a^2 - b^2) - (b*atan(x/b))/(a^2 - b^2)]
    @test_int [x/((x-1)*(x^2+1)), x, 5, atan(x)/2 + (1/2)*log(1 - x) - (1/4)*log(1 + x^2)]
    @test_int [x/(1+x^3), x, 6, -(atan((1 - 2*x)/sqrt(3))/sqrt(3)) - (1/3)*log(1 + x) + (1/6)*log(1 - x + x^2)]
    @test_int [x^3/((x-1)^2*(x^3+1)), x, 3, 1/(2*(1 - x)) + (3/4)*log(1 - x) - (1/12)*log(1 + x) - (1/3)*log(1 - x + x^2)]
    @test_int [1/(1+x^4), x, 9, -(atan(1 - sqrt(2)*x)/(2*sqrt(2))) + atan(1 + sqrt(2)*x)/(2*sqrt(2)) - log(1 - sqrt(2)*x + x^2)/(4*sqrt(2)) + log(1 + sqrt(2)*x + x^2)/(4*sqrt(2))]
    @test_int [x^2/(1+x^4), x, 9, -(atan(1 - sqrt(2)*x)/(2*sqrt(2))) + atan(1 + sqrt(2)*x)/(2*sqrt(2)) + log(1 - sqrt(2)*x + x^2)/(4*sqrt(2)) - log(1 + sqrt(2)*x + x^2)/(4*sqrt(2))]
    @test_int [1/(1+x^2+x^4), x, 9, -(atan((1 - 2*x)/sqrt(3))/(2*sqrt(3))) + atan((1 + 2*x)/sqrt(3))/(2*sqrt(3)) - (1/4)*log(1 - x + x^2) + (1/4)*log(1 + x + x^2)]

    #= ::Section::Closed:: =#
    #=Examples*involving*a+b*x=#

    @test_int [(a+b*x)^p, x, 1, (a + b*x)^(1 + p)/(b*(1 + p))]
    @test_int [x*(a+b*x)^p, x, 2, -((a*(a + b*x)^(1 + p))/(b^2*(1 + p))) + (a + b*x)^(2 + p)/(b^2*(2 + p))]
    @test_int [x^2*(a+b*x)^p, x, 2, (a^2*(a + b*x)^(1 + p))/(b^3*(1 + p)) - (2*a*(a + b*x)^(2 + p))/(b^3*(2 + p)) + (a + b*x)^(3 + p)/(b^3*(3 + p))]
    @test_int [1/(a+b*x), x, 1, log(a + b*x)/b]
    @test_int [1/(a+b*x)^2, x, 1, -(1/(b*(a + b*x)))]
    @test_int [x/(a+b*x), x, 2, x/b - (a*log(a + b*x))/b^2]
    @test_int [x^2/(a+b*x), x, 2, -((a*x)/b^2) + x^2/(2*b) + (a^2*log(a + b*x))/b^3]
    @test_int [1/(x*(a+b*x)), x, 3, log(x)/a - log(a + b*x)/a]
    @test_int [1/(x^2*(a+b*x)), x, 2, -(1/(a*x)) - (b*log(x))/a^2 + (b*log(a + b*x))/a^2]
    @test_int [1/(x*(a+b*x))^2, x, 2, -(1/(a^2*x)) - b/(a^2*(a + b*x)) - (2*b*log(x))/a^3 + (2*b*log(a + b*x))/a^3]
    @test_int [1/(c^2+x^2), x, 1, atan(x/c)/c]
    @test_int [1/(c^2-x^2), x, 1, atanh(x/c)/c]

    #= ::Section::Closed:: =#
    #=More*complicated*rational*function*examples=#

    #= Mostly*examples*contributed*by*David*M. Dahm, who*also*developed*the*code*to*integrate*them. =#

    @test_int [1/(2*x^3-1), x, 6, -(atan((1 + 2*2^(1/3)*x)/sqrt(3))/(2^(1/3)*sqrt(3))) + log(1 - 2^(1/3)*x)/(3*2^(1/3)) - log(1 + 2^(1/3)*x + 2^(2/3)*x^2)/(6*2^(1/3))]
    @test_int [1/(x^3-2), x, 6, -(atan((1 + 2^(2/3)*x)/sqrt(3))/(2^(2/3)*sqrt(3))) + log(2^(1/3) - x)/(3*2^(2/3)) - log(2^(2/3) + 2^(1/3)*x + x^2)/(6*2^(2/3))]
    @test_int [1/(a*x^3-b), x, 6, -(atan((b^(1/3) + 2*a^(1/3)*x)/(sqrt(3)*b^(1/3)))/(sqrt(3)*a^(1/3)*b^(2/3))) + log(b^(1/3) - a^(1/3)*x)/(3*a^(1/3)*b^(2/3)) - log(b^(2/3) + a^(1/3)*b^(1/3)*x + a^(2/3)*x^2)/(6*a^(1/3)*b^(2/3))]
    @test_int [1/(x^4-2), x, 3, -(atan(x/2^(1/4))/(2*2^(3/4))) - atanh(x/2^(1/4))/(2*2^(3/4))]
    @test_int [1/(5*x^4-1), x, 3, -(atan(5^(1/4)*x)/(2*5^(1/4))) - atanh(5^(1/4)*x)/(2*5^(1/4))]
    @test_int [1/(3*x^4+7), x, 9, If(13<9, -(atan(1 - (3/7)^(1/4)*sqrt(2)*x)/(2*sqrt(2)*3^(1/4)*7^(3/4))) + atan(1 + (3/7)^(1/4)*sqrt(2)*x)/(2*sqrt(2)*3^(1/4)*7^(3/4)) - log(sqrt(21) - sqrt(2)*3^(3/4)*7^(1/4)*x + 3*x^2)/(4*sqrt(2)*3^(1/4)*7^(3/4)) + log(sqrt(21) + sqrt(2)*3^(3/4)*7^(1/4)*x + 3*x^2)/(4*sqrt(2)*3^(1/4)*7^(3/4)), -(atan(1 - (3/7)^(1/4)*sqrt(2)*x)/(2*sqrt(2)*3^(1/4)*7^(3/4))) + atan(1 + (3/7)^(1/4)*sqrt(2)*x)/(2*sqrt(2)*3^(1/4)*7^(3/4)) - log(sqrt(21) - sqrt(2)*3^(3/4)*7^(1/4)*x + 3*x^2)/(4*sqrt(2)*3^(1/4)*7^(3/4)) + log(sqrt(21) + sqrt(2)*3^(3/4)*7^(1/4)*x + 3*x^2)/(4*sqrt(2)*3^(1/4)*7^(3/4)))]
    @test_int [1/(x^4+3*x^2-1), x, 3, (-sqrt(2/(13*(3 + sqrt(13)))))*atan(sqrt(2/(3 + sqrt(13)))*x) - sqrt((1/26)*(3 + sqrt(13)))*atanh(sqrt(2/(-3 + sqrt(13)))*x)]
    @test_int [1/(x^4-3*x^2-1), x, 3, (-sqrt((1/26)*(3 + sqrt(13))))*atan(sqrt(2/(-3 + sqrt(13)))*x) - sqrt(2/(13*(3 + sqrt(13))))*atanh(sqrt(2/(3 + sqrt(13)))*x)]
    @test_int [1/(x^4-3*x^2+1), x, 3, (-sqrt(2/(5*(3 + sqrt(5)))))*atanh(sqrt(2/(3 + sqrt(5)))*x) + sqrt((1/10)*(3 + sqrt(5)))*atanh(sqrt((1/2)*(3 + sqrt(5)))*x)]
    @test_int [1/(x^4-4*x^2+1), x, 3, atanh(x/sqrt(2 - sqrt(3)))/(2*sqrt(3*(2 - sqrt(3)))) - atanh(x/sqrt(2 + sqrt(3)))/(2*sqrt(3*(2 + sqrt(3))))]
    @test_int [1/(x^4+4*x^2+1), x, 3, atan(x/sqrt(2 - sqrt(3)))/(2*sqrt(3*(2 - sqrt(3)))) - atan(x/sqrt(2 + sqrt(3)))/(2*sqrt(3*(2 + sqrt(3))))]
    @test_int [1/(x^4+x^2+2), x, 9, (-(1/2))*sqrt((1/14)*(-1 + 2*sqrt(2)))*atan((sqrt(-1 + 2*sqrt(2)) - 2*x)/sqrt(1 + 2*sqrt(2))) + (1/2)*sqrt((1/14)*(-1 + 2*sqrt(2)))*atan((sqrt(-1 + 2*sqrt(2)) + 2*x)/sqrt(1 + 2*sqrt(2))) - log(sqrt(2) - sqrt(-1 + 2*sqrt(2))*x + x^2)/(4*sqrt(2*(-1 + 2*sqrt(2)))) + log(sqrt(2) + sqrt(-1 + 2*sqrt(2))*x + x^2)/(4*sqrt(2*(-1 + 2*sqrt(2))))]
    @test_int [1/(x^4-x^2+2), x, 9, (-(1/2))*sqrt((1/14)*(1 + 2*sqrt(2)))*atan((sqrt(1 + 2*sqrt(2)) - 2*x)/sqrt(-1 + 2*sqrt(2))) + (1/2)*sqrt((1/14)*(1 + 2*sqrt(2)))*atan((sqrt(1 + 2*sqrt(2)) + 2*x)/sqrt(-1 + 2*sqrt(2))) - log(sqrt(2) - sqrt(1 + 2*sqrt(2))*x + x^2)/(4*sqrt(2*(1 + 2*sqrt(2)))) + log(sqrt(2) + sqrt(1 + 2*sqrt(2))*x + x^2)/(4*sqrt(2*(1 + 2*sqrt(2))))]
    @test_int [1/(x^6-1), x, 10, atan((1 - 2*x)/sqrt(3))/(2*sqrt(3)) - atan((1 + 2*x)/sqrt(3))/(2*sqrt(3)) - atanh(x)/3 + (1/12)*log(1 - x + x^2) - (1/12)*log(1 + x + x^2)]
    @test_int [1/(x^6-2), x, 10, atan(1/sqrt(3) - (2^(5/6)*x)/sqrt(3))/(2*2^(5/6)*sqrt(3)) - atan(1/sqrt(3) + (2^(5/6)*x)/sqrt(3))/(2*2^(5/6)*sqrt(3)) - atanh(x/2^(1/6))/(3*2^(5/6)) + log(2^(1/3) - 2^(1/6)*x + x^2)/(12*2^(5/6)) - log(2^(1/3) + 2^(1/6)*x + x^2)/(12*2^(5/6))]
    @test_int [1/(x^6+2), x, 10, atan(x/2^(1/6))/(3*2^(5/6)) - atan(sqrt(3) - 2^(5/6)*x)/(6*2^(5/6)) + atan(sqrt(3) + 2^(5/6)*x)/(6*2^(5/6)) - log(2^(1/3) - 2^(1/6)*sqrt(3)*x + x^2)/(4*2^(5/6)*sqrt(3)) + log(2^(1/3) + 2^(1/6)*sqrt(3)*x + x^2)/(4*2^(5/6)*sqrt(3))]
    @test_int [1/(x^8+1), x, 19, -(atan((sqrt(2 - sqrt(2)) - 2*x)/sqrt(2 + sqrt(2)))/(4*sqrt(2*(2 - sqrt(2))))) - atan((sqrt(2 + sqrt(2)) - 2*x)/sqrt(2 - sqrt(2)))/(4*sqrt(2*(2 + sqrt(2)))) + atan((sqrt(2 - sqrt(2)) + 2*x)/sqrt(2 + sqrt(2)))/(4*sqrt(2*(2 - sqrt(2)))) + atan((sqrt(2 + sqrt(2)) + 2*x)/sqrt(2 - sqrt(2)))/(4*sqrt(2*(2 + sqrt(2)))) - (1/16)*sqrt(2 - sqrt(2))*log(1 - sqrt(2 - sqrt(2))*x + x^2) + (1/16)*sqrt(2 - sqrt(2))*log(1 + sqrt(2 - sqrt(2))*x + x^2) - (1/16)*sqrt(2 + sqrt(2))*log(1 - sqrt(2 + sqrt(2))*x + x^2) + (1/16)*sqrt(2 + sqrt(2))*log(1 + sqrt(2 + sqrt(2))*x + x^2)]
    @test_int [1/(x^8-1), x, 13, -(atan(x)/4) + atan(1 - sqrt(2)*x)/(4*sqrt(2)) - atan(1 + sqrt(2)*x)/(4*sqrt(2)) - atanh(x)/4 + log(1 - sqrt(2)*x + x^2)/(8*sqrt(2)) - log(1 + sqrt(2)*x + x^2)/(8*sqrt(2))]
    @test_int [1/(x^8-x^4+1), x, 19, -(atan((sqrt(2 - sqrt(3)) - 2*x)/sqrt(2 + sqrt(3)))/(2*sqrt(6))) - atan((sqrt(2 + sqrt(3)) - 2*x)/sqrt(2 - sqrt(3)))/(2*sqrt(6)) + atan((sqrt(2 - sqrt(3)) + 2*x)/sqrt(2 + sqrt(3)))/(2*sqrt(6)) + atan((sqrt(2 + sqrt(3)) + 2*x)/sqrt(2 - sqrt(3)))/(2*sqrt(6)) - log(1 - sqrt(2 - sqrt(3))*x + x^2)/(4*sqrt(6)) + log(1 + sqrt(2 - sqrt(3))*x + x^2)/(4*sqrt(6)) - log(1 - sqrt(2 + sqrt(3))*x + x^2)/(4*sqrt(6)) + log(1 + sqrt(2 + sqrt(3))*x + x^2)/(4*sqrt(6))]
    @test_int [x^7/(x^12+1), x, 7, -(atan((1 - 2*x^4)/sqrt(3))/(4*sqrt(3))) - (1/12)*log(1 + x^4) + (1/24)*log(1 - x^4 + x^8)]

    #= ::Section::Closed:: =#
    #=Examples*involving*logarithms=#

    @test_int [log(x), x, 1, -x + x*log(x)]
    @test_int [x*log(x), x, 1, -(x^2/4) + (1/2)*x^2*log(x)]
    @test_int [x^2*log(x), x, 1, -(x^3/9) + (1/3)*x^3*log(x)]
    @test_int [x^p*log(x), x, 1, -(x^(1 + p)/(1 + p)^2) + (x^(1 + p)*log(x))/(1 + p)]
    @test_int [(log(x))^2, x, 2, 2*x - 2*x*log(x) + x*log(x)^2]
    @test_int [x^9*log(x)^11, x, 11, -((6237*x^10)/156250000) + (6237*x^10*log(x))/15625000 - (6237*x^10*log(x)^2)/3125000 + (2079*x^10*log(x)^3)/312500 - (2079*x^10*log(x)^4)/125000 + (2079*x^10*log(x)^5)/62500 - (693*x^10*log(x)^6)/12500 + (99*x^10*log(x)^7)/1250 - (99*x^10*log(x)^8)/1000 + (11/100)*x^10*log(x)^9 - (11/100)*x^10*log(x)^10 + (1/10)*x^10*log(x)^11]
    @test_int [log(x)^2/x, x, 2, log(x)^3/3]
    @test_int [1/log(x), x, 1, LogIntegral(x)]
    @test_int [1/log(x+1), x, 2, LogIntegral(1 + x)]
    @test_int [1/(x*log(x)), x, 2, log(log(x))]
    @test_int [1/(x*log(x))^2, x, 3, -ExpIntegralEi(-log(x)) - 1/(x*log(x))]
    @test_int [(log(x))^p/x, x, 2, log(x)^(1 + p)/(1 + p)]
    @test_int [log(x)*(a*x+b), x, 2, (-b)*x - (a*x^2)/4 + b*x*log(x) + (1/2)*a*x^2*log(x)]
    @test_int [(a*x+b)^2*log(x), x, 4, (-b^2)*x - (1/2)*a*b*x^2 - (a^2*x^3)/9 - (b^3*log(x))/(3*a) + ((b + a*x)^3*log(x))/(3*a)]
    @test_int [log(x)/(a*x+b)^2, x, 2, (x*log(x))/(b*(b + a*x)) - log(b + a*x)/(a*b)]
    @test_int [x*log(a*x+b), x, 3, (b*x)/(2*a) - x^2/4 - (b^2*log(b + a*x))/(2*a^2) + (1/2)*x^2*log(b + a*x)]
    @test_int [x^2*log(a*x+b), x, 3, -((b^2*x)/(3*a^2)) + (b*x^2)/(6*a) - x^3/9 + (b^3*log(b + a*x))/(3*a^3) + (1/3)*x^3*log(b + a*x)]
    @test_int [log(x^2+a^2), x, 3, -2*x + 2*a*atan(x/a) + x*log(a^2 + x^2)]
    @test_int [x*log(x^2+a^2), x, 3, -(x^2/2) + (1/2)*(a^2 + x^2)*log(a^2 + x^2)]
    @test_int [x^2*log(x^2+a^2), x, 4, (2*a^2*x)/3 - (2*x^3)/9 - (2/3)*a^3*atan(x/a) + (1/3)*x^3*log(a^2 + x^2)]
    @test_int [x^4*log(x^2+a^2), x, 4, -((2*a^4*x)/5) + (2*a^2*x^3)/15 - (2*x^5)/25 + (2/5)*a^5*atan(x/a) + (1/5)*x^5*log(a^2 + x^2)]
    @test_int [log(x^2-a^2), x, 3, -2*x + 2*a*atanh(x/a) + x*log(-a^2 + x^2)]
    @test_int [log(log(log(log(x)))), x, 0, CannotIntegrate(log(log(log(log(x)))), x)]

    #= ::Section::Closed:: =#
    #=Examples*involving*circular*functions=#

    @test_int [sin(x), x, 1, -cos(x)]
    @test_int [cos(x), x, 1, sin(x)]
    @test_int [tan(x), x, 1, -log(cos(x))]
    @test_int [1/tan(x), x, 1, log(sin(x))]
    @test_int [1/(1+tan(x))^2, x, 2, (1/2)*log(cos(x) + sin(x)) - 1/(2*(1 + tan(x)))]
    @test_int [1/cos(x), x, 1, atanh(sin(x))]
    @test_int [1/sin(x), x, 1, -atanh(cos(x))]
    @test_int [sin(x)^2, x, 2, x/2 - (1/2)*cos(x)*sin(x)]
    @test_int [x^3*sin(x^2), x, 3, (-(1/2))*x^2*cos(x^2) + sin(x^2)/2]
    @test_int [sin(x)^3, x, 2, -cos(x) + cos(x)^3/3]
    @test_int [sin(x)^p, x, 1, (cos(x)*HypergeometricFunctions._₂F₁(1/2, (1 + p)/2, (3 + p)/2, sin(x)^2)*sin(x)^(1 + p))/((1 + p)*sqrt(cos(x)^2))]
    @test_int [(sin(x)^2+1)^2*cos(x), x, 3, sin(x) + (2*sin(x)^3)/3 + sin(x)^5/5]
    @test_int [cos(x)^2, x, 2, x/2 + (1/2)*cos(x)*sin(x)]
    @test_int [cos(x)^3, x, 2, sin(x) - sin(x)^3/3]
    @test_int [1/cos(x)^2, x, 2, tan(x)]
    @test_int [sin(x)*sin(2*x), x, 1, sin(x)/2 - (1/6)*sin(3*x)]
    @test_int [x*sin(x), x, 2, (-x)*cos(x) + sin(x)]
    @test_int [x^2*sin(x), x, 3, 2*cos(x) - x^2*cos(x) + 2*x*sin(x)]
    @test_int [x*sin(x)^2, x, 2, x^2/4 - (1/2)*x*cos(x)*sin(x) + sin(x)^2/4]
    @test_int [x^2*sin(x)^2, x, 4, -(x/4) + x^3/6 + (1/4)*cos(x)*sin(x) - (1/2)*x^2*cos(x)*sin(x) + (1/2)*x*sin(x)^2]
    @test_int [x*sin(x)^3, x, 3, (-(2/3))*x*cos(x) + (2*sin(x))/3 - (1/3)*x*cos(x)*sin(x)^2 + sin(x)^3/9]
    @test_int [x*cos(x), x, 2, cos(x) + x*sin(x)]
    @test_int [x^2*cos(x), x, 3, 2*x*cos(x) - 2*sin(x) + x^2*sin(x)]
    @test_int [x*cos(x)^2, x, 2, x^2/4 + cos(x)^2/4 + (1/2)*x*cos(x)*sin(x)]
    @test_int [x^2*cos(x)^2, x, 4, -(x/4) + x^3/6 + (1/2)*x*cos(x)^2 - (1/4)*cos(x)*sin(x) + (1/2)*x^2*cos(x)*sin(x)]
    @test_int [x*cos(x)^3, x, 3, (2*cos(x))/3 + cos(x)^3/9 + (2/3)*x*sin(x) + (1/3)*x*cos(x)^2*sin(x)]
    @test_int [sin(x)/x, x, 1, SinIntegral(x)]
    @test_int [cos(x)/x, x, 1, CosIntegral(x)]
    @test_int [sin(x)/x^2, x, 2, CosIntegral(x) - sin(x)/x]
    @test_int [sin(x)^2/x, x, 3, (-(1/2))*CosIntegral(2*x) + log(x)/2]
    @test_int [tan(x)^3, x, 2, log(cos(x)) + tan(x)^2/2]

    @test_int [sin(a+b*x), x, 1, -(cos(a + b*x)/b)]
    @test_int [cos(a+b*x), x, 1, sin(a + b*x)/b]
    @test_int [tan(a+b*x), x, 1, -(log(cos(a + b*x))/b)]
    @test_int [1/tan(a+b*x), x, 1, log(sin(a + b*x))/b]
    @test_int [1/sin(a+b*x), x, 1, -(atanh(cos(a + b*x))/b)]
    @test_int [1/cos(a+b*x), x, 1, atanh(sin(a + b*x))/b]
    @test_int [sin(a+b*x)^2, x, 2, x/2 - (cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [sin(a+b*x)^3, x, 2, -(cos(a + b*x)/b) + cos(a + b*x)^3/(3*b)]
    @test_int [cos(a+b*x)^2, x, 2, x/2 + (cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [cos(a+b*x)^3, x, 2, sin(a + b*x)/b - sin(a + b*x)^3/(3*b)]
    @test_int [1/cos(a+b*x)^2, x, 2, tan(a + b*x)/b]

    @test_int [1/(1+cos(x)), x, 1, sin(x)/(1 + cos(x))]
    @test_int [1/(1-cos(x)), x, 1, -(sin(x)/(1 - cos(x)))]
    @test_int [1/(1+sin(x)), x, 1, -(cos(x)/(1 + sin(x)))]
    @test_int [1/(1-sin(x)), x, 1, cos(x)/(1 - sin(x))]
    @test_int [1/(a+b*sin(x)), x, 3, (2*atan((b + a*tan(x/2))/sqrt(a^2 - b^2)))/sqrt(a^2 - b^2)]
    @test_int [1/(a+b*sin(x)+cos(x)), x, 3, -((2*atanh((b - (1 - a)*tan(x/2))/sqrt(1 - a^2 + b^2)))/sqrt(1 - a^2 + b^2))]
    @test_int [x^2*sin(a+b*x)^2, x, 4, -(x/(4*b^2)) + x^3/6 + (cos(a + b*x)*sin(a + b*x))/(4*b^3) - (x^2*cos(a + b*x)*sin(a + b*x))/(2*b) + (x*sin(a + b*x)^2)/(2*b^2)]
    @test_int [cos(x)*cos(2*x), x, 1, sin(x)/2 + (1/6)*sin(3*x)]
    @test_int [x^2*cos(a+b*x)^2, x, 4, -(x/(4*b^2)) + x^3/6 + (x*cos(a + b*x)^2)/(2*b^2) - (cos(a + b*x)*sin(a + b*x))/(4*b^3) + (x^2*cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [1/tan(x)^3, x, 2, (-(1/2))*cot(x)^2 - log(sin(x))]
    @test_int [x^3*tan(x)^4, x, 17, -(x^2/2) + (4*im*x^3)/3 + x^4/4 - 4*x^2*log(1 + ℯ^(2*im*x)) + log(cos(x)) + 4*im*x*PolyLog(2, -ℯ^(2*im*x)) - 2*PolyLog(3, -ℯ^(2*im*x)) + x*tan(x) - x^3*tan(x) - (1/2)*x^2*tan(x)^2 + (1/3)*x^3*tan(x)^3]
    @test_int [x^3*tan(x)^6, x, 34, (19*x^2)/20 - (23*im*x^3)/15 - x^4/4 + (23/5)*x^2*log(1 + ℯ^(2*im*x)) - 2*log(cos(x)) - (23/5)*im*x*PolyLog(2, -ℯ^(2*im*x)) + (23/10)*PolyLog(3, -ℯ^(2*im*x)) - (19/10)*x*tan(x) + x^3*tan(x) - tan(x)^2/20 + (4/5)*x^2*tan(x)^2 + (1/10)*x*tan(x)^3 - (1/3)*x^3*tan(x)^3 - (3/20)*x^2*tan(x)^4 + (1/5)*x^3*tan(x)^5]
    @test_int [x*tan(x)^2, x, 3, -(x^2/2) + log(cos(x)) + x*tan(x)]
    @test_int [sin(2*x)*cos(3*x), x, 1, cos(x)/2 - (1/10)*cos(5*x)]
    @test_int [sin(x)^2*cos(x)^2, x, 3, x/8 + (1/8)*cos(x)*sin(x) - (1/4)*cos(x)^3*sin(x)]
    @test_int [1/(sin(x)^2*cos(x)^2), x, 3, -cot(x) + tan(x)]

    @test_int [d^x*sin(x), x, 1, -((d^x*cos(x))/(1 + log(d)^2)) + (d^x*log(d)*sin(x))/(1 + log(d)^2)]
    @test_int [d^x*cos(x), x, 1, (d^x*cos(x)*log(d))/(1 + log(d)^2) + (d^x*sin(x))/(1 + log(d)^2)]
    @test_int [x*d^x*sin(x), x, 4, (2*d^x*cos(x)*log(d))/(1 + log(d)^2)^2 - (d^x*x*cos(x))/(1 + log(d)^2) + (d^x*sin(x))/(1 + log(d)^2)^2 - (d^x*log(d)^2*sin(x))/(1 + log(d)^2)^2 + (d^x*x*log(d)*sin(x))/(1 + log(d)^2)]
    @test_int [x*d^x*cos(x), x, 4, (d^x*cos(x))/(1 + log(d)^2)^2 - (d^x*cos(x)*log(d)^2)/(1 + log(d)^2)^2 + (d^x*x*cos(x)*log(d))/(1 + log(d)^2) - (2*d^x*log(d)*sin(x))/(1 + log(d)^2)^2 + (d^x*x*sin(x))/(1 + log(d)^2)]
    @test_int [x^2*d^x*sin(x), x, 11, (2*d^x*cos(x))/(1 + log(d)^2)^3 - (6*d^x*cos(x)*log(d)^2)/(1 + log(d)^2)^3 + (4*d^x*x*cos(x)*log(d))/(1 + log(d)^2)^2 - (d^x*x^2*cos(x))/(1 + log(d)^2) - (6*d^x*log(d)*sin(x))/(1 + log(d)^2)^3 + (2*d^x*log(d)^3*sin(x))/(1 + log(d)^2)^3 + (2*d^x*x*sin(x))/(1 + log(d)^2)^2 - (2*d^x*x*log(d)^2*sin(x))/(1 + log(d)^2)^2 + (d^x*x^2*log(d)*sin(x))/(1 + log(d)^2)]
    @test_int [x^2*d^x*cos(x), x, 11, -((6*d^x*cos(x)*log(d))/(1 + log(d)^2)^3) + (2*d^x*cos(x)*log(d)^3)/(1 + log(d)^2)^3 + (2*d^x*x*cos(x))/(1 + log(d)^2)^2 - (2*d^x*x*cos(x)*log(d)^2)/(1 + log(d)^2)^2 + (d^x*x^2*cos(x)*log(d))/(1 + log(d)^2) - (2*d^x*sin(x))/(1 + log(d)^2)^3 + (6*d^x*log(d)^2*sin(x))/(1 + log(d)^2)^3 - (4*d^x*x*log(d)*sin(x))/(1 + log(d)^2)^2 + (d^x*x^2*sin(x))/(1 + log(d)^2)]
    @test_int [x^3*d^x*sin(x), x, 25, -((24*d^x*cos(x)*log(d))/(1 + log(d)^2)^4) + (24*d^x*cos(x)*log(d)^3)/(1 + log(d)^2)^4 + (6*d^x*x*cos(x))/(1 + log(d)^2)^3 - (18*d^x*x*cos(x)*log(d)^2)/(1 + log(d)^2)^3 + (6*d^x*x^2*cos(x)*log(d))/(1 + log(d)^2)^2 - (d^x*x^3*cos(x))/(1 + log(d)^2) - (6*d^x*sin(x))/(1 + log(d)^2)^4 + (36*d^x*log(d)^2*sin(x))/(1 + log(d)^2)^4 - (6*d^x*log(d)^4*sin(x))/(1 + log(d)^2)^4 - (18*d^x*x*log(d)*sin(x))/(1 + log(d)^2)^3 + (6*d^x*x*log(d)^3*sin(x))/(1 + log(d)^2)^3 + (3*d^x*x^2*sin(x))/(1 + log(d)^2)^2 - (3*d^x*x^2*log(d)^2*sin(x))/(1 + log(d)^2)^2 + (d^x*x^3*log(d)*sin(x))/(1 + log(d)^2)]
    @test_int [x^3*d^x*cos(x), x, 25, -((6*d^x*cos(x))/(1 + log(d)^2)^4) + (36*d^x*cos(x)*log(d)^2)/(1 + log(d)^2)^4 - (6*d^x*cos(x)*log(d)^4)/(1 + log(d)^2)^4 - (18*d^x*x*cos(x)*log(d))/(1 + log(d)^2)^3 + (6*d^x*x*cos(x)*log(d)^3)/(1 + log(d)^2)^3 + (3*d^x*x^2*cos(x))/(1 + log(d)^2)^2 - (3*d^x*x^2*cos(x)*log(d)^2)/(1 + log(d)^2)^2 + (d^x*x^3*cos(x)*log(d))/(1 + log(d)^2) + (24*d^x*log(d)*sin(x))/(1 + log(d)^2)^4 - (24*d^x*log(d)^3*sin(x))/(1 + log(d)^2)^4 - (6*d^x*x*sin(x))/(1 + log(d)^2)^3 + (18*d^x*x*log(d)^2*sin(x))/(1 + log(d)^2)^3 - (6*d^x*x^2*log(d)*sin(x))/(1 + log(d)^2)^2 + (d^x*x^3*sin(x))/(1 + log(d)^2)]

    @test_int [sin(x)*sin(2*x)*sin(3*x), x, 5, (-(1/8))*cos(2*x) - (1/16)*cos(4*x) + (1/24)*cos(6*x)]
    @test_int [cos(x)*cos(2*x)*cos(3*x), x, 5, x/4 + (1/8)*sin(2*x) + (1/16)*sin(4*x) + (1/24)*sin(6*x)]
    @test_int [sin(k*x)^3*x^2, x, 6, (14*cos(k*x))/(9*k^3) - (2*x^2*cos(k*x))/(3*k) - (2*cos(k*x)^3)/(27*k^3) + (4*x*sin(k*x))/(3*k^2) - (x^2*cos(k*x)*sin(k*x)^2)/(3*k) + (2*x*sin(k*x)^3)/(9*k^2)]
    @test_int [x*cos(k/sin(x))*cos(x)/sin(x)^2, x, 0, CannotIntegrate(x*cos(k*csc(x))*cot(x)*csc(x), x)]

    #= Mixed*angles*and*half*angles. =#
    @test_int [cos(x)/(sin(x)*tan(x/2)), x, 4, -x - cot(x/2)]
    @test_int [sin(a*x)/(b+c*sin(a*x))^2, x, 5, -((2*c*atan((c + b*tan((a*x)/2))/sqrt(b^2 - c^2)))/(a*(b^2 - c^2)^(3/2))) - (b*cos(a*x))/(a*(b^2 - c^2)*(b + c*sin(a*x)))]

    #= ::Section::Closed:: =#
    #=Examples*involving*trig*functions*of*logarithms=#

    @test_int [sin(log(x)), x, 1, (-(1/2))*x*cos(log(x)) + (1/2)*x*sin(log(x))]
    @test_int [cos(log(x)), x, 1, (1/2)*x*cos(log(x)) + (1/2)*x*sin(log(x))]

    #= ::Section::Closed:: =#
    #=Examples*involving*exponentials=#

    @test_int [ℯ^x, x, 1, ℯ^x]
    @test_int [a^x, x, 1, a^x/log(a)]
    @test_int [ℯ^(a*x), x, 1, ℯ^(a*x)/a]
    @test_int [ℯ^(a*x)/x, x, 1, ExpIntegralEi(a*x)]
    @test_int [1/(a+b*ℯ^(m*x)), x, 4, x/a - log(a + b*ℯ^(m*x))/(a*m)]
    @test_int [ℯ^(2*x)/(1+ℯ^x), x, 3, ℯ^x - log(1 + ℯ^x)]
    @test_int [ℯ^(2*x)*ℯ^(a*x), x, 2, ℯ^((2 + a)*x)/(2 + a)]
    @test_int [1/(a*ℯ^(m*x)+b*ℯ^(-m*x)), x, 2, atan((sqrt(a)*ℯ^(m*x))/sqrt(b))/(sqrt(a)*sqrt(b)*m)]
    @test_int [x*ℯ^(a*x), x, 2, -(ℯ^(a*x)/a^2) + (ℯ^(a*x)*x)/a]
    @test_int [x^20*ℯ^x, x, 21, 2432902008176640000*ℯ^x - 2432902008176640000*ℯ^x*x + 1216451004088320000*ℯ^x*x^2 - 405483668029440000*ℯ^x*x^3 + 101370917007360000*ℯ^x*x^4 - 20274183401472000*ℯ^x*x^5 + 3379030566912000*ℯ^x*x^6 - 482718652416000*ℯ^x*x^7 + 60339831552000*ℯ^x*x^8 - 6704425728000*ℯ^x*x^9 + 670442572800*ℯ^x*x^10 - 60949324800*ℯ^x*x^11 + 5079110400*ℯ^x*x^12 - 390700800*ℯ^x*x^13 + 27907200*ℯ^x*x^14 - 1860480*ℯ^x*x^15 + 116280*ℯ^x*x^16 - 6840*ℯ^x*x^17 + 380*ℯ^x*x^18 - 20*ℯ^x*x^19 + ℯ^x*x^20]
    @test_int [a^x/b^x, x, 2, a^x/(b^x*(log(a) - log(b)))]
    @test_int [a^x*b^x, x, 2, (a^x*b^x)/(log(a) + log(b))]
    @test_int [a^x/x^2, x, 2, -(a^x/x) + ExpIntegralEi(x*log(a))*log(a)]
    @test_int [x*a^x/(1+b*x)^2, x, 5, a^x/(b^2*(1 + b*x)) + ExpIntegralEi(((1 + b*x)*log(a))/b)/(a^b^(-1)*b^2) - (ExpIntegralEi(((1 + b*x)*log(a))/b)*log(a))/(a^b^(-1)*b^3)]
    @test_int [x*ℯ^(a*x)/(1+a*x)^2, x, 1, ℯ^(a*x)/(a^2*(1 + a*x))]
    @test_int [x*k^(x^2), x, 1, k^x^2/(2*log(k))]
    @test_int [ℯ^(x^2), x, 1, (1/2)*sqrt(pi)*erfi(x)]
    @test_int [x*ℯ^(x^2), x, 1, ℯ^x^2/2]
    @test_int [(x+1)*ℯ^(1/x)/x^4, x, 7, -ℯ^(1/x) - ℯ^(1/x)/x^2 + ℯ^(1/x)/x]
    @test_int [(2*x^3+x)*(ℯ^(x^2))^2*ℯ^(1-x*ℯ^(x^2))/(1-x*ℯ^(x^2))^2, x, -3, -(ℯ^(1 - ℯ^x^2*x)/(-1 + ℯ^x^2*x))]
    @test_int [ℯ^(ℯ^(ℯ^(ℯ^x))), x, 1, CannotIntegrate(ℯ^ℯ^ℯ^ℯ^x, x)]

    #= Examples*involving*exponentials*and*logarithms. =#
    @test_int [ℯ^x*log(x), x, 2, -ExpIntegralEi(x) + ℯ^x*log(x)]
    @test_int [x*ℯ^x*log(x), x, 5, -ℯ^x + ExpIntegralEi(x) - ℯ^x*log(x) + ℯ^x*x*log(x)]
    @test_int [ℯ^(2*x)*log(ℯ^x), x, 3, -(ℯ^(2*x)/4) + (1/2)*ℯ^(2*x)*log(ℯ^x)]

    #= ::Section::Closed:: =#
    #=Examples*involving*square*roots=#

    @test_int [sqrt(2)*x^2 + 2*x, x, 1, x^2 + (sqrt(2)*x^3)/3]
    @test_int [log(x)/sqrt(a*x+b), x, 4, -((4*sqrt(b + a*x))/a) + (4*sqrt(b)*atanh(sqrt(b + a*x)/sqrt(b)))/a + (2*sqrt(b + a*x)*log(x))/a]
    @test_int [sqrt(a+b*x)*sqrt(c+d*x), x, 5, ((b*c - a*d)*sqrt(a + b*x)*sqrt(c + d*x))/(4*b*d) + ((a + b*x)^(3/2)*sqrt(c + d*x))/(2*b) - ((b*c - a*d)^2*atanh((sqrt(d)*sqrt(a + b*x))/(sqrt(b)*sqrt(c + d*x))))/(4*b^(3/2)*d^(3/2))]
    @test_int [sqrt(a+b*x), x, 1, (2*(a + b*x)^(3/2))/(3*b)]
    @test_int [x*sqrt(a+b*x), x, 2, -((2*a*(a + b*x)^(3/2))/(3*b^2)) + (2*(a + b*x)^(5/2))/(5*b^2)]
    @test_int [x^2*sqrt(a+b*x), x, 2, (2*a^2*(a + b*x)^(3/2))/(3*b^3) - (4*a*(a + b*x)^(5/2))/(5*b^3) + (2*(a + b*x)^(7/2))/(7*b^3)]
    @test_int [sqrt(a+b*x)/x, x, 3, 2*sqrt(a + b*x) - 2*sqrt(a)*atanh(sqrt(a + b*x)/sqrt(a))]
    @test_int [sqrt(a+b*x)/x^2, x, 3, -(sqrt(a + b*x)/x) - (b*atanh(sqrt(a + b*x)/sqrt(a)))/sqrt(a)]
    @test_int [1/sqrt(a+b*x), x, 1, (2*sqrt(a + b*x))/b]
    @test_int [x/sqrt(a+b*x), x, 2, -((2*a*sqrt(a + b*x))/b^2) + (2*(a + b*x)^(3/2))/(3*b^2)]
    @test_int [x^2/sqrt(a+b*x), x, 2, (2*a^2*sqrt(a + b*x))/b^3 - (4*a*(a + b*x)^(3/2))/(3*b^3) + (2*(a + b*x)^(5/2))/(5*b^3)]
    @test_int [1/(x*sqrt(a+b*x)), x, 2, -((2*atanh(sqrt(a + b*x)/sqrt(a)))/sqrt(a))]
    @test_int [1/(x^2*sqrt(a+b*x)), x, 3, -(sqrt(a + b*x)/(a*x)) + (b*atanh(sqrt(a + b*x)/sqrt(a)))/a^(3/2)]
    @test_int [sqrt(a+b*x)^p, x, 1, (2*(a + b*x)^((2 + p)/2))/(b*(2 + p))]
    @test_int [x*sqrt(a+b*x)^p, x, 2, -((2*a*(a + b*x)^((2 + p)/2))/(b^2*(2 + p))) + (2*(a + b*x)^((4 + p)/2))/(b^2*(4 + p))]
    @test_int [atan((-sqrt(2)+2*x)/sqrt(2)), x, 6, atan(1 - sqrt(2)*x)/sqrt(2) - x*atan(1 - sqrt(2)*x) - log(1 - sqrt(2)*x + x^2)/(2*sqrt(2))]
    @test_int [1/sqrt(x^2-1), x, 2, atanh(x/sqrt(-1 + x^2))]
    @test_int [sqrt(x+1)*sqrt(x), x, 4, (1/4)*sqrt(x)*sqrt(1 + x) + (1/2)*x^(3/2)*sqrt(1 + x) - asinh(sqrt(x))/4]

    @test_int [sin(sqrt(x)), x, 3, -2*sqrt(x)*cos(sqrt(x)) + 2*sin(sqrt(x))]
    @test_int [x*sqrt(1 - x^2)^(-9/4), x, 1, 4/(1 - x^2)^(1/8)]
    @test_int [x/sqrt(1 - x^4), x, 2, asin(x^2)/2]
    @test_int [1/(x*sqrt(1 + x^4)), x, 3, (-(1/2))*atanh(sqrt(1 + x^4))]
    @test_int [x/sqrt(1 + x^2 + x^4), x, 3, (1/2)*asinh((1 + 2*x^2)/sqrt(3))]
    @test_int [1/(x*sqrt(x^2 - 1 - x^4)), x, 3, (-(1/2))*atan((2 - x^2)/(2*sqrt(-1 + x^2 - x^4)))]

    #= Examples*generated*by*differentiating*various*functions. =#
    @test_int [(1 + x)/((1 - x)^2*sqrt(1 + x^2)), x, 1, sqrt(1 + x^2)/(1 - x)]
    @test_int [1/sqrt(1 + x^2), x, 1, asinh(x)]
    @test_int [(sqrt(x)*sqrt(1 + x) + sqrt(x)*sqrt(2 + x) + sqrt(1 + x)*sqrt(2 + x))/(2*sqrt(x)*sqrt(1 + x)*sqrt(2 + x)), x, 3, sqrt(x) + sqrt(1 + x) + sqrt(2 + x)]
    @test_int [(-2*sqrt(1 + x^3) + 5*x^4*sqrt(1 + x^3) - 3*x^2*sqrt(1 - 2*x + x^5))/(2*sqrt(1 + x^3)*sqrt(1 - 2*x + x^5)), x, 5, -sqrt(1 + x^3) + sqrt(1 - 2*x + x^5)]

    #= ::Section::Closed:: =#
    #=Examples*from*James*Davenport's*thesis=#

    @test_int [1/sqrt(x^2-1)+10/sqrt(x^2-4), x, 5, 10*atanh(x/sqrt(-4 + x^2)) + atanh(x/sqrt(-1 + x^2))]
    @test_int [sqrt(x+sqrt(x^2+a^2))/x, x, 6, 2*sqrt(x + sqrt(a^2 + x^2)) - 2*sqrt(a)*atan(sqrt(x + sqrt(a^2 + x^2))/sqrt(a)) - 2*sqrt(a)*atanh(sqrt(x + sqrt(a^2 + x^2))/sqrt(a))]

    #= Another*such*example*from*James*Davenport's*thesis*(p. 146). =#
    #= It*contains*a*point*of*order*3, which*is*found*by*use*of*Mazur's =#
    #= bound*on*the*torsion*of*elliptic*curves*over*the*rationals; =#
    @test_int [(3*x^2)/(2*(1 + x^3 + sqrt(1 + x^3))), x, 4, log(1 + sqrt(1 + x^3))]

    #= ::Section::Closed:: =#
    #=Examples*quoted*by*Joel*Moses=#

    @test_int [1/sqrt(2*h*r^2-alpha^2), r, 2, atanh((sqrt(2)*sqrt(h)*r)/sqrt(-alpha^2 + 2*h*r^2))/(sqrt(2)*sqrt(h))]
    @test_int [1/(r*sqrt(2*h*r^2-alpha^2-epsilon^2)), r, 3, atan(sqrt(-alpha^2 - epsilon^2 + 2*h*r^2)/sqrt(alpha^2 + epsilon^2))/sqrt(alpha^2 + epsilon^2)]
    @test_int [1/(r*sqrt(2*h*r^2-alpha^2-2*k*r)), r, 2, -(atan((alpha^2 + k*r)/(alpha*sqrt(-alpha^2 - 2*k*r + 2*h*r^2)))/alpha)]
    @test_int [1/(r*sqrt(2*h*r^2-alpha^2-epsilon^2-2*k*r)), r, 2, -(atan((alpha^2 + epsilon^2 + k*r)/(sqrt(alpha^2 + epsilon^2)*sqrt(-alpha^2 - epsilon^2 - 2*k*r + 2*h*r^2)))/sqrt(alpha^2 + epsilon^2))]
    @test_int [r/sqrt(2*e*r^2-alpha^2), r, 1, sqrt(-alpha^2 + 2*e*r^2)/(2*e)]
    @test_int [r/sqrt(2*e*r^2-alpha^2-epsilon^2), r, 1, sqrt(-alpha^2 - epsilon^2 + 2*e*r^2)/(2*e)]
    @test_int [r/sqrt(2*e*r^2-alpha^2-2*k*r^4), r, 3, -(atan((e - 2*k*r^2)/(sqrt(2)*sqrt(k)*sqrt(-alpha^2 + 2*e*r^2 - 2*k*r^4)))/(2*sqrt(2)*sqrt(k)))]
    @test_int [r/sqrt(2*e*r^2-alpha^2-2*k*r), r, 3, sqrt(-alpha^2 - 2*k*r + 2*e*r^2)/(2*e) - (k*atanh((k - 2*e*r)/(sqrt(2)*sqrt(e)*sqrt(-alpha^2 - 2*k*r + 2*e*r^2))))/(2*sqrt(2)*e^(3/2))]
    @test_int [1/(r*sqrt(2*h*r^2-alpha^2-2*k*r^4)), r, 3, -(atan((alpha^2 - h*r^2)/(alpha*sqrt(-alpha^2 + 2*h*r^2 - 2*k*r^4)))/(2*alpha))]
    @test_int [1/(r*sqrt(2*h*r^2-alpha^2-epsilon^2-2*k*r^4)), r, 3, -(atan((alpha^2 + epsilon^2 - h*r^2)/(sqrt(alpha^2 + epsilon^2)*sqrt(-alpha^2 - epsilon^2 + 2*h*r^2 - 2*k*r^4)))/(2*sqrt(alpha^2 + epsilon^2)))]

    #= ::Section::Closed:: =#
    #=Examples*from*Novosibirsk=#

    #= Many*of*these*integrals*used*to*require*Steve*Harrington's*code*to*evaluate. =#
    #= They*originated*in*Novosibirsk*as*examples*of*using*Analytik. =#
    #= There*are*still*a*few*examples*that*could*be*evaluated*using*better*heuristics. =#

    @test_int [a*sin(3*x+5)^2*cos(3*x+5), x, 3, (1/9)*a*sin(5 + 3*x)^3]
    @test_int [log(x^2)/x^3, x, 1, -(1/(2*x^2)) - log(x^2)/(2*x^2)]
    @test_int [x*sin(x+a), x, 2, (-x)*cos(a + x) + sin(a + x)]
    @test_int [(log(x)*(1-x)-1)/(ℯ^x*log(x)^2), x, 1, x/(ℯ^x*log(x))]
    @test_int [x^3/(a*x^2+b), x, 3, x^2/(2*a) - (b*log(b + a*x^2))/(2*a^2)]
    @test_int [x^(1/2)*(x+1)^(-7/2), x, 2, (2*x^(3/2))/(5*(1 + x)^(5/2)) + (4*x^(3/2))/(15*(1 + x)^(3/2))]
    @test_int [x^(-1)*(x+1)^(-1), x, 3, log(x) - log(1 + x)]
    @test_int [x^(-1/2)*(2*x-1)^(-1), x, 2, (-sqrt(2))*atanh(sqrt(2)*sqrt(x))]
    @test_int [(x^2+1)*x^(1/2), x, 2, (2*x^(3/2))/3 + (2*x^(7/2))/7]
    @test_int [x^(-1)*(x-a)^(1/3), x, 5, 3*(-a + x)^(1/3) + sqrt(3)*a^(1/3)*atan((a^(1/3) - 2*(-a + x)^(1/3))/(sqrt(3)*a^(1/3))) + (1/2)*a^(1/3)*log(x) - (3/2)*a^(1/3)*log(a^(1/3) + (-a + x)^(1/3))]
    @test_int [x*sinh(x), x, 2, x*cosh(x) - sinh(x)]
    @test_int [x*cosh(x), x, 2, -cosh(x) + x*sinh(x)]
    @test_int [sinh(2*x)/cosh(2*x), x, 1, (1/2)*log(cosh(2*x))]
    @test_int [(im*eps*sinh(x)-1)/(eps*im*cosh(x)+im*a-x), x, 1, log(a + im*x + eps*cosh(x))]
    @test_int [sin(2*x+3)*cos(x)^2, x, 4, (-(1/4))*cos(3 + 2*x) - (1/16)*cos(3 + 4*x) + (1/4)*x*sin(3)]
    @test_int [x*atan(x), x, 3, -(x/2) + atan(x)/2 + (1/2)*x^2*atan(x)]
    @test_int [x*acot(x), x, 3, x/2 + (1/2)*x^2*acot(x) - atan(x)/2]
    @test_int [x*log(x^2+a), x, 3, -(x^2/2) + (1/2)*(a + x^2)*log(a + x^2)]
    @test_int [sin(x+a)*cos(x), x, 3, (-(1/4))*cos(a + 2*x) + (1/2)*x*sin(a)]
    @test_int [cos(x+a)*sin(x), x, 3, (-(1/4))*cos(a + 2*x) - (1/2)*x*sin(a)]
    @test_int [(1+sin(x))^(1/2), x, 1, -((2*cos(x))/sqrt(1 + sin(x)))]
    @test_int [(1-sin(x))^(1/2), x, 1, (2*cos(x))/sqrt(1 - sin(x))]
    @test_int [(1+cos(x))^(1/2), x, 1, (2*sin(x))/sqrt(1 + cos(x))]
    @test_int [(1-cos(x))^(1/2), x, 1, -((2*sin(x))/sqrt(1 - cos(x)))]
    @test_int [1/(x^(1/2)-(x-1)^(1/2)), x, 3, (2/3)*(-1 + x)^(3/2) + (2*x^(3/2))/3]
    @test_int [1/(1-(x+1)^(1/2)), x, 4, -2*sqrt(1 + x) - 2*log(1 - sqrt(1 + x))]
    @test_int [x/(x^4+36)^(1/2), x, 2, (1/2)*asinh(x^2/6)]
    @test_int [1/(x^(1/3)+x^(1/2)), x, 4, 6*x^(1/6) - 3*x^(1/3) + 2*sqrt(x) - 6*log(1 + x^(1/6))]
    @test_int [log(2+3*x^2), x, 3, -2*x + 2*sqrt(2/3)*atan(sqrt(3/2)*x) + x*log(2 + 3*x^2)]
    @test_int [cot(x), x, 1, log(sin(x))]
    @test_int [cot(x)^4, x, 3, x + cot(x) - cot(x)^3/3]
    @test_int [tanh(x), x, 1, log(cosh(x))]
    @test_int [coth(x), x, 1, log(sinh(x))]
    @test_int [b^x, x, 1, b^x/log(b)]
    @test_int [(x^4+x^(-4)+2)^(1/2), x, 4, -((x*sqrt(2 + 1/x^4 + x^4))/(1 + x^4)) + (x^5*sqrt(2 + 1/x^4 + x^4))/(3*(1 + x^4))]
    @test_int [(2*x+1)/(3*x+2), x, 2, (2*x)/3 - (1/9)*log(2 + 3*x)]
    @test_int [x*log(x+(x^2+1)^(1/2)), x, 3, (-(1/4))*x*sqrt(1 + x^2) + asinh(x)/4 + (1/2)*x^2*log(x + sqrt(1 + x^2))]
    @test_int [x*(ℯ^x*sin(x)+1)^2, x, 14, -((3*ℯ^(2*x))/32) + (1/8)*ℯ^(2*x)*x + x^2/2 + ℯ^x*cos(x) - ℯ^x*x*cos(x) - (1/32)*ℯ^(2*x)*cos(2*x) + ℯ^x*x*sin(x) + (1/16)*ℯ^(2*x)*cos(x)*sin(x) - (1/4)*ℯ^(2*x)*x*cos(x)*sin(x) - (1/16)*ℯ^(2*x)*sin(x)^2 + (1/4)*ℯ^(2*x)*x*sin(x)^2 + (1/32)*ℯ^(2*x)*sin(2*x)]
    @test_int [x*ℯ^x*cos(x), x, 4, (1/2)*ℯ^x*x*cos(x) - (1/2)*ℯ^x*sin(x) + (1/2)*ℯ^x*x*sin(x)]

    #= ::Section::Closed:: =#
    #=Examples*from*Herbert*Stoyan=#

    @test_int [1/(x-3)^4, x, 1, 1/(3*(3 - x)^3)]
    @test_int [x/(x^3-1), x, 6, atan((1 + 2*x)/sqrt(3))/sqrt(3) + (1/3)*log(1 - x) - (1/6)*log(1 + x + x^2)]
    @test_int [x/(x^4-1), x, 2, (-(1/2))*atanh(x^2)]
    @test_int [log(x)*(x^3+1)/(x^4+2), x, 10, (1/8)*(2 + im*(-2)^(1/4))*log(x)*log(1 - ((1 + im)*x)/2^(3/4)) + (1/16)*(4 + (1 - im)*2^(3/4))*log(x)*log(1 + ((1 + im)*x)/2^(3/4)) + (1/8)*(2 + (-2)^(1/4))*log(x)*log(1 - ((-1)^(3/4)*x)/2^(1/4)) + (1/8)*(2 - (-2)^(1/4))*log(x)*log(1 + ((-1)^(3/4)*x)/2^(1/4)) + (1/16)*(4 + (1 - im)*2^(3/4))*PolyLog(2, -(((1 + im)*x)/2^(3/4))) + (1/8)*(2 + im*(-2)^(1/4))*PolyLog(2, ((1 + im)*x)/2^(3/4)) + (1/8)*(2 - (-2)^(1/4))*PolyLog(2, -(((-1)^(3/4)*x)/2^(1/4))) + (1/8)*(2 + (-2)^(1/4))*PolyLog(2, ((-1)^(3/4)*x)/2^(1/4))]
    @test_int [log(x)+log(x+1)+log(x+2), x, 6, -3*x + x*log(x) + (1 + x)*log(1 + x) + (2 + x)*log(2 + x)]
    @test_int [1/(x^3+5), x, 6, -(atan((5^(1/3) - 2*x)/(sqrt(3)*5^(1/3)))/(sqrt(3)*5^(2/3))) + log(5^(1/3) + x)/(3*5^(2/3)) - log(5^(2/3) - 5^(1/3)*x + x^2)/(6*5^(2/3))]
    @test_int [1/sqrt(1+x^2), x, 1, asinh(x)]
    @test_int [sqrt(x^2+3), x, 2, (1/2)*x*sqrt(3 + x^2) + (3/2)*asinh(x/sqrt(3))]
    @test_int [x/(x+1)^2, x, 2, 1/(1 + x) + log(1 + x)]

    #= ::Section::Closed:: =#
    #=Examples*from*Moses' SIN*program=#

    @test_int [asin(x), x, 2, sqrt(1 - x^2) + x*asin(x)]
    @test_int [x^2*asin(x), x, 4, sqrt(1 - x^2)/3 - (1/9)*(1 - x^2)^(3/2) + (1/3)*x^3*asin(x)]
    @test_int [sec(x)^2/(1+sec(x)^2-3*tan(x)), x, 4, -log(cos(x) - sin(x)) + log(2*cos(x) - sin(x))]
    @test_int [1/sec(x)^2, x, 2, x/2 + (1/2)*cos(x)*sin(x)]
    @test_int [(5*x^2-3*x-2)/(x^2*(x-2)), x, 2, -(1/x) + 3*log(2 - x) + 2*log(x)]
    @test_int [1/(4*x^2+9)^(1/2), x, 1, (1/2)*asinh((2*x)/3)]
    @test_int [(x^2+4)^(-1/2), x, 1, asinh(x/2)]
    @test_int [1/(9*x^2-12*x+10), x, 2, -(atan((2 - 3*x)/sqrt(6))/(3*sqrt(6)))]
    @test_int [1/(x^8-2*x^7+2*x^6-2*x^5+x^4), x, 3, 1/(2*(1 - x)) - 1/(3*x^3) - 1/x^2 - 2/x - (5/2)*log(1 - x) + 2*log(x) + (1/4)*log(1 + x^2)]
    @test_int [(a*x^3+b*x^2+c*x+d)/((x+1)*x*(x-3)), x, 2, a*x + (1/12)*(27*a + 9*b + 3*c + d)*log(3 - x) - (1/3)*d*log(x) - (1/4)*(a - b + c - d)*log(1 + x)]
    @test_int [1/(2-log(x^2+1))^5, x, 0, Unintegrable(1/(2 - log(1 + x^2))^5, x)]

    #= ::Section::Closed:: =#
    #=Miscellaneous*examples=#

    #= The*next*integral*appeared*in*Risch's*1968*paper. =#

    @test_int [2*x*ℯ^(x^2)*log(x)+ℯ^(x^2)/x+(log(x)-2)/(log(x)^2+x)^2+((2/x)*log(x)+(1/x)+1)/(log(x)^2+x), x, 9, ℯ^x^2*log(x) - log(x)/(x + log(x)^2) + log(x + log(x)^2)]

    #= The*following*integral*would*not*evaluate*in*REDUCE*3.3.0. =#

    @test_int [exp(x*z+x/2)*sin(pi*z)^4*x^4, z, 4, (24*ℯ^(x/2 + x*z)*pi^4*x^3)/(64*pi^4 + 20*pi^2*x^2 + x^4) - (24*ℯ^(x/2 + x*z)*pi^3*x^4*cos(pi*z)*sin(pi*z))/(64*pi^4 + 20*pi^2*x^2 + x^4) + (12*ℯ^(x/2 + x*z)*pi^2*x^5*sin(pi*z)^2)/(64*pi^4 + 20*pi^2*x^2 + x^4) - (4*ℯ^(x/2 + x*z)*pi*x^4*cos(pi*z)*sin(pi*z)^3)/(16*pi^2 + x^2) + (ℯ^(x/2 + x*z)*x^5*sin(pi*z)^4)/(16*pi^2 + x^2)]

    #= Examples*involving*the*error*function. =#

    @test_int [erf(x), x, 1, 1/(ℯ^x^2*sqrt(pi)) + x*erf(x)]
    @test_int [erf(x+a), x, 1, 1/(ℯ^(a + x)^2*sqrt(pi)) + (a + x)*erf(a + x)]

    #= Some*interesting*integrals*of*algebraic*functions; =#
    #= The*Chebyshev*integral. =#

    @test_int [(2*x^6+4*x^5+7*x^4-3*x^3-x*x-8*x-8)/((2*x^2-1)^2*sqrt(x^4+4*x^3+2*x^2+1)), x, -10, ((1 + 2*x)*sqrt(1 + 2*x^2 + 4*x^3 + x^4))/(2*(-1 + 2*x^2)) - atanh((x*(2 + x)*(7 - x + 27*x^2 + 33*x^3))/((2 + 37*x^2 + 31*x^3)*sqrt(1 + 2*x^2 + 4*x^3 + x^4)))]

    #= This*integral*came*from*Dr. G.S. Joyce*of*Imperial*College*London. =#

    @test_int [(1+2*y)*sqrt(1-5*y-5*y^2)/(y*(1+y)*(2+y)*sqrt(1-y-y^2)), y, -2, (-(1/4))*atanh(((1 - 3*y)*sqrt(1 - 5*y - 5*y^2))/((1 - 5*y)*sqrt(1 - y - y^2))) - (1/2)*atanh(((4 + 3*y)*sqrt(1 - 5*y - 5*y^2))/((6 + 5*y)*sqrt(1 - y - y^2))) + (9/4)*atanh(((11 + 7*y)*sqrt(1 - 5*y - 5*y^2))/(3*(7 + 5*y)*sqrt(1 - y - y^2)))]

    #= This*one*has*a*simple*result. =#

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

    #= This*used*to*reveal*bugs*in*the*integrator*which*have*been*fixed. =#

    @test_int [sqrt(-4*sqrt(2) + 9)*x - sqrt(x^4 + 2*x^2 + 4*x + 1)*sqrt(2), x, -1, (1/2)*sqrt(9 - 4*sqrt(2))*x^2 - sqrt(2)*((-(1/3))*sqrt(1 + 4*x + 2*x^2 + x^4) + (1/3)*(1 + x)*sqrt(1 + 4*x + 2*x^2 + x^4) + (4*im*(-13 + 3*sqrt(33))^(1/3)*sqrt(1 + 4*x + 2*x^2 + x^4))/(4*2^(2/3)*(-im + sqrt(3)) - 2*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3) + 6*im*(-13 + 3*sqrt(33))^(1/3)*x) - (8*2^(2/3)*sqrt(3/(-13 + 3*sqrt(33) + 4*(-26 + 6*sqrt(33))^(1/3)))*sqrt((im*(-19899 + 3445*sqrt(33) + (-26 + 6*sqrt(33))^(2/3)*(-2574 + 466*sqrt(33)) + (-26 + 6*sqrt(33))^(1/3)*(-19899 + 3445*sqrt(33)) + (59697 - 10335*sqrt(33))*x))/((-39 - 13*im*sqrt(3) + 9*im*sqrt(11) + 9*sqrt(33) + 4*im*(3*im + sqrt(3))*(-26 + 6*sqrt(33))^(1/3))*(26 - 6*sqrt(33) + (-13 + 13*im*sqrt(3) - 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-4 - 4*im*sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x)))*sqrt(1 + 4*x + 2*x^2 + x^4)*Elliptic.E(asin(sqrt(26 - 6*sqrt(33) + (-13 - 13*im*sqrt(3) + 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + 4*im*(im + sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x)/(sqrt((39 + 13*im*sqrt(3) - 9*im*sqrt(11) - 9*sqrt(33) + 4*(3 - im*sqrt(3))*(-26 + 6*sqrt(33))^(1/3))/(39 - 13*im*sqrt(3) + 9*im*sqrt(11) - 9*sqrt(33) + 4*(3 + im*sqrt(3))*(-26 + 6*sqrt(33))^(1/3)))*sqrt(26 - 6*sqrt(33) + (-13 + 13*im*sqrt(3) - 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-4 - 4*im*sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x))), (4*(21 + 7*im*sqrt(3) - 3*im*sqrt(11) - 3*sqrt(33)) + (3 - im*sqrt(3) - 3*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3))/(4*(21 - 7*im*sqrt(3) + 3*im*sqrt(11) - 3*sqrt(33)) + (3 + im*sqrt(3) + 3*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3))))/((4*2^(2/3) - (-13 + 3*sqrt(33))^(1/3) - 2^(1/3)*(-13 + 3*sqrt(33))^(2/3) + 3*(-13 + 3*sqrt(33))^(1/3)*x)*sqrt((im*(1 + x))/((104 - 24*sqrt(33) + (-13 - 13*im*sqrt(3) + 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + 4*im*(im + sqrt(3))*(-26 + 6*sqrt(33))^(2/3))*(26 - 6*sqrt(33) + (-13 + 13*im*sqrt(3) - 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-4 - 4*im*sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x)))*sqrt(26 - 6*sqrt(33) + (-13 + 13*im*sqrt(3) - 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-4 - 4*im*sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x)*sqrt(26 - 6*sqrt(33) + (-13 - 13*im*sqrt(3) + 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + 4*im*(im + sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x)) + ((2^(1/3)*(13 - 13*im*sqrt(3) + 9*im*sqrt(11) - 3*sqrt(33)) + 4*2^(2/3)*(1 + im*sqrt(3))*(-13 + 3*sqrt(33))^(1/3) + 20*(-13 + 3*sqrt(33))^(2/3))*(4*2^(2/3)*(im + sqrt(3)) + 8*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(-im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3))*sqrt((52 - 12*sqrt(33) - 2^(1/3)*(-13 + 3*sqrt(33))^(4/3) + 4*(-26 + 6*sqrt(33))^(2/3))/(-13 + 3*sqrt(33) + 4*(-26 + 6*sqrt(33))^(1/3)))*sqrt((1/(1 + x))*(-8*im*(-13 + 3*sqrt(33)) + (-43*im - 13*sqrt(3) + 9*sqrt(11) + 5*im*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (2*im + 4*sqrt(3) - 2*im*sqrt(33))*(-26 + 6*sqrt(33))^(2/3) + (8*im*(-13 + 3*sqrt(33)) + (13*im - 13*sqrt(3) + 9*sqrt(11) - 3*im*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + 4*(im + sqrt(3))*(-26 + 6*sqrt(33))^(2/3))*x))*sqrt(1 + 4*x + 2*x^2 + x^4)*Elliptic.F(asin((sqrt(52 - 12*sqrt(33) - 2^(1/3)*(-13 + 3*sqrt(33))^(4/3) + 4*(-26 + 6*sqrt(33))^(2/3))*sqrt(26 - 6*sqrt(33) + (-13 - 13*im*sqrt(3) + 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + 4*im*(im + sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x))/(2^(1/6)*sqrt(3)*(-13 + 3*sqrt(33))^(2/3)*sqrt(39 + 13*im*sqrt(3) - 9*im*sqrt(11) - 9*sqrt(33) + 4*(3 - im*sqrt(3))*(-26 + 6*sqrt(33))^(1/3))*sqrt(1 + x))), (4*(21*im - 7*sqrt(3) + 3*sqrt(11) - 3*im*sqrt(33)) + (3*im + sqrt(3) + 3*sqrt(11) + 3*im*sqrt(33))*(-26 + 6*sqrt(33))^(1/3))/(-56*sqrt(3) + 24*sqrt(11) + 2*(sqrt(3) + 3*sqrt(11))*(-26 + 6*sqrt(33))^(1/3))))/(3*2^(2/3)*3^(3/4)*(-13 + 3*sqrt(33))^(1/3)*sqrt(39 + 13*im*sqrt(3) - 9*im*sqrt(11) - 9*sqrt(33) + 4*(3 - im*sqrt(3))*(-26 + 6*sqrt(33))^(1/3))*sqrt(1 + x)*(4*2^(2/3)*(-im + sqrt(3)) - 2*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3) + 6*im*(-13 + 3*sqrt(33))^(1/3)*x)*sqrt(26 - 6*sqrt(33) + (-13 - 13*im*sqrt(3) + 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + 4*im*(im + sqrt(3))*(-26 + 6*sqrt(33))^(2/3) + 6*(-13 + 3*sqrt(33))*x)*sqrt((8*(-13 + 3*sqrt(33)) - (5 - 3*im*sqrt(3) + 3*im*sqrt(11) + sqrt(33))*(-26 + 6*sqrt(33))^(2/3) + (-26 + 6*sqrt(33))^(1/3)*(-41 + 15*im*sqrt(3) - 3*im*sqrt(11) + 7*sqrt(33)) + (104 - 24*sqrt(33) + (-13 - 13*im*sqrt(3) + 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + 4*im*(im + sqrt(3))*(-26 + 6*sqrt(33))^(2/3))*x)/((-39 - 13*im*sqrt(3) + 9*im*sqrt(11) + 9*sqrt(33) + 4*im*(3*im + sqrt(3))*(-26 + 6*sqrt(33))^(1/3))*(1 + x)))) + ((4*2^(2/3) + 2*(-13 + 3*sqrt(33))^(1/3) - 2^(1/3)*(-13 + 3*sqrt(33))^(2/3))*(4*2^(2/3)*(im + sqrt(3)) - 4*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(-im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3))*(4*2^(2/3)*(-im + sqrt(3)) + 4*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3))*sqrt((-39 + 13*im*sqrt(3) - 9*im*sqrt(11) + 9*sqrt(33) - 4*im*(-3*im + sqrt(3))*(-26 + 6*sqrt(33))^(1/3))/(104 - 24*sqrt(33) + (-13 + 13*im*sqrt(3) - 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-4 - 4*im*sqrt(3))*(-26 + 6*sqrt(33))^(2/3)))*sqrt(1 + x)*sqrt((104 - 24*sqrt(33) + 2*(1 + 14*im*sqrt(3) - 6*im*sqrt(11) + sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-7 - im*sqrt(3) - 3*im*sqrt(11) + sqrt(33))*(-26 + 6*sqrt(33))^(2/3) + 2*(-52 + 12*sqrt(33) + 2^(1/3)*(-13 + 3*sqrt(33))^(4/3) - 4*(-26 + 6*sqrt(33))^(2/3))*x)/((-39 + 13*im*sqrt(3) - 9*im*sqrt(11) + 9*sqrt(33) - 4*im*(-3*im + sqrt(3))*(-26 + 6*sqrt(33))^(1/3))*(1 + x)))*sqrt((104 - 24*sqrt(33) + 2*(1 - 14*im*sqrt(3) + 6*im*sqrt(11) + sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-7 + im*sqrt(3) + 3*im*sqrt(11) + sqrt(33))*(-26 + 6*sqrt(33))^(2/3) + 2*(-52 + 12*sqrt(33) + 2^(1/3)*(-13 + 3*sqrt(33))^(4/3) - 4*(-26 + 6*sqrt(33))^(2/3))*x)/((-39 - 13*im*sqrt(3) + 9*im*sqrt(11) + 9*sqrt(33) + 4*im*(3*im + sqrt(3))*(-26 + 6*sqrt(33))^(1/3))*(1 + x)))*sqrt(1 + 4*x + 2*x^2 + x^4)*Elliptic.Pi((2^(1/3)*(4*2^(1/3)*(-3*im + sqrt(3)) + (3*im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3)))/(4*2^(2/3)*(-im + sqrt(3)) - 8*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3)), asin(sqrt(13 - 3*sqrt(33) - 2^(1/3)*(-13 + 3*sqrt(33))^(4/3) + 4*(-26 + 6*sqrt(33))^(2/3) + (-39 + 9*sqrt(33))*x)/(2^(1/6)*sqrt(3)*(-13 + 3*sqrt(33))^(2/3)*sqrt((-39 + 13*im*sqrt(3) - 9*im*sqrt(11) + 9*sqrt(33) - 4*im*(-3*im + sqrt(3))*(-26 + 6*sqrt(33))^(1/3))/(104 - 24*sqrt(33) + (-13 + 13*im*sqrt(3) - 9*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3) + (-4 - 4*im*sqrt(3))*(-26 + 6*sqrt(33))^(2/3)))*sqrt(1 + x))), (4*(21 - 7*im*sqrt(3) + 3*im*sqrt(11) - 3*sqrt(33)) + (3 + im*sqrt(3) + 3*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3))/(4*(21 + 7*im*sqrt(3) - 3*im*sqrt(11) - 3*sqrt(33)) + (3 - im*sqrt(3) - 3*im*sqrt(11) + 3*sqrt(33))*(-26 + 6*sqrt(33))^(1/3))))/(2^(1/6)*sqrt(3)*(4*2^(2/3)*(im + sqrt(3)) + 2*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(-im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3) - 6*im*(-13 + 3*sqrt(33))^(1/3)*x)*(4*2^(2/3)*(-im + sqrt(3)) - 2*im*(-13 + 3*sqrt(33))^(1/3) + 2^(1/3)*(im + sqrt(3))*(-13 + 3*sqrt(33))^(2/3) + 6*im*(-13 + 3*sqrt(33))^(1/3)*x)*sqrt(13 - 3*sqrt(33) - 2^(1/3)*(-13 + 3*sqrt(33))^(4/3) + 4*(-26 + 6*sqrt(33))^(2/3) + (-39 + 9*sqrt(33))*x)))]

    #= It*is*interesting*to*see*how*much*of*this*one*can*be*done; =#

    @test_int [(12*log(x/mc^2)*x^2*pi^2*mc^3*(-8*x-12*mc^2+3*mc) + pi^2*(12*x^4*mc+3*x^4+176*x^3*mc^3-24*x^3*mc^2 - 144*x^2*mc^5-48*x*mc^7+24*x*mc^6+4*mc^9-3*mc^8))/(384*ℯ^(x/y)*x^2), x, 20, ((3 - 4*mc)*mc^8*pi^2)/(ℯ^(x/y)*(384*x)) + ((3/8)*mc^5*pi^2*y)/ℯ^(x/y) + ((1/48)*(3 - 22*mc)*mc^2*pi^2*x*y)/ℯ^(x/y) - ((1/128)*(1 + 4*mc)*pi^2*x^2*y)/ℯ^(x/y) + ((1/48)*(3 - 22*mc)*mc^2*pi^2*y^2)/ℯ^(x/y) + ((1/4)*mc^3*pi^2*y^2)/ℯ^(x/y) - ((1/64)*(1 + 4*mc)*pi^2*x*y^2)/ℯ^(x/y) - ((1/64)*(1 + 4*mc)*pi^2*y^3)/ℯ^(x/y) + (1/16)*(1 - 2*mc)*mc^6*pi^2*ExpIntegralEi(-(x/y)) + ((3 - 4*mc)*mc^8*pi^2*ExpIntegralEi(-(x/y)))/(384*y) + (1/32)*mc^3*pi^2*(3*mc - 12*mc^2 - 8*y)*y*ExpIntegralEi(-(x/y)) - ((1/32)*mc^3*pi^2*(3*(1 - 4*mc)*mc - 8*x)*y*log(x/mc^2))/ℯ^(x/y) + ((1/4)*mc^3*pi^2*y^2*log(x/mc^2))/ℯ^(x/y)]

    #= The*following*integrals*reveal*deficiencies*in*the*current*integrator; =#

    @test_int [sin(2*x)/cos(x), x, 2, -2*cos(x)]

    #= This*example, which*appeared*in*Tobey's*thesis, needs*factorization =#
    #= over*algebraic*fields. It*currently*gives*an*ugly*answer*and*so*has =#
    #= been*suppressed; =#

    @test_int [(7*x^13+10*x^8+4*x^7-7*x^6-4*x^3-4*x^2+3*x+3)/(x^14-2*x^8-2*x^7-2*x^4-4*x^3-x^2+2*x+1), x, -5, (1/2)*((1 + sqrt(2))*log(1 + x + sqrt(2)*x + sqrt(2)*x^2 - x^7) - (-1 + sqrt(2))*log(-1 + (-1 + sqrt(2))*x + sqrt(2)*x^2 + x^7))]
    end

max_stars_repo_path: test/0 Independent test suites/Hearn Problems.jl


text:
"""
    sumr(;binfile::AbstractString="", outdir::AbstractString=".", pseudocount::Number=1.0)

Extract the summary information of data matrix.

Input Arguments
---------
- `binfile` is a Julia Binary file generated by `csv2bin` function.
- `outdir` is specified the directory you want to save the result.
- `pseudocount` is specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.

Output Files
---------
- `Sample_NoCounts.csv` : Sum of counts in each column.
- `Feature_Means.csv` : Mean in each row.
- `Feature_LogMeans.csv` : Log10(Mean+pseudocount) in each row.
- `Feature_FTTMeans.csv` : FTT(Mean+pseudocount) in each row.
- `Feature_Vars.csv` : Sample variance in each row.
- `Feature_LogVars.csv` : Log10(Var+pseudocount) in each row.
- `Feature_FTTVars.csv` : FTT(Var+pseudocount) in each row.
- `Feature_CV2s.csv` : Coefficient of Variation in each row.
- `Feature_NoZeros.csv` : Number of zero-elements in each row.
"""

# Column-wise statistics
function nocounts(binfile::AbstractString)
    N, M = nm(binfile)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    nc = zeros(UInt32, M)
    progress = Progress(N)
    open(binfile) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            # Update
            nc = nc .+ x
            # Progress Bar
            next!(progress)
        end
        close(stream)
    end
    return nc
end

# Row-wise statistics
function stats(binfile::AbstractString, pseudocount::Number, nc::AbstractArray)
    N, M = nm(binfile)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    x_div = zeros(UInt32, M)
    x_cpm = zeros(UInt32, M)
    x_cpt = zeros(UInt32, M)
    x_cpmed = zeros(UInt32, M)

    m = zeros(N)
    lm = zeros(N)
    fttm = zeros(N)

    v = zeros(N)
    lv = zeros(N)
    fttv = zeros(N)

    cpmm = zeros(N)
    lcpmm = zeros(N)
    fttcpmm = zeros(N)

    cpmv = zeros(N)
    lcpmv = zeros(N)
    fttcpmv = zeros(N)

    cptm = zeros(N)
    lcptm = zeros(N)
    fttcptm = zeros(N)

    cptv = zeros(N)
    lcptv = zeros(N)
    fttcptv = zeros(N)

    cpmedm = zeros(N)
    lcpmedm = zeros(N)
    fttcpmedm = zeros(N)

    cpmedv = zeros(N)
    lcpmedv = zeros(N)
    fttcpmedv = zeros(N)

    c = zeros(N)
    nz = zeros(N)
    progress = Progress(N)
    open(binfile) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            # Update
            m[n] = mean(x)
            lm[n] = mean(log10.(x .+ pseudocount))
            fttm[n] = mean(sqrt.(x) .+ sqrt.(x .+ 1))

            v[n] = var(x)
            lv[n] = var(log10.(x .+ pseudocount))
            fttv[n] = var(sqrt.(x) .+ sqrt.(x .+ 1))

            x_div = x ./ nc

            x_cpm = 1e6 .* x_div
            cpmm[n] = mean(x_cpm)
            lcpmm[n] = mean(log10.(x_cpm .+ pseudocount))
            fttcpmm[n] = mean(sqrt.(x_cpm) .+ sqrt.(x_cpm .+ 1))

            cpmv[n] = var(x_cpm)
            lcpmv[n] = var(log10.(x_cpm .+ pseudocount))
            fttcpmv[n] = var(sqrt.(x_cpm) .+ sqrt.(x_cpm .+ 1))

            x_cpt = 1e4 .* x_div
            cptm[n] = mean(x_cpt)
            lcptm[n] = mean(log10.(x_cpt .+ pseudocount))
            fttcptm[n] = mean(sqrt.(x_cpt) .+ sqrt.(x_cpt .+ 1))

            cptv[n] = var(x_cpt)
            lcptv[n] = var(log10.(x_cpt .+ pseudocount))
            fttcptv[n] = var(sqrt.(x_cpt) .+ sqrt.(x_cpt .+ 1))

            x_cpmed = median(nc) .* x_div
            cpmedm[n] = mean(x_cpmed)
            lcpmedm[n] = mean(log10.(x_cpmed .+ pseudocount))
            fttcpmedm[n] = mean(sqrt.(x_cpmed) .+ sqrt.(x_cpmed .+ 1))

            cpmedv[n] = var(x_cpmed)
            lcpmedv[n] = var(log10.(x_cpmed .+ pseudocount))
            fttcpmedv[n] = var(sqrt.(x_cpmed) .+ sqrt.(x_cpmed .+ 1))

            c[n] = v[n] / m[n]^2
            for mm = 1:M
                if x[mm] != 0
                    nz[n] = nz[n] + 1
                end
            end
            # Progress Bar
            next!(progress)
        end
        close(stream)
    end
    return m, lm, fttm, v, lv, fttv, cpmm, lcpmm, fttcpmm, cpmv, lcpmv, fttcpmv, cptm, lcptm, fttcptm, cptv, lcptv, fttcptv, cpmedm, lcpmedm, fttcpmedm, cpmedv, lcpmedv, fttcpmedv, c, nz
end

function sumr(;binfile::AbstractString="", outdir::AbstractString=".", pseudocount::Number=1.0)
    # 1 / 2 : Column-wise statistics
    println("1 / 2 : Column-wise statistics are calculated...")
    Sample_NoCounts = nocounts(binfile)

    # 2 / 2 : Row-wise statistics
    println("2 / 2 : Row-wise statistics are calculated...")
    Feature_Means, Feature_LogMeans, Feature_FTTMeans, Feature_Vars, Feature_LogVars, Feature_FTTVars, Feature_CPMMeans, Feature_LogCPMMeans, Feature_FTTCPMMeans, Feature_CPMVars, Feature_LogCPMVars, Feature_FTTCPMVars, Feature_CPTMeans, Feature_LogCPTMeans, Feature_FTTCPTMeans, Feature_CPTVars, Feature_LogCPTVars, Feature_FTTCPTVars, Feature_CPMEDMeans, Feature_LogCPMEDMeans, Feature_FTTCPMEDMeans, Feature_CPMEDVars, Feature_LogCPMEDVars, Feature_FTTCPMEDVars, Feature_CV2s, Feature_NoZeros = stats(binfile, pseudocount, Sample_NoCounts)

    # Save
    writecsv(joinpath(outdir, "Sample_NoCounts.csv"), Sample_NoCounts)

    writecsv(joinpath(outdir, "Feature_Means.csv"), Feature_Means)
    writecsv(joinpath(outdir, "Feature_LogMeans.csv"), Feature_LogMeans)
    writecsv(joinpath(outdir, "Feature_FTTMeans.csv"), Feature_FTTMeans)

    writecsv(joinpath(outdir, "Feature_CPMMeans.csv"), Feature_CPMMeans)
    writecsv(joinpath(outdir, "Feature_LogCPMMeans.csv"), Feature_LogCPMMeans)
    writecsv(joinpath(outdir, "Feature_FTTCPMMeans.csv"), Feature_FTTCPMMeans)

    writecsv(joinpath(outdir, "Feature_CPTMeans.csv"), Feature_CPTMeans)
    writecsv(joinpath(outdir, "Feature_LogCPTMeans.csv"), Feature_LogCPTMeans)
    writecsv(joinpath(outdir, "Feature_FTTCPTMeans.csv"), Feature_FTTCPTMeans)

    writecsv(joinpath(outdir, "Feature_CPMEDMeans.csv"), Feature_CPMEDMeans)
    writecsv(joinpath(outdir, "Feature_LogCPMEDMeans.csv"), Feature_LogCPMEDMeans)
    writecsv(joinpath(outdir, "Feature_FTTCPMEDMeans.csv"), Feature_FTTCPMEDMeans)

    writecsv(joinpath(outdir, "Feature_Vars.csv"), Feature_Vars)
    writecsv(joinpath(outdir, "Feature_LogVars.csv"), Feature_LogVars)
    writecsv(joinpath(outdir, "Feature_FTTVars.csv"), Feature_FTTVars)

    writecsv(joinpath(outdir, "Feature_CPMVars.csv"), Feature_CPMVars)
    writecsv(joinpath(outdir, "Feature_LogCPMVars.csv"), Feature_LogCPMVars)
    writecsv(joinpath(outdir, "Feature_FTTCPMVars.csv"), Feature_FTTCPMVars)

    writecsv(joinpath(outdir, "Feature_CPTVars.csv"), Feature_CPTVars)
    writecsv(joinpath(outdir, "Feature_LogCPTVars.csv"), Feature_LogCPTVars)
    writecsv(joinpath(outdir, "Feature_FTTCPTVars.csv"), Feature_FTTCPTVars)

    writecsv(joinpath(outdir, "Feature_CPMEDVars.csv"), Feature_CPMEDVars)
    writecsv(joinpath(outdir, "Feature_LogCPMEDVars.csv"), Feature_LogCPMEDVars)
    writecsv(joinpath(outdir, "Feature_FTTCPMEDVars.csv"), Feature_FTTCPMEDVars)

    writecsv(joinpath(outdir, "Feature_CV2s.csv"), Feature_CV2s)
    writecsv(joinpath(outdir, "Feature_NoZeros.csv"), Feature_NoZeros)
end


text:
# function : current state s, current action a, state space S maps to next state s2
function simplechain_transition(s,a,S)
    if a == 1  # go left
        s2 = s-1 # decrease state value by 1
    elseif a == 2 # go right
        s2 = s+1 # increase state value by 1
    elseif a == -1
        s2 = s
    else
        error("action must be {1,2}, but was : ", action)
    end
    if s2 > maximum(S)  # if s2 is larger than the maximum state
        s2 = minimum(S) # set it to the start state (smallest state)
    elseif s2 < minimum(S) # if s2 is smaller than smallest state
        s2 = minimum(S) # make it the smallest state
    end
    return s2  # return next state
end

function simplechain_perturb_action(a::Int, failchance)
    ϵ = rand()
    if ϵ ≤ 1.0 - failchance
        return a
    elseif ϵ ≤ 1.0 - failchance / 2
        return Int(-1)
    else
        return a==1 ? 2 : 1
    end
end

function create_simple_chain(num_states::Int; stochastic=false, failchance=0.1)
    S = 1:num_states
    A = 1:2

    if !stochastic
        p = (s,a)->simplechain_transition(s,a,S)  # mask out S in above function
    else
        p = (s,a)->simplechain_transition(s,simplechain_perturb_action(a,failchance), S)
    end
    d0 = ()->minimum(S)
    r = (s,a,s′)-> -1.0
    γ = (s,a,s′)-> (s==num_states && s′==1) ? 0.0 : 1.0
    meta = Dict{Symbol,Any}()
    meta[:minreward] = -1.0
    meta[:maxreward] = -1.0
    meta[:minreturn] = -Inf
    meta[:maxreturn] = -Float64(num_states)
    meta[:stochastic] = stochastic
    meta[:minhorizon] = num_states
    meta[:maxhorizon] = Inf
    meta[:discounted] = false
    m = MDP(S,A,p,r,γ,d0,meta, ()->nothing)
    return m
end

function finite_horizon_transition(m0, s, a, maxT)
    t = s[1]
    t += 1
    s′ = m0.p(s[2],a)
    if t > maxT || m0.γ(s[2], a, s′)==0.0
        t = 1
        s′ = m0.d0()
    end

    return (t, s′)
end

function create_simple_chain_finitetime(num_states::Int; stochastic=false, failchance=0.1, droptime=true)
    maxT = num_states * 20
    S = (1:maxT, 1:num_states)
    A = 1:2
    m0 = create_simple_chain(num_states; stochastic=stochastic, failchance=failchance)
    p = (s,a)->finite_horizon_transition(m0, s, a, maxT)
    γ = (s,a,s′)-> (s[2]==num_states && s′[2]==1) || (s[1] ≥ s′[1]) ? 0.0 : 1.0
    r = (s,a,s′)->m0.r(s[2],a,s′[2])
    d0 = ()->(1,m0.d0())
    meta = Dict{Symbol,Any}()
    meta[:minreward] = -1.0
    meta[:maxreward] = -1.0
    meta[:minreturn] = -Float64(maxT)
    meta[:maxreturn] = -Float64(num_states)
    meta[:stochastic] = stochastic
    meta[:minhorizon] = num_states
    meta[:maxhorizon] = maxT
    meta[:discounted] = false

    if droptime
        X = S[2]
        obs = s->s[2]
        m = POMDP(S,A,X,p,obs,r,γ,d0,meta,()->nothing)
    else
        m = MDP(S,A,p,r,γ,d0,meta,()->nothing)
    end

    return m
end


text:
# Source

IteratorInterfaceExtensions.isiterable(x::DataFrames.DataFrame) = true
TableTraits.isiterabletable(x::DataFrames.DataFrame) = true

function TableTraits.getiterator(df::DataFrames.DataFrame)
    return TableTraitsUtils.create_tableiterator(getfield(df, :columns), names(df))
end

# # Sink

function _DataFrame(x)
    cols, names = create_columns_from_iterabletable(x, na_representation=:missing)

    return DataFrames.DataFrame(cols, names)
end

DataFrames.DataFrame(x::AbstractVector{T}) where {T<:NamedTuple} = _DataFrame(x)

function DataFrames.DataFrame(x)
    if TableTraits.isiterabletable(x)
        return _DataFrame(x)
    else
        return convert(DataFrames.DataFrame, x)
    end
end


text:
interface_robot_chess(robot_decart) =
    (
        isborder = robot_decart.isborder,
        putmarker! = () ->  if iseven(sum(robot_decart.coordinates()))
                                robot_decart.putmarker!()
                            end,
        move! =  robot_decart.move!,
    )

function chess_mark!(r::Robot)
    robot = interface_robot(r)
    decart = interface_decart((x=0,y=0))

    robot_decart = interface_robot_decart(robot, decart)
    robot_chess = interface_robot_chess(robot_decart)
    line = interface_line(robot_chess.move!)
    trajectories = interface_trajectories(robot_chess)

    trajectories.snake!(Ost, Nord) do side
        robot_chess.putmarker!()
        line.movements!(robot_chess.putmarker!, side)
    end
end

max_stars_repo_name: Anpilogova/-

max_issues_repo_name: Anpilogova/-

max_forks_repo_name: Anpilogova/-


text:
export alias_elimination

function flatten(sys::ODESystem)
    if isempty(sys.systems)
        return sys
    else
        return ODESystem(equations(sys),
                         independent_variable(sys),
                         observed=observed(sys))
    end
end

using SymbolicUtils: Rewriters

function fixpoint_sub(x, dict)
    y = substitute(x, dict)
    while !isequal(x, y)
        y = x
        x = substitute(y, dict)
    end

    return x
end

function substitute_aliases(diffeqs, dict)
    lhss(diffeqs) .~ fixpoint_sub.(rhss(diffeqs), (dict,))
end

isvar(s::Sym) = !isparameter(s)
isvar(s::Term) = isvar(s.op)
isvar(s::Any) = false

function filterexpr(f, s)
    vs = []
    Rewriters.Prewalk(Rewriters.Chain([@rule((~x::f) => push!(vs, ~x))]))(s)
    vs
end

function make_lhs_0(eq)
    if eq.lhs isa Number && iszero(eq.lhs)
        return eq
    else
        0 ~ eq.lhs - eq.rhs
    end
end

function alias_elimination(sys::ODESystem)
    eqs = vcat(equations(sys), observed(sys))

    # make all algebraic equations have 0 on LHS
    eqs = map(eqs) do eq
        if eq.lhs isa Term && eq.lhs.op isa Differential
            eq
        else
            make_lhs_0(eq)
        end
    end

    newstates = map(eqs) do eq
        if eq.lhs isa Term && eq.lhs.op isa Differential
            filterexpr(isvar, eq.lhs)
        else
            []
        end
    end |> Iterators.flatten |> collect |> unique

    all_vars = map(eqs) do eq
        filterexpr(isvar, eq.rhs)
    end |> Iterators.flatten |> collect |> unique

    alg_idxs = findall(x->!(x.lhs isa Term) && iszero(x.lhs), eqs)

    eliminate = setdiff(all_vars, newstates)

    outputs = solve_for(eqs[alg_idxs], eliminate)

    diffeqs = eqs[setdiff(1:length(eqs), alg_idxs)]

    diffeqs′ = substitute_aliases(diffeqs, Dict(eliminate .=> outputs))

    ODESystem(diffeqs′, sys.iv, newstates, parameters(sys), observed=eliminate .~ outputs)
end


text:
### Solver options ###
opts = SolverOptions()
opts.square_root = false
opts.verbose=false
opts.constraint_tolerance = 1e-5
opts.cost_tolerance = 1e-6
model, obj = TrajectoryOptimization.Dynamics.pendulum
solver = Solver(model,obj,dt=0.1,opts=opts)

# -Initial state and control trajectories
X_interp = line_trajectory(solver)
U = ones(solver.model.m,solver.N-1)

results, stats = solve(solver,X_interp,U)

@test norm(results.X[end] - solver.obj.xf) < 1e-5
@test max_violation(results) < 1e-5
@test stats["iterations"] > stats["iterations (infeasible)"]
@test stats["outer loop iterations"] > stats["outer loop iterations (infeasible)"]
# Constraints
u_min = -2
u_max = 2
x_min = [-10;-10]
x_max = [10; 10]
obj_c = ConstrainedObjective(obj, u_min=u_min, u_max=u_max, x_min=x_min, x_max=x_max)

solver_con = Solver(model,obj_c,dt=0.1,opts=opts)

# Linear interpolation for state trajectory
X_interp = line_trajectory(solver_con)
U = ones(solver_con.model.m,solver_con.N-1)

results_con, stats_con = solve(solver_con,X_interp,U)
@test norm(results_con.X[end] - solver_con.obj.xf) < 1e-5
@test max_violation(results_con) < 1e-5
@test stats_con["iterations"] > stats_con["iterations (infeasible)"]
@test stats_con["outer loop iterations"] > stats_con["outer loop iterations (infeasible)"]


text:
@testset "upsample_nearest, integer scale via reshape" begin
    x = reshape(Float32[1. 2.; 3. 4.], (2,2,1,1))
    @test upsample_nearest(x, (3,3))[1,:] == [1,1,1, 2,2,2]

    y = upsample_nearest(x, (2,3))
    @test size(y) == (4,6,1,1)
    ∇upsample_nearest(y, (2,3)) == [6 12; 18 24]

    gradtest(x -> upsample_nearest(x, (2,3)), rand(2,2,1,1))

    y2 = upsample_nearest(x, size=(4,6))
    @test y ≈ y2
    gradtest(x -> upsample_nearest(x, size=(4,6)), rand(2,2,1,1))

    @test_throws ArgumentError ∇upsample_nearest(y, (2,4))
    @test_throws ArgumentError upsample_nearest(x, (1,2,3,4,5))
    @test_throws ArgumentError upsample_nearest(x, size=(3,4))
end

@testset "upsample_bilinear 2d" begin
    x = Float32[1 2; 3 4][:,:,:,:]
    x = cat(x,x; dims=3)
    x = cat(x,x; dims=4)

    # this output matches the one of pytorch v1.5.0
    # nn.UpsamplingBilinear2d(scale_factor=(3,2), align_corners=True)
    # for above x
    y_true = Float32[ 1//1  4//3   5//3   2//1;
                      7//5 26//15 31//15 12//5;
                      9//5 32//15 37//15 14//5;
                     11//5 38//15 43//15 16//5;
                     13//5 44//15 49//15 18//5;
                      3//1 10//3  11//3   4//1][:,:,:,:]
    y_true = cat(y_true,y_true; dims=3)
    y_true = cat(y_true,y_true; dims=4)

    y = upsample_bilinear(x, (3, 2))
    @test size(y) == size(y_true)
    @test eltype(y) == Float32
    @test y ≈ y_true

    gradtest(x->upsample_bilinear(x, (3, 2)), x, atol=1e-3) # works to higher precision for Float64

    # additional grad check, also compliant with pytorch
    o = ones(Float32,6,4,2,1)
    grad_true = 6*ones(Float32,2,2,2,1)
    @test ∇upsample_bilinear(o; size = (2,2)) ≈ grad_true

    y_true_2 = Rational{Int}[1//1  5//4  6//4  7//4 2//1;
                             3//2  7//4  8//4  9//4 5//2;
                             4//2  9//4 10//4 11//4 6//2;
                             5//2 11//4 12//4 13//4 7//2;
                             3//1 13//4 14//4 15//4 4//1][:,:,:,:]

     # check for real-valued single-number argument and type stability for rationals
    upsample_bilinear(x, 2.5) == y_true_2

    # check Integer support for forward pass
    # grads are always assumed to be floats, so no extension there
    x = UInt8[1 3; 3 5][:,:,:,:]
    y_true_int = UInt8[1 2 3; 2 3 4; 3 4 5][:,:,:,:]
    y = upsample_bilinear(x, 1.5)

    @test eltype(y) == UInt8
    @test y == y_true_int
end

@testset "pixel_shuffle" begin
    x = reshape(1:16, (2, 2, 4, 1))
    # [:, :, 1, 1] =
    #     1  3
    #     2  4
    # [:, :, 2, 1] =
    #     5  7
    #     6  8
    # [:, :, 3, 1] =
    #     9  11
    #     10  12
    # [:, :, 4, 1] =
    #     13  15
    #     14  16

    y_true = [1  9 3 11
              5 13 7 15
              2 10 4 12
              6 14 8 16][:,:,:,:]

    y = pixel_shuffle(x, 2)
    @test size(y) == size(y_true)
    @test y_true == y

    x = reshape(1:32, (2, 2, 8, 1))
    y_true = zeros(Int, 4, 4, 2, 1)
    y_true[:,:,1,1] .= [ 1   9  3  11
                         5  13  7  15
                         2  10  4  12
                         6  14  8  16 ]

    y_true[:,:,2,1] .= [ 17  25  19  27
                         21  29  23  31
                         18  26  20  28
                         22  30  24  32]

    y = pixel_shuffle(x, 2)
    @test size(y) == size(y_true)
    @test y_true == y

    x = reshape(1:4*3*27*2, (4,3,27,2))
    y = pixel_shuffle(x, 3)
    @test size(y) == (12, 9, 3, 2)
    # batch dimension is preserved
    x1 = x[:,:,:,[1]]
    x2 = x[:,:,:,[2]]
    y1 = pixel_shuffle(x1, 3)
    y2 = pixel_shuffle(x2, 3)
    @test cat(y1, y2, dims=4) == y

    for d in [1, 2, 3]
        r = rand(1:5)
        n = rand(1:5)
        c = rand(1:5)
        insize = rand(1:5, d)
        x = rand(insize..., r^d*c, n)

        y = pixel_shuffle(x, r)
        @test size(y) == ((r .* insize)..., c, n)

        gradtest(x -> pixel_shuffle(x, r), x)
    end
end


text:
module Quantica

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) Quantica

using Requires

function __init__()
      @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("plot_makie.jl")
      @require VegaLite = "112f6efa-9a02-5b7d-90c0-432ed331239a" include("plot_vegalite.jl")
end

using StaticArrays, NearestNeighbors, SparseArrays, LinearAlgebra, OffsetArrays,
      ProgressMeter, LinearMaps, Random

using ExprTools

using SparseArrays: getcolptr, AbstractSparseMatrix

export sublat, bravais, lattice, dims, sites, supercell, unitcell,
       hopping, onsite, @onsite!, @hopping!, parameters,
       hamiltonian, parametric, bloch, bloch!, optimize!, similarmatrix,
       flatten, wrap, transform!, combine,
       spectrum, bandstructure, marchingmesh, defaultmethod, bands, vertices,
       energies, states,
       momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM

export LatticePresets, RegionPresets, HamiltonianPresets

export LinearAlgebraPackage, ArpackPackage, KrylovKitPackage

export @SMatrix, @SVector, SMatrix, SVector, SA

export ishermitian, I

export SparseMatrixCSC

const NameType = Symbol
const nametype = Symbol

const TOOMANYITERS = 10^8

include("iterators.jl")
include("tools.jl")
include("presets.jl")
include("lattice.jl")
include("model.jl")
include("hamiltonian.jl")
include("parametric.jl")
include("mesh.jl")
include("diagonalizer.jl")
include("bandstructure.jl")
include("KPM.jl")
include("convert.jl")

end


text:
function compute_kinetics_vector(x,m,t)::Array{Float64,1}

    # initialize -
    κ = m["κ"]
    S = m["S"]
    number_of_rates = m["number_of_rates"]
    rate_vector = zeros(number_of_rates)

    # main loop -
    tmp_set = Set{Float64}()
    for j ∈ 1:number_of_rates

        # TODO: we should cache the index vectors for each rate in the model dictionary -
        idx_reactants = findall(x->x<0.0, S[:,j])
        for index ∈ idx_reactants

            σ = S[index,j]
            term = x[index]^(-σ)
            push!(tmp_set,term)
        end

        # compute the rate -
        rate_vector[j] = κ[j]*prod(tmp_set)

        # clean-up
        empty!(tmp_set)
    end

    # return -
    return rate_vector
end


text:
using StatsBase
using StatsPlots
using Random

"""
    compute_running_proportion(p_heads::Float64, N::Int)

Computes the running proportion of heads when 0 < p_heads < 1 is the probability of getting an head
and N is the number of trials.
"""
function compute_running_proportion(p_heads::Float64, N::Int)
    @assert p_heads >= 0 && p_heads <= 1
    flip_sequence = sample(0:1,pweights([1-p_heads,p_heads]), N)
    r = cumsum(flip_sequence)
    n = 1:N
    run_prop = r ./ n
    return run_prop
end

"""
    plot_running_proportion(p_heads::Float64, N::Int)

Reproduces figure 4.1 with custom probability of heads and number of flips.
"""
function plot_running_proportion(p_heads::Float64, N::Int)
    run_prop = compute_running_proportion(p_heads, N)
    n = 1:N
    plot(n,run_prop,
        title="Running Proportion of Heads",
        legend=false,
        marker = :dot,
        xscale = :log10)
    plot!(ylims=(0,1))
    # Add annotation with end proportion
    annotation = "End Proportion = " * string(run_prop[end])
    annotate!(N/4,0.7,annotation)
    # Add reference line (true probability)
    plot!(n, fill(p_heads, N), linestyle = :dot)
    xlabel!("Flip Number")
    ylabel!("Proportion Heads")
end


text:
"""
$(SIGNATURES)

Circular path for a swept mesh.
Specify the radius by evaluating it ahead of time.
For example, define `my_circle(θ) = f_swept_circle(0.1, θ)`, and then pass `my_circle` to [`create_swept_mesh`](@ref).
"""
function f_swept_circle(r::Float64, θ::Float64)
    n̂_1 = SVector{3,Float64}( cos(θ), sin(θ), 0.0)
    n̂_2 = SVector{3,Float64}(-sin(θ), cos(θ), 0.0)
    return r * n̂_1, n̂_1, n̂_2
end

"""
$(SIGNATURES)

Straight path for a swept mesh.
Pass this to [`create_swept_mesh`](@ref) to make a basic swept mesh.
"""
function f_swept_triv(θ::Float64)
	n̂_1 = SVector(0.0, 0.0, -1.0)
	n̂_2 = SVector(0.0, 1.0,  0.0)
	return n̂_2 * θ, n̂_1, n̂_2
end

function add_rot_sym_segment!(eM, fun_gen, θ, ϕ, rad, is_open::NTuple{2,Bool})
    # x̂ is in the "radial" direction
    # ŷ is along the path

    p1, x̂1, ŷ1 = fun_gen(θ[1])
    p2, x̂2, ŷ2 = fun_gen(θ[2])
    p3 = (p1 + p2) * 0.5
    p4 = p1 + AngleAxis(ϕ[1], ŷ1...) * x̂1 * rad[1]
    p6 = p1 + AngleAxis(ϕ[2], ŷ1...) * x̂1 * rad[1]
    p5 = p2 + AngleAxis(ϕ[1], ŷ2...) * x̂2 * rad[2]
    p7 = p2 + AngleAxis(ϕ[2], ŷ2...) * x̂2 * rad[2]
    n_offset = length(eM.point)
    i1 = SVector{4,Int64}(1,3,4,6) .+ n_offset
    i2 = SVector{4,Int64}(3,2,5,7) .+ n_offset
    i3 = SVector{4,Int64}(3,4,6,7) .+ n_offset
    i4 = SVector{4,Int64}(4,3,5,7) .+ n_offset
    i5 = SVector{3,Int64}(4,6,7) .+ n_offset
    i6 = SVector{3,Int64}(4,7,5) .+ n_offset
    append!(eM.point, [p1, p2, p3, p4, p5, p6, p7])
    append!(eM.tri, [i5, i6])
    append!(eM.tet, [i1, i2, i3, i4])
    ϵ = zeros(7)
    ϵ[1:3] .= 1.0
    if is_open[1]
        ϵ[1] = 0.0
        push!(eM.tri, SVector{3,Int64}(1,6,4) .+ n_offset)
    end
    if is_open[2]
        ϵ[2] = 0.0
        push!(eM.tri, SVector{3,Int64}(2,5,7) .+ n_offset)
    end
    append!(eM.ϵ, ϵ)
    return nothing
end

"""
$(SIGNATURES)

Creates a mesh by sweeping a 3D path for the given inputs:
* `fun_gen`: function that takes a single input arc length and outputs 1.) the position on path, 2.) a single direction normal to path and 3.) direction along the path.
* `lr`: arc length locations of nodes on path
* `rad`: thickness of swept mesh
* `n_side`: number of sides of swept path
* `is_open`: if the path starts and ends at the same point (e.g. ring) set this to false
* `rot_half`: if the sides of the path are "off" set this to false
"""
function create_swept_mesh(fun_gen::Function, lr::Union{LinRange,Vector{Float64}}, rad::Union{Float64,Vector{Float64}}, n_side::Int64=4,
		is_open::Bool=true; rot_half::Bool=true)

	l_lr = length(lr)
	if isa(rad, Float64)
		rad = zeros(l_lr) .+ rad
	else
		(l_lr == length(rad)) || error("the length of lr and length of rad must be the same")
	end
    eM = eMesh{Tri,Tet}()
    Δ_ϕ = 2 * pi / n_side
    rad = rad ./ cos(Δ_ϕ / 2)
	n_θ = l_lr - 1
    for k_θ = 1:n_θ
        θ = (lr[k_θ], lr[k_θ + 1])
		rad_k = (rad[k_θ], rad[k_θ + 1])
        for k_ϕ = 1:n_side
            ϕ_0 = Δ_ϕ * (k_ϕ  - 0.5 * rot_half)
            ϕ_1 = ϕ_0 + Δ_ϕ
            ϕ = (ϕ_0, ϕ_1)
            if is_open
                is_open_ = (k_θ==1, k_θ==n_θ)
            else
                is_open_ = (false, false)
            end
            add_rot_sym_segment!(eM, fun_gen, θ, ϕ, rad_k, is_open_)
        end
    end
	remove_degenerate!(eM)
    mesh_repair!(eM)
    return eM
end

function f_swept_helix(θ::Float64, coil_sep::Float64)
    delta_z = 1 / (2 * pi) * coil_sep
    r = SVector{3,Float64}(cos(θ), sin(θ), θ * delta_z)
    dir_1 = normalize(SVector{3,Float64}(-sin(θ), cos(θ), delta_z))
    dir_2 = SVector{3,Float64}( cos(θ), sin(θ), 0.0)
    return r, dir_1, dir_2
end

f_swept_circle(θ::Float64) = f_swept_helix(θ, 0.0)