text:
@testset "4.4.2.1 (a+b cot)^m (c+d cot)^n" begin
    (A, B, a, b, c, d, e, f, m, n, x, ) = @variables A B a b c d e f m n x

    #= ::Package:: =#

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^n*(a+b*cot(e+f*x))^m=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^n*(a+im*a*cot(e+f*x))^m=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+im*a*cot(e+f*x))^m*(d*cot(e+f*x))^n*with*m*symbolic=#

    @test_int [(a + a*im*cot(c + d*x))^n, x, 2, (im*(a + im*a*cot(c + d*x))^n*HypergeometricFunctions._₂F₁(1, n, 1 + n, (1/2)*(1 + im*cot(c + d*x))))/(2*d*n)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^n*(a+a*cot(e+f*x))^m=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^(n/2)*(a+a*cot(e+f*x))^m=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [(a + a*cot(c + d*x))*(e*cot(c + d*x))^(5/2), x, 5, -((sqrt(2)*a*e^(5/2)*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/d) + (2*a*e^2*sqrt(e*cot(c + d*x)))/d - (2*a*e*(e*cot(c + d*x))^(3/2))/(3*d) - (2*a*(e*cot(c + d*x))^(5/2))/(5*d)]
    @test_int [(a + a*cot(c + d*x))*(e*cot(c + d*x))^(3/2), x, 4, -((sqrt(2)*a*e^(3/2)*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/d) - (2*a*e*sqrt(e*cot(c + d*x)))/d - (2*a*(e*cot(c + d*x))^(3/2))/(3*d)]
    @test_int [(a + a*cot(c + d*x))*(e*cot(c + d*x))^(1/2), x, 3, (sqrt(2)*a*sqrt(e)*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/d - (2*a*sqrt(e*cot(c + d*x)))/d]
    @test_int [(a + a*cot(c + d*x))/(e*cot(c + d*x))^(1/2), x, 2, (sqrt(2)*a*atan((sqrt(e)*(1 - cot(c + d*x)))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*sqrt(e))]
    @test_int [(a + a*cot(c + d*x))/(e*cot(c + d*x))^(3/2), x, 3, -((sqrt(2)*a*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*e^(3/2))) + (2*a)/(d*e*sqrt(e*cot(c + d*x)))]
    @test_int [(a + a*cot(c + d*x))/(e*cot(c + d*x))^(5/2), x, 4, -((sqrt(2)*a*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*e^(5/2))) + (2*a)/(3*d*e*(e*cot(c + d*x))^(3/2)) + (2*a)/(d*e^2*sqrt(e*cot(c + d*x)))]

    @test_int [(a + a*cot(c + d*x))^2*(e*cot(c + d*x))^(5/2), x, 16, (sqrt(2)*a^2*e^(5/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/d - (sqrt(2)*a^2*e^(5/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/d + (4*a^2*e^2*sqrt(e*cot(c + d*x)))/d - (4*a^2*(e*cot(c + d*x))^(5/2))/(5*d) - (2*a^2*(e*cot(c + d*x))^(7/2))/(7*d*e) + (a^2*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d) - (a^2*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d)]
    @test_int [(a + a*cot(c + d*x))^2*(e*cot(c + d*x))^(3/2), x, 15, -((sqrt(2)*a^2*e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/d) + (sqrt(2)*a^2*e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/d - (4*a^2*(e*cot(c + d*x))^(3/2))/(3*d) - (2*a^2*(e*cot(c + d*x))^(5/2))/(5*d*e) + (a^2*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d) - (a^2*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d)]
    @test_int [(a + a*cot(c + d*x))^2*(e*cot(c + d*x))^(1/2), x, 15, -((sqrt(2)*a^2*sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/d) + (sqrt(2)*a^2*sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/d - (4*a^2*sqrt(e*cot(c + d*x)))/d - (2*a^2*(e*cot(c + d*x))^(3/2))/(3*d*e) - (a^2*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d) + (a^2*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d)]
    @test_int [(a + a*cot(c + d*x))^2/(e*cot(c + d*x))^(1/2), x, 14, (sqrt(2)*a^2*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*sqrt(e)) - (sqrt(2)*a^2*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*sqrt(e)) - (2*a^2*sqrt(e*cot(c + d*x)))/(d*e) - (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*sqrt(e)) + (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*sqrt(e))]
    @test_int [(a + a*cot(c + d*x))^2/(e*cot(c + d*x))^(3/2), x, 13, (sqrt(2)*a^2*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*e^(3/2)) - (sqrt(2)*a^2*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*e^(3/2)) + (2*a^2)/(d*e*sqrt(e*cot(c + d*x))) + (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*e^(3/2)) - (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*e^(3/2))]
    @test_int [(a + a*cot(c + d*x))^2/(e*cot(c + d*x))^(5/2), x, 14, -((sqrt(2)*a^2*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*e^(5/2))) + (sqrt(2)*a^2*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*e^(5/2)) + (2*a^2)/(3*d*e*(e*cot(c + d*x))^(3/2)) + (4*a^2)/(d*e^2*sqrt(e*cot(c + d*x))) + (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*e^(5/2)) - (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*e^(5/2))]
    @test_int [(a + a*cot(c + d*x))^2/(e*cot(c + d*x))^(7/2), x, 14, -((sqrt(2)*a^2*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*e^(7/2))) + (sqrt(2)*a^2*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(d*e^(7/2)) + (2*a^2)/(5*d*e*(e*cot(c + d*x))^(5/2)) + (4*a^2)/(3*d*e^2*(e*cot(c + d*x))^(3/2)) - (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*e^(7/2)) + (a^2*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*d*e^(7/2))]

    @test_int [(a + a*cot(c + d*x))^3*(e*cot(c + d*x))^(5/2), x, 7, (2*sqrt(2)*a^3*e^(5/2)*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/d + (4*a^3*e^2*sqrt(e*cot(c + d*x)))/d + (4*a^3*e*(e*cot(c + d*x))^(3/2))/(3*d) - (4*a^3*(e*cot(c + d*x))^(5/2))/(5*d) - (40*a^3*(e*cot(c + d*x))^(7/2))/(63*d*e) - (2*(e*cot(c + d*x))^(7/2)*(a^3 + a^3*cot(c + d*x)))/(9*d*e)]
    @test_int [(a + a*cot(c + d*x))^3*(e*cot(c + d*x))^(3/2), x, 6, -((2*sqrt(2)*a^3*e^(3/2)*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/d) + (4*a^3*e*sqrt(e*cot(c + d*x)))/d - (4*a^3*(e*cot(c + d*x))^(3/2))/(3*d) - (32*a^3*(e*cot(c + d*x))^(5/2))/(35*d*e) - (2*(e*cot(c + d*x))^(5/2)*(a^3 + a^3*cot(c + d*x)))/(7*d*e)]
    @test_int [(a + a*cot(c + d*x))^3*(e*cot(c + d*x))^(1/2), x, 5, -((2*sqrt(2)*a^3*sqrt(e)*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/d) - (4*a^3*sqrt(e*cot(c + d*x)))/d - (8*a^3*(e*cot(c + d*x))^(3/2))/(5*d*e) - (2*(e*cot(c + d*x))^(3/2)*(a^3 + a^3*cot(c + d*x)))/(5*d*e)]
    @test_int [(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(1/2), x, 4, (2*sqrt(2)*a^3*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*sqrt(e)) - (16*a^3*sqrt(e*cot(c + d*x)))/(3*d*e) - (2*sqrt(e*cot(c + d*x))*(a^3 + a^3*cot(c + d*x)))/(3*d*e)]
    @test_int [(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(3/2), x, 4, (2*sqrt(2)*a^3*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*e^(3/2)) - (4*a^3*sqrt(e*cot(c + d*x)))/(d*e^2) + (2*(a^3 + a^3*cot(c + d*x)))/(d*e*sqrt(e*cot(c + d*x)))]
    @test_int [(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(5/2), x, 4, -((2*sqrt(2)*a^3*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*e^(5/2))) + (16*a^3)/(3*d*e^2*sqrt(e*cot(c + d*x))) + (2*(a^3 + a^3*cot(c + d*x)))/(3*d*e*(e*cot(c + d*x))^(3/2))]
    @test_int [(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(7/2), x, 5, -((2*sqrt(2)*a^3*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*e^(7/2))) + (8*a^3)/(5*d*e^2*(e*cot(c + d*x))^(3/2)) + (4*a^3)/(d*e^3*sqrt(e*cot(c + d*x))) + (2*(a^3 + a^3*cot(c + d*x)))/(5*d*e*(e*cot(c + d*x))^(5/2))]
    @test_int [(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(9/2), x, 6, (2*sqrt(2)*a^3*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(d*e^(9/2)) + (32*a^3)/(35*d*e^2*(e*cot(c + d*x))^(5/2)) + (4*a^3)/(3*d*e^3*(e*cot(c + d*x))^(3/2)) - (4*a^3)/(d*e^4*sqrt(e*cot(c + d*x))) + (2*(a^3 + a^3*cot(c + d*x)))/(7*d*e*(e*cot(c + d*x))^(7/2))]

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [1/(a + a*cot(c + d*x))*(e*cot(c + d*x))^(5/2), x, 7, (e^(5/2)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(a*d) - (e^(5/2)*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(sqrt(2)*a*d) - (2*e^2*sqrt(e*cot(c + d*x)))/(a*d)]
    @test_int [1/(a + a*cot(c + d*x))*(e*cot(c + d*x))^(3/2), x, 6, -((e^(3/2)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(a*d)) + (e^(3/2)*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(sqrt(2)*a*d)]
    @test_int [1/(a + a*cot(c + d*x))*(e*cot(c + d*x))^(1/2), x, 6, (sqrt(e)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(a*d) + (sqrt(e)*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(sqrt(2)*a*d)]
    @test_int [1/(a + a*cot(c + d*x))/(e*cot(c + d*x))^(1/2), x, 6, -(atan(sqrt(e*cot(c + d*x))/sqrt(e))/(a*d*sqrt(e))) - atanh((sqrt(e)*(1 + cot(c + d*x)))/(sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*a*d*sqrt(e))]
    @test_int [1/(a + a*cot(c + d*x))/(e*cot(c + d*x))^(3/2), x, 7, atan(sqrt(e*cot(c + d*x))/sqrt(e))/(a*d*e^(3/2)) - atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*a*d*e^(3/2)) + 2/(a*d*e*sqrt(e*cot(c + d*x)))]
    @test_int [1/(a + a*cot(c + d*x))/(e*cot(c + d*x))^(5/2), x, 10, -(atan(sqrt(e*cot(c + d*x))/sqrt(e))/(a*d*e^(5/2))) + atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x))))/(sqrt(2)*a*d*e^(5/2)) + 2/(3*a*d*e*(e*cot(c + d*x))^(3/2)) - 2/(a*d*e^2*sqrt(e*cot(c + d*x)))]

    @test_int [1/(a + a*cot(c + d*x))^2*(e*cot(c + d*x))^(5/2), x, 17, -((3*e^(5/2)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(2*a^2*d)) - (e^(5/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(2*sqrt(2)*a^2*d) + (e^(5/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(2*sqrt(2)*a^2*d) + (e^2*sqrt(e*cot(c + d*x)))/(2*d*(a^2 + a^2*cot(c + d*x))) - (e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(4*sqrt(2)*a^2*d) + (e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(4*sqrt(2)*a^2*d)]
    @test_int [1/(a + a*cot(c + d*x))^2*(e*cot(c + d*x))^(3/2), x, 18, (e^(3/2)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(2*a^2*d) + (e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(2*sqrt(2)*a^2*d) - (e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(2*sqrt(2)*a^2*d) - (e*sqrt(e*cot(c + d*x)))/(2*d*(a^2 + a^2*cot(c + d*x))) - (e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(4*sqrt(2)*a^2*d) + (e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(4*sqrt(2)*a^2*d)]
    @test_int [1/(a + a*cot(c + d*x))^2*(e*cot(c + d*x))^(1/2), x, 17, (sqrt(e)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(2*a^2*d) + (sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(2*sqrt(2)*a^2*d) - (sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(2*sqrt(2)*a^2*d) + sqrt(e*cot(c + d*x))/(2*d*(a^2 + a^2*cot(c + d*x))) + (sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(4*sqrt(2)*a^2*d) - (sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(4*sqrt(2)*a^2*d)]
    @test_int [1/(a + a*cot(c + d*x))^2/(e*cot(c + d*x))^(1/2), x, 18, -((3*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(2*a^2*d*sqrt(e))) - atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e))/(2*sqrt(2)*a^2*d*sqrt(e)) + atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e))/(2*sqrt(2)*a^2*d*sqrt(e)) - sqrt(e*cot(c + d*x))/(2*d*e*(a^2 + a^2*cot(c + d*x))) + log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x)))/(4*sqrt(2)*a^2*d*sqrt(e)) - log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x)))/(4*sqrt(2)*a^2*d*sqrt(e))]
    @test_int [1/(a + a*cot(c + d*x))^2/(e*cot(c + d*x))^(3/2), x, 18, (5*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(2*a^2*d*e^(3/2)) - atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e))/(2*sqrt(2)*a^2*d*e^(3/2)) + atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e))/(2*sqrt(2)*a^2*d*e^(3/2)) + 5/(2*a^2*d*e*sqrt(e*cot(c + d*x))) - 1/(2*d*e*sqrt(e*cot(c + d*x))*(a^2 + a^2*cot(c + d*x))) - log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x)))/(4*sqrt(2)*a^2*d*e^(3/2)) + log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x)))/(4*sqrt(2)*a^2*d*e^(3/2))]
    @test_int [1/(a + a*cot(c + d*x))^2/(e*cot(c + d*x))^(5/2), x, 20, -((7*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(2*a^2*d*e^(5/2))) + atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e))/(2*sqrt(2)*a^2*d*e^(5/2)) - atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e))/(2*sqrt(2)*a^2*d*e^(5/2)) + 7/(6*a^2*d*e*(e*cot(c + d*x))^(3/2)) - 9/(2*a^2*d*e^2*sqrt(e*cot(c + d*x))) - 1/(2*d*e*(e*cot(c + d*x))^(3/2)*(a^2 + a^2*cot(c + d*x))) - log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x)))/(4*sqrt(2)*a^2*d*e^(5/2)) + log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x)))/(4*sqrt(2)*a^2*d*e^(5/2))]

    @test_int [1/(a + a*cot(c + d*x))^3*(e*cot(c + d*x))^(5/2), x, 8, -((e^(5/2)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(8*a^3*d)) + (e^(5/2)*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(2*sqrt(2)*a^3*d) - (5*e^2*sqrt(e*cot(c + d*x)))/(8*a^3*d*(1 + cot(c + d*x))) + (e^2*sqrt(e*cot(c + d*x)))/(4*a*d*(a + a*cot(c + d*x))^2)]
    @test_int [1/(a + a*cot(c + d*x))^3*(e*cot(c + d*x))^(3/2), x, 8, (5*e^(3/2)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(8*a^3*d) + (e^(3/2)*atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(2*sqrt(2)*a^3*d) - (e*sqrt(e*cot(c + d*x)))/(4*a*d*(a + a*cot(c + d*x))^2) + (e*sqrt(e*cot(c + d*x)))/(8*d*(a^3 + a^3*cot(c + d*x)))]
    @test_int [1/(a + a*cot(c + d*x))^3*(e*cot(c + d*x))^(1/2), x, 8, -((sqrt(e)*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(8*a^3*d)) - (sqrt(e)*atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x)))))/(2*sqrt(2)*a^3*d) + sqrt(e*cot(c + d*x))/(4*a*d*(a + a*cot(c + d*x))^2) + (3*sqrt(e*cot(c + d*x)))/(8*d*(a^3 + a^3*cot(c + d*x)))]
    @test_int [1/(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(1/2), x, 8, -((11*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(8*a^3*d*sqrt(e))) - atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*a^3*d*sqrt(e)) - (7*sqrt(e*cot(c + d*x)))/(8*a^3*d*e*(1 + cot(c + d*x))) - sqrt(e*cot(c + d*x))/(4*a*d*e*(a + a*cot(c + d*x))^2)]
    @test_int [1/(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(3/2), x, 9, (31*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(8*a^3*d*e^(3/2)) + atanh((sqrt(e) + sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*a^3*d*e^(3/2)) + 27/(8*a^3*d*e*sqrt(e*cot(c + d*x))) - 9/(8*a^3*d*e*sqrt(e*cot(c + d*x))*(1 + cot(c + d*x))) - 1/(4*a*d*e*sqrt(e*cot(c + d*x))*(a + a*cot(c + d*x))^2)]
    @test_int [1/(a + a*cot(c + d*x))^3/(e*cot(c + d*x))^(5/2), x, 10, -((59*atan(sqrt(e*cot(c + d*x))/sqrt(e)))/(8*a^3*d*e^(5/2))) + atan((sqrt(e) - sqrt(e)*cot(c + d*x))/(sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*a^3*d*e^(5/2)) + 55/(24*a^3*d*e*(e*cot(c + d*x))^(3/2)) - 63/(8*a^3*d*e^2*sqrt(e*cot(c + d*x))) - 11/(8*a^3*d*e*(e*cot(c + d*x))^(3/2)*(1 + cot(c + d*x))) - 1/(4*a*d*e*(e*cot(c + d*x))^(3/2)*(a + a*cot(c + d*x))^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^n*(a+a*cot(e+f*x))^(m/2)=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [cot(x)^2*sqrt(1 + cot(x)), x, 12, (-sqrt((1/2)*(1 + sqrt(2))))*atan((sqrt(2*(1 + sqrt(2))) - 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) + sqrt((1/2)*(1 + sqrt(2)))*atan((sqrt(2*(1 + sqrt(2))) + 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) - (2/3)*(1 + cot(x))^(3/2) + log(1 + sqrt(2) + cot(x) - sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(2*sqrt(2*(1 + sqrt(2)))) - log(1 + sqrt(2) + cot(x) + sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(2*sqrt(2*(1 + sqrt(2))))]
    @test_int [cot(x)^1*sqrt(1 + cot(x)), x, 6, sqrt((1/2)*(-1 + sqrt(2)))*atan((4 - 3*sqrt(2) + (2 - sqrt(2))*cot(x))/(2*sqrt(-7 + 5*sqrt(2))*sqrt(1 + cot(x)))) + sqrt((1/2)*(1 + sqrt(2)))*atanh((4 + 3*sqrt(2) + (2 + sqrt(2))*cot(x))/(2*sqrt(7 + 5*sqrt(2))*sqrt(1 + cot(x)))) - 2*sqrt(1 + cot(x))]

    @test_int [cot(x)^2*(1 + cot(x))^(3/2), x, 8, (-sqrt(-1 + sqrt(2)))*atan((3 - 2*sqrt(2) + (1 - sqrt(2))*cot(x))/(sqrt(2*(-7 + 5*sqrt(2)))*sqrt(1 + cot(x)))) - sqrt(1 + sqrt(2))*atanh((3 + 2*sqrt(2) + (1 + sqrt(2))*cot(x))/(sqrt(2*(7 + 5*sqrt(2)))*sqrt(1 + cot(x)))) + 2*sqrt(1 + cot(x)) - (2/5)*(1 + cot(x))^(5/2)]
    @test_int [cot(x)^1*(1 + cot(x))^(3/2), x, 14, (-sqrt(1 + sqrt(2)))*atan((sqrt(2*(1 + sqrt(2))) - 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) + sqrt(1 + sqrt(2))*atan((sqrt(2*(1 + sqrt(2))) + 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) - 2*sqrt(1 + cot(x)) - (2/3)*(1 + cot(x))^(3/2) - log(1 + sqrt(2) + cot(x) - sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(2*sqrt(1 + sqrt(2))) + log(1 + sqrt(2) + cot(x) + sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(2*sqrt(1 + sqrt(2)))]

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [cot(x)^2/sqrt(1 + cot(x)), x, 12, (-(1/2))*sqrt(1 + sqrt(2))*atan((sqrt(2*(1 + sqrt(2))) - 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) + (1/2)*sqrt(1 + sqrt(2))*atan((sqrt(2*(1 + sqrt(2))) + 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) - 2*sqrt(1 + cot(x)) - log(1 + sqrt(2) + cot(x) - sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(4*sqrt(1 + sqrt(2))) + log(1 + sqrt(2) + cot(x) + sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(4*sqrt(1 + sqrt(2)))]
    @test_int [cot(x)^1/sqrt(1 + cot(x)), x, 5, (1/2)*sqrt(-1 + sqrt(2))*atan((3 - 2*sqrt(2) + (1 - sqrt(2))*cot(x))/(sqrt(2*(-7 + 5*sqrt(2)))*sqrt(1 + cot(x)))) + (1/2)*sqrt(1 + sqrt(2))*atanh((3 + 2*sqrt(2) + (1 + sqrt(2))*cot(x))/(sqrt(2*(7 + 5*sqrt(2)))*sqrt(1 + cot(x))))]

    @test_int [cot(x)^2/(1 + cot(x))^(3/2), x, 6, (1/2)*sqrt((1/2)*(-1 + sqrt(2)))*atan((4 - 3*sqrt(2) + (2 - sqrt(2))*cot(x))/(2*sqrt(-7 + 5*sqrt(2))*sqrt(1 + cot(x)))) + (1/2)*sqrt((1/2)*(1 + sqrt(2)))*atanh((4 + 3*sqrt(2) + (2 + sqrt(2))*cot(x))/(2*sqrt(7 + 5*sqrt(2))*sqrt(1 + cot(x)))) + 1/sqrt(1 + cot(x))]
    @test_int [cot(x)^1/(1 + cot(x))^(3/2), x, 13, (1/2)*sqrt((1/2)*(1 + sqrt(2)))*atan((sqrt(2*(1 + sqrt(2))) - 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) - (1/2)*sqrt((1/2)*(1 + sqrt(2)))*atan((sqrt(2*(1 + sqrt(2))) + 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) - 1/sqrt(1 + cot(x)) - log(1 + sqrt(2) + cot(x) - sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(4*sqrt(2*(1 + sqrt(2)))) + log(1 + sqrt(2) + cot(x) + sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(4*sqrt(2*(1 + sqrt(2))))]

    @test_int [cot(x)^2/(1 + cot(x))^(5/2), x, 8, (1/4)*sqrt(-1 + sqrt(2))*atan((3 - 2*sqrt(2) + (1 - sqrt(2))*cot(x))/(sqrt(2*(-7 + 5*sqrt(2)))*sqrt(1 + cot(x)))) + (1/4)*sqrt(1 + sqrt(2))*atanh((3 + 2*sqrt(2) + (1 + sqrt(2))*cot(x))/(sqrt(2*(7 + 5*sqrt(2)))*sqrt(1 + cot(x)))) + 1/(3*(1 + cot(x))^(3/2)) - 1/sqrt(1 + cot(x))]
    @test_int [cot(x)^1/(1 + cot(x))^(5/2), x, 13, (1/4)*sqrt(1 + sqrt(2))*atan((sqrt(2*(1 + sqrt(2))) - 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) - (1/4)*sqrt(1 + sqrt(2))*atan((sqrt(2*(1 + sqrt(2))) + 2*sqrt(1 + cot(x)))/sqrt(2*(-1 + sqrt(2)))) - 1/(3*(1 + cot(x))^(3/2)) + log(1 + sqrt(2) + cot(x) - sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(8*sqrt(1 + sqrt(2))) - log(1 + sqrt(2) + cot(x) + sqrt(2*(1 + sqrt(2)))*sqrt(1 + cot(x)))/(8*sqrt(1 + sqrt(2)))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^n*(a+b*cot(e+f*x))^m=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^(n/2)*(a+b*cot(e+f*x))^m=#

    #= ::Subsubsection::Closed:: =#
    #=m>0=#

    @test_int [(a + b*cot(c + d*x))*(e*cot(c + d*x))^(3/2), x, 12, -(((a + b)*e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d)) + ((a + b)*e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - (2*a*e*sqrt(e*cot(c + d*x)))/d - (2*b*(e*cot(c + d*x))^(3/2))/(3*d) - ((a - b)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d) + ((a - b)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d)]
    @test_int [(a + b*cot(c + d*x))*(e*cot(c + d*x))^(1/2), x, 11, ((a - b)*sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - ((a - b)*sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - (2*b*sqrt(e*cot(c + d*x)))/d - ((a + b)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d) + ((a + b)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d)]
    @test_int [(a + b*cot(c + d*x))/(e*cot(c + d*x))^(1/2), x, 10, ((a + b)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*sqrt(e)) - ((a + b)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*sqrt(e)) + ((a - b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*sqrt(e)) - ((a - b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*sqrt(e))]
    @test_int [(a + b*cot(c + d*x))/(e*cot(c + d*x))^(3/2), x, 11, -(((a - b)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(3/2))) + ((a - b)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(3/2)) + (2*a)/(d*e*sqrt(e*cot(c + d*x))) + ((a + b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(3/2)) - ((a + b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(3/2))]
    @test_int [(a + b*cot(c + d*x))/(e*cot(c + d*x))^(5/2), x, 12, -(((a + b)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(5/2))) + ((a + b)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(5/2)) + (2*a)/(3*d*e*(e*cot(c + d*x))^(3/2)) + (2*b)/(d*e^2*sqrt(e*cot(c + d*x))) - ((a - b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(5/2)) + ((a - b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(5/2))]

    @test_int [(a + b*cot(c + d*x))^2*(e*cot(c + d*x))^(3/2), x, 13, -(((a^2 + 2*a*b - b^2)*e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d)) + ((a^2 + 2*a*b - b^2)*e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - (2*(a^2 - b^2)*e*sqrt(e*cot(c + d*x)))/d - (4*a*b*(e*cot(c + d*x))^(3/2))/(3*d) - (2*b^2*(e*cot(c + d*x))^(5/2))/(5*d*e) - ((a^2 - 2*a*b - b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d) + ((a^2 - 2*a*b - b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d)]
    @test_int [(a + b*cot(c + d*x))^2*(e*cot(c + d*x))^(1/2), x, 12, ((a^2 - 2*a*b - b^2)*sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - ((a^2 - 2*a*b - b^2)*sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - (4*a*b*sqrt(e*cot(c + d*x)))/d - (2*b^2*(e*cot(c + d*x))^(3/2))/(3*d*e) - ((a^2 + 2*a*b - b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d) + ((a^2 + 2*a*b - b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d)]
    @test_int [(a + b*cot(c + d*x))^2/(e*cot(c + d*x))^(1/2), x, 11, ((a^2 + 2*a*b - b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*sqrt(e)) - ((a^2 + 2*a*b - b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*sqrt(e)) - (2*b^2*sqrt(e*cot(c + d*x)))/(d*e) + ((a^2 - 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*sqrt(e)) - ((a^2 - 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*sqrt(e))]
    @test_int [(a + b*cot(c + d*x))^2/(e*cot(c + d*x))^(3/2), x, 11, -(((a^2 - 2*a*b - b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(3/2))) + ((a^2 - 2*a*b - b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(3/2)) + (2*a^2)/(d*e*sqrt(e*cot(c + d*x))) + ((a^2 + 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(3/2)) - ((a^2 + 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(3/2))]
    @test_int [(a + b*cot(c + d*x))^2/(e*cot(c + d*x))^(5/2), x, 12, -(((a^2 + 2*a*b - b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(5/2))) + ((a^2 + 2*a*b - b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(5/2)) + (2*a^2)/(3*d*e*(e*cot(c + d*x))^(3/2)) + (4*a*b)/(d*e^2*sqrt(e*cot(c + d*x))) - ((a^2 - 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(5/2)) + ((a^2 - 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(5/2))]
    @test_int [(a + b*cot(c + d*x))^2/(e*cot(c + d*x))^(7/2), x, 13, ((a^2 - 2*a*b - b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(7/2)) - ((a^2 - 2*a*b - b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(7/2)) + (2*a^2)/(5*d*e*(e*cot(c + d*x))^(5/2)) + (4*a*b)/(3*d*e^2*(e*cot(c + d*x))^(3/2)) - (2*(a^2 - b^2))/(d*e^3*sqrt(e*cot(c + d*x))) - ((a^2 + 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(7/2)) + ((a^2 + 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(7/2))]

    @test_int [(a + b*cot(c + d*x))^3*(e*cot(c + d*x))^(3/2), x, 14, -(((a - b)*(a^2 + 4*a*b + b^2)*e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d)) + ((a - b)*(a^2 + 4*a*b + b^2)*e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - (2*a*(a^2 - 3*b^2)*e*sqrt(e*cot(c + d*x)))/d - (2*b*(3*a^2 - b^2)*(e*cot(c + d*x))^(3/2))/(3*d) - (32*a*b^2*(e*cot(c + d*x))^(5/2))/(35*d*e) - (2*b^2*(e*cot(c + d*x))^(5/2)*(a + b*cot(c + d*x)))/(7*d*e) - ((a + b)*(a^2 - 4*a*b + b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d) + ((a + b)*(a^2 - 4*a*b + b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d)]
    @test_int [(a + b*cot(c + d*x))^3*(e*cot(c + d*x))^(1/2), x, 13, ((a + b)*(a^2 - 4*a*b + b^2)*sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - ((a + b)*(a^2 - 4*a*b + b^2)*sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d) - (2*b*(3*a^2 - b^2)*sqrt(e*cot(c + d*x)))/d - (8*a*b^2*(e*cot(c + d*x))^(3/2))/(5*d*e) - (2*b^2*(e*cot(c + d*x))^(3/2)*(a + b*cot(c + d*x)))/(5*d*e) - ((a - b)*(a^2 + 4*a*b + b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d) + ((a - b)*(a^2 + 4*a*b + b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d)]
    @test_int [(a + b*cot(c + d*x))^3/(e*cot(c + d*x))^(1/2), x, 12, ((a - b)*(a^2 + 4*a*b + b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*sqrt(e)) - ((a - b)*(a^2 + 4*a*b + b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*sqrt(e)) - (16*a*b^2*sqrt(e*cot(c + d*x)))/(3*d*e) - (2*b^2*sqrt(e*cot(c + d*x))*(a + b*cot(c + d*x)))/(3*d*e) + ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*sqrt(e)) - ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*sqrt(e))]
    @test_int [(a + b*cot(c + d*x))^3/(e*cot(c + d*x))^(3/2), x, 12, -(((a + b)*(a^2 - 4*a*b + b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(3/2))) + ((a + b)*(a^2 - 4*a*b + b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(3/2)) - (2*b*(a^2 + b^2)*sqrt(e*cot(c + d*x)))/(d*e^2) + (2*a^2*(a + b*cot(c + d*x)))/(d*e*sqrt(e*cot(c + d*x))) + ((a - b)*(a^2 + 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(3/2)) - ((a - b)*(a^2 + 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(3/2))]
    @test_int [(a + b*cot(c + d*x))^3/(e*cot(c + d*x))^(5/2), x, 12, -(((a - b)*(a^2 + 4*a*b + b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(5/2))) + ((a - b)*(a^2 + 4*a*b + b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(5/2)) + (16*a^2*b)/(3*d*e^2*sqrt(e*cot(c + d*x))) + (2*a^2*(a + b*cot(c + d*x)))/(3*d*e*(e*cot(c + d*x))^(3/2)) - ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(5/2)) + ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(5/2))]
    @test_int [(a + b*cot(c + d*x))^3/(e*cot(c + d*x))^(7/2), x, 13, ((a + b)*(a^2 - 4*a*b + b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(7/2)) - ((a + b)*(a^2 - 4*a*b + b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(7/2)) + (8*a^2*b)/(5*d*e^2*(e*cot(c + d*x))^(3/2)) - (2*a*(a^2 - 3*b^2))/(d*e^3*sqrt(e*cot(c + d*x))) + (2*a^2*(a + b*cot(c + d*x)))/(5*d*e*(e*cot(c + d*x))^(5/2)) - ((a - b)*(a^2 + 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(7/2)) + ((a - b)*(a^2 + 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(7/2))]
    @test_int [(a + b*cot(c + d*x))^3/(e*cot(c + d*x))^(9/2), x, 14, ((a - b)*(a^2 + 4*a*b + b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(9/2)) - ((a - b)*(a^2 + 4*a*b + b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*d*e^(9/2)) + (32*a^2*b)/(35*d*e^2*(e*cot(c + d*x))^(5/2)) - (2*a*(a^2 - 3*b^2))/(3*d*e^3*(e*cot(c + d*x))^(3/2)) - (2*b*(3*a^2 - b^2))/(d*e^4*sqrt(e*cot(c + d*x))) + (2*a^2*(a + b*cot(c + d*x)))/(7*d*e*(e*cot(c + d*x))^(7/2)) + ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(9/2)) - ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*d*e^(9/2))]

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [1/(a + b*cot(c + d*x))*(e*cot(c + d*x))^(5/2), x, 15, (2*a^(5/2)*e^(5/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(b^(3/2)*(a^2 + b^2)*d) - ((a + b)*e^(5/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d) + ((a + b)*e^(5/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d) - (2*e^2*sqrt(e*cot(c + d*x)))/(b*d) + ((a - b)*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d) - ((a - b)*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d)]
    @test_int [1/(a + b*cot(c + d*x))*(e*cot(c + d*x))^(3/2), x, 14, -((2*a^(3/2)*e^(3/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(sqrt(b)*(a^2 + b^2)*d)) - ((a - b)*e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d) + ((a - b)*e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d) - ((a + b)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d) + ((a + b)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d)]
    @test_int [1/(a + b*cot(c + d*x))*(e*cot(c + d*x))^(1/2), x, 14, (2*sqrt(a)*sqrt(b)*sqrt(e)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/((a^2 + b^2)*d) + ((a + b)*sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d) - ((a + b)*sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d) - ((a - b)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d) + ((a - b)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d)]
    @test_int [1/(a + b*cot(c + d*x))/(e*cot(c + d*x))^(1/2), x, 14, -((2*b^(3/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(sqrt(a)*(a^2 + b^2)*d*sqrt(e))) + ((a - b)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d*sqrt(e)) - ((a - b)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d*sqrt(e)) + ((a + b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d*sqrt(e)) - ((a + b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d*sqrt(e))]
    @test_int [1/(a + b*cot(c + d*x))/(e*cot(c + d*x))^(3/2), x, 15, (2*b^(5/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(a^(3/2)*(a^2 + b^2)*d*e^(3/2)) - ((a + b)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d*e^(3/2)) + ((a + b)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d*e^(3/2)) + 2/(a*d*e*sqrt(e*cot(c + d*x))) + ((a - b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d*e^(3/2)) - ((a - b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d*e^(3/2))]
    @test_int [1/(a + b*cot(c + d*x))/(e*cot(c + d*x))^(5/2), x, 16, -((2*b^(7/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(a^(5/2)*(a^2 + b^2)*d*e^(5/2))) - ((a - b)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d*e^(5/2)) + ((a - b)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)*d*e^(5/2)) + 2/(3*a*d*e*(e*cot(c + d*x))^(3/2)) - (2*b)/(a^2*d*e^2*sqrt(e*cot(c + d*x))) - ((a + b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d*e^(5/2)) + ((a + b)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)*d*e^(5/2))]

    @test_int [1/(a + b*cot(c + d*x))^2*(e*cot(c + d*x))^(7/2), x, 16, (a^(5/2)*(3*a^2 + 7*b^2)*e^(7/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(b^(5/2)*(a^2 + b^2)^2*d) + ((a^2 - 2*a*b - b^2)*e^(7/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) - ((a^2 - 2*a*b - b^2)*e^(7/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) - ((3*a^2 + 2*b^2)*e^3*sqrt(e*cot(c + d*x)))/(b^2*(a^2 + b^2)*d) + (a^2*e^2*(e*cot(c + d*x))^(3/2))/(b*(a^2 + b^2)*d*(a + b*cot(c + d*x))) + ((a^2 + 2*a*b - b^2)*e^(7/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d) - ((a^2 + 2*a*b - b^2)*e^(7/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d)]
    @test_int [1/(a + b*cot(c + d*x))^2*(e*cot(c + d*x))^(5/2), x, 15, -((a^(3/2)*(a^2 + 5*b^2)*e^(5/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(b^(3/2)*(a^2 + b^2)^2*d)) - ((a^2 + 2*a*b - b^2)*e^(5/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) + ((a^2 + 2*a*b - b^2)*e^(5/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) + (a^2*e^2*sqrt(e*cot(c + d*x)))/(b*(a^2 + b^2)*d*(a + b*cot(c + d*x))) + ((a^2 - 2*a*b - b^2)*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d) - ((a^2 - 2*a*b - b^2)*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d)]
    @test_int [1/(a + b*cot(c + d*x))^2*(e*cot(c + d*x))^(3/2), x, 15, -((sqrt(a)*(a^2 - 3*b^2)*e^(3/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(sqrt(b)*(a^2 + b^2)^2*d)) - ((a^2 - 2*a*b - b^2)*e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) + ((a^2 - 2*a*b - b^2)*e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) - (a*e*sqrt(e*cot(c + d*x)))/((a^2 + b^2)*d*(a + b*cot(c + d*x))) - ((a^2 + 2*a*b - b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d) + ((a^2 + 2*a*b - b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d)]
    @test_int [1/(a + b*cot(c + d*x))^2*(e*cot(c + d*x))^(1/2), x, 15, (sqrt(b)*(3*a^2 - b^2)*sqrt(e)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(sqrt(a)*(a^2 + b^2)^2*d) + ((a^2 + 2*a*b - b^2)*sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) - ((a^2 + 2*a*b - b^2)*sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d) + (b*sqrt(e*cot(c + d*x)))/((a^2 + b^2)*d*(a + b*cot(c + d*x))) - ((a^2 - 2*a*b - b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d) + ((a^2 - 2*a*b - b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d)]
    @test_int [1/(a + b*cot(c + d*x))^2/(e*cot(c + d*x))^(1/2), x, 15, -((b^(3/2)*(5*a^2 + b^2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(a^(3/2)*(a^2 + b^2)^2*d*sqrt(e))) + ((a^2 - 2*a*b - b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d*sqrt(e)) - ((a^2 - 2*a*b - b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d*sqrt(e)) - (b^2*sqrt(e*cot(c + d*x)))/(a*(a^2 + b^2)*d*e*(a + b*cot(c + d*x))) + ((a^2 + 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d*sqrt(e)) - ((a^2 + 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d*sqrt(e))]
    @test_int [1/(a + b*cot(c + d*x))^2/(e*cot(c + d*x))^(3/2), x, 16, (b^(5/2)*(7*a^2 + 3*b^2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(a^(5/2)*(a^2 + b^2)^2*d*e^(3/2)) - ((a^2 + 2*a*b - b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d*e^(3/2)) + ((a^2 + 2*a*b - b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^2*d*e^(3/2)) + (2*a^2 + 3*b^2)/(a^2*(a^2 + b^2)*d*e*sqrt(e*cot(c + d*x))) - b^2/(a*(a^2 + b^2)*d*e*sqrt(e*cot(c + d*x))*(a + b*cot(c + d*x))) + ((a^2 - 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d*e^(3/2)) - ((a^2 - 2*a*b - b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^2*d*e^(3/2))]

    @test_int [1/(a + b*cot(c + d*x))^3*(e*cot(c + d*x))^(9/2), x, 17, (a^(5/2)*(15*a^4 + 46*a^2*b^2 + 63*b^4)*e^(9/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(4*b^(7/2)*(a^2 + b^2)^3*d) + ((a - b)*(a^2 + 4*a*b + b^2)*e^(9/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) - ((a - b)*(a^2 + 4*a*b + b^2)*e^(9/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) - ((15*a^4 + 31*a^2*b^2 + 8*b^4)*e^4*sqrt(e*cot(c + d*x)))/(4*b^3*(a^2 + b^2)^2*d) + (a^2*e^2*(e*cot(c + d*x))^(5/2))/(2*b*(a^2 + b^2)*d*(a + b*cot(c + d*x))^2) + (a^2*(5*a^2 + 13*b^2)*e^3*(e*cot(c + d*x))^(3/2))/(4*b^2*(a^2 + b^2)^2*d*(a + b*cot(c + d*x))) - ((a + b)*(a^2 - 4*a*b + b^2)*e^(9/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d) + ((a + b)*(a^2 - 4*a*b + b^2)*e^(9/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d)]
    @test_int [1/(a + b*cot(c + d*x))^3*(e*cot(c + d*x))^(7/2), x, 16, -((a^(3/2)*(3*a^4 + 6*a^2*b^2 + 35*b^4)*e^(7/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(4*b^(5/2)*(a^2 + b^2)^3*d)) + ((a + b)*(a^2 - 4*a*b + b^2)*e^(7/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) - ((a + b)*(a^2 - 4*a*b + b^2)*e^(7/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) + (a^2*e^2*(e*cot(c + d*x))^(3/2))/(2*b*(a^2 + b^2)*d*(a + b*cot(c + d*x))^2) + (a^2*(3*a^2 + 11*b^2)*e^3*sqrt(e*cot(c + d*x)))/(4*b^2*(a^2 + b^2)^2*d*(a + b*cot(c + d*x))) + ((a - b)*(a^2 + 4*a*b + b^2)*e^(7/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d) - ((a - b)*(a^2 + 4*a*b + b^2)*e^(7/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d)]
    @test_int [1/(a + b*cot(c + d*x))^3*(e*cot(c + d*x))^(5/2), x, 16, -((sqrt(a)*(a^4 + 18*a^2*b^2 - 15*b^4)*e^(5/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(4*b^(3/2)*(a^2 + b^2)^3*d)) - ((a - b)*(a^2 + 4*a*b + b^2)*e^(5/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) + ((a - b)*(a^2 + 4*a*b + b^2)*e^(5/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) + (a^2*e^2*sqrt(e*cot(c + d*x)))/(2*b*(a^2 + b^2)*d*(a + b*cot(c + d*x))^2) - (a*(a^2 + 9*b^2)*e^2*sqrt(e*cot(c + d*x)))/(4*b*(a^2 + b^2)^2*d*(a + b*cot(c + d*x))) + ((a + b)*(a^2 - 4*a*b + b^2)*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d) - ((a + b)*(a^2 - 4*a*b + b^2)*e^(5/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d)]
    @test_int [1/(a + b*cot(c + d*x))^3*(e*cot(c + d*x))^(3/2), x, 16, -(((3*a^4 - 26*a^2*b^2 + 3*b^4)*e^(3/2)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(4*sqrt(a)*sqrt(b)*(a^2 + b^2)^3*d)) - ((a + b)*(a^2 - 4*a*b + b^2)*e^(3/2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) + ((a + b)*(a^2 - 4*a*b + b^2)*e^(3/2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) - (a*e*sqrt(e*cot(c + d*x)))/(2*(a^2 + b^2)*d*(a + b*cot(c + d*x))^2) - ((3*a^2 - 5*b^2)*e*sqrt(e*cot(c + d*x)))/(4*(a^2 + b^2)^2*d*(a + b*cot(c + d*x))) - ((a - b)*(a^2 + 4*a*b + b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d) + ((a - b)*(a^2 + 4*a*b + b^2)*e^(3/2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d)]
    @test_int [1/(a + b*cot(c + d*x))^3*(e*cot(c + d*x))^(1/2), x, 16, (sqrt(b)*(15*a^4 - 18*a^2*b^2 - b^4)*sqrt(e)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(4*a^(3/2)*(a^2 + b^2)^3*d) + ((a - b)*(a^2 + 4*a*b + b^2)*sqrt(e)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) - ((a - b)*(a^2 + 4*a*b + b^2)*sqrt(e)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d) + (b*sqrt(e*cot(c + d*x)))/(2*(a^2 + b^2)*d*(a + b*cot(c + d*x))^2) + (b*(7*a^2 - b^2)*sqrt(e*cot(c + d*x)))/(4*a*(a^2 + b^2)^2*d*(a + b*cot(c + d*x))) - ((a + b)*(a^2 - 4*a*b + b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d) + ((a + b)*(a^2 - 4*a*b + b^2)*sqrt(e)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d)]
    @test_int [1/(a + b*cot(c + d*x))^3/(e*cot(c + d*x))^(1/2), x, 16, -((b^(3/2)*(35*a^4 + 6*a^2*b^2 + 3*b^4)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(4*a^(5/2)*(a^2 + b^2)^3*d*sqrt(e))) + ((a + b)*(a^2 - 4*a*b + b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d*sqrt(e)) - ((a + b)*(a^2 - 4*a*b + b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d*sqrt(e)) - (b^2*sqrt(e*cot(c + d*x)))/(2*a*(a^2 + b^2)*d*e*(a + b*cot(c + d*x))^2) - (b^2*(11*a^2 + 3*b^2)*sqrt(e*cot(c + d*x)))/(4*a^2*(a^2 + b^2)^2*d*e*(a + b*cot(c + d*x))) + ((a - b)*(a^2 + 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d*sqrt(e)) - ((a - b)*(a^2 + 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d*sqrt(e))]
    @test_int [1/(a + b*cot(c + d*x))^3/(e*cot(c + d*x))^(3/2), x, 17, (b^(5/2)*(63*a^4 + 46*a^2*b^2 + 15*b^4)*atan((sqrt(b)*sqrt(e*cot(c + d*x)))/(sqrt(a)*sqrt(e))))/(4*a^(7/2)*(a^2 + b^2)^3*d*e^(3/2)) - ((a - b)*(a^2 + 4*a*b + b^2)*atan(1 - (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d*e^(3/2)) + ((a - b)*(a^2 + 4*a*b + b^2)*atan(1 + (sqrt(2)*sqrt(e*cot(c + d*x)))/sqrt(e)))/(sqrt(2)*(a^2 + b^2)^3*d*e^(3/2)) + (8*a^4 + 31*a^2*b^2 + 15*b^4)/(4*a^3*(a^2 + b^2)^2*d*e*sqrt(e*cot(c + d*x))) - b^2/(2*a*(a^2 + b^2)*d*e*sqrt(e*cot(c + d*x))*(a + b*cot(c + d*x))^2) - (b^2*(13*a^2 + 5*b^2))/(4*a^2*(a^2 + b^2)^2*d*e*sqrt(e*cot(c + d*x))*(a + b*cot(c + d*x))) + ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) - sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d*e^(3/2)) - ((a + b)*(a^2 - 4*a*b + b^2)*log(sqrt(e) + sqrt(e)*cot(c + d*x) + sqrt(2)*sqrt(e*cot(c + d*x))))/(2*sqrt(2)*(a^2 + b^2)^3*d*e^(3/2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*cot(e+f*x))^n*(a+b*cot(e+f*x))^m*with*m*symbolic=#

    @test_int [(a + b*cot(c + d*x))^n, x, 5, -((b*(a + b*cot(c + d*x))^(1 + n)*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (a + b*cot(c + d*x))/(a - sqrt(-b^2))))/(2*sqrt(-b^2)*(a - sqrt(-b^2))*d*(1 + n))) + (b*(a + b*cot(c + d*x))^(1 + n)*HypergeometricFunctions._₂F₁(1, 1 + n, 2 + n, (a + b*cot(c + d*x))/(a + sqrt(-b^2))))/(2*sqrt(-b^2)*(a + sqrt(-b^2))*d*(1 + n))]

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*(d*tan(e+f*x))^n*(a+b*cot(e+f*x))^m=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(d*tan(e+f*x))^n*(a+im*a*cot(e+f*x))^m=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(d*tan(e+f*x))^n*(a+a*cot(e+f*x))^m=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*tan(e+f*x))^n*(a+b*cot(e+f*x))^m=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*tan(e+f*x))^n*(a+b*cot(e+f*x))^m*with*m*symbolic=#

    @test_int [(a + b*cot(e + f*x))^m*(d*tan(e + f*x))^n, x, 8, -((AppellF1(1 - n, -m, 1, 2 - n, -((b*cot(e + f*x))/a), (-im)*cot(e + f*x))*cot(e + f*x)*(a + b*cot(e + f*x))^m*(d*tan(e + f*x))^n)/((1 + (b*cot(e + f*x))/a)^m*(2*f*(1 - n)))) - (AppellF1(1 - n, -m, 1, 2 - n, -((b*cot(e + f*x))/a), im*cot(e + f*x))*cot(e + f*x)*(a + b*cot(e + f*x))^m*(d*tan(e + f*x))^n)/((1 + (b*cot(e + f*x))/a)^m*(2*f*(1 - n)))]

    #= ::Title::Closed:: =#
    #=Integrands*of*the*form*(a+b*cot(e+f*x))^m*(c+d*cot(e+f*x))^n=#

    #= ::Section:: =#
    #=Integrands*of*the*form*(a+a*im*cot(e+f*x))^m*(c-c*im*cot(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+a*im*cot(e+f*x))^m*(c+d*cot(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+a*im*cot(e+f*x))^m*(c+d*cot(e+f*x))^(n/2)=#

    @test_int [(1 + im*cot(c + d*x))/sqrt(a + b*cot(c + d*x)), x, 3, (2*im*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/(sqrt(a - im*b)*d)]
    @test_int [(1 - im*cot(c + d*x))/sqrt(a + b*cot(c + d*x)), x, 3, -((2*im*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/(sqrt(a + im*b)*d))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(a+b*cot(e+f*x))^m*(c+d*cot(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*cot(e+f*x))^m*(c+d*cot(e+f*x))^n=#

    #= ::Subsubsection:: =#
    #=m>0=#

    #= ::Subsubsection::Closed:: =#
    #=m<0=#

    @test_int [(A + B*cot(c + d*x))/(a + b*cot(c + d*x)), x, 2, ((a*A + b*B)*x)/(a^2 + b^2) - ((A*b - a*B)*log(b*cos(c + d*x) + a*sin(c + d*x)))/((a^2 + b^2)*d)]

    @test_int [(A + B*cot(c + d*x))/(a + b*cot(c + d*x))^2, x, 3, ((a^2*A - A*b^2 + 2*a*b*B)*x)/(a^2 + b^2)^2 + (A*b - a*B)/((a^2 + b^2)*d*(a + b*cot(c + d*x))) - ((2*a*A*b - a^2*B + b^2*B)*log(b*cos(c + d*x) + a*sin(c + d*x)))/((a^2 + b^2)^2*d)]

    @test_int [(A + B*cot(c + d*x))/(a + b*cot(c + d*x))^3, x, 4, ((a^3*A - 3*a*A*b^2 + 3*a^2*b*B - b^3*B)*x)/(a^2 + b^2)^3 + (A*b - a*B)/(2*(a^2 + b^2)*d*(a + b*cot(c + d*x))^2) + (2*a*A*b - a^2*B + b^2*B)/((a^2 + b^2)^2*d*(a + b*cot(c + d*x))) - ((3*a^2*A*b - A*b^3 - a^3*B + 3*a*b^2*B)*log(b*cos(c + d*x) + a*sin(c + d*x)))/((a^2 + b^2)^3*d)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(a+b*cot(e+f*x))^m*(c+d*cot(e+f*x))^(n/2)=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [(A + B*cot(c + d*x))*(a + b*cot(c + d*x))^(5/2), x, 10, ((a - im*b)^(5/2)*(im*A + B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/d - ((a + im*b)^(5/2)*(im*A - B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/d - (2*(2*a*A*b + a^2*B - b^2*B)*sqrt(a + b*cot(c + d*x)))/d - (2*(A*b + a*B)*(a + b*cot(c + d*x))^(3/2))/(3*d) - (2*B*(a + b*cot(c + d*x))^(5/2))/(5*d)]
    @test_int [(A + B*cot(c + d*x))*(a + b*cot(c + d*x))^(3/2), x, 9, ((a - im*b)^(3/2)*(im*A + B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/d - ((a + im*b)^(3/2)*(im*A - B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/d - (2*(A*b + a*B)*sqrt(a + b*cot(c + d*x)))/d - (2*B*(a + b*cot(c + d*x))^(3/2))/(3*d)]
    @test_int [(A + B*cot(c + d*x))*(a + b*cot(c + d*x))^(1/2), x, 8, (sqrt(a - im*b)*(im*A + B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/d - (sqrt(a + im*b)*(im*A - B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/d - (2*B*sqrt(a + b*cot(c + d*x)))/d]

    @test_int [(-a + b*cot(c + d*x))*(a + b*cot(c + d*x))^(5/2), x, 10, -(((im*a - b)*(a - im*b)^(5/2)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/d) + ((a + im*b)^(5/2)*(im*a + b)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/d + (2*b*(a^2 + b^2)*sqrt(a + b*cot(c + d*x)))/d - (2*b*(a + b*cot(c + d*x))^(5/2))/(5*d)]
    @test_int [(-a + b*cot(c + d*x))*(a + b*cot(c + d*x))^(3/2), x, 13, (b*(a^2 + b^2)*atanh((sqrt(a + sqrt(a^2 + b^2)) - sqrt(2)*sqrt(a + b*cot(c + d*x)))/sqrt(a - sqrt(a^2 + b^2))))/(sqrt(2)*sqrt(a - sqrt(a^2 + b^2))*d) - (b*(a^2 + b^2)*atanh((sqrt(a + sqrt(a^2 + b^2)) + sqrt(2)*sqrt(a + b*cot(c + d*x)))/sqrt(a - sqrt(a^2 + b^2))))/(sqrt(2)*sqrt(a - sqrt(a^2 + b^2))*d) - (2*b*(a + b*cot(c + d*x))^(3/2))/(3*d) + (b*(a^2 + b^2)*log(a + sqrt(a^2 + b^2) + b*cot(c + d*x) - sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*sqrt(a + b*cot(c + d*x))))/(2*sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*d) - (b*(a^2 + b^2)*log(a + sqrt(a^2 + b^2) + b*cot(c + d*x) + sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*sqrt(a + b*cot(c + d*x))))/(2*sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*d)]
    @test_int [(-a + b*cot(c + d*x))*(a + b*cot(c + d*x))^(1/2), x, 13, (b*sqrt(a^2 + b^2)*atanh((sqrt(a + sqrt(a^2 + b^2)) - sqrt(2)*sqrt(a + b*cot(c + d*x)))/sqrt(a - sqrt(a^2 + b^2))))/(sqrt(2)*sqrt(a - sqrt(a^2 + b^2))*d) - (b*sqrt(a^2 + b^2)*atanh((sqrt(a + sqrt(a^2 + b^2)) + sqrt(2)*sqrt(a + b*cot(c + d*x)))/sqrt(a - sqrt(a^2 + b^2))))/(sqrt(2)*sqrt(a - sqrt(a^2 + b^2))*d) - (2*b*sqrt(a + b*cot(c + d*x)))/d - (b*sqrt(a^2 + b^2)*log(a + sqrt(a^2 + b^2) + b*cot(c + d*x) - sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*sqrt(a + b*cot(c + d*x))))/(2*sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*d) + (b*sqrt(a^2 + b^2)*log(a + sqrt(a^2 + b^2) + b*cot(c + d*x) + sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*sqrt(a + b*cot(c + d*x))))/(2*sqrt(2)*sqrt(a + sqrt(a^2 + b^2))*d)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [(A + B*cot(c + d*x))/(a + b*cot(c + d*x))^(1/2), x, 7, ((im*A + B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/(sqrt(a - im*b)*d) - ((im*A - B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/(sqrt(a + im*b)*d)]
    @test_int [(A + B*cot(c + d*x))/(a + b*cot(c + d*x))^(3/2), x, 8, ((im*A + B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/((a - im*b)^(3/2)*d) - ((im*A - B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/((a + im*b)^(3/2)*d) + (2*(A*b - a*B))/((a^2 + b^2)*d*sqrt(a + b*cot(c + d*x)))]
    @test_int [(A + B*cot(c + d*x))/(a + b*cot(c + d*x))^(5/2), x, 9, ((im*A + B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/((a - im*b)^(5/2)*d) - ((im*A - B)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/((a + im*b)^(5/2)*d) + (2*(A*b - a*B))/(3*(a^2 + b^2)*d*(a + b*cot(c + d*x))^(3/2)) + (2*(2*a*A*b - a^2*B + b^2*B))/((a^2 + b^2)^2*d*sqrt(a + b*cot(c + d*x)))]

    @test_int [(-a + b*cot(c + d*x))/(a + b*cot(c + d*x))^(1/2), x, 7, -(((im*a - b)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/(sqrt(a - im*b)*d)) + ((im*a + b)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/(sqrt(a + im*b)*d)]
    @test_int [(-a + b*cot(c + d*x))/(a + b*cot(c + d*x))^(3/2), x, 8, -(((im*a - b)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/((a - im*b)^(3/2)*d)) + ((im*a + b)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/((a + im*b)^(3/2)*d) - (4*a*b)/((a^2 + b^2)*d*sqrt(a + b*cot(c + d*x)))]
    @test_int [(-a + b*cot(c + d*x))/(a + b*cot(c + d*x))^(5/2), x, 9, -(((im*a - b)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a - im*b)))/((a - im*b)^(5/2)*d)) + ((im*a + b)*atanh(sqrt(a + b*cot(c + d*x))/sqrt(a + im*b)))/((a + im*b)^(5/2)*d) - (4*a*b)/(3*(a^2 + b^2)*d*(a + b*cot(c + d*x))^(3/2)) - (2*b*(3*a^2 - b^2))/((a^2 + b^2)^2*d*sqrt(a + b*cot(c + d*x)))]
    end

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.2.1 (a+b cot)^m (c+d cot)^n.jl


text:
module MusicTransformer

using Flux
using Transformers

include("layers.jl")
include("models.jl")
include("musicencoders.jl")
include("generate.jl")
include("pretrain.jl")

end

ext: jl


text:
module CellMLToolkit

using MathML

using ModelingToolkit
using SymbolicUtils: FnType, Sym, operation, arguments
using EzXML
using Memoize

include("structures.jl")
include("accessors.jl")
include("components.jl")
include("import.jl")

function read_cellml(path::AbstractString, tspan)
    @warn "read_cellml is deprecated, please use CellModel"
    ml = CellModel(path)
    ODEProblem(ml, tspan)
end

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

export CellModel, ODEProblem
export read_cellml, parse_cellml
export list_params, list_states
export readxml, getxml, getsys
export update_list!

getsys(ml::CellModel) = ml.sys

"""
    constructs a CellModel struct for the CellML model defined in path
"""
function CellModel(path::AbstractString)
    doc = load_cellml(path)
    CellModel(doc, process_components(doc))
end

list_params(ml::CellModel) = find_sys_p(ml.doc, ml.sys)
list_states(ml::CellModel) = find_sys_u0(ml.doc, ml.sys)

import ModelingToolkit.ODEProblem

"""
    ODEProblem constructs an ODEProblem from a CellModel
"""
function ODEProblem(ml::CellModel, tspan;
        jac=false, level=1, p=last.(list_params(ml)), u0=last.(list_states(ml)))
    ODEProblem(ml.sys, u0, tspan, p; jac=jac)
end

function update_list!(l, sym, val)
    i = findfirst(isequal(sym), Symbol.(first.(l)))
    if i != nothing
        l[i] = (first(l[i]) => val)
    else
        @warn "symbol $sym not found"
    end
end

end # module


text:
struct SampleBased <: ConversionTrait end

function convert_arguments(::SampleBased, args::NTuple{N,AbstractVector{<:Number}}) where {N}
    return args
end

function convert_arguments(P::SampleBased, positions::Vararg{AbstractVector})
    return convert_arguments(P, positions)
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module MCObservables

import Base: show, <<, push!, mean, var, count, isempty, merge, merge!, zero, zeros, sum
export MCObservable, ScalarObservable, VectorObservable
export mean, var, stderror, confidence_interval
export show, dump_plot
export merge, merge!

export confidence_interval

using Distributions

abstract MCObservable
abstract ScalarObservable <: MCObservable
abstract VectorObservable <: MCObservable

isempty(obs::MCObservable) = count(obs) == 0

zero{Obs<:MCObservable}(::Type{Obs}) = Obs()
zero(obs::MCObservable) = zero(typeof(obs))
function zeros{Obs<:MCObservable}(::Type{Obs}, dim...)
  reshape(Obs[Obs() for i in 1:prod(dim)], dim)
end

function show(io::IO, obs::MCObservable)
  if !isempty(obs)
    print(io, mean(obs), " +/- ", stderror(obs))
  else
    print(io, "No Entries")
  end
end

function dump_plot(io::IO, obs::MCObservable)
  print(io, mean(obs), " ", stderror(obs))
end

const confidence_rate_1sigma = 0.5erf(0.5sqrt(2.0))

include("util.jl")
include("observableset.jl")
include("parsesigma.jl")
include("tiny.jl")
include("tinyvector.jl")
include("simple.jl")
include("simplevector.jl")
include("binning.jl")
include("binningvector.jl")
include("jackknife.jl")
include("jackknifevector.jl")

## these three definitions are needed to resolve ambiguousness with Base.<<(Any,Integer)
<<(obs::MCObservable, x::Int32) = push!(obs,x)
<<(obs::MCObservable, x::Int64) = push!(obs,x)
<<(obs::MCObservable, x::Integer) = push!(obs,x)

<<(obs::MCObservable, x) = push!(obs,x)

end # module


text:
apply(ψ::Vector{<:Complex}, m::AbstractMatrix) = m*ψ

"""Tailored apply for XGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,XGate}, N::Int)
    i = cg.iwire[1]
    χ = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = χ[:, [2, 1], :]
    return reshape(χ, :)
end

"""Tailored apply for YGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,YGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= -im.*ψs[:, 2, :]
    χ[:, 2, :] .=  im.*ψs[:, 1, :]
    return reshape(χ, :)
end

"""Tailored apply for ZGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,ZGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .=   ψs[:, 1, :]
    χ[:, 2, :] .= .-ψs[:, 2, :]
    return reshape(χ, :)
end

"""Tailored apply for HadamardGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,HadamardGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= (ψs[:, 1, :] .+ ψs[:, 2, :]) ./ sqrt(2)
    χ[:, 2, :] .= (ψs[:, 1, :] .- ψs[:, 2, :]) ./ sqrt(2)
    return reshape(χ, :)
end

"""Tailored apply for SGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,SGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= ψs[:, 1, :]
    χ[:, 2, :] .= im .* ψs[:, 2, :]
    return reshape(χ, :)
end

"""Tailored apply for SdagGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,SdagGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= ψs[:, 1, :]
    χ[:, 2, :] .= -im .* ψs[:, 2, :]
    return reshape(χ, :)
end

"""Tailored apply for TGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,TGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= ψs[:, 1, :]
    χ[:, 2, :] .= Base.exp(im*π/4) .* ψs[:, 2, :]
    return reshape(χ, :)
end

"""Tailored apply for TdagGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,TdagGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= ψs[:, 1, :]
    χ[:, 2, :] .= Base.exp(-im*π/4) .* ψs[:, 2, :]
    return reshape(χ, :)
end

"""Tailored apply for RzGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,RzGate}, N::Int)
    i = cg.iwire[1]
    θ = cg.gate.θ[]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= Base.exp(-im*θ/2).*ψs[:, 1, :]
    χ[:, 2, :] .= Base.exp( im*θ/2).*ψs[:, 2, :]
    return reshape(χ, :)
end

"""Tailored apply for PhaseShiftGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,PhaseShiftGate}, N::Int)
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    χ[:, 1, :] .= ψs[:, 1, :]
    χ[:, 2, :] .= Base.exp(im*cg.gate.ϕ[]) .* ψs[:, 2, :]
    return reshape(χ, :)
end

"""Tailored apply for a general single qubit gate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{1,T}, N::Int) where {T<:AbstractGate}
    i = cg.iwire[1]
    ψs = reshape(ψ, 2^(i-1), 2, 2^(N-i))
    χ = similar(ψs)
    U = matrix(cg.gate)
    χ[:, 1, :] .= U[1, 1] .* ψs[:, 1, :] .+ U[1, 2] .* ψs[:, 2, :]
    χ[:, 2, :] .= U[2, 1] .* ψs[:, 1, :] .+ U[2, 2] .* ψs[:, 2, :]
    return reshape(χ, :)
end
#

"""Tailored apply for SwapGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{2,SwapGate}, N::Int)
    i, j = cg.iwire
    i, j = i < j ? (i, j) : (j, i) # sort them

    ψs = reshape(ψ, 2^(i-1), 2, 2^(j-i-1), 2, 2^(N-j))
    χ = similar(ψs)

    χ[:, 1, :, 1, :] .= ψs[:, 1, :, 1, :]
    χ[:, 2, :, 2, :] .= ψs[:, 2, :, 2, :]
    χ[:, 1, :, 2, :] .= ψs[:, 2, :, 1, :]
    χ[:, 2, :, 1, :] .= ψs[:, 1, :, 2, :]

    return reshape(χ, :)
end

"""Tailored apply for a general ControlledGate"""
function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{M,ControlledGate{G}}, N::Int) where {M,G}
    χ = copy(ψ) # TODO: copy only what's needed
    χ = reshape(χ, fill(2, N)...)
    T = target_wires(cg.gate)
    itarget  = cg.iwire[1:T]
    icontrol = cg.iwire[T+1:M]

    new_iwire = Tuple(i - sum(icontrol .< i) for i in itarget)

    slice_target = Tuple(i in icontrol ? 2 : Colon() for i in 1:N)
    new_cg = CircuitGate(new_iwire, cg.gate.U)

    χ[slice_target...] = reshape(apply(reshape(χ[slice_target...], :), new_cg),
                                 fill(2, N-M+T)...)
    return reshape(χ, :)
end

"""
    apply(ψ::Vector{<:Complex}, cg::CircuitGate{M,G}) where {M,G}

Apply a [`CircuitGate`](@ref) to a quantum state vector `ψ`.
# Examples

```jldoctest
julia> cg = circuit_gate(1, HadamardGate());
julia> ψ = [1 0];
julia> apply(ψ, cg)
2-element Array{Complex{Float64},1}:
 0.7071067811865475 + 0.0im
 0.7071067811865475 + 0.0im
```
"""
function apply(ψ::Vector{<:Complex}, cg::CircuitGate{M,G}) where {M,G}
    l = length(ψ)::Int
    N = Qaintessent.intlog2(l)
    l == 2^N || error("Vector length must be a power of 2")
    req_wires(cg) <= N || error("CircuitGate requires a minimum of $(req_wires(cg)) qubits, input vector `ψ` has $N qubits")
    _apply(ψ, cg, N)
end

function _apply(ψ::Vector{<:Complex}, cg::CircuitGate{M,G}, N::Int) where {M,G}
    U = matrix(cg.gate)
    ψs = reshape(ψ, fill(2, N)...)
    χ = similar(ψs)
    it = binary_digits(M, 0)
    m = binary_digits(M, 0)
    for i in 1:2^M
        binary_digits!(it, i - 1)
        # cannot use .= here since broadcasting fails for scalar numbers
        χ[sliced_index(it, cg.iwire, N)...] = sum(U[i, j] .* ψs[sliced_index(binary_digits!(m, j - 1), cg.iwire, N)...] for j in 1:2^M)
    end
    return reshape(χ, :)
end

"""
    apply(ψ::Vector{<:Complex}, cgs::Vector{<:CircuitGate})

Apply a sequence of [`CircuitGate`](@ref)(s) to a quantum state vector `ψ`.

# Examples

```jldoctest
julia> cgs = [circuit_gate(1, HadamardGate()),
                circuit_gate(1, X),
                circuit_gate(1, Y)];
julia> ψ = [1 0];
julia> apply(ψ, cgs)
2-element Array{Complex{Float64},1}:
 0.0 - 0.7071067811865475im
 0.0 + 0.7071067811865475im
```
"""
function apply(ψ::Vector{<:Complex}, cgs::Vector{<:CircuitGate})
    length(cgs) != 0 || error("Vector of length 0 cannot be applied")
    l = length(ψ)
    N = Qaintessent.intlog2(l)
    l == 2^N || error("Vector length must be a power of 2")
    req = maximum(req_wires.(cgs))
    req <= N || error("CircuitGates require a minimum of $req qubits, input vector `ψ` has $N qubits")
    for cg in cgs
        ψ = _apply(ψ, cg, N)
    end
    return ψ
end

"""
    apply(ψ::Vector{<:Complex}, m::Moment)

returns state vector of `N` qubits after applying a `Moment{N}` object to a quantum state vector of `N` qubits `ψ`
"""
function apply(ψ::Vector{<:Complex}, m::Moment)
    Nmoment = req_wires(m)
    l = length(ψ)
    N = Qaintessent.intlog2(l)
    l == 2^N || error("Vector length must be a power of 2")
    Nmoment <= N || error("Moment affecting $Nmoment qubits applied to $N qubits")
    _apply(ψ, m, N)
end

function _apply(ψ::Vector{<:Complex}, m::Moment, N::Int)
    for gate in m
        ψ = _apply(ψ, gate, N)
    end
    return ψ
end

function apply(ψ::Vector{<:Complex}, m::Vector{Moment})
    length(m) != 0 || error("Vector of length 0 cannot be applied")
    Nmoment = maximum(req_wires.(m))
    l = length(ψ)
    N = Qaintessent.intlog2(l)
    l == 2^N || error("Vector length must be a power of 2")
    Nmoment <= N || error("Moment affecting $Nmoment qubits applied to $N qubits")
    for moment in m
        ψ = _apply(ψ, moment, N)
    end
    return ψ
end

"""
    apply(ψ::Vector{<:Complex}, m::MeasurementOperator{M,G}) where {M,G<:AbstractGate}

returns state vector of `N` qubits after applying a `Moment{N}` object to a quantum state vector of `N` qubits `ψ`
"""

function apply(ψ::Vector{<:Complex}, m::MeasurementOperator{M,G}) where {M,G<:AbstractGate}
    Nmoment = num_wires(m)
    l = length(ψ)::Int
    N = Qaintessent.intlog2(l)
    l == 2^N || error("Vector length must be a power of 2")
    Nmoment <= N || error("MeasurementOperator affecting $Nmoment qubits applied to $N qubits")
    c = circuit_gate((m.iwire...), m.operator)
    _apply(ψ, c, N)
end

function apply(ψ::Vector{<:Complex}, m::MeasurementOperator{M,G}) where {M,G<:AbstractMatrix}
    Nmoment = num_wires(m)
    l = length(ψ)::Int
    N = Qaintessent.intlog2(l)
    l == 2^N || error("Vector length must be a power of 2")
    Nmoment <= N || error("MeasurementOperator affecting $Nmoment qubits applied to $N qubits")
    apply(ψ, m.operator)
end

"""
    apply(ψ::Vector{<:Complex}, c::Circuit{N}) where {N}

returns list of expectation values from measurement operators in `c.meas` after applying circuit gates in `c.cgc` on state vector of `N` qubits `ψ`
"""
function apply(ψ::Vector{<:Complex}, c::Circuit{N}) where {N}
    length(ψ)::Int == 2^N || error("Size of vector `ψ` must match Circuit size of $(2^N)")
    length(c) != 0 || error("Circuit does not contain any gates")
    length(c.meas) != 0 || error("Circuit does not contain any measurement operators")
    ψl = copy(ψ)
    for moment in c.moments
        ψl = _apply(ψl, moment, N)
    end
    ψr = apply.((ψl,), c.meas)
    return real.(dot.((ψl,), ψr))
end

(c::Circuit)(ψ) = apply(ψ, c)


text:
# R formula to julia Formula

# a special wrapper for rcopy specialized for formulas
function rcopy_formula(s::Ptr{S}) where S<:Sxp
    r = rcopy(s)
    if isa(r, Number) && isinteger(r)
        Int(r)
    else
        r
    end
end

function rcopy(::Type{Expr}, l::Ptr{LangSxp})
    op = rcopy(Symbol, l[1])
    if op == :(:)
        op = :&
    end
    args = cdr(l)
    if op == Symbol("(")
        f = rcopy(Expr, l[2])
    else
        cargs = []
        for (t, s) in pairs(args)
            if isNull(t)
                push!(cargs, rcopy_formula(s))
            else
                push!(cargs, Expr(:(=), rcopy_formula(t), rcopy_formula(s)))
            end
        end
        f = Expr(:call, op, cargs...)
    end
    # unwind these opeators
    if op in (:+, :*, :&) && isa(f.args[2], Expr) && f.args[2].args[1] == op
        f = Expr(:call, op, f.args[2].args[2:end]..., f.args[3])
    end
    f
end

function rcopy(::Type{Formula}, l::Ptr{LangSxp})
    ex_orig = rcopy(Expr, l)
    ex = parse!(copy(ex_orig))
    lhs, rhs = ex.args[2:3]
    Formula(ex_orig, ex, lhs, rhs)
end

# julia Formula to R formula

# formula
function sexp_formula(e::Expr)
    e.head == :call || error("invalid formula object")
    op = e.args[1]
    if op == :&
        op = :(:)
    end
    if length(e.args) > 3 && op in (:+, :*, :(:))
        rlang_p(op, sexp_formula(Expr(e.head, e.args[1:end-1]...)), sexp_formula(e.args[end]))
    else
        rlang_p(op, map(sexp_formula, e.args[2:end])...)
    end
end
sexp_formula(e::Symbol) = sexp(SymSxp, e)
sexp_formula(n::Integer) = sexp(RClass{:numeric}, Float64(n))
sexp_formula(n::Number) = sexp(n)

# R formula objects
function sexp(::Type{RClass{:formula}}, f::Formula)
    s = protect(sexp_formula(f.ex_orig == :() ? f.ex : f.ex_orig))
    try
        setattrib!(s, Const.ClassSymbol, "formula")
        setattrib!(s, ".Environment", Const.GlobalEnv)
    finally
        unprotect(1)
    end
    s
end


text:
#
# types.jl --
#
# Type definitions for local filters.
#
#------------------------------------------------------------------------------
#
# This file is part of the `LocalFilters.jl` package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2020, Éric Thiébaut.
#

"""
All neighborhoods are instances of a type derived from `Neighborhood`.
"""
abstract type Neighborhood{N} end

"""

    RectangularBox(start, stop)

yields a neighborhood which is a rectangular (Cartesian) box defined by the
bounds of the multi-dimensional indices in the box.

Another possibility is to specify the dimensions of the box and the offsets of
its central element:

    RectangularBox(dims, offs)

with `dims` a `N`-tuple of dimensions and `offs` either a `N`-tuple of indices
of an instance of `CartesianIndex{N}`.

A `RectangularBox` can also be defined by the index ranges along all the
dimensions.  For example:

    RectangularBox(-3:3, 0, -2:1)
    RectangularBox((-3:3, 0, -2:1))

both yield a 3-dimensional `RectangularBox` of size `7×1×4` and whose first
index varies on `-3:3`, its second index is `0` while its third index varies on
`-2:1`.

Finally, a `RectangularBox` can be defined as:

    RectangularBox(R)

where `R` is an instance of `CartesianIndices`.

"""
struct RectangularBox{N} <: Neighborhood{N}
    # For a RectangularBox, the most used operation is getting its limits, so
    # this is what we store.
    start::CartesianIndex{N}
    stop::CartesianIndex{N}
end

"""

A `Kernel` can be used to define a weighted neighborhood (for weighted local
average or for convolution) or a structuring element (for mathematical
morphology).  It is a rectangular array of coefficients over a, possibly
off-centered, rectangular neighborhood.  In general, it is sufficient to
specify `::Kernel{T,N}` in the signature of methods, with `T` the type of the
coefficients and `N` the number of dimensions (the third parameter `A` of the
type is to fully qualify the type of the array of coefficients).

A kernel is built as:

    B = Kernel([T,] C, start=defaultstart(C))

where `C` is the array of coefficients (which can be retrieved by `coefs(B)`)
and `start` the initial `CartesianIndex` for indexing the kernel (which can be
retrieved by `initialindex(B)`).  The `start` parameter let the caller choose
an arbitrary origin for the kernel coefficients; when a filter is applied, the
following mapping is assumed:

    B[k] ≡ C[k + off]

where `off = initialindex(C) - initialindex(B)`.

If `start` is omitted, its value is set so that the *origin* (whose index is
`zero(CartesianIndex{N})` with `N` the number of dimensions) of the kernel
indices is at the geometric center of the array of coefficients (see
[`LocalFilters.defaultstart`](@ref)).  Optional argument `T` is to impose the
type of the coefficients.

To convert the element type of the coefficients of an existing kernel, do:

    Kernel(T, K)

which yields a kernel whose coefficients are those of the kernel `K`
converted to type `T`.

It is also possible to convert instances of [`RectangularBox`](@ref) into a
kernel with boolean coefficients by calling:

    Kernel(B)

where `B` is the neighborhood to convert into an instance of `Kernel`.

"""
struct Kernel{T,N,A<:AbstractArray{T,N}} <: Neighborhood{N}
    coefs::A
    offset::CartesianIndex{N}
    start::CartesianIndex{N}
    stop::CartesianIndex{N}
    function Kernel{T,N,A}(C::A, start::CartesianIndex{N}
                           ) where {T,N,A<:AbstractArray{T,N}}
        # In local filtering operations, the kernel `B` will be used as
        #
        #     A[i-k]⋄B[k]    (∀ k ∈ [kmin, kmax])
        #
        # with `A` the source, `i` the destination index, `⋄` a given binary
        # operation and `k` (here, indices and offsets are multi-dimensional).
        # The offset between the indices in the kernel `B` and those in the
        # array of coefficients `C` is such that:
        #
        #     kmin + off = jmin = initialindex(C)
        #     kmax + off = jmax = finalindex(C)
        #
        # Hence:
        #
        #     kmin = start
        #     off = jmin - kmin
        #     kmax = jmax - off
        #
        # gives the [limits](@ref) and [offset](@ref) for index `k` in the
        # kernel `B`.
        @assert eltype(C) === T
        @assert ndims(C) == N
        jmin, jmax = limits(C)
        off = jmin - start
        return new{T,N,A}(C, off, start, jmax - off)
    end
end

"""

`IndexInterval` is an union of the types of any argument suitable to specify an
interval of indices along a dimension (that is, an integer or an integer valued
unit range).

"""
const IndexInterval = Union{Integer,AbstractUnitRange{<:Integer}}

"""

`Dimensions{N}` is an `N`-tuple of integers, that is the type of an argument
suitable for specifying a list of dimensions.

"""
const Dimensions{N} = NTuple{N,Integer}

"""

`UnitIndexRange` is any integer valued range with unit step, that is the type
of an argument suitable for specifying a range of contiguous indices.

"""
const UnitIndexRange = AbstractUnitRange{<:Integer}

"""

`UnitIndexRanges{N}` is an `N`-tuple of `UnitIndexRange`, that is the type of
an argument suitable for specifying a list of Cartesian indices.

"""
const UnitIndexRanges{N} = NTuple{N,UnitIndexRange}

"""

`CartesianRegion{N}` is an union of the types of anything suitable to define a
Cartesian region in `N` dimensions.  That is, an interval of Cartesian indices
in `N` dimensions.  Methods [`initialindex`](@ref), [`finalindex`](@ref),
[`limits`](@ref) and [`cartesianregion`](@ref) can be applied to anything
whose type belongs to `CartesianRegion`.

"""
const CartesianRegion{N} = Union{NTuple{2,CartesianIndex{N}},
                                 UnitIndexRanges{N},
                                 Neighborhood{N},
                                 CartesianIndices{N}}


text:
using Knet
libknet8handle = Libdl.dlopen(Knet.libknet8)

SIZE1 = 1000
SIZE2 = 100
ITER = 100000
mat32 = KnetArray(rand(Float32,SIZE1,SIZE2))
mat32b = KnetArray(rand(Float32,SIZE1,SIZE2))
col32 = KnetArray(rand(Float32,SIZE1))
row32 = KnetArray(rand(Float32,1,SIZE2))
out32 = similar(mat32)
mat64 = KnetArray(rand(Float64,SIZE1,SIZE2))
mat64b = KnetArray(rand(Float64,SIZE1,SIZE2))
col64 = KnetArray(rand(Float64,SIZE1))
row64 = KnetArray(rand(Float64,1,SIZE2))
out64 = similar(mat64)
f32 = f64 = nothing

function cuda12test(fname, jname=fname, o...)
    println(fname)
    fcpu = eval(parse(jname))
    global f32 = Libdl.dlsym(libknet8handle, fname*"_32_12")
    @time cuda12rep(f32,mat32,mat32b,out32)
    isapprox(Array(out32),fcpu(Array(mat32),Array(mat32b))) || warn("$fname 32 mat mat")
    @time cuda12rep(f32,mat32,col32,out32)
    isapprox(Array(out32),fcpu(Array(mat32),Array(col32))) || warn("$fname 32 mat col")
    @time cuda12rep(f32,mat32,row32,out32)
    isapprox(Array(out32),fcpu(Array(mat32),Array(row32))) || warn("$fname 32 mat row")
    @time cuda12rep(f32,row32,col32,out32)
    isapprox(Array(out32),fcpu(Array(row32),Array(col32))) || warn("$fname 32 row col")
    global f64 = Libdl.dlsym(libknet8handle, fname*"_64_12")
    @time cuda12rep(f64,mat64,mat64b,out64)
    isapprox(Array(out64),fcpu(Array(mat64),Array(mat64b))) || warn("$fname 64 mat mat")
    @time cuda12rep(f64,mat64,col64,out64)
    isapprox(Array(out64),fcpu(Array(mat64),Array(col64))) || warn("$fname 64 mat col")
    @time cuda12rep(f64,mat64,row64,out64)
    isapprox(Array(out64),fcpu(Array(mat64),Array(row64))) || warn("$fname 64 mat row")
    @time cuda12rep(f64,row64,col64,out64)
    isapprox(Array(out64),fcpu(Array(row64),Array(col64))) || warn("$fname 64 row col")
end

function cuda12rep{T}(f,x::KnetArray{T},y::KnetArray{T},z::KnetArray{T})
    n = Cint(length(z))
    if size(x,1)==1; sx=size(z,1); else; sx=1; end
    if size(y,1)==1; sy=size(z,1); else; sy=1; end
    nx = length(x); ny = length(y)
    for i=1:ITER
        ccall(f,Void,(Cint,Ptr{T},Cint,Cint,Ptr{T},Cint,Cint,Ptr{T}),n,x,sx,nx,y,sy,ny,z)
    end
    Knet.@cuda(cudart,cudaDeviceSynchronize,())
    Knet.@cuda(cudart,cudaGetLastError,())
end

for f in Knet.cuda12
    isa(f,Tuple) || (f=(f,))
    cuda12test(f...)
end


text:
using Zygote, Test
using Zygote: bufferfrom
using Base.Threads: @spawn

function threads1(x)
  ch = Channel(0)
  @sync begin
    @spawn put!(ch, x^2)
    take!(ch)
  end
end

@test gradient(threads1, 5) == (10,)

function threads2(xs)
  n = length(xs)
  chunks = view(xs, 1:n÷2), view(xs, n÷2+1:n)
  p = bufferfrom([0.0, 0.0])
  @sync begin
    for i = 1:2
      @spawn begin
        s = zero(eltype(chunks[i]))
        for j = 1:length(chunks[i])
          s += chunks[i][j]
        end
        p[i] = s
      end
    end
  end
  return p[1]+p[2]
end

@test gradient(threads2, [1, 2, 3, 4]) == ([1, 1, 1, 1],)


text:
module RecurrenceAnalysis

using Distances, Statistics, LinearAlgebra, SparseArrays, DelayEmbeddings, StaticArrays
import Base.Meta.parse

const METRICS = Dict(
    "euclidean"=>Euclidean(),
    "max"=>Chebyshev(),
    "inf"=>Chebyshev(),
    "cityblock"=>Cityblock(),
    "manhattan"=>Cityblock(),
    "taxicab"=>Cityblock(),
    "min"=>Cityblock()
)
const DEFAULT_METRIC = Euclidean()
getmetric(m::Metric) = m
function getmetric(normtype::AbstractString)
    normtype = lowercase(normtype)
    !haskey(METRICS,normtype) && error("incorrect norm type. Accepted values are \""
        *join(keys(METRICS),"\", \"", "\" or \"") * "\".")
    METRICS[normtype]
end

export RecurrenceMatrix, CrossRecurrenceMatrix, JointRecurrenceMatrix,
       AbstractRecurrenceMatrix, WithinRange, NeighborNumber, FAN

export embed,
       reconstruct,
       Dataset,
       distancematrix,
       textrecurrenceplot,
       recurrenceplot,
       recurrencerate,
       recurrencestructures,
       dl_average, dl_max, dl_entropy,
       vl_average, vl_max, vl_entropy,
       rt_average, rt_max, rt_entropy,
       determinism,
       divergence,
       trend,
       laminarity,
       trappingtime,
       meanrecurrencetime,
       nmprt,
       rqa,
       sorteddistances,
       skeletonize,
       @windowed,
       rna,
       # deprecated:
       transitivity

include("matrices/distance_matrix.jl")
include("matrices/matrices.jl")
include("matrices/plot.jl")
include("matrices/skeletonization.jl")
include("rqa/histograms.jl")
include("rqa/rqa.jl")
include("rqa/radius.jl")
include("rna/graphs.jl")
include("rna/rna.jl")
include("rqa/windowed.jl")
include("deprecate.jl")

end


text:
using UtilitiesForMRI, Test

@testset "UtilitiesForMRI.jl" begin
    include("./test_kspace_sampling.jl")
    include("./test_motionpars_utils.jl")
    include("./test_nfft.jl")
    include("./test_optimizers.jl")
    include("./test_rotations.jl")
    include("./test_translations.jl")
    include("./test_updownscaling.jl")
end

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
using IterTools

cd(dirname(@__FILE__))
lines = open("three.txt", "r") do f
    lines::Array{Array{Int64,1},1} = []
    for l in eachline(f)
        line = [(c == '#' ?  1 : 0) for c in l]
        push!(lines, line)
    end
    lines
end

function make_slide(right::Int64, down::Int64, m::Int64)
    function slide((x, y)::Tuple{Int64, Int64})::Tuple{Int64, Int64}
        x̂ = (x + right - 1) % m + 1 # OK, here it is really annoying that indices in Julia start at 1 ...
        ŷ = y + down
        x̂, ŷ
    end
end

route(x::Int64, y::Int64, r::Int64, d::Int64, m::Int64, l::Int64) =
    takewhile((_,y)::Tuple{Int64, Int64} -> y <= l, iterated(make_slide(r, d, m), (x, y)))

function make_tree(lines::Array{Array{Int64,1},1})
    tree((x, y)::Tuple{Int64, Int64})::Bool = lines[y][x]
end

trees(x::Int64, y::Int64, right::Int64, down::Int64, lines::Array{Array{Int64,1},1}) =
    sum(map(make_tree(lines), route(x, y, right, down, length(lines[1]), length(lines))))

t1 = trees(1, 1, 1, 1, lines)
t2 = trees(1, 1, 3, 1, lines)
t3 = trees(1, 1, 5, 1, lines)
t4 = trees(1, 1, 7, 1, lines)
t5 = trees(1, 1, 1, 2, lines)

println(t2)
println(t1*t2*t3*t4*t5)


text:
#=
    Genetic algorithm is a search algorithm, which a computational
    metaphor on how natural evolution happen. It is useful to solve
    a wide range of problems that are hard to solve using brute-force.
=#

using Pkg
# Pkg.add("Test")
using Test
using Random

# To make everything reproducible
Random.seed!(42)

# Create a gene value
function createGene()
    return floor(Int, rand() * 10)
end

# Create an individual
function createIndividual()
    return [createGene() for i in 1:5]
end

# Create population
function createPopulation()
    return [createIndividual() for i in 1:10]
end

# Genetic operation: Crossover
function crossover(ind1, ind2)
    size = length(ind1)
    randomIndex = floor(Int64, rand() * size)
    return crossoverAtIndex(ind1, ind2, randomIndex)
end

function crossoverAtIndex(ind1, ind2, index)
    return vcat(ind1[1:index], ind2[index + 1:length(ind2)])
end

@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 2) == [1, 2, 30, 40]
@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 1) == [1, 20, 30, 40]
@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 0) == [10, 20, 30, 40]
@test crossoverAtIndex([1, 2, 3, 4], [10, 20, 30, 40], 4) == [1, 2, 3, 4]

# Genetic operation: Mutation
function mutate(ind)
    size = length(ind)
    randomIndex = floor(Int64, rand() * size + 1)
    return mutateAtIndex(ind, randomIndex)
end

function mutateAtIndex(ind, index)
    t = copy(ind)
    t[index] = createGene()
    return t
end

# Return true if the two individuals are the same, expect with a mutation at a particular index
function checkIfSimilar(ind1, ind2)
    pairs = collect(zip(ind1, ind2))
    # the predicate is the same than t -> first(t) == last(t)
    allIdentical = filter(((a, b),) -> a == b, pairs)
    allDifferent = filter(((a, b),) -> a != b, pairs)
    oneDiff = length(allDifferent) == 1
    allTheRestTheSame = length(allIdentical) == (length(ind1) - 1)
    return oneDiff && allTheRestTheSame
end
@test checkIfSimilar([1, 2, 3, 4], mutateAtIndex([1, 2, 3, 4], 2))
@test checkIfSimilar([1, 2, 3, 4], mutateAtIndex([1, 2, 3, 4], 1))
@test checkIfSimilar([1, 2, 3, 4], mutateAtIndex([1, 2, 3, 4], 4))

# Fitness function
function fitness(ind)
    solution = 1:5
    allDifferent = filter(((a,b),) -> a != b, collect(zip(solution, ind)))
    return length(allDifferent)
end

@test fitness([3, 2, 1, 4, 5]) == 2
@test fitness(reverse(collect(1:5))) == 4
@test fitness(1:5) == 0

function bestFitnessOf(population)
    return maximum(map(fitness, population))
end

# Pick best individual
function pickBestIndividual(population)
    indexMaxFitness = argmin(map(fitness, population))
    bestIndividual = population[indexMaxFitness]
    return bestIndividual
end

@test pickBestIndividual([[3, 2, 1, 4, 5], reverse(collect(1:5)), 1:5]) == 1:5
@test pickBestIndividual([[3, 2, 1, 4, 5], reverse(collect(1:5))]) == [3, 2, 1, 4, 5]

# Selection using the tournament of size k
# Assume that the smallest fitness is the best
function selectIndividual(population; k=5)
    populationSize = length(population)
    selectedIndividuals = []
    for tmp in 1:k
        anInd = population[floor(Int64, rand()*populationSize) + 1]
        push!(selectedIndividuals, anInd)
    end
    return pickBestIndividual(selectedIndividuals)
end

@test selectIndividual([[3, 2, 1, 4, 5], reverse(collect(1:5)), 1:5]) == 1:5
@test selectIndividual([[3, 2, 1, 4, 5], reverse(collect(1:5))]) == [3, 2, 1, 4, 5]

# We are now ready to run the algorithm
function run(;maxNumberOfIterations=10, probMutation=0.2)
    population = createPopulation()
    numberOfIndividuals = length(population)
    fitnesses = []
    for iteration in 1:maxNumberOfIterations
        newPopulation = []
        for it in 1:numberOfIndividuals
            ind1 = selectIndividual(population)
            ind2 = selectIndividual(population)
            newIndividual = crossover(ind1, ind2)
            if (rand() <= probMutation)
                newIndividual = mutate(newIndividual)
            end
            push!(newPopulation, newIndividual)
        end
        population = newPopulation
        push!(fitnesses, bestFitnessOf(population))
    end
    return pickBestIndividual(population), fitnesses
end

run(maxNumberOfIterations=40)


text:
using Artifacts

const gmsh_path = artifact"gmsh4.7.0"

const gmsh_dirs      = readdir(gmsh_path,join=true)

@assert length(gmsh_dirs)==1 "there should be only one directory for the untared gmsh file. Got $dirs"

include(joinpath(first(gmsh_dirs),"lib","gmsh.jl"))

export gmsh # gmsh module

"""
macro gmsh(ex)

Initialize `gmsh` through `gmsh.initilize(), execute `ex`, the close `gmsh`
through `gmsh.finalize()`.
"""
macro gmsh(ex)
    return quote
        gmsh.initialize()
        $(esc(ex))
        gmsh.finalize()
    end
end

"""
    read_geo(fname::String;dim=3)

Read a `.geo` file and generate a [`GenericMesh`](@ref) together with a
[`Domain`](@ref) for it using `gmsh`.
"""
function read_geo(fname;dim=3,h=nothing,order=nothing)
    assert_extension(fname, ".geo")
    gmsh.initialize()
    gmsh.open(fname)
    h     === nothing || _gmsh_set_meshsize(h)
    order === nothing || _gmsh_set_meshorder(order)
    gmsh.model.mesh.generate(dim)
    Ω = _initialize_domain(dim)
    M = _initialize_mesh(Ω)
    if dim == 2
        M = convert_to_2d(M)
    end
    gmsh.finalize()
    return Ω,M
end

"""
    read_msh(fname::String)

Similar to [`read_geo`](@ref), but the mesh is simply read from the input file
instead of generated.
"""
function read_msh(fname;dim=3)
    assert_extension(fname, ".msh")
    gmsh.initialize()
    gmsh.open(fname)
    Ω = _initialize_domain(dim)
    M = _initialize_mesh(Ω)
    gmsh.finalize()
    return Ω, M
end

"""
    _initialize_domain(d)

Construct a `Domain` from the current `gmsh` model, starting from entities of dimension `d`.

This is a helper function, and should not be called by itself since it assumes that `gmsh` has been initialized.
"""
function _initialize_domain(dim)
    Ω = Domain() # Create empty domain
    dim_tags = gmsh.model.getEntities(dim)
    for (_, tag) in dim_tags
        if haskey(ENTITIES,(dim,tag))
            ent = ENTITIES[(dim,tag)]
        else
            ent = ElementaryEntity(dim, tag)
            _fill_entity_boundary!(ent)
        end
        push!(Ω, ent)
    end
    return Ω
end

"""
    _fill_entity_boundary!

Use the `gmsh` API to add the boundary of an `ElementaryEntity`.

This is a helper function, and should not be called by itself since it assumes that `gmsh` has been initialized.
"""
function _fill_entity_boundary!(ent)
    combine  = true # FIXME: what should we use here?
    oriented = false
    dim_tags = gmsh.model.getBoundary((ent.dim, ent.tag),combine,oriented)
    for (d, t) in dim_tags
        if haskey(ENTITIES,(d,t))
            bnd = ENTITIES[(d,t)]
        else
            bnd = ElementaryEntity(d, t)
            _fill_entity_boundary!(bnd)
        end
        push!(ent.boundary, bnd)
    end
    return ent
end

"""
    _initialize_mesh(Ω::Domain)

Performs all the GMSH API calls to extract the information necessary to
construct the mesh.
"""
function _initialize_mesh(Ω::Domain)
    tags, coord, _ = gmsh.model.mesh.getNodes()
    nodes = reinterpret(SVector{3,Float64}, coord) |> collect
    # map gmsh type tags to actual internal types
    etypes = [_type_tag_to_etype(e) for e in gmsh.model.mesh.getElementTypes()]
    # Recursively populating the dictionaries
    elements = OrderedDict{DataType,Matrix{Int}}()
    ent2tags = OrderedDict{AbstractEntity,OrderedDict{DataType,Vector{Int}}}()
    elements, ent2tag = _domain_to_mesh!(elements, ent2tags, Ω)
    return GenericMesh{3,Float64}(;nodes, etypes, elements, ent2tags)
end

"""
    _domain_to_mesh!(elements, ent2tag, Ω::Domain)

Recursively populating the dictionaries `elements` and `ent2tag`.
"""
function _domain_to_mesh!(elements, ent2tag, Ω::Domain)
    isempty(Ω) && (return elements, ent2tag)
    for ω in Ω
        _ent_to_mesh!(elements, ent2tag, ω)
    end
    Γ = skeleton(Ω)
    _domain_to_mesh!(elements, ent2tag, Γ)
end

"""
    _ent_to_mesh!(elements, ent2tag, ω::ElementaryEntity)

For each element type used to mesh `ω`:
- push into `elements::OrderedDict` the pair `etype=>ntags`;
- push into `ent2tag::OrderedDict` the pair `etype=>etags`;

where:
- `etype::DataType` determines the type of the element (see
    [`_type_tag_to_etype`](@ref));
- `ntags::Matrix{Int}` gives the indices of the nodes defining those
    elements;
- `etags::Vector{Int}` gives the indices of those elements in `elements`.
"""
function _ent_to_mesh!(elements, ent2tag, ω::ElementaryEntity)
    ω in keys(ent2tag) && (return elements, ent2tag)
    etypes_to_etags = OrderedDict{DataType,Vector{Int}}()
    # Loop on GMSH element types (integer)
    type_tags, _, ntagss = gmsh.model.mesh.getElements(key(ω)...)
    for (type_tag, ntags) in zip(type_tags, ntagss)
        _, _, _, Np, _ = gmsh.model.mesh.getElementProperties(type_tag)
        ntags = reshape(ntags, Int(Np), :)
        etype = _type_tag_to_etype(type_tag)
        if etype in keys(elements)
            etag = size(elements[etype], 2) .+ collect(1:size(ntags,2))
            ntags = hcat(elements[etype], ntags)
        else
            etag = collect(1:size(ntags, 2))
        end
        push!(elements, etype => ntags)
        push!(etypes_to_etags, etype => etag)
    end
    push!(ent2tag, ω => etypes_to_etags)
    return elements, ent2tag
end

"""
    gmsh_disk(;rx=0.5,ry=0.5,center=(0,0,0)) -> Ω, M

Use `gmsh` API to generate a disk and return a [`Domain`](@ref) `Ω` and a
[`GenericMesh`](@ref) of the disk.
"""
function gmsh_disk(;rx=0.5,ry=0.5,center=(0.,0.,0.),dim=2,h=min(rx,ry)/10,order=1)
    gmsh.initialize()
    _gmsh_set_meshsize(h)
    _gmsh_set_meshorder(order)
    gmsh.model.occ.addDisk(center...,rx,ry)
    gmsh.model.occ.synchronize()
    gmsh.model.mesh.generate(dim)
    Ω = _initialize_domain(2)
    M = _initialize_mesh(Ω)
    M = convert_to_2d(M)
    gmsh.finalize()
    return Ω,M
end

"""
    gmsh_sphere(;radius=0.5,center=(0,0,0),dim=3,h=radius/10,order=1) -> Ω, M

Use `gmsh` API to generate a sphere and return `Ω::Domain` and `M::GenericMesh`.
Only entities of dimension `≤ dim` are meshed.
"""
function gmsh_sphere(;radius=0.5,center=(0., 0., 0.),dim=3,h=radius/10,order=1,recombine=false)
    gmsh.initialize()
    _gmsh_set_meshsize(h)
    _gmsh_set_meshorder(order)
    gmsh.model.occ.addSphere(center..., radius)
    gmsh.model.occ.synchronize()
    gmsh.model.mesh.generate(dim)
    recombine && gmsh.model.mesh.recombine()
    Ω = _initialize_domain(3)
    M = _initialize_mesh(Ω)
    gmsh.finalize()
    return Ω, M
end

"""
    gmsh_box(;origin=(0,0,0),widths=(0,0,0)) -> Ω, M

Use `gmsh` API to generate an axis aligned box. Return `Ω::Domain` and `M::GenericMesh`.
"""
function gmsh_box(;origin=(0., 0., 0.),widths=(1., 1., 1.),h=0.1)
    @gmsh begin
        _gmsh_set_meshsize(h)
        gmsh.model.occ.addBox(origin..., widths...)
        gmsh.model.occ.synchronize()
        gmsh.model.mesh.generate()
        Ω = _initialize_domain(3)
        M = _initialize_mesh(Ω)
    end
    return Ω, M
end

"""
    gmsh_rectangle(;origin,widths)
"""
function gmsh_rectangle(;origin=(0.,0.,0.),dx=1,dy=1,dim=2,h=0.1)
    @gmsh begin
        _gmsh_set_meshsize(h)
        gmsh.model.occ.addRectangle(origin...,dx,dy)
        gmsh.model.occ.synchronize()
        gmsh.model.mesh.generate(dim)
        Ω = _initialize_domain(dim)
        M = _initialize_mesh(Ω)
        if dim == 2
            M = convert_to_2d(M)
        end
    end
    return Ω,M
end

function _gmsh_set_meshsize(hmax, hmin=hmax)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", hmin)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", hmax)
end

function _gmsh_set_meshorder(order)
    gmsh.option.setNumber("Mesh.ElementOrder", order)
end

function _gmsh_set_verbosity(i)
    gmsh.option.setNumber("General.Verbosity",i)
end

function _gmsh_summary(model)
    gmsh.model.setCurrent(model)
    @printf("List of entities in model `%s`: \n", model)
    @printf("|%10s|%10s|%10s|\n","name","dimension","tag")
    ents = gmsh.model.getEntities()
    # pgroups = gmsh.model.getPhysicalGroups()
    for ent in ents
        name = gmsh.model.getEntityName(ent...)
        dim, tag = ent
        @printf("|%10s|%10d|%10d|\n", name, dim, tag)
    end
    println()
end

function _gmsh_summary()
    models = gmsh.model.list()
    for model in models
        _gmsh_summary(model)
    end
end

"""
    _type_tag_to_etype(tag)

Mapping of `gmsh` element types, encoded as an integer, to the internal
equivalent of those. This function assumes `gmsh` has been initilized.
"""
function _type_tag_to_etype(tag)
    T = Float64
    name,dim,order,num_nodes,ref_nodes,num_primary_nodes  = gmsh.model.mesh.getElementProperties(tag)
    num_nodes = Int(num_nodes) #convert to Int64
    if occursin("Point",name)
        etype = SVector{3,T}
    elseif occursin("Line",name)
    etype = LagrangeLine{num_nodes,SVector{3,T}}
    elseif occursin("Triangle",name)
        etype = LagrangeTriangle{num_nodes,SVector{3,T}}
    elseif occursin("Quadrilateral",name)
        etype = LagrangeRectangle{num_nodes,SVector{3,T}}
    elseif occursin("Tetrahedron",name)
        etype = LagrangeTetrahedron{num_nodes,SVector{3,T}}
    else
        error("gmsh element of family $name does not an internal equivalent")
    end
    return etype
end

"""
    _etype_to_type_tag(etype)

Mapping of internal element types, to the integer tag of `gmsh` elements. This
function assumes `gmsh` has been initialized.
"""
function _etype_to_type_tag(el::LagrangeElement)
    etype = typeof(el)
    tag = 1
    while true
        E   = _type_tag_to_etype(tag)
        E === etype && (return tag)
        tag = tag + 1
    end
end

# use gmsh API to extract the reference nodes
# TODO: we should just code these internally at some point instead of relying on
# gmsh.
@generated function reference_nodes(el::LagrangeElement{ReferenceLine})
    gmsh.initialize()
    p           = Geometry.order(el)
    family_name = "line"
    tag = gmsh.model.mesh.getElementType(family_name,p)
    name,dim,order,num_nodes,ref_nodes,num_primary_nodes  = gmsh.model.mesh.getElementProperties(tag)
    num_nodes = Int(num_nodes)
    gmsh.finalize()
    ref_nodes = SVector{num_nodes}(ref_nodes)
    return :($ref_nodes)
end

@generated function reference_nodes(el::LagrangeElement{ReferenceTriangle})
    gmsh.initialize()
    p           = Geometry.order(el)
    family_name = "triangle"
    tag = gmsh.model.mesh.getElementType(family_name,p)
    name,dim,order,num_nodes,ref_nodes,num_primary_nodes  = gmsh.model.mesh.getElementProperties(tag)
    num_nodes = Int(num_nodes)
    ref_nodes = reshape(ref_nodes,2,num_nodes)
    gmsh.finalize()
    sref_nodes = svector(num_nodes) do i
        SVector{2}(ref_nodes[:,i])
    end
    return :($sref_nodes)
end

"""
    struct GmshParametricEntity{M} <: AbstractEntity

An attempt to wrap a `gmsh` entity into an internal object for which a
parametric access is available. Like a `ParametricEntity`, the element is given
by the image of `el(x)` for `x ∈ domain(el)`, where `el::GmshParametricEntity`.

The interface implements a somewhat *hacky* way to use `gmsh` as a `CAD` reader.
Some important points to keep in mind (which makes this *hacky*):
- Each call `el(x)` entails calling a function on the `gmsh` library (linked
  dynamically). In particular the call signature for the `gmsh`'s `getValue`
  method requires passing a `Vector` for the input point, and therefore some
  unnecessary allocations are incurred. Overall, this means **this interface is inneficient**.
- Only untrimmed surfaces should be considered; thus **this interface is
  limited** to somewhat simple surfaces.
- **No guarantee is given on the quality of the parametrization**. For instance, a
  sphere may be parametrized using spherical coordinates, for which a geometric
  singularity exists at the poles. For complex surfaces, it is typically much
  better to simply mesh the surface using `gmsh` and then import the file.

Conclusion: although I expect this functionality to become more useful in the
future for e.g. doing isogeometric analysis, for the moment **it should be used
with extreme caution**.
"""
struct GmshParametricEntity{M} <: AbstractEntity
    tag::Int
    domain::HyperRectangle{M,Float64}
    function GmshParametricEntity{M}(tag,domain) where {M}
        dim = M
        ent = new{M}(dim,domain)
        _global_add_entity!(ent)
        return ent
    # TODO: throw a warning if the surface is trimmed
    end
end

function GmshParametricEntity(dim::Int,tag::Int,model=gmsh.model.getCurrent())
    low_corner,high_corner = gmsh.model.getParametrizationBounds(dim,tag)
    rec = HyperRectangle(low_corner,high_corner)
    return GmshParametricEntity{dim}(tag,rec)
end
GmshParametricEntity(dim::Integer,tag::Integer,args...;kwargs...) = GmshParametricEntity(Int(dim),Int(tag),args...;kwargs...)

function (par::GmshParametricEntity{N})(x) where {N}
    if N === 1
        return gmsh.model.getValue(N,par.tag,x)
    elseif N===2
        return gmsh.model.getValue(N,par.tag,[x[1],x[2]])
    else
        error("got N=$N, values must be 1 or 2")
    end
end

function jacobian(psurf::GmshParametricEntity{N},s::SVector) where {N}
    if N==1
        jac = gmsh.model.getDerivative(N,psurf.tag,s)
        return SMatrix{3,1}(jac)
    elseif N==2
        jac = gmsh.model.getDerivative(N,psurf.tag,[s[1],s[2]])
        return SMatrix{3,2}(jac)
    else
        error("got N=$N, values must be 1 or 2")
    end
end

# struct GmshParametricBody{M} <: AbstractParametricBody{3,M,Float64}
#     parts::Vector{GmshParametricEntity{M}}
# end

# function GmshParametricBody(dim,tag,model=gmsh.model.getCurrent())
#     dimtags = gmsh.model.getBoundary((dim,tag))
#     body    = GmshParametricBody{2}([])
#     for dimtag in dimtags
#         patch = GmshParametricEntity(Int(dimtag[1]),Int(dimtag[2]),model)
#         push!(body.parts,patch)
#     end
#     return body
# end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp