text:
using LinearAlgebra, DiffEqOperators, Random, Test, BandedMatrices, SparseArrays

@show "2D"

@testset "2D Multiplication with no boundary points and dx = 1.0" begin

    # Test (Lxx + Lyy)*M, dx = 1.0, no coefficient
    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N)

    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(0.1i)+sin(0.1j)
        end
    end

    Lxx = CenteredDifference{1}(2,2,1.0,N)
    Lyy = CenteredDifference{2}(2,2,1.0,N)
    A = Lxx + Lyy

    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lxx*M)[1:N,2:N+1] +(Lyy*M)[2:N+1,1:N])

    # Test a single axis, multiple operators: (Lx + Lxx)*M, dx = 1.0
    Lx = CenteredDifference{1}(1,2,1.0,N)
    A = Lx + Lxx

    M_temp = zeros(100,102)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx*M)+(Lxx*M))

    # Test a single axis, multiple operators: (Ly + Lyy)*M, dx = 1.0, no coefficient
    Ly = CenteredDifference{2}(1,2,1.0,N)
    A = Ly + Lyy

    M_temp = zeros(102,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly*M)+(Lyy*M))

    # Test multiple operators on both axis: (Lx + Ly + Lxx + Lyy)*M, no coefficient
    A = Lx + Ly + Lxx + Lyy
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx*M)[1:N,2:N+1] +(Ly*M)[2:N+1,1:N] + (Lxx*M)[1:N,2:N+1] +(Lyy*M)[2:N+1,1:N])
end

@testset "2D Multiplication with identical bpc and dx = 1.0" begin

    # Test (Lxxxx + Lyyyy)*M, dx = 1.0, no coefficient, two boundary points on each axis
    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N)

    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(0.1i)+sin(0.1j)
        end
    end

    Lx4 = CenteredDifference{1}(4,4,1.0,N)
    Ly4 = CenteredDifference{2}(4,4,1.0,N)
    A = Lx4 + Ly4

    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

    # Test a single axis, multiple operators: (Lxxx + Lxxxx)*M, dx = 1.0
    Lx3 = CenteredDifference{1}(3,4,1.0,N)
    A = Lx3 + Lx4

    M_temp = zeros(100,102)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)+(Lx4*M))

    # Test a single axis, multiple operators: (Lyyy + Lyyyy)*M, dx = 1.0, no coefficient
    Ly3 = CenteredDifference{2}(3,4,1.0,N)
    A = Ly3 + Ly4

    M_temp = zeros(102,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)+(Ly4*M))

    # Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N] + (Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

    # Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Lx3 + Ly3
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N])

end

@testset "2D Multiplication with differing bpc and dx = 1.0" begin

    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N+2)
    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(0.1i)+sin(0.1j)
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,1.0,N)
    # Lx3 has 1 boundary point
    Lx3 = CenteredDifference{1}(3,3,1.0,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,1.0,N)

    # Test a single axis, multiple operators: (Lxx+Lxxxx)*M, dx = 1.0
    A = Lx2+Lx4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx4*M))

    # Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M, dx = 1.0
    A += Lx3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))

    # Ly2 has 0 boundary points
    Ly2 = CenteredDifference{2}(2,2,1.0,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,1.0,N)
    # Ly4 has 2 boundary points
    Ly4 = CenteredDifference{2}(4,4,1.0,N)
    M_temp = zeros(N+2,N)

    # Test a single axis, multiple operators: (Lyy+Lyyyy)*M, dx = 1.0
    A = Ly2+Ly4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly4*M))

    # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M, dx = 1.0
    A += Ly3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))

    # Test multiple operators on both axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N] + (Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

end

@testset "2D Multiplication with identical bpc and non-trivial dx = dy = 0.1" begin

    # Test (Lxxxx + Lyyyy)*M, dx = 0.1, dy = 0.01, no coefficient, two boundary points on each axis
    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N)

    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(0.1i)+sin(0.1j)
        end
    end

    Lx4 = CenteredDifference{1}(4,4,0.1,N)
    Ly4 = CenteredDifference{2}(4,4,0.1,N)
    A = Lx4 + Ly4

    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

    # Test a single axis, multiple operators: (Lxxx + Lxxxx)*M, dx = 0.1
    Lx3 = CenteredDifference{1}(3,4,0.1,N)
    A = Lx3 + Lx4

    M_temp = zeros(100,102)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)+(Lx4*M))

    # Test a single axis, multiple operators: (Lyyy + Lyyyy)*M, dx = 0.01, no coefficient
    Ly3 = CenteredDifference{2}(3,4,0.1,N)
    A = Ly3 + Ly4

    M_temp = zeros(102,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)+(Ly4*M))

    # Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient dx =
    A = Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N] + (Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

    # Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Lx3 + Ly3
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N])

end

@testset "2D Multiplication with identical bpc and non-trivial dx = 0.1, dy = 0.25" begin

    # Test (Lxxxx + Lyyyy)*M, dx = 0.1, dy = 0.01, no coefficient, two boundary points on each axis
    dx = 0.1
    dy = 0.25
    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N)

    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(dx*i)+sin(dy*j)
        end
    end

    Lx4 = CenteredDifference{1}(4,4,dx,N)
    Ly4 = CenteredDifference{2}(4,4,dy,N)
    A = Lx4 + Ly4

    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

    # Test a single axis, multiple operators: (Lxxx + Lxxxx)*M, dx = 0.1
    Lx3 = CenteredDifference{1}(3,4,dx,N)
    A = Lx3 + Lx4

    M_temp = zeros(100,102)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)+(Lx4*M))

    # Test a single axis, multiple operators: (Lyyy + Lyyyy)*M, dx = 0.01, no coefficient
    Ly3 = CenteredDifference{2}(3,4,dy,N)
    A = Ly3 + Ly4

    M_temp = zeros(102,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)+(Ly4*M))

    # Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient dx =
    A = Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N] + (Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

    # Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Lx3 + Ly3
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N])

end

@testset "2D Multiplication with differing bpc and non-trivial dx = dy = 0.1" begin

    dx = 0.1
    dy = 0.1
    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N+2)
    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(dx*i)+sin(dy*j)
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Test a single axis, multiple operators: (Lxx+Lxxxx)*M, dx = 1.0
    A = Lx2+Lx4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx4*M))

    # Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M, dx = 1.0
    A += Lx3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))

    # Ly2 has 0 boundary points
    Ly2 = CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = CenteredDifference{2}(4,4,dy,N)
    M_temp = zeros(N+2,N)

    # Test a single axis, multiple operators: (Lyy+Lyyyy)*M, dx = 1.0
    A = Ly2+Ly4
    mul!(M_temp, A, M)
    @test isapprox(M_temp, ((Ly2*M) + (Ly4*M)), atol = 1e-7)

    # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M, dx = 1.0
    A += Ly3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))

    # Test multiple operators on both axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N] + (Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

end

@testset "2D Multiplication with differing bpc and non-trivial dx = 0.1, dy = 0.25" begin

    dx = 0.1
    dy = 0.25
    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N+2)
    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(dx*i)+sin(dy*j)
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Test a single axis, multiple operators: (Lxx+Lxxxx)*M, dx = 1.0
    A = Lx2+Lx4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx4*M)) || maximum(abs.(M_temp - ((Lx2*M) + (Lx4*M)))) < 1e-10

    # Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M, dx = 1.0
    A += Lx3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))

    # Ly2 has 0 boundary points
    Ly2 = CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = CenteredDifference{2}(4,4,dy,N)
    M_temp = zeros(N+2,N)

    # Test a single axis, multiple operators: (Lyy+Lyyyy)*M, dx = 1.0
    A = Ly2+Ly4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly4*M))

    # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M, dx = 1.0
    A += Ly3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))

    # Test multiple operators on both axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N] + (Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])

end

# THis testset uses the last testset which has a several non-trivial cases,
# and additionally tests coefficient handling. All operators are handled by the
# fast 2D/3D dispatch.
@testset "2D coefficient handling" begin

    dx = 0.1
    dy = 0.25
    N = 100
    M = zeros(N+2,N+2)
    M_temp = zeros(N,N+2)
    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(dx*i)+sin(dy*j)
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = 5.5*CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = 0.5*CenteredDifference{1}(4,4,dx,N)

    # Test a single axis, multiple operators: (Lxx+Lxxxx)*M, dx = 1.0
    A = Lx2+Lx4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx4*M))

    # Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M, dx = 1.0
    A += Lx3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))

    # Ly2 has 0 boundary points
    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = 2.0*CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)
    M_temp = zeros(N+2,N)

    # Test a single axis, multiple operators: (Lyy+Lyyyy)*M, dx = 1.0
    A = Ly2+Ly4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly4*M))

    # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M, dx = 1.0
    A += Ly3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))

    # Test multiple operators on both axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(100,100)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Lx3*M)[1:N,2:N+1] +(Ly3*M)[2:N+1,1:N] + (Lx4*M)[1:N,2:N+1] +(Ly4*M)[2:N+1,1:N])
end

@show "irregular"

@testset "x and y are both irregular grids" begin

    N = 100
    dx = cumsum(rand(N+2))
    dy = cumsum(rand(N+2))
    M = zeros(N+2,N+2)

    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(dx[i])+sin(dy[j])
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Ly2 has 0 boundary points
    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)

    # Test composition of all first-dimension operators
    A = Lx2+Lx3+Lx4
    M_temp = zeros(N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Lx2*M + Lx3*M + Lx4*M)

    # Test composition of all second-dimension operators
    A = Ly2+Ly3+Ly4
    M_temp = zeros(N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Ly2*M + Ly3*M + Ly4*M)

    # Test composition of all operators
    A = Lx2+Lx3+Lx4+Ly2+Ly3+Ly4
    M_temp = zeros(N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Lx3*M)[1:N,2:N+1]+(Lx4*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Ly3*M)[2:N+1,1:N]+(Ly4*M)[2:N+1,1:N])

end

@testset "regular x grid (dx=0.25) and irregular y grid" begin

    N = 100
    dx = 0.25
    dx_vec = 0.25:0.25:(N+2)*0.25
    dy = cumsum(rand(N+2))
    M = zeros(N+2,N+2)

    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(dx*i)+sin(dy[j])
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Ly2 has 0 boundary points
    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)

    # Test that composition of all x-operators works
    A = Lx2 + Lx3 + Lx4
    M_temp = zeros(N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Lx2*M + Lx3*M + Lx4*M)

    # Test that composition of all y-operators works
    A = Ly2 + Ly3 + Ly4
    M_temp = zeros(N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Ly2*M + Ly3*M + Ly4*M)

    # Test that composition of both x and y operators works
    A = Lx2 + Ly2 + Lx3 + Ly3 + Ly4 + Lx4
    M_temp = zeros(N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Lx3*M)[1:N,2:N+1]+(Lx4*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Ly3*M)[2:N+1,1:N]+(Ly4*M)[2:N+1,1:N])

    # Last case where we now have some `irregular-grid` operators operating on the
    # regular-spaced axis x

    # These operators are operating on the regular grid x, but are constructed as though
    # they were irregular grid operators. Hence we test if we can seperate irregular and
    # regular gird operators on the same axis
    # Lx2 has 0 boundary points
    _Lx2 = 4.532*CenteredDifference{1}(2,2,dx_vec,N)
    # Lx3 has 1 boundary point
    _Lx3 = 0.235*CenteredDifference{1}(3,3,dx_vec,N)
    # Lx4 has 2 boundary points
    _Lx4 = CenteredDifference{1}(4,4,dx_vec,N)

    A += _Lx2 + _Lx3 + _Lx4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((_Lx2*M)[1:N,2:N+1]+(_Lx3*M)[1:N,2:N+1]+(_Lx4*M)[1:N,2:N+1]+(Lx2*M)[1:N,2:N+1]+(Lx3*M)[1:N,2:N+1]+(Lx4*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Ly3*M)[2:N+1,1:N]+(Ly4*M)[2:N+1,1:N])

end

@testset "irregular x grid and regular y grid (dy = 0.25)" begin

    N = 100
    dy = 0.25
    dy_vec = 0.25:0.25:(N+2)*0.25
    dx = cumsum(rand(N+2))
    M = zeros(N+2,N+2)

    for i in 1:N+2
        for j in 1:N+2
            M[i,j] = cos(dx[i])+sin(dy*j)
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Ly2 has 0 boundary points
    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)

    # Test that composition of all x-operators works
    A = Lx2 + Lx3 + Lx4
    M_temp = zeros(N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Lx2*M + Lx3*M + Lx4*M)

    # Test that composition of all y-operators works
    A = Ly2 + Ly3 + Ly4
    M_temp = zeros(N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Ly2*M + Ly3*M + Ly4*M)

    # Test that composition of both x and y operators works
    A = Lx2 + Ly2 + Lx3 + Ly3 + Ly4 + Lx4
    M_temp = zeros(N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Lx3*M)[1:N,2:N+1]+(Lx4*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Ly3*M)[2:N+1,1:N]+(Ly4*M)[2:N+1,1:N])

    # Last case where we now have some `irregular-grid` operators operating on the
    # regular-spaced axis y

    # These operators are operating on the regular grid y, but are constructed as though
    # they were irregular grid operators. Hence we test if we can seperate irregular and
    # regular gird operators on the same axis
    # Ly2 has 0 boundary points
    _Ly2 = CenteredDifference{2}(2,2,dy_vec,N)
    # Ly3 has 1 boundary point
    _Ly3 = CenteredDifference{2}(3,3,dy_vec,N)
    # Ly4 has 2 boundary points
    _Ly4 = 12.1*CenteredDifference{2}(4,4,dy_vec,N)

    A += _Ly2 + _Ly3 + _Ly4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1]+(Lx3*M)[1:N,2:N+1]+(Lx4*M)[1:N,2:N+1]+(Ly2*M)[2:N+1,1:N]+(Ly3*M)[2:N+1,1:N]+(Ly4*M)[2:N+1,1:N]+(_Ly2*M)[2:N+1,1:N]+(_Ly3*M)[2:N+1,1:N]+(_Ly4*M)[2:N+1,1:N])

end

################################################################################
# 3D Multiplication Tests
################################################################################

@show "3D"

@testset "3D Multiplication with no boundary points and dx = dy = dz = 1.0" begin

    # Test (Lxx + Lyy + Lzz)*M, dx = dy = dz = 1.0, no coefficient
    N = 100
    M = zeros(N+2,N+2,N+2)
    M_temp = zeros(N,N,N)

    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(0.1i)+sin(0.1j) + exp(0.01k)
            end
        end
    end

    Lxx = CenteredDifference{1}(2,2,1.0,N)
    Lyy = CenteredDifference{2}(2,2,1.0,N)
    Lzz = CenteredDifference{3}(2,2,1.0,N)
    A = Lxx + Lyy + Lzz

    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lxx*M)[1:N,2:N+1,2:N+1] + (Lyy*M)[2:N+1,1:N,2:N+1] + (Lzz*M)[2:N+1,2:N+1,1:N])

    # Test a single axis, multiple operators: (Lx + Lxx)*M, dx = 1.0
    Lx = CenteredDifference{1}(1,2,1.0,N)
    A = Lx + Lxx

    M_temp = zeros(N,N+2,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx*M)+(Lxx*M))

    # Test a single axis, multiple operators: (Ly + Lyy)*M, dy = 1.0, no coefficient
    Ly = CenteredDifference{2}(1,2,1.0,N)
    A = Ly + Lyy

    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly*M)+(Lyy*M))

    # Test a single axis, multiple operators: (Lz + Lzz)*M, dz = 1.0, no coefficient
    Lz = CenteredDifference{3}(1,2,1.0,N)
    A = Lz + Lzz

    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lz*M)+(Lzz*M))

    # Test multiple operators on both axis: (Lx + Ly + Lxx + Lyy)*M, no coefficient
    A = Lx + Ly + Lxx + Lyy
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx*M)[1:N,2:N+1,:] +(Ly*M)[2:N+1,1:N,:] + (Lxx*M)[1:N,2:N+1,:] +(Lyy*M)[2:N+1,1:N,:])

    # Test multiple operators on both axis: (Lx + Lxx + Lz + Lzz)*M, no coefficient
    A = Lx + Lxx + Lz + Lzz
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx*M)[1:N,:,2:N+1] + (Lxx*M)[1:N,:,2:N+1]  + (Lz*M)[2:N+1,:,1:N] +(Lzz*M)[2:N+1,:,1:N])

    # Test multiple operators on both axis: (Ly + Lyy + Lz + Lzz)*M, no coefficient
    A = Ly + Lyy + Lz + Lzz
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly*M)[:,1:N,2:N+1] + (Lyy*M)[:,1:N,2:N+1] + (Lz*M)[:,2:N+1,1:N] +(Lzz*M)[:,2:N+1,1:N])

    # Test multiple operators on both axis: (Lx + Ly + Lxx + Lyy + Lz + Lzz)*M, no coefficient
    A = Lx + Ly + Lxx + Lyy + Lz + Lzz
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx*M)[1:N,2:N+1,2:N+1] +(Ly*M)[2:N+1,1:N,2:N+1] + (Lxx*M)[1:N,2:N+1,2:N+1] +(Lyy*M)[2:N+1,1:N,2:N+1] + (Lz*M)[2:N+1,2:N+1,1:N] +(Lzz*M)[2:N+1,2:N+1,1:N])

end

@testset "3D Multiplication with identical bpc and dx = dy = dz = 1.0" begin

    N = 100
    M = zeros(N+2,N+2,N+2)
    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(0.1i)+sin(0.1j) + exp(0.01k)
            end
        end
    end

    # Test a single axis, multiple operators: (Lxxx + Lxxxx)*M, dx = dy = dz = 1.0
    # Lx3 and Lx4 have the same number of boundary points
    Lx3 = CenteredDifference{1}(3,4,1.0,N)
    Lx4 = CenteredDifference{1}(4,4,1.0,N)
    A = Lx3 + Lx4

    M_temp = zeros(N,N+2,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)+(Lx4*M))

    # Test a single axis, multiple operators: (Lyyy + Lyyyy)*M, dx = dy = dz = 1.0, no coefficient
    # Ly3 and Ly4 have the same number of boundary points
    Ly3 = CenteredDifference{2}(3,4,1.0,N)
    Ly4 = CenteredDifference{2}(4,4,1.0,N)
    A = Ly3 + Ly4

    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)+(Ly4*M))

    # Test a single axis, multiple operators: (Lzzz + Lzzzz)*M, dx = dy = dz = 1.0, no coefficient
    # Lz3 and Lz4 have the same number of boundary points
    Lz3 = CenteredDifference{3}(3,4,1.0,N)
    Lz4 = CenteredDifference{3}(4,4,1.0,N)
    A = Lz3 + Lz4

    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lz3*M)+(Lz4*M))

    # Test (Lxxxx + Lyyyy)*M, dx = 1.0, no coefficient, two boundary points on each axis
    M_temp = zeros(N,N,N+2)
    A = Lx4 + Ly4
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])

    # Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Lx3 + Ly3
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:])

    # Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])

    # Test (Lxxxx + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
    M_temp = zeros(N,N+2,N)
    A = Lx4 + Lz4
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test (Lxxx + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Lx3 + Lz3
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N])

    # Test multiple operators on both axis: (Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
    A = Lx3 + Lz3 + Lx4 + Lz4
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test (Lyyyy + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
    M_temp = zeros(N+2,N,N)
    A = Ly4 + Lz4
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test (Lyyy + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Ly3 + Lz3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N])

    # Test multiple operators on both axis: (Lyyy + Lzzz + Lyyyy + Lzzzz)*M, no coefficient
    A = Ly3 + Lz3 + Ly4 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test a single operator on each axis: (Lx3 + Ly3 + Lz3)*M, no coefficient
    A = Lx3 + Ly3 + Lz3
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N])

    # Test a single operator on each axis: (Lx4 + Ly4 + Lz4)*M, no coefficient
    A = Lx4 + Ly4 + Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])

    # Test multiple operators on each axis: (Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4)*M, no coefficient
    A = Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])

end

@show "3D BPC"

@testset "3D Multiplication with differing bpc and dx = dy = dz = 1.0" begin

    N = 100
    M = zeros(N+2,N+2,N+2)
    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(0.1i)+sin(0.1j) + exp(0.01k)
            end
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,1.0,N)
    # Lx3 has 1 boundary point
    Lx3 = CenteredDifference{1}(3,3,1.0,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,1.0,N)

    # Test a single axis, multiple operators: (Lxx+Lxxxx)*M, dx = 1.0
    A = Lx2+Lx4
    M_temp = zeros(N,N+2,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx4*M))

    # Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M, dx = 1.0
    A += Lx3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))

    # Ly2 has 0 boundary points
    Ly2 = CenteredDifference{2}(2,2,1.0,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,1.0,N)
    # Ly4 has 2 boundary points
    Ly4 = CenteredDifference{2}(4,4,1.0,N)

    # Test a single axis, multiple operators: (Lyy+Lyyyy)*M, dx = 1.0
    A = Ly2+Ly4
    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly4*M))

    # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M, dx = 1.0
    A += Ly3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))

    # Lz2 has 0 boundary points
    Lz2 = CenteredDifference{3}(2,2,1.0,N)
    # Lz3 has 1 boundary point
    Lz3 = CenteredDifference{3}(3,3,1.0,N)
    # Lz4 has 2 boundary points
    Lz4 = CenteredDifference{3}(4,4,1.0,N)

    # Test a single axis, multiple operators: (Lzy+Lzzzz)*M, dz = 1.0
    A = Lz2+Lz4
    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lz2*M) + (Lz4*M))

    # Test a single axis, multiple operators: (Lzz++Lzzz+Lzzzz)*M, dz = 1.0
    A += Lz3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lz2*M) + (Lz3*M) + (Lz4*M))

    # Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
    A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
    A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test operators on all three axis (Lxx + Lyy + Lzz + Lxxx + Lyyy + Lzzz + Lxxxx + Lyyyy + Lzzzz)*M, no coefficient
    A = Lx2 + Ly2 + Lz2 + Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]

end

@testset "3D Multiplication with identical bpc and dx = 0.5, dy = 0.25, dz = 0.05, no coefficients" begin

    N = 100
    M = zeros(N+2,N+2,N+2)
    dx = 0.5
    dy = 0.25
    dz = 0.05
    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(dx*i)+sin(dy*j) + exp(dz*k)
            end
        end
    end

    # Test a single axis, multiple operators: (Lxxx + Lxxxx)*M
    # Lx3 and Lx4 have the same number of boundary points
    Lx3 = CenteredDifference{1}(3,4,dx,N)
    Lx4 = CenteredDifference{1}(4,4,dx,N)
    A = Lx3 + Lx4

    M_temp = zeros(N,N+2,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)+(Lx4*M))

    # Test a single axis, multiple operators: (Lyyy + Lyyyy)*M
    # Ly3 and Ly4 have the same number of boundary points
    Ly3 = CenteredDifference{2}(3,4,dy,N)
    Ly4 = CenteredDifference{2}(4,4,dy,N)
    A = Ly3 + Ly4

    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)+(Ly4*M))

    # Test a single axis, multiple operators: (Lzzz + Lzzzz)*M
    # Lz3 and Lz4 have the same number of boundary points
    Lz3 = CenteredDifference{3}(3,4,dz,N)
    Lz4 = CenteredDifference{3}(4,4,dz,N)
    A = Lz3 + Lz4

    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lz3*M)+(Lz4*M))

    # Test (Lxxxx + Lyyyy)*M, two boundary points on each axis
    M_temp = zeros(N,N,N+2)
    A = Lx4 + Ly4
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])

    # Test (Lxxx + Lyyy)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Lx3 + Ly3
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:])

    # Test multiple operators on both axis: (Lxxx + Lyyy + Lxxxx + Lyyyy)*M
    A = Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])

    # Test (Lxxxx + Lzzzz)*M, two boundary points on each axis
    M_temp = zeros(N,N+2,N)
    A = Lx4 + Lz4
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test (Lxxx + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Lx3 + Lz3
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N])

    # Test multiple operators on both axis: (Lxxx + Lzzz + Lxxxx + Lzzzz)*M
    A = Lx3 + Lz3 + Lx4 + Lz4
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test (Lyyyy + Lzzzz)*M, dx = 1.0, no coefficient, two boundary points on each axis
    M_temp = zeros(N+2,N,N)
    A = Ly4 + Lz4
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test (Lyyy + Lzzz)*M, no coefficient. These operators have non-symmetric interior stencils
    A = Ly3 + Lz3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N])

    # Test multiple operators on both axis: (Lyyy + Lzzz + Lyyyy + Lzzzz)*M
    A = Ly3 + Lz3 + Ly4 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test a single operator on each axis: (Lx3 + Ly3 + Lz3)*M
    A = Lx3 + Ly3 + Lz3
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N])

    # Test a single operator on each axis: (Lx4 + Ly4 + Lz4)*M
    A = Lx4 + Ly4 + Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])

    # Test multiple operators on each axis: (Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4)*M
    A = Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx3*M)[1:N,2:N+1,2:N+1] + (Ly3*M)[2:N+1,1:N,2:N+1] +(Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] + (Ly4*M)[2:N+1,1:N,2:N+1] +(Lz4*M)[2:N+1,2:N+1,1:N])

