text:
function blochmatrix1()

    A = BlochMatrix()
    B = BlochMatrix(0.0, 0, 0, 0, 0, 0, 0, 0, 0)

    show(devnull, "text/plain", A)
    @test eltype(A) == eltype(B)
    return A == B

end

function blochmatrix2()

    A = BlochSim.BlochDynamicsMatrix()
    B = BlochSim.BlochDynamicsMatrix(0f0, 0, 0)

    @test eltype(A) != eltype(B)
    return A == B

end

function blochmatrix3()

    A = FreePrecessionMatrix()
    B = FreePrecessionMatrix(false, 0, 0)
    C = BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9)
    copyto!(C, A)

    show(devnull, "text/plain", A)
    @test eltype(A) != eltype(B) == eltype(C)
    return A == B == C

end

function sum1()

    M1 = MagnetizationMC((1, 2, 3), (1, 2, 3))
    M2 = MagnetizationMC((0.5, 0.5, 0.5), (0.1, 0.1, 0.1))
    M = M1 + M2
    correct = MagnetizationMC((1.5, 2.5, 3.5), (1.1, 2.1, 3.1))

    return M == correct

end

function subtract1()

    M1 = Magnetization(4, 3, 1)
    M = M1 - M1
    correct = Magnetization()

    return M == correct

end

function subtract2()

    M1 = MagnetizationMC((1, 2, 3), (1, 2, 3))
    M2 = MagnetizationMC((0.5, 0.5, 0.5), (0.1, 0.1, 0.1))
    M = M1 - M2
    correct = MagnetizationMC((0.5, 1.5, 2.5), (0.9, 1.9, 2.9))

    return M == correct

end

function subtract3()

    A = BlochMatrix(1, 0, 0, 0, 1, 0, 0, 0, 1)
    C = I - A
    correct = BlochMatrix()

    return C == correct

end

function times1()

    M1 = MagnetizationMC((1, 1, 1), (2, 2, 2))
    M = M1 * 2
    correct = MagnetizationMC((2, 2, 2), (4, 4, 4))

    return M == correct

end

function times2()

    A = FreePrecessionMatrix(1, 1, 1)
    M = Magnetization(1, 1, 1)
    B = A * M
    correct = Magnetization(2, 0, 1)

    return B == correct

end

function times3()

    A = BlochMcConnellMatrix{Int}(2)
    for i = 1:2, j = 1:2
        A.A[i][j].a11 = 1
        A.A[i][j].a21 = 2
        A.A[i][j].a31 = 3
        A.A[i][j].a12 = 4
        A.A[i][j].a22 = 5
        A.A[i][j].a32 = 6
        A.A[i][j].a13 = 7
        A.A[i][j].a23 = 8
        A.A[i][j].a33 = 9
    end
    M = MagnetizationMC((0, 0, 1), (0, 0, 1))
    B = A * M
    correct = Matrix(A) * Vector(M)

    show(devnull, "text/plain", A)
    return Vector(B) == correct

end

function times4()

    A = ExcitationMatrix(BlochMatrix(2, 0, 0, 0, 2, 0, 0, 0, 2))
    M = MagnetizationMC((1, 1, 1), (2, 2, 3))
    B = A * M
    correct = MagnetizationMC((2, 2, 2), (4, 4, 6))

    show(devnull, "text/plain", A)
    return B == correct

end

function times5()

    M = idealspoiling * MagnetizationMC((1, 2, 1), (0.1, 0.2, 0.3))
    correct = MagnetizationMC((0, 0, 1), (0, 0, 0.3))

    return M == correct

end

function times6()

    A = BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9)
    B = FreePrecessionMatrix(1, 1, 1)
    C = A * B
    correct = Matrix(A) * Matrix(B)

    return Matrix(C) == correct

end

function times7()

    A = ExcitationMatrix(BlochMatrix(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9))
    C = A * idealspoiling
    correct = BlochMatrix(0, 0, 0, 0, 0, 0, 0.7, 0.8, 0.9)

    return C == correct

end

function times8()

    A = ExcitationMatrix(BlochMatrix(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9))
    B = BlochMcConnellMatrix(3)
    for i = 1:3, j = 1:3
        B.A[i][j].a11 = 1
        B.A[i][j].a21 = 2
        B.A[i][j].a31 = 3
        B.A[i][j].a12 = 4
        B.A[i][j].a22 = 5
        B.A[i][j].a32 = 6
        B.A[i][j].a13 = 7
        B.A[i][j].a23 = 8
        B.A[i][j].a33 = 9
    end
    C = A * B
    correct = kron(I(3), Matrix(A.A)) * Matrix(B)

    return Matrix(C) ≈ correct

end

function times9()

    A = idealspoiling
    B = BlochMcConnellMatrix(3)
    for i = 1:3, j = 1:3
        B.A[i][j].a11 = 1
        B.A[i][j].a21 = 2
        B.A[i][j].a31 = 3
        B.A[i][j].a12 = 4
        B.A[i][j].a22 = 5
        B.A[i][j].a32 = 6
        B.A[i][j].a13 = 7
        B.A[i][j].a23 = 8
        B.A[i][j].a33 = 9
    end
    C = A * B
    correct = kron(I(3), [0 0 0; 0 0 0; 0 0 1]) * Matrix(B)

    return Matrix(C) == correct

end

function rdiv1()

    M = MagnetizationMC((2.0, 4.0, 8.0), (2.0, 4.0, 6.0)) / 2
    correct = MagnetizationMC((1.0, 2.0, 4.0), (1.0, 2.0, 3.0))

    return M == correct

end

function ldiv1()

    M = Magnetization(1, 2, 3)
    A = BlochMatrix(1, 0, 0, 0, 1, 0, 0, 0, 1)
    B = A \ M
    correct = Magnetization(1, 2, 3)

    return B == correct

end

function ldiv2()

    M = MagnetizationMC((1, 2, 3), (1, 2, 3))
    A = BlochMcConnellMatrix(((
        BlochMatrix(1, 0, 0, 0, 1, 0, 0, 0, 1),
        BlochMatrix{Int}()), (
        BlochMatrix{Int}(),
        BlochMatrix(1, 0, 0, 0, 1, 0, 0, 0, 1)))
    )
    B = A \ M
    correct = MagnetizationMC((1, 2, 3), (1, 2, 3))

    return B == correct

end

function add1()

    M1 = MagnetizationMC((1, 1, 1), (2, 2, 2))
    M2 = MagnetizationMC((1, 2, 3), (1, 2, 3))
    V = Vector{Int}(undef, 6)
    add!(V, M1, M2)
    correct = [2, 3, 4, 3, 4, 5]

    return V == correct

end

function subtract1()

    M1 = MagnetizationMC((1, 1, 1), (2, 2, 2))
    M2 = MagnetizationMC((1, 2, 3), (1, 2, 3))
    subtract!(M2, M1)
    correct = MagnetizationMC((0, 1, 2), (-1, 0, 1))

    return M2 == correct

end

function mul1()

    M2 = Magnetization()
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    M1 = Magnetization(1, 2, 3)
    mul!(M2, A, M1)
    correct = Matrix(A) * Vector(M1)

    return Vector(M2) == correct

end

function mul2()

    A = BlochSim.BlochDynamicsMatrix()
    A.R1 = 1
    A.R2 = 2
    A.Δω = 3
    t = 10
    correct = Matrix(A) * t
    mul!(A, t)

    return Matrix(A) == correct

end

function mul3()

    A = BlochSim.ExchangeDynamicsMatrix()
    A.r = 2
    t = 10
    correct = Matrix(A) * t
    mul!(A, t)

    return Matrix(A) == correct

end

function mul4()

    T11 = 1
    T12 = 0.4
    T21 = 0.08
    T22 = 0.02
    τ12 = 0.05
    τ21 = 0.1

    A = BlochSim.BlochMcConnellDynamicsMatrix(2)
    A.A[1].R2 = -1 / T21 - 1 / τ12
    A.A[1].Δω = 2π
    A.A[1].R1 = -1 / T11 - 1 / τ12
    A.A[2].R2 = -1 / T22 - 1 / τ21
    A.A[2].Δω = 2π
    A.A[2].R1 = -1 / T12 - 1 / τ21
    A.E[1].r = 1 / τ12
    A.E[2].r = 1 / τ21
    t = 10
    correct = Matrix(A) * t
    mul!(A, t)

    return Matrix(A) == correct

end

function mul5()

    C = BlochMatrix()
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    t = 10
    mul!(C, A, t)
    correct = Matrix(A) * t

    return Matrix(C) == correct

end

function mul6()

    C = BlochMcConnellMatrix(2)
    A = BlochMcConnellMatrix(2)
    for i = 1:2, j = 1:2
        A.A[i][j].a11 = 1
        A.A[i][j].a21 = 2
        A.A[i][j].a31 = 3
        A.A[i][j].a12 = 4
        A.A[i][j].a22 = 5
        A.A[i][j].a32 = 6
        A.A[i][j].a13 = 7
        A.A[i][j].a23 = 8
        A.A[i][j].a33 = 9
    end
    t = 10
    mul!(C, A, t)
    correct = Matrix(A) * t

    return Matrix(C) == correct

end

function mul7()

    C = BlochMatrix()
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    B = BlochMatrix()
    B.a11 = 10
    B.a21 = 20
    B.a31 = 30
    B.a12 = 40
    B.a22 = 50
    B.a32 = 60
    B.a13 = 70
    B.a23 = 80
    B.a33 = 90
    correct = A * B
    mul!(C, A, B)

    return C == correct

end

function mul8()

    C = BlochMatrix()
    A = BlochSim.BlochDynamicsMatrix()
    A.R1 = 1
    A.R2 = 2
    A.Δω = 3
    B = BlochMatrix()
    B.a11 = 10
    B.a21 = 20
    B.a31 = 30
    B.a12 = 40
    B.a22 = 50
    B.a32 = 60
    B.a13 = 70
    B.a23 = 80
    B.a33 = 90
    correct = Matrix(A) * Matrix(B)
    mul!(C, A, B)

    return Matrix(C) == correct

end

function mul9()

    C = BlochMatrix()
    A = BlochSim.ExchangeDynamicsMatrix()
    A.r = 1
    B = BlochMatrix()
    B.a11 = 10
    B.a21 = 20
    B.a31 = 30
    B.a12 = 40
    B.a22 = 50
    B.a32 = 60
    B.a13 = 70
    B.a23 = 80
    B.a33 = 90
    correct = Matrix(A) * Matrix(B)
    mul!(C, A, B)

    return Matrix(C) == correct

end

function mul10()

    C = BlochMatrix()
    A = BlochSim.BlochDynamicsMatrix()
    A.R1 = 1
    A.R2 = 2
    A.Δω = 3
    B = BlochSim.BlochDynamicsMatrix()
    B.R1 = 10
    B.R2 = 20
    B.Δω = 30
    correct = Matrix(A) * Matrix(B)
    mul!(C, A, B)

    return Matrix(C) == correct

end

function mul11()

    C = BlochMatrix()
    A = BlochSim.ExchangeDynamicsMatrix()
    A.r = 2
    B = BlochSim.ExchangeDynamicsMatrix()
    B.r = 10
    correct = Matrix(A) * Matrix(B)
    mul!(C, A, B)

    return Matrix(C) == correct

end

function mul12()

    C = BlochMatrix()
    A = BlochSim.BlochDynamicsMatrix()
    A.R1 = 1
    A.R2 = 2
    A.Δω = 3
    B = BlochSim.ExchangeDynamicsMatrix()
    B.r = 10
    correct = Matrix(A) * Matrix(B)
    mul!(C, A, B)

    return Matrix(C) == correct

end

function mul13()

    C = BlochMatrix()
    A = BlochSim.BlochDynamicsMatrix()
    A.R1 = 1
    A.R2 = 2
    A.Δω = 3
    B = BlochSim.ExchangeDynamicsMatrix()
    B.r = 10
    correct = Matrix(B) * Matrix(A)
    mul!(C, B, A)

    return Matrix(C) == correct

end

function mul14()

    C = BlochMcConnellMatrix(2)
    A = BlochMcConnellMatrix(2)
    for i = 1:2, j = 1:2
        A.A[i][j].a11 = 1 + i + 10j
        A.A[i][j].a21 = 2 + i + 10j
        A.A[i][j].a31 = 3 + i + 10j
        A.A[i][j].a12 = 4 + i + 10j
        A.A[i][j].a22 = 5 + i + 10j
        A.A[i][j].a32 = 6 + i + 10j
        A.A[i][j].a13 = 7 + i + 10j
        A.A[i][j].a23 = 8 + i + 10j
        A.A[i][j].a33 = 9 + i + 10j
    end
    B = BlochMcConnellMatrix(2)
    for i = 1:2, j = 1:2
        B.A[i][j].a11 = 10 + i + 10j
        B.A[i][j].a21 = 20 + i + 10j
        B.A[i][j].a31 = 30 + i + 10j
        B.A[i][j].a12 = 40 + i + 10j
        B.A[i][j].a22 = 50 + i + 10j
        B.A[i][j].a32 = 60 + i + 10j
        B.A[i][j].a13 = 70 + i + 10j
        B.A[i][j].a23 = 80 + i + 10j
        B.A[i][j].a33 = 90 + i + 10j
    end
    mul!(C, A, B)
    correct = A * B

    return C == correct

end

function mul15()

    # Use Int because Float64 had slight floating point errors
    C = BlochMcConnellMatrix{Int}(2)
    A = BlochSim.BlochMcConnellDynamicsMatrix{Int}(2)
    A.A[1].R2 = 1
    A.A[1].Δω = 2
    A.A[1].R1 = 3
    A.A[2].R2 = 4
    A.A[2].Δω = 5
    A.A[2].R1 = 6
    A.E[1].r = 7
    A.E[2].r = 8
    B = BlochMcConnellMatrix{Int}(2)
    for i = 1:2, j = 1:2
        B.A[i][j].a11 = 10 + i + 10j
        B.A[i][j].a21 = 20 + i + 10j
        B.A[i][j].a31 = 30 + i + 10j
        B.A[i][j].a12 = 40 + i + 10j
        B.A[i][j].a22 = 50 + i + 10j
        B.A[i][j].a32 = 60 + i + 10j
        B.A[i][j].a13 = 70 + i + 10j
        B.A[i][j].a23 = 80 + i + 10j
        B.A[i][j].a33 = 90 + i + 10j
    end
    mul!(C, A, B)
    correct = Matrix(A) * Matrix(B)

    return Matrix(C) == correct

end

function mul16()

    T11 = 1
    T12 = 0.4
    T21 = 0.08
    T22 = 0.02
    τ12 = 0.05
    τ21 = 0.1

    C = BlochMcConnellMatrix(2)
    A = BlochSim.BlochMcConnellDynamicsMatrix(2)
    A.A[1].R2 = -1 / T21 - 1 / τ12
    A.A[1].Δω = 2π
    A.A[1].R1 = -1 / T11 - 1 / τ12
    A.A[2].R2 = -1 / T22 - 1 / τ21
    A.A[2].Δω = 2π
    A.A[2].R1 = -1 / T12 - 1 / τ21
    A.E[1].r = 1 / τ12
    A.E[2].r = 1 / τ21
    mul!(C, A, A)
    correct = Matrix(A) * Matrix(A)

    return Matrix(C) == correct

end

function mul17()

    M = MagnetizationMC((1, 2, 3), (4, 5, 6))
    mul!(M, 2)
    correct = MagnetizationMC((2, 4, 6), (8, 10, 12))

    return M == correct

end

function mul18()

    A = ExcitationMatrix(BlochMatrix(1, 1, 1, 2, 2, 2, 3, 3, 3))
    M = Magnetization(1, 2, 3)
    V = Vector{Int}(undef, 3)
    mul!(V, A, M)
    correct = [14, 14, 14]

    return V == correct

end

function mul19()

    C = BlochMatrix()
    A = FreePrecessionMatrix(1, 1, 1)
    B = ExcitationMatrix(BlochMatrix(1, 2, 3, 4, 5, 6, 7, 8, 9))
    mul!(C, A, B)
    correct = BlochMatrix(3, 1, 3, 9, 1, 6, 15, 1, 9)

    return C == correct

end

function div1()

    M = MagnetizationMC((1, 2, 3), (4, 5, 6.0))
    div!(M, 2)
    correct = MagnetizationMC((0.5, 1, 1.5), (2, 2.5, 3))

    return M == correct

end

function muladd1()

    M2 = Magnetization(1, 1, 1)
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    M1 = Magnetization(1, 2, 3)
    correct = Matrix(A) * Vector(M1) + Vector(M2)
    muladd!(M2, A, M1)

    return Vector(M2) == correct

end

function muladd2()

    C = BlochMatrix()
    C.a11 = 10
    C.a21 = 20
    C.a31 = 30
    C.a12 = 40
    C.a22 = 50
    C.a32 = 60
    C.a13 = 70
    C.a23 = 80
    C.a33 = 90
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    t = 10
    correct = Matrix(A) * t + Matrix(C)
    muladd!(C, A, t)

    return Matrix(C) == correct

