text:
@testset "6.5.7 (d hyper)^m (a+b (c sech)^n)^p" begin
    (a, b, c, d, x, ) = @variables a b c d x

    #= ::Package:: =#

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

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

    @test_int [(a + b*sech(c + d*x)^2)*sinh(c + d*x)^4, x, 5, (3/8)*(a - 4*b)*x - ((5*a - 4*b)*cosh(c + d*x)*sinh(c + d*x))/(8*d) + (a*cosh(c + d*x)^3*sinh(c + d*x))/(4*d) + (b*tanh(c + d*x))/d]
    @test_int [(a + b*sech(c + d*x)^2)*sinh(c + d*x)^3, x, 3, -(((a - b)*cosh(c + d*x))/d) + (a*cosh(c + d*x)^3)/(3*d) + (b*sech(c + d*x))/d]
    @test_int [(a + b*sech(c + d*x)^2)*sinh(c + d*x)^2, x, 4, (-(1/2))*(a - 2*b)*x + (a*cosh(c + d*x)*sinh(c + d*x))/(2*d) - (b*tanh(c + d*x))/d]
    @test_int [(a + b*sech(c + d*x)^2)*sinh(c + d*x), x, 3, (a*cosh(c + d*x))/d - (b*sech(c + d*x))/d]
    @test_int [csch(c + d*x)*(a + b*sech(c + d*x)^2), x, 3, -(((a + b)*atanh(cosh(c + d*x)))/d) + (b*sech(c + d*x))/d]
    @test_int [csch(c + d*x)^2*(a + b*sech(c + d*x)^2), x, 3, -(((a + b)*coth(c + d*x))/d) - (b*tanh(c + d*x))/d]
    @test_int [csch(c + d*x)^3*(a + b*sech(c + d*x)^2), x, 4, ((a + 3*b)*atanh(cosh(c + d*x)))/(2*d) - ((a + b)*coth(c + d*x)*csch(c + d*x))/(2*d) - (b*sech(c + d*x))/d]
    @test_int [csch(c + d*x)^4*(a + b*sech(c + d*x)^2), x, 3, ((a + 2*b)*coth(c + d*x))/d - ((a + b)*coth(c + d*x)^3)/(3*d) + (b*tanh(c + d*x))/d]

    @test_int [(a + b*sech(c + d*x)^2)^2*sinh(c + d*x)^4, x, 6, (1/8)*(3*a^2 - 24*a*b + 8*b^2)*x - (a*(a - 8*b)*cosh(c + d*x)*sinh(c + d*x))/(8*d) - ((a^2 - 8*a*b + 4*b^2)*tanh(c + d*x))/(4*d) + (a^2*sinh(c + d*x)^4*tanh(c + d*x))/(4*d) - (b^2*tanh(c + d*x)^3)/(3*d)]
    @test_int [(a + b*sech(c + d*x)^2)^2*sinh(c + d*x)^3, x, 3, -((a*(a - 2*b)*cosh(c + d*x))/d) + (a^2*cosh(c + d*x)^3)/(3*d) + ((2*a - b)*b*sech(c + d*x))/d + (b^2*sech(c + d*x)^3)/(3*d)]
    @test_int [(a + b*sech(c + d*x)^2)^2*sinh(c + d*x)^2, x, 5, (-(1/2))*a*(a - 4*b)*x + (a*(a - 4*b)*tanh(c + d*x))/(2*d) + (a^2*sinh(c + d*x)^2*tanh(c + d*x))/(2*d) + (b^2*tanh(c + d*x)^3)/(3*d)]
    @test_int [(a + b*sech(c + d*x)^2)^2*sinh(c + d*x), x, 3, (a^2*cosh(c + d*x))/d - (2*a*b*sech(c + d*x))/d - (b^2*sech(c + d*x)^3)/(3*d)]
    @test_int [csch(c + d*x)*(a + b*sech(c + d*x)^2)^2, x, 4, -(((a + b)^2*atanh(cosh(c + d*x)))/d) + (b*(2*a + b)*sech(c + d*x))/d + (b^2*sech(c + d*x)^3)/(3*d)]
    @test_int [csch(c + d*x)^2*(a + b*sech(c + d*x)^2)^2, x, 3, -(((a + b)^2*coth(c + d*x))/d) - (2*b*(a + b)*tanh(c + d*x))/d + (b^2*tanh(c + d*x)^3)/(3*d)]
    @test_int [csch(c + d*x)^3*(a + b*sech(c + d*x)^2)^2, x, 5, ((a + b)*(a + 5*b)*atanh(cosh(c + d*x)))/(2*d) - ((3*a^2 + 6*a*b + 5*b^2)*coth(c + d*x)*csch(c + d*x))/(6*d) - (b*(6*a + 5*b)*sech(c + d*x))/(3*d) + (b^2*csch(c + d*x)^2*sech(c + d*x)^3)/(3*d)]
    @test_int [csch(c + d*x)^4*(a + b*sech(c + d*x)^2)^2, x, 3, ((a + b)*(a + 3*b)*coth(c + d*x))/d - ((a + b)^2*coth(c + d*x)^3)/(3*d) + (b*(2*a + 3*b)*tanh(c + d*x))/d - (b^2*tanh(c + d*x)^3)/(3*d)]

    @test_int [(a + b*sech(c + d*x)^2)^3*sinh(c + d*x)^4, x, 6, (3/8)*a*(a^2 - 12*a*b + 8*b^2)*x - (3*a*(a^2 - 12*a*b + 8*b^2)*tanh(c + d*x))/(8*d) + (b*(6*a^2 - 23*a*b - 8*b^2)*tanh(c + d*x)^3)/(8*d) - (3*(5*a - 16*b)*b^2*tanh(c + d*x)^5)/(40*d) - (3*(a - 2*b)*sinh(c + d*x)^2*tanh(c + d*x)*(a + b - b*tanh(c + d*x)^2)^2)/(8*d) + (cosh(c + d*x)*sinh(c + d*x)^3*(a + b - b*tanh(c + d*x)^2)^3)/(4*d)]
    @test_int [(a + b*sech(c + d*x)^2)^3*sinh(c + d*x)^3, x, 3, -((a^2*(a - 3*b)*cosh(c + d*x))/d) + (a^3*cosh(c + d*x)^3)/(3*d) + (3*a*(a - b)*b*sech(c + d*x))/d + ((3*a - b)*b^2*sech(c + d*x)^3)/(3*d) + (b^3*sech(c + d*x)^5)/(5*d)]
    @test_int [(a + b*sech(c + d*x)^2)^3*sinh(c + d*x)^2, x, 6, -(a^2*(a - 6*b)*x)/2 + a^3/(4*d*(1 - tanh(c + d*x))) - (3*a^2*b*tanh(c + d*x))/d + (b^2*(3*a + b)*tanh(c + d*x)^3)/(3*d) - (b^3*tanh(c + d*x)^5)/(5*d) - a^3/(4*d*(1 + tanh(c + d*x))), (-(1/2))*a^2*(a - 6*b)*x - (b*(81*a^2 - 28*a*b - 4*b^2)*tanh(c + d*x))/(30*d) - ((33*a - 2*b)*b*tanh(c + d*x)*(a + b - b*tanh(c + d*x)^2))/(30*d) - (7*b*tanh(c + d*x)*(a + b - b*tanh(c + d*x)^2)^2)/(10*d) + (cosh(c + d*x)*sinh(c + d*x)*(a + b - b*tanh(c + d*x)^2)^3)/(2*d)]
    @test_int [(a + b*sech(c + d*x)^2)^3*sinh(c + d*x), x, 3, (a^3*cosh(c + d*x))/d - (3*a^2*b*sech(c + d*x))/d - (a*b^2*sech(c + d*x)^3)/d - (b^3*sech(c + d*x)^5)/(5*d)]
    @test_int [csch(c + d*x)*(a + b*sech(c + d*x)^2)^3, x, 4, -(((a + b)^3*atanh(cosh(c + d*x)))/d) + (b*(3*a^2 + 3*a*b + b^2)*sech(c + d*x))/d + (b^2*(3*a + b)*sech(c + d*x)^3)/(3*d) + (b^3*sech(c + d*x)^5)/(5*d)]
    @test_int [csch(c + d*x)^2*(a + b*sech(c + d*x)^2)^3, x, 3, -(((a + b)^3*coth(c + d*x))/d) - (3*b*(a + b)^2*tanh(c + d*x))/d + (b^2*(a + b)*tanh(c + d*x)^3)/d - (b^3*tanh(c + d*x)^5)/(5*d)]
    @test_int [csch(c + d*x)^3*(a + b*sech(c + d*x)^2)^3, x, 5, ((a + b)^2*(a + 7*b)*atanh(cosh(c + d*x)))/(2*d) - ((a + b)^2*(a + 7*b)*sech(c + d*x))/(2*d) - (b*(6*a^2 + 15*a*b + 7*b^2)*sech(c + d*x)^3)/(6*d) - (b^2*(5*a + 7*b)*sech(c + d*x)^5)/(10*d) - ((a + b)*(b + a*cosh(c + d*x)^2)^2*csch(c + d*x)^2*sech(c + d*x)^5)/(2*d)]
    @test_int [csch(c + d*x)^4*(a + b*sech(c + d*x)^2)^3, x, 3, ((a + b)^2*(a + 4*b)*coth(c + d*x))/d - ((a + b)^3*coth(c + d*x)^3)/(3*d) + (3*b*(a + b)*(a + 2*b)*tanh(c + d*x))/d - (b^2*(3*a + 4*b)*tanh(c + d*x)^3)/(3*d) + (b^3*tanh(c + d*x)^5)/(5*d)]

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

    @test_int [sinh(c + d*x)^4/(a + b*sech(c + d*x)^2), x, 6, ((3*a^2 + 12*a*b + 8*b^2)*x)/(8*a^3) - (sqrt(b)*(a + b)^(3/2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a^3*d) - ((5*a + 4*b)*cosh(c + d*x)*sinh(c + d*x))/(8*a^2*d) + (cosh(c + d*x)^3*sinh(c + d*x))/(4*a*d)]
    @test_int [sinh(c + d*x)^3/(a + b*sech(c + d*x)^2), x, 4, (sqrt(b)*(a + b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(a^(5/2)*d) - ((a + b)*cosh(c + d*x))/(a^2*d) + cosh(c + d*x)^3/(3*a*d)]
    @test_int [sinh(c + d*x)^2/(a + b*sech(c + d*x)^2), x, 5, -(((a + 2*b)*x)/(2*a^2)) + (sqrt(b)*sqrt(a + b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a^2*d) + (cosh(c + d*x)*sinh(c + d*x))/(2*a*d)]
    @test_int [sinh(c + d*x)/(a + b*sech(c + d*x)^2), x, 3, -((sqrt(b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(a^(3/2)*d)) + cosh(c + d*x)/(a*d)]
    @test_int [csch(c + d*x)/(a + b*sech(c + d*x)^2), x, 4, (sqrt(b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(sqrt(a)*(a + b)*d) - atanh(cosh(c + d*x))/((a + b)*d)]
    @test_int [csch(c + d*x)^2/(a + b*sech(c + d*x)^2), x, 3, (sqrt(b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/((a + b)^(3/2)*d) - coth(c + d*x)/((a + b)*d)]
    @test_int [csch(c + d*x)^3/(a + b*sech(c + d*x)^2), x, 5, -((sqrt(a)*sqrt(b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/((a + b)^2*d)) + ((a - b)*atanh(cosh(c + d*x)))/(2*(a + b)^2*d) - (coth(c + d*x)*csch(c + d*x))/(2*(a + b)*d)]
    @test_int [csch(c + d*x)^4/(a + b*sech(c + d*x)^2), x, 4, -((a*sqrt(b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/((a + b)^(5/2)*d)) + (a*coth(c + d*x))/((a + b)^2*d) - coth(c + d*x)^3/(3*(a + b)*d)]

    @test_int [sinh(c + d*x)^4/(a + b*sech(c + d*x)^2)^2, x, 7, (3*(a^2 + 8*a*b + 8*b^2)*x)/(8*a^4) - (3*sqrt(b)*sqrt(a + b)*(a + 2*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^4*d) - ((5*a + 6*b)*cosh(c + d*x)*sinh(c + d*x))/(8*a^2*d*(a + b - b*tanh(c + d*x)^2)) + (cosh(c + d*x)^3*sinh(c + d*x))/(4*a*d*(a + b - b*tanh(c + d*x)^2)) - (3*b*(3*a + 4*b)*tanh(c + d*x))/(8*a^3*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sinh(c + d*x)^3/(a + b*sech(c + d*x)^2)^2, x, 5, (sqrt(b)*(3*a + 5*b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(2*a^(7/2)*d) - ((a + 2*b)*cosh(c + d*x))/(a^3*d) + cosh(c + d*x)^3/(3*a^2*d) - (b*(a + b)*cosh(c + d*x))/(2*a^3*d*(b + a*cosh(c + d*x)^2))]
    @test_int [sinh(c + d*x)^2/(a + b*sech(c + d*x)^2)^2, x, 6, -(((a + 4*b)*x)/(2*a^3)) + (sqrt(b)*(3*a + 4*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^3*sqrt(a + b)*d) + (cosh(c + d*x)*sinh(c + d*x))/(2*a*d*(a + b - b*tanh(c + d*x)^2)) + (b*tanh(c + d*x))/(a^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sinh(c + d*x)/(a + b*sech(c + d*x)^2)^2, x, 4, (-3*sqrt(b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(2*a^(5/2)*d) + (3*cosh(c + d*x))/(2*a^2*d) - cosh(c + d*x)^3/(2*a*d*(b + a*cosh(c + d*x)^2))]
    @test_int [csch(c + d*x)/(a + b*sech(c + d*x)^2)^2, x, 5, (sqrt(b)*(3*a + b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(2*a^(3/2)*(a + b)^2*d) - atanh(cosh(c + d*x))/((a + b)^2*d) - (b*cosh(c + d*x))/(2*a*(a + b)*d*(b + a*cosh(c + d*x)^2))]
    @test_int [csch(c + d*x)^2/(a + b*sech(c + d*x)^2)^2, x, 4, (3*sqrt(b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*(a + b)^(5/2)*d) - (3*coth(c + d*x))/(2*(a + b)^2*d) + coth(c + d*x)/(2*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [csch(c + d*x)^3/(a + b*sech(c + d*x)^2)^2, x, 6, -(((3*a - b)*sqrt(b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(2*sqrt(a)*(a + b)^3*d)) + ((a - 3*b)*atanh(cosh(c + d*x)))/(2*(a + b)^3*d) - ((a - b)*cosh(c + d*x))/(2*(a + b)^2*d*(b + a*cosh(c + d*x)^2)) - (coth(c + d*x)*csch(c + d*x))/(2*(a + b)*d*(b + a*cosh(c + d*x)^2))]
    @test_int [csch(c + d*x)^4/(a + b*sech(c + d*x)^2)^2, x, 5, -(((3*a - 2*b)*sqrt(b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*(a + b)^(7/2)*d)) + ((a - b)*coth(c + d*x))/((a + b)^3*d) - coth(c + d*x)^3/(3*(a + b)^2*d) - (a*b*tanh(c + d*x))/(2*(a + b)^3*d*(a + b - b*tanh(c + d*x)^2))]

    @test_int [sinh(c + d*x)^4/(a + b*sech(c + d*x)^2)^3, x, 8, (3*(a^2 + 12*a*b + 16*b^2)*x)/(8*a^5) - (3*sqrt(b)*(5*a^2 + 20*a*b + 16*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^5*sqrt(a + b)*d) - ((5*a + 8*b)*cosh(c + d*x)*sinh(c + d*x))/(8*a^2*d*(a + b - b*tanh(c + d*x)^2)^2) + (cosh(c + d*x)^3*sinh(c + d*x))/(4*a*d*(a + b - b*tanh(c + d*x)^2)^2) - (b*(7*a + 12*b)*tanh(c + d*x))/(8*a^3*d*(a + b - b*tanh(c + d*x)^2)^2) - (3*b*(a + 2*b)*tanh(c + d*x))/(2*a^4*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sinh(c + d*x)^3/(a + b*sech(c + d*x)^2)^3, x, 6, (5*sqrt(b)*(3*a + 7*b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(8*a^(9/2)*d) - ((a + 3*b)*cosh(c + d*x))/(a^4*d) + cosh(c + d*x)^3/(3*a^3*d) + (b^2*(a + b)*cosh(c + d*x))/(4*a^4*d*(b + a*cosh(c + d*x)^2)^2) - (b*(9*a + 13*b)*cosh(c + d*x))/(8*a^4*d*(b + a*cosh(c + d*x)^2))]
    @test_int [sinh(c + d*x)^2/(a + b*sech(c + d*x)^2)^3, x, 7, -(((a + 6*b)*x)/(2*a^4)) + (sqrt(b)*(15*a^2 + 40*a*b + 24*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^4*(a + b)^(3/2)*d) + (cosh(c + d*x)*sinh(c + d*x))/(2*a*d*(a + b - b*tanh(c + d*x)^2)^2) + (3*b*tanh(c + d*x))/(4*a^2*d*(a + b - b*tanh(c + d*x)^2)^2) + (b*(11*a + 12*b)*tanh(c + d*x))/(8*a^3*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sinh(c + d*x)/(a + b*sech(c + d*x)^2)^3, x, 5, (-15*sqrt(b)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(8*a^(7/2)*d) + (15*cosh(c + d*x))/(8*a^3*d) - cosh(c + d*x)^5/(4*a*d*(b + a*cosh(c + d*x)^2)^2) - (5*cosh(c + d*x)^3)/(8*a^2*d*(b + a*cosh(c + d*x)^2))]
    @test_int [csch(c + d*x)/(a + b*sech(c + d*x)^2)^3, x, 6, (sqrt(b)*(15*a^2 + 10*a*b + 3*b^2)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(8*a^(5/2)*(a + b)^3*d) - atanh(cosh(c + d*x))/((a + b)^3*d) - (b*cosh(c + d*x)^3)/(4*a*(a + b)*d*(b + a*cosh(c + d*x)^2)^2) - (b*(7*a + 3*b)*cosh(c + d*x))/(8*a^2*(a + b)^2*d*(b + a*cosh(c + d*x)^2))]
    @test_int [csch(c + d*x)^2/(a + b*sech(c + d*x)^2)^3, x, 5, (15*sqrt(b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*(a + b)^(7/2)*d) - (15*coth(c + d*x))/(8*(a + b)^3*d) + coth(c + d*x)/(4*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) + (5*coth(c + d*x))/(8*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [csch(c + d*x)^3/(a + b*sech(c + d*x)^2)^3, x, 7, -((sqrt(b)*(15*a^2 - 10*a*b - b^2)*atan((sqrt(a)*cosh(c + d*x))/sqrt(b)))/(8*a^(3/2)*(a + b)^4*d)) + ((a - 5*b)*atanh(cosh(c + d*x)))/(2*(a + b)^4*d) + ((2*a - b)*b*cosh(c + d*x))/(4*a*(a + b)^2*d*(b + a*cosh(c + d*x)^2)^2) - ((4*a^2 - 9*a*b - b^2)*cosh(c + d*x))/(8*a*(a + b)^3*d*(b + a*cosh(c + d*x)^2)) - (cosh(c + d*x)*coth(c + d*x)^2)/(2*(a + b)*d*(b + a*cosh(c + d*x)^2)^2)]
    @test_int [csch(c + d*x)^4/(a + b*sech(c + d*x)^2)^3, x, 6, -((5*(3*a - 4*b)*sqrt(b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*(a + b)^(9/2)*d)) + ((a - 2*b)*coth(c + d*x))/((a + b)^4*d) - coth(c + d*x)^3/(3*(a + b)^3*d) - (a*b*tanh(c + d*x))/(4*(a + b)^3*d*(a + b - b*tanh(c + d*x)^2)^2) - ((7*a - 4*b)*b*tanh(c + d*x))/(8*(a + b)^4*d*(a + b - b*tanh(c + d*x)^2))]

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

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

    @test_int [cosh(c + d*x)^4*(a + b*sech(c + d*x)^2), x, 3, ((3*a + 4*b)*x)/8 + ((3*a + 4*b)*cosh(c + d*x)*sinh(c + d*x))/(8*d) + (a*cosh(c + d*x)^3*sinh(c + d*x))/(4*d)]
    @test_int [cosh(c + d*x)^3*(a + b*sech(c + d*x)^2), x, 3, ((a + b)*sinh(c + d*x))/d + (a*sinh(c + d*x)^3)/(3*d)]
    @test_int [cosh(c + d*x)^2*(a + b*sech(c + d*x)^2), x, 2, ((a + 2*b)*x)/2 + (a*cosh(c + d*x)*sinh(c + d*x))/(2*d)]
    @test_int [cosh(c + d*x)^1*(a + b*sech(c + d*x)^2), x, 2, (b*atan(sinh(c + d*x)))/d + (a*sinh(c + d*x))/d]
    @test_int [sech(c + d*x)^1*(a + b*sech(c + d*x)^2), x, 2, ((2*a + b)*atan(sinh(c + d*x)))/(2*d) + (b*sech(c + d*x)*tanh(c + d*x))/(2*d)]
    @test_int [sech(c + d*x)^2*(a + b*sech(c + d*x)^2), x, 3, ((a + b)*tanh(c + d*x))/d - (b*tanh(c + d*x)^3)/(3*d), ((3*a + 2*b)*tanh(c + d*x))/(3*d) + (b*sech(c + d*x)^2*tanh(c + d*x))/(3*d)]
    @test_int [sech(c + d*x)^3*(a + b*sech(c + d*x)^2), x, 3, ((4*a + 3*b)*atan(sinh(c + d*x)))/(8*d) + ((4*a + 3*b)*sech(c + d*x)*tanh(c + d*x))/(8*d) + (b*sech(c + d*x)^3*tanh(c + d*x))/(4*d)]
    @test_int [sech(c + d*x)^4*(a + b*sech(c + d*x)^2), x, 3, ((a + b)*tanh(c + d*x))/d - ((a + 2*b)*tanh(c + d*x)^3)/(3*d) + (b*tanh(c + d*x)^5)/(5*d), ((5*a + 4*b)*tanh(c + d*x))/(5*d) + (b*sech(c + d*x)^4*tanh(c + d*x))/(5*d) - ((5*a + 4*b)*tanh(c + d*x)^3)/(15*d)]

    @test_int [cosh(c + d*x)^4*(a + b*sech(c + d*x)^2)^2, x, 4, (1/8)*(3*a^2 + 8*a*b + 8*b^2)*x + (3*a*(a + 2*b)*cosh(c + d*x)*sinh(c + d*x))/(8*d) + (a*cosh(c + d*x)^3*sinh(c + d*x)*(a + b - b*tanh(c + d*x)^2))/(4*d)]
    @test_int [cosh(c + d*x)^3*(a + b*sech(c + d*x)^2)^2, x, 4, (b^2*atan(sinh(c + d*x)))/d + (a*(a + 2*b)*sinh(c + d*x))/d + (a^2*sinh(c + d*x)^3)/(3*d)]
    @test_int [cosh(c + d*x)^2*(a + b*sech(c + d*x)^2)^2, x, 5, (1/2)*a*(a + 4*b)*x + (a^2*cosh(c + d*x)*sinh(c + d*x))/(2*d) + (b^2*tanh(c + d*x))/d]
    @test_int [cosh(c + d*x)^1*(a + b*sech(c + d*x)^2)^2, x, 5, (b*(4*a + b)*atan(sinh(c + d*x)))/(2*d) + (a^2*sinh(c + d*x))/d + (b^2*sech(c + d*x)*tanh(c + d*x))/(2*d)]
    @test_int [sech(c + d*x)^1*(a + b*sech(c + d*x)^2)^2, x, 4, ((8*a^2 + 8*a*b + 3*b^2)*atan(sinh(c + d*x)))/(8*d) + (3*b*(2*a + b)*sech(c + d*x)*tanh(c + d*x))/(8*d) + (b*sech(c + d*x)^3*(a + b + a*sinh(c + d*x)^2)*tanh(c + d*x))/(4*d)]
    @test_int [sech(c + d*x)^2*(a + b*sech(c + d*x)^2)^2, x, 3, ((a + b)^2*tanh(c + d*x))/d - (2*b*(a + b)*tanh(c + d*x)^3)/(3*d) + (b^2*tanh(c + d*x)^5)/(5*d)]
    @test_int [sech(c + d*x)^3*(a + b*sech(c + d*x)^2)^2, x, 5, ((8*a^2 + 12*a*b + 5*b^2)*atan(sinh(c + d*x)))/(16*d) + ((8*a^2 + 12*a*b + 5*b^2)*sech(c + d*x)*tanh(c + d*x))/(16*d) + (b*(8*a + 5*b)*sech(c + d*x)^3*tanh(c + d*x))/(24*d) + (b*sech(c + d*x)^5*(a + b + a*sinh(c + d*x)^2)*tanh(c + d*x))/(6*d)]
    @test_int [sech(c + d*x)^4*(a + b*sech(c + d*x)^2)^2, x, 3, ((a + b)^2*tanh(c + d*x))/d - ((a + b)*(a + 3*b)*tanh(c + d*x)^3)/(3*d) + (b*(2*a + 3*b)*tanh(c + d*x)^5)/(5*d) - (b^2*tanh(c + d*x)^7)/(7*d)]

    @test_int [cosh(c + d*x)^4*(a + b*sech(c + d*x)^2)^3, x, 6, (3/8)*a*(a^2 + 4*a*b + 8*b^2)*x + (3*a^2*(a + 4*b)*cosh(c + d*x)*sinh(c + d*x))/(8*d) + (a^3*cosh(c + d*x)^3*sinh(c + d*x))/(4*d) + (b^3*tanh(c + d*x))/d]
    @test_int [cosh(c + d*x)^3*(a + b*sech(c + d*x)^2)^3, x, 5, (b^2*(6*a + b)*atan(sinh(c + d*x)))/(2*d) + (a^2*(a + 3*b)*sinh(c + d*x))/d + (a^3*sinh(c + d*x)^3)/(3*d) + (b^3*sech(c + d*x)*tanh(c + d*x))/(2*d)]
    @test_int [cosh(c + d*x)^2*(a + b*sech(c + d*x)^2)^3, x, 5, (1/2)*a^2*(a + 6*b)*x + (a^3*cosh(c + d*x)*sinh(c + d*x))/(2*d) + (b^2*(3*a + b)*tanh(c + d*x))/d - (b^3*tanh(c + d*x)^3)/(3*d)]
    @test_int [cosh(c + d*x)^1*(a + b*sech(c + d*x)^2)^3, x, 6, (3*b*(8*a^2 + 4*a*b + b^2)*atan(sinh(c + d*x)))/(8*d) + (a^3*sinh(c + d*x))/d + (3*b^2*(4*a + b)*sech(c + d*x)*tanh(c + d*x))/(8*d) + (b^3*sech(c + d*x)^3*tanh(c + d*x))/(4*d)]
    @test_int [sech(c + d*x)^1*(a + b*sech(c + d*x)^2)^3, x, 5, ((2*a + b)*(8*a^2 + 8*a*b + 5*b^2)*atan(sinh(c + d*x)))/(16*d) + (b*(44*a^2 + 44*a*b + 15*b^2)*sech(c + d*x)*tanh(c + d*x))/(48*d) + (5*b*(2*a + b)*sech(c + d*x)^3*(a + b + a*sinh(c + d*x)^2)*tanh(c + d*x))/(24*d) + (b*sech(c + d*x)^5*(a + b + a*sinh(c + d*x)^2)^2*tanh(c + d*x))/(6*d)]
    @test_int [sech(c + d*x)^2*(a + b*sech(c + d*x)^2)^3, x, 3, ((a + b)^3*tanh(c + d*x))/d - (b*(a + b)^2*tanh(c + d*x)^3)/d + (3*b^2*(a + b)*tanh(c + d*x)^5)/(5*d) - (b^3*tanh(c + d*x)^7)/(7*d)]
    @test_int [sech(c + d*x)^3*(a + b*sech(c + d*x)^2)^3, x, 6, ((64*a^3 + 144*a^2*b + 120*a*b^2 + 35*b^3)*atan(sinh(c + d*x)))/(128*d) + ((64*a^3 + 144*a^2*b + 120*a*b^2 + 35*b^3)*sech(c + d*x)*tanh(c + d*x))/(128*d) + (b*(72*a^2 + 92*a*b + 35*b^2)*sech(c + d*x)^3*tanh(c + d*x))/(192*d) + (b*(12*a + 7*b)*sech(c + d*x)^5*(a + b + a*sinh(c + d*x)^2)*tanh(c + d*x))/(48*d) + (b*sech(c + d*x)^7*(a + b + a*sinh(c + d*x)^2)^2*tanh(c + d*x))/(8*d)]
    @test_int [sech(c + d*x)^4*(a + b*sech(c + d*x)^2)^3, x, 3, ((a + b)^3*tanh(c + d*x))/d - ((a + b)^2*(a + 4*b)*tanh(c + d*x)^3)/(3*d) + (3*b*(a + b)*(a + 2*b)*tanh(c + d*x)^5)/(5*d) - (b^2*(3*a + 4*b)*tanh(c + d*x)^7)/(7*d) + (b^3*tanh(c + d*x)^9)/(9*d)]

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

    @test_int [cosh(c + d*x)^4/(a + b*sech(c + d*x)^2), x, 6, ((3*a^2 - 4*a*b + 8*b^2)*x)/(8*a^3) - (b^(5/2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a^3*sqrt(a + b)*d) + ((3*a - 4*b)*cosh(c + d*x)*sinh(c + d*x))/(8*a^2*d) + (cosh(c + d*x)^3*sinh(c + d*x))/(4*a*d)]
    @test_int [cosh(c + d*x)^3/(a + b*sech(c + d*x)^2), x, 4, (b^2*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(a^(5/2)*sqrt(a + b)*d) + ((a - b)*sinh(c + d*x))/(a^2*d) + sinh(c + d*x)^3/(3*a*d)]
    @test_int [cosh(c + d*x)^2/(a + b*sech(c + d*x)^2), x, 5, ((a - 2*b)*x)/(2*a^2) + (b^(3/2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a^2*sqrt(a + b)*d) + (cosh(c + d*x)*sinh(c + d*x))/(2*a*d)]
    @test_int [cosh(c + d*x)/(a + b*sech(c + d*x)^2), x, 3, -((b*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(a^(3/2)*sqrt(a + b)*d)) + sinh(c + d*x)/(a*d)]
    @test_int [sech(c + d*x)/(a + b*sech(c + d*x)^2), x, 2, atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b))/(sqrt(a)*sqrt(a + b)*d)]
    @test_int [sech(c + d*x)^2/(a + b*sech(c + d*x)^2), x, 2, atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b))/(sqrt(b)*sqrt(a + b)*d)]
    @test_int [sech(c + d*x)^3/(a + b*sech(c + d*x)^2), x, 4, atan(sinh(c + d*x))/(b*d) - (sqrt(a)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(b*sqrt(a + b)*d)]
    @test_int [sech(c + d*x)^4/(a + b*sech(c + d*x)^2), x, 3, -((a*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(b^(3/2)*sqrt(a + b)*d)) + tanh(c + d*x)/(b*d)]
    @test_int [sech(c + d*x)^5/(a + b*sech(c + d*x)^2), x, 5, -(((2*a - b)*atan(sinh(c + d*x)))/(2*b^2*d)) + (a^(3/2)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(b^2*sqrt(a + b)*d) + (sech(c + d*x)*tanh(c + d*x))/(2*b*d)]
    @test_int [sech(c + d*x)^6/(a + b*sech(c + d*x)^2), x, 4, (a^2*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(b^(5/2)*sqrt(a + b)*d) - ((a - b)*tanh(c + d*x))/(b^2*d) - tanh(c + d*x)^3/(3*b*d)]

    @test_int [cosh(c + d*x)^3/(a + b*sech(c + d*x)^2)^2, x, 5, (b^2*(6*a + 5*b)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(2*a^(7/2)*(a + b)^(3/2)*d) + ((a - 2*b)*sinh(c + d*x))/(a^3*d) + sinh(c + d*x)^3/(3*a^2*d) - (b^3*sinh(c + d*x))/(2*a^3*(a + b)*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [cosh(c + d*x)^2/(a + b*sech(c + d*x)^2)^2, x, 6, ((a - 4*b)*x)/(2*a^3) + (b^(3/2)*(5*a + 4*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^3*(a + b)^(3/2)*d) + (cosh(c + d*x)*sinh(c + d*x))/(2*a*d*(a + b - b*tanh(c + d*x)^2)) + (b*(a + 2*b)*tanh(c + d*x))/(2*a^2*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [cosh(c + d*x)/(a + b*sech(c + d*x)^2)^2, x, 5, -((b*(4*a + 3*b)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(2*a^(5/2)*(a + b)^(3/2)*d)) + sinh(c + d*x)/(a^2*d) + (b^2*sinh(c + d*x))/(2*a^2*(a + b)*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)/(a + b*sech(c + d*x)^2)^2, x, 3, ((2*a + b)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(2*a^(3/2)*(a + b)^(3/2)*d) - (b*sinh(c + d*x))/(2*a*(a + b)*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)^2/(a + b*sech(c + d*x)^2)^2, x, 3, atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b))/(2*sqrt(b)*(a + b)^(3/2)*d) + tanh(c + d*x)/(2*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sech(c + d*x)^3/(a + b*sech(c + d*x)^2)^2, x, 3, atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b))/(2*sqrt(a)*(a + b)^(3/2)*d) + sinh(c + d*x)/(2*(a + b)*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)^4/(a + b*sech(c + d*x)^2)^2, x, 3, ((a + 2*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*b^(3/2)*(a + b)^(3/2)*d) - (a*tanh(c + d*x))/(2*b*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sech(c + d*x)^5/(a + b*sech(c + d*x)^2)^2, x, 5, atan(sinh(c + d*x))/(b^2*d) - (sqrt(a)*(2*a + 3*b)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(2*b^2*(a + b)^(3/2)*d) - (a*sinh(c + d*x))/(2*b*(a + b)*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)^6/(a + b*sech(c + d*x)^2)^2, x, 5, -((a*(3*a + 4*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*b^(5/2)*(a + b)^(3/2)*d)) + tanh(c + d*x)/(b^2*d) + (a^2*tanh(c + d*x))/(2*b^2*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sech(c + d*x)^7/(a + b*sech(c + d*x)^2)^2, x, 6, -(((4*a - b)*atan(sinh(c + d*x)))/(2*b^3*d)) + (a^(3/2)*(4*a + 5*b)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(2*b^3*(a + b)^(3/2)*d) + (a*(2*a + b)*sinh(c + d*x))/(2*b^2*(a + b)*d*(a + b + a*sinh(c + d*x)^2)) + (sech(c + d*x)*tanh(c + d*x))/(2*b*d*(a + b + a*sinh(c + d*x)^2))]

    @test_int [cosh(c + d*x)^2/(a + b*sech(c + d*x)^2)^3, x, 7, ((a - 6*b)*x)/(2*a^4) + (b^(3/2)*(35*a^2 + 56*a*b + 24*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^4*(a + b)^(5/2)*d) + (cosh(c + d*x)*sinh(c + d*x))/(2*a*d*(a + b - b*tanh(c + d*x)^2)^2) + (b*(2*a + 3*b)*tanh(c + d*x))/(4*a^2*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) + (b*(4*a + 3*b)*(a + 4*b)*tanh(c + d*x))/(8*a^3*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [cosh(c + d*x)^1/(a + b*sech(c + d*x)^2)^3, x, 6, -((3*b*(4*(a + b)^2 + (2*a + b)^2)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(8*a^(7/2)*(a + b)^(5/2)*d)) + sinh(c + d*x)/(a^3*d) - (b^3*sinh(c + d*x))/(4*a^3*(a + b)*d*(a + b + a*sinh(c + d*x)^2)^2) + (3*b^2*(4*a + 3*b)*sinh(c + d*x))/(8*a^3*(a + b)^2*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)^1/(a + b*sech(c + d*x)^2)^3, x, 4, ((8*a^2 + 8*a*b + 3*b^2)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(8*a^(5/2)*(a + b)^(5/2)*d) - (b*cosh(c + d*x)^2*sinh(c + d*x))/(4*a*(a + b)*d*(a + b + a*sinh(c + d*x)^2)^2) - (3*b*(2*a + b)*sinh(c + d*x))/(8*a^2*(a + b)^2*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)^2/(a + b*sech(c + d*x)^2)^3, x, 4, (3*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*sqrt(b)*(a + b)^(5/2)*d) + tanh(c + d*x)/(4*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) + (3*tanh(c + d*x))/(8*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sech(c + d*x)^3/(a + b*sech(c + d*x)^2)^3, x, 4, ((4*a + b)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(8*a^(3/2)*(a + b)^(5/2)*d) - (b*sinh(c + d*x))/(4*a*(a + b)*d*(a + b + a*sinh(c + d*x)^2)^2) + ((4*a + b)*sinh(c + d*x))/(8*a*(a + b)^2*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)^4/(a + b*sech(c + d*x)^2)^3, x, 4, ((a + 4*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*b^(3/2)*(a + b)^(5/2)*d) - (a*tanh(c + d*x))/(4*b*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) + ((a + 4*b)*tanh(c + d*x))/(8*b*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sech(c + d*x)^5/(a + b*sech(c + d*x)^2)^3, x, 4, (3*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(8*sqrt(a)*(a + b)^(5/2)*d) + sinh(c + d*x)/(4*(a + b)*d*(a + b + a*sinh(c + d*x)^2)^2) + (3*sinh(c + d*x))/(8*(a + b)^2*d*(a + b + a*sinh(c + d*x)^2))]
    @test_int [sech(c + d*x)^6/(a + b*sech(c + d*x)^2)^3, x, 4, ((3*a^2 + 8*a*b + 8*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*b^(5/2)*(a + b)^(5/2)*d) - (a*sech(c + d*x)^2*tanh(c + d*x))/(4*b*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) - (3*a*(a + 2*b)*tanh(c + d*x))/(8*b^2*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [sech(c + d*x)^7/(a + b*sech(c + d*x)^2)^3, x, 6, atan(sinh(c + d*x))/(b^3*d) - (sqrt(a)*(8*a^2 + 20*a*b + 15*b^2)*atan((sqrt(a)*sinh(c + d*x))/sqrt(a + b)))/(8*b^3*(a + b)^(5/2)*d) - (a*sinh(c + d*x))/(4*b*(a + b)*d*(a + b + a*sinh(c + d*x)^2)^2) - (a*(4*a + 7*b)*sinh(c + d*x))/(8*b^2*(a + b)^2*d*(a + b + a*sinh(c + d*x)^2))]

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

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^2)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^2)^p=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [(a + b*sech(c + d*x)^2)*tanh(c + d*x)^4, x, 4, a*x - (a*tanh(c + d*x))/d - (a*tanh(c + d*x)^3)/(3*d) + (b*tanh(c + d*x)^5)/(5*d)]
    @test_int [(a + b*sech(c + d*x)^2)*tanh(c + d*x)^3, x, 4, (a*log(cosh(c + d*x)))/d + ((a - b)*sech(c + d*x)^2)/(2*d) + (b*sech(c + d*x)^4)/(4*d)]
    @test_int [(a + b*sech(c + d*x)^2)*tanh(c + d*x)^2, x, 4, a*x - (a*tanh(c + d*x))/d + (b*tanh(c + d*x)^3)/(3*d)]
    @test_int [(a + b*sech(c + d*x)^2)*tanh(c + d*x)^1, x, 3, (a*log(cosh(c + d*x)))/d - (b*sech(c + d*x)^2)/(2*d)]
    @test_int [(a + b*sech(c + d*x)^2)*tanh(c + d*x)^0, x, 3, a*x + (b*tanh(c + d*x))/d]
    @test_int [coth(c + d*x)^1*(a + b*sech(c + d*x)^2), x, 4, -((b*log(cosh(c + d*x)))/d) + ((a + b)*log(sinh(c + d*x)))/d]
    @test_int [coth(c + d*x)^2*(a + b*sech(c + d*x)^2), x, 4, a*x - ((a + b)*coth(c + d*x))/d]
    @test_int [coth(c + d*x)^3*(a + b*sech(c + d*x)^2), x, 4, -(((a + b)*csch(c + d*x)^2)/(2*d)) + (a*log(sinh(c + d*x)))/d]
    @test_int [coth(c + d*x)^4*(a + b*sech(c + d*x)^2), x, 4, a*x - (a*coth(c + d*x))/d - ((a + b)*coth(c + d*x)^3)/(3*d)]
    @test_int [coth(c + d*x)^5*(a + b*sech(c + d*x)^2), x, 4, -(((2*a + b)*csch(c + d*x)^2)/(2*d)) - ((a + b)*csch(c + d*x)^4)/(4*d) + (a*log(sinh(c + d*x)))/d]

    @test_int [(a + b*sech(c + d*x)^2)^2*tanh(c + d*x)^4, x, 4, a^2*x - (a^2*tanh(c + d*x))/d - (a^2*tanh(c + d*x)^3)/(3*d) + (b*(2*a + b)*tanh(c + d*x)^5)/(5*d) - (b^2*tanh(c + d*x)^7)/(7*d)]
    @test_int [(a + b*sech(c + d*x)^2)^2*tanh(c + d*x)^3, x, 4, (a^2*log(cosh(c + d*x)))/d + (a*(a - 2*b)*sech(c + d*x)^2)/(2*d) + ((2*a - b)*b*sech(c + d*x)^4)/(4*d) + (b^2*sech(c + d*x)^6)/(6*d)]
    @test_int [(a + b*sech(c + d*x)^2)^2*tanh(c + d*x)^2, x, 4, a^2*x - (a^2*tanh(c + d*x))/d + (b*(2*a + b)*tanh(c + d*x)^3)/(3*d) - (b^2*tanh(c + d*x)^5)/(5*d)]
    @test_int [(a + b*sech(c + d*x)^2)^2*tanh(c + d*x)^1, x, 4, (a^2*log(cosh(c + d*x)))/d - (a*b*sech(c + d*x)^2)/d - (b^2*sech(c + d*x)^4)/(4*d)]
    @test_int [(a + b*sech(c + d*x)^2)^2*tanh(c + d*x)^0, x, 4, a^2*x + (b*(2*a + b)*tanh(c + d*x))/d - (b^2*tanh(c + d*x)^3)/(3*d)]
    @test_int [coth(c + d*x)^1*(a + b*sech(c + d*x)^2)^2, x, 4, -((b*(2*a + b)*log(cosh(c + d*x)))/d) + ((a + b)^2*log(sinh(c + d*x)))/d + (b^2*sech(c + d*x)^2)/(2*d)]
    @test_int [coth(c + d*x)^2*(a + b*sech(c + d*x)^2)^2, x, 4, a^2*x - ((a + b)^2*coth(c + d*x))/d - (b^2*tanh(c + d*x))/d]
    @test_int [coth(c + d*x)^3*(a + b*sech(c + d*x)^2)^2, x, 4, -(((a + b)^2*csch(c + d*x)^2)/(2*d)) + (b^2*log(cosh(c + d*x)))/d + ((a^2 - b^2)*log(sinh(c + d*x)))/d]
    @test_int [coth(c + d*x)^4*(a + b*sech(c + d*x)^2)^2, x, 4, a^2*x - ((a^2 - b^2)*coth(c + d*x))/d - ((a + b)^2*coth(c + d*x)^3)/(3*d)]
    @test_int [coth(c + d*x)^5*(a + b*sech(c + d*x)^2)^2, x, 4, -((a*(a + b)*csch(c + d*x)^2)/d) - ((a + b)^2*csch(c + d*x)^4)/(4*d) + (a^2*log(sinh(c + d*x)))/d]
    @test_int [coth(c + d*x)^6*(a + b*sech(c + d*x)^2)^2, x, 4, a^2*x - (a^2*coth(c + d*x))/d - ((a^2 - b^2)*coth(c + d*x)^3)/(3*d) - ((a + b)^2*coth(c + d*x)^5)/(5*d)]
    @test_int [coth(c + d*x)^7*(a + b*sech(c + d*x)^2)^2, x, 5, -((a*(a + b)*csch(c + d*x)^2)/d) - ((a + b)^2*csch(c + d*x)^4)/(4*d) - ((b + a*cosh(c + d*x)^2)^3*csch(c + d*x)^6)/(6*(a + b)*d) + (a^2*log(sinh(c + d*x)))/d]

    @test_int [(a + b*sech(c + d*x)^2)^3*tanh(c + d*x)^4, x, 4, a^3*x - (a^3*tanh(c + d*x))/d - (a^3*tanh(c + d*x)^3)/(3*d) + (b*(3*a^2 + 3*a*b + b^2)*tanh(c + d*x)^5)/(5*d) - (b^2*(3*a + 2*b)*tanh(c + d*x)^7)/(7*d) + (b^3*tanh(c + d*x)^9)/(9*d)]
    @test_int [(a + b*sech(c + d*x)^2)^3*tanh(c + d*x)^3, x, 5, (a^3*log(cosh(c + d*x)))/d - (3*a^2*b*sech(c + d*x)^2)/(2*d) - (3*a*b^2*sech(c + d*x)^4)/(4*d) - (b^3*sech(c + d*x)^6)/(6*d) + ((b + a*cosh(c + d*x)^2)^4*sech(c + d*x)^8)/(8*b*d)]
    @test_int [(a + b*sech(c + d*x)^2)^3*tanh(c + d*x)^2, x, 4, a^3*x - (a^3*tanh(c + d*x))/d + (b*(3*a^2 + 3*a*b + b^2)*tanh(c + d*x)^3)/(3*d) - (b^2*(3*a + 2*b)*tanh(c + d*x)^5)/(5*d) + (b^3*tanh(c + d*x)^7)/(7*d)]
    @test_int [(a + b*sech(c + d*x)^2)^3*tanh(c + d*x)^1, x, 4, (a^3*log(cosh(c + d*x)))/d - (3*a^2*b*sech(c + d*x)^2)/(2*d) - (3*a*b^2*sech(c + d*x)^4)/(4*d) - (b^3*sech(c + d*x)^6)/(6*d)]
    @test_int [(a + b*sech(c + d*x)^2)^3*tanh(c + d*x)^0, x, 4, a^3*x + (b*(3*a^2 + 3*a*b + b^2)*tanh(c + d*x))/d - (b^2*(3*a + 2*b)*tanh(c + d*x)^3)/(3*d) + (b^3*tanh(c + d*x)^5)/(5*d)]
    @test_int [coth(c + d*x)^1*(a + b*sech(c + d*x)^2)^3, x, 4, -((b*(3*a^2 + 3*a*b + b^2)*log(cosh(c + d*x)))/d) + ((a + b)^3*log(sinh(c + d*x)))/d + (b^2*(3*a + b)*sech(c + d*x)^2)/(2*d) + (b^3*sech(c + d*x)^4)/(4*d)]
    @test_int [coth(c + d*x)^2*(a + b*sech(c + d*x)^2)^3, x, 4, a^3*x - ((a + b)^3*coth(c + d*x))/d - (b^2*(3*a + 2*b)*tanh(c + d*x))/d + (b^3*tanh(c + d*x)^3)/(3*d)]
    @test_int [coth(c + d*x)^3*(a + b*sech(c + d*x)^2)^3, x, 4, -(((a + b)^3*csch(c + d*x)^2)/(2*d)) + (b^2*(3*a + 2*b)*log(cosh(c + d*x)))/d + ((a - 2*b)*(a + b)^2*log(sinh(c + d*x)))/d - (b^3*sech(c + d*x)^2)/(2*d)]
    @test_int [coth(c + d*x)^4*(a + b*sech(c + d*x)^2)^3, x, 4, a^3*x - ((a - 2*b)*(a + b)^2*coth(c + d*x))/d - ((a + b)^3*coth(c + d*x)^3)/(3*d) + (b^3*tanh(c + d*x))/d]
    @test_int [coth(c + d*x)^5*(a + b*sech(c + d*x)^2)^3, x, 4, -(((2*a - b)*(a + b)^2*csch(c + d*x)^2)/(2*d)) - ((a + b)^3*csch(c + d*x)^4)/(4*d) - (b^3*log(cosh(c + d*x)))/d + ((a^3 + b^3)*log(sinh(c + d*x)))/d]
    @test_int [coth(c + d*x)^6*(a + b*sech(c + d*x)^2)^3, x, 4, a^3*x - ((a^3 + b^3)*coth(c + d*x))/d - ((a - 2*b)*(a + b)^2*coth(c + d*x)^3)/(3*d) - ((a + b)^3*coth(c + d*x)^5)/(5*d)]
    @test_int [coth(c + d*x)^7*(a + b*sech(c + d*x)^2)^3, x, 4, -((3*a^2*(a + b)*csch(c + d*x)^2)/(2*d)) - (3*a*(a + b)^2*csch(c + d*x)^4)/(4*d) - ((a + b)^3*csch(c + d*x)^6)/(6*d) + (a^3*log(sinh(c + d*x)))/d]

    @test_int [(a + b*sech(c + d*x)^2)^4, x, 4, a^4*x + (b*(2*a + b)*(2*a^2 + 2*a*b + b^2)*tanh(c + d*x))/d - (b^2*(6*a^2 + 8*a*b + 3*b^2)*tanh(c + d*x)^3)/(3*d) + (b^3*(4*a + 3*b)*tanh(c + d*x)^5)/(5*d) - (b^4*tanh(c + d*x)^7)/(7*d)]

    @test_int [(a + b*sech(c + d*x)^2)^5, x, 4, a^5*x + (b*(5*a^4 + 10*a^3*b + 10*a^2*b^2 + 5*a*b^3 + b^4)*tanh(c + d*x))/d - (b^2*(10*a^3 + 20*a^2*b + 15*a*b^2 + 4*b^3)*tanh(c + d*x)^3)/(3*d) + (b^3*(10*a^2 + 15*a*b + 6*b^2)*tanh(c + d*x)^5)/(5*d) - (b^4*(5*a + 4*b)*tanh(c + d*x)^7)/(7*d) + (b^5*tanh(c + d*x)^9)/(9*d)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [tanh(c + d*x)^5/(a + b*sech(c + d*x)^2), x, 4, -(((a + 2*b)*log(cosh(c + d*x)))/(b^2*d)) + ((a + b)^2*log(b + a*cosh(c + d*x)^2))/(2*a*b^2*d) - sech(c + d*x)^2/(2*b*d)]
    @test_int [tanh(c + d*x)^4/(a + b*sech(c + d*x)^2), x, 6, x/a - ((a + b)^(3/2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a*b^(3/2)*d) + tanh(c + d*x)/(b*d)]
    @test_int [tanh(c + d*x)^3/(a + b*sech(c + d*x)^2), x, 4, -(log(cosh(c + d*x))/(b*d)) + ((a + b)*log(b + a*cosh(c + d*x)^2))/(2*a*b*d)]
    @test_int [tanh(c + d*x)^2/(a + b*sech(c + d*x)^2), x, 5, x/a - (sqrt(a + b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a*sqrt(b)*d)]
    @test_int [tanh(c + d*x)^1/(a + b*sech(c + d*x)^2), x, 2, log(b + a*cosh(c + d*x)^2)/(2*a*d)]
    @test_int [tanh(c + d*x)^0/(a + b*sech(c + d*x)^2), x, 3, x/a - (sqrt(b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a*sqrt(a + b)*d)]
    @test_int [coth(c + d*x)^1/(a + b*sech(c + d*x)^2), x, 4, (b*log(b + a*cosh(c + d*x)^2))/(2*a*(a + b)*d) + log(sinh(c + d*x))/((a + b)*d)]
    @test_int [coth(c + d*x)^2/(a + b*sech(c + d*x)^2), x, 6, x/a - (b^(3/2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a*(a + b)^(3/2)*d) - coth(c + d*x)/((a + b)*d)]
    @test_int [coth(c + d*x)^3/(a + b*sech(c + d*x)^2), x, 4, -(csch(c + d*x)^2/(2*(a + b)*d)) + (b^2*log(b + a*cosh(c + d*x)^2))/(2*a*(a + b)^2*d) + ((a + 2*b)*log(sinh(c + d*x)))/((a + b)^2*d)]
    @test_int [coth(c + d*x)^4/(a + b*sech(c + d*x)^2), x, 7, x/a - (b^(5/2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(a*(a + b)^(5/2)*d) - ((a + 2*b)*coth(c + d*x))/((a + b)^2*d) - coth(c + d*x)^3/(3*(a + b)*d)]

    @test_int [tanh(c + d*x)^5/(a + b*sech(c + d*x)^2)^2, x, 4, (a + b)^2/(2*a^2*b*d*(b + a*cosh(c + d*x)^2)) + log(cosh(c + d*x))/(b^2*d) + ((a^(-2) - b^(-2))*log(b + a*cosh(c + d*x)^2))/(2*d)]
    @test_int [tanh(c + d*x)^4/(a + b*sech(c + d*x)^2)^2, x, 6, x/a^2 + ((a - 2*b)*sqrt(a + b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^2*b^(3/2)*d) - ((a + b)*tanh(c + d*x))/(2*a*b*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [tanh(c + d*x)^3/(a + b*sech(c + d*x)^2)^2, x, 4, (a + b)/(2*a^2*d*(b + a*cosh(c + d*x)^2)) + log(b + a*cosh(c + d*x)^2)/(2*a^2*d)]
    @test_int [tanh(c + d*x)^2/(a + b*sech(c + d*x)^2)^2, x, 6, x/a^2 - ((a + 2*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^2*sqrt(b)*sqrt(a + b)*d) - tanh(c + d*x)/(2*a*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [tanh(c + d*x)^1/(a + b*sech(c + d*x)^2)^2, x, 4, b/(2*a^2*d*(b + a*cosh(c + d*x)^2)) + log(b + a*cosh(c + d*x)^2)/(2*a^2*d)]
    @test_int [tanh(c + d*x)^0/(a + b*sech(c + d*x)^2)^2, x, 5, x/a^2 - (sqrt(b)*(3*a + 2*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^2*(a + b)^(3/2)*d) - (b*tanh(c + d*x))/(2*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [coth(c + d*x)^1/(a + b*sech(c + d*x)^2)^2, x, 4, b^2/(2*a^2*(a + b)*d*(b + a*cosh(c + d*x)^2)) + (b*(2*a + b)*log(b + a*cosh(c + d*x)^2))/(2*a^2*(a + b)^2*d) + log(sinh(c + d*x))/((a + b)^2*d)]
    @test_int [coth(c + d*x)^2/(a + b*sech(c + d*x)^2)^2, x, 7, x/a^2 - (b^(3/2)*(5*a + 2*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^2*(a + b)^(5/2)*d) - ((2*a - b)*coth(c + d*x))/(2*a*(a + b)^2*d) - (b*coth(c + d*x))/(2*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [coth(c + d*x)^3/(a + b*sech(c + d*x)^2)^2, x, 4, b^3/(2*a^2*(a + b)^2*d*(b + a*cosh(c + d*x)^2)) - csch(c + d*x)^2/(2*(a + b)^2*d) + (b^2*(3*a + b)*log(b + a*cosh(c + d*x)^2))/(2*a^2*(a + b)^3*d) + ((a + 3*b)*log(sinh(c + d*x)))/((a + b)^3*d)]
    @test_int [coth(c + d*x)^4/(a + b*sech(c + d*x)^2)^2, x, 8, x/a^2 - (b^(5/2)*(7*a + 2*b)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(2*a^2*(a + b)^(7/2)*d) - ((2*a^2 + 6*a*b - b^2)*coth(c + d*x))/(2*a*(a + b)^3*d) - ((2*a - 3*b)*coth(c + d*x)^3)/(6*a*(a + b)^2*d) - (b*coth(c + d*x)^3)/(2*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]

    @test_int [tanh(c + d*x)^6/(a + b*sech(c + d*x)^2)^3, x, 7, x/a^3 - (sqrt(a + b)*(3*a^2 - 4*a*b + 8*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^3*b^(5/2)*d) - ((a + b)*tanh(c + d*x)^3)/(4*a*b*d*(a + b - b*tanh(c + d*x)^2)^2) + ((3*a - 4*b)*(a + b)*tanh(c + d*x))/(8*a^2*b^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [tanh(c + d*x)^5/(a + b*sech(c + d*x)^2)^3, x, 4, -(a + b)^2/(4*a^3*d*(b + a*cosh(c + d*x)^2)^2) + (a + b)/(a^3*d*(b + a*cosh(c + d*x)^2)) + log(b + a*cosh(c + d*x)^2)/(2*a^3*d)]
    @test_int [tanh(c + d*x)^4/(a + b*sech(c + d*x)^2)^3, x, 7, x/a^3 + ((a^2 - 4*a*b - 8*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^3*b^(3/2)*sqrt(a + b)*d) - ((a + b)*tanh(c + d*x))/(4*a*b*d*(a + b - b*tanh(c + d*x)^2)^2) + ((a - 4*b)*tanh(c + d*x))/(8*a^2*b*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [tanh(c + d*x)^3/(a + b*sech(c + d*x)^2)^3, x, 4, -(b*(a + b))/(4*a^3*d*(b + a*cosh(c + d*x)^2)^2) + (a + 2*b)/(2*a^3*d*(b + a*cosh(c + d*x)^2)) + log(b + a*cosh(c + d*x)^2)/(2*a^3*d)]
    @test_int [tanh(c + d*x)^2/(a + b*sech(c + d*x)^2)^3, x, 7, x/a^3 - ((3*a^2 + 12*a*b + 8*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^3*sqrt(b)*(a + b)^(3/2)*d) - tanh(c + d*x)/(4*a*d*(a + b - b*tanh(c + d*x)^2)^2) - ((3*a + 4*b)*tanh(c + d*x))/(8*a^2*(a + b)*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [tanh(c + d*x)^1/(a + b*sech(c + d*x)^2)^3, x, 4, -b^2/(4*a^3*d*(b + a*cosh(c + d*x)^2)^2) + b/(a^3*d*(b + a*cosh(c + d*x)^2)) + log(b + a*cosh(c + d*x)^2)/(2*a^3*d)]
    @test_int [tanh(c + d*x)^0/(a + b*sech(c + d*x)^2)^3, x, 6, x/a^3 - (sqrt(b)*(15*a^2 + 20*a*b + 8*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^3*(a + b)^(5/2)*d) - (b*tanh(c + d*x))/(4*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) - (b*(7*a + 4*b)*tanh(c + d*x))/(8*a^2*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [coth(c + d*x)^1/(a + b*sech(c + d*x)^2)^3, x, 4, -(b^3/(4*a^3*(a + b)*d*(b + a*cosh(c + d*x)^2)^2)) + (b^2*(3*a + 2*b))/(2*a^3*(a + b)^2*d*(b + a*cosh(c + d*x)^2)) + (b*(3*a^2 + 3*a*b + b^2)*log(b + a*cosh(c + d*x)^2))/(2*a^3*(a + b)^3*d) + log(sinh(c + d*x))/((a + b)^3*d)]
    @test_int [coth(c + d*x)^2/(a + b*sech(c + d*x)^2)^3, x, 8, x/a^3 - (b^(3/2)*(35*a^2 + 28*a*b + 8*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^3*(a + b)^(7/2)*d) - ((8*a^2 - 11*a*b - 4*b^2)*coth(c + d*x))/(8*a^2*(a + b)^3*d) - (b*coth(c + d*x))/(4*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) - (b*(9*a + 4*b)*coth(c + d*x))/(8*a^2*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]
    @test_int [coth(c + d*x)^3/(a + b*sech(c + d*x)^2)^3, x, 4, -(b^4/(4*a^3*(a + b)^2*d*(b + a*cosh(c + d*x)^2)^2)) + (b^3*(2*a + b))/(a^3*(a + b)^3*d*(b + a*cosh(c + d*x)^2)) - csch(c + d*x)^2/(2*(a + b)^3*d) + (b^2*(6*a^2 + 4*a*b + b^2)*log(b + a*cosh(c + d*x)^2))/(2*a^3*(a + b)^4*d) + ((a + 4*b)*log(sinh(c + d*x)))/((a + b)^4*d)]
    @test_int [coth(c + d*x)^4/(a + b*sech(c + d*x)^2)^3, x, 9, x/a^3 - (b^(5/2)*(63*a^2 + 36*a*b + 8*b^2)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(8*a^3*(a + b)^(9/2)*d) - ((8*a^3 + 32*a^2*b - 15*a*b^2 - 4*b^3)*coth(c + d*x))/(8*a^2*(a + b)^4*d) - ((8*a^2 - 39*a*b - 12*b^2)*coth(c + d*x)^3)/(24*a^2*(a + b)^3*d) - (b*coth(c + d*x)^3)/(4*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^2) - (b*(11*a + 4*b)*coth(c + d*x)^3)/(8*a^2*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2))]

    @test_int [1/(a + b*sech(c + d*x)^2)^4, x, 7, x/a^4 - (sqrt(b)*(35*a^3 + 70*a^2*b + 56*a*b^2 + 16*b^3)*atanh((sqrt(b)*tanh(c + d*x))/sqrt(a + b)))/(16*a^4*(a + b)^(7/2)*d) - (b*tanh(c + d*x))/(6*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^3) - (b*(11*a + 6*b)*tanh(c + d*x))/(24*a^2*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2)^2) - (b*(19*a^2 + 22*a*b + 8*b^2)*tanh(c + d*x))/(16*a^3*(a + b)^3*d*(a + b - b*tanh(c + d*x)^2))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^2)^(p/2)*when*a+b=0=#

    @test_int [(1 - sech(x)^2)^(3/2), x, 4, coth(x)*log(cosh(x))*sqrt(tanh(x)^2) - (1/2)*coth(x)*(tanh(x)^2)^(3/2)]
    @test_int [sqrt(1 - sech(x)^2), x, 3, coth(x)*log(cosh(x))*sqrt(tanh(x)^2)]
    @test_int [1/sqrt(1 - sech(x)^2), x, 3, (log(sinh(x))*tanh(x))/sqrt(tanh(x)^2)]

    @test_int [(-1 + sech(x)^2)^(3/2), x, 4, (-coth(x))*log(cosh(x))*sqrt(-tanh(x)^2) + (1/2)*tanh(x)*sqrt(-tanh(x)^2)]
    @test_int [sqrt(-1 + sech(x)^2), x, 3, coth(x)*log(cosh(x))*sqrt(-tanh(x)^2)]
    @test_int [1/sqrt(-1 + sech(x)^2), x, 3, (log(sinh(x))*tanh(x))/sqrt(-tanh(x)^2)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^2)^(p/2)=#

    #= ::Subsubsection::Closed:: =#
    #=p>0=#

    @test_int [sqrt(a + b*sech(x)^2)*tanh(x)^5, x, 7, sqrt(a)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - sqrt(a + b*sech(x)^2) + ((a + 2*b)*(a + b*sech(x)^2)^(3/2))/(3*b^2) - (a + b*sech(x)^2)^(5/2)/(5*b^2)]
    @test_int [sqrt(a + b*sech(x)^2)*tanh(x)^4, x, 9, -((a^2 + 6*a*b - 3*b^2)*atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2)))/(8*b^(3/2)) + sqrt(a)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) + ((a - 3*b)*tanh(x)*sqrt(a + b - b*tanh(x)^2))/(8*b) - (tanh(x)^3*sqrt(a + b - b*tanh(x)^2))/4]
    @test_int [sqrt(a + b*sech(x)^2)*tanh(x)^3, x, 6, sqrt(a)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - sqrt(a + b*sech(x)^2) + (a + b*sech(x)^2)^(3/2)/(3*b)]
    @test_int [sqrt(a + b*sech(x)^2)*tanh(x)^2, x, 8, -((a - b)*atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2)))/(2*sqrt(b)) + sqrt(a)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) - (tanh(x)*sqrt(a + b - b*tanh(x)^2))/2]
    @test_int [sqrt(a + b*sech(x)^2)*tanh(x)^1, x, 5, sqrt(a)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - sqrt(a + b*sech(x)^2)]
    @test_int [sqrt(a + b*sech(x)^2)*tanh(x)^0, x, 6, sqrt(b)*atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) + sqrt(a)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))]
    @test_int [coth(x)^1*sqrt(a + b*sech(x)^2), x, 7, sqrt(a)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - sqrt(a + b)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b))]
    @test_int [coth(x)^2*sqrt(a + b*sech(x)^2), x, 6, sqrt(a)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) - coth(x)*sqrt(a + b - b*tanh(x)^2)]
    @test_int [coth(x)^3*sqrt(a + b*sech(x)^2), x, 8, sqrt(a)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - ((2*a + b)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b)))/(2*sqrt(a + b)) - (1/2)*coth(x)^2*sqrt(a + b*sech(x)^2)]
    @test_int [coth(x)^4*sqrt(a + b*sech(x)^2), x, 7, sqrt(a)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) - ((3*a + 2*b)*coth(x)*sqrt(a + b - b*tanh(x)^2))/(3*(a + b)) - (coth(x)^3*sqrt(a + b - b*tanh(x)^2))/3]
    @test_int [coth(x)^5*sqrt(a + b*sech(x)^2), x, 9, sqrt(a)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - ((8*a^2 + 12*a*b + 3*b^2)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b)))/(8*(a + b)^(3/2)) - ((4*a + 3*b)*coth(x)^2*sqrt(a + b*sech(x)^2))/(8*(a + b)) - (1/4)*coth(x)^4*sqrt(a + b*sech(x)^2)]

    @test_int [(a + b*sech(x)^2)^(3/2)*tanh(x)^3, x, 7, a^(3/2)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - a*sqrt(a + b*sech(x)^2) - (a + b*sech(x)^2)^(3/2)/3 + (a + b*sech(x)^2)^(5/2)/(5*b)]
    @test_int [(a + b*sech(x)^2)^(3/2)*tanh(x)^2, x, 9, -((3*a^2 - 6*a*b - b^2)*atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2)))/(8*sqrt(b)) + a^(3/2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) - ((5*a + b)*tanh(x)*sqrt(a + b - b*tanh(x)^2))/8 + (b*tanh(x)^3*sqrt(a + b - b*tanh(x)^2))/4]
    @test_int [(a + b*sech(x)^2)^(3/2)*tanh(x)^1, x, 6, a^(3/2)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - a*sqrt(a + b*sech(x)^2) - (a + b*sech(x)^2)^(3/2)/3]
    @test_int [(a + b*sech(x)^2)^(3/2)*tanh(x)^0, x, 7, (sqrt(b)*(3*a + b)*atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2)))/2 + a^(3/2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) + (b*tanh(x)*sqrt(a + b - b*tanh(x)^2))/2]
    @test_int [coth(x)^1*(a + b*sech(x)^2)^(3/2), x, 8, a^(3/2)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a)) - (a + b)^(3/2)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b)) + b*sqrt(a + b*sech(x)^2)]
    @test_int [coth(x)^2*(a + b*sech(x)^2)^(3/2), x, 8, -(b^(3/2)*atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2))) + a^(3/2)*atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2)) - (a + b)*coth(x)*sqrt(a + b - b*tanh(x)^2)]

    @test_int [(a + b*sech(c + d*x)^2)^(5/2), x, 8, (sqrt(b)*(15*a^2 + 10*a*b + 3*b^2)*atan((sqrt(b)*tanh(c + d*x))/sqrt(a + b - b*tanh(c + d*x)^2)))/(8*d) + (a^(5/2)*atanh((sqrt(a)*tanh(c + d*x))/sqrt(a + b - b*tanh(c + d*x)^2)))/d + (b*(7*a + 3*b)*tanh(c + d*x)*sqrt(a + b - b*tanh(c + d*x)^2))/(8*d) + (b*tanh(c + d*x)*(a + b - b*tanh(c + d*x)^2)^(3/2))/(4*d)]

    #= ::Subsubsection::Closed:: =#
    #=p<0=#

    @test_int [tanh(x)^5/sqrt(a + b*sech(x)^2), x, 6, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/sqrt(a) + ((a + 2*b)*sqrt(a + b*sech(x)^2))/b^2 - (a + b*sech(x)^2)^(3/2)/(3*b^2)]
    @test_int [tanh(x)^4/sqrt(a + b*sech(x)^2), x, 8, -((a + 3*b)*atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2)))/(2*b^(3/2)) + atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/sqrt(a) + (tanh(x)*sqrt(a + b - b*tanh(x)^2))/(2*b)]
    @test_int [tanh(x)^3/sqrt(a + b*sech(x)^2), x, 5, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/sqrt(a) + sqrt(a + b*sech(x)^2)/b]
    @test_int [tanh(x)^2/sqrt(a + b*sech(x)^2), x, 7, -(atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/sqrt(b)) + atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/sqrt(a)]
    @test_int [tanh(x)^1/sqrt(a + b*sech(x)^2), x, 4, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/sqrt(a)]
    @test_int [tanh(x)^0/sqrt(a + b*sech(x)^2), x, 3, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/sqrt(a)]
    @test_int [coth(x)^1/sqrt(a + b*sech(x)^2), x, 7, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/sqrt(a) - atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b))/sqrt(a + b)]
    @test_int [coth(x)^2/sqrt(a + b*sech(x)^2), x, 6, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/sqrt(a) - (coth(x)*sqrt(a + b - b*tanh(x)^2))/(a + b)]
    @test_int [coth(x)^3/sqrt(a + b*sech(x)^2), x, 8, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/sqrt(a) - ((2*a + 3*b)*atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b)))/(2*(a + b)^(3/2)) - (coth(x)^2*sqrt(a + b*sech(x)^2))/(2*(a + b))]

    @test_int [tanh(x)^5/(a + b*sech(x)^2)^(3/2), x, 6, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(3/2) - (a + b)^2/(a*b^2*sqrt(a + b*sech(x)^2)) - sqrt(a + b*sech(x)^2)/b^2]
    @test_int [tanh(x)^4/(a + b*sech(x)^2)^(3/2), x, 8, atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/b^(3/2) + atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(3/2) - ((a + b)*tanh(x))/(a*b*sqrt(a + b - b*tanh(x)^2))]
    @test_int [tanh(x)^3/(a + b*sech(x)^2)^(3/2), x, 5, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(3/2) - (a + b)/(a*b*sqrt(a + b*sech(x)^2))]
    @test_int [tanh(x)^2/(a + b*sech(x)^2)^(3/2), x, 5, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(3/2) - tanh(x)/(a*sqrt(a + b - b*tanh(x)^2))]
    @test_int [tanh(x)^1/(a + b*sech(x)^2)^(3/2), x, 5, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(3/2) - 1/(a*sqrt(a + b*sech(x)^2))]
    @test_int [tanh(x)^0/(a + b*sech(x)^2)^(3/2), x, 4, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(3/2) - (b*tanh(x))/(a*(a + b)*sqrt(a + b - b*tanh(x)^2))]
    @test_int [coth(x)^1/(a + b*sech(x)^2)^(3/2), x, 8, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(3/2) - atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b))/(a + b)^(3/2) - b/(a*(a + b)*sqrt(a + b*sech(x)^2))]
    @test_int [coth(x)^2/(a + b*sech(x)^2)^(3/2), x, 7, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(3/2) - (b*coth(x))/(a*(a + b)*sqrt(a + b - b*tanh(x)^2)) - ((a - b)*coth(x)*sqrt(a + b - b*tanh(x)^2))/(a*(a + b)^2)]

    @test_int [tanh(x)^6/(a + b*sech(x)^2)^(5/2), x, 9, -(atan((sqrt(b)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/b^(5/2)) + atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(5/2) - ((a + b)*tanh(x)^3)/(3*a*b*(a + b - b*tanh(x)^2)^(3/2)) - ((a^(-2) - b^(-2))*tanh(x))/sqrt(a + b - b*tanh(x)^2)]
    @test_int [tanh(x)^5/(a + b*sech(x)^2)^(5/2), x, 6, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(5/2) - (a + b)^2/(3*a*b^2*(a + b*sech(x)^2)^(3/2)) - (1/a^2 - 1/b^2)/sqrt(a + b*sech(x)^2)]
    @test_int [tanh(x)^4/(a + b*sech(x)^2)^(5/2), x, 7, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(5/2) - ((a + b)*tanh(x))/(3*a*b*(a + b - b*tanh(x)^2)^(3/2)) + ((a - 3*b)*tanh(x))/(3*a^2*b*sqrt(a + b - b*tanh(x)^2))]
    @test_int [tanh(x)^3/(a + b*sech(x)^2)^(5/2), x, 6, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(5/2) - (a + b)/(3*a*b*(a + b*sech(x)^2)^(3/2)) - 1/(a^2*sqrt(a + b*sech(x)^2))]
    @test_int [tanh(x)^2/(a + b*sech(x)^2)^(5/2), x, 7, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(5/2) - tanh(x)/(3*a*(a + b - b*tanh(x)^2)^(3/2)) - ((2*a + 3*b)*tanh(x))/(3*a^2*(a + b)*sqrt(a + b - b*tanh(x)^2))]
    @test_int [tanh(x)^1/(a + b*sech(x)^2)^(5/2), x, 6, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(5/2) - 1/(3*a*(a + b*sech(x)^2)^(3/2)) - 1/(a^2*sqrt(a + b*sech(x)^2))]
    @test_int [tanh(x)^0/(a + b*sech(x)^2)^(5/2), x, 6, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(5/2) - (b*tanh(x))/(3*a*(a + b)*(a + b - b*tanh(x)^2)^(3/2)) - (b*(5*a + 3*b)*tanh(x))/(3*a^2*(a + b)^2*sqrt(a + b - b*tanh(x)^2))]
    @test_int [coth(x)^1/(a + b*sech(x)^2)^(5/2), x, 9, atanh(sqrt(a + b*sech(x)^2)/sqrt(a))/a^(5/2) - atanh(sqrt(a + b*sech(x)^2)/sqrt(a + b))/(a + b)^(5/2) - b/(3*a*(a + b)*(a + b*sech(x)^2)^(3/2)) - (b*(2*a + b))/(a^2*(a + b)^2*sqrt(a + b*sech(x)^2))]
    @test_int [coth(x)^2/(a + b*sech(x)^2)^(5/2), x, 8, atanh((sqrt(a)*tanh(x))/sqrt(a + b - b*tanh(x)^2))/a^(5/2) - (b*coth(x))/(3*a*(a + b)*(a + b - b*tanh(x)^2)^(3/2)) - (b*(7*a + 3*b)*coth(x))/(3*a^2*(a + b)^2*sqrt(a + b - b*tanh(x)^2)) - ((a - 3*b)*(3*a + b)*coth(x)*sqrt(a + b - b*tanh(x)^2))/(3*a^2*(a + b)^3)]

    @test_int [1/(a + b*sech(c + d*x)^2)^(7/2), x, 7, atanh((sqrt(a)*tanh(c + d*x))/sqrt(a + b - b*tanh(c + d*x)^2))/(a^(7/2)*d) - (b*tanh(c + d*x))/(5*a*(a + b)*d*(a + b - b*tanh(c + d*x)^2)^(5/2)) - (b*(9*a + 5*b)*tanh(c + d*x))/(15*a^2*(a + b)^2*d*(a + b - b*tanh(c + d*x)^2)^(3/2)) - (b*(33*a^2 + 40*a*b + 15*b^2)*tanh(c + d*x))/(15*a^3*(a + b)^3*d*sqrt(a + b - b*tanh(c + d*x)^2))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^2)^p*when*p*symbolic=#

    #= ::Section:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^3)^p=#

    #= ::Section:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^4)^p=#

    #= ::Section:: =#
    #=Integrands*of*the*form*tanh(e+f*x)^m*(a+b*sech(e+f*x)^n)^p=#
    end

max_stars_repo_path: test/6 Hyperbolic functions/6.5 Hyperbolic secant/6.5.7 (d hyper)^m (a+b (c sech)^n)^p.jl


text:
"""
    ast_is(ast::VPtr, what::Symbol)::Bool

Helper for quickly recognizing kinds of ASTs
"""
ast_is(ast::VPtr, what::Symbol)::Bool = ccall(sbml(what), Cint, (VPtr,), ast) != 0

"""
    parse_math_children(ast::VPtr)::Vector{Math}

Recursively parse all children of an AST node.
"""
parse_math_children(ast::VPtr)::Vector{Math} = [
    parse_math(ccall(sbml(:ASTNode_getChild), VPtr, (VPtr, Cuint), ast, i - 1)) for
    i = 1:ccall(sbml(:ASTNode_getNumChildren), Cuint, (VPtr,), ast)
]

# Mapping of AST node type value subset to relational operations. Depends on
# `ASTNodeType.h` (also see below the case with AST_NAME_TIME)
const relational_opers = Dict{Int32,String}(
    308 => "eq",
    309 => "geq",
    310 => "gt",
    311 => "leq",
    312 => "lt",
    313 => "neq",
)

function relational_oper(t::Int)
    haskey(relational_opers, t) ||
        throw(DomainError(t, "Unknown ASTNodeType value for relational operator"))
    relational_opers[t]
end

"""
    parse_math(ast::VPtr)::Math

This attempts to parse out a decent Julia-esque ([`Math`](@ref) AST from a
pointer to `ASTNode_t`.
"""
function parse_math(ast::VPtr)::Math
    if ast_is(ast, :ASTNode_isName)
        if ccall(sbml(:ASTNode_getType), Cint, (VPtr,), ast) == 262
            # This is a special case checking for the value of "simulation
            # time" as defined by SBML. The constant `262` is the value of the
            # enum AST_NAME_TIME in `libsbml/src/sbml/math/ASTNodeType.h`,
            # needs to be kept up to date with the library (otherwise this
            # breaks).
            return MathTime(get_string(ast, :ASTNode_getName))
        else
            return MathIdent(get_string(ast, :ASTNode_getName))
        end
    elseif ast_is(ast, :ASTNode_isConstant)
        return MathConst(get_string(ast, :ASTNode_getName))
    elseif ast_is(ast, :ASTNode_isInteger)
        return MathVal(ccall(sbml(:ASTNode_getInteger), Cint, (VPtr,), ast))
    elseif ast_is(ast, :ASTNode_isRational)
        return MathVal(
            ccall(sbml(:ASTNode_getNumerator), Cint, (VPtr,), ast) //
            ccall(sbml(:ASTNode_getDenominator), Cint, (VPtr,), ast),
        )
    elseif ast_is(ast, :ASTNode_isReal)
        return MathVal(ccall(sbml(:ASTNode_getReal), Cdouble, (VPtr,), ast))
    elseif ast_is(ast, :ASTNode_isFunction)
        return MathApply(get_string(ast, :ASTNode_getName), parse_math_children(ast))
    elseif ast_is(ast, :ASTNode_isOperator)
        return MathApply(
            string(Char(ccall(sbml(:ASTNode_getCharacter), Cchar, (VPtr,), ast))),
            parse_math_children(ast),
        )
    elseif ast_is(ast, :ASTNode_isRelational)
        return MathApply(
            relational_oper(Int(ccall(sbml(:ASTNode_getType), Cint, (VPtr,), ast))),
            parse_math_children(ast),
        )
    elseif ast_is(ast, :ASTNode_isLogical)
        return MathApply(get_string(ast, :ASTNode_getName), parse_math_children(ast))
    elseif ast_is(ast, :ASTNode_isLambda)
        children = parse_math_children(ast)
        if !isempty(children)
            body = pop!(children)
            return MathLambda(broadcast((x::MathIdent) -> x.id, children), body)
        else
            @warn "invalid function definition found"
            return MathIdent("?invalid?")
        end
    else
        @warn "unsupported math element found"
        return MathIdent("?unsupported?")
    end
end

ext: jl


text:
using Documenter, PopGen

makedocs(;
    modules=[PopGen],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/pdimens/PopGen.jl/blob/{commit}{path}#L{line}",
    sitename="PopGen",
    authors="Pavel Dimens, Jason Selwyn",
    assets=String[],
)

deploydocs(;
    repo="github.com/pdimens/PopGen.jl",
)


text:
# Trace a path back through the network

@enum LegType transit transfer access egress

struct Leg
    start_time::Int32
    end_time::Int32
    origin_stop::Union{Missing, Int64}
    destination_stop::Union{Missing, Int64}
    type::LegType
    route::Union{Missing, Int64}
end

# trace the transit portion of the path
function trace_path(res::RaptorResult, stop::Int64)::Vector{Leg}
    legs = Vector{Leg}()

    current_stop = stop
    for round in size(res.times_at_stops_each_round, 1):-1:2
        time_after_round = res.times_at_stops_each_round[round, current_stop]
        time_before_round = res.times_at_stops_each_round[round - 1, current_stop]

        if time_after_round == time_before_round
            # not updated this round
            continue
        end

        prev_stop = res.prev_stop[round, current_stop]

        # okay, something happened this round. figure out if it was a transfer.
        local leg
        if round % 2 == 1
            # transfer rounds are round 3, 5, etc.
            # transfers are assumed to start the instant you get off the previous vehicle
            prev_time = res.times_at_stops_each_round[round - 1, prev_stop]
            leg = Leg(prev_time, time_after_round, prev_stop, current_stop, transfer, missing)
        else
            prev_route = res.prev_route[round, current_stop]
            prev_time = res.prev_boardtime[round, current_stop]
            leg = Leg(prev_time, time_after_round, prev_stop, current_stop, transit, prev_route)
        end

        push!(legs, leg)

        # prepare for previous iteration (we're going backwards)
        current_stop = prev_stop
    end

    reverse!(legs)

    return legs
end

function trace_path(res::StreetRaptorResult, destination::Int64)::Vector{Leg}
    # get the transit path
    dest_stop = res.egress_stop_for_destination[destination]
    legs = trace_path(res.raptor_result, dest_stop)

    # add the egress
    dest_time = res.times_at_destinations[destination]
    final_stop_arr_time = res.raptor_result.times_at_stops_each_round[size(res.raptor_result.times_at_stops_each_round, 1), dest_stop]
    egress_leg = Leg(final_stop_arr_time, dest_time, dest_stop, missing, egress, missing)
    push!(legs, egress_leg)

    # add the access
    initial_board_stop = legs[1].origin_stop
    # first round is access times
    arrival_time_at_initial_board_stop = res.raptor_result.times_at_stops_each_round[1, initial_board_stop]
    access_leg = Leg(res.request.departure_time, arrival_time_at_initial_board_stop, missing, initial_board_stop, access, missing)

    pushfirst!(legs, access_leg)

    return legs
end


text:
#
# This example shows the use of processes and channels with Simulate.jl
#
# pb, 2019-11-09
#

using Simulate, Printf, Random

function simple(input::Channel, output::Channel, name, id, op)
    token = take!(input)         # take something from the input
    now!(SF(println, @sprintf("%5.2f: %s %d took token %d", tau(), name, id, token)))
    d = delay!(rand())           # after a delay
    put!(output, op(token, id))  # put it out with some op applied
end

reset!(𝐶)
Random.seed!(123)

ch1 = Channel(32)  # create two channels
ch2 = Channel(32)

for i in 1:2:8    # create and register 8 SimProcesses SP
    process!(SP(i, simple, ch1, ch2, "foo", i, +))
    process!(SP(i+1, simple, ch2, ch1, "bar", i+1, *))
end

put!(ch1, 1)  # put first token into channel 1
run!(𝐶, 10)   # and run for 10 time units


text:
using ForwardDiff
using LinearAlgebra
using Plots
using Statistics

n = 1

# Params
mass = 1
g = -9.81
ℓ = 0.5

# Lagrangian
L(q,q̇) = 1/2*mass*ℓ^2*q̇^2 - mass*g*ℓ*(1-cos(q))

# Discrete Lagrangian (midpoint)
Ld(q,q_,h) = h*L((q+q_)/2, (q_-q)/h)
Ld2(q,q_,h) = h*mass*ℓ*( 1/2*ℓ*((q_-q)/h)^2 - g*(1-cos((q+q_)/2)) )
# Ld2(q,q_,h) = g*(1-cos((q+q_)/2))

q = rand()
q_ = rand()
h = 0.01
Ld(q,q_,h) - Ld2(q,q_,h)

# Discrete Euler Lagrange
D1(q,q_,h) = h*mass*ℓ*( -ℓ/h^2*(q_-q) - g/2*sin((q+q_)/2) )
D2(_q,q,h) = mass*ℓ^2/h*(q-_q) - mass*ℓ*h*g/2*sin((_q+q)/2)

Ld_aug(Q) = Ld2(Q[1],Q[2],h)

Q = [q,q_]
u = [0.0]
ForwardDiff.gradient(Ld_aug,Q)
D1(q,q_,h)
D2(q,q_,h)

function dyn_discrete(q1,q2,u; eps=1e-10)
    q3 = 2*q2-q1  # Guess using Euler

    # Discrete Euler-Lagrange
    P(q1,q2,q3,u) = D2(q1,q2,h) + D1(q2,q3,h) + u[1]
    D2_Ld = ForwardDiff.gradient(Ld_aug,[q1;q2])[n+1:end]

    for i = 1:10
        hessLd = ForwardDiff.hessian(Ld_aug,[q2;q3])
        gradLd = ForwardDiff.gradient(Ld_aug,[q2;q3])
        D1_Ld = gradLd[1:n]
        val = D2_Ld + D1_Ld + u
        q3 = q3 - hessLd[1:n,n+1:end]\val
        if norm(val,Inf) < eps
            break
        end
    end
    return q3
end
q1 = q
q2 = q + q_*h
q3 = dyn_discrete([q1],[q2],u)

# Continuous
L_aug(Q) = L(Q[1],Q[2])
gradL = ForwardDiff.gradient(L_aug,Q)
hessL = ForwardDiff.hessian(L_aug,Q)
function dyn(Q,u)
    gradL = ForwardDiff.gradient(L_aug,Q)
    hessL = ForwardDiff.hessian(L_aug,Q)
    Ldd = hessL[n+1:end,n+1:end]
    qdd = Ldd\(u + gradL[1:n] - hessL[1:n,n+1:end]*Q[n+1:end])
    return [Q[2]; qdd]
end
dyn(Q,u)

# Analytical
function dyn0(Q,u)
    return [Q[2]; u[1] - g/ℓ*sin(Q[1])]
end
dyn0(Q,u)
dyn0(Q,u) ≈ dyn(Q,u)

function rk3(f::Function, dt::Float64)
        # Runge-Kutta 3 (zero order hold)
    fd(x,u,dt=dt) = begin
        k1 = f(x,u)*dt
        k2 = f(x+k1/2,u)*dt
        k3 = f(x - k1 + 2*k2,u)*dt
        return x + (k1 + 4k2 + k3)/6
    end
end

dyn_discrete_rk3 = rk3(dyn,h)
dyn_discrete_rk3(Q,u)

function simulate1(Q0,N)
    Qs = zeros(length(Q0),N)
    Qs[:,1] = Q0
    for k = 2:N
        Qs[:,k] = dyn_discrete_rk3(Qs[:,k-1],[0])
    end
    return Qs
end

function simulate2(Q0,N)
    Q0 = [Q0[1:n];Q0[1:n]+Q[n+1:end]*h]
    Qs = zeros(n,N)
    Qs[:,1] = Q0[1:n]
    Qs[:,2] = Q0[n+1:end]
    for k = 3:N
        Qs[:,k] = dyn_discrete(Qs[:,k-2],Qs[:,k-1],[0])
    end
    return Qs
end

Q0 = [deg2rad(170),0]
res1 = simulate1(Q0,501)
plot(rad2deg.(res1[1,3:end]),label=:continuous)

res2 = simulate2(Q0,501)
plot!(rad2deg.(res2[1,:]),label=:discrete)


text:
import LinearAlgebra

@testset "U" begin
    @testset "ucase" begin
        @test SPICE._ucase("test") == uppercase("test")
    end
    @testset "ucrss" begin
        v1 = randn(3)
        v2 = randn(3)
        @test SPICE._ucrss(v1, v2) ≈ LinearAlgebra.normalize(LinearAlgebra.cross(v1, v2))
    end
    @testset "uddf" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :spk))
            et = str2et("JAN 1 2009")

            function udfunc(et_in)
                _, new_et = spkpos("MERCURY", et_in, "J2000", "LT+S", "MOON")
                new_et
            end

            deriv = uddf(udfunc, et, 1.0)
            @test deriv ≈ -0.000135670940 atol=1e-10
        finally
            kclear()
        end
    end
    @testset "union" begin
        # SPICEINT_CELL
        one = SpiceIntCell(8)
        two = SpiceIntCell(8)
        push!(one, 1, 2, 3)
        push!(two, 2, 3, 4)
        out = union(one, two)
        @test out == [1, 2, 3, 4]
        # SPICECHAR_CELL
        one = SpiceCharCell(8, 8)
        two = SpiceCharCell(8, 8)
        push!(one, "1", "2", "3")
        push!(two, "2", "3", "4")
        out = union(one, two)
        @test out == ["1", "2", "3", "4"]
        # SPICEDOUBLE_CELL
        one = SpiceDoubleCell(8)
        two = SpiceDoubleCell(8)
        push!(one, 1.0, 2.0, 3.0)
        push!(two, 2.0, 3.0, 4.0)
        out = union(one, two)
        @test out == [1.0, 2.0, 3.0, 4.0]
    end
    @testset "unitim" begin
        try
            furnsh(path(CORE, :lsk))
            et = str2et("Dec 19 2003")
            actual = unitim(et, :ET, :JED)
            expected = 2452992.5007428653
            @test actual ≈ expected
            @test_throws SpiceError unitim(et, :TEST, :JED)
            @test_throws SpiceError unitim(et, :ET, :TEST)
        finally
            kclear()
        end
    end
    @testset "unload" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck))
            @test ktotal(:ALL) == 2
            unload(path(CORE, :pck))
            @test ktotal(:ALL) == 1
        finally
            kclear()
        end
    end
    @testset "unorm" begin
        v1 = randn(3)
        @test all(SPICE._unorm(v1) .≈ (LinearAlgebra.normalize(v1),
                                       LinearAlgebra.norm(v1)))
        v1 = randn(8)
        @test all(SPICE._unormg(v1) .≈ (LinearAlgebra.normalize(v1),
                                        LinearAlgebra.norm(v1)))
    end
    @testset "utc2et" begin
        try
            furnsh(path(CORE, :lsk))
            utcstr = "December 1, 2004 15:04:11"
            output = utc2et(utcstr)
            @test output ≈ 155185515.1831043
            @test_throws SpiceError utc2et("")
            @test_throws SpiceError utc2et("1 1 1")
        finally
            kclear()
        end
    end
end


text:
# batch-wise matrix multiplication
# wrapper for batched_gemm!
export batched_mul, batched_transpose, batched_adjoint

include("./batchedadjtrans.jl")

"""
    batched_mul(A, B) -> C

Batched matrix multiplication. Result has `C[:,:,k] == A[:,:,k] * B[:,:,k]` for all `k`.
"""
function batched_mul(A::AbstractArray{T1, 3}, B::AbstractArray{T2, 3}) where {T1, T2}
    axes(A, 3) == axes(B, 3) || throw(DimensionMismatch("batch size mismatch"))
    T = promote_type(T1, T2)
    C = similar(A, T, (axes(A, 1), axes(B, 2), axes(A, 3)))
    batched_mul!(C, A, B)
end

"""
    batched_mul!(C, A, B) -> C

In-place batched matrix multiplication,
equivalent to `mul!(C[:,:,k], A[:,:,k], B[:,:,k])` for all `k`.
"""
function batched_mul! end

_unbatch(A) = A
_unbatch(A::BatchedAdjOrTrans) = A.parent

# batched_gemm!

const _GemmFloat = Union{Float64, Float32, ComplexF64, ComplexF32}

_BATCHED_GEMM_LIST = [
    (:(StridedArray{T, 3}), 'N'),
    (:(BatchedTranspose{T, <:StridedArray{T, 3}}), 'T'),
    (:(BatchedAdjoint{T, <:StridedArray{T, 3}}), 'C')
]

for (TA, transA) in _BATCHED_GEMM_LIST, (TB, transB) in _BATCHED_GEMM_LIST
    @eval function batched_mul!(C::StridedArray{T, 3}, A::$TA, B::$TB) where {T<:_GemmFloat}
        batched_gemm!($transA, $transB, one(T), _unbatch(A), _unbatch(B), zero(T), C)
        C
    end
end

# fallback

_BATCHED_LIST = [
    (:(AbstractArray{<:Any, 3}), :identity),
    (:(BatchedTranspose{<:Any, <:AbstractArray{<:Any, 3}}), :transpose),
    (:(BatchedAdjoint{<:Any, <:AbstractArray{<:Any, 3}}), :adjoint)
]
for (TA, fA) in _BATCHED_LIST, (TB, fB) in _BATCHED_LIST
    @eval function batched_mul!(C::AbstractArray{<:Any, 3}, A::$TA, B::$TB)
        axes(A, 3) == axes(B, 3) == axes(C, 3) || throw(DimensionMismatch("batch size mismatch"))
        @debug "calling fallback method for batched_mul!" typeof(A) typeof(B) typeof(C)
        A′, B′ = _unbatch(A), _unbatch(B)
        @inbounds for k in axes(C, 3)
            @views mul!(C[:,:,k], $fA(A′[:,:,k]), $fB(B′[:,:,k]))
        end
        C
    end
end


text:
"""
    SqExponentialKernel()

Squared exponential kernel.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the squared exponential kernel is defined as
```math
k(x, x') = \\exp\\bigg(- \\frac{\\|x - x'\\|_2^2}{2}\\bigg).
```

See also: [`GammaExponentialKernel`](@ref)
"""
struct SqExponentialKernel <: SimpleKernel end

kappa(::SqExponentialKernel, d²::Real) = exp(-d² / 2)

metric(::SqExponentialKernel) = SqEuclidean()

iskroncompatible(::SqExponentialKernel) = true

Base.show(io::IO, ::SqExponentialKernel) = print(io, "Squared Exponential Kernel")

## Aliases ##

"""
    RBFKernel()

Alias of [`SqExponentialKernel`](@ref).
"""
const RBFKernel = SqExponentialKernel

"""
    GaussianKernel()

Alias of [`SqExponentialKernel`](@ref).
"""
const GaussianKernel = SqExponentialKernel

"""
    SEKernel()

Alias of [`SqExponentialKernel`](@ref).
"""
const SEKernel = SqExponentialKernel

"""
    ExponentialKernel()

Exponential kernel.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the exponential kernel is defined as
```math
k(x, x') = \\exp\\big(- \\|x - x'\\|_2\\big).
```

See also: [`GammaExponentialKernel`](@ref)
"""
struct ExponentialKernel <: SimpleKernel end

kappa(::ExponentialKernel, d::Real) = exp(-d)

metric(::ExponentialKernel) = Euclidean()

iskroncompatible(::ExponentialKernel) = true

Base.show(io::IO, ::ExponentialKernel) = print(io, "Exponential Kernel")

## Aliases ##

"""
    LaplacianKernel()

Alias of [`ExponentialKernel`](@ref).
"""
const LaplacianKernel = ExponentialKernel

"""
    Matern12Kernel()

Alias of [`ExponentialKernel`](@ref).
"""
const Matern12Kernel = ExponentialKernel

"""
    GammaExponentialKernel(; γ::Real=2.0)

γ-exponential kernel with parameter `γ`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the γ-exponential kernel[^RW] with parameter
``\\gamma \\in (0, 2]`` is defined as
```math
k(x, x'; \\gamma) = \\exp\\big(- \\|x - x'\\|_2^{\\gamma}\\big).
```

See also: [`ExponentialKernel`](@ref), [`SqExponentialKernel`](@ref)

[^RW]: C. E. Rasmussen & C. K. I. Williams (2006). Gaussian Processes for Machine Learning.
"""
struct GammaExponentialKernel{Tγ<:Real} <: SimpleKernel
    γ::Vector{Tγ}
    function GammaExponentialKernel(; gamma::Real=2.0, γ::Real=gamma)
        @check_args(GammaExponentialKernel, γ, zero(γ) < γ ≤ 2, "γ ∈ (0, 2]")
        return new{typeof(γ)}([γ])
    end
end

@functor GammaExponentialKernel

kappa(κ::GammaExponentialKernel, d::Real) = exp(-d^first(κ.γ))

metric(::GammaExponentialKernel) = Euclidean()

iskroncompatible(::GammaExponentialKernel) = true

function Base.show(io::IO, κ::GammaExponentialKernel)
    return print(io, "Gamma Exponential Kernel (γ = ", first(κ.γ), ")")
end


text:
value(x) = x
cuify(x) = error("To use LinSolveGPUFactorize, you must do `using CuArrays`")
promote_u0(u0,p,t0) = u0
promote_tspan(u0,p,tspan,prob,kwargs) = tspan
get_tmp(x) = nothing
isdistribution(u0) = false

function SciMLBase.tmap(args...)
  error("Zygote must be added to differentiate Zygote? If you see this error, report it.")
end

function __init__()
  @require ApproxFun="28f2ccd6-bb30-5033-b560-165f7b14dc2f" begin
    eval_u0(u0::ApproxFun.Fun) = false
  end

  @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin

    promote_u0(u0::AbstractArray{<:Measurements.Measurement},p::AbstractArray{<:Measurements.Measurement},t0) = u0
    promote_u0(u0,p::AbstractArray{<:Measurements.Measurement},t0) = eltype(p).(u0)

    value(x::Type{Measurements.Measurement{T}}) where {T} = T
    value(x::Measurements.Measurement) = Measurements.value(x)

    @inline fastpow(x::Measurements.Measurement, y::Measurements.Measurement) = x^y

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Measurements.Measurement,t) = abs(Measurements.value(u))
  end

  @require MonteCarloMeasurements="0987c9cc-fe09-11e8-30f0-b96dd679fdca" begin

    promote_u0(u0::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},p::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},t0) = u0
    promote_u0(u0,p::AbstractArray{<:MonteCarloMeasurements.AbstractParticles},t0) = eltype(p).(u0)

    value(x::Type{MonteCarloMeasurements.AbstractParticles{T,N}}) where {T,N} = T
    value(x::MonteCarloMeasurements.AbstractParticles) = mean(x.particles)

    @inline fastpow(x::MonteCarloMeasurements.AbstractParticles, y::MonteCarloMeasurements.AbstractParticles) = x^y

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N}) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(value.(t)))))
    end
    @inline ODE_DEFAULT_NORM(u::MonteCarloMeasurements.AbstractParticles,t) = abs(value(u))
  end

  @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin
    # Support adaptive errors should be errorless for exponentiation
    value(x::Type{Unitful.AbstractQuantity{T,D,U}}) where {T,D,U} = T
    value(x::Unitful.AbstractQuantity) = x.val
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Unitful.AbstractQuantity,t) = abs(value(u))
    @inline UNITLESS_ABS2(x::Unitful.AbstractQuantity) = real(abs2(x)/oneunit(x)*oneunit(x))
  end

  @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    include("tracker.jl")
  end

  # Piracy, should get upstreamed
  @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
    cuify(x::AbstractArray) = CUDA.CuArray(x)
    function LinearAlgebra.ldiv!(x::CUDA.CuArray,_qr::CUDA.CUSOLVER.CuQR,b::CUDA.CuArray)
      _x = UpperTriangular(_qr.R) \ (_qr.Q' * reshape(b,length(b),1))
      x .= vec(_x)
      CUDA.unsafe_free!(_x)
      return x
    end
    # make `\` work
    LinearAlgebra.ldiv!(F::CUDA.CUSOLVER.CuQR, b::CUDA.CuArray) = (x = similar(b); ldiv!(x, F, b); x)
    default_factorize(A::CUDA.CuArray) = qr(A)

    ODE_DEFAULT_NORM(u::CUDA.CuArray,t) = sqrt(real(sum(abs2,u))/length(u))

    @inline function ODE_DEFAULT_NORM(u::CUDA.CuArray{<:ForwardDiff.Dual},t)
      sqrt(sum(abs2,value.(u)) / length(u))
    end

    @inline function ODE_DEFAULT_NORM(u::CUDA.CuArray{<:ForwardDiff.Dual},t::ForwardDiff.Dual)
      sqrt(sum(abs2,u) / length(u))
    end
  end

  @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
    include("reversediff.jl")
  end

  @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
    function ∇tmap(cx, f, args...)
      ys_and_backs = SciMLBase.tmap((args...) -> Zygote._pullback(cx, f, args...), args...)
      if isempty(ys_and_backs)
        ys_and_backs, _ -> (NoTangent(),NoTangent())
      else
        ys, backs = Zygote.unzip(ys_and_backs)
        function ∇tmap_internal(Δ)
          Δf_and_args_zipped = SciMLBase.tmap((f, δ) -> f(δ), backs, Δ)
          Δf_and_args = Zygote.unzip(Δf_and_args_zipped)
          Δf = reduce(Zygote.accum, Δf_and_args[1])
          (Δf, Δf_and_args[2:end]...)
        end
        ys,∇tmap_internal
      end
    end

    function ∇responsible_map(cx, f, args...)
      ys_and_backs = SciMLBase.responsible_map((args...) -> Zygote._pullback(cx, f, args...), args...)
      if isempty(ys_and_backs)
        ys_and_backs, _ -> (NoTangent(),NoTangent())
      else
        ys, backs = Zygote.unzip(ys_and_backs)
        ys, function ∇responsible_map_internal(Δ)
          # Apply pullbacks in reverse order. Needed for correctness if `f` is stateful.
          Δf_and_args_zipped = SciMLBase.responsible_map((f, δ) -> f(δ), Zygote._tryreverse(SciMLBase.responsible_map, backs, Δ)...)
          Δf_and_args = Zygote.unzip(Zygote._tryreverse(SciMLBase.responsible_map, Δf_and_args_zipped))
          Δf = reduce(Zygote.accum, Δf_and_args[1])
          (Δf, Δf_and_args[2:end]...)
        end
      end
    end
  end

  @require GeneralizedGenerated="6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb" begin
    SciMLBase.numargs(::GeneralizedGenerated.RuntimeFn{Args}) where Args = GeneralizedGenerated.from_type(Args) |> length
  end

  @require Pardiso="46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" begin
    mutable struct MKLPardisoFactorize
      A
      ps::Pardiso.MKLPardisoSolver
      verbose::Bool
      firsttime::Bool
    end
    MKLPardisoFactorize(;verbose=false) = MKLPardisoFactorize(nothing,Pardiso.MKLPardisoSolver(),verbose,true)
    function (p::MKLPardisoFactorize)(x,A,b,update_matrix=false;kwargs...)
      if p.firsttime
        Pardiso.set_phase!(p.ps, Pardiso.ANALYSIS)
        Pardiso.pardiso(p.ps, x, A, b)
        p.firsttime = false
      end

      if update_matrix
        Pardiso.set_phase!(p.ps, Pardiso.NUM_FACT)
        Pardiso.pardiso(p.ps, x, A, b)
        p.A = A
      end

      Pardiso.set_phase!(p.ps, Pardiso.SOLVE_ITERATIVE_REFINE)
      Pardiso.pardiso(p.ps, x, A, b)
    end
    function (p::MKLPardisoFactorize)(::Type{Val{:init}},f,u0_prototype)
      if eltype(u0_prototype) <: Complex
        mattype = Pardiso.COMPLEX_NONSYM
      else
        mattype = Pardiso.REAL_NONSYM
      end
      Pardiso.set_matrixtype!(p.ps, mattype)
      if p.verbose
          Pardiso.set_msglvl!(p.ps, Pardiso.MESSAGE_LEVEL_ON)
      end
      MKLPardisoFactorize(nothing,p.ps,p.verbose,true)
    end

    mutable struct PardisoFactorize
      A
      ps::Pardiso.PardisoSolver
      verbose::Bool
      firsttime::Bool
    end
    PardisoFactorize(;verbose=false) = PardisoFactorize(nothing,Pardiso.PardisoSolver(),verbose,true)
    function (p::PardisoFactorize)(x,A,b,update_matrix=false;kwargs...)
      if p.firsttime
        Pardiso.set_phase!(p.ps, Pardiso.ANALYSIS)
        Pardiso.pardiso(p.ps, x, A, b)
        p.firsttime = false
      end

      if update_matrix
        Pardiso.set_phase!(p.ps, Pardiso.NUM_FACT)
        Pardiso.pardiso(p.ps, x, A, b)
        p.A = A
      end

      Pardiso.set_phase!(p.ps, Pardiso.SOLVE_ITERATIVE_REFINE)
      Pardiso.pardiso(p.ps, x, A, b)
    end
    function (p::PardisoFactorize)(::Type{Val{:init}},f,u0_prototype)
      if eltype(u0_prototype) <: Complex
        mattype = Pardiso.COMPLEX_NONSYM
      else
        mattype = Pardiso.REAL_NONSYM
      end
      Pardiso.set_matrixtype!(p.ps, mattype)
      if p.verbose
          Pardiso.set_msglvl!(p.ps, Pardiso.MESSAGE_LEVEL_ON)
      end
      PardisoFactorize(nothing,p.ps,p.verbose,true)
    end

    mutable struct PardisoIterate
      A
      ps::Pardiso.PardisoSolver
      verbose::Bool
      firsttime::Bool
    end
    PardisoIterate(;verbose=false) = PardisoIterate(nothing,Pardiso.PardisoSolver(),verbose,true)
    function (p::PardisoIterate)(x,A,b,update_matrix=false;kwargs...)
      if p.firsttime
        Pardiso.set_phase!(p.ps, Pardiso.ANALYSIS)
        Pardiso.pardiso(p.ps, x, A, b)
        p.firsttime = false
      end

      #=
      # Pardiso iterative solver doesn't need a separate factorization phase
      # It will auto-update as it needs to.
      if update_matrix
        Pardiso.set_phase!(p.ps, Pardiso.NUM_FACT)
        Pardiso.pardiso(p.ps, x, A, b)
        p.A = A
      end
      =#

      Pardiso.set_phase!(ps, Pardiso.NUM_FACT_SOLVE_REFINE)
      Pardiso.pardiso(ps, X, A, B)
    end
    function (p::PardisoIterate)(::Type{Val{:init}},f,u0_prototype)
      if eltype(u0_prototype) <: Complex
        mattype = Pardiso.COMPLEX_NONSYM
      else
        mattype = Pardiso.REAL_NONSYM
      end
      Pardiso.set_matrixtype!(p.ps, mattype)
      Pardiso.set_solver!(ps, Pardiso.ITERATIVE_SOLVER)
      if p.verbose
          Pardiso.set_msglvl!(p.ps, Pardiso.MESSAGE_LEVEL_ON)
      end
      PardisoIterate(nothing,p.ps,p.verbose,true)
    end
  end
end


text:
"""
`module  JAC.ImpactExcitation`
    ... a submodel of JAC that contains all methods for computing electron impact excitation cross sections and collision strengths.
"""
module ImpactExcitation

    using Printf, ..AngularMomentum, ..Basics, ..Defaults, ..ManyElectron, ..Radial, ..Nuclear, ..TableStrings

    """
    `struct  ImpactExcitationSettings`  ... defines a type for the details and parameters of computing electron-impact excitation lines.

        + electronEnergies        ::Array{Float64,1}             ... List of impact-energies of the incoming elecgtrons.
        + includeBreit            ::Bool                         ... True if the Breit interaction is to be included, and false otherwise.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True if particular lines are selected for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
        + maxKappa                ::Int64                        ... Maximum kappa value of partial waves to be included.
        + energyShift             ::Float64                      ... An overall energy shift for all transitions |i> --> |f>.
    """
    struct Settings
        electronEnergies          ::Array{Float64,1}
        includeBreit              ::Bool
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
        maxKappa                  ::Int64
        energyShift               ::Float64
    end

    """
    `ImpactExcitation.Settings()`  ... constructor for the default values of electron-impact excitation line computations.
    """
    function Settings()
       Settings( Float64[], false, false, false, Tuple{Int64,Int64}[], 0, 0.)
    end

    # `Base.show(io::IO, settings::ImpactExcitation.Settings)`  ... prepares a proper printout of the variable settings::ImpactExcitation.Settings.
    function Base.show(io::IO, settings::ImpactExcitation.Settings)
        println(io, "electronEnergies:           $(settings.electronEnergies)  ")
        println(io, "includeBreit:               $(settings.includeBreit)  ")
        println(io, "printBeforeComputation:     $(settings.printBeforeComputation)  ")
        println(io, "selectLines:                $(settings.selectLines)  ")
        println(io, "selectedLines:              $(settings.selectedLines)  ")
        println(io, "maxKappa:                   $(settings.maxKappa)  ")
        println(io, "energyShift:                $(settings.energyShift)  ")
    end

    """
    `struct  ImpactExcitation.Channel`
        ... defines a type for a electron-impact excitaiton channel to help characterize the incoming and outgoing (continuum) states of
            many electron-states with a single free electron

        + initialKappa     ::Int64              ... partial-wave of the incoming free electron
        + finalKappa       ::Int64              ... partial-wave of the outgoing free electron
        + symmetry         ::LevelSymmetry      ... total angular momentum and parity of the scattering state
        + initialPhase     ::Float64            ... phase of the incoming partial wave
        + finalPhase       ::Float64            ... phase of the outgoing partial wave
        + amplitude        ::Complex{Float64}   ... Collision amplitude associated with the given channel.
    """
    struct  Channel
        initialKappa       ::Int64
        finalKappa         ::Int64
        symmetry           ::LevelSymmetry
        initialPhase       ::Float64
        finalPhase         ::Float64
        amplitude          ::Complex{Float64}
    end

    # `Base.show(io::IO, channel::ImpactExcitation.Channel)`  ... prepares a proper printout of the variable channel::ImpactExcitation.Channel.
    function Base.show(io::IO, channel::ImpactExcitation.Channel)
        println(io, "initialKappa:       $(channel.initialKappa)  ")
        println(io, "finalKappa:         $(channel.finalKappa)  ")
        println(io, "symmetry:           $(channel.symmetry)  ")
        println(io, "initialPhase:       $(channel.initialPhase)  ")
        println(io, "finalPhase:         $(channel.finalPhase)  ")
        println(io, "amplitude:          $(channel.amplitude)  ")
    end

    """
    `struct  ImpactExcitation.Line`
        ... defines a type for a electron-impact excitation line that may include the definition of channels and their corresponding amplitudes.

        + initialLevel           ::Level         ... initial- (bound-state) level
        + finalLevel             ::Level         ... final- (bound-state) level
        + initialElectronEnergy  ::Float64       ... energy of the incoming (initial-state) free-electron
        + finalElectronEnergy    ::Float64       ... energy of the outgoing (final-state) free-electron
        + crossSection           ::Float64       ... total cross section of this line
        + hasChannels            ::Bool          ... Determines whether the individual scattering (sub-) channels are defined in terms of their
                                                     kappa's, phases and the total angular momentum/parity as well as the amplitude, or not.
        + channels               ::Array{ImpactExcitation.Channel,1}  ... List of ImpactExcitation channels of this line.
    """
    struct  Line
        initialLevel             ::Level
        finalLevel               ::Level
        initialElectronEnergy    ::Float64
        finalElectronEnergy      ::Float64
        crossSection             ::Float64
        hasChannels              ::Bool
        channels                 ::Array{ImpactExcitation.Channel,1}
    end

    """
    `ImpactExcitation.Line()`  ... 'empty' constructor for an electron-impact excitation line between a specified initial and final level.
    """
    function Line()
        Line(Level(), Level(), 0., 0., 0., false, ImpactExcitation.Channel[] )
    end

    """
    `ImpactExcitation.Line(initialLevel::Level, finalLevel::Level, crossSection::Float64)`
        ... constructor for an electron-impact excitation line between a specified initial and final level.
    """
    function Line(initialLevel::Level, finalLevel::Level, crossSection::Float64)
        Line(initialLevel, finalLevel, 0., 0., crossSection, false, ImpactExcitation.Channel[] )
    end

    # `Base.show(io::IO, line::ImpactExcitation.Line)`  ... prepares a proper printout of the variable line::ImpactExcitation.Line.
    function Base.show(io::IO, line::ImpactExcitation.Line)
        println(io, "initialLevel:            $(line.initialLevel)  ")
        println(io, "finalLevel:              $(line.finalLevel)  ")
        println(io, "initialElectronEnergy:   $(line.initialElectronEnergy)  ")
        println(io, "finalElectronEnergy:     $(line.finalElectronEnergy)  ")
        println(io, "crossSection:            $(line.crossSection)  ")
        println(io, "hasChannels:             $(line.hasChannels)  ")
        println(io, "channels:                $(line.channels)  ")
    end

    """
    `ImpactExcitation.computeAmplitudesProperties(line::ImpactExcitation.Line, grid::Radial.Grid, settings::ImpactExcitation.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::ImpactExcitation.Line is returned for which the amplitudes and
            properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::ImpactExcitation.Line, grid::Radial.Grid, settings::ImpactExcitation.Settings)
        global JAC_counter
        newChannels = ImpactExcitation.Channel[]
        for channel in line.channels
            # Generate two continuum orbitals
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("ImpactExcitation.computeAmplitudesProperties-aa: warning ... no coninuum orbitals are generated.") end
            initialPhase = 0.;    finalPhase = 0.
            # Define a proper continuum basis from the initialLevel.basis and finalLevel.basis with the two continuum orbital
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("ImpactExcitation.computeAmplitudesProperties-ab: warning ... no continuum bases are generated.") end
            # Compute the transition matrix for the two constructed continuum bases
            JAC_counter = JAC_counter + 1
            if   JAC_counter < 20   println("ImpactExcitation.computeAmplitudesProperties-ac: warning ... no transition matrix is computed.") end
            # matrix    = ImpactExcitation.computeMatrix(channel.multipole, channel.gauge, line.omega, line.finalLevel.basis,
            #                                                line.initialLevel.basis, grid, settings)
            # amplitude = line.finalLevel.mc * matrix * line.initialLevel.mc
            amplitude = 1.0
            push!( newChannels, ImpactExcitation.Channel( channel.initialKappa, channel.finalKappa, channel.symmetry,
                                                          initialPhase, finalPhase, amplitude) )
        end
        # Calculate the electron-impact excitation strength and cross section
        JAC_counter = JAC_counter + 1
        if   JAC_counter < 20   println("ImpactExcitation.computeAmplitudesProperties-ba: warning ... cs and strength set to -1.") end
        crossSection = -1.
        line = ImpactExcitation.Line( line.initialLevel, line.finalLevel, line.initialElectronEnergy, line.finalElectronEnergy,
                                      crossSection, true, newChannels)
        return( line )
    end

    """
    `ImpactExcitation.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                   settings::ImpactExcitation.Settings; output=true)`
        ... to compute the electron-impact excitation transition amplitudes and all properties as requested by the given settings.
            A list of lines::Array{ImpactExcitation.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::ImpactExcitation.Settings;
                           output=true)
        println("")
        printstyled("ImpactExcitation.computePathways(): The computation of electron-impact excitation cross sections starts now ... \n", color=:light_green)
        printstyled("--------------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        lines = ImpactExcitation.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBeforeComputation    ImpactExcitation.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = ImpactExcitation.Line[]
        for  line in lines
            newLine = ImpactExcitation.computeAmplitudesProperties(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        ImpactExcitation.displayResults(lines)
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `ImpactExcitation.computeMatrix(finalBasis::Basis, initialBasis::Basis, settings::ImpactExcitation.Settings)`
        ... to compute the transition matrix  (<finalContinuumCSF_r|| V(e-e) ||initialContinuumCSF_s>)  between the CSF_r from the
            finalContinuumBasis and the CSF_s from the initialContinuumBasis. A (non-quadratic) matrix::Array{Float64,2} with
            dimensions [length(finalContinuumBasis.csfs) x length(initialContinuumBasis.csfs)] is returned. Note that this transition
            matrix is typically specific to just one Eimex channel due to the different energies, partial waves and overall symmetry
            of the scattering states. **Not yet implemented !**
    """
    function computeMatrix(finalBasis::Basis, initialBasis::Basis, settings::ImpactExcitation.Settings)
        error("Not yet implemented.")
    end

    """
    `ImpactExcitation.determineChannels(finalLevel::Level, initialLevel::Level, settings::ImpactExcitation.Settings)`
        ... to determine a list of electron-impact excitation Channels for a transitions from the initial to the final level and by
            taking into account the particular settings of for this computation; an Array{ImpactExcitation.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, settings::ImpactExcitation.Settings)
        channels = ImpactExcitation.Channel[];
        symi = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        #
        for  inKappa = -(settings.maxKappa+1):settings.maxKappa
            if  inKappa == 0    continue    end
            symtList = AngularMomentum.allowedTotalSymmetries(symi, inKappa)
            for  symt in symtList
                outKappaList = AngularMomentum.allowedKappaSymmetries(symt, symf)
                for  outKappa in outKappaList
                    push!(channels, ImpactExcitation.Channel( inKappa, outKappa, symt, 0., 0.,Complex(0.)) )
                end
            end
        end
        return( channels )
    end

    """
    `ImpactExcitation.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::ImpactExcitation.Settings)`
        ... to determine a list of ImpactExcitation.Line's for transitions between levels from the initial- and final-state multiplets,
            and by taking into account the particular selections and settings for this computation; an Array{ImpactExcitation.Line,1} is
            returned. Apart from the level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::ImpactExcitation.Settings)
        if    settings.selectLines    selectLines   = true;   selectedLines = Basics.determine("selected lines", settings.selectedLines)
        else                          selectLines   = false
        end

        lines = ImpactExcitation.Line[]
        for  i = 1:length(initialMultiplet.levels)
            for  f = 1:length(finalMultiplet.levels)
                if  selectLines  &&  !(haskey(selectedLines, (i,f)) )    continue   end
                for  en in settings.electronEnergies
                    initialElectronEnergy  = en
                    finalElectronEnergy    = en - (initialMultiplet.levels[i].energy - finalMultiplet.levels[f].energy) + settings.energyShift
                    if  finalElectronEnergy < 0    continue   end

                    channels = ImpactExcitation.determineChannels(finalMultiplet.levels[f], initialMultiplet.levels[i], settings)
                    push!( lines, ImpactExcitation.Line(initialMultiplet.levels[i], finalMultiplet.levels[f],
                                                        initialElectronEnergy, finalElectronEnergy, 0., true, channels) )
                end
            end
        end
        return( lines )
    end

    """
    `ImpactExcitation.displayLines(lines::Array{ImpactExcitation.Line,1})`
        ... to display a list of lines and channels that have been selected due to the prior settings. A neat table of all
            selected transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{ImpactExcitation.Line,1})
        println(" ")
        println("  Selected electron-impact ionization lines:")
        println(" ")
        println("  ", TableStrings.hLine(180))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                                sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=3);                                sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(12, "f--Energy--i"; na=4)
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=5)
        sa = sa * TableStrings.center(12, "Energy e_in"; na=3);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(12, "Energy e_out"; na=4);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(57, "List of partial waves and total symmetries"; na=4)
        sb = sb * TableStrings.flushleft(57, "partial-in [total J^P] partial-out        "; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(180))
        #
        NoLines = 0
        for  line in lines
            NoLines = NoLines + 1
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            en = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", en))                          * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.initialElectronEnergy))  * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.finalElectronEnergy))    * "    "
            kappaInOutSymmetryList = Tuple{Int64,Int64,LevelSymmetry}[]
            for  i in 1:length(line.channels)
                push!( kappaInOutSymmetryList, (line.channels[i].initialKappa, line.channels[i].finalKappa, line.channels[i].symmetry) )
            end
            wa = TableStrings.kappaKappaSymmetryTupels(95, kappaInOutSymmetryList)
            sb = sa * wa[1];    println( sb )
            for  i = 2:length(wa)
                NoLines = NoLines + 1
                sb = TableStrings.hBlank( length(sa) ) * wa[i];    println( sb )
            end
            # Avoid long table printouts
            if NoLines > 100  println("\n  ImpactExcitation.displayLines():  A maximum of 100 lines are printed in this table. \n")
               break
            end
        end
        println("  ", TableStrings.hLine(180))
        #
        return( nothing )
    end

    """
    `ImpactExcitation.displayResults(lines::Array{ImpactExcitation.Line,1})`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayResults(lines::Array{ImpactExcitation.Line,1})
        println(" ")
        println("  Electron-impact excitation cross sections:")
        println(" ")
        println("  ", TableStrings.hLine(113))
        sa = "  ";   sb = "  "
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                                sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=3);                                sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(12, "f--Energy--i"; na=4)
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=5)
        sa = sa * TableStrings.center(12, "Energy e_in"; na=3);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(12, "Energy e_out"; na=3);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(15, "Cross section"; na=3)
        sb = sb * TableStrings.center(15, TableStrings.inUnits("cross section"); na=3)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(113))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            en = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", en))                          * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.initialElectronEnergy))  * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.finalElectronEnergy))    * "     "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection))    * "    "
            println(sa)
        end
        println("  ", TableStrings.hLine(113))
        #
        return( nothing )
    end

end # module