end

@testset "3D Multiplication with differing bpc and dx = 0.5, dy = 0.25, dz = 0.05, no coefficients" begin

    N = 100
    M = zeros(N+2,N+2,N+2)
    dx = 0.5
    dy = 0.25
    dz = 0.05
    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(dx*i)+sin(dy*j) + exp(dz*k)
            end
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Test a single axis, multiple operators: (Lxx+Lxxxx)*M
    A = Lx2+Lx4
    M_temp = zeros(N,N+2,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx4*M))

    # Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M
    A += Lx3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))

    # Ly2 has 0 boundary points
    Ly2 = CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = CenteredDifference{2}(4,4,dy,N)

    # Test a single axis, multiple operators: (Lyy+Lyyyy)*M
    A = Ly2+Ly4
    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly4*M))

    # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M
    A += Ly3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))

    # Lz2 has 0 boundary points
    Lz2 = CenteredDifference{3}(2,2,dz,N)
    # Lz3 has 1 boundary point
    Lz3 = CenteredDifference{3}(3,3,dz,N)
    # Lz4 has 2 boundary points
    Lz4 = CenteredDifference{3}(4,4,dz,N)

    # Test a single axis, multiple operators: (Lzy+Lzzzz)*M
    A = Lz2+Lz4
    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lz2*M) + (Lz4*M))

    # Test a single axis, multiple operators: (Lzz++Lzzz+Lzzzz)*M
    A += Lz3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lz2*M) + (Lz3*M) + (Lz4*M))

    # Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
    A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
    A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test operators on all three axis (Lxx + Lyy + Lzz + Lxxx + Lyyy + Lzzz + Lxxxx + Lyyyy + Lzzzz)*M, no coefficient
    A = Lx2 + Ly2 + Lz2 + Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]

end

@show "3D with coefficients"

# This is the same test set as the above test set with the addition of coefficients
@testset "3D Multiplication with coefficients" begin

    N = 100
    M = zeros(N+2,N+2,N+2)
    dx = 0.5
    dy = 0.25
    dz = 0.05
    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(dx*i)+sin(dy*j) + exp(dz*k)
            end
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = 1.234*CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 2.456*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Test a single axis, multiple operators: (Lxx+Lxxxx)*M
    A = Lx2+Lx4
    M_temp = zeros(N,N+2,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx4*M))

    # Test a single axis, multiple operators: (Lxx++Lxxx+Lxxxx)*M
    A += Lx3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M) + (Lx3*M) + (Lx4*M))

    # Ly2 has 0 boundary points
    Ly2 = CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = 4.014*CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 1.49*CenteredDifference{2}(4,4,dy,N)

    # Test a single axis, multiple operators: (Lyy+Lyyyy)*M
    A = Ly2+Ly4
    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly4*M))

    # Test a single axis, multiple operators: (Lyy++Lyyy+Lyyyy)*M
    A += Ly3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M) + (Ly3*M) + (Ly4*M))

    # Lz2 has 0 boundary points
    Lz2 = 1.546*CenteredDifference{3}(2,2,dz,N)
    # Lz3 has 1 boundary point
    Lz3 = CenteredDifference{3}(3,3,dz,N)
    # Lz4 has 2 boundary points
    Lz4 = 0.55*CenteredDifference{3}(4,4,dz,N)

    # Test a single axis, multiple operators: (Lzy+Lzzzz)*M
    A = Lz2+Lz4
    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lz2*M) + (Lz4*M))

    # Test a single axis, multiple operators: (Lzz++Lzzz+Lzzzz)*M
    A += Lz3
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lz2*M) + (Lz3*M) + (Lz4*M))

    # Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:])

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
    A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M
    A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test operators on all three axis (Lxx + Lyy + Lzz + Lxxx + Lyyy + Lzzz + Lxxxx + Lyyyy + Lzzzz)*M, no coefficient
    A = Lx2 + Ly2 + Lz2 + Lx3 + Ly3 + Lz3 + Lx4 + Ly4 + Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]

end

@show "3D irregular"

@testset "x, y, and z are all irregular grids" begin

    N = 100
    dx = cumsum(rand(N+2))
    dy = cumsum(rand(N+2))
    dz = cumsum(rand(N+2))
    M = zeros(N+2,N+2,N+2)

    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(dx[i])+sin(dy[j]) + exp(dz[k])
            end
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Ly2 has 0 boundary points
    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)

    # Lz2 has 0 boundary points
    Lz2 = CenteredDifference{3}(2,2,dz,N)
    # Lz3 has 1 boundary point
    Lz3 = CenteredDifference{3}(3,3,dz,N)
    # Lz4 has 2 boundary points
    Lz4 = CenteredDifference{3}(4,4,dz,N)

    # Test composition of all first-dimension operators
    A = Lx2+Lx3+Lx4
    M_temp = zeros(N,N+2,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Lx2*M + Lx3*M + Lx4*M)

    # Test composition of all second-dimension operators
    A = Ly2+Ly3+Ly4
    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Ly2*M + Ly3*M + Ly4*M)

    # Test composition of all third-dimension operators
    A = Lz2+Lz3+Lz4
    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ (Lz2*M + Lz3*M + Lz4*M)

    # Test multiple operators on two axis: (Lxx + Lyy + Lxxx + Lyyy + Lxxxx + Lyyyy)*M, no coefficient
    A = Lx2 + Ly2 + Lx3 + Ly3 + Lx4 + Ly4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)

    # Numerical errors accumulating for this test case, more so than the other tests
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Lx3*M)[1:N,2:N+1,:] +(Ly3*M)[2:N+1,1:N,:] + (Lx4*M)[1:N,2:N+1,:] +(Ly4*M)[2:N+1,1:N,:]) rtol=1e-3

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
    A = Lx2 + Lz2 + Lx3 + Lz3 + Lx4 + Lz4
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lx3*M)[1:N,:,2:N+1] +(Lz3*M)[2:N+1,:,1:N] + (Lx4*M)[1:N,:,2:N+1] +(Lz4*M)[2:N+1,:,1:N])

    # Test multiple operators on two axis: (Lxx + Lzz + Lxxx + Lzzz + Lxxxx + Lzzzz)*M, no coefficient
    A = Ly2 + Lz2 + Ly3 + Lz3 + Ly4 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)

    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Ly3*M)[:,1:N,2:N+1] +(Lz3*M)[:,2:N+1,1:N] + (Ly4*M)[:,1:N,2:N+1] +(Lz4*M)[:,2:N+1,1:N])

    # Test composition of all operators
    A = Lx2+Lx3+Lx4+Ly2+Ly3+Ly4+Lz2+Lz3+Lz4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1])+(Lz4*M)[2:N+1,2:N+1,1:N]

end

@testset "irregular x grid and regular y and z grids (dy = 0.25, dz = 0.05)" begin

    N = 100
    dy = 0.25
    dz = 0.05
    dy_vec = 0.25:0.25:(N+2)*0.25
    dz_vec = 0.05:0.05:(N+2)*0.05
    dx = cumsum(rand(N+2))
    M = zeros(N+2,N+2,N+2)

    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(dx[i])+sin(dy*j)+exp(dz*k)
            end
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Ly2 has 0 boundary points
    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)

    # Lz2 has 0 boundary points
    Lz2 = 1.4*CenteredDifference{3}(2,2,dz,N)
    # Lz3 has 1 boundary point
    Lz3 = CenteredDifference{3}(3,3,dz,N)
    # Lz4 has 2 boundary points
    Lz4 = CenteredDifference{3}(4,4,dz,N)

    # Test that composition of both x and y operators works
    A = Lx2 + Ly2 + Lx3 + Ly3 + Ly4 + Lx4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Lx3*M)[1:N,2:N+1,:]+(Lx4*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Ly3*M)[2:N+1,1:N,:]+(Ly4*M)[2:N+1,1:N,:])

    # Test that composition of both x and z operators works
    A = Lx2 + Lz2 + Lx3 + Lz3 + Lz4 + Lx4
    M_temp = zeros(N,N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,:,2:N+1]+(Lx3*M)[1:N,:,2:N+1]+(Lx4*M)[1:N,:,2:N+1]+(Lz2*M)[2:N+1,:,1:N]+(Lz3*M)[2:N+1,:,1:N]+(Lz4*M)[2:N+1,:,1:N])

    # Test that composition of x, y, and z operators works
    A = Lx2 + Lz2 + Lx3 + Lz3 + Lz4 + Lx4 + Ly2 + Ly3 + Ly4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N])

    # Last case where we now have some `irregular-grid` operators operating on the
    # regular-spaced axis y

    # These operators are operating on the regular grid y, but are constructed as though
    # they were irregular grid operators. Hence we test if we can seperate irregular and
    # regular gird operators on the same axis
    # Ly2 has 0 boundary points
    _Ly2 = CenteredDifference{2}(2,2,dy_vec,N)
    # Lz3 has 1 boundary point
    _Lz3 = CenteredDifference{3}(3,3,dz_vec,N)
    # Ly4 has 2 boundary points
    _Ly4 = 12.1*CenteredDifference{2}(4,4,dy_vec,N)

    A += _Ly2 + _Lz3 + _Ly4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N] + (_Ly2*M)[2:N+1,1:N,2:N+1]+(_Lz3*M)[2:N+1,2:N+1,1:N]+(_Ly4*M)[2:N+1,1:N,2:N+1])

end

@testset "irregular y grid and regular x and z gris (dx = 0.25, dz = 0.05)" begin

    N = 100
    dx = 0.25
    dz = 0.05
    dx_vec = 0.25:0.25:(N+2)*0.25
    dz_vec = 0.05:0.05:(N+2)*0.05
    dy = cumsum(rand(N+2))
    M = zeros(N+2,N+2,N+2)

    for i in 1:N+2
        for j in 1:N+2
            for k in 1:N+2
                M[i,j,k] = cos(dx*i)+sin(dy[j])+exp(dz*k)
            end
        end
    end

    # Lx2 has 0 boundary points
    Lx2 = CenteredDifference{1}(2,2,dx,N)
    # Lx3 has 1 boundary point
    Lx3 = 1.45*CenteredDifference{1}(3,3,dx,N)
    # Lx4 has 2 boundary points
    Lx4 = CenteredDifference{1}(4,4,dx,N)

    # Ly2 has 0 boundary points
    Ly2 = 8.14*CenteredDifference{2}(2,2,dy,N)
    # Ly3 has 1 boundary point
    Ly3 = CenteredDifference{2}(3,3,dy,N)
    # Ly4 has 2 boundary points
    Ly4 = 4.567*CenteredDifference{2}(4,4,dy,N)

    # Lz2 has 0 boundary points
    Lz2 = 1.4*CenteredDifference{3}(2,2,dz,N)
    # Lz3 has 1 boundary point
    Lz3 = CenteredDifference{3}(3,3,dz,N)
    # Lz4 has 2 boundary points
    Lz4 = CenteredDifference{3}(4,4,dz,N)

    # Test that composition of both x and y operators works
    A = Lx2 + Ly2 + Lx3 + Ly3 + Ly4 + Lx4
    M_temp = zeros(N,N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,:]+(Lx3*M)[1:N,2:N+1,:]+(Lx4*M)[1:N,2:N+1,:]+(Ly2*M)[2:N+1,1:N,:]+(Ly3*M)[2:N+1,1:N,:]+(Ly4*M)[2:N+1,1:N,:])

    # Test that composition of both y and z operators works
    A = Ly2 + Ly3 + Ly4 + Lz2 + Lz3 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    # Need to figure out why this test is exploding
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])

    ############################################################################
    # Tests to isolate the above problem
    ############################################################################

    A = Ly2 + Ly3 + Ly4
    M_temp = zeros(N+2,N,N+2)
    mul!(M_temp, A, M)
    @test M_temp ≈ Ly2*M + Ly3*M + Ly4*M

    ### Test the addition of Lz2, Lz3, Lz4 seperately

    A = Ly2 + Ly3 + Ly4 + Lz2
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N])

    A = Ly2 + Ly3 + Ly4 + Lz3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz3*M)[:,2:N+1,1:N])

    A = Ly2 + Ly3 + Ly4 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz4*M)[:,2:N+1,1:N])

    ###

    A = Ly2 + Ly3 + Ly4 + Lz2 + Lz3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N])

    A = Ly2 + Ly3 + Ly4 + Lz2 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])

    A = Ly2 + Ly3 + Ly4 + Lz3 + Lz4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])

    ###

    A = Ly2 + Ly3 + Ly4 + Lz2 + Lz2
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz2*M)[:,2:N+1,1:N])
    # It appears that multiple z operators with some y operators is causing the issues

    ###

    A = Lz2 + Lz3 + Lz4
    M_temp = zeros(N+2,N+2,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ Lz2*M + Lz3*M + Lz4*M

    A = Lz2 + Lz3 + Lz4 + Ly2
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])

    A = Lz2 + Lz3 + Lz4 + Ly3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
    # It seems that the y paddign could be the issue

    A = Lz2 + Lz3 + Lz4 + Ly4
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
    # It seems that the y paddign could be the issue

    A = Lz2 + Lz3 + Lz4 + Ly4 + Ly3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
    # It seems that the y paddign could be the issue

    A = Lz2 + Lz3 + Lz4 + Ly4 + Ly2
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Ly4*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N]+(Lz4*M)[:,2:N+1,1:N])
    # It seems that the y paddign could be the issue

    ###

    A = Lz2 + Lz3 +Ly3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N])

    A = Lz3 +Ly3
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly3*M)[:,1:N,2:N+1]+(Lz3*M)[:,2:N+1,1:N])

    A = Lz2 + Lz3 +Ly2
    M_temp = zeros(N+2,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Ly2*M)[:,1:N,2:N+1]+(Lz2*M)[:,2:N+1,1:N]+(Lz3*M)[:,2:N+1,1:N])

    # It seems that the padding of y is forcing multiple z operators to fail

    ############################################################################
    ############################################################################

    # Test that composition of x, y, and z operators works
    A = Lx2 + Lz2 + Lx3 + Lz3 + Lz4 + Lx4 + Ly2 + Ly3 + Ly4
    M_temp = zeros(N,N,N)
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N])

    # Last case where we now have some `irregular-grid` operators operating on the
    # regular-spaced axis x and z

    # These operators are operating on the regular grid x and z, but are constructed as though
    # they were irregular grid operators. Hence we test if we can seperate irregular and
    # regular gird operators on the same axis
    # Lx2 has 0 boundary points
    _Lx2 = CenteredDifference{1}(2,2,dx_vec,N)
    # Lz3 has 1 boundary point
    _Lz3 = CenteredDifference{3}(3,3,dz_vec,N)
    # Lx4 has 2 boundary points
    _Lx4 = 12.1*CenteredDifference{1}(4,4,dx_vec,N)

    A += _Lx2 + _Lz3 + _Lx4
    mul!(M_temp, A, M)
    @test M_temp ≈ ((Lx2*M)[1:N,2:N+1,2:N+1]+(Ly2*M)[2:N+1,1:N,2:N+1]+(Lz2*M)[2:N+1,2:N+1,1:N]+(Lx3*M)[1:N,2:N+1,2:N+1] +(Ly3*M)[2:N+1,1:N,2:N+1]
     + (Lz3*M)[2:N+1,2:N+1,1:N] + (Lx4*M)[1:N,2:N+1,2:N+1] +(Ly4*M)[2:N+1,1:N,2:N+1]+(Lz4*M)[2:N+1,2:N+1,1:N] + (_Lx2*M)[1:N,2:N+1,2:N+1]+(_Lz3*M)[2:N+1,2:N+1,1:N]+(_Lx4*M)[1:N,2:N+1,2:N+1])