end

function muladd3()

    C = BlochMcConnellMatrix(3)
    for i = 1:3, j = 1:3
        C.A[i][j].a11 = 10
        C.A[i][j].a21 = 20
        C.A[i][j].a31 = 30
        C.A[i][j].a12 = 40
        C.A[i][j].a22 = 50
        C.A[i][j].a32 = 60
        C.A[i][j].a13 = 70
        C.A[i][j].a23 = 80
        C.A[i][j].a33 = 90
    end
    A = BlochMcConnellMatrix(3)
    for i = 1:3, j = 1:3
        A.A[i][j].a11 = 1
        A.A[i][j].a21 = 2
        A.A[i][j].a31 = 3
        A.A[i][j].a12 = 4
        A.A[i][j].a22 = 5
        A.A[i][j].a32 = 6
        A.A[i][j].a13 = 7
        A.A[i][j].a23 = 8
        A.A[i][j].a33 = 9
    end
    t = 10
    correct = Matrix(A) * t + Matrix(C)
    muladd!(C, A, t)

    return Matrix(C) == correct

end

function muladd4()

    C = BlochMatrix()
    C.a11 = 10
    C.a21 = 20
    C.a31 = 30
    C.a12 = 40
    C.a22 = 50
    C.a32 = 60
    C.a13 = 70
    C.a23 = 80
    C.a33 = 90
    t = 10
    correct = I(3) * t + Matrix(C)
    muladd!(C, I, t)

    return Matrix(C) == correct

end

function muladd5()

    C = BlochMcConnellMatrix(3)
    for i = 1:3, j = 1:3
        C.A[i][j].a11 = 10
        C.A[i][j].a21 = 20
        C.A[i][j].a31 = 30
        C.A[i][j].a12 = 40
        C.A[i][j].a22 = 50
        C.A[i][j].a32 = 60
        C.A[i][j].a13 = 70
        C.A[i][j].a23 = 80
        C.A[i][j].a33 = 90
    end
    t = 10
    correct = I(9) * t + Matrix(C)
    muladd!(C, I, t)

    return Matrix(C) == correct

end

function muladd6()

    C = BlochMatrix()
    C.a11 = 1
    C.a21 = 1
    C.a31 = 1
    C.a12 = 1
    C.a22 = 1
    C.a32 = 1
    C.a13 = 1
    C.a23 = 1
    C.a33 = 1
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    B = BlochMatrix()
    B.a11 = 10
    B.a21 = 20
    B.a31 = 30
    B.a12 = 40
    B.a22 = 50
    B.a32 = 60
    B.a13 = 70
    B.a23 = 80
    B.a33 = 90
    correct = Matrix(A) * Matrix(B) + Matrix(C)
    muladd!(C, A, B)

    return Matrix(C) == correct

end

function muladd7()

    C = BlochMatrix()
    C.a11 = 1
    C.a21 = 1
    C.a31 = 1
    C.a12 = 1
    C.a22 = 1
    C.a32 = 1
    C.a13 = 1
    C.a23 = 1
    C.a33 = 1
    A = BlochSim.BlochDynamicsMatrix()
    A.R2 = 1
    A.Δω = 2
    A.R1 = 3
    B = BlochMatrix()
    B.a11 = 10
    B.a21 = 20
    B.a31 = 30
    B.a12 = 40
    B.a22 = 50
    B.a32 = 60
    B.a13 = 70
    B.a23 = 80
    B.a33 = 90
    correct = Matrix(A) * Matrix(B) + Matrix(C)
    muladd!(C, A, B)

    return Matrix(C) == correct

end

function muladd8()

    C = BlochMatrix()
    C.a11 = 1
    C.a21 = 1
    C.a31 = 1
    C.a12 = 1
    C.a22 = 1
    C.a32 = 1
    C.a13 = 1
    C.a23 = 1
    C.a33 = 1
    A = BlochSim.ExchangeDynamicsMatrix()
    A.r = 1
    B = BlochMatrix()
    B.a11 = 10
    B.a21 = 20
    B.a31 = 30
    B.a12 = 40
    B.a22 = 50
    B.a32 = 60
    B.a13 = 70
    B.a23 = 80
    B.a33 = 90
    correct = Matrix(A) * Matrix(B) + Matrix(C)
    muladd!(C, A, B)

    return Matrix(C) == correct

end

function muladd9()

    C = BlochMatrix()
    C.a11 = 1
    C.a21 = 1
    C.a31 = 1
    C.a12 = 1
    C.a22 = 1
    C.a32 = 1
    C.a13 = 1
    C.a23 = 1
    C.a33 = 1
    A = BlochSim.BlochDynamicsMatrix()
    A.R2 = 1
    A.Δω = 2
    A.R1 = 3
    B = BlochSim.BlochDynamicsMatrix()
    B.R2 = 10
    B.Δω = 20
    B.R1 = 30
    correct = Matrix(A) * Matrix(B) + Matrix(C)
    muladd!(C, A, B)

    return Matrix(C) == correct

end

function muladd10()

    C = BlochMatrix()
    C.a11 = 1
    C.a21 = 1
    C.a31 = 1
    C.a12 = 1
    C.a22 = 1
    C.a32 = 1
    C.a13 = 1
    C.a23 = 1
    C.a33 = 1
    A = BlochSim.ExchangeDynamicsMatrix()
    A.r = 1
    B = BlochSim.ExchangeDynamicsMatrix()
    B.r = 10
    correct = Matrix(A) * Matrix(B) + Matrix(C)
    muladd!(C, A, B)

    return Matrix(C) == correct

end

function muladd11()

    C = BlochMatrix()
    C.a11 = 1
    C.a21 = 1
    C.a31 = 1
    C.a12 = 1
    C.a22 = 1
    C.a32 = 1
    C.a13 = 1
    C.a23 = 1
    C.a33 = 1
    A = BlochSim.BlochDynamicsMatrix()
    A.R2 = 1
    A.Δω = 2
    A.R1 = 3
    B = BlochSim.ExchangeDynamicsMatrix()
    B.r = 10
    correct = Matrix(A) * Matrix(B) + Matrix(C)
    muladd!(C, A, B)

    return Matrix(C) == correct

end

function muladd12()

    C = BlochMatrix()
    C.a11 = 1
    C.a21 = 1
    C.a31 = 1
    C.a12 = 1
    C.a22 = 1
    C.a32 = 1
    C.a13 = 1
    C.a23 = 1
    C.a33 = 1
    A = BlochSim.BlochDynamicsMatrix()
    A.R2 = 1
    A.Δω = 2
    A.R1 = 3
    B = BlochSim.ExchangeDynamicsMatrix()
    B.r = 10
    correct = Matrix(B) * Matrix(A) + Matrix(C)
    muladd!(C, B, A)

    return Matrix(C) == correct

end

function subtractmul1()

    M2 = Magnetization()
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    M1 = Magnetization(10.0, 20.0, 30.0)
    subtractmul!(M2, I, A, M1)
    correct = (I - Matrix(A)) * Vector(M1)

    return Vector(M2) == correct

end

function subtractmul2()

    M2 = MagnetizationMC(2)
    A = BlochMcConnellMatrix(2)
    for i = 1:2, j = 1:2
        A.A[i][j].a11 = 1 + i + 10j
        A.A[i][j].a21 = 2 + i + 10j
        A.A[i][j].a31 = 3 + i + 10j
        A.A[i][j].a12 = 4 + i + 10j
        A.A[i][j].a22 = 5 + i + 10j
        A.A[i][j].a32 = 6 + i + 10j
        A.A[i][j].a13 = 7 + i + 10j
        A.A[i][j].a23 = 8 + i + 10j
        A.A[i][j].a33 = 9 + i + 10j
    end
    M1 = MagnetizationMC((10.0, 20.0, 30.0), (1.0, 2.0, 3.0))
    subtractmul!(M2, I, A, M1)
    correct = (I - Matrix(A)) * Vector(M1)

    return Vector(M2) == correct

end

function subtractmuladd1()

    M2 = Magnetization(1.0, 1.0, 1.0)
    A = BlochMatrix()
    A.a11 = 1
    A.a21 = 2
    A.a31 = 3
    A.a12 = 4
    A.a22 = 5
    A.a32 = 6
    A.a13 = 7
    A.a23 = 8
    A.a33 = 9
    M1 = Magnetization(10.0, 20.0, 30.0)
    subtractmul!(M2, I, A, M1)
    correct = (I - Matrix(A)) * Vector(M1)

    return Vector(M2) == correct

end

function getblock1()

    A = BlochSim.BlochMcConnellDynamicsMatrix(2)
    A.A[1].R2 = 1
    A.A[1].Δω = 2
    A.A[1].R1 = 3
    A.A[2].R2 = 10
    A.A[2].Δω = 20
    A.A[2].R1 = 30
    A.E[1].r = 100
    A.E[2].r = 200
    result = [Matrix(BlochSim.getblock(A, 1, 1)) Matrix(BlochSim.getblock(A, 1, 2));
              Matrix(BlochSim.getblock(A, 2, 1)) Matrix(BlochSim.getblock(A, 2, 2))]
    correct = Matrix(A)

    return result == correct

end

function getblock2()

    A = BlochMcConnellMatrix(2)
    for i = 1:2, j = 1:2
        A.A[i][j].a11 = 10 + i + 10j
        A.A[i][j].a21 = 20 + i + 10j
        A.A[i][j].a31 = 30 + i + 10j
        A.A[i][j].a12 = 40 + i + 10j
        A.A[i][j].a22 = 50 + i + 10j
        A.A[i][j].a32 = 60 + i + 10j
        A.A[i][j].a13 = 70 + i + 10j
        A.A[i][j].a23 = 80 + i + 10j
        A.A[i][j].a33 = 90 + i + 10j
    end
    result = [Matrix(BlochSim.getblock(A, 1, 1)) Matrix(BlochSim.getblock(A, 1, 2));
              Matrix(BlochSim.getblock(A, 2, 1)) Matrix(BlochSim.getblock(A, 2, 2))]
    correct = Matrix(A)

    return result == correct

end

function absolutesum1()

    A = BlochSim.BlochDynamicsMatrix()
    A.R1 = 1
    A.R2 = 2
    A.Δω = 3
    result = BlochSim.absolutesum(A)
    correct = norm(Matrix(A), 1)

    return result == correct

end

function absolutesum2()

    A = BlochSim.ExchangeDynamicsMatrix()
    A.r = 3
    result = BlochSim.absolutesum(A)
    correct = norm(Matrix(A), 1)

    return result == correct

end

function absolutesum3()

    A = BlochSim.BlochMcConnellDynamicsMatrix{Int}(2)
    A.A[1].R2 = 1
    A.A[1].Δω = 2
    A.A[1].R1 = 3
    A.A[2].R2 = 4
    A.A[2].Δω = 5
    A.A[2].R1 = 6
    A.E[1].r = 7
    A.E[2].r = 8
    result = BlochSim.absolutesum(A)
    correct = norm(Matrix(A), 1)

    return result == correct

end

@testset "Bloch Matrices" begin

    @test blochmatrix1()
    @test blochmatrix2()
    @test blochmatrix3()

end

@testset "Bloch Matrix Math" begin

    @testset "+" begin

        @test sum1()

    end

    @testset "-" begin

        @test subtract1()
        @test subtract2()
        @test subtract3()

    end

    @testset "*" begin

        @test times1()
        @test times2()
        @test times3()
        @test times4()
        @test times5()
        @test times6()
        @test times7()
        @test times8()
        @test times9()

    end

    @testset "/" begin

        @test rdiv1()

    end

    @testset "\\" begin

        @test ldiv1()
        @test ldiv2()

    end

    @testset "add!" begin

        @test add1()

    end

    @testset "subtract!" begin

        @test subtract1()

    end

    @testset "mul!" begin

        @test mul1()
        @test mul2()
        @test mul3()
        @test mul4()
        @test mul5()
        @test mul6()
        @test mul7()
        @test mul8()
        @test mul9()
        @test mul10()
        @test mul11()
        @test mul12()
        @test mul13()
        @test mul14()
        @test mul15()
        @test mul16()
        @test mul17()
        @test mul18()
        @test mul19()

    end

    @testset "div!" begin

        @test div1()

    end

    @testset "muladd!" begin

        @test muladd1()
        @test muladd2()
        @test muladd3()
        @test muladd4()
        @test muladd5()
        @test muladd6()
        @test muladd7()
        @test muladd8()
        @test muladd9()
        @test muladd10()
        @test muladd11()
        @test muladd12()

    end

    @testset "subtractmul!" begin

        @test subtractmul1()
        @test subtractmul2()

    end

    @testset "subtractmuladd!" begin

        @test subtractmuladd1()

    end

    @testset "getblock" begin

        @test getblock1()
        @test getblock2()

    end

    @testset "absolutesum" begin

        @test absolutesum1()
        @test absolutesum2()
        @test absolutesum3()

    end

end


text:
using LazySets, Plots

include("SparsePolynomialZonotope.jl")

# ==========
# Slide 4
# ==========
Z = Zonotope([0, 0.], [2 1 2; 0 2 2.])
fig = plot(Z)

# equivalent:
# Z = Zonotope([0, 0.], [[2, 0], [1, 2.], [2, 2.]])

# plot generators
plot!(fig, [0, 2.], [0, 0.], arrow=:arrow, c=:black, lw=2.0, lab="")
plot!(fig, [0, 1.], [0, 2.], arrow=:arrow, c=:black, lw=2.0, lab="")
plot!(fig, [0, 2.], [0, 2.], arrow=:arrow, c=:black, lw=2.0, lab="")
savefig(fig, "slide_4.png")

# ==========
# Slide 5
# ==========
SP = SparsePolynomialZonotope([0, 0.], [2 1 2; 0 2 2.], [1 0 3; 0 1 1.])
#fig = plot(SP)

# ==========
# Slide 10
# ==========

# ==========
# Slide 12
# ==========

# ==========
# Slide 13
# ==========

ext: jl


text:
"""
    SFAM.jl

Description:
    Options, structures, and logic for the Simplified Fuzzy ARTMAP (SFAM) module.

References:
    [1] G. A. Carpenter, S. Grossberg, N. Markuzon, J. H. Reynolds, and D. B. Rosen, “Fuzzy ARTMAP: A Neural Network Architecture for Incremental Supervised Learning of Analog Multidimensional Maps,” IEEE Trans. Neural Networks, vol. 3, no. 5, pp. 698-713, 1992, doi: 10.1109/72.159059.
"""

"""
    opts_SFAM(;kwargs)

Implements a Simple Fuzzy ARTMAP learner's options.

# Keyword Arguments
- `rho::Float`: vigilance value, [0, 1], default 0.75.
- `alpha::Float`: choice parameter, alpha > 0, default 1e-7.
- `epsilon::Float`: match tracking parameter, (0, 1), default 1e-3
- `beta::Float`: learning parameter, (0, 1], default 1.0.
- `uncommitted::Bool`: uncommitted node flag, default true.
- `display::Bool`: display flag, default true.
- `max_epoch::Int`: maximum number of epochs during training, default 1.
"""
@with_kw mutable struct opts_SFAM <: ARTOpts @deftype Float
    # Vigilance parameter: [0, 1]
    rho = 0.75; @assert rho >= 0.0 && rho <= 1.0
    # Choice parameter: alpha > 0
    alpha = 1e-7; @assert alpha > 0.0
    # Match tracking parameter
    epsilon = 1e-3; @assert epsilon > 0.0 && epsilon < 1.0
    # Learning parameter: (0, 1]
    beta = 1.0; @assert beta > 0.0 && beta <= 1.0
    # Uncommitted node flag
    uncommitted::Bool = true
    # Display flag
    display::Bool = true
    # Maximum number of epochs during training
    max_epochs::Int = 1
end # opts_SFAM()

"""
    SFAM <: ARTMAP

Simple Fuzzy ARTMAP struct.

For module options, see [`AdaptiveResonance.opts_SFAM`](@ref).

# Option Parameters
- `opts::opts_SFAM`: Simplified Fuzzy ARTMAP options struct.
- `config::DataConfig`: data configuration struct.

# Working Parameters
- `W::RealMatrix`: category weight matrix.
- `labels::IntegerVector`: incremental list of labels corresponding to each F2 node, self-prescribed or supervised.
- `n_categories::Int`: number of category weights (F2 nodes).
- `epoch::Int`: current training epoch.

# References
1. G. A. Carpenter, S. Grossberg, N. Markuzon, J. H. Reynolds, and D. B. Rosen, “Fuzzy ARTMAP: A Neural Network Architecture for Incremental Supervised Learning of Analog Multidimensional Maps,” IEEE Trans. Neural Networks, vol. 3, no. 5, pp. 698-713, 1992, doi: 10.1109/72.159059.
"""
mutable struct SFAM <: ARTMAP
    opts::opts_SFAM
    config::DataConfig
    W::RealMatrix
    labels::IntegerVector
    n_categories::Int
    epoch::Int