end


text:
__precompile__()
module EC

using DelimitedFiles,DataFrames

include("./data.jl")
include("./io.jl")

export
    ec_grab,
    ec_list,
    cycle,
    total_cycles,
    get_measurement,
    get_params,
    get_curr_range,
    get_eq_anodickathodiccurrents,
    get_ind_eq_anodickathodicpotentials
end

ext: jl

max_stars_repo_name: lambcount/EC

max_issues_repo_name: lambcount/EC

max_forks_repo_name: lambcount/EC


text:
"""
Constructs a `DataFrames.DataFrame` from `layer`.

Where each `field` is a column, and each `feature` is a row. The name of
the geometry column (if any) defaults to `geometry0`.
"""
function geodataframe(layer::AG.FeatureLayer)
    sink = DataStreams.Data.stream!(AG.Source(layer), DataFrames.DataFrame)
    DataStreams.Data.close!(sink)
end

"""
Returns the dataset in `filename` as a `DataFrames.DataFrame`.

For details, see http://www.gdal.org/ogr_sql.html.

### Arguments
* `filename`: the dataset to run the query on
* (optional) `query`: the SQL SELECT query to be run on the dataset/layer

### Keyword Arguments
* `layer::Int` (default=`0`): the layer in the dataset to retrieve
* `dialect::String` (default=`SQLite`): the SQL dialect for `query` (if any)
"""
function read(filename::AbstractString; layer::Int=0)
    AG.registerdrivers() do
        AG.read(filename) do dataset
            geodataframe(AG.getlayer(dataset, layer))
        end
    end
end

function read(
        filename::AbstractString,
        query::AbstractString;
        dialect::String = "SQLite",
        layer::Int=0
    )
    AG.registerdrivers() do
        AG.read(filename) do dataset
            AG.executesql(dataset, query, dialect=dialect) do results
                geodataframe(results)
            end
        end
    end
end


text:
module ComputerVisionWithJuliaBook
using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, "..")))

include("CVUI.jl")

function julia_main()::Cint
    # do something based on ARGS?
    CVUI.initUI()
    return 0 # if things finished successfully
end
#comment when building native app
#julia_main()
end

max_stars_repo_name: bemcho/ComputerVisionWithJuliaBook

max_issues_repo_name: bemcho/ComputerVisionWithJuliaBook

max_forks_repo_name: bemcho/ComputerVisionWithJuliaBook


text:
"""
DuckDB data chunk
"""
mutable struct DataChunk
    handle::duckdb_data_chunk

    function DataChunk(handle::duckdb_data_chunk, destroy::Bool)
        result = new(handle)
        if destroy
            finalizer(_destroy_data_chunk, result)
        end
        return result
    end
end

function get_column_count(chunk::DataChunk)
    return duckdb_data_chunk_get_column_count(chunk.handle)
end

function get_size(chunk::DataChunk)
    return duckdb_data_chunk_get_size(chunk.handle)
end

function set_size(chunk::DataChunk, size::Int64)
    return duckdb_data_chunk_set_size(chunk.handle, size)
end

function get_vector(chunk::DataChunk, col_idx::Int64)::Vec
    if col_idx < 1 || col_idx > get_column_count(chunk)
        throw(
            InvalidInputException(
                string(
                    "get_array column index ",
                    col_idx,
                    " out of range, expected value between 1 and ",
                    get_column_count(chunk)
                )
            )
        )
    end
    return Vec(duckdb_data_chunk_get_vector(chunk.handle, col_idx))
end

function get_array(chunk::DataChunk, col_idx::Int64, ::Type{T})::Vector{T} where {T}
    return get_array(get_vector(chunk, col_idx), T)
end

function get_validity(chunk::DataChunk, col_idx::Int64)::ValidityMask
    return get_validity(get_vector(chunk, col_idx))
end

function all_valid(chunk::DataChunk, col_idx::Int64)
    return all_valid(get_vector(chunk, col_idx))
end

# this is only required when we own the data chunk
function _destroy_data_chunk(chunk::DataChunk)
    if chunk.handle != C_NULL
        duckdb_destroy_data_chunk(chunk.handle)
    end
    return chunk.handle = C_NULL
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


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