end # SFAM <: ARTMAP

"""
    SFAM()

Implements a Simple Fuzzy ARTMAP learner.

# Examples
```julia-repl
julia> SFAM()
SFAM
    opts: opts_SFAM
    ...
```
"""
function SFAM()
    opts = opts_SFAM()
    SFAM(opts)
end # SFAM()

"""
    SFAM(;kwargs...)

Implements a Simple Fuzzy ARTMAP learner with keyword arguments.

# Examples
```julia-repl
julia> SFAM()
SFAM
    opts: opts_SFAM
    ...
```
"""
function SFAM(;kwargs...)
    opts = opts_SFAM(;kwargs...)
    SFAM(opts)
end # SFAM(;kwargs...)

"""
    SFAM(opts)

Implements a Simple Fuzzy ARTMAP learner with specified options.

# Examples
```julia-repl
julia> opts = opts_SFAM()
julia> SFAM(opts)
SFAM
    opts: opts_SFAM
    ...
```
"""
function SFAM(opts::opts_SFAM)
    SFAM(
        opts,                           # opts_SFAM
        DataConfig(),                   # config
        Array{Float}(undef, 0, 0),      # W
        Array{Int}(undef, 0),           # labels
        0,                              # n_categories
        0                               # epoch
    )
end # SFAM(opts::opts_SFAM)

# SFAM incremental training method
function train!(art::SFAM, x::RealVector, y::Integer ; preprocessed::Bool=false)
    # Run the sequential initialization procedure
    sample = init_train!(x, art, preprocessed)

    # Initialization
    if !(y in art.labels)
        # Initialize W and labels
        if isempty(art.W)
            art.W = Array{Float64}(undef, art.config.dim_comp, 1)
            art.W[:, 1] = sample
        else
            art.W = [art.W sample]
        end
        push!(art.labels, y)
        art.n_categories += 1
    else
        # Baseline vigilance parameter
        rho_baseline = art.opts.rho

        # Compute activation function
        T = zeros(art.n_categories)
        for jx in 1:art.n_categories
            T[jx] = activation(art, sample, art.W[:, jx])
        end

        # Sort activation function values in descending order
        index = sortperm(T, rev=true)
        mismatch_flag = true
        for jx in 1:art.n_categories
            # Compute match function
            M = art_match(art, sample, art.W[:, index[jx]])
            # Current winner
            if M >= rho_baseline
                if y == art.labels[index[jx]]
                    # Learn
                    @debug "Learning"
                    art.W[:, index[jx]] = learn(art, sample, art.W[:, index[jx]])
                    mismatch_flag = false
                    break
                else
                    # Match tracking
                    @debug "Match tracking"
                    rho_baseline = M + art.opts.epsilon
                end
            end
        end

        # If we triggered a mismatch
        if mismatch_flag
            # Create new weight vector
            @debug "Mismatch"
            art.W = hcat(art.W, sample)
            push!(art.labels, y)
            art.n_categories += 1
        end
    end

    # ARTMAP guarantees correct training classification, so just return the label
    return y
end

# SFAM incremental classification method
function classify(art::SFAM, x::RealVector ; preprocessed::Bool=false, get_bmu::Bool=false)
    # Run the sequential initialization procedure
    sample = init_classify!(x, art, preprocessed)

    # Compute activation function
    T = zeros(art.n_categories)
    for jx in 1:art.n_categories
        T[jx] = activation(art, sample, art.W[:, jx])
    end

    # Sort activation function values in descending order
    index = sortperm(T, rev=true)
    mismatch_flag = true
    for jx in 1:art.n_categories
        # Compute match function
        M = art_match(art, sample, art.W[:, index[jx]])
        # Current winner
        if M >= art.opts.rho
            y_hat = art.labels[index[jx]]
            mismatch_flag = false
            break
        end
    end

    # If we did not find a resonant category
    if mismatch_flag
        @debug "Mismatch"
        # Report either the best matching unit or the mismatch label -1
        y_hat = get_bmu ? art.labels[index[1]] : -1
    end

    return y_hat
end

"""
    stopping_conditions(art::SFAM)

Stopping conditions for Simple Fuzzy ARTMAP, checked at the end of every epoch.
"""
function stopping_conditions(art::SFAM)
    # Compute the stopping condition, return a bool
    return art.epoch >= art.opts.max_epochs
end # stopping_conditions(art::SFAM)

"""
    learn(art::SFAM, x::RealVector, W::RealVector)

Returns a single updated weight for the Simple Fuzzy ARTMAP module for weight
vector W and sample x.
"""
function learn(art::SFAM, x::RealVector, W::RealVector)
    # Update W
    return art.opts.beta .* element_min(x, W) .+ W .* (1 - art.opts.beta)
end # learn(art::SFAM, x::RealVector, W::RealVector)

"""
    activation(art::SFAM, x::RealVector, W::RealVector)

Returns the activation value of the Simple Fuzzy ARTMAP module with weight W
and sample x.
"""
function activation(art::SFAM, x::RealVector, W::RealVector)
    # Compute T and return
    return norm(element_min(x, W), 1) / (art.opts.alpha + norm(W, 1))
end # activation(art::SFAM, x::RealVector, W::RealVector)

"""
    art_match(art::SFAM, x::RealVector, W::RealVector)

Returns the match function for the Simple Fuzzy ARTMAP module with weight W and
sample x.
"""
function art_match(art::SFAM, x::RealVector, W::RealVector)
    # Compute M and return
    return norm(element_min(x, W), 1) / art.config.dim
end # art_match(art::SFAM, x::RealVector, W::RealVector)


text:
abstract type AbstractTransform end

function logjacobian(g::AbstractTransform, λ::Tuple, y)
    return sum(x -> log∂(g, λ, x), y)
end

struct DoNothing <: AbstractTransform end

transform(::DoNothing, y) = y
log∂(::DoNothing, _) = 0

struct BoxCox <: AbstractTransform end

transform(::BoxCox, λ::Number, y) = λ == 0 ? log(y) : expm1(log(y) * λ) / λ
log∂(::BoxCox, λ::Number, y) = (λ - 1) * log(y)

struct YeoJohnson <: AbstractTransform end

function transform(::YeoJohnson, λ::Number, y)
    λ′ = y >= 0 ? λ : 2 - λ
    return sign(y) * (λ′ == 0 ? log1p(abs(y)) : expm1(log1p(abs(y)) * λ′) / λ′)
end
log∂(::YeoJohnson, λ::Number, y) = sign(y) * (λ - 1) *  log(abs(y))

struct SumTanh <: AbstractTransform end

function transform(::SumTanh, α::Vector, β::Vector, c::Vector, y)
    return sum(α .* tanh.(β .* (y .+ c)))
end
function log∂(::SumTanh, α::Vector, β::Vector, c::Vector, y)
    return sum(α .* β .* sech.(β .* (y .+ c)) .^ 2)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
include("graph/node.jl")

export  RepresentationGraph

# Representation graphs are used to indicate whether particular stages
# of an individual are synchronised with the genome.
type RepresentationGraph
  nodes::Dict{AbstractString, RepresentationGraphNode}

  # Produces the representation graph for a given species.
  RepresentationGraph(s::Species) =
    RepresentationGraph([RepresentationGraphNode(stv) for (stl, stv) in s.stages])

  # Constructs a representation graph from a list of provided
  # representation graph nodes.
  RepresentationGraph(nodes::Vector{RepresentationGraphNode}) = begin

    # Construct a map storing each node within the graph by its label.
    node_map = Dict{AbstractString, RepresentationGraphNode}()
    for n in nodes
      node_map[n.label] = n
    end

    # Determine the list of sources for each node within the graph.
    for (label, node) in node_map
      if !isroot(node)
        node.sources = [node_map[node.parent]]
      end
    end

    # Determine the list of destinations and Lamarckian links for
    # each node within the graph.
    for (label, node) in node_map
      if !isroot(node)
        push!(node_map[node.parent].destinations, node)
        if node.lamarckian
          push!(node_map[node.parent].sources, node)
          push!(node.destinations, node_map[node.parent])
        end
      end
    end

    # For each node, determine the list of nodes which it can ultimately reach,
    # and which nodes can ultimately reach it.
    for node in nodes
      node.can_reach = find_reachable(node)
      node.can_reach_from = find_reachable_from(node)
    end

    # Construct the representation graph using the prepared map of
    # nodes within the graph.
    return new(node_map)

  end
end

function describe(g::RepresentationGraph)
  for n in nodes(g)
    if n.clean
      status = "CLEAN"
    else
      status = "DIRTY"
    end
    println("$(n.label) <$(status)>")
  end
end

function path_to_clean(from::RepresentationGraphNode)
  path = RepresentationGraphNode[from]
  paths = Vector{RepresentationGraphNode}[path]
  while !path[1].clean
    path = pop!(paths)

    # Create a new branch for each unvisited node.
    for src in path[1].sources
      if src != path[1]
        push!(paths, [[src], path])
      end
    end
  end
  return [n.label for n in path]
end

# Returns a list of all the nodes which can be reached from a given node.
function find_reachable(node::RepresentationGraphNode)
  target = node
  q = [node]
  reachable = RepresentationGraphNode[]
  buffer = RepresentationGraphNode[]
  while !isempty(q)
    node = shift!(q)
    buffer = filter(d -> d != target && !in(d, reachable), node.destinations)
    reachable = vcat(reachable, buffer)
    q = vcat(q, buffer)
  end
  return reachable
end

# Returns a list of all the nodes from which this node can be reached.
function find_reachable_from(node::RepresentationGraphNode)
  target = node
  q = [node]
  reachable_from = RepresentationGraphNode[]
  buffer = RepresentationGraphNode[]
  while !isempty(q)
    node = shift!(q)
    buffer = filter(d -> d != target && !in(d, reachable_from), node.sources)
    reachable_from = vcat(reachable_from, buffer)
    q = vcat(q, buffer)
  end
  return reachable_from
end

# Returns a list of the nodes within a given representation graph.
nodes(g::RepresentationGraph) = collect(values(g.nodes))

# Initialises the representation graph for a given species by marking all
# nodes as dirty, except for the canonical genotype.
function initialize!(g::RepresentationGraph)
  for node in nodes(g)
    node.clean = isroot(node)
  end
  return g
end

# Produces a deep clone of a given representation graph.
clone(g::RepresentationGraph) =
  RepresentationGraph([copy(n) for n in nodes(g)])

# Returns all clean nodes within a given representation graph.
find_clean(g::RepresentationGraph) = filter(n -> n.clean, nodes(g))

# Calculates the chain of synchronisation operations required to sync a given
# stage within the representation graph. Marks all nodes along the chain as
# clean thereafter, before touching the edited node and dirtying any nodes
# reachable from the edited node.
function partial_sync!(g::RepresentationGraph, s::AbstractString)
  # Calculate the shortest path to synchronising the requested stage.
  path = path_to_clean(g.nodes[s])

  # Mark all the nodes along the path as clean.
  for label in path
    g.nodes[label].clean = true
  end

  # Mark all reachable nodes as dirty.
  for n in g.nodes[s].can_reach
    g.nodes[n.label].clean = false
  end

  # Return the chain of operations (as a vector of pairs) required to synchronise
  # the representation graph for the desired operation.
  return Tuple{AbstractString, AbstractString}[
    (path[i], path[i + 1]) for i in 1:(length(path) - 1)]
end

# Calculates the chain of synchronisation operations required to fully
# synchronise all the stages of a given representation graph. Marks all
# nodes within the graph as clean thereafter.
function full_sync!(g::RepresentationGraph)
  ops = Tuple{AbstractString, AbstractString}[]
  q = find_clean(g)
  while !isempty(q)
    n1 = pop!(q)
    for n2 in n1.destinations
      if !n2.clean
        n2.clean = true
        push!(q, n2)
        push!(ops, (n1.label, n2.label))
      end
    end
  end
  return ops
end


text:
mutable struct WasmtimeLinker
    wasmtime_linker_ptr::Ptr{wasmtime_linker_t}

    WasmtimeLinker(ptr::Ptr{wasmtime_linker_t}) = finalizer(wasmtime_linker_delete, new(ptr))
end

WasmtimeLinker(engine) = wasmtime_linker_new(engine) |> WasmtimeLinker

Base.unsafe_convert(::Type{Ptr{wasmtime_linker_t}}, linker::WasmtimeLinker) = linker.wasmtime_linker_ptr
Base.show(io::IO, ::WasmtimeLinker) = write(io, "WasmtimeLinker()")

define_wasi!(linker) = (@wt_check wasmtime_linker_define_wasi(linker));


text:
# Test list parsing
@test @julog(cons(x, cons(y, cons(z, cend())))) == @julog(list[x, y, z])
@test @julog(cons(x, cons(y, cons(z, W)))) == @julog(list[x, y, z | W])
@test @julog(list[x|[y|[z]]]) == @julog(list[x, y, z])

# Test list functionality with some basic list operations
clauses = @julog [
    member(X, [X | Y]) <<= true,
    member(X, [Y | YS]) <<= member(X, YS),
    append([], L, L) <<= true,
    append([X | XS], YS, [X | ZS]) <<= append(XS, YS, ZS),
    reverse(X, Y) <<= reverse(X, [], Y),
    reverse([], YS, YS) <<= true,
    reverse([X | XS], Accu, YS) <<= reverse(XS, [X | Accu], YS)
]

# Test list membership
@test resolve(@julog(member(banana, [avocado, banana, coconut])), clauses)[1] == true
@test resolve(@julog(member(durian, [avocado, banana, coconut])), clauses)[1] == false

# Test list appending
@test resolve(@julog(append([h,e,l,l], [o], [h,e,l,l,o])), clauses)[1] == true
@test resolve(@julog(append([o], [h,e,l,l], [h,e,l,l,o])), clauses)[1] == false

# Test list reversal
@test resolve(@julog(reverse([r,e,g,a,l], [l,a,g,e,r])), clauses)[1] == true
# A palindrome!
@test resolve(@julog(reverse([d,e,l,e,v,e,l,e,d], [d,e,l,e,v,e,l,e,d])), clauses)[1] == true
# Not a palindrome (but yes an ambigram)
@test resolve(@julog(reverse([p,a,s,s,e,d], [p,a,s,s,e,d])), clauses)[1] == false


text:
struct SubspaceDogleg <: TRSPSolver
end


text:
# Types
# =====

const BIGBED_MAGIC = 0x8789F2EB
const BIGWIG_MAGIC = 0x888FFC26
const BIGBED_BTREE_MAGIC = 0x78CA8C91
const BIGBED_RTREE_MAGIC = 0x2468ACE0

const BIGBED_MAX_ZOOM_LEVELS = 10
const BIGBED_CURRENT_VERSION = 4

const BIGWIG_DATATYPE_BEDGRAPH  = UInt8(1)
const BIGWIG_DATATYPE_VARSTEP   = UInt8(2)
const BIGWIG_DATATYPE_FIXEDSTEP = UInt8(3)

# See Supplemental Table 5
immutable BigBedHeader
    magic::UInt32
    version::UInt16
    zoom_levels::UInt16
    chromosome_tree_offset::UInt64
    full_data_offset::UInt64
    full_index_offset::UInt64
    field_count::UInt16
    defined_field_count::UInt16
    auto_sql_offset::UInt64
    total_summary_offset::UInt64
    uncompress_buf_size::UInt32
    reserved::UInt64
end

function Base.read(io::IO, ::Type{BigBedHeader})
    return BigBedHeader(
        read(io, UInt32), read(io, UInt16), read(io, UInt16),
        read(io, UInt64), read(io, UInt64), read(io, UInt64),
        read(io, UInt16), read(io, UInt16), read(io, UInt64),
        read(io, UInt64), read(io, UInt32), read(io, UInt64))
end

function Base.write(io::IO, header::BigBedHeader)
    write(io, header.magic)
    write(io, header.version)
    write(io, header.zoom_levels)
    write(io, header.chromosome_tree_offset)
    write(io, header.full_data_offset)
    write(io, header.full_index_offset)
    write(io, header.field_count)
    write(io, header.defined_field_count)
    write(io, header.auto_sql_offset)
    write(io, header.total_summary_offset)
    write(io, header.uncompress_buf_size)
    write(io, header.reserved)
end

# See Supplemental Table 6
immutable BigBedZoomHeader
    reduction_level::UInt32
    reserved::UInt32
    data_offset::UInt64
    index_offset::UInt64
end

function Base.read(io::IO, ::Type{BigBedZoomHeader})
    return BigBedZoomHeader(
        read(io, UInt32), read(io, UInt32), read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, header::BigBedZoomHeader)
    write(io, header.reduction_level)
    write(io, header.reserved)
    write(io, header.data_offset)
    write(io, header.index_offset)
end

# Supplemental Table 7
immutable BigBedTotalSummary
    bases_covered::UInt64
    min_val::Float64
    max_val::Float64
    sum_data::Float64
    sum_squares::Float64
end

function Base.read(io::IO, ::Type{BigBedTotalSummary})
    return BigBedTotalSummary(
        read(io, UInt64), read(io, UInt64), read(io, UInt64),
        read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, summary::BigBedTotalSummary)
    write(io, summary.bases_covered)
    write(io, summary.min_val)
    write(io, summary.max_val)
    write(io, summary.sum_data)
    write(io, summary.sum_squares)
end

# Supplemental Table 19
immutable BigBedZoomData
    chrom_id::UInt32
    chrom_start::UInt32
    chrom_end::UInt32
    valid_count::UInt32
    min_val::Float32
    max_val::Float32
    sum_data::Float32
    sum_squares::Float32
end

function Base.read(io::IO, ::Type{BigBedZoomData})
    return BigBedZoomData(
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, Float32), read(io, Float32), read(Float32), read(Float32))
end

function Base.write(io::IO, data::BigBedZoomData)
    write(io, data.chrom_id)
    write(io, data.chrom_start)
    write(io, data.chrom_end)
    write(io, data.valid_count)
    write(io, data.min_val)
    write(io, data.max_val)
    write(io, data.sum_data)
    write(io, data.sum_squares)
end

# Supplemental Table 8
immutable BigBedBTreeHeader
    magic::UInt32
    block_size::UInt32
    key_size::UInt32
    val_size::UInt32
    item_count::UInt64
    reserved::UInt64
end

function Base.read(io::IO, ::Type{BigBedBTreeHeader})
    return BigBedBTreeHeader(
        read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, UInt32), read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, header::BigBedBTreeHeader)
    write(io, header.magic, header.block_size, header.key_size,
          header.val_size, header.item_count, header.reserved)
end

# Supplementary Table 9
immutable BigBedBTreeNode
    isleaf::UInt8
    reserved::UInt8
    count::UInt16
end

function Base.read(io::IO, ::Type{BigBedBTreeNode})
    return BigBedBTreeNode(read(io, UInt8), read(io, UInt8), read(io, UInt16))
end

function Base.write(io::IO, node::BigBedBTreeNode)
    write(io, node.isleaf)
    write(io, node.reserved)
    write(io, node.count)
end

# Supplemental Table 10
type BigBedBTreeLeafNode
    key::Vector{UInt8}
    chrom_id::UInt32
    chrom_size::UInt32
end

function BigBedBTreeLeafNode(keysize::Integer)
    BigBedBTreeLeafNode(Array(UInt8, keysize), 0, 0)
end

function Base.read!(io::IO, node::BigBedBTreeLeafNode)
    nb = readbytes!(io, node.key, length(node.key))
    if nb < length(node.key)
        error("Unexpected end of input.")
    end
    node.chrom_id = read(io, UInt32)
    node.chrom_size = read(io, UInt32)
    return nb + 8
end

# Supplemental Table 11
type BigBedBTreeInternalNode
    key::Vector{UInt8}
    child_offset::UInt64
end

function BigBedBTreeInternalNode(keysize::Integer)
    BigBedBTreeInternalNode(Array(UInt8, keysize), 0)
end

function Base.read!(io::IO, node::BigBedBTreeInternalNode)
    if readbytes!(io, node.key, length(node.key)) < length(node.key)
        error("Unexpected end of input.")
    end
    node.child_offset = read(io, UInt64)
end

# Supplemental Table 13
immutable BigWigSectionHeader
    chrom_id::UInt32
    chrom_start::UInt32
    chrom_end::UInt32
    item_step::UInt32
    item_span::UInt32
    data_type::UInt8
    reserved::UInt8
    item_count::UInt16
end

function Base.write(io::IO, header::BigWigSectionHeader)
    write(io, header.chrom_id)
    write(io, header.chrom_start)
    write(io, header.chrom_end)
    write(io, header.item_step)
    write(io, header.item_span)
    write(io, header.data_type)
    write(io, header.reserved)
    write(io, header.item_count)
end

immutable BedGraphItem
    chrom_start::UInt32
    chrom_end::UInt32
    val::Float32
end

function Base.write(io::IO, item::BedGraphItem)
    write(io, item.chrom_start)
    write(io, item.chrom_end)
    write(io, item.val)
end

# Supplemental Table 14
immutable BigBedRTreeHeader
    magic::UInt32
    block_size::UInt32
    item_count::UInt64
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    end_file_offset::UInt64
    items_per_slot::UInt32
    reserved::UInt32
end

function Base.read(io::IO, ::Type{BigBedRTreeHeader})
    return BigBedRTreeHeader(
        read(io, UInt32), read(io, UInt32), read(io, UInt64), read(io, UInt32),
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt64),
        read(io, UInt32), read(io, UInt32))
end

function Base.write(io::IO, header::BigBedRTreeHeader)
    write(io, header.magic)
    write(io, header.block_size)
    write(io, header.item_count)
    write(io, header.start_chrom_ix)
    write(io, header.start_base)
    write(io, header.end_chrom_ix)
    write(io, header.end_base)
    write(io, header.end_file_offset)
    write(io, header.items_per_slot)
    write(io, header.reserved)
end

# Supplemental Table 15
immutable BigBedRTreeNode
    isleaf::UInt8
    reserved::UInt8
    count::UInt16
end

function Base.read(io::IO, ::Type{BigBedRTreeNode})
    return BigBedRTreeNode(read(io, UInt8), read(io, UInt8), read(io, UInt16))
end

function Base.write(io::IO, node::BigBedRTreeNode)
    write(io, node.isleaf)
    write(io, node.reserved)
    write(io, node.count)
end

# Supplemental Table 16
immutable BigBedRTreeLeafNode
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    data_offset::UInt64
    data_size::UInt64
end

function Base.read(io::IO, ::Type{BigBedRTreeLeafNode})
    return BigBedRTreeLeafNode(
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, UInt64), read(io, UInt64))
end

function Base.write(io::IO, node::BigBedRTreeLeafNode)
    write(io, node.start_chrom_ix)
    write(io, node.start_base)
    write(io, node.end_chrom_ix)
    write(io, node.end_base)
    write(io, node.data_offset)
    write(io, node.data_size)
end

# Supplemental Table 17
immutable BigBedRTreeInternalNode
    start_chrom_ix::UInt32
    start_base::UInt32
    end_chrom_ix::UInt32
    end_base::UInt32
    data_offset::UInt64
end

function Base.read(io::IO, ::Type{BigBedRTreeInternalNode})
    return BigBedRTreeInternalNode(
        read(io, UInt32), read(io, UInt32), read(io, UInt32), read(io, UInt32),
        read(io, UInt64))
end

function Base.write(io::IO, node::BigBedRTreeInternalNode)
    write(io, node.start_chrom_ix)
    write(io, node.start_base)
    write(io, node.end_chrom_ix)
    write(io, node.end_base)
    write(io, node.data_offset)
end


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