@test integrate(x^100*SpecialFunctions.gamma(0,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(101,a*x)+1//101*x^101*SpecialFunctions.gamma(0,a*x))
@test integrate(x^2*SpecialFunctions.gamma(0,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(3,a*x)+1//3*x^3*SpecialFunctions.gamma(0,a*x))
@test integrate(x*SpecialFunctions.gamma(0,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(0,a*x)+-1//2*a^-2*SpecialFunctions.gamma(2,a*x))
@test integrate(SpecialFunctions.gamma(0,a*x), x) == :(x*SpecialFunctions.gamma(0,a*x)+-1*a^-1*exp(-1*a*x))
@test integrate(x^-2*SpecialFunctions.gamma(0,a*x), x) == :(a*SpecialFunctions.gamma(-1,a*x)+-1*x^-1*SpecialFunctions.gamma(0,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(0,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(-2,a*x)+-1//2*x^-2*SpecialFunctions.gamma(0,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(0,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(0,a*x)+1//3*a^3*SpecialFunctions.gamma(-3,a*x))
@test integrate(x^3*SpecialFunctions.gamma(1,a*x), x) == :(-6*a^-4*exp(-1*a*x)+-1*a^-1*x^3*exp(-1*a*x)+-6*x*a^-3*exp(-1*a*x)+-3*a^-2*x^2*exp(-1*a*x))
@test integrate(x^2*SpecialFunctions.gamma(1,a*x), x) == :(-2*a^-3*exp(-1*a*x)+-1*a^-1*x^2*exp(-1*a*x)+-2*x*a^-2*exp(-1*a*x))
@test integrate(x*SpecialFunctions.gamma(1,a*x), x) == :(-1*a^-2*exp(-1*a*x)+-1*x*a^-1*exp(-1*a*x))
@test integrate(SpecialFunctions.gamma(1,a*x), x) == :(-1*a^-1*exp(-1*a*x))
@test integrate(x^-1*SpecialFunctions.gamma(1,a*x), x) == :(expintegralei(-1*a*x))
@test integrate(x^-2*SpecialFunctions.gamma(1,a*x), x) == :(-1*a*expintegralei(-1*a*x)+-1*x^-1*exp(-1*a*x))
@test integrate(x^-3*SpecialFunctions.gamma(1,a*x), x) == :((1/2)*a^2*expintegralei(-1*a*x)+-1//2*x^-2*exp(-1*a*x)+(1/2)*a*x^-1*exp(-1*a*x))
@test integrate(x^-4*SpecialFunctions.gamma(1,a*x), x) == :(-1//3*x^-3*exp(-1*a*x)+-1//6*a^3*expintegralei(-1*a*x)+-1//6*a^2*x^-1*exp(-1*a*x)+1//6*a*x^-2*exp(-1*a*x))
@test integrate(x^100*SpecialFunctions.gamma(2,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(103,a*x)+1//101*x^101*SpecialFunctions.gamma(2,a*x))
@test integrate(x^2*SpecialFunctions.gamma(2,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(5,a*x)+1//3*x^3*SpecialFunctions.gamma(2,a*x))
@test integrate(x*SpecialFunctions.gamma(2,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(2,a*x)+-1//2*a^-2*SpecialFunctions.gamma(4,a*x))
@test integrate(SpecialFunctions.gamma(2,a*x), x) == :(x*SpecialFunctions.gamma(2,a*x)+-1*a^-1*SpecialFunctions.gamma(3,a*x))
@test integrate(x^-1*SpecialFunctions.gamma(2,a*x), x) == :(-1*exp(-1*a*x)+expintegralei(-1*a*x))
@test integrate(x^-2*SpecialFunctions.gamma(2,a*x), x) == :(a*exp(-1*a*x)+-1*x^-1*SpecialFunctions.gamma(2,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(2,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(0,a*x)+-1//2*x^-2*SpecialFunctions.gamma(2,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(2,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(2,a*x)+1//3*a^3*SpecialFunctions.gamma(-1,a*x))
@test integrate(x^100*SpecialFunctions.gamma(3,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(104,a*x)+1//101*x^101*SpecialFunctions.gamma(3,a*x))
@test integrate(x^2*SpecialFunctions.gamma(3,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(6,a*x)+1//3*x^3*SpecialFunctions.gamma(3,a*x))
@test integrate(x*SpecialFunctions.gamma(3,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(3,a*x)+-1//2*a^-2*SpecialFunctions.gamma(5,a*x))
@test integrate(SpecialFunctions.gamma(3,a*x), x) == :(x*SpecialFunctions.gamma(3,a*x)+-1*a^-1*SpecialFunctions.gamma(4,a*x))
@test integrate(x^-1*SpecialFunctions.gamma(3,a*x), x) == :(-1*SpecialFunctions.gamma(2,a*x)+-2*exp(-1*a*x)+2*expintegralei(-1*a*x))
@test integrate(x^-2*SpecialFunctions.gamma(3,a*x), x) == :(a*SpecialFunctions.gamma(2,a*x)+-1*x^-1*SpecialFunctions.gamma(3,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(3,a*x), x) == :((1/2)*a^2*exp(-1*a*x)+-1//2*x^-2*SpecialFunctions.gamma(3,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(3,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(3,a*x)+1//3*a^3*SpecialFunctions.gamma(0,a*x))
@test integrate(x^100*SpecialFunctions.gamma(-1,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(100,a*x)+1//101*x^101*SpecialFunctions.gamma(-1,a*x))
@test integrate(x^3*SpecialFunctions.gamma(-1,a*x), x) == :(-1//4*a^-4*SpecialFunctions.gamma(3,a*x)+1//4*x^4*SpecialFunctions.gamma(-1,a*x))
@test integrate(x^2*SpecialFunctions.gamma(-1,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(2,a*x)+1//3*x^3*SpecialFunctions.gamma(-1,a*x))
@test integrate(x*SpecialFunctions.gamma(-1,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(-1,a*x)+-1//2*a^-2*exp(-1*a*x))
@test integrate(SpecialFunctions.gamma(-1,a*x), x) == :(x*SpecialFunctions.gamma(-1,a*x)+-1*a^-1*SpecialFunctions.gamma(0,a*x))
@test integrate(x^-2*SpecialFunctions.gamma(-1,a*x), x) == :(a*SpecialFunctions.gamma(-2,a*x)+-1*x^-1*SpecialFunctions.gamma(-1,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(-1,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(-3,a*x)+-1//2*x^-2*SpecialFunctions.gamma(-1,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(-1,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(-1,a*x)+1//3*a^3*SpecialFunctions.gamma(-4,a*x))
@test integrate(x^100*SpecialFunctions.gamma(-2,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(99,a*x)+1//101*x^101*SpecialFunctions.gamma(-2,a*x))
@test integrate(x^3*SpecialFunctions.gamma(-2,a*x), x) == :(-1//4*a^-4*SpecialFunctions.gamma(2,a*x)+1//4*x^4*SpecialFunctions.gamma(-2,a*x))
@test integrate(x^2*SpecialFunctions.gamma(-2,a*x), x) == :(-1//3*a^-3*exp(-1*a*x)+1//3*x^3*SpecialFunctions.gamma(-2,a*x))
@test integrate(x*SpecialFunctions.gamma(-2,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(-2,a*x)+-1//2*a^-2*SpecialFunctions.gamma(0,a*x))
@test integrate(SpecialFunctions.gamma(-2,a*x), x) == :(x*SpecialFunctions.gamma(-2,a*x)+-1*a^-1*SpecialFunctions.gamma(-1,a*x))
@test integrate(x^-2*SpecialFunctions.gamma(-2,a*x), x) == :(a*SpecialFunctions.gamma(-3,a*x)+-1*x^-1*SpecialFunctions.gamma(-2,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(-2,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(-4,a*x)+-1//2*x^-2*SpecialFunctions.gamma(-2,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(-2,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(-2,a*x)+1//3*a^3*SpecialFunctions.gamma(-5,a*x))
@test integrate(x^100*SpecialFunctions.gamma(-3,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(98,a*x)+1//101*x^101*SpecialFunctions.gamma(-3,a*x))
@test integrate(x^3*SpecialFunctions.gamma(-3,a*x), x) == :(-1//4*a^-4*exp(-1*a*x)+1//4*x^4*SpecialFunctions.gamma(-3,a*x))
@test integrate(x^2*SpecialFunctions.gamma(-3,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(0,a*x)+1//3*x^3*SpecialFunctions.gamma(-3,a*x))
@test integrate(x*SpecialFunctions.gamma(-3,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(-3,a*x)+-1//2*a^-2*SpecialFunctions.gamma(-1,a*x))
@test integrate(SpecialFunctions.gamma(-3,a*x), x) == :(x*SpecialFunctions.gamma(-3,a*x)+-1*a^-1*SpecialFunctions.gamma(-2,a*x))
@test integrate(x^-2*SpecialFunctions.gamma(-3,a*x), x) == :(a*SpecialFunctions.gamma(-4,a*x)+-1*x^-1*SpecialFunctions.gamma(-3,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(-3,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(-5,a*x)+-1//2*x^-2*SpecialFunctions.gamma(-3,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(-3,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(-3,a*x)+1//3*a^3*SpecialFunctions.gamma(-6,a*x))
@test integrate(x^100*SpecialFunctions.gamma(1/2,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(203//2,a*x)+1//101*x^101*SpecialFunctions.gamma(1/2,a*x))
@test integrate(x^2*SpecialFunctions.gamma(1/2,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(7//2,a*x)+1//3*x^3*SpecialFunctions.gamma(1/2,a*x))
@test integrate(x*SpecialFunctions.gamma(1/2,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(1/2,a*x)+-1//2*a^-2*SpecialFunctions.gamma(5//2,a*x))
@test integrate(SpecialFunctions.gamma(1/2,a*x), x) == :(x*SpecialFunctions.gamma(1/2,a*x)+-1*a^-1*SpecialFunctions.gamma(3//2,a*x))
@test integrate(x^-2*SpecialFunctions.gamma(1/2,a*x), x) == :(a*SpecialFunctions.gamma(-1//2,a*x)+-1*x^-1*SpecialFunctions.gamma(1/2,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(1/2,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(-3//2,a*x)+-1//2*x^-2*SpecialFunctions.gamma(1/2,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(1/2,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(1/2,a*x)+1//3*a^3*SpecialFunctions.gamma(-5//2,a*x))
@test integrate(x^100*SpecialFunctions.gamma(3//2,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(205//2,a*x)+1//101*x^101*SpecialFunctions.gamma(3//2,a*x))
@test integrate(x^2*SpecialFunctions.gamma(3//2,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(9//2,a*x)+1//3*x^3*SpecialFunctions.gamma(3//2,a*x))
@test integrate(x*SpecialFunctions.gamma(3//2,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(3//2,a*x)+-1//2*a^-2*SpecialFunctions.gamma(7//2,a*x))
@test integrate(SpecialFunctions.gamma(3//2,a*x), x) == :(x*SpecialFunctions.gamma(3//2,a*x)+-1*a^-1*SpecialFunctions.gamma(5//2,a*x))
@test integrate(x^-2*SpecialFunctions.gamma(3//2,a*x), x) == :(a*SpecialFunctions.gamma(1/2,a*x)+-1*x^-1*SpecialFunctions.gamma(3//2,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(3//2,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(-1//2,a*x)+-1//2*x^-2*SpecialFunctions.gamma(3//2,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(3//2,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(3//2,a*x)+1//3*a^3*SpecialFunctions.gamma(-3//2,a*x))
@test integrate((d*x)^m*SpecialFunctions.gamma(3,b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(3,b*x)+-1*b^-1*(b*x)^(-1m)*(d*x)^m*(1+m)^-1*SpecialFunctions.gamma(4+m,b*x))
@test integrate((d*x)^m*SpecialFunctions.gamma(2,b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(2,b*x)+-1*b^-1*(b*x)^(-1m)*(d*x)^m*(1+m)^-1*SpecialFunctions.gamma(3+m,b*x))
@test integrate((d*x)^m*SpecialFunctions.gamma(1,b*x), x) == :(-1*b^-1*(b*x)^(-1m)*(d*x)^m*SpecialFunctions.gamma(1+m,b*x))
@test integrate((d*x)^m*SpecialFunctions.gamma(0,b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(0,b*x)+-1*b^-1*(b*x)^(-1m)*(d*x)^m*(1+m)^-1*SpecialFunctions.gamma(1+m,b*x))
@test integrate((d*x)^m*SpecialFunctions.gamma(-1,b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(-1,b*x)+-1*b^-1*(b*x)^(-1m)*(d*x)^m*(1+m)^-1*SpecialFunctions.gamma(m,b*x))
@test integrate((d*x)^m*SpecialFunctions.gamma(-2,b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(-2,b*x)+-1*b^-1*(b*x)^(-1m)*(d*x)^m*(1+m)^-1*SpecialFunctions.gamma(-1+m,b*x))
@test integrate(x^m*SpecialFunctions.gamma(n,x), x) == :(-1*(1+m)^-1*SpecialFunctions.gamma(1+m+n,x)+x^(1+m)*(1+m)^-1*SpecialFunctions.gamma(n,x))
@test integrate(x^m*SpecialFunctions.gamma(n,b*x), x) == :(x^(1+m)*(1+m)^-1*SpecialFunctions.gamma(n,b*x)+-1*b^-1*x^m*(b*x)^(-1m)*(1+m)^-1*SpecialFunctions.gamma(1+m+n,b*x))
@test integrate((d*x)^m*SpecialFunctions.gamma(n,x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(n,x)+-1*x^(-1m)*(d*x)^m*(1+m)^-1*SpecialFunctions.gamma(1+m+n,x))
@test integrate((d*x)^m*SpecialFunctions.gamma(n,b*x), x) == :(d^-1*(d*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(n,b*x)+-1*b^-1*(b*x)^(-1m)*(d*x)^m*(1+m)^-1*SpecialFunctions.gamma(1+m+n,b*x))
@test integrate(x^100*SpecialFunctions.gamma(n,a*x), x) == :(-1//101*a^-101*SpecialFunctions.gamma(101+n,a*x)+1//101*x^101*SpecialFunctions.gamma(n,a*x))
@test integrate(x^2*SpecialFunctions.gamma(n,a*x), x) == :(-1//3*a^-3*SpecialFunctions.gamma(3+n,a*x)+1//3*x^3*SpecialFunctions.gamma(n,a*x))
@test integrate(x*SpecialFunctions.gamma(n,a*x), x) == :((1/2)*x^2*SpecialFunctions.gamma(n,a*x)+-1//2*a^-2*SpecialFunctions.gamma(2+n,a*x))
@test integrate(SpecialFunctions.gamma(n,a*x), x) == :(x*SpecialFunctions.gamma(n,a*x)+-1*a^-1*SpecialFunctions.gamma(1+n,a*x))
@test integrate(x^-2*SpecialFunctions.gamma(n,a*x), x) == :(a*SpecialFunctions.gamma(-1+n,a*x)+-1*x^-1*SpecialFunctions.gamma(n,a*x))
@test integrate(x^-3*SpecialFunctions.gamma(n,a*x), x) == :((1/2)*a^2*SpecialFunctions.gamma(-2+n,a*x)+-1//2*x^-2*SpecialFunctions.gamma(n,a*x))
@test integrate(x^-4*SpecialFunctions.gamma(n,a*x), x) == :(-1//3*x^-3*SpecialFunctions.gamma(n,a*x)+1//3*a^3*SpecialFunctions.gamma(-3+n,a*x))
@test integrate(x^2*SpecialFunctions.gamma(n,2x), x) == :(-1//24*SpecialFunctions.gamma(3+n,2x)+1//3*x^3*SpecialFunctions.gamma(n,2x))
@test integrate(x*SpecialFunctions.gamma(n,2x), x) == :(-1//8*SpecialFunctions.gamma(2+n,2x)+(1/2)*x^2*SpecialFunctions.gamma(n,2x))
@test integrate(SpecialFunctions.gamma(n,2x), x) == :(-1//2*SpecialFunctions.gamma(1+n,2x)+x*SpecialFunctions.gamma(n,2x))
@test integrate(x^-2*SpecialFunctions.gamma(n,2x), x) == :(2*SpecialFunctions.gamma(-1+n,2x)+-1*x^-1*SpecialFunctions.gamma(n,2x))
@test integrate(x^-3*SpecialFunctions.gamma(n,2x), x) == :(2*SpecialFunctions.gamma(-2+n,2x)+-1//2*x^-2*SpecialFunctions.gamma(n,2x))
@test integrate(x^-4*SpecialFunctions.gamma(n,2x), x) == :(8//3*SpecialFunctions.gamma(-3+n,2x)+-1//3*x^-3*SpecialFunctions.gamma(n,2x))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(0,a+b*x), x) == :(-1//4*b^-4*(b*c+-1*a*d)^3*exp(-1a+-1*b*x)+1//4*d^-1*(c+d*x)^4*SpecialFunctions.gamma(0,a+b*x)+-1//4*b^-4*d^-1*(b*c+-1*a*d)^4*SpecialFunctions.gamma(0,a+b*x)+-1//4*b^-4*d^3*SpecialFunctions.gamma(4,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//4*d*b^-4*(b*c+-1*a*d)^2*SpecialFunctions.gamma(2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//4*b^-4*d^2*(b*c+-1*a*d)*SpecialFunctions.gamma(3,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(0,a+b*x), x) == :(-1//3*b^-3*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+1//3*d^-1*(c+d*x)^3*SpecialFunctions.gamma(0,a+b*x)+-1//3*b^-3*d^-1*(b*c+-1*a*d)^3*SpecialFunctions.gamma(0,a+b*x)+-1//3*b^-3*d^2*SpecialFunctions.gamma(3,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//3*d*b^-3*(b*c+-1*a*d)*SpecialFunctions.gamma(2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)*SpecialFunctions.gamma(0,a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*SpecialFunctions.gamma(0,a+b*x)+-1//2*b^-2*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-1//2*d*b^-2*SpecialFunctions.gamma(2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*SpecialFunctions.gamma(0,a+b*x))
@test integrate(SpecialFunctions.gamma(0,a+b*x), x) == :(-1*b^-1*exp(-1a+-1*b*x)+b^-1*(a+b*x)*SpecialFunctions.gamma(0,a+b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.gamma(0,a+b*x), x) == :(-1*d^-1*(c+d*x)^-1*SpecialFunctions.gamma(0,a+b*x)+b*d^-1*(b*c+-1*a*d)^-1*SpecialFunctions.gamma(0,a+b*x)+-1*b*d^-1*(b*c+-1*a*d)^-1*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-3*SpecialFunctions.gamma(0,a+b*x), x) == :(-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.gamma(0,a+b*x)+(1/2)*b^2*d^-1*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(0,a+b*x)+-1//2*b^2*d^-1*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//2*b^2*d^-2*(b*c+-1*a*d)^-1*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-4*SpecialFunctions.gamma(0,a+b*x), x) == :(-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.gamma(0,a+b*x)+1//3*b^3*d^-1*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(0,a+b*x)+-1//3*b^3*d^-1*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//3*b^3*d^-3*(b*c+-1*a*d)^-1*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//3*b^3*d^-2*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^4*SpecialFunctions.gamma(1,a+b*x), x) == :(-1*b^-1*(c+d*x)^4*exp(-1a+-1*b*x)+-24*b^-5*d^4*exp(-1a+-1*b*x)+-24*b^-4*d^3*(c+d*x)*exp(-1a+-1*b*x)+-12*b^-3*d^2*(c+d*x)^2*exp(-1a+-1*b*x)+-4*d*b^-2*(c+d*x)^3*exp(-1a+-1*b*x))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(1,a+b*x), x) == :(-1*b^-1*(c+d*x)^3*exp(-1a+-1*b*x)+-6*b^-4*d^3*exp(-1a+-1*b*x)+-6*b^-3*d^2*(c+d*x)*exp(-1a+-1*b*x)+-3*d*b^-2*(c+d*x)^2*exp(-1a+-1*b*x))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(1,a+b*x), x) == :(-1*b^-1*(c+d*x)^2*exp(-1a+-1*b*x)+-2*b^-3*d^2*exp(-1a+-1*b*x)+-2*d*b^-2*(c+d*x)*exp(-1a+-1*b*x))
@test integrate((c+d*x)*SpecialFunctions.gamma(1,a+b*x), x) == :(-1*d*b^-2*exp(-1a+-1*b*x)+-1*b^-1*(c+d*x)*exp(-1a+-1*b*x))
@test integrate(SpecialFunctions.gamma(1,a+b*x), x) == :(-1*b^-1*exp(-1a+-1*b*x))
@test integrate((c+d*x)^-1*SpecialFunctions.gamma(1,a+b*x), x) == :(d^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-2*SpecialFunctions.gamma(1,a+b*x), x) == :(-1*d^-1*(c+d*x)^-1*exp(-1a+-1*b*x)+-1*b*d^-2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-3*SpecialFunctions.gamma(1,a+b*x), x) == :(-1//2*d^-1*(c+d*x)^-2*exp(-1a+-1*b*x)+(1/2)*b*d^-2*(c+d*x)^-1*exp(-1a+-1*b*x)+(1/2)*b^2*d^-3*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-4*SpecialFunctions.gamma(1,a+b*x), x) == :(-1//3*d^-1*(c+d*x)^-3*exp(-1a+-1*b*x)+-1//6*b^2*d^-3*(c+d*x)^-1*exp(-1a+-1*b*x)+-1//6*b^3*d^-4*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+1//6*b*d^-2*(c+d*x)^-2*exp(-1a+-1*b*x))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(2,a+b*x), x) == :(1//4*d^-1*(c+d*x)^4*SpecialFunctions.gamma(2,a+b*x)+-1//4*b^-4*d^3*SpecialFunctions.gamma(6,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+1//4*b^-4*d^2*(b*c+-1*a*d)*SpecialFunctions.gamma(5,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(2,a+b*x), x) == :(1//3*d^-1*(c+d*x)^3*SpecialFunctions.gamma(2,a+b*x)+-1//3*b^-3*d^2*SpecialFunctions.gamma(5,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+1//3*d*b^-3*(b*c+-1*a*d)*SpecialFunctions.gamma(4,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)*SpecialFunctions.gamma(2,a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*SpecialFunctions.gamma(2,a+b*x)+(1/2)*b^-2*(b*c+-1*a*d)*SpecialFunctions.gamma(3,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//2*d*b^-2*SpecialFunctions.gamma(4,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate(SpecialFunctions.gamma(2,a+b*x), x) == :(-1*b^-1*SpecialFunctions.gamma(3,a+b*x)+b^-1*(a+b*x)*SpecialFunctions.gamma(2,a+b*x))
@test integrate((c+d*x)^-1*SpecialFunctions.gamma(2,a+b*x), x) == :(-1*d^-1*exp(-1a+-1*b*x)+d^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1*d^-2*(b*c+-1*a*d)*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-2*SpecialFunctions.gamma(2,a+b*x), x) == :(b*d^-2*exp(-1a+-1*b*x)+-1*d^-1*(c+d*x)^-1*SpecialFunctions.gamma(2,a+b*x)+-1*b*d^-3*(b*c+-1*a*d)*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-3*SpecialFunctions.gamma(2,a+b*x), x) == :(-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.gamma(2,a+b*x)+(1/2)*b^2*d^-3*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//2*b^2*d^-4*(b*c+-1*a*d)*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-4*SpecialFunctions.gamma(2,a+b*x), x) == :(-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.gamma(2,a+b*x)+1//3*b^3*d^-4*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//3*b^3*d^-5*(b*c+-1*a*d)*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-5*SpecialFunctions.gamma(2,a+b*x), x) == :(-1//4*d^-1*(c+d*x)^-4*SpecialFunctions.gamma(2,a+b*x)+1//4*b^4*d^-5*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//4*b^4*d^-6*(b*c+-1*a*d)*SpecialFunctions.gamma(-3,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(3,a+b*x), x) == :(1//4*d^-1*(c+d*x)^4*SpecialFunctions.gamma(3,a+b*x)+-1//4*b^-4*d^3*SpecialFunctions.gamma(7,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+(1/2)*b^-4*d^2*(b*c+-1*a*d)*SpecialFunctions.gamma(6,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//4*d*b^-4*(b*c+-1*a*d)^2*SpecialFunctions.gamma(5,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(3,a+b*x), x) == :(1//3*d^-1*(c+d*x)^3*SpecialFunctions.gamma(3,a+b*x)+-1//3*b^-3*d^2*SpecialFunctions.gamma(6,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//3*b^-3*(b*c+-1*a*d)^2*SpecialFunctions.gamma(4,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+2//3*d*b^-3*(b*c+-1*a*d)*SpecialFunctions.gamma(5,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)*SpecialFunctions.gamma(3,a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*SpecialFunctions.gamma(3,a+b*x)+-1*b^-2*(b*c+-1*a*d)*SpecialFunctions.gamma(4,a+b*x)+-1//2*d*b^-2*SpecialFunctions.gamma(5,a+b*x)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*SpecialFunctions.gamma(3,a+b*x))
@test integrate(SpecialFunctions.gamma(3,a+b*x), x) == :(-1*b^-1*SpecialFunctions.gamma(4,a+b*x)+b^-1*(a+b*x)*SpecialFunctions.gamma(3,a+b*x))
@test integrate((c+d*x)^-1*SpecialFunctions.gamma(3,a+b*x), x) == :(-3*d^-1*exp(-1a+-1*b*x)+d^-2*(b*c+-1*a*d)*exp(-1a+-1*b*x)+-1*d^-1*(a+b*x)*exp(-1a+-1*b*x)+2*d^-1*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+d^-3*(b*c+-1*a*d)^2*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x))+-1*d^-2*(-2*a*d+2*b*c)*exp(-1a+b*c*d^-1)*expintegralei(-1*b*d^-1*(c+d*x)))
@test integrate((c+d*x)^-2*SpecialFunctions.gamma(3,a+b*x), x) == :(b*d^-2*SpecialFunctions.gamma(2,a+b*x)+-1*d^-1*(c+d*x)^-1*SpecialFunctions.gamma(3,a+b*x)+-1*b*d^-3*(b*c+-1*a*d)*exp(-1a+-1*b*x)+b*d^-4*(b*c+-1*a*d)^2*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-3*SpecialFunctions.gamma(3,a+b*x), x) == :((1/2)*b^2*d^-3*exp(-1a+-1*b*x)+-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.gamma(3,a+b*x)+(1/2)*b^2*d^-5*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1*b^2*d^-4*(b*c+-1*a*d)*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-4*SpecialFunctions.gamma(3,a+b*x), x) == :(-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.gamma(3,a+b*x)+1//3*b^3*d^-4*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-2//3*b^3*d^-5*(b*c+-1*a*d)*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+1//3*b^3*d^-6*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-5*SpecialFunctions.gamma(3,a+b*x), x) == :(-1//4*d^-1*(c+d*x)^-4*SpecialFunctions.gamma(3,a+b*x)+1//4*b^4*d^-5*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//2*b^4*d^-6*(b*c+-1*a*d)*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+1//4*b^4*d^-7*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-3,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(-1,a+b*x), x) == :(-1*b^-4*(b*c+-1*a*d)^3*SpecialFunctions.gamma(0,a+b*x)+-1//4*b^-4*d^3*SpecialFunctions.gamma(3,a+b*x)+1//4*d^-1*(c+d*x)^4*SpecialFunctions.gamma(-1,a+b*x)+-1*b^-4*d^2*(b*c+-1*a*d)*SpecialFunctions.gamma(2,a+b*x)+-3//2*d*b^-4*(b*c+-1*a*d)^2*exp(-1a+-1*b*x)+-1//4*b^-4*d^-1*(b*c+-1*a*d)^4*SpecialFunctions.gamma(-1,a+b*x))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(-1,a+b*x), x) == :(-1*b^-3*(b*c+-1*a*d)^2*SpecialFunctions.gamma(0,a+b*x)+1//3*d^-1*(c+d*x)^3*SpecialFunctions.gamma(-1,a+b*x)+-1//3*d*b^-3*(-2*a*d+3*b*c)*exp(-1a+-1*b*x)+-1//3*b^-3*d^-1*(b*c+-1*a*d)^3*SpecialFunctions.gamma(-1,a+b*x)+-1//3*b^-3*d^2*SpecialFunctions.gamma(2,b*x)*exp(-1a))
@test integrate((c+d*x)*SpecialFunctions.gamma(-1,a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*SpecialFunctions.gamma(-1,a+b*x)+-1*b^-2*(b*c+-1*a*d)*SpecialFunctions.gamma(0,a+b*x)+-1//2*d*b^-2*exp(-1a+-1*b*x)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-1,a+b*x))
@test integrate(SpecialFunctions.gamma(-1,a+b*x), x) == :(-1*b^-1*SpecialFunctions.gamma(0,a+b*x)+b^-1*(a+b*x)*SpecialFunctions.gamma(-1,a+b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.gamma(-1,a+b*x), x) == :(-1*b*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(0,a+b*x)+-1*d^-1*(c+d*x)^-1*SpecialFunctions.gamma(-1,a+b*x)+b*d^-1*(b*c+-1*a*d)^-1*SpecialFunctions.gamma(-1,a+b*x)+b*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-3*SpecialFunctions.gamma(-1,a+b*x), x) == :(-1*b^2*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(0,a+b*x)+-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.gamma(-1,a+b*x)+b^2*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+(1/2)*b^2*d^-1*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-1,a+b*x)+(1/2)*b^2*d^-1*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-4*SpecialFunctions.gamma(-1,a+b*x), x) == :(-1*b^3*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(0,a+b*x)+-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.gamma(-1,a+b*x)+b^3*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+1//3*b^3*d^-1*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-1,a+b*x)+1//3*b^3*d^-2*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+2//3*b^3*d^-1*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(-2,a+b*x), x) == :(-1*b^-4*(b*c+-1*a*d)^3*SpecialFunctions.gamma(-1,a+b*x)+1//4*d^-1*(c+d*x)^4*SpecialFunctions.gamma(-2,a+b*x)+-3//2*d*b^-4*(b*c+-1*a*d)^2*SpecialFunctions.gamma(0,a+b*x)+-1//4*b^-4*d^-1*(b*c+-1*a*d)^4*SpecialFunctions.gamma(-2,a+b*x)+-1//4*b^-4*d^2*(-3*a*d+4*b*c)*exp(-1a+-1*b*x)+-1//4*b^-4*d^3*SpecialFunctions.gamma(2,b*x)*exp(-1a))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(-2,a+b*x), x) == :(-1*b^-3*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-1,a+b*x)+-1//3*b^-3*d^2*exp(-1a+-1*b*x)+1//3*d^-1*(c+d*x)^3*SpecialFunctions.gamma(-2,a+b*x)+-1*d*b^-3*(b*c+-1*a*d)*SpecialFunctions.gamma(0,a+b*x)+-1//3*b^-3*d^-1*(b*c+-1*a*d)^3*SpecialFunctions.gamma(-2,a+b*x))
@test integrate((c+d*x)*SpecialFunctions.gamma(-2,a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*SpecialFunctions.gamma(-2,a+b*x)+-1*b^-2*(b*c+-1*a*d)*SpecialFunctions.gamma(-1,a+b*x)+-1//2*d*b^-2*SpecialFunctions.gamma(0,a+b*x)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-2,a+b*x))
@test integrate(SpecialFunctions.gamma(-2,a+b*x), x) == :(-1*b^-1*SpecialFunctions.gamma(-1,a+b*x)+b^-1*(a+b*x)*SpecialFunctions.gamma(-2,a+b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.gamma(-2,a+b*x), x) == :(-1*b*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-1,a+b*x)+-1*d^-1*(c+d*x)^-1*SpecialFunctions.gamma(-2,a+b*x)+b*d*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(0,a+b*x)+b*d^-1*(b*c+-1*a*d)^-1*SpecialFunctions.gamma(-2,a+b*x)+-1*b*d*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-3*SpecialFunctions.gamma(-2,a+b*x), x) == :(-1*b^2*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-1,a+b*x)+-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.gamma(-2,a+b*x)+(1/2)*b^2*d^-1*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-2,a+b*x)+-1//2*b^2*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+3//2*d*b^2*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(0,a+b*x)+-3//2*d*b^2*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-4*SpecialFunctions.gamma(-2,a+b*x), x) == :(-1*b^3*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(-1,a+b*x)+-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.gamma(-2,a+b*x)+-1*b^3*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+2*d*b^3*(b*c+-1*a*d)^-5*SpecialFunctions.gamma(0,a+b*x)+1//3*b^3*d^-1*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-2,a+b*x)+-2*d*b^3*(b*c+-1*a*d)^-5*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1//3*b^3*d^-1*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(-3,a+b*x), x) == :(-1*b^-4*(b*c+-1*a*d)^3*SpecialFunctions.gamma(-2,a+b*x)+-1//4*b^-4*d^3*exp(-1a+-1*b*x)+1//4*d^-1*(c+d*x)^4*SpecialFunctions.gamma(-3,a+b*x)+-1*b^-4*d^2*(b*c+-1*a*d)*SpecialFunctions.gamma(0,a+b*x)+-3//2*d*b^-4*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-1,a+b*x)+-1//4*b^-4*d^-1*(b*c+-1*a*d)^4*SpecialFunctions.gamma(-3,a+b*x))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(-3,a+b*x), x) == :(-1*b^-3*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-2,a+b*x)+-1//3*b^-3*d^2*SpecialFunctions.gamma(0,a+b*x)+1//3*d^-1*(c+d*x)^3*SpecialFunctions.gamma(-3,a+b*x)+-1*d*b^-3*(b*c+-1*a*d)*SpecialFunctions.gamma(-1,a+b*x)+-1//3*b^-3*d^-1*(b*c+-1*a*d)^3*SpecialFunctions.gamma(-3,a+b*x))
@test integrate((c+d*x)*SpecialFunctions.gamma(-3,a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*SpecialFunctions.gamma(-3,a+b*x)+-1*b^-2*(b*c+-1*a*d)*SpecialFunctions.gamma(-2,a+b*x)+-1//2*d*b^-2*SpecialFunctions.gamma(-1,a+b*x)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*SpecialFunctions.gamma(-3,a+b*x))
@test integrate(SpecialFunctions.gamma(-3,a+b*x), x) == :(-1*b^-1*SpecialFunctions.gamma(-2,a+b*x)+b^-1*(a+b*x)*SpecialFunctions.gamma(-3,a+b*x))
@test integrate((c+d*x)^-2*SpecialFunctions.gamma(-3,a+b*x), x) == :(-1*b*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-2,a+b*x)+-1*d^-1*(c+d*x)^-1*SpecialFunctions.gamma(-3,a+b*x)+b*d*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-1,a+b*x)+b*d^-1*(b*c+-1*a*d)^-1*SpecialFunctions.gamma(-3,a+b*x)+-1*b*d^2*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(0,a+b*x)+b*d^2*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-3*SpecialFunctions.gamma(-3,a+b*x), x) == :(-1*b^2*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-2,a+b*x)+-1//2*d^-1*(c+d*x)^-2*SpecialFunctions.gamma(-3,a+b*x)+(1/2)*b^2*d^-1*(b*c+-1*a*d)^-2*SpecialFunctions.gamma(-3,a+b*x)+-2*b^2*d^2*(b*c+-1*a*d)^-5*SpecialFunctions.gamma(0,a+b*x)+3//2*d*b^2*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(-1,a+b*x)+(1/2)*d*b^2*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+2*b^2*d^2*(b*c+-1*a*d)^-5*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^-4*SpecialFunctions.gamma(-3,a+b*x), x) == :(-1*b^3*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(-2,a+b*x)+-1//3*d^-1*(c+d*x)^-3*SpecialFunctions.gamma(-3,a+b*x)+2*d*b^3*(b*c+-1*a*d)^-5*SpecialFunctions.gamma(-1,a+b*x)+-10//3*b^3*d^2*(b*c+-1*a*d)^-6*SpecialFunctions.gamma(0,a+b*x)+1//3*b^3*d^-1*(b*c+-1*a*d)^-3*SpecialFunctions.gamma(-3,a+b*x)+1//3*b^3*(b*c+-1*a*d)^-4*SpecialFunctions.gamma(-2,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+4//3*d*b^3*(b*c+-1*a*d)^-5*SpecialFunctions.gamma(-1,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+10//3*b^3*d^2*(b*c+-1*a*d)^-6*SpecialFunctions.gamma(0,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate(x^5//2*SpecialFunctions.gamma(2,a+b*x), x) == :(2//7*x^7//2*SpecialFunctions.gamma(2,a+b*x)+-2//7*b^-3*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(11//2,b*x)*exp(-1a)+-2//7*a*b^-3*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(9//2,b*x)*exp(-1a))
@test integrate(x^3//2*SpecialFunctions.gamma(2,a+b*x), x) == :(2//5*x^5//2*SpecialFunctions.gamma(2,a+b*x)+-2//5*b^-2*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(9//2,b*x)*exp(-1a)+-2//5*a*b^-2*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(7//2,b*x)*exp(-1a))
@test integrate(x^(1/2)*SpecialFunctions.gamma(2,a+b*x), x) == :(2//3*x^3//2*SpecialFunctions.gamma(2,a+b*x)+-2//3*b^-1*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(7//2,b*x)*exp(-1a)+-2//3*a*b^-1*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(5//2,b*x)*exp(-1a))
@test integrate(x^-1//2*SpecialFunctions.gamma(2,a+b*x), x) == :(2*x^(1/2)*SpecialFunctions.gamma(2,a+b*x)+-2*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(5//2,b*x)*exp(-1a)+-2*a*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(3//2,b*x)*exp(-1a))
@test integrate(x^-3//2*SpecialFunctions.gamma(2,a+b*x), x) == :(-2*x^-1//2*SpecialFunctions.gamma(2,a+b*x)+2*a*x^-1//2*(b*x)^(1/2)*SpecialFunctions.gamma(1/2,b*x)*exp(-1a)+2*b*x^(1/2)*(b*x)^-1//2*SpecialFunctions.gamma(3//2,b*x)*exp(-1a))
@test integrate(x^-5//2*SpecialFunctions.gamma(2,a+b*x), x) == :(-2//3*x^-3//2*SpecialFunctions.gamma(2,a+b*x)+2//3*b*x^-1//2*(b*x)^(1/2)*SpecialFunctions.gamma(1/2,b*x)*exp(-1a)+2//3*a*b*x^-1//2*(b*x)^(1/2)*SpecialFunctions.gamma(-1//2,b*x)*exp(-1a))
@test integrate(x^-7//2*SpecialFunctions.gamma(2,a+b*x), x) == :(-2//5*x^-5//2*SpecialFunctions.gamma(2,a+b*x)+2//5*b^2*x^-1//2*(b*x)^(1/2)*SpecialFunctions.gamma(-1//2,b*x)*exp(-1a)+2//5*a*b^2*x^-1//2*(b*x)^(1/2)*SpecialFunctions.gamma(-3//2,b*x)*exp(-1a))
@test integrate(x^-9//2*SpecialFunctions.gamma(2,a+b*x), x) == :(-2//7*x^-7//2*SpecialFunctions.gamma(2,a+b*x)+2//7*b^3*x^-1//2*(b*x)^(1/2)*SpecialFunctions.gamma(-3//2,b*x)*exp(-1a)+2//7*a*b^3*x^-1//2*(b*x)^(1/2)*SpecialFunctions.gamma(-5//2,b*x)*exp(-1a))
@test integrate(x^4//3*SpecialFunctions.gamma(2,a+b*x), x) == :(3//7*x^7//3*SpecialFunctions.gamma(2,a+b*x)+-3//7*b^-2*x^1//3*(b*x)^-1//3*SpecialFunctions.gamma(13//3,b*x)*exp(-1a)+-3//7*a*b^-2*x^1//3*(b*x)^-1//3*SpecialFunctions.gamma(10//3,b*x)*exp(-1a))
@test integrate(x^2//3*SpecialFunctions.gamma(2,a+b*x), x) == :(3//5*x^5//3*SpecialFunctions.gamma(2,a+b*x)+-3//5*b^-1*x^2//3*(b*x)^-2//3*SpecialFunctions.gamma(11//3,b*x)*exp(-1a)+-3//5*a*b^-1*x^2//3*(b*x)^-2//3*SpecialFunctions.gamma(8//3,b*x)*exp(-1a))
@test integrate(x^1//3*SpecialFunctions.gamma(2,a+b*x), x) == :(3//4*x^4//3*SpecialFunctions.gamma(2,a+b*x)+-3//4*b^-1*x^1//3*(b*x)^-1//3*SpecialFunctions.gamma(10//3,b*x)*exp(-1a)+-3//4*a*b^-1*x^1//3*(b*x)^-1//3*SpecialFunctions.gamma(7//3,b*x)*exp(-1a))
@test integrate(x^-1//3*SpecialFunctions.gamma(2,a+b*x), x) == :(3//2*x^2//3*SpecialFunctions.gamma(2,a+b*x)+-3//2*x^2//3*(b*x)^-2//3*SpecialFunctions.gamma(8//3,b*x)*exp(-1a)+-3//2*a*x^2//3*(b*x)^-2//3*SpecialFunctions.gamma(5//3,b*x)*exp(-1a))
@test integrate(x^-2//3*SpecialFunctions.gamma(2,a+b*x), x) == :(3*x^1//3*SpecialFunctions.gamma(2,a+b*x)+-3*x^1//3*(b*x)^-1//3*SpecialFunctions.gamma(7//3,b*x)*exp(-1a)+-3*a*x^1//3*(b*x)^-1//3*SpecialFunctions.gamma(4//3,b*x)*exp(-1a))
@test integrate(x^-4//3*SpecialFunctions.gamma(2,a+b*x), x) == :(-3*x^-1//3*SpecialFunctions.gamma(2,a+b*x)+3*a*x^-1//3*(b*x)^1//3*SpecialFunctions.gamma(2//3,b*x)*exp(-1a)+3*b*x^2//3*(b*x)^-2//3*SpecialFunctions.gamma(5//3,b*x)*exp(-1a))
@test integrate((c+d*x)^m*SpecialFunctions.gamma(3,a+b*x), x) == :(d^-1*(1+m)^-1*(c+d*x)^(1+m)*SpecialFunctions.gamma(3,a+b*x)+-1*b^-1*(b*d^-1*(c+d*x))^(-1m)*(1+m)^-1*(c+d*x)^m*SpecialFunctions.gamma(4+m,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+b^-1*d^-1*(b*d^-1*(c+d*x))^(-1m)*(1+m)^-1*(c+d*x)^m*(-2*a*d+2*b*c)*SpecialFunctions.gamma(3+m,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+-1*b^-1*d^-2*(b*d^-1*(c+d*x))^(-1m)*(1+m)^-1*(c+d*x)^m*(b*c+-1*a*d)^2*SpecialFunctions.gamma(2+m,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^m*SpecialFunctions.gamma(2,a+b*x), x) == :(d^-1*(1+m)^-1*(c+d*x)^(1+m)*SpecialFunctions.gamma(2,a+b*x)+-1*b^-1*(b*d^-1*(c+d*x))^(-1m)*(1+m)^-1*(c+d*x)^m*SpecialFunctions.gamma(3+m,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1)+b^-1*d^-1*(b*d^-1*(c+d*x))^(-1m)*(1+m)^-1*(c+d*x)^m*(b*c+-1*a*d)*SpecialFunctions.gamma(2+m,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^m*SpecialFunctions.gamma(1,a+b*x), x) == :(-1*b^-1*(b*d^-1*(c+d*x))^(-1m)*(c+d*x)^m*SpecialFunctions.gamma(1+m,b*d^-1*(c+d*x))*exp(-1a+b*c*d^-1))
@test integrate((c+d*x)^4*SpecialFunctions.gamma(n,a+b*x), x) == :(-1*b^-5*(b*c+-1*a*d)^4*SpecialFunctions.gamma(1+n,a+b*x)+-1//5*b^-5*d^4*SpecialFunctions.gamma(5+n,a+b*x)+1//5*d^-1*(c+d*x)^5*SpecialFunctions.gamma(n,a+b*x)+-1*b^-5*d^3*(b*c+-1*a*d)*SpecialFunctions.gamma(4+n,a+b*x)+-2*d*b^-5*(b*c+-1*a*d)^3*SpecialFunctions.gamma(2+n,a+b*x)+-2*b^-5*d^2*(b*c+-1*a*d)^2*SpecialFunctions.gamma(3+n,a+b*x)+-1//5*b^-5*d^-1*(b*c+-1*a*d)^5*SpecialFunctions.gamma(n,a+b*x))
@test integrate((c+d*x)^3*SpecialFunctions.gamma(n,a+b*x), x) == :(-1*b^-4*(b*c+-1*a*d)^3*SpecialFunctions.gamma(1+n,a+b*x)+-1//4*b^-4*d^3*SpecialFunctions.gamma(4+n,a+b*x)+1//4*d^-1*(c+d*x)^4*SpecialFunctions.gamma(n,a+b*x)+-1*b^-4*d^2*(b*c+-1*a*d)*SpecialFunctions.gamma(3+n,a+b*x)+-3//2*d*b^-4*(b*c+-1*a*d)^2*SpecialFunctions.gamma(2+n,a+b*x)+-1//4*b^-4*d^-1*(b*c+-1*a*d)^4*SpecialFunctions.gamma(n,a+b*x))
@test integrate((c+d*x)^2*SpecialFunctions.gamma(n,a+b*x), x) == :(-1*b^-3*(b*c+-1*a*d)^2*SpecialFunctions.gamma(1+n,a+b*x)+-1//3*b^-3*d^2*SpecialFunctions.gamma(3+n,a+b*x)+1//3*d^-1*(c+d*x)^3*SpecialFunctions.gamma(n,a+b*x)+-1*d*b^-3*(b*c+-1*a*d)*SpecialFunctions.gamma(2+n,a+b*x)+-1//3*b^-3*d^-1*(b*c+-1*a*d)^3*SpecialFunctions.gamma(n,a+b*x))
@test integrate((c+d*x)*SpecialFunctions.gamma(n,a+b*x), x) == :((1/2)*d^-1*(c+d*x)^2*SpecialFunctions.gamma(n,a+b*x)+-1*b^-2*(b*c+-1*a*d)*SpecialFunctions.gamma(1+n,a+b*x)+-1//2*d*b^-2*SpecialFunctions.gamma(2+n,a+b*x)+-1//2*b^-2*d^-1*(b*c+-1*a*d)^2*SpecialFunctions.gamma(n,a+b*x))
@test integrate(SpecialFunctions.gamma(n,a+b*x), x) == :(-1*b^-1*SpecialFunctions.gamma(1+n,a+b*x)+b^-1*(a+b*x)*SpecialFunctions.gamma(n,a+b*x))
@test integrate(x^2*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n))), x) == :(1//3*x^3*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n)))+-1//3*x^3*(c*x^n)^(-3*n^-1)*(d*(a+b*log(c*x^n)))^p*(-1*b^-1*n^-1*(3+-1*b*d*n)*(a+b*log(c*x^n)))^(-1p)*SpecialFunctions.gamma(p,-1*b^-1*n^-1*(3+-1*b*d*n)*(a+b*log(c*x^n)))*exp(-3*a*b^-1*n^-1))
@test integrate(x*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n))), x) == :((1/2)*x^2*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n)))+-1//2*x^2*(c*x^n)^(-2*n^-1)*(d*(a+b*log(c*x^n)))^p*(-1*b^-1*n^-1*(2+-1*b*d*n)*(a+b*log(c*x^n)))^(-1p)*SpecialFunctions.gamma(p,-1*b^-1*n^-1*(2+-1*b*d*n)*(a+b*log(c*x^n)))*exp(-2*a*b^-1*n^-1))
@test integrate(SpecialFunctions.gamma(p,d*(a+b*log(c*x^n))), x) == :(x*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n)))+-1*x*(c*x^n)^(-1*n^-1)*(d*(a+b*log(c*x^n)))^p*(-1*b^-1*n^-1*(1+-1*b*d*n)*(a+b*log(c*x^n)))^(-1p)*SpecialFunctions.gamma(p,-1*b^-1*n^-1*(1+-1*b*d*n)*(a+b*log(c*x^n)))*exp(-1*a*b^-1*n^-1))
@test integrate(x^-1*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n))), x) == :(b^-1*n^-1*(a+b*log(c*x^n))*SpecialFunctions.gamma(p,a*d+b*d*log(c*x^n))+-1*b^-1*d^-1*n^-1*SpecialFunctions.gamma(1+p,a*d+b*d*log(c*x^n)))
@test integrate(x^-2*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n))), x) == :(-1*x^-1*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n)))+x^-1*(c*x^n)^(n^-1)*(d*(a+b*log(c*x^n)))^p*(b^-1*n^-1*(1+b*d*n)*(a+b*log(c*x^n)))^(-1p)*SpecialFunctions.gamma(p,b^-1*n^-1*(1+b*d*n)*(a+b*log(c*x^n)))*exp(a*b^-1*n^-1))
@test integrate(x^-3*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n))), x) == :(-1//2*x^-2*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n)))+(1/2)*x^-2*(c*x^n)^(2*n^-1)*(d*(a+b*log(c*x^n)))^p*(b^-1*n^-1*(2+b*d*n)*(a+b*log(c*x^n)))^(-1p)*SpecialFunctions.gamma(p,b^-1*n^-1*(2+b*d*n)*(a+b*log(c*x^n)))*exp(2*a*b^-1*n^-1))
@test integrate((e*x)^m*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n))), x) == :(e^-1*(e*x)^(1+m)*(1+m)^-1*SpecialFunctions.gamma(p,d*(a+b*log(c*x^n)))+-1*e^-1*(c*x^n)^(-1*n^-1*(1+m))*(d*(a+b*log(c*x^n)))^p*(e*x)^(1+m)*(-1*b^-1*n^-1*(a+b*log(c*x^n))*(1+m+-1*b*d*n))^(-1p)*(1+m)^-1*SpecialFunctions.gamma(p,-1*b^-1*n^-1*(a+b*log(c*x^n))*(1+m+-1*b*d*n))*exp(-1*a*b^-1*n^-1*(1+m)))
@test integrate((c+d*x)^3*loggamma(a+b*x), x) == :(b^-1*(c+d*x)^3*SpecialFunctions.polygamma(-2,x)+-6*b^-4*d^3*SpecialFunctions.polygamma(-5,x)+-3*d*b^-2*(c+d*x)^2*SpecialFunctions.polygamma(-3,x)+6*b^-3*d^2*(c+d*x)*SpecialFunctions.polygamma(-4,x))
@test integrate((c+d*x)^2*loggamma(a+b*x), x) == :(b^-1*(c+d*x)^2*SpecialFunctions.polygamma(-2,x)+2*b^-3*d^2*SpecialFunctions.polygamma(-4,x)+-2*d*b^-2*(c+d*x)*SpecialFunctions.polygamma(-3,x))
@test integrate((c+d*x)*loggamma(a+b*x), x) == :(b^-1*(c+d*x)*SpecialFunctions.polygamma(-2,x)+-1*d*b^-2*SpecialFunctions.polygamma(-3,x))
@test integrate(loggamma(a+b*x), x) == :(b^-1*SpecialFunctions.polygamma(-2,x))
@test integrate((c+d*x)^2*log(SpecialFunctions.gamma(a+b*x)), x) == :(b^-1*(c+d*x)^2*SpecialFunctions.polygamma(-2,x)+2*b^-3*d^2*SpecialFunctions.polygamma(-4,x)+-1//3*d^-1*(c+d*x)^3*loggamma(a+b*x)+1//3*d^-1*(c+d*x)^3*log(SpecialFunctions.gamma(a+b*x))+-2*d*b^-2*(c+d*x)*SpecialFunctions.polygamma(-3,x))
@test integrate((c+d*x)*log(SpecialFunctions.gamma(a+b*x)), x) == :(b^-1*(c+d*x)*SpecialFunctions.polygamma(-2,x)+(1/2)*d^-1*(c+d*x)^2*log(SpecialFunctions.gamma(a+b*x))+-1*d*b^-2*SpecialFunctions.polygamma(-3,x)+-1//2*d^-1*(c+d*x)^2*loggamma(a+b*x))
@test integrate(log(SpecialFunctions.gamma(a+b*x)), x) == :(x*log(SpecialFunctions.gamma(a+b*x))+b^-1*SpecialFunctions.polygamma(-2,x)+-1*x*loggamma(a+b*x))
@test integrate((c+d*x)^3*SpecialFunctions.polygamma(n,x), x) == :(b^-1*(c+d*x)^3*SpecialFunctions.polygamma(-1+n,x)+-6*b^-4*d^3*SpecialFunctions.polygamma(-4+n,x)+-3*d*b^-2*(c+d*x)^2*SpecialFunctions.polygamma(-2+n,x)+6*b^-3*d^2*(c+d*x)*SpecialFunctions.polygamma(-3+n,x))
@test integrate((c+d*x)^2*SpecialFunctions.polygamma(n,x), x) == :(b^-1*(c+d*x)^2*SpecialFunctions.polygamma(-1+n,x)+2*b^-3*d^2*SpecialFunctions.polygamma(-3+n,x)+-2*d*b^-2*(c+d*x)*SpecialFunctions.polygamma(-2+n,x))
@test integrate((c+d*x)*SpecialFunctions.polygamma(n,x), x) == :(b^-1*(c+d*x)*SpecialFunctions.polygamma(-1+n,x)+-1*d*b^-2*SpecialFunctions.polygamma(-2+n,x))
@test integrate(SpecialFunctions.polygamma(n,x), x) == :(b^-1*SpecialFunctions.polygamma(-1+n,x))
@test integrate(x^2*SpecialFunctions.polygamma(1,x), x) == :(2*b^-3*SpecialFunctions.polygamma(-2,x)+b^-1*x^2*SpecialFunctions.polygamma(0,x)+-2*x*b^-2*loggamma(a+b*x))
@test integrate(x^-2*SpecialFunctions.polygamma(1,x)+-1*b*x^-1*SpecialFunctions.polygamma(2,x), x) == :(-1*x^-1*SpecialFunctions.polygamma(1,x))
@test integrate(x^-2*SpecialFunctions.polygamma(n,x)+-1*b*x^-1*SpecialFunctions.polygamma(1+n,x), x) == :(-1*x^-1*SpecialFunctions.polygamma(n,x))
@test integrate(SpecialFunctions.gamma(a+b*x)^n*SpecialFunctions.polygamma(0,x), x) == :(b^-1*n^-1*SpecialFunctions.gamma(a+b*x)^n)
@test integrate(factorial(a+b*x)^n*SpecialFunctions.polygamma(0,x), x) == :(b^-1*n^-1*factorial(a+b*x)^n)

max_stars_repo_path: test/8 Special functions/8.6 Gamma functions.jl


text:
@testset "environments" begin

    @testset "test API" begin

        for env in [CartPoleEnv(), WrappedEnv(CloneStatePreprocessor(), CartPoleEnv())]
            reset!(env)
            action_space = get_action_space(env)
            observation_space = get_observation_space(env)

            obs = observe(env)

            for _ in 1:1000
                if get_terminal(obs)
                    reset!(env)
                    obs = observe(env)
                end
                @test get_state(obs) ∈ observation_space
                action = rand(action_space)
                env(action)
                obs = observe(env)
            end
        end

    end

end