@test integrate(x^3*(a+b*x^2)^(1/2)*(A+B*x), x) == :(-1//120*a*b^-2*(a+b*x^2)^3//2*(16A+15*B*x)+1//5*A*b^-1*x^2*(a+b*x^2)^3//2+1//6*B*b^-1*x^3*(a+b*x^2)^3//2+1//16*B*a^3*b^-5//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//16*B*x*a^2*b^-2*(a+b*x^2)^(1/2))
@test integrate(x^2*(a+b*x^2)^(1/2)*(A+B*x), x) == :(-1//60*b^-2*(a+b*x^2)^3//2*(8*B*a+-15*A*b*x)+-1//8*A*a^2*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//5*B*b^-1*x^2*(a+b*x^2)^3//2+-1//8*A*a*x*b^-1*(a+b*x^2)^(1/2))
@test integrate(x*(a+b*x^2)^(1/2)*(A+B*x), x) == :(1//12*b^-1*(a+b*x^2)^3//2*(4A+3*B*x)+-1//8*B*a^2*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1//8*B*a*x*b^-1*(a+b*x^2)^(1/2))
@test integrate((a+b*x^2)^(1/2)*(A+B*x), x) == :((1/2)*A*x*(a+b*x^2)^(1/2)+1//3*B*b^-1*(a+b*x^2)^3//2+(1/2)*A*a*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^-1*(a+b*x^2)^(1/2)*(A+B*x), x) == :((1/2)*(a+b*x^2)^(1/2)*(2A+B*x)+-1*A*a^(1/2)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+(1/2)*B*a*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^-2*(a+b*x^2)^(1/2)*(A+B*x), x) == :(A*b^(1/2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1*B*a^(1/2)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1*x^-1*(a+b*x^2)^(1/2)*(A+-1*B*x))
@test integrate(x^-3*(a+b*x^2)^(1/2)*(A+B*x), x) == :(B*b^(1/2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1//2*x^-2*(a+b*x^2)^(1/2)*(A+2*B*x)+-1//2*A*b*a^-1//2*arctanh(a^-1//2*(a+b*x^2)^(1/2)))
@test integrate(x^3*(a+b*x^2)^3//2*(A+B*x), x) == :(-1//560*a*b^-2*(a+b*x^2)^5//2*(32A+35*B*x)+1//7*A*b^-1*x^2*(a+b*x^2)^5//2+1//8*B*b^-1*x^3*(a+b*x^2)^5//2+3//128*B*a^4*b^-5//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//64*B*x*a^2*b^-2*(a+b*x^2)^3//2+3//128*B*x*a^3*b^-2*(a+b*x^2)^(1/2))
@test integrate(x^2*(a+b*x^2)^3//2*(A+B*x), x) == :(-1//210*b^-2*(a+b*x^2)^5//2*(12*B*a+-35*A*b*x)+-1//16*A*a^3*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//7*B*b^-1*x^2*(a+b*x^2)^5//2+-1//16*A*x*a^2*b^-1*(a+b*x^2)^(1/2)+-1//24*A*a*x*b^-1*(a+b*x^2)^3//2)
@test integrate(x*(a+b*x^2)^3//2*(A+B*x), x) == :(1//30*b^-1*(a+b*x^2)^5//2*(6A+5*B*x)+-1//16*B*a^3*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1//16*B*x*a^2*b^-1*(a+b*x^2)^(1/2)+-1//24*B*a*x*b^-1*(a+b*x^2)^3//2)
@test integrate((a+b*x^2)^3//2*(A+B*x), x) == :(1//4*A*x*(a+b*x^2)^3//2+1//5*B*b^-1*(a+b*x^2)^5//2+3//8*A*a*x*(a+b*x^2)^(1/2)+3//8*A*a^2*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^-1*(a+b*x^2)^3//2*(A+B*x), x) == :(1//12*(a+b*x^2)^3//2*(4A+3*B*x)+-1*A*a^3//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+1//8*a*(a+b*x^2)^(1/2)*(8A+3*B*x)+3//8*B*a^2*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^-2*(a+b*x^2)^3//2*(A+B*x), x) == :((1/2)*(a+b*x^2)^(1/2)*(2*B*a+3*A*b*x)+-1*B*a^3//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1//3*x^-1*(a+b*x^2)^3//2*(3A+-1*B*x)+3//2*A*a*b^(1/2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^-3*(a+b*x^2)^3//2*(A+B*x), x) == :((1/2)*x^-1*(a+b*x^2)^(1/2)*(-3*B*a+3*A*b*x)+-1//2*x^-2*(a+b*x^2)^3//2*(A+-1*B*x)+-3//2*A*b*a^(1/2)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+3//2*B*a*b^(1/2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^3*(a+b*x^2)^5//2*(A+B*x), x) == :(-1//5040*a*b^-2*(a+b*x^2)^7//2*(160A+189*B*x)+1//9*A*b^-1*x^2*(a+b*x^2)^7//2+1//10*B*b^-1*x^3*(a+b*x^2)^7//2+3//256*B*a^5*b^-5//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//128*B*x*a^3*b^-2*(a+b*x^2)^3//2+1//160*B*x*a^2*b^-2*(a+b*x^2)^5//2+3//256*B*x*a^4*b^-2*(a+b*x^2)^(1/2))
@test integrate(x^2*(a+b*x^2)^5//2*(A+B*x), x) == :(-1//504*b^-2*(a+b*x^2)^7//2*(16*B*a+-63*A*b*x)+-5//128*A*a^4*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//9*B*b^-1*x^2*(a+b*x^2)^7//2+-5//128*A*x*a^3*b^-1*(a+b*x^2)^(1/2)+-5//192*A*x*a^2*b^-1*(a+b*x^2)^3//2+-1//48*A*a*x*b^-1*(a+b*x^2)^5//2)
@test integrate(x*(a+b*x^2)^5//2*(A+B*x), x) == :(1//56*b^-1*(a+b*x^2)^7//2*(8A+7*B*x)+-5//128*B*a^4*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-5//128*B*x*a^3*b^-1*(a+b*x^2)^(1/2)+-5//192*B*x*a^2*b^-1*(a+b*x^2)^3//2+-1//48*B*a*x*b^-1*(a+b*x^2)^5//2)
@test integrate((a+b*x^2)^5//2*(A+B*x), x) == :(1//6*A*x*(a+b*x^2)^5//2+1//7*B*b^-1*(a+b*x^2)^7//2+5//16*A*x*a^2*(a+b*x^2)^(1/2)+5//16*A*a^3*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+5//24*A*a*x*(a+b*x^2)^3//2)
@test integrate(x^-1*(a+b*x^2)^5//2*(A+B*x), x) == :(1//30*(a+b*x^2)^5//2*(6A+5*B*x)+-1*A*a^5//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+1//16*a^2*(a+b*x^2)^(1/2)*(16A+5*B*x)+1//24*a*(a+b*x^2)^3//2*(8A+5*B*x)+5//16*B*a^3*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^-2*(a+b*x^2)^5//2*(A+B*x), x) == :(1//12*(a+b*x^2)^3//2*(4*B*a+15*A*b*x)+-1*B*a^5//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1//5*x^-1*(a+b*x^2)^5//2*(5A+-1*B*x)+1//8*a*(a+b*x^2)^(1/2)*(8*B*a+15*A*b*x)+15//8*A*a^2*b^(1/2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^-3*(a+b*x^2)^5//2*(A+B*x), x) == :(-1//4*x^-2*(a+b*x^2)^5//2*(2A+-1*B*x)+-1//12*x^-1*(a+b*x^2)^3//2*(15*B*a+-10*A*b*x)+-5//2*A*b*a^3//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+5//8*a*b*(a+b*x^2)^(1/2)*(4A+3*B*x)+15//8*B*a^2*b^(1/2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^3*(a+b*x^2)^-1//2*(A+B*x), x) == :(-1//24*a*b^-2*(a+b*x^2)^(1/2)*(16A+9*B*x)+1//3*A*b^-1*x^2*(a+b*x^2)^(1/2)+1//4*B*b^-1*x^3*(a+b*x^2)^(1/2)+3//8*B*a^2*b^-5//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^2*(a+b*x^2)^-1//2*(A+B*x), x) == :(-1//6*b^-2*(a+b*x^2)^(1/2)*(4*B*a+-3*A*b*x)+-1//2*A*a*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//3*B*b^-1*x^2*(a+b*x^2)^(1/2))
@test integrate(x*(a+b*x^2)^-1//2*(A+B*x), x) == :((1/2)*b^-1*(a+b*x^2)^(1/2)*(2A+B*x)+-1//2*B*a*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate((a+b*x^2)^-1//2*(A+B*x), x) == :(A*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+B*b^-1*(a+b*x^2)^(1/2))
@test integrate(x^-1*(a+b*x^2)^-1//2*(A+B*x), x) == :(B*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1*A*a^-1//2*arctanh(a^-1//2*(a+b*x^2)^(1/2)))
@test integrate(x^-2*(a+b*x^2)^-1//2*(A+B*x), x) == :(-1*B*a^-1//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1*A*a^-1*x^-1*(a+b*x^2)^(1/2))
@test integrate(x^-3*(a+b*x^2)^-1//2*(A+B*x), x) == :((1/2)*A*b*a^-3//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1*B*a^-1*x^-1*(a+b*x^2)^(1/2)+-1//2*A*a^-1*x^-2*(a+b*x^2)^(1/2))
@test integrate(x^3*(a+b*x^2)^-3//2*(A+B*x), x) == :((1/2)*b^-2*(a+b*x^2)^(1/2)*(4A+3*B*x)+-1*b^-1*x^2*(a+b*x^2)^-1//2*(A+B*x)+-3//2*B*a*b^-5//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2))
@test integrate(x^2*(a+b*x^2)^-3//2*(A+B*x), x) == :(A*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+2*B*b^-2*(a+b*x^2)^(1/2)+-1*x*b^-1*(a+b*x^2)^-1//2*(A+B*x))
@test integrate(x*(a+b*x^2)^-3//2*(A+B*x), x) == :(B*b^-3//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1*b^-1*(a+b*x^2)^-1//2*(A+B*x))
@test integrate((a+b*x^2)^-3//2*(A+B*x), x) == :(-1*a^-1*b^-1*(a+b*x^2)^-1//2*(B*a+-1*A*b*x))
@test integrate(x^-1*(a+b*x^2)^-3//2*(A+B*x), x) == :(a^-1*(a+b*x^2)^-1//2*(A+B*x)+-1*A*a^-3//2*arctanh(a^-1//2*(a+b*x^2)^(1/2)))
@test integrate(x^-2*(a+b*x^2)^-3//2*(A+B*x), x) == :(-1*B*a^-3//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+a^-1*x^-1*(a+b*x^2)^-1//2*(A+B*x)+-2*A*a^-2*x^-1*(a+b*x^2)^(1/2))
@test integrate(x^-3*(a+b*x^2)^-3//2*(A+B*x), x) == :(a^-1*x^-2*(a+b*x^2)^-1//2*(A+B*x)+-2*B*a^-2*x^-1*(a+b*x^2)^(1/2)+-3//2*A*a^-2*x^-2*(a+b*x^2)^(1/2)+3//2*A*b*a^-5//2*arctanh(a^-1//2*(a+b*x^2)^(1/2)))
@test integrate(x^3*(a+b*x^2)^-5//2*(A+B*x), x) == :(B*b^-5//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1//3*b^-2*(a+b*x^2)^-1//2*(2A+3*B*x)+-1//3*b^-1*x^2*(a+b*x^2)^-3//2*(A+B*x))
@test integrate(x^2*(a+b*x^2)^-5//2*(A+B*x), x) == :(-2//3*B*b^-2*(a+b*x^2)^-1//2+-1//3*a^-1*b^-1*x^2*(a+b*x^2)^-3//2*(B*a+-1*A*b*x))
@test integrate(x*(a+b*x^2)^-5//2*(A+B*x), x) == :(1//3*b^-1*(a+b*x^2)^-3//2*(-1A+-1*B*x)+1//3*B*x*a^-1*b^-1*(a+b*x^2)^-1//2)
@test integrate((a+b*x^2)^-5//2*(A+B*x), x) == :(1//3*a^-1*b^-1*(a+b*x^2)^-3//2*(-1*B*a+A*b*x)+2//3*A*x*a^-2*(a+b*x^2)^-1//2)
@test integrate(x^-1*(a+b*x^2)^-5//2*(A+B*x), x) == :(-1*A*a^-5//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+1//3*a^-1*(a+b*x^2)^-3//2*(A+B*x)+1//3*a^-2*(a+b*x^2)^-1//2*(3A+2*B*x))
@test integrate(x^-2*(a+b*x^2)^-5//2*(A+B*x), x) == :(-1*B*a^-5//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-8//3*A*a^-3*x^-1*(a+b*x^2)^(1/2)+1//3*a^-1*x^-1*(a+b*x^2)^-3//2*(A+B*x)+1//3*a^-2*x^-1*(a+b*x^2)^-1//2*(4A+3*B*x))
@test integrate(x^-3*(a+b*x^2)^-5//2*(A+B*x), x) == :(-8//3*B*a^-3*x^-1*(a+b*x^2)^(1/2)+-5//2*A*a^-3*x^-2*(a+b*x^2)^(1/2)+1//3*a^-1*x^-2*(a+b*x^2)^-3//2*(A+B*x)+1//3*a^-2*x^-2*(a+b*x^2)^-1//2*(5A+4*B*x)+5//2*A*b*a^-7//2*arctanh(a^-1//2*(a+b*x^2)^(1/2)))
@test integrate(x*(1+-1*x^2)^-1//2*(1+-1x), x) == :(-1//2*arcsin(x)+(1+-1*x^2)^(1/2)*(-1+(1/2)*x))
@test integrate((1+-1*x^2)^-1//2*(x+-1*x^2), x) == :(-1//2*arcsin(x)+(1+-1*x^2)^(1/2)*(-1+(1/2)*x))
@test integrate((-3+x^2)^-1*(3+x^2), x) == :(x+-2*3^(1/2)*arctanh(1//3*x*3^(1/2)))
@test integrate((1+x^2)^-1*(-1+x^2), x) == :(x+-2*arctan(x))
@test integrate(x^7*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(B*b^-9//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1//35*a^-1*b^-5*(a+b*x^2)^(1/2)*(-128*C*a+16*A*b)+-1//7*a^-1*b^-1*x^7*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a))+-1//35*x*a^-1*b^-4*(a+b*x^2)^-1//2*(-1*x*(-64*C*a+8*A*b)+35*B*a)+-1//35*a^-1*b^-2*x^5*(a+b*x^2)^-5//2*(-1*x*(A*b+-8*C*a)+7*B*a)+-1//105*a^-1*b^-3*x^3*(a+b*x^2)^-3//2*(-1*x*(-48*C*a+6*A*b)+35*B*a))
@test integrate(x^6*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(C*b^-9//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1//35*b^-4*(a+b*x^2)^-1//2*(16B+35*C*x)+-1//35*b^-2*x^4*(a+b*x^2)^-5//2*(6B+7*C*x)+-1//105*b^-3*x^2*(a+b*x^2)^-3//2*(24B+35*C*x)+-1//7*a^-1*b^-1*x^6*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a)))
@test integrate(x^5*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(1//105*b^-4*(a+b*x^2)^-3//2*(4*A*b+24*C*a)+-1//35*a^-1*b^-4*(a+b*x^2)^-1//2*(4*A*b+24*C*a)+-1//7*a^-1*b^-1*x^5*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a))+-1//35*a^-1*b^-2*x^4*(a+b*x^2)^-5//2*(A*b+6*C*a+-5*B*b*x))
@test integrate(x^4*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(-1//105*a^-1*b^-3*(a+b*x^2)^-3//2*(x*(6*A*b+15*C*a)+8*B*a)+-1//7*a^-1*b^-1*x^4*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a))+-1//35*a^-1*b^-2*x^2*(a+b*x^2)^-5//2*(x*(2*A*b+5*C*a)+4*B*a)+1//35*x*a^-2*b^-3*(a+b*x^2)^-1//2*(2*A*b+5*C*a))
@test integrate(x^3*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(-1//105*a^-1*b^-3*(a+b*x^2)^-3//2*(6*A*b+8*C*a+-3*B*b*x)+-1//7*a^-1*b^-1*x^3*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a))+-1//35*x*a^-1*b^-2*(a+b*x^2)^-5//2*(x*(3*A*b+4*C*a)+3*B*a)+2//35*B*x*a^-2*b^-2*(a+b*x^2)^-1//2)
@test integrate(x^2*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(-1//35*a^-1*b^-2*(a+b*x^2)^-5//2*(x*(3*C*a+4*A*b)+2*B*a)+-1//7*a^-1*b^-1*x^2*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a))+1//105*x*a^-3*b^-2*(a+b*x^2)^-1//2*(6*C*a+8*A*b)+1//105*x*a^-2*b^-2*(a+b*x^2)^-3//2*(3*C*a+4*A*b))
@test integrate(x*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(-1//35*a^-1*b^-2*(a+b*x^2)^-5//2*(2*C*a+5*A*b+-1*B*b*x)+-1//7*x*a^-1*b^-1*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a))+4//105*B*x*a^-2*b^-1*(a+b*x^2)^-3//2+8//105*B*x*a^-3*b^-1*(a+b*x^2)^-1//2)
@test integrate((a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(-1//7*a^-1*b^-1*(a+b*x^2)^-7//2*(B*a+-1*x*(A*b+-1*C*a))+1//35*x*a^-2*b^-1*(a+b*x^2)^-5//2*(C*a+6*A*b)+1//105*x*a^-4*b^-1*(a+b*x^2)^-1//2*(8*C*a+48*A*b)+1//105*x*a^-3*b^-1*(a+b*x^2)^-3//2*(4*C*a+24*A*b))
@test integrate(x^-1*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(-1*A*a^-9//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+1//35*a^-4*(a+b*x^2)^-1//2*(35A+16*B*x)+1//35*a^-2*(a+b*x^2)^-5//2*(7A+6*B*x)+1//105*a^-3*(a+b*x^2)^-3//2*(35A+24*B*x)+1//7*a^-1*b^-1*(a+b*x^2)^-7//2*(A*b+-1*C*a+B*b*x))
@test integrate(x^-2*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :(-1*B*a^-9//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+1//7*a^-1*(a+b*x^2)^-7//2*(B+-1*x*(-1C+A*b*a^-1))+1//35*a^-4*(a+b*x^2)^-1//2*(35B+-1*x*(-16C+93*A*b*a^-1))+1//35*a^-2*(a+b*x^2)^-5//2*(7B+-1*x*(-6C+13*A*b*a^-1))+1//105*a^-3*(a+b*x^2)^-3//2*(35B+-1*x*(-24C+87*A*b*a^-1))+-1*A*a^-5*x^-1*(a+b*x^2)^(1/2))
@test integrate(x^-3*(a+b*x^2)^-9//2*(A+B*x+C*x^2), x) == :((1/2)*a^-11//2*(-2*C*a+9*A*b)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1//7*a^-2*(a+b*x^2)^-7//2*(a*(-1C+A*b*a^-1)+B*b*x)+-1//35*a^-5*(a+b*x^2)^-1//2*(-35*C*a+140*A*b+93*B*b*x)+-1//35*a^-3*(a+b*x^2)^-5//2*(-7*C*a+14*A*b+13*B*b*x)+-1//105*a^-4*(a+b*x^2)^-3//2*(-35*C*a+105*A*b+87*B*b*x)+-1*B*a^-5*x^-1*(a+b*x^2)^(1/2)+-1//2*A*a^-5*x^-2*(a+b*x^2)^(1/2))
@test integrate(x^3*(a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :(1//6*x^6*(A*b+C*a)+1//7*x^7*(B*b+D*a)+1//4*A*a*x^4+1//5*B*a*x^5+1//8*C*b*x^8+1//9*D*b*x^9)
@test integrate(x^2*(a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :(1//5*x^5*(A*b+C*a)+1//6*x^6*(B*b+D*a)+1//3*A*a*x^3+1//4*B*a*x^4+1//7*C*b*x^7+1//8*D*b*x^8)
@test integrate(x*(a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :(1//4*x^4*(A*b+C*a)+1//5*x^5*(B*b+D*a)+(1/2)*A*a*x^2+1//3*B*a*x^3+1//6*C*b*x^6+1//7*D*b*x^7)
@test integrate((a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :(1//3*x^3*(A*b+C*a)+1//4*x^4*(B*b+D*a)+A*a*x+(1/2)*B*a*x^2+1//5*C*b*x^5+1//6*D*b*x^6)
@test integrate(x^-1*(a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*x^2*(A*b+C*a)+1//3*x^3*(B*b+D*a)+A*a*log(x)+B*a*x+1//4*C*b*x^4+1//5*D*b*x^5)
@test integrate(x^-2*(a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :(x*(A*b+C*a)+(1/2)*x^2*(B*b+D*a)+B*a*log(x)+-1*A*a*x^-1+1//3*C*b*x^3+1//4*D*b*x^4)
@test integrate(x^-3*(a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :(x*(B*b+D*a)+(A*b+C*a)*log(x)+(1/2)*C*b*x^2+-1*B*a*x^-1+-1//2*A*a*x^-2+1//3*D*b*x^3)
@test integrate(x^-4*(a+b*x^2)*(A+B*x+C*x^2+D*x^3), x) == :((B*b+D*a)*log(x)+-1*x^-1*(A*b+C*a)+C*b*x+(1/2)*D*b*x^2+-1//2*B*a*x^-2+-1//3*A*a*x^-3)
@test integrate(x^3*(a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(1//4*A*a^2*x^4+1//5*B*a^2*x^5+1//6*a*x^6*(C*a+2*A*b)+1//7*a*x^7*(D*a+2*B*b)+1//8*b*x^8*(A*b+2*C*a)+1//9*b*x^9*(B*b+2*D*a)+1//10*C*b^2*x^10+1//11*D*b^2*x^11)
@test integrate(x^2*(a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(1//3*A*a^2*x^3+1//4*B*a^2*x^4+1//5*a*x^5*(C*a+2*A*b)+1//6*a*x^6*(D*a+2*B*b)+1//7*b*x^7*(A*b+2*C*a)+1//8*b*x^8*(B*b+2*D*a)+1//9*C*b^2*x^9+1//10*D*b^2*x^10)
@test integrate(x*(a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(1//3*B*a^2*x^3+1//4*C*a^2*x^4+1//5*a*x^5*(D*a+2*B*b)+1//6*A*b^-1*(a+b*x^2)^3+1//7*b*x^7*(B*b+2*D*a)+1//8*C*b^2*x^8+1//9*D*b^2*x^9+1//3*C*a*b*x^6)
@test integrate((a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(A*x*a^2+1//3*a*x^3*(C*a+2*A*b)+1//4*D*a^2*x^4+1//5*b*x^5*(A*b+2*C*a)+1//6*B*b^-1*(a+b*x^2)^3+1//7*C*b^2*x^7+1//8*D*b^2*x^8+1//3*D*a*b*x^6)
@test integrate(x^-1*(a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(A*a^2*log(x)+B*x*a^2+1//3*a*x^3*(D*a+2*B*b)+1//4*A*b^2*x^4+1//5*b*x^5*(B*b+2*D*a)+1//6*C*b^-1*(a+b*x^2)^3+1//7*D*b^2*x^7+A*a*b*x^2)
@test integrate(x^-2*(a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(B*a^2*log(x)+a*x*(C*a+2*A*b)+-1*A*a^2*x^-1+1//3*b*x^3*(A*b+2*C*a)+1//4*B*b^2*x^4+1//5*C*b^2*x^5+1//6*D*b^-1*(a+b*x^2)^3+B*a*b*x^2)
@test integrate(x^-3*(a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(a*x*(D*a+2*B*b)+a*(C*a+2*A*b)*log(x)+(1/2)*b*x^2*(A*b+2*C*a)+-1*B*a^2*x^-1+-1//2*A*a^2*x^-2+1//3*b*x^3*(B*b+2*D*a)+1//4*C*b^2*x^4+1//5*D*b^2*x^5)
@test integrate(x^-4*(a+b*x^2)^2*(A+B*x+C*x^2+D*x^3), x) == :(a*(D*a+2*B*b)*log(x)+b*x*(A*b+2*C*a)+(1/2)*b*x^2*(B*b+2*D*a)+-1*a*x^-1*(C*a+2*A*b)+-1//2*B*a^2*x^-2+-1//3*A*a^2*x^-3+1//3*C*b^2*x^3+1//4*D*b^2*x^4)
@test integrate(x^3*(a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(1//4*A*a^3*x^4+1//5*B*a^3*x^5+1//6*a^2*x^6*(C*a+3*A*b)+1//7*a^2*x^7*(D*a+3*B*b)+1//10*b^2*x^10*(A*b+3*C*a)+1//11*b^2*x^11*(B*b+3*D*a)+1//12*C*b^3*x^12+1//13*D*b^3*x^13+1//3*a*b*x^9*(B*b+D*a)+3//8*a*b*x^8*(A*b+C*a))
@test integrate(x^2*(a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(1//3*A*a^3*x^3+1//4*B*a^3*x^4+1//5*a^2*x^5*(C*a+3*A*b)+1//6*a^2*x^6*(D*a+3*B*b)+1//9*b^2*x^9*(A*b+3*C*a)+1//10*b^2*x^10*(B*b+3*D*a)+1//11*C*b^3*x^11+1//12*D*b^3*x^12+3//7*a*b*x^7*(A*b+C*a)+3//8*a*b*x^8*(B*b+D*a))
@test integrate(x*(a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(1//3*B*a^3*x^3+1//4*C*a^3*x^4+1//5*a^2*x^5*(D*a+3*B*b)+1//8*A*b^-1*(a+b*x^2)^4+1//9*b^2*x^9*(B*b+3*D*a)+1//10*C*b^3*x^10+1//11*D*b^3*x^11+(1/2)*C*b*a^2*x^6+3//7*a*b*x^7*(B*b+D*a)+3//8*C*a*b^2*x^8)
@test integrate((a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(A*x*a^3+1//3*a^2*x^3*(C*a+3*A*b)+1//4*D*a^3*x^4+1//7*b^2*x^7*(A*b+3*C*a)+1//8*B*b^-1*(a+b*x^2)^4+1//9*C*b^3*x^9+1//10*D*b^3*x^10+(1/2)*D*b*a^2*x^6+3//5*a*b*x^5*(A*b+C*a)+3//8*D*a*b^2*x^8)
@test integrate(x^-1*(a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(A*a^3*log(x)+B*x*a^3+1//3*a^2*x^3*(D*a+3*B*b)+1//6*A*b^3*x^6+1//7*b^2*x^7*(B*b+3*D*a)+1//8*C*b^-1*(a+b*x^2)^4+1//9*D*b^3*x^9+3//2*A*b*a^2*x^2+3//4*A*a*b^2*x^4+3//5*a*b*x^5*(B*b+D*a))
@test integrate(x^-2*(a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(B*a^3*log(x)+x*a^2*(C*a+3*A*b)+-1*A*a^3*x^-1+1//5*b^2*x^5*(A*b+3*C*a)+1//6*B*b^3*x^6+1//7*C*b^3*x^7+1//8*D*b^-1*(a+b*x^2)^4+a*b*x^3*(A*b+C*a)+3//2*B*b*a^2*x^2+3//4*B*a*b^2*x^4)
@test integrate(x^-3*(a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(x*a^2*(D*a+3*B*b)+a^2*(C*a+3*A*b)*log(x)+-1*B*a^3*x^-1+-1//2*A*a^3*x^-2+1//4*b^2*x^4*(A*b+3*C*a)+1//5*b^2*x^5*(B*b+3*D*a)+1//6*C*b^3*x^6+1//7*D*b^3*x^7+a*b*x^3*(B*b+D*a)+3//2*a*b*x^2*(A*b+C*a))
@test integrate(x^-4*(a+b*x^2)^3*(A+B*x+C*x^2+D*x^3), x) == :(a^2*(D*a+3*B*b)*log(x)+-1*a^2*x^-1*(C*a+3*A*b)+-1//2*B*a^3*x^-2+-1//3*A*a^3*x^-3+1//3*b^2*x^3*(A*b+3*C*a)+1//4*b^2*x^4*(B*b+3*D*a)+1//5*C*b^3*x^5+1//6*D*b^3*x^6+3*a*b*x*(A*b+C*a)+3//2*a*b*x^2*(B*b+D*a))
@test integrate(x^4*(a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :(1//3*b^-2*x^3*(A*b+-1*C*a)+1//4*b^-2*x^4*(B*b+-1*D*a)+1//5*C*b^-1*x^5+1//6*D*b^-1*x^6+a^3//2*b^-7//2*(A*b+-1*C*a)*arctan(x*a^-1//2*b^(1/2))+(1/2)*a^2*b^-4*(B*b+-1*D*a)*log(a+b*x^2)+-1*a*x*b^-3*(A*b+-1*C*a)+-1//2*a*b^-3*x^2*(B*b+-1*D*a))
@test integrate(x^3*(a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*b^-2*x^2*(A*b+-1*C*a)+1//3*b^-2*x^3*(B*b+-1*D*a)+1//4*C*b^-1*x^4+1//5*D*b^-1*x^5+a^3//2*b^-7//2*(B*b+-1*D*a)*arctan(x*a^-1//2*b^(1/2))+-1*a*x*b^-3*(B*b+-1*D*a)+-1//2*a*b^-3*(A*b+-1*C*a)*log(a+b*x^2))
@test integrate(x^2*(a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :(x*b^-2*(A*b+-1*C*a)+(1/2)*b^-2*x^2*(B*b+-1*D*a)+1//3*C*b^-1*x^3+1//4*D*b^-1*x^4+-1*a^(1/2)*b^-5//2*(A*b+-1*C*a)*arctan(x*a^-1//2*b^(1/2))+-1//2*a*b^-3*(B*b+-1*D*a)*log(a+b*x^2))
@test integrate(x*(a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :(x*b^-2*(B*b+-1*D*a)+(1/2)*C*b^-1*x^2+(1/2)*b^-2*(A*b+-1*C*a)*log(a+b*x^2)+1//3*D*b^-1*x^3+-1*a^(1/2)*b^-5//2*(B*b+-1*D*a)*arctan(x*a^-1//2*b^(1/2)))
@test integrate((a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :(C*x*b^-1+(1/2)*D*b^-1*x^2+(1/2)*b^-2*(B*b+-1*D*a)*log(a+b*x^2)+a^-1//2*b^-3//2*(A*b+-1*C*a)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-1*(a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :(A*a^-1*log(x)+D*x*b^-1+a^-1//2*b^-3//2*(B*b+-1*D*a)*arctan(x*a^-1//2*b^(1/2))+-1//2*a^-1*b^-1*(A*b+-1*C*a)*log(a+b*x^2))
@test integrate(x^-2*(a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :(B*a^-1*log(x)+-1*A*a^-1*x^-1+-1*a^-3//2*b^-1//2*(A*b+-1*C*a)*arctan(x*a^-1//2*b^(1/2))+-1//2*a^-1*b^-1*(B*b+-1*D*a)*log(a+b*x^2))
@test integrate(x^-3*(a+b*x^2)^-1*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*a^-2*(A*b+-1*C*a)*log(a+b*x^2)+-1*B*a^-1*x^-1+-1*a^-2*(A*b+-1*C*a)*log(x)+-1//2*A*a^-1*x^-2+-1*a^-3//2*b^-1//2*(B*b+-1*D*a)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^4*(a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*x*b^-3*(-5*C*a+3*A*b)+(1/2)*b^-3*x^2*(-3*D*a+2*B*b)+1//4*D*b^-2*x^4+-1//2*a*b^-4*(-3*D*a+2*B*b)*log(a+b*x^2)+-1//2*a^(1/2)*b^-7//2*(-5*C*a+3*A*b)*arctan(x*a^-1//2*b^(1/2))+-1//6*a^-1*b^-2*x^3*(-5*C*a+3*A*b)+-1//2*a^-1*b^-1*x^4*(a+b*x^2)^-1*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a)))
@test integrate(x^3*(a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*x*b^-3*(-5*D*a+3*B*b)+(1/2)*b^-3*(A*b+-2*C*a)*log(a+b*x^2)+1//3*D*b^-2*x^3+-1//2*a^(1/2)*b^-7//2*(-5*D*a+3*B*b)*arctan(x*a^-1//2*b^(1/2))+-1//2*a^-1*b^-2*x^2*(A*b+-2*C*a)+-1//2*a^-1*b^-1*x^3*(a+b*x^2)^-1*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a)))
@test integrate(x^2*(a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*D*b^-2*x^2+(1/2)*b^-3*(B*b+-2*D*a)*log(a+b*x^2)+(1/2)*a^-1//2*b^-5//2*(A*b+-3*C*a)*arctan(x*a^-1//2*b^(1/2))+-1//2*x*a^-1*b^-2*(A*b+-3*C*a)+-1//2*a^-1*b^-1*x^2*(a+b*x^2)^-1*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a)))
@test integrate(x*(a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :(D*x*b^-2+(1/2)*C*b^-2*log(a+b*x^2)+(1/2)*a^-1//2*b^-5//2*(B*b+-3*D*a)*arctan(x*a^-1//2*b^(1/2))+-1//2*x*a^-1*b^-1*(a+b*x^2)^-1*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a)))
@test integrate((a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*D*b^-2*log(a+b*x^2)+(1/2)*a^-1*b^-1*(a+b*x^2)^-1*(x*(A*b+-1*C*a)+-1*a*(B+-1*D*a*b^-1))+(1/2)*a^-3//2*b^-3//2*(A*b+C*a)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-1*(a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :(A*a^-2*log(x)+-1//2*A*a^-2*log(a+b*x^2)+(1/2)*a^-1*b^-1*(a+b*x^2)^-1*(A*b+x*(B*b+-1*D*a)+-1*C*a)+(1/2)*a^-3//2*b^-3//2*(B*b+D*a)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-2*(a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :(B*a^-2*log(x)+-1*A*a^-2*x^-1+-1//2*B*a^-2*log(a+b*x^2)+(1/2)*a^-1*b^-1*(a+b*x^2)^-1*(B*b+-1*D*a+-1*b*x*(-1C+A*b*a^-1))+-1//2*a^-5//2*b^-1//2*(-1*C*a+3*A*b)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-3*(a+b*x^2)^-2*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*a^-3*(-1*C*a+2*A*b)*log(a+b*x^2)+-1*B*a^-2*x^-1+-1*a^-3*(-1*C*a+2*A*b)*log(x)+-1//2*A*a^-2*x^-2+-1//2*a^-1*(a+b*x^2)^-1*(-1C+x*(-1D+B*b*a^-1)+A*b*a^-1)+-1//2*a^-5//2*b^-1//2*(-1*D*a+3*B*b)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^4*(a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*b^-4*(B*b+-3*D*a)*log(a+b*x^2)+-1//2*a^-1*b^-3*x^2*(B*b+-3*D*a)+1//8*x*a^-1*b^-3*(-3*A*b+15*C*a)+1//8*a^-1//2*b^-7//2*(-15*C*a+3*A*b)*arctan(x*a^-1//2*b^(1/2))+-1//4*a^-1*b^-1*x^4*(a+b*x^2)^-2*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a))+1//8*a^-1*b^-2*x^3*(a+b*x^2)^-1*(A*b+x*(-8*D*a+4*B*b)+-5*C*a))
@test integrate(x^3*(a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*C*b^-3*log(a+b*x^2)+1//8*x*a^-1*b^-3*(-3*B*b+15*D*a)+1//8*a^-1//2*b^-7//2*(-15*D*a+3*B*b)*arctan(x*a^-1//2*b^(1/2))+-1//4*a^-1*b^-1*x^3*(a+b*x^2)^-2*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a))+-1//8*a^-1*b^-2*x^2*(a+b*x^2)^-1*(-1*x*(-7*D*a+3*B*b)+4*C*a))
@test integrate(x^2*(a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*D*b^-3*log(a+b*x^2)+1//8*a^-3//2*b^-5//2*(A*b+3*C*a)*arctan(x*a^-1//2*b^(1/2))+-1//4*a^-1*b^-1*x^2*(a+b*x^2)^-2*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a))+-1//8*x*a^-1*b^-2*(a+b*x^2)^-1*(A*b+-1*x*(-6*D*a+2*B*b)+3*C*a))
@test integrate(x*(a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :(-1//8*a^-1*b^-2*(a+b*x^2)^-1*(-1*x*(B*b+-5*D*a)+2*A*b+2*C*a)+1//8*a^-3//2*b^-5//2*(B*b+3*D*a)*arctan(x*a^-1//2*b^(1/2))+-1//4*x*a^-1*b^-1*(a+b*x^2)^-2*(a*(B+-1*D*a*b^-1)+-1*x*(A*b+-1*C*a)))
@test integrate((a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :(-1//8*a^-2*b^-2*(a+b*x^2)^-1*(4*D*a^2+-1*b*x*(C*a+3*A*b))+1//4*a^-1*b^-1*(a+b*x^2)^-2*(x*(A*b+-1*C*a)+-1*a*(B+-1*D*a*b^-1))+1//8*a^-5//2*b^-3//2*(C*a+3*A*b)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-1*(a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :(A*a^-3*log(x)+-1//2*A*a^-3*log(a+b*x^2)+1//4*a^-1*b^-1*(a+b*x^2)^-2*(A*b+x*(B*b+-1*D*a)+-1*C*a)+1//8*a^-2*b^-1*(a+b*x^2)^-1*(x*(D*a+3*B*b)+4*A*b)+1//8*a^-5//2*b^-3//2*(D*a+3*B*b)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-2*(a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :(B*a^-3*log(x)+-1*A*a^-3*x^-1+-1//2*B*a^-3*log(a+b*x^2)+1//8*a^-2*(a+b*x^2)^-1*(4B+-1*x*(-3C+7*A*b*a^-1))+-1//8*a^-7//2*b^-1//2*(-3*C*a+15*A*b)*arctan(x*a^-1//2*b^(1/2))+1//4*a^-1*b^-1*(a+b*x^2)^-2*(B*b+-1*D*a+-1*b*x*(-1C+A*b*a^-1)))
@test integrate(x^-3*(a+b*x^2)^-3*(A+B*x+C*x^2+D*x^3), x) == :((1/2)*a^-4*(-1*C*a+3*A*b)*log(a+b*x^2)+-1*B*a^-3*x^-1+-1*a^-4*(-1*C*a+3*A*b)*log(x)+-1//2*A*a^-3*x^-2+-1//4*a^-1*(a+b*x^2)^-2*(-1C+x*(-1D+B*b*a^-1)+A*b*a^-1)+-1//8*a^-3*(a+b*x^2)^-1*(x*(-3*D*a+7*B*b)+-4*C*a+8*A*b)+-1//8*a^-7//2*b^-1//2*(-3*D*a+15*B*b)*arctan(x*a^-1//2*b^(1/2)))
@test integrate((5+x^2)^-2*(-1x+4*x^3), x) == :(2*log(5+x^2)+21*(10+2*x^2)^-1)
@test integrate((-2+x^2)^-1//2*(x^3+-1x), x) == :((-2+x^2)^(1/2)+1//3*(-2+x^2)^3//2)
@test integrate((1+2*x^2)^-1*(-1*x^2+2*x^4), x) == :(-1x+1//3*x^3+(1/2)*2^(1/2)*arctan(x*2^(1/2)))
@test integrate((1+x^2)^-1*(x^3+x^4), x) == :((1/2)*x^2+-1x+-1//2*log(1+x^2)+1//3*x^3+arctan(x))
@test integrate(x^6*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(1//5*b^-4*x^5*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//7*b^-3*x^7*(d*b^2+f*a^2+-1*a*b*e)+1//9*b^-2*x^9*(b*e+-1*a*f)+1//11*f*b^-1*x^11+x*a^2*b^-6*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1*a^5//2*b^-13//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2))+-1//3*a*b^-5*x^3*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2))
@test integrate(x^4*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(1//3*b^-4*x^3*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//5*b^-3*x^5*(d*b^2+f*a^2+-1*a*b*e)+1//7*b^-2*x^7*(b*e+-1*a*f)+1//9*f*b^-1*x^9+a^3//2*b^-11//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2))+-1*a*x*b^-5*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2))
@test integrate(x^2*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(x*b^-4*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//3*b^-3*x^3*(d*b^2+f*a^2+-1*a*b*e)+1//5*b^-2*x^5*(b*e+-1*a*f)+1//7*f*b^-1*x^7+-1*a^(1/2)*b^-9//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate((a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(x*b^-3*(d*b^2+f*a^2+-1*a*b*e)+1//3*b^-2*x^3*(b*e+-1*a*f)+1//5*f*b^-1*x^5+a^-1//2*b^-7//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-2*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(x*b^-2*(b*e+-1*a*f)+-1*c*a^-1*x^-1+1//3*f*b^-1*x^3+-1*a^-3//2*b^-5//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-4*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(f*x*b^-1+a^-2*x^-1*(b*c+-1*a*d)+-1//3*c*a^-1*x^-3+a^-5//2*b^-3//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-6*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(-1*a^-3*x^-1*(c*b^2+e*a^2+-1*a*b*d)+-1//5*c*a^-1*x^-5+1//3*a^-2*x^-3*(b*c+-1*a*d)+-1*a^-7//2*b^-1//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-8*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(a^-4*x^-1*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1//3*a^-3*x^-3*(c*b^2+e*a^2+-1*a*b*d)+-1//7*c*a^-1*x^-7+1//5*a^-2*x^-5*(b*c+-1*a*d)+a^-9//2*b^(1/2)*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-10*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//5*a^-3*x^-5*(c*b^2+e*a^2+-1*a*b*d)+-1//9*c*a^-1*x^-9+1//3*a^-4*x^-3*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//7*a^-2*x^-7*(b*c+-1*a*d)+-1*b*a^-5*x^-1*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1*a^-11//2*b^3//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-12*(a+b*x^2)^-1*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//7*a^-3*x^-7*(c*b^2+e*a^2+-1*a*b*d)+-1//11*c*a^-1*x^-11+1//5*a^-4*x^-5*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//9*a^-2*x^-9*(b*c+-1*a*d)+a^-6*b^2*x^-1*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+a^-13//2*b^5//2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)*arctan(x*a^-1//2*b^(1/2))+-1//3*b*a^-5*x^-3*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2))
@test integrate(x^6*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//6*b^-5*x^3*(-11*f*a^3+5*c*b^3+-7*a*d*b^2+9*b*e*a^2)+1//7*b^-3*x^7*(b*e+-2*a*f)+1//9*f*b^-2*x^9+(1/2)*a^-1*x^7*(a+b*x^2)^-1*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+(1/2)*a^3//2*b^-13//2*(-11*f*a^3+5*c*b^3+-7*a*d*b^2+9*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+-1//2*a*x*b^-6*(-11*f*a^3+5*c*b^3+-7*a*d*b^2+9*b*e*a^2)+-1//10*a^-1*b^-4*x^5*(-11*f*a^3+5*c*b^3+-7*a*d*b^2+9*b*e*a^2))
@test integrate(x^4*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :((1/2)*x*b^-5*(-9*f*a^3+3*c*b^3+-5*a*d*b^2+7*b*e*a^2)+1//5*b^-3*x^5*(b*e+-2*a*f)+1//7*f*b^-2*x^7+(1/2)*a^-1*x^5*(a+b*x^2)^-1*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+-1//2*a^(1/2)*b^-11//2*(-9*f*a^3+3*c*b^3+-5*a*d*b^2+7*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+-1//6*a^-1*b^-4*x^3*(-9*f*a^3+3*c*b^3+-5*a*d*b^2+7*b*e*a^2))
@test integrate(x^2*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//3*b^-3*x^3*(b*e+-2*a*f)+1//5*f*b^-2*x^5+(1/2)*a^-1*x^3*(a+b*x^2)^-1*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+(1/2)*a^-1//2*b^-9//2*(c*b^3+-7*f*a^3+-3*a*d*b^2+5*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+-1//2*x*a^-1*b^-4*(c*b^3+-7*f*a^3+-3*a*d*b^2+5*b*e*a^2))
@test integrate((a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(x*b^-3*(b*e+-2*a*f)+1//3*f*b^-2*x^3+(1/2)*x*a^-1*(a+b*x^2)^-1*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+(1/2)*a^-3//2*b^-7//2*(c*b^3+5*f*a^3+a*d*b^2+-3*b*e*a^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-2*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(f*x*b^-2+-1*c*a^-2*x^-1+-1//2*x*a^-1*(a+b*x^2)^-1*(-1d+a*e*b^-1+b*c*a^-1+-1*f*a^2*b^-2)+-1//2*a^-5//2*b^-5//2*(3*c*b^3+3*f*a^3+-1*a*d*b^2+-1*b*e*a^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-4*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(a^-3*x^-1*(-1*a*d+2*b*c)+-1//3*c*a^-2*x^-3+(1/2)*x*a^-1*(a+b*x^2)^-1*(e+c*a^-2*b^2+-1*a*f*b^-1+-1*b*d*a^-1)+(1/2)*a^-7//2*b^-3//2*(f*a^3+5*c*b^3+b*e*a^2+-3*a*d*b^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-6*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(-1*a^-4*x^-1*(e*a^2+3*c*b^2+-2*a*b*d)+-1//5*c*a^-2*x^-5+1//3*a^-3*x^-3*(-1*a*d+2*b*c)+-1//2*x*a^-4*(a+b*x^2)^-1*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1//2*a^-9//2*b^-1//2*(-1*f*a^3+7*c*b^3+-5*a*d*b^2+3*b*e*a^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-8*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(a^-5*x^-1*(-1*f*a^3+4*c*b^3+-3*a*d*b^2+2*b*e*a^2)+-1//3*a^-4*x^-3*(e*a^2+3*c*b^2+-2*a*b*d)+-1//7*c*a^-2*x^-7+1//5*a^-3*x^-5*(-1*a*d+2*b*c)+(1/2)*a^-11//2*b^(1/2)*(-3*f*a^3+9*c*b^3+-7*a*d*b^2+5*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+(1/2)*b*x*a^-5*(a+b*x^2)^-1*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2))
@test integrate(x^-10*(a+b*x^2)^-2*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//5*a^-4*x^-5*(e*a^2+3*c*b^2+-2*a*b*d)+-1//9*c*a^-2*x^-9+1//3*a^-5*x^-3*(-1*f*a^3+4*c*b^3+-3*a*d*b^2+2*b*e*a^2)+1//7*a^-3*x^-7*(-1*a*d+2*b*c)+-1*b*a^-6*x^-1*(-2*f*a^3+5*c*b^3+-4*a*d*b^2+3*b*e*a^2)+-1//2*a^-13//2*b^3//2*(-5*f*a^3+11*c*b^3+-9*a*d*b^2+7*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+-1//2*x*a^-6*b^2*(a+b*x^2)^-1*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2))
@test integrate(x^8*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(1//6*b^-6*x^3*(-23*f*a^3+5*c*b^3+-9*a*d*b^2+15*b*e*a^2)+1//7*b^-4*x^7*(b*e+-3*a*f)+1//9*f*b^-3*x^9+-1//4*a*x*b^-7*(-63*f*a^3+15*c*b^3+-27*a*d*b^2+43*b*e*a^2)+-1//20*a^-1*b^-5*x^5*(-29*f*a^3+5*c*b^3+-9*a*d*b^2+17*b*e*a^2)+1//4*a^-1*x^9*(a+b*x^2)^-2*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+1//8*a^3//2*b^-15//2*(-143*f*a^3+35*c*b^3+-63*a*d*b^2+99*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+-1//8*x*a^2*b^-7*(a+b*x^2)^-1*(-17*f*a^3+5*c*b^3+-9*a*d*b^2+13*b*e*a^2))
@test integrate(x^6*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :((1/2)*x*b^-6*(-21*f*a^3+3*c*b^3+-7*a*d*b^2+13*b*e*a^2)+1//5*b^-4*x^5*(b*e+-3*a*f)+1//7*f*b^-3*x^7+-1//8*a^(1/2)*b^-13//2*(-99*f*a^3+15*c*b^3+-35*a*d*b^2+63*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+-1//12*a^-1*b^-5*x^3*(-27*f*a^3+3*c*b^3+-7*a*d*b^2+15*b*e*a^2)+1//4*a^-1*x^7*(a+b*x^2)^-2*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+1//8*a*x*b^-6*(a+b*x^2)^-1*(-15*f*a^3+3*c*b^3+-7*a*d*b^2+11*b*e*a^2))
@test integrate(x^4*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(1//3*b^-4*x^3*(b*e+-3*a*f)+1//5*f*b^-3*x^5+-1//4*x*a^-1*b^-5*(c*b^3+-25*f*a^3+-5*a*d*b^2+13*b*e*a^2)+-1//8*x*b^-5*(a+b*x^2)^-1*(c*b^3+-13*f*a^3+-5*a*d*b^2+9*b*e*a^2)+1//4*a^-1*x^5*(a+b*x^2)^-2*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+1//8*a^-1//2*b^-11//2*(-63*f*a^3+3*c*b^3+-15*a*d*b^2+35*b*e*a^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^2*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(x*b^-4*(b*e+-3*a*f)+1//3*f*b^-3*x^3+1//4*a^-1*x^3*(a+b*x^2)^-2*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+1//8*a^-3//2*b^-9//2*(c*b^3+35*f*a^3+-15*b*e*a^2+3*a*d*b^2)*arctan(x*a^-1//2*b^(1/2))+-1//8*x*a^-1*b^-4*(a+b*x^2)^-1*(c*b^3+11*f*a^3+-7*b*e*a^2+3*a*d*b^2))
@test integrate((a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(f*x*b^-3+1//4*x*a^-1*(a+b*x^2)^-2*(c+-1*a*b^-3*(d*b^2+f*a^2+-1*a*b*e))+1//8*a^-5//2*b^-7//2*(-15*f*a^3+3*c*b^3+a*d*b^2+3*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+1//8*x*a^-2*b^-3*(a+b*x^2)^-1*(3*c*b^3+9*f*a^3+a*d*b^2+-5*b*e*a^2))
@test integrate(x^-2*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(-1*c*a^-3*x^-1+-1//4*x*a^-1*(a+b*x^2)^-2*(-1d+a*e*b^-1+b*c*a^-1+-1*f*a^2*b^-2)+-1//8*a^-7//2*b^-5//2*(-3*f*a^3+15*c*b^3+-1*b*e*a^2+-3*a*d*b^2)*arctan(x*a^-1//2*b^(1/2))+-1//8*x*a^-3*b^-2*(a+b*x^2)^-1*(5*f*a^3+7*c*b^3+-1*b*e*a^2+-3*a*d*b^2))
@test integrate(x^-4*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(a^-4*x^-1*(-1*a*d+3*b*c)+-1//3*c*a^-3*x^-3+1//4*x*a^-1*(a+b*x^2)^-2*(e+c*a^-2*b^2+-1*a*f*b^-1+-1*b*d*a^-1)+1//8*a^-9//2*b^-3//2*(f*a^3+35*c*b^3+-15*a*d*b^2+3*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+1//8*x*a^-4*b^-1*(a+b*x^2)^-1*(f*a^3+11*c*b^3+-7*a*d*b^2+3*b*e*a^2))
@test integrate(x^-6*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(-1*a^-5*x^-1*(e*a^2+6*c*b^2+-3*a*b*d)+-1//5*c*a^-3*x^-5+1//3*a^-4*x^-3*(-1*a*d+3*b*c)+-1//4*x*a^-4*(a+b*x^2)^-2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1//8*x*a^-5*(a+b*x^2)^-1*(-3*f*a^3+15*c*b^3+-11*a*d*b^2+7*b*e*a^2)+-1//8*a^-11//2*b^-1//2*(-3*f*a^3+63*c*b^3+-35*a*d*b^2+15*b*e*a^2)*arctan(x*a^-1//2*b^(1/2)))
@test integrate(x^-8*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(a^-6*x^-1*(-1*f*a^3+10*c*b^3+-6*a*d*b^2+3*b*e*a^2)+-1//3*a^-5*x^-3*(e*a^2+6*c*b^2+-3*a*b*d)+-1//7*c*a^-3*x^-7+1//5*a^-4*x^-5*(-1*a*d+3*b*c)+1//8*a^-13//2*b^(1/2)*(-15*f*a^3+99*c*b^3+-63*a*d*b^2+35*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+1//4*b*x*a^-5*(a+b*x^2)^-2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//8*b*x*a^-6*(a+b*x^2)^-1*(-7*f*a^3+19*c*b^3+-15*a*d*b^2+11*b*e*a^2))
@test integrate(x^-10*(a+b*x^2)^-3*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//5*a^-5*x^-5*(e*a^2+6*c*b^2+-3*a*b*d)+-1//9*c*a^-3*x^-9+1//3*a^-6*x^-3*(-1*f*a^3+10*c*b^3+-6*a*d*b^2+3*b*e*a^2)+1//7*a^-4*x^-7*(-1*a*d+3*b*c)+-1*b*a^-7*x^-1*(-3*f*a^3+15*c*b^3+-10*a*d*b^2+6*b*e*a^2)+-1//8*a^-15//2*b^3//2*(-35*f*a^3+143*c*b^3+-99*a*d*b^2+63*b*e*a^2)*arctan(x*a^-1//2*b^(1/2))+-1//4*x*a^-6*b^2*(a+b*x^2)^-2*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1//8*x*a^-7*b^2*(a+b*x^2)^-1*(-11*f*a^3+23*c*b^3+-19*a*d*b^2+15*b*e*a^2))
@test integrate(x^5*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//5*b^-6*(a+b*x^2)^5//2*(c*b^3+-10*f*a^3+-3*a*d*b^2+6*b*e*a^2)+1//7*b^-6*(a+b*x^2)^7//2*(d*b^2+10*f*a^2+-4*a*b*e)+1//9*b^-6*(a+b*x^2)^9//2*(b*e+-5*a*f)+1//11*f*b^-6*(a+b*x^2)^11//2+a^2*b^-6*(a+b*x^2)^(1/2)*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1//3*a*b^-6*(a+b*x^2)^3//2*(-5*f*a^3+2*c*b^3+-3*a*d*b^2+4*b*e*a^2))
@test integrate(x^3*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//3*b^-5*(a+b*x^2)^3//2*(c*b^3+-4*f*a^3+-2*a*d*b^2+3*b*e*a^2)+1//5*b^-5*(a+b*x^2)^5//2*(d*b^2+6*f*a^2+-3*a*b*e)+1//7*b^-5*(a+b*x^2)^7//2*(b*e+-4*a*f)+1//9*f*b^-5*(a+b*x^2)^9//2+-1*a*b^-5*(a+b*x^2)^(1/2)*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2))
@test integrate(x*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(b^-4*(a+b*x^2)^(1/2)*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//3*b^-4*(a+b*x^2)^3//2*(d*b^2+3*f*a^2+-2*a*b*e)+1//5*b^-4*(a+b*x^2)^5//2*(b*e+-3*a*f)+1//7*f*b^-4*(a+b*x^2)^7//2)
@test integrate(x^-1*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(b^-3*(a+b*x^2)^(1/2)*(d*b^2+f*a^2+-1*a*b*e)+-1*c*a^-1//2*arctanh(a^-1//2*(a+b*x^2)^(1/2))+1//3*b^-3*(a+b*x^2)^3//2*(b*e+-2*a*f)+1//5*f*b^-3*(a+b*x^2)^5//2)
@test integrate(x^-3*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(b^-2*(a+b*x^2)^(1/2)*(b*e+-1*a*f)+(1/2)*a^-3//2*(b*c+-2*a*d)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+1//3*f*b^-2*(a+b*x^2)^3//2+-1//2*c*a^-1*x^-2*(a+b*x^2)^(1/2))
@test integrate(x^-5*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(f*b^-1*(a+b*x^2)^(1/2)+-1//8*a^-5//2*(3*c*b^2+8*e*a^2+-4*a*b*d)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1//4*c*a^-1*x^-4*(a+b*x^2)^(1/2)+1//8*a^-2*x^-2*(a+b*x^2)^(1/2)*(-4*a*d+3*b*c))
@test integrate(x^-7*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//16*a^-7//2*(-16*f*a^3+5*c*b^3+-6*a*d*b^2+8*b*e*a^2)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1//6*c*a^-1*x^-6*(a+b*x^2)^(1/2)+-1//16*a^-3*x^-2*(a+b*x^2)^(1/2)*(5*c*b^2+8*e*a^2+-6*a*b*d)+1//24*a^-2*x^-4*(a+b*x^2)^(1/2)*(-6*a*d+5*b*c))
@test integrate(x^-9*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//8*c*a^-1*x^-8*(a+b*x^2)^(1/2)+-1//128*b*a^-9//2*(-64*f*a^3+35*c*b^3+-40*a*d*b^2+48*b*e*a^2)*arctanh(a^-1//2*(a+b*x^2)^(1/2))+-1//192*a^-3*x^-4*(a+b*x^2)^(1/2)*(35*c*b^2+48*e*a^2+-40*a*b*d)+1//48*a^-2*x^-6*(a+b*x^2)^(1/2)*(-8*a*d+7*b*c)+1//128*a^-4*x^-2*(a+b*x^2)^(1/2)*(-64*f*a^3+35*c*b^3+-40*a*d*b^2+48*b*e*a^2))
@test integrate(x^4*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//10*f*b^-1*x^9*(a+b*x^2)^(1/2)+1//80*b^-2*x^7*(a+b*x^2)^(1/2)*(-9*a*f+10*b*e)+1//256*a^2*b^-11//2*(-63*f*a^3+96*c*b^3+-80*a*d*b^2+70*b*e*a^2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//384*b^-4*x^3*(a+b*x^2)^(1/2)*(-63*f*a^3+96*c*b^3+-80*a*d*b^2+70*b*e*a^2)+1//480*b^-3*x^5*(a+b*x^2)^(1/2)*(63*f*a^2+80*d*b^2+-70*a*b*e)+-1//256*a*x*b^-5*(a+b*x^2)^(1/2)*(-63*f*a^3+96*c*b^3+-80*a*d*b^2+70*b*e*a^2))
@test integrate(x^2*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//128*a*b^-9//2*(-35*f*a^3+64*c*b^3+-48*a*d*b^2+40*b*e*a^2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//8*f*b^-1*x^7*(a+b*x^2)^(1/2)+1//48*b^-2*x^5*(a+b*x^2)^(1/2)*(-7*a*f+8*b*e)+1//128*x*b^-4*(a+b*x^2)^(1/2)*(-35*f*a^3+64*c*b^3+-48*a*d*b^2+40*b*e*a^2)+1//192*b^-3*x^3*(a+b*x^2)^(1/2)*(35*f*a^2+48*d*b^2+-40*a*b*e))
@test integrate((a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//16*b^-7//2*(-5*f*a^3+16*c*b^3+-8*a*d*b^2+6*b*e*a^2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//6*f*b^-1*x^5*(a+b*x^2)^(1/2)+1//16*x*b^-3*(a+b*x^2)^(1/2)*(5*f*a^2+8*d*b^2+-6*a*b*e)+1//24*b^-2*x^3*(a+b*x^2)^(1/2)*(-5*a*f+6*b*e))
@test integrate(x^-2*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(1//8*b^-5//2*(3*f*a^2+8*d*b^2+-4*a*b*e)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1*c*a^-1*x^-1*(a+b*x^2)^(1/2)+1//4*f*b^-1*x^3*(a+b*x^2)^(1/2)+1//8*x*b^-2*(a+b*x^2)^(1/2)*(-3*a*f+4*b*e))
@test integrate(x^-4*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :((1/2)*b^-3//2*(-1*a*f+2*b*e)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+(1/2)*f*x*b^-1*(a+b*x^2)^(1/2)+-1//3*c*a^-1*x^-3*(a+b*x^2)^(1/2)+1//3*a^-2*x^-1*(a+b*x^2)^(1/2)*(-3*a*d+2*b*c))
@test integrate(x^-6*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(f*b^-1//2*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+-1//5*c*a^-1*x^-5*(a+b*x^2)^(1/2)+-1//15*a^-3*x^-1*(a+b*x^2)^(1/2)*(8*c*b^2+15*e*a^2+-10*a*b*d)+1//15*a^-2*x^-3*(a+b*x^2)^(1/2)*(-5*a*d+4*b*c))
@test integrate(x^-8*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//7*c*a^-1*x^-7*(a+b*x^2)^(1/2)+-1//105*a^-3*x^-3*(a+b*x^2)^(1/2)*(24*c*b^2+35*e*a^2+-28*a*b*d)+1//35*a^-2*x^-5*(a+b*x^2)^(1/2)*(-7*a*d+6*b*c)+1//105*a^-4*x^-1*(a+b*x^2)^(1/2)*(-105*f*a^3+48*c*b^3+-56*a*d*b^2+70*b*e*a^2))
@test integrate(x^-10*(a+b*x^2)^-1//2*(c+d*x^2+e*x^4+f*x^6), x) == :(-1//9*c*a^-1*x^-9*(a+b*x^2)^(1/2)+-1//105*a^-3*x^-5*(a+b*x^2)^(1/2)*(16*c*b^2+21*e*a^2+-18*a*b*d)+1//63*a^-2*x^-7*(a+b*x^2)^(1/2)*(-9*a*d+8*b*c)+1//315*a^-4*x^-3*(a+b*x^2)^(1/2)*(-105*f*a^3+64*c*b^3+-72*a*d*b^2+84*b*e*a^2)+-2//315*b*a^-5*x^-1*(a+b*x^2)^(1/2)*(-105*f*a^3+64*c*b^3+-72*a*d*b^2+84*b*e*a^2))
@test integrate(x^8*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(1//16*b^-15//2*(-429*D*a^3+16*A*b^3+-72*B*a*b^2+198*C*b*a^2)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//6*D*b^-3*x^9*(a+b*x^2)^-3//2+1//7*a^-1*x^9*(a+b*x^2)^-7//2*(A+-1*a*b^-3*(B*b^2+D*a^2+-1*C*a*b))+-1//16*x*a^-1*b^-7*(a+b*x^2)^(1/2)*(-3*a*(24*B*b^2+143*D*a^2+-66*C*a*b)+16*A*b^3)+-1//30*a^-2*b^-5*x^5*(a+b*x^2)^-1//2*(-3*a*(24*B*b^2+143*D*a^2+-66*C*a*b)+16*A*b^3)+-1//35*a^-2*b^-3*x^9*(a+b*x^2)^-5//2*(-1*a*(9*B*b^2+23*D*a^2+-16*C*a*b)+2*A*b^3)+-1//210*a^-2*b^-4*x^7*(a+b*x^2)^-3//2*(-3*a*(24*B*b^2+143*D*a^2+-66*C*a*b)+16*A*b^3)+1//24*a^-2*b^-6*x^3*(a+b*x^2)^(1/2)*(-3*a*(24*B*b^2+143*D*a^2+-66*C*a*b)+16*A*b^3))
@test integrate(x^6*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(1//8*b^-13//2*(8*B*b^2+99*D*a^2+-36*C*a*b)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+1//4*D*b^-3*x^7*(a+b*x^2)^-3//2+1//7*a^-1*x^7*(a+b*x^2)^-7//2*(A+-1*a*b^-3*(B*b^2+D*a^2+-1*C*a*b))+-1//8*x*a^-1*b^-6*(a+b*x^2)^(1/2)*(8*B*b^2+99*D*a^2+-36*C*a*b)+1//5*a^-1*b^-3*x^7*(a+b*x^2)^-5//2*(B*b^2+3*D*a^2+-2*C*a*b)+1//12*a^-1*b^-5*x^3*(a+b*x^2)^-1//2*(8*B*b^2+99*D*a^2+-36*C*a*b)+1//60*a^-1*b^-4*x^5*(a+b*x^2)^-3//2*(8*B*b^2+99*D*a^2+-36*C*a*b))
@test integrate(x^4*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :((1/2)*b^-11//2*(-9*D*a+2*C*b)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+(1/2)*D*x*b^-5*(a+b*x^2)^(1/2)+-1//3*x*b^-5*(a+b*x^2)^-1//2*(-15*D*a+4*C*b)+1//7*a^-1*x^5*(a+b*x^2)^-7//2*(A+-1*a*b^-3*(B*b^2+D*a^2+-1*C*a*b))+1//3*a*x*b^-5*(a+b*x^2)^-3//2*(C*b+-3*D*a)+1//35*a^-2*b^-3*x^5*(a+b*x^2)^-5//2*(a*(5*B*b^2+19*D*a^2+-12*C*a*b)+2*A*b^3))
@test integrate((a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(A*x*a^-1*(a+b*x^2)^-7//2+1//3*a^-2*x^3*(a+b*x^2)^-7//2*(B*a+6*A*b)+1//15*a^-3*x^5*(a+b*x^2)^-7//2*(a*(3*C*a+4*B*b)+24*A*b^2)+1//105*a^-4*x^7*(a+b*x^2)^-7//2*(a*(8*B*b^2+15*D*a^2+6*C*a*b)+48*A*b^3))
@test integrate(x^-2*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(-1*A*a^-1*x^-1*(a+b*x^2)^-7//2+-1*x*a^-2*(a+b*x^2)^-7//2*(-1*B*a+8*A*b)+-1//3*a^-3*x^3*(a+b*x^2)^-7//2*(-1*a*(C*a+6*B*b)+48*A*b^2)+-1//15*a^-4*x^5*(a+b*x^2)^-7//2*(-3*D*a^3+4*b*(-1*a*(C*a+6*B*b)+48*A*b^2))+-2//105*b*a^-5*x^7*(a+b*x^2)^-7//2*(-3*D*a^3+4*b*(-1*a*(C*a+6*B*b)+48*A*b^2)))
@test integrate(x^-4*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(-1//3*A*a^-1*x^-3*(a+b*x^2)^-7//2+1//3*x*a^-3*(a+b*x^2)^-7//2*(-3*a*(-1*C*a+8*B*b)+80*A*b^2)+1//3*a^-4*x^3*(a+b*x^2)^-7//2*(-1*a*(-1*D*a^2+48*B*b^2+-6*C*a*b)+160*A*b^3)+1//3*a^-2*x^-1*(a+b*x^2)^-7//2*(-3*B*a+10*A*b)+4//15*b*a^-5*x^5*(a+b*x^2)^-7//2*(-1*a*(-1*D*a^2+48*B*b^2+-6*C*a*b)+160*A*b^3)+8//105*a^-6*b^2*x^7*(a+b*x^2)^-7//2*(-1*a*(-1*D*a^2+48*B*b^2+-6*C*a*b)+160*A*b^3))
@test integrate(x^-6*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(-1//3*a^-3*x^-1*(a+b*x^2)^-7//2*(-1*a*(-3*C*a+10*B*b)+24*A*b^2)+-1//5*A*a^-1*x^-5*(a+b*x^2)^-7//2+-1//21*x*a^-4*(a+b*x^2)^-7//2*(-1*a*(3*D*a^2+80*B*b^2+-24*C*a*b)+192*A*b^3)+-1//35*x*a^-5*(a+b*x^2)^-5//2*(-2*a*(3*D*a^2+80*B*b^2+-24*C*a*b)+384*A*b^3)+-1//105*x*a^-7*(a+b*x^2)^-1//2*(-16*a*(3*D*a^2+80*B*b^2+-24*C*a*b)+3072*A*b^3)+-1//105*x*a^-6*(a+b*x^2)^-3//2*(-8*a*(3*D*a^2+80*B*b^2+-24*C*a*b)+1536*A*b^3)+1//15*a^-2*x^-3*(a+b*x^2)^-7//2*(-5*B*a+12*A*b))
@test integrate(x^-8*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(-1//7*A*a^-1*x^-7*(a+b*x^2)^-7//2+-1//15*a^-3*x^-3*(a+b*x^2)^-7//2*(-1*a*(-5*C*a+12*B*b)+24*A*b^2)+1//3*a^-4*x^-1*(a+b*x^2)^-7//2*(-1*a*(3*D*a^2+24*B*b^2+-10*C*a*b)+48*A*b^3)+1//5*a^-2*x^-5*(a+b*x^2)^-7//2*(-1*B*a+2*A*b)+8//21*b*x*a^-5*(a+b*x^2)^-7//2*(-1*a*(3*D*a^2+24*B*b^2+-10*C*a*b)+48*A*b^3)+16//35*b*x*a^-6*(a+b*x^2)^-5//2*(-1*a*(3*D*a^2+24*B*b^2+-10*C*a*b)+48*A*b^3)+64//105*b*x*a^-7*(a+b*x^2)^-3//2*(-1*a*(3*D*a^2+24*B*b^2+-10*C*a*b)+48*A*b^3)+128//105*b*x*a^-8*(a+b*x^2)^-1//2*(-1*a*(3*D*a^2+24*B*b^2+-10*C*a*b)+48*A*b^3))
@test integrate(x^-10*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6), x) == :(-1//9*A*a^-1*x^-9*(a+b*x^2)^-7//2+-1//45*a^-3*x^-5*(a+b*x^2)^-7//2*(-9*a*(-1*C*a+2*B*b)+32*A*b^2)+1//45*a^-4*x^-3*(a+b*x^2)^-7//2*(-3*a*(5*D*a^2+24*B*b^2+-12*C*a*b)+128*A*b^3)+1//63*a^-2*x^-7*(a+b*x^2)^-7//2*(-9*B*a+16*A*b)+-256//315*x*a^-9*b^2*(a+b*x^2)^-1//2*(-3*a*(5*D*a^2+24*B*b^2+-12*C*a*b)+128*A*b^3)+-128//315*x*a^-8*b^2*(a+b*x^2)^-3//2*(-3*a*(5*D*a^2+24*B*b^2+-12*C*a*b)+128*A*b^3)+-32//105*x*a^-7*b^2*(a+b*x^2)^-5//2*(-3*a*(5*D*a^2+24*B*b^2+-12*C*a*b)+128*A*b^3)+-16//63*x*a^-6*b^2*(a+b*x^2)^-7//2*(-3*a*(5*D*a^2+24*B*b^2+-12*C*a*b)+128*A*b^3)+-2//9*b*a^-5*x^-1*(a+b*x^2)^-7//2*(-3*a*(5*D*a^2+24*B*b^2+-12*C*a*b)+128*A*b^3))
@test integrate((a+b*x^2)^-1//2*(c*x^5+d*x^7+e*x^9+f*x^11), x) == :(1//5*b^-6*(a+b*x^2)^5//2*(c*b^3+-10*f*a^3+-3*a*d*b^2+6*b*e*a^2)+1//7*b^-6*(a+b*x^2)^7//2*(d*b^2+10*f*a^2+-4*a*b*e)+1//9*b^-6*(a+b*x^2)^9//2*(b*e+-5*a*f)+1//11*f*b^-6*(a+b*x^2)^11//2+a^2*b^-6*(a+b*x^2)^(1/2)*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+-1//3*a*b^-6*(a+b*x^2)^3//2*(-5*f*a^3+2*c*b^3+-3*a*d*b^2+4*b*e*a^2))
@test integrate((a+b*x^2)^-1//2*(c*x^3+d*x^5+e*x^7+f*x^9), x) == :(1//3*b^-5*(a+b*x^2)^3//2*(c*b^3+-4*f*a^3+-2*a*d*b^2+3*b*e*a^2)+1//5*b^-5*(a+b*x^2)^5//2*(d*b^2+6*f*a^2+-3*a*b*e)+1//7*b^-5*(a+b*x^2)^7//2*(b*e+-4*a*f)+1//9*f*b^-5*(a+b*x^2)^9//2+-1*a*b^-5*(a+b*x^2)^(1/2)*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2))
@test integrate((a+b*x^2)^-1//2*(c*x+d*x^3+e*x^5+f*x^7), x) == :(b^-4*(a+b*x^2)^(1/2)*(c*b^3+-1*f*a^3+b*e*a^2+-1*a*d*b^2)+1//3*b^-4*(a+b*x^2)^3//2*(d*b^2+3*f*a^2+-2*a*b*e)+1//5*b^-4*(a+b*x^2)^5//2*(b*e+-3*a*f)+1//7*f*b^-4*(a+b*x^2)^7//2)
@test integrate(x^2*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6+F*x^8), x) == :((1/2)*b^-11//2*(-9*F*a+2*D*b)*arctanh(x*b^(1/2)*(a+b*x^2)^-1//2)+(1/2)*F*x*b^-5*(a+b*x^2)^(1/2)+-1*x*b^-5*(a+b*x^2)^-1//2*(D*b+-4*F*a)+1//7*a^-1*b^-4*x^3*(a+b*x^2)^-7//2*(A*b^4+-1*a*(B*b^3+-1*F*a^3+D*b*a^2+-1*C*a*b^2))+1//35*a^-2*b^-4*x^3*(a+b*x^2)^-5//2*(a*(-24*F*a^3+3*B*b^3+-10*C*a*b^2+17*D*b*a^2)+4*A*b^4)+1//105*a^-3*b^-4*x^3*(a+b*x^2)^-3//2*(a*(6*B*b^3+162*F*a^3+-71*D*b*a^2+15*C*a*b^2)+8*A*b^4))
@test integrate(x^-2*(a+b*x^2)^-9//2*(A+B*x^2+C*x^4+D*x^6+F*x^8), x) == :(-1*A*a^-1*x^-1*(a+b*x^2)^-7//2+-1*x*a^-2*(a+b*x^2)^-7//2*(-1*B*a+8*A*b)+-1//3*a^-3*x^3*(a+b*x^2)^-7//2*(-1*a*(C*a+6*B*b)+48*A*b^2)+-1//15*a^-4*x^5*(a+b*x^2)^-7//2*(-1*a*(3*D*a^2+24*B*b^2+4*C*a*b)+192*A*b^3)+-1//105*a^-5*x^7*(a+b*x^2)^-7//2*(-1*a*(15*F*a^3+48*B*b^3+6*D*b*a^2+8*C*a*b^2)+384*A*b^4))

max_stars_repo_path: test/1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.8 P(x) (c x)^m (a+b x^2)^p.jl


text:
using MovingWeightedLeastSquares
using Base.Random
using Base.Test
using Base.Iterators
using DynamicPolynomials, MultivariatePolynomials

include("sample-generator.jl")
include("cell-linked-list.jl")
include("1d-approximation.jl")
include("single-array-input.jl")
include("1d-derivative.jl")
include("2d-approximation.jl")
include("2d-derivative.jl")


text:
function ispangram(input::AbstractString)

end


text:
#
# Functionality that should be eventually included in Base.LibGit2
#

const NULL_OID = LibGit2.Oid()

"""
    checkout_safe!(repo::GitRepo, rev::AbstractString) -> nothing

Checks out a revision from a Git repository. The `rev` can be a branch, tag, SHA, or partial
SHA. Meant to be a temporary workaround to deal with limitations of `LibGit2.checkout!`
"""
function checkout_safe!(repo::LibGit2.GitRepo, rev::AbstractString)
    # Avoid using `checkout!` unless the repo has no uncommitted changes
    length(LibGit2.GitStatus(repo)) == 0 || error("Repo has uncommited changes")

    # TODO: `revparseid` issues:
    # - Doesn't notice ambiguous revisions. i.e. a tag with the same name as a branch.
    #   Issue could be from the `libgit2` library.
    # - Returns a value when it cannot find the specified revision
    oid = LibGit2.revparseid(repo, rev)
    oid != NULL_OID || error("Unable to determine SHA of \"$rev\"")

    # TODO: `checkout!` seems broken in multiple ways including:
    # - Only uses full SHAs. Can not use tags, branches, or partial SHAs
    # - Defaults to using a `force=true` which wipes out any uncommited changes
    # - Using `force=false` seems to always fail
    LibGit2.checkout!(repo, string(oid), force=true)
end

function checkout_safe!(repo::LibGit2.GitRepo, oid::Base.LibGit2.Oid)
    length(LibGit2.GitStatus(repo)) == 0 || error("Repo has uncommited changes")
    LibGit2.checkout!(repo, string(oid), force=true)
end


text:
"""
    ReachSolution{SN, RSN<:ReachSet{SN}} <: AbstractSolution

Type that wraps the solution of a reachability problem as a sequence of lazy
sets, and a dictionary of options.

### Fields

- `Xk`       -- the list of [`ReachSet`](@ref)s
- `options`  -- the dictionary of options
"""
struct ReachSolution{SN, RSN<:ReachSet{SN}} <: AbstractSolution
    Xk::Vector{RSN}
    options::AbstractOptions
end

# constructor with no options
ReachSolution(Xk::Vector{RSN}) where {SN, RSN<:ReachSet{SN}} =
    ReachSolution{SN, RSN<:ReachSet{SN}}(Xk, Options())


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "PSDD and PlainVtree Learner Test" begin
    data = dataset(twenty_datasets("nltcs"); do_shuffle=false, batch_size=-1);
    train_data = train(data);
    clt = learn_chow_liu_tree(WXData(train_data));
    #clt = parse_clt("circuits/test.clt")
    vtree = learn_vtree_from_clt(clt; vtree_mode="balanced"); # or "linear"
    @test vtree isa PlainVtree

    mktempdir() do tmp
        save(vtree, "$tmp/test.vtree.dot")
        psdd = compile_psdd_from_clt(clt, vtree);
        @test psdd isa ProbΔ
        save_as_dot(psdd, "$tmp/test.psdd.dot")
    end

end