text:
@testset "5.3.7 Inverse tangent functions" begin
    (a, b, c, d, e, f, m, n, x, ) = @variables a b c d e f m n x

    #= ::Package:: =#

    #= ::Title::Closed:: =#
    #=Integrands*involving*inverse*tangents*of*algebraic*functions=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*u*atan(a+b*x^n)=#

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

    @test_int [x^(n-1)*atan(a + b*x^n), x, 4, ((a + b*x^n)*atan(a + b*x^n))/(b*n) - log(1 + (a + b*x^n)^2)/(2*b*n)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(f*x)^m*(a+b*atan(c*x/sqrt(d+e*x^2)))=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*(a+b*atan(c*x/sqrt(d+e*x^2)))*when*e=c^2=#

    @test_int [x^5*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 6, (5*d^2*x*sqrt(d + e*x^2))/(96*(-e)^(5/2)) + (5*d*x^3*sqrt(d + e*x^2))/(144*(-e)^(3/2)) + (x^5*sqrt(d + e*x^2))/(36*sqrt(-e)) + (x^6*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/6 + (5*d^3*sqrt(-e)*atanh((sqrt(e)*x)/sqrt(d + e*x^2)))/(96*e^(7/2))]
    @test_int [x^3*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 5, (3*d*x*sqrt(d + e*x^2))/(32*(-e)^(3/2)) + (x^3*sqrt(d + e*x^2))/(16*sqrt(-e)) + (x^4*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/4 - (3*d^2*sqrt(-e)*atanh((sqrt(e)*x)/sqrt(d + e*x^2)))/(32*e^(5/2))]
    @test_int [x^1*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 4, (x*sqrt(d + e*x^2))/(4*sqrt(-e)) + (x^2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/2 + (d*sqrt(-e)*atanh((sqrt(e)*x)/sqrt(d + e*x^2)))/(4*e^(3/2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^1, x, 8, -((sqrt(d)*sqrt(-e)*sqrt(1 + (e*x^2)/d)*asinh((sqrt(e)*x)/sqrt(d))^2)/(2*sqrt(e)*sqrt(d + e*x^2))) + (sqrt(d)*sqrt(-e)*sqrt(1 + (e*x^2)/d)*asinh((sqrt(e)*x)/sqrt(d))*log(1 - ℯ^(2*asinh((sqrt(e)*x)/sqrt(d)))))/(sqrt(e)*sqrt(d + e*x^2)) - (sqrt(d)*sqrt(-e)*sqrt(1 + (e*x^2)/d)*asinh((sqrt(e)*x)/sqrt(d))*log(x))/(sqrt(e)*sqrt(d + e*x^2)) + atan((sqrt(-e)*x)/sqrt(d + e*x^2))*log(x) + (sqrt(d)*sqrt(-e)*sqrt(1 + (e*x^2)/d)*PolyLog(2, ℯ^(2*asinh((sqrt(e)*x)/sqrt(d)))))/(2*sqrt(e)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^3, x, 2, -(sqrt(-e)*sqrt(d + e*x^2))/(2*d*x) - atan((sqrt(-e)*x)/sqrt(d + e*x^2))/(2*x^2)]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^5, x, 3, -(sqrt(-e)*sqrt(d + e*x^2))/(12*d*x^3) - ((-e)^(3/2)*sqrt(d + e*x^2))/(6*d^2*x) - atan((sqrt(-e)*x)/sqrt(d + e*x^2))/(4*x^4)]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^7, x, 4, -(sqrt(-e)*sqrt(d + e*x^2))/(30*d*x^5) - (2*(-e)^(3/2)*sqrt(d + e*x^2))/(45*d^2*x^3) - (4*(-e)^(5/2)*sqrt(d + e*x^2))/(45*d^3*x) - atan((sqrt(-e)*x)/sqrt(d + e*x^2))/(6*x^6)]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^9, x, 5, -(sqrt(-e)*sqrt(d + e*x^2))/(56*d*x^7) - (3*(-e)^(3/2)*sqrt(d + e*x^2))/(140*d^2*x^5) - ((-e)^(5/2)*sqrt(d + e*x^2))/(35*d^3*x^3) - (2*(-e)^(7/2)*sqrt(d + e*x^2))/(35*d^4*x) - atan((sqrt(-e)*x)/sqrt(d + e*x^2))/(8*x^8)]

    @test_int [x^6*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 4, (d^3*sqrt(d + e*x^2))/(7*(-e)^(7/2)) - (d^2*(d + e*x^2)^(3/2))/(7*(-e)^(7/2)) + (3*d*(d + e*x^2)^(5/2))/(35*(-e)^(7/2)) - (d + e*x^2)^(7/2)/(49*(-e)^(7/2)) + (x^7*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/7]
    @test_int [x^4*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 4, (d^2*sqrt(d + e*x^2))/(5*(-e)^(5/2)) - (2*d*(d + e*x^2)^(3/2))/(15*(-e)^(5/2)) + (d + e*x^2)^(5/2)/(25*(-e)^(5/2)) + (x^5*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/5]
    @test_int [x^2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 4, (d*sqrt(d + e*x^2))/(3*(-e)^(3/2)) - (d + e*x^2)^(3/2)/(9*(-e)^(3/2)) + (x^3*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/3]
    @test_int [x^0*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 2, sqrt(d + e*x^2)/sqrt(-e) + x*atan((sqrt(-e)*x)/sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^2, x, 4, -(atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x) - (sqrt(-e)*atanh(sqrt(d + e*x^2)/sqrt(d)))/sqrt(d)]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^4, x, 5, -(sqrt(-e)*sqrt(d + e*x^2))/(6*d*x^2) - atan((sqrt(-e)*x)/sqrt(d + e*x^2))/(3*x^3) - ((-e)^(3/2)*atanh(sqrt(d + e*x^2)/sqrt(d)))/(6*d^(3/2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^6, x, 6, -(sqrt(-e)*sqrt(d + e*x^2))/(20*d*x^4) - (3*(-e)^(3/2)*sqrt(d + e*x^2))/(40*d^2*x^2) - atan((sqrt(-e)*x)/sqrt(d + e*x^2))/(5*x^5) - (3*(-e)^(5/2)*atanh(sqrt(d + e*x^2)/sqrt(d)))/(40*d^(5/2))]

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

    @test_int [x^(9/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 6, (60*d^2*sqrt(x)*sqrt(d + e*x^2))/(847*(-e)^(5/2)) + (36*d*x^(5/2)*sqrt(d + e*x^2))/(847*(-e)^(3/2)) + (4*x^(9/2)*sqrt(d + e*x^2))/(121*sqrt(-e)) + (2*x^(11/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/11 + (30*d^(11/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(847*e^(13/4)*sqrt(d + e*x^2))]
    @test_int [x^(5/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 5, (20*d*sqrt(x)*sqrt(d + e*x^2))/(147*(-e)^(3/2)) + (4*x^(5/2)*sqrt(d + e*x^2))/(49*sqrt(-e)) + (2*x^(7/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/7 - (10*d^(7/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(147*e^(9/4)*sqrt(d + e*x^2))]
    @test_int [x^(1/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 4, (4*sqrt(x)*sqrt(d + e*x^2))/(9*sqrt(-e)) + (2*x^(3/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/3 + (2*d^(3/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(9*e^(5/4)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^(3/2), x, 3, (-2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/sqrt(x) + (2*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(d^(1/4)*e^(1/4)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^(7/2), x, 4, (-4*sqrt(-e)*sqrt(d + e*x^2))/(15*d*x^(3/2)) - (2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/(5*x^(5/2)) - (2*sqrt(-e)*e^(3/4)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(15*d^(5/4)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^(11/2), x, 5, (-4*sqrt(-e)*sqrt(d + e*x^2))/(63*d*x^(7/2)) - (20*(-e)^(3/2)*sqrt(d + e*x^2))/(189*d^2*x^(3/2)) - (2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/(9*x^(9/2)) + (10*sqrt(-e)*e^(7/4)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(189*d^(9/4)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^(15/2), x, 6, (-4*sqrt(-e)*sqrt(d + e*x^2))/(143*d*x^(11/2)) - (36*(-e)^(3/2)*sqrt(d + e*x^2))/(1001*d^2*x^(7/2)) - (60*(-e)^(5/2)*sqrt(d + e*x^2))/(1001*d^3*x^(3/2)) - (2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/(13*x^(13/2)) - (30*sqrt(-e)*e^(11/4)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(1001*d^(13/4)*sqrt(d + e*x^2))]

    @test_int [x^(7/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 7, (28*d*x^(3/2)*sqrt(d + e*x^2))/(405*(-e)^(3/2)) + (4*x^(7/2)*sqrt(d + e*x^2))/(81*sqrt(-e)) - (28*d^2*sqrt(-e)*sqrt(x)*sqrt(d + e*x^2))/(135*e^(5/2)*(sqrt(d) + sqrt(e)*x)) + (2*x^(9/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/9 + (28*d^(9/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.E(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(135*e^(11/4)*sqrt(d + e*x^2)) - (14*d^(9/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(135*e^(11/4)*sqrt(d + e*x^2))]
    @test_int [x^(3/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)), x, 6, (4*x^(3/2)*sqrt(d + e*x^2))/(25*sqrt(-e)) + (12*d*sqrt(-e)*sqrt(x)*sqrt(d + e*x^2))/(25*e^(3/2)*(sqrt(d) + sqrt(e)*x)) + (2*x^(5/2)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/5 - (12*d^(5/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.E(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(25*e^(7/4)*sqrt(d + e*x^2)) + (6*d^(5/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(25*e^(7/4)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^(1/2), x, 5, (-4*sqrt(-e)*sqrt(x)*sqrt(d + e*x^2))/(sqrt(e)*(sqrt(d) + sqrt(e)*x)) + 2*sqrt(x)*atan((sqrt(-e)*x)/sqrt(d + e*x^2)) + (4*d^(1/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.E(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(e^(3/4)*sqrt(d + e*x^2)) - (2*d^(1/4)*sqrt(-e)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(e^(3/4)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^(5/2), x, 6, (-4*sqrt(-e)*sqrt(d + e*x^2))/(3*d*sqrt(x)) + (4*sqrt(-e^2)*sqrt(x)*sqrt(d + e*x^2))/(3*d*(sqrt(d) + sqrt(e)*x)) - (2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/(3*x^(3/2)) - (4*sqrt(-e)*e^(1/4)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.E(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(3*d^(3/4)*sqrt(d + e*x^2)) + (2*sqrt(-e)*e^(1/4)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(3*d^(3/4)*sqrt(d + e*x^2))]
    @test_int [atan((sqrt(-e)*x)/sqrt(d + e*x^2))/x^(9/2), x, 7, (-4*sqrt(-e)*sqrt(d + e*x^2))/(35*d*x^(5/2)) - (12*(-e)^(3/2)*sqrt(d + e*x^2))/(35*d^2*sqrt(x)) - (12*sqrt(-e)*e^(3/2)*sqrt(x)*sqrt(d + e*x^2))/(35*d^2*(sqrt(d) + sqrt(e)*x)) - (2*atan((sqrt(-e)*x)/sqrt(d + e*x^2)))/(7*x^(7/2)) + (12*sqrt(-e)*e^(5/4)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.E(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(35*d^(7/4)*sqrt(d + e*x^2)) - (6*sqrt(-e)*e^(5/4)*(sqrt(d) + sqrt(e)*x)*sqrt((d + e*x^2)/(sqrt(d) + sqrt(e)*x)^2)*Elliptic.F(2*atan((e^(1/4)*sqrt(x))/d^(1/4)), 1/2))/(35*d^(7/4)*sqrt(d + e*x^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*u*atan(a+b*x+c*x^2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(a+b*x+c*x^2)=#

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

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*u^m*(a+b*atan(sqrt(1-c*x)/sqrt(1+c*x)))^n=#

    @test_int [(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^n/(1 - c^2*x^2), x, 0, Unintegrable((a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^n/(1 - c^2*x^2), x)]

    @test_int [(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^3/(1 - c^2*x^2), x, 9, -((2*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^3*atanh(1 - 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/c) + (3*im*b*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^2*PolyLog(2, 1 - 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(2*c) - (3*im*b*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^2*PolyLog(2, -1 + 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(2*c) + (3*b^2*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))*PolyLog(3, 1 - 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(2*c) - (3*b^2*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))*PolyLog(3, -1 + 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(2*c) - (3*im*b^3*PolyLog(4, 1 - 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(4*c) + (3*im*b^3*PolyLog(4, -1 + 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(4*c)]
    @test_int [(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^2/(1 - c^2*x^2), x, 7, -((2*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^2*atanh(1 - 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/c) + (im*b*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))*PolyLog(2, 1 - 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/c - (im*b*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))*PolyLog(2, -1 + 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/c + (b^2*PolyLog(3, 1 - 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(2*c) - (b^2*PolyLog(3, -1 + 2/(1 + (im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(2*c)]
    @test_int [(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^1/(1 - c^2*x^2), x, 4, -((a*log(sqrt(1 - c*x)/sqrt(1 + c*x)))/c) - (im*b*PolyLog(2, -((im*sqrt(1 - c*x))/sqrt(1 + c*x))))/(2*c) + (im*b*PolyLog(2, (im*sqrt(1 - c*x))/sqrt(1 + c*x)))/(2*c)]
    @test_int [1/((a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^1*(1 - c^2*x^2)), x, 0, Unintegrable(1/((1 - c^2*x^2)*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))), x)]
    @test_int [1/((a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^2*(1 - c^2*x^2)), x, 0, Unintegrable(1/((1 - c^2*x^2)*(a + b*atan(sqrt(1 - c*x)/sqrt(1 + c*x)))^2), x)]

    #= ::Title::Closed:: =#
    #=Integrands*involving*inverse*tangents*of*trig*functions=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(Trig(a+b*x))=#

    @test_int [x^m*atan(tan(a + b*x)), x, 2, -((b*x^(2 + m))/(2 + 3*m + m^2)) + (x^(1 + m)*atan(tan(a + b*x)))/(1 + m)]

    @test_int [x^2*atan(tan(a + b*x)), x, 2, -((b*x^4)/12) + (1/3)*x^3*atan(tan(a + b*x))]
    @test_int [x^1*atan(tan(a + b*x)), x, 2, -((b*x^3)/6) + (1/2)*x^2*atan(tan(a + b*x))]
    @test_int [x^0*atan(tan(a + b*x)), x, 2, atan(tan(a + b*x))^2/(2*b)]
    @test_int [atan(tan(a + b*x))/x^1, x, 2, b*x - (b*x - atan(tan(a + b*x)))*log(x)]

    @test_int [x^m*atan(cot(a + b*x)), x, 2, (b*x^(2 + m))/(2 + 3*m + m^2) + (x^(1 + m)*atan(cot(a + b*x)))/(1 + m)]

    @test_int [x^2*atan(cot(a + b*x)), x, 2, (b*x^4)/12 + (1/3)*x^3*atan(cot(a + b*x))]
    @test_int [x^1*atan(cot(a + b*x)), x, 2, (b*x^3)/6 + (1/2)*x^2*atan(cot(a + b*x))]
    @test_int [x^0*atan(cot(a + b*x)), x, 2, -(atan(cot(a + b*x))^2/(2*b))]
    @test_int [atan(cot(a + b*x))/x^1, x, 2, (-b)*x + (b*x + atan(cot(a + b*x)))*log(x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(c+d*Trig(a+b*x))=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(c+d*tan(a+b*x))=#

    @test_int [atan(tan(a + b*x)), x, 2, atan(tan(a + b*x))^2/(2*b)]

    @test_int [x^2*atan(c + d*tan(a + b*x)), x, 11, (1/3)*x^3*atan(c + d*tan(a + b*x)) + (1/6)*im*x^3*log(1 + ((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d)) - (1/6)*im*x^3*log(1 + ((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d))) + (x^2*PolyLog(2, -(((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d))))/(4*b) - (x^2*PolyLog(2, -(((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d)))))/(4*b) + (im*x*PolyLog(3, -(((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d))))/(4*b^2) - (im*x*PolyLog(3, -(((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d)))))/(4*b^2) - PolyLog(4, -(((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d)))/(8*b^3) + PolyLog(4, -(((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d))))/(8*b^3)]
    @test_int [x^1*atan(c + d*tan(a + b*x)), x, 9, (1/2)*x^2*atan(c + d*tan(a + b*x)) + (1/4)*im*x^2*log(1 + ((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d)) - (1/4)*im*x^2*log(1 + ((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d))) + (x*PolyLog(2, -(((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d))))/(4*b) - (x*PolyLog(2, -(((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d)))))/(4*b) + (im*PolyLog(3, -(((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d))))/(8*b^2) - (im*PolyLog(3, -(((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d)))))/(8*b^2)]
    @test_int [x^0*atan(c + d*tan(a + b*x)), x, 7, x*atan(c + d*tan(a + b*x)) + (1/2)*im*x*log(1 + ((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d)) - (1/2)*im*x*log(1 + ((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d))) + PolyLog(2, -(((1 + im*c + d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c - d)))/(4*b) - PolyLog(2, -(((c + im*(1 - d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 + d))))/(4*b)]
    @test_int [atan(c + d*tan(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + d*tan(a + b*x))/x, x)]

    @test_int [x^2*atan(c + (1 + im*c)*tan(a + b*x)), x, 7, -((b*x^4)/12) + (1/3)*x^3*atan(c + (1 + im*c)*tan(a + b*x)) - (1/6)*im*x^3*log(1 - im*c*ℯ^(2*im*a + 2*im*b*x)) - (x^2*PolyLog(2, im*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) - (im*x*PolyLog(3, im*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b^2) + PolyLog(4, im*c*ℯ^(2*im*a + 2*im*b*x))/(8*b^3)]
    @test_int [x^1*atan(c + (1 + im*c)*tan(a + b*x)), x, 6, -((b*x^3)/6) + (1/2)*x^2*atan(c + (1 + im*c)*tan(a + b*x)) - (1/4)*im*x^2*log(1 - im*c*ℯ^(2*im*a + 2*im*b*x)) - (x*PolyLog(2, im*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) - (im*PolyLog(3, im*c*ℯ^(2*im*a + 2*im*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c + (1 + im*c)*tan(a + b*x)), x, 5, -((b*x^2)/2) + x*atan(c + (1 + im*c)*tan(a + b*x)) - (1/2)*im*x*log(1 - im*c*ℯ^(2*im*a + 2*im*b*x)) - PolyLog(2, im*c*ℯ^(2*im*a + 2*im*b*x))/(4*b)]
    @test_int [atan(c + (1 + im*c)*tan(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + (1 + im*c)*tan(a + b*x))/x, x)]

    @test_int [x^2*atan(c + (-1 + im*c)*tan(a + b*x)), x, 7, (b*x^4)/12 + (1/3)*x^3*atan(c - (1 - im*c)*tan(a + b*x)) + (1/6)*im*x^3*log(1 + im*c*ℯ^(2*im*a + 2*im*b*x)) + (x^2*PolyLog(2, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) + (im*x*PolyLog(3, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b^2) - PolyLog(4, (-im)*c*ℯ^(2*im*a + 2*im*b*x))/(8*b^3)]
    @test_int [x^1*atan(c + (-1 + im*c)*tan(a + b*x)), x, 6, (b*x^3)/6 + (1/2)*x^2*atan(c - (1 - im*c)*tan(a + b*x)) + (1/4)*im*x^2*log(1 + im*c*ℯ^(2*im*a + 2*im*b*x)) + (x*PolyLog(2, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) + (im*PolyLog(3, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c + (-1 + im*c)*tan(a + b*x)), x, 5, (b*x^2)/2 + x*atan(c - (1 - im*c)*tan(a + b*x)) + (1/2)*im*x*log(1 + im*c*ℯ^(2*im*a + 2*im*b*x)) + PolyLog(2, (-im)*c*ℯ^(2*im*a + 2*im*b*x))/(4*b)]
    @test_int [atan(c + (-1 + im*c)*tan(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + (-1 + im*c)*tan(a + b*x))/x, x)]

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

    @test_int [atan(cot(a + b*x)), x, 2, -(atan(cot(a + b*x))^2/(2*b))]

    @test_int [x^2*atan(c + d*cot(a + b*x)), x, 11, (1/3)*x^3*atan(c + d*cot(a + b*x)) + (1/6)*im*x^3*log(1 - ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d)) - (1/6)*im*x^3*log(1 - ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d))) + (x^2*PolyLog(2, ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d)))/(4*b) - (x^2*PolyLog(2, ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d))))/(4*b) + (im*x*PolyLog(3, ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d)))/(4*b^2) - (im*x*PolyLog(3, ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d))))/(4*b^2) - PolyLog(4, ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d))/(8*b^3) + PolyLog(4, ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d)))/(8*b^3)]
    @test_int [x^1*atan(c + d*cot(a + b*x)), x, 9, (1/2)*x^2*atan(c + d*cot(a + b*x)) + (1/4)*im*x^2*log(1 - ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d)) - (1/4)*im*x^2*log(1 - ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d))) + (x*PolyLog(2, ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d)))/(4*b) - (x*PolyLog(2, ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d))))/(4*b) + (im*PolyLog(3, ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d)))/(8*b^2) - (im*PolyLog(3, ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d))))/(8*b^2)]
    @test_int [x^0*atan(c + d*cot(a + b*x)), x, 7, x*atan(c + d*cot(a + b*x)) + (1/2)*im*x*log(1 - ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d)) - (1/2)*im*x*log(1 - ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d))) + PolyLog(2, ((1 + im*c - d)*ℯ^(2*im*a + 2*im*b*x))/(1 + im*c + d))/(4*b) - PolyLog(2, ((c + im*(1 + d))*ℯ^(2*im*a + 2*im*b*x))/(c + im*(1 - d)))/(4*b)]
    @test_int [atan(c + d*cot(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + d*cot(a + b*x))/x, x)]

    @test_int [x^2*atan(c + (1 - im*c)*cot(a + b*x)), x, 7, (b*x^4)/12 + (1/3)*x^3*atan(c + (1 - im*c)*cot(a + b*x)) + (1/6)*im*x^3*log(1 - im*c*ℯ^(2*im*a + 2*im*b*x)) + (x^2*PolyLog(2, im*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) + (im*x*PolyLog(3, im*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b^2) - PolyLog(4, im*c*ℯ^(2*im*a + 2*im*b*x))/(8*b^3)]
    @test_int [x^1*atan(c + (1 - im*c)*cot(a + b*x)), x, 6, (b*x^3)/6 + (1/2)*x^2*atan(c + (1 - im*c)*cot(a + b*x)) + (1/4)*im*x^2*log(1 - im*c*ℯ^(2*im*a + 2*im*b*x)) + (x*PolyLog(2, im*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) + (im*PolyLog(3, im*c*ℯ^(2*im*a + 2*im*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c + (1 - im*c)*cot(a + b*x)), x, 5, (b*x^2)/2 + x*atan(c + (1 - im*c)*cot(a + b*x)) + (1/2)*im*x*log(1 - im*c*ℯ^(2*im*a + 2*im*b*x)) + PolyLog(2, im*c*ℯ^(2*im*a + 2*im*b*x))/(4*b)]
    @test_int [atan(c + (1 - im*c)*cot(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + (1 - im*c)*cot(a + b*x))/x, x)]

    @test_int [x^2*atan(c + (-1 - im*c)*cot(a + b*x)), x, 7, -((b*x^4)/12) + (1/3)*x^3*atan(c - (1 + im*c)*cot(a + b*x)) - (1/6)*im*x^3*log(1 + im*c*ℯ^(2*im*a + 2*im*b*x)) - (x^2*PolyLog(2, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) - (im*x*PolyLog(3, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b^2) + PolyLog(4, (-im)*c*ℯ^(2*im*a + 2*im*b*x))/(8*b^3)]
    @test_int [x^1*atan(c + (-1 - im*c)*cot(a + b*x)), x, 6, -((b*x^3)/6) + (1/2)*x^2*atan(c - (1 + im*c)*cot(a + b*x)) - (1/4)*im*x^2*log(1 + im*c*ℯ^(2*im*a + 2*im*b*x)) - (x*PolyLog(2, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(4*b) - (im*PolyLog(3, (-im)*c*ℯ^(2*im*a + 2*im*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c + (-1 - im*c)*cot(a + b*x)), x, 5, -((b*x^2)/2) + x*atan(c - (1 + im*c)*cot(a + b*x)) - (1/2)*im*x*log(1 + im*c*ℯ^(2*im*a + 2*im*b*x)) - PolyLog(2, (-im)*c*ℯ^(2*im*a + 2*im*b*x))/(4*b)]
    @test_int [atan(c + (-1 - im*c)*cot(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + (-1 - im*c)*cot(a + b*x))/x, x)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(c+d*Hyper(a+b*x))=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(c+d*sinh(a+b*x))=#

    @test_int [atan(sinh(x)), x, 6, -2*x*atan(ℯ^x) + x*atan(sinh(x)) + im*PolyLog(2, (-im)*ℯ^x) - im*PolyLog(2, im*ℯ^x)]
    @test_int [x*atan(sinh(x)), x, 8, (-x^2)*atan(ℯ^x) + (1/2)*x^2*atan(sinh(x)) + im*x*PolyLog(2, (-im)*ℯ^x) - im*x*PolyLog(2, im*ℯ^x) - im*PolyLog(3, (-im)*ℯ^x) + im*PolyLog(3, im*ℯ^x)]
    @test_int [x^2*atan(sinh(x)), x, 10, (-(2/3))*x^3*atan(ℯ^x) + (1/3)*x^3*atan(sinh(x)) + im*x^2*PolyLog(2, (-im)*ℯ^x) - im*x^2*PolyLog(2, im*ℯ^x) - 2*im*x*PolyLog(3, (-im)*ℯ^x) + 2*im*x*PolyLog(3, im*ℯ^x) + 2*im*PolyLog(4, (-im)*ℯ^x) - 2*im*PolyLog(4, im*ℯ^x)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*x^m*atan(c+d*cosh(a+b*x))=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(c+d*tanh(a+b*x))=#

    @test_int [(e + f*x)^3*atan(tanh(a + b*x)), x, 12, -(((e + f*x)^4*atan(ℯ^(2*a + 2*b*x)))/(4*f)) + ((e + f*x)^4*atan(tanh(a + b*x)))/(4*f) + (im*(e + f*x)^3*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) - (im*(e + f*x)^3*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b) - (3*im*f*(e + f*x)^2*PolyLog(3, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^2) + (3*im*f*(e + f*x)^2*PolyLog(3, im*ℯ^(2*a + 2*b*x)))/(8*b^2) + (3*im*f^2*(e + f*x)*PolyLog(4, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^3) - (3*im*f^2*(e + f*x)*PolyLog(4, im*ℯ^(2*a + 2*b*x)))/(8*b^3) - (3*im*f^3*PolyLog(5, (-im)*ℯ^(2*a + 2*b*x)))/(16*b^4) + (3*im*f^3*PolyLog(5, im*ℯ^(2*a + 2*b*x)))/(16*b^4)]
    @test_int [(e + f*x)^2*atan(tanh(a + b*x)), x, 10, -(((e + f*x)^3*atan(ℯ^(2*a + 2*b*x)))/(3*f)) + ((e + f*x)^3*atan(tanh(a + b*x)))/(3*f) + (im*(e + f*x)^2*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) - (im*(e + f*x)^2*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b) - (im*f*(e + f*x)*PolyLog(3, (-im)*ℯ^(2*a + 2*b*x)))/(4*b^2) + (im*f*(e + f*x)*PolyLog(3, im*ℯ^(2*a + 2*b*x)))/(4*b^2) + (im*f^2*PolyLog(4, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^3) - (im*f^2*PolyLog(4, im*ℯ^(2*a + 2*b*x)))/(8*b^3)]
    @test_int [(e + f*x)^1*atan(tanh(a + b*x)), x, 8, -(((e + f*x)^2*atan(ℯ^(2*a + 2*b*x)))/(2*f)) + ((e + f*x)^2*atan(tanh(a + b*x)))/(2*f) + (im*(e + f*x)*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) - (im*(e + f*x)*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b) - (im*f*PolyLog(3, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^2) + (im*f*PolyLog(3, im*ℯ^(2*a + 2*b*x)))/(8*b^2)]
    @test_int [(e + f*x)^0*atan(tanh(a + b*x)), x, 6, (-x)*atan(ℯ^(2*a + 2*b*x)) + x*atan(tanh(a + b*x)) + (im*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) - (im*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b)]
    @test_int [atan(tanh(a + b*x))/(e + f*x)^1, x, 0, CannotIntegrate(atan(tanh(a + b*x))/(e + f*x), x)]

    @test_int [x^2*atan(c + d*tanh(a + b*x)), x, 11, (1/3)*x^3*atan(c + d*tanh(a + b*x)) + (1/6)*im*x^3*log(1 + ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)) - (1/6)*im*x^3*log(1 + ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)) + (im*x^2*PolyLog(2, -(((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d))))/(4*b) - (im*x^2*PolyLog(2, -(((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d))))/(4*b) - (im*x*PolyLog(3, -(((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d))))/(4*b^2) + (im*x*PolyLog(3, -(((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d))))/(4*b^2) + (im*PolyLog(4, -(((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d))))/(8*b^3) - (im*PolyLog(4, -(((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d))))/(8*b^3)]
    @test_int [x^1*atan(c + d*tanh(a + b*x)), x, 9, (1/2)*x^2*atan(c + d*tanh(a + b*x)) + (1/4)*im*x^2*log(1 + ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)) - (1/4)*im*x^2*log(1 + ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)) + (im*x*PolyLog(2, -(((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d))))/(4*b) - (im*x*PolyLog(2, -(((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d))))/(4*b) - (im*PolyLog(3, -(((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d))))/(8*b^2) + (im*PolyLog(3, -(((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d))))/(8*b^2)]
    @test_int [x^0*atan(c + d*tanh(a + b*x)), x, 7, x*atan(c + d*tanh(a + b*x)) + (1/2)*im*x*log(1 + ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)) - (1/2)*im*x*log(1 + ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)) + (im*PolyLog(2, -(((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d))))/(4*b) - (im*PolyLog(2, -(((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d))))/(4*b)]
    @test_int [atan(c + d*tanh(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + d*tanh(a + b*x))/x, x)]

    @test_int [x^2*atan(c + (im + c)*tanh(a + b*x)), x, 7, (-(1/12))*im*b*x^4 + (1/3)*x^3*atan(c + (im + c)*tanh(a + b*x)) + (1/6)*im*x^3*log(1 + im*c*ℯ^(2*a + 2*b*x)) + (im*x^2*PolyLog(2, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b) - (im*x*PolyLog(3, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b^2) + (im*PolyLog(4, (-im)*c*ℯ^(2*a + 2*b*x)))/(8*b^3)]
    @test_int [x^1*atan(c + (im + c)*tanh(a + b*x)), x, 6, (-(1/6))*im*b*x^3 + (1/2)*x^2*atan(c + (im + c)*tanh(a + b*x)) + (1/4)*im*x^2*log(1 + im*c*ℯ^(2*a + 2*b*x)) + (im*x*PolyLog(2, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b) - (im*PolyLog(3, (-im)*c*ℯ^(2*a + 2*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c + (im + c)*tanh(a + b*x)), x, 5, (-(1/2))*im*b*x^2 + x*atan(c + (im + c)*tanh(a + b*x)) + (1/2)*im*x*log(1 + im*c*ℯ^(2*a + 2*b*x)) + (im*PolyLog(2, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b)]
    @test_int [atan(c + (im + c)*tanh(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + (im + c)*tanh(a + b*x))/x, x)]

    @test_int [x^2*atan(c - (im - c)*tanh(a + b*x)), x, 7, (1/12)*im*b*x^4 + (1/3)*x^3*atan(c - (im - c)*tanh(a + b*x)) - (1/6)*im*x^3*log(1 - im*c*ℯ^(2*a + 2*b*x)) - (im*x^2*PolyLog(2, im*c*ℯ^(2*a + 2*b*x)))/(4*b) + (im*x*PolyLog(3, im*c*ℯ^(2*a + 2*b*x)))/(4*b^2) - (im*PolyLog(4, im*c*ℯ^(2*a + 2*b*x)))/(8*b^3)]
    @test_int [x^1*atan(c - (im - c)*tanh(a + b*x)), x, 6, (1/6)*im*b*x^3 + (1/2)*x^2*atan(c - (im - c)*tanh(a + b*x)) - (1/4)*im*x^2*log(1 - im*c*ℯ^(2*a + 2*b*x)) - (im*x*PolyLog(2, im*c*ℯ^(2*a + 2*b*x)))/(4*b) + (im*PolyLog(3, im*c*ℯ^(2*a + 2*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c - (im - c)*tanh(a + b*x)), x, 5, (1/2)*im*b*x^2 + x*atan(c - (im - c)*tanh(a + b*x)) - (1/2)*im*x*log(1 - im*c*ℯ^(2*a + 2*b*x)) - (im*PolyLog(2, im*c*ℯ^(2*a + 2*b*x)))/(4*b)]
    @test_int [atan(c - (im - c)*tanh(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c - (im - c)*tanh(a + b*x))/x, x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*atan(c+d*coth(a+b*x))=#

    @test_int [(e + f*x)^3*atan(coth(a + b*x)), x, 12, ((e + f*x)^4*atan(ℯ^(2*a + 2*b*x)))/(4*f) + ((e + f*x)^4*atan(coth(a + b*x)))/(4*f) - (im*(e + f*x)^3*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) + (im*(e + f*x)^3*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b) + (3*im*f*(e + f*x)^2*PolyLog(3, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^2) - (3*im*f*(e + f*x)^2*PolyLog(3, im*ℯ^(2*a + 2*b*x)))/(8*b^2) - (3*im*f^2*(e + f*x)*PolyLog(4, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^3) + (3*im*f^2*(e + f*x)*PolyLog(4, im*ℯ^(2*a + 2*b*x)))/(8*b^3) + (3*im*f^3*PolyLog(5, (-im)*ℯ^(2*a + 2*b*x)))/(16*b^4) - (3*im*f^3*PolyLog(5, im*ℯ^(2*a + 2*b*x)))/(16*b^4)]
    @test_int [(e + f*x)^2*atan(coth(a + b*x)), x, 10, ((e + f*x)^3*atan(ℯ^(2*a + 2*b*x)))/(3*f) + ((e + f*x)^3*atan(coth(a + b*x)))/(3*f) - (im*(e + f*x)^2*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) + (im*(e + f*x)^2*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b) + (im*f*(e + f*x)*PolyLog(3, (-im)*ℯ^(2*a + 2*b*x)))/(4*b^2) - (im*f*(e + f*x)*PolyLog(3, im*ℯ^(2*a + 2*b*x)))/(4*b^2) - (im*f^2*PolyLog(4, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^3) + (im*f^2*PolyLog(4, im*ℯ^(2*a + 2*b*x)))/(8*b^3)]
    @test_int [(e + f*x)^1*atan(coth(a + b*x)), x, 8, ((e + f*x)^2*atan(ℯ^(2*a + 2*b*x)))/(2*f) + ((e + f*x)^2*atan(coth(a + b*x)))/(2*f) - (im*(e + f*x)*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) + (im*(e + f*x)*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b) + (im*f*PolyLog(3, (-im)*ℯ^(2*a + 2*b*x)))/(8*b^2) - (im*f*PolyLog(3, im*ℯ^(2*a + 2*b*x)))/(8*b^2)]
    @test_int [(e + f*x)^0*atan(coth(a + b*x)), x, 6, x*atan(ℯ^(2*a + 2*b*x)) + x*atan(coth(a + b*x)) - (im*PolyLog(2, (-im)*ℯ^(2*a + 2*b*x)))/(4*b) + (im*PolyLog(2, im*ℯ^(2*a + 2*b*x)))/(4*b)]
    @test_int [atan(coth(a + b*x))/(e + f*x)^1, x, 0, CannotIntegrate(atan(coth(a + b*x))/(e + f*x), x)]

    @test_int [x^2*atan(c + d*coth(a + b*x)), x, 11, (1/3)*x^3*atan(c + d*coth(a + b*x)) + (1/6)*im*x^3*log(1 - ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)) - (1/6)*im*x^3*log(1 - ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)) + (im*x^2*PolyLog(2, ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)))/(4*b) - (im*x^2*PolyLog(2, ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)))/(4*b) - (im*x*PolyLog(3, ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)))/(4*b^2) + (im*x*PolyLog(3, ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)))/(4*b^2) + (im*PolyLog(4, ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)))/(8*b^3) - (im*PolyLog(4, ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)))/(8*b^3)]
    @test_int [x^1*atan(c + d*coth(a + b*x)), x, 9, (1/2)*x^2*atan(c + d*coth(a + b*x)) + (1/4)*im*x^2*log(1 - ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)) - (1/4)*im*x^2*log(1 - ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)) + (im*x*PolyLog(2, ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)))/(4*b) - (im*x*PolyLog(2, ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)))/(4*b) - (im*PolyLog(3, ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)))/(8*b^2) + (im*PolyLog(3, ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)))/(8*b^2)]
    @test_int [x^0*atan(c + d*coth(a + b*x)), x, 7, x*atan(c + d*coth(a + b*x)) + (1/2)*im*x*log(1 - ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)) - (1/2)*im*x*log(1 - ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)) + (im*PolyLog(2, ((im - c - d)*ℯ^(2*a + 2*b*x))/(im - c + d)))/(4*b) - (im*PolyLog(2, ((im + c + d)*ℯ^(2*a + 2*b*x))/(im + c - d)))/(4*b)]
    @test_int [atan(c + d*coth(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + d*coth(a + b*x))/x, x)]

    @test_int [x^2*atan(c + (im + c)*coth(a + b*x)), x, 7, (-(1/12))*im*b*x^4 + (1/3)*x^3*atan(c + (im + c)*coth(a + b*x)) + (1/6)*im*x^3*log(1 - im*c*ℯ^(2*a + 2*b*x)) + (im*x^2*PolyLog(2, im*c*ℯ^(2*a + 2*b*x)))/(4*b) - (im*x*PolyLog(3, im*c*ℯ^(2*a + 2*b*x)))/(4*b^2) + (im*PolyLog(4, im*c*ℯ^(2*a + 2*b*x)))/(8*b^3)]
    @test_int [x^1*atan(c + (im + c)*coth(a + b*x)), x, 6, (-(1/6))*im*b*x^3 + (1/2)*x^2*atan(c + (im + c)*coth(a + b*x)) + (1/4)*im*x^2*log(1 - im*c*ℯ^(2*a + 2*b*x)) + (im*x*PolyLog(2, im*c*ℯ^(2*a + 2*b*x)))/(4*b) - (im*PolyLog(3, im*c*ℯ^(2*a + 2*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c + (im + c)*coth(a + b*x)), x, 5, (-(1/2))*im*b*x^2 + x*atan(c + (im + c)*coth(a + b*x)) + (1/2)*im*x*log(1 - im*c*ℯ^(2*a + 2*b*x)) + (im*PolyLog(2, im*c*ℯ^(2*a + 2*b*x)))/(4*b)]
    @test_int [atan(c + (im + c)*coth(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c + (im + c)*coth(a + b*x))/x, x)]

    @test_int [x^2*atan(c - (im - c)*coth(a + b*x)), x, 7, (1/12)*im*b*x^4 + (1/3)*x^3*atan(c - (im - c)*coth(a + b*x)) - (1/6)*im*x^3*log(1 + im*c*ℯ^(2*a + 2*b*x)) - (im*x^2*PolyLog(2, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b) + (im*x*PolyLog(3, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b^2) - (im*PolyLog(4, (-im)*c*ℯ^(2*a + 2*b*x)))/(8*b^3)]
    @test_int [x^1*atan(c - (im - c)*coth(a + b*x)), x, 6, (1/6)*im*b*x^3 + (1/2)*x^2*atan(c - (im - c)*coth(a + b*x)) - (1/4)*im*x^2*log(1 + im*c*ℯ^(2*a + 2*b*x)) - (im*x*PolyLog(2, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b) + (im*PolyLog(3, (-im)*c*ℯ^(2*a + 2*b*x)))/(8*b^2)]
    @test_int [x^0*atan(c - (im - c)*coth(a + b*x)), x, 5, (1/2)*im*b*x^2 + x*atan(c - (im - c)*coth(a + b*x)) - (1/2)*im*x*log(1 + im*c*ℯ^(2*a + 2*b*x)) - (im*PolyLog(2, (-im)*c*ℯ^(2*a + 2*b*x)))/(4*b)]
    @test_int [atan(c - (im - c)*coth(a + b*x))/x^1, x, 0, CannotIntegrate(atan(c - (im - c)*coth(a + b*x))/x, x)]

    #= ::Title::Closed:: =#
    #=Integrands*involving*inverse*tangents*of*exponentials=#

    @test_int [atan(ℯ^x), x, 4, (1/2)*im*PolyLog(2, (-im)*ℯ^x) - (1/2)*im*PolyLog(2, im*ℯ^x)]
    @test_int [x*atan(ℯ^x), x, 7, (1/2)*im*x*PolyLog(2, (-im)*ℯ^x) - (1/2)*im*x*PolyLog(2, im*ℯ^x) - (1/2)*im*PolyLog(3, (-im)*ℯ^x) + (1/2)*im*PolyLog(3, im*ℯ^x)]
    @test_int [x^2*atan(ℯ^x), x, 9, (1/2)*im*x^2*PolyLog(2, (-im)*ℯ^x) - (1/2)*im*x^2*PolyLog(2, im*ℯ^x) - im*x*PolyLog(3, (-im)*ℯ^x) + im*x*PolyLog(3, im*ℯ^x) + im*PolyLog(4, (-im)*ℯ^x) - im*PolyLog(4, im*ℯ^x)]

    @test_int [atan(ℯ^(a + b*x)), x, 4, (im*PolyLog(2, (-im)*ℯ^(a + b*x)))/(2*b) - (im*PolyLog(2, im*ℯ^(a + b*x)))/(2*b)]
    @test_int [x*atan(ℯ^(a + b*x)), x, 7, (im*x*PolyLog(2, (-im)*ℯ^(a + b*x)))/(2*b) - (im*x*PolyLog(2, im*ℯ^(a + b*x)))/(2*b) - (im*PolyLog(3, (-im)*ℯ^(a + b*x)))/(2*b^2) + (im*PolyLog(3, im*ℯ^(a + b*x)))/(2*b^2)]
    @test_int [x^2*atan(ℯ^(a + b*x)), x, 9, (im*x^2*PolyLog(2, (-im)*ℯ^(a + b*x)))/(2*b) - (im*x^2*PolyLog(2, im*ℯ^(a + b*x)))/(2*b) - (im*x*PolyLog(3, (-im)*ℯ^(a + b*x)))/b^2 + (im*x*PolyLog(3, im*ℯ^(a + b*x)))/b^2 + (im*PolyLog(4, (-im)*ℯ^(a + b*x)))/b^3 - (im*PolyLog(4, im*ℯ^(a + b*x)))/b^3]

    @test_int [atan(a + b*f^(c + d*x)), x, 6, -((atan(a + b*f^(c + d*x))*log(2/(1 - im*(a + b*f^(c + d*x)))))/(d*log(f))) + (atan(a + b*f^(c + d*x))*log((2*b*f^(c + d*x))/((im - a)*(1 - im*(a + b*f^(c + d*x))))))/(d*log(f)) + (im*PolyLog(2, 1 - 2/(1 - im*(a + b*f^(c + d*x)))))/(2*d*log(f)) - (im*PolyLog(2, 1 - (2*b*f^(c + d*x))/((im - a)*(1 - im*(a + b*f^(c + d*x))))))/(2*d*log(f))]
    @test_int [x*atan(a + b*f^(c + d*x)), x, 9, (1/2)*x^2*atan(a + b*f^(c + d*x)) - (1/4)*im*x^2*log(1 - (im*b*f^(c + d*x))/(1 - im*a)) + (1/4)*im*x^2*log(1 + (im*b*f^(c + d*x))/(1 + im*a)) - (im*x*PolyLog(2, (im*b*f^(c + d*x))/(1 - im*a)))/(2*d*log(f)) + (im*x*PolyLog(2, -((im*b*f^(c + d*x))/(1 + im*a))))/(2*d*log(f)) + (im*PolyLog(3, (im*b*f^(c + d*x))/(1 - im*a)))/(2*d^2*log(f)^2) - (im*PolyLog(3, -((im*b*f^(c + d*x))/(1 + im*a))))/(2*d^2*log(f)^2), (1/4)*im*x^2*log(1 - im*a - im*b*f^(c + d*x)) - (1/4)*im*x^2*log(1 + im*a + im*b*f^(c + d*x)) + (1/4)*im*x^2*log(1 - (b*f^(c + d*x))/(im - a)) - (1/4)*im*x^2*log(1 + (b*f^(c + d*x))/(im + a)) + (im*x*PolyLog(2, (b*f^(c + d*x))/(im - a)))/(2*d*log(f)) - (im*x*PolyLog(2, -((b*f^(c + d*x))/(im + a))))/(2*d*log(f)) - (im*PolyLog(3, (b*f^(c + d*x))/(im - a)))/(2*d^2*log(f)^2) + (im*PolyLog(3, -((b*f^(c + d*x))/(im + a))))/(2*d^2*log(f)^2)]
    @test_int [x^2*atan(a + b*f^(c + d*x)), x, 11, (1/3)*x^3*atan(a + b*f^(c + d*x)) - (1/6)*im*x^3*log(1 - (im*b*f^(c + d*x))/(1 - im*a)) + (1/6)*im*x^3*log(1 + (im*b*f^(c + d*x))/(1 + im*a)) - (im*x^2*PolyLog(2, (im*b*f^(c + d*x))/(1 - im*a)))/(2*d*log(f)) + (im*x^2*PolyLog(2, -((im*b*f^(c + d*x))/(1 + im*a))))/(2*d*log(f)) + (im*x*PolyLog(3, (im*b*f^(c + d*x))/(1 - im*a)))/(d^2*log(f)^2) - (im*x*PolyLog(3, -((im*b*f^(c + d*x))/(1 + im*a))))/(d^2*log(f)^2) - (im*PolyLog(4, (im*b*f^(c + d*x))/(1 - im*a)))/(d^3*log(f)^3) + (im*PolyLog(4, -((im*b*f^(c + d*x))/(1 + im*a))))/(d^3*log(f)^3), (1/6)*im*x^3*log(1 - im*a - im*b*f^(c + d*x)) - (1/6)*im*x^3*log(1 + im*a + im*b*f^(c + d*x)) + (1/6)*im*x^3*log(1 - (b*f^(c + d*x))/(im - a)) - (1/6)*im*x^3*log(1 + (b*f^(c + d*x))/(im + a)) + (im*x^2*PolyLog(2, (b*f^(c + d*x))/(im - a)))/(2*d*log(f)) - (im*x^2*PolyLog(2, -((b*f^(c + d*x))/(im + a))))/(2*d*log(f)) - (im*x*PolyLog(3, (b*f^(c + d*x))/(im - a)))/(d^2*log(f)^2) + (im*x*PolyLog(3, -((b*f^(c + d*x))/(im + a))))/(d^2*log(f)^2) + (im*PolyLog(4, (b*f^(c + d*x))/(im - a)))/(d^3*log(f)^3) - (im*PolyLog(4, -((b*f^(c + d*x))/(im + a))))/(d^3*log(f)^3)]

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

    #= ::Title::Closed:: =#
    #=Miscellaneous*integrands*involving*inverse*tangents=#

    #= ::Section::Closed:: =#
    #=Problems*from*Calculus*textbooks=#

    #= ::Subsubsection::Closed:: =#
    #=Edwards*and*Penney*Calculus=#

    @test_int [atan(x)/(-1 + x)^3, x, 5, 1/(4*(1 - x)) - atan(x)/(2*(1 - x)^2) - (1/4)*log(1 - x) + (1/8)*log(1 + x^2)]
    @test_int [atan(1 + 2*x)/(4 + 3*x)^3, x, 9, -(1/(34*(4 + 3*x))) + (8/867)*atan(1 + 2*x) - atan(1 + 2*x)/(6*(4 + 3*x)^2) + (5/289)*log(4 + 3*x) - (5/578)*log(1 + 2*x + 2*x^2)]

    #= ::Subsubsection::Closed:: =#
    #=Thomas*Calculus, 8th*Edition=#

    @test_int [atan(sqrt(1 + x)), x, 4, -sqrt(1 + x) + 2*atan(sqrt(1 + x)) + x*atan(sqrt(1 + x))]

    #= ::Section::Closed:: =#
    #=Miscellaneous*integrands*involving*inverse*tangents=#

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

    @test_int [(x + x^3 + (1 + x)^2*atan(x))/((1 + x)^2*(1 + x^2)), x, 5, 1/(1 + x) + atan(x)^2/2 + log(1 + x)]

    @test_int [x^3*atan(sqrt(x + 1) - sqrt(x)), x, 9, -(sqrt(x)/8) + x^(3/2)/24 - x^(5/2)/40 + x^(7/2)/56 + (pi*x^4)/16 + atan(sqrt(x))/8 - (1/8)*x^4*atan(sqrt(x))]
    @test_int [x^2*atan(sqrt(x + 1) - sqrt(x)), x, 8, sqrt(x)/6 - x^(3/2)/18 + x^(5/2)/30 + (pi*x^3)/12 - atan(sqrt(x))/6 - (1/6)*x^3*atan(sqrt(x))]
    @test_int [x^1*atan(sqrt(x + 1) - sqrt(x)), x, 7, -(sqrt(x)/4) + x^(3/2)/12 + (pi*x^2)/8 + atan(sqrt(x))/4 - (1/4)*x^2*atan(sqrt(x))]
    @test_int [x^0*atan(sqrt(x + 1) - sqrt(x)), x, 6, sqrt(x)/2 + (pi*x)/4 - atan(sqrt(x))/2 - (1/2)*x*atan(sqrt(x))]
    @test_int [atan(sqrt(x + 1) - sqrt(x))/x^1, x, 6, (1/4)*pi*log(x) - (1/2)*im*PolyLog(2, (-im)*sqrt(x)) + (1/2)*im*PolyLog(2, im*sqrt(x))]
    @test_int [atan(sqrt(x + 1) - sqrt(x))/x^2, x, 6, -(pi/(4*x)) + 1/(2*sqrt(x)) + atan(sqrt(x))/2 + atan(sqrt(x))/(2*x)]
    @test_int [atan(sqrt(x + 1) - sqrt(x))/x^3, x, 7, -(pi/(8*x^2)) + 1/(12*x^(3/2)) - 1/(4*sqrt(x)) - atan(sqrt(x))/4 + atan(sqrt(x))/(4*x^2)]
    @test_int [atan(sqrt(x + 1) - sqrt(x))/x^4, x, 8, -(pi/(12*x^3)) + 1/(30*x^(5/2)) - 1/(18*x^(3/2)) + 1/(6*sqrt(x)) + atan(sqrt(x))/6 + atan(sqrt(x))/(6*x^3)]

    @test_int [atan(c*x/sqrt(a - c^2*x^2))^m/sqrt(d - c^2*d/a*x^2), x, 2, (sqrt(a - c^2*x^2)*atan((c*x)/sqrt(a - c^2*x^2))^(1 + m))/(c*(1 + m)*sqrt(d - (c^2*d*x^2)/a))]

    @test_int [atan(c*x/sqrt(a - c^2*x^2))^2/sqrt(d - c^2*d/a*x^2), x, 2, (sqrt(a - c^2*x^2)*atan((c*x)/sqrt(a - c^2*x^2))^3)/(3*c*sqrt(d - (c^2*d*x^2)/a))]
    @test_int [atan(c*x/sqrt(a - c^2*x^2))^1/sqrt(d - c^2*d/a*x^2), x, 2, (sqrt(a - c^2*x^2)*atan((c*x)/sqrt(a - c^2*x^2))^2)/(2*c*sqrt(d - (c^2*d*x^2)/a))]
    @test_int [1/atan(c*x/sqrt(a - c^2*x^2))^1/sqrt(d - c^2*d/a*x^2), x, 2, (sqrt(a - c^2*x^2)*log(atan((c*x)/sqrt(a - c^2*x^2))))/(c*sqrt(d - (c^2*d*x^2)/a))]
    @test_int [1/atan(c*x/sqrt(a - c^2*x^2))^2/sqrt(d - c^2*d/a*x^2), x, 2, -(sqrt(a - c^2*x^2)/(c*sqrt(d - (c^2*d*x^2)/a)*atan((c*x)/sqrt(a - c^2*x^2))))]
    @test_int [1/atan(c*x/sqrt(a - c^2*x^2))^3/sqrt(d - c^2*d/a*x^2), x, 2, -(sqrt(a - c^2*x^2)/(2*c*sqrt(d - (c^2*d*x^2)/a)*atan((c*x)/sqrt(a - c^2*x^2))^2))]

    @test_int [atan(e*x/sqrt(-a*e^2/b - e^2*x^2))^m/sqrt(a + b*x^2), x, 2, (sqrt(-((a*e^2)/b) - e^2*x^2)*atan((e*x)/sqrt(-((a*e^2)/b) - e^2*x^2))^(1 + m))/(e*(1 + m)*sqrt(a + b*x^2))]

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

    @test_int [atan(c*(a + b*x))*log(d*(a + b*x))/(a + b*x), x, 9, (im*log(d*(a + b*x))*PolyLog(2, (-im)*c*(a + b*x)))/(2*b) - (im*log(d*(a + b*x))*PolyLog(2, im*c*(a + b*x)))/(2*b) - (im*PolyLog(3, (-im)*c*(a + b*x)))/(2*b) + (im*PolyLog(3, im*c*(a + b*x)))/(2*b)]

    @test_int [ℯ^(c*(a + b*x))*atan(sinh(a*c + b*c*x)), x, 5, (ℯ^(a*c + b*c*x)*atan(sinh(c*(a + b*x))))/(b*c) - log(1 + ℯ^(2*c*(a + b*x)))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))*atan(cosh(a*c + b*c*x)), x, 8, (ℯ^(a*c + b*c*x)*atan(cosh(c*(a + b*x))))/(b*c) - ((1 - sqrt(2))*log(3 - 2*sqrt(2) + ℯ^(2*c*(a + b*x))))/(2*b*c) - ((1 + sqrt(2))*log(3 + 2*sqrt(2) + ℯ^(2*c*(a + b*x))))/(2*b*c)]
    @test_int [ℯ^(c*(a + b*x))*atan(tanh(a*c + b*c*x)), x, 13, atan(1 - sqrt(2)*ℯ^(a*c + b*c*x))/(sqrt(2)*b*c) - atan(1 + sqrt(2)*ℯ^(a*c + b*c*x))/(sqrt(2)*b*c) + (ℯ^(a*c + b*c*x)*atan(tanh(c*(a + b*x))))/(b*c) - log(1 + ℯ^(2*c*(a + b*x)) - sqrt(2)*ℯ^(a*c + b*c*x))/(2*sqrt(2)*b*c) + log(1 + ℯ^(2*c*(a + b*x)) + sqrt(2)*ℯ^(a*c + b*c*x))/(2*sqrt(2)*b*c)]
    @test_int [ℯ^(c*(a + b*x))*atan(coth(a*c + b*c*x)), x, 13, -(atan(1 - sqrt(2)*ℯ^(a*c + b*c*x))/(sqrt(2)*b*c)) + atan(1 + sqrt(2)*ℯ^(a*c + b*c*x))/(sqrt(2)*b*c) + (ℯ^(a*c + b*c*x)*atan(coth(c*(a + b*x))))/(b*c) + log(1 + ℯ^(2*c*(a + b*x)) - sqrt(2)*ℯ^(a*c + b*c*x))/(2*sqrt(2)*b*c) - log(1 + ℯ^(2*c*(a + b*x)) + sqrt(2)*ℯ^(a*c + b*c*x))/(2*sqrt(2)*b*c)]
    @test_int [ℯ^(c*(a + b*x))*atan(sech(a*c + b*c*x)), x, 8, (ℯ^(a*c + b*c*x)*atan(sech(c*(a + b*x))))/(b*c) + ((1 - sqrt(2))*log(3 - 2*sqrt(2) + ℯ^(2*c*(a + b*x))))/(2*b*c) + ((1 + sqrt(2))*log(3 + 2*sqrt(2) + ℯ^(2*c*(a + b*x))))/(2*b*c)]
    @test_int [ℯ^(c*(a + b*x))*atan(csch(a*c + b*c*x)), x, 5, (ℯ^(a*c + b*c*x)*atan(csch(c*(a + b*x))))/(b*c) + log(1 + ℯ^(2*c*(a + b*x)))/(b*c)]

    @test_int [((a + b*atan(c*x^n))*(d + e*log(f*x^m)))/x, x, 13, a*d*log(x) + (a*e*log(f*x^m)^2)/(2*m) + (im*b*d*PolyLog(2, (-im)*c*x^n))/(2*n) + (im*b*e*log(f*x^m)*PolyLog(2, (-im)*c*x^n))/(2*n) - (im*b*d*PolyLog(2, im*c*x^n))/(2*n) - (im*b*e*log(f*x^m)*PolyLog(2, im*c*x^n))/(2*n) - (im*b*e*m*PolyLog(3, (-im)*c*x^n))/(2*n^2) + (im*b*e*m*PolyLog(3, im*c*x^n))/(2*n^2)]
    end

max_stars_repo_path: test/5 Inverse trig functions/5.3 Inverse tangent/5.3.7 Inverse tangent functions.jl


text:
module Synapses

using ..Lif, ..LifConfig

export BaseSynapse, Synapse, PlasticSynapse, update!, send!, receive!, reset!, show, isconnected

import Base: show, getindex, setindex!

# =========================================================================== #

abstract type BaseSynapse end
abstract type DynamicSynapse <: BaseSynapse end

# =========================================================================== #
# macro pcheck(var)
#     quote
#         @assert(!isnan($var) && !isinf($var))
#     end
# end
# =========================================================================== #
# enumeration of concrete synpase types
# =========================================================================== #
@enum SynapseType Static DC Depressing Facilitating Tsodyks
# =========================================================================== #
# struct like container for extra synapse parameters
# =========================================================================== #
mutable struct SynapseParameters
    weight::Float
    delay::Time
    tau::Time
    aux::Dict{Symbol, Any}
end
SynapseParameters(weight, delay, tau) = SynapseParameters(weight, delay, tau, Dict{Symbol, Any}())
# --------------------------------------------------------------------------- #
getindex(s::SynapseParameters, k::Symbol) = return s.aux[k]
# --------------------------------------------------------------------------- #
setindex!(s::SynapseParameters, x::Any, k::Symbol) = (s.aux[k] = x)
# --------------------------------------------------------------------------- #
function show(io::IO, s::SynapseParameters)
    #TODO, FIXME: this formatting turns out quite poor...
    fn = [:weight, :delay, :tau]
    str = ""
    for f in fn
        str *= string(f) * ": " * string(getfield(s, f)) * "\n"
    end
    for f in keys(s.aux)
        str *= string(f) * ": " * string(s[f]) * "\n"
    end
    print(io, str[1:end-1])
end
# =========================================================================== #
# interal synapse factory function
# =========================================================================== #
function build_synapse!(s::T, weight, delay, tau; args...) where T<:BaseSynapse
    s.weight = weight
    s.queue = Array{Float,1}()
    s.onset = 0.0 #+Inf
    s.p = SynapseParameters(weight, delay, tau, Dict(args))

    return s
end
# =========================================================================== #
# external synapse factory function
# =========================================================================== #
function Synapse(typ::SynapseType, args...)
    local out::BaseSynapse
    if typ == Static
        out = StaticSynapse(args...)
    elseif typ == DC
        out = DCSynapse(args...)
    elseif typ == Depressing
        out = DepressingSynapse(args...)
    elseif typ == Facilitating
        out = FacilitatingSynapse(args...)
    elseif typ == Tsodyks
        out = TsodyksSynapse(args...)
    else
        error("Invalid synapse type")
    end

    return out
end
# ---------------------------------------------------------------------------- #
function Synapse()
    return StaticSynapse()
end
# =========================================================================== #
# static synapse type
# =========================================================================== #
mutable struct StaticSynapse <: BaseSynapse
    weight::Float         #current weight
    queue::Array{Float,1} #spike queue
    onset::Float          #current spike onset time
    p::SynapseParameters    #extra parameters

    function StaticSynapse(weight, delay, tau)
        return build_synapse!(new(), weight, delay, tau)
    end
end
StaticSynapse() = StaticSynapse(0.0, +Inf, 0.0)
# =========================================================================== #
# DC synapse type - gap junction
# =========================================================================== #
mutable struct DCSynapse <: BaseSynapse
    weight::Float         #current weight
    queue::Array{Float,1} #gap junction queue
    onset::Float          #current spike onset time
    p::SynapseParameters    #extra parameters

    function DCSynapse(weight, delay)
        return build_synapse!(new(), weight, delay, 0.0)
    end
end
DCSynapse() = DCSynapse(0.0, 0.0)
# =========================================================================== #
# Depressing synapse type
# =========================================================================== #
mutable struct DepressingSynapse <: DynamicSynapse
    weight::Float         #current weight
    queue::Array{Float,1} #spike queue
    onset::Float          #current spike onset time
    p::SynapseParameters    #extra parameters

    function DepressingSynapse(weight, delay, tau, dx, taux)
        return build_synapse!(new(), weight, delay, tau, x=1.0, dx=dx, taux=taux)
    end
end
DepressingSynapse() = DepressingSynapse(0.0, +Inf, 0.0, 0.0, 0.0)
# =========================================================================== #
# Facilitating synapse type
# =========================================================================== #
mutable struct FacilitatingSynapse <: DynamicSynapse
    weight::Float         #current weight
    queue::Array{Float,1} #spike queue
    onset::Float          #current spike onset time
    p::SynapseParameters    #extra parameters

    function FacilitatingSynapse(weight, delay, tau, dx, taux)
        return build_synapse!(new(), weight, delay, tau, x=0.0, dx=dx, taux=taux)
    end
end
FacilitatingSynapse() = FacilitatingSynapse(0.0, +Inf, 0.0, 0.0, 0.0)
# =========================================================================== #
# Tsodyks synapse type
# =========================================================================== #
mutable struct TsodyksSynapse <: DynamicSynapse
    weight::Float         #current weight
    queue::Array{Float,1} #spike queue
    onset::Float          #current spike onset time
    p::SynapseParameters    #extra parameters

    function TsodyksSynapse(weight, delay, tau, U, taud, tauf)
        return build_synapse!(new(), weight, delay, tau, U=U, taud=taud,
            tauf=tauf, x=1.0, y=0.0, u=U)
    end
end
TsodyksSynapse() = TsodyksSynapse(0.0, +Inf, 0.0, 0.0, 0.0, 0.0)
# =========================================================================== #
# BaseSynapse methods
# =========================================================================== #
function send!(pre::N, s::T, spk::Bool, tnow::Time) where {T<:BaseSynapse, N<:BaseNeuron}
    if spk
        #send a spike to the synapse
        push!(s.queue, tnow + s.p.delay)
    end
end
# --------------------------------------------------------------------------- #
function receive!(post::N, s::T, dt::Time, tnow::Time) where {T<:BaseSynapse, N<:BaseNeuron}
    #receive 'input' from the current state of the synapse
    if !isempty(s.queue) && (tnow >= s.queue[1])
        s.onset = popfirst!(s.queue)
        spike = true
    else
        spike = false
    end

    # NOTE: this update *MUST* be performed before calculation of the current
    # <out> below
    update!(s, dt, spike, tnow)

    # alpha synapse
    # t = (tnow - s.onset) / s.p.tau
    # out = s.weight * t * exp(1.0-t)

    #
    # out = s.weight * (1.0 + 1.2*exp(-t/15.0) - exp(-t/20.0))
    # biexp(t, a, b) = (a*b/(a-b)) * (exp(-t.*b) - exp(-t.*a))
    # biexp2(t, τa, τb) = (1.0/(τb-τa)) * (exp(-t./τb) - exp(-t./τa))
    #   where τa = 1/a and τb = 1/b

    # @pcheck(s.weight); @pcheck(tnow); @pcheck(s.onset); @pcheck(s.p.tau)

    # exponential synapse
    if s.onset > 0.0
        out = s.weight * exp(-(tnow-s.onset)/s.p.tau)
        # out = s.weight
    else
        out = 0.0
    end

    # update!(s, dt, spike, tnow)

    return out
end
# --------------------------------------------------------------------------- #
function reset!(s::T) where {T<:BaseSynapse}
    s.onset = 0.0
    s.queue = Array{Float,1}()
    s.weight = s.p.weight
end
# --------------------------------------------------------------------------- #
function update!(s::T, dt::Time, spike::Bool, tnow::Time) where {T<:BaseSynapse}
    #nothing to do, synaptic weight is static
end
# --------------------------------------------------------------------------- #
function isconnected(s::T) where {T<:BaseSynapse}
    if s.p.weight == 0.0 || isinf(s.p.delay)
        b = false
    else
        b = true
    end
    return b
end
# --------------------------------------------------------------------------- #
function show(io::IO, s::T) where {T<:BaseSynapse}
    #TODO, FIXME: this formatting is quite poor...
    show(io, s.p.weight)
end
# =========================================================================== #
# DynamicSynapse methods
# =========================================================================== #
function update!(s::DepressingSynapse, dt::Time, spike::Bool, tnow::Time)
    if spike
        s.weight = s.p.weight * s.p[:x]
        s.p[:x] *= 1.0 - s.p[:dx]
    end
    s.p[:x] += dt * ((1.0 - s.p[:x]) / s.p[:taux])
    # s.p[:x] *= (1.0 - (s.p[:dx] * spike))
end
# =========================================================================== #
function update!(s::FacilitatingSynapse, dt::Time, spike::Bool, tnow::Time)
    # NOTE: for a facilitating synapse, the weight should represent the maximum
    # output that the synapse can produce (i.e. an asymptotic limit on <x>), so
    # we then the release probability to begin growing <x> (which asymptotes at
    # 1)
    error("FacilitatingSynapse update() function is not complete")
    if spike
        s.weight = s.p.weight * (1.0 + s.p[:x])
        s.p[:x] += s.p[:dx]
    end

    s.p[:x] += dt * ((-s.p[:x]) / s.p[:taux])
end
# =========================================================================== #
# function update!(s::FDSynapse, dt::Time, spike::Bool, tnow::Time)
#     if spike
#         s.weight = s.p.weight * (1.0 + s.p[:x])
#     end
#     # (1.0/(τb-τa)) * (exp(-t./τb) - exp(-t./τa))
#     s.p[:x] =
# end
# =========================================================================== #
function update2!(s::TsodyksSynapse, dt::Time, spike::Bool, tnow::Time)
    # this is just the forward Euler approximation of update!
    if spike
        s.weight = s.p.weight * s.p[:u] * s.p[:x]

        # NOTE: it appears that this update can occur either before or after the
        # integration, so long as it is triggered by a spike
        # NOTE: <x> depends on <u> and *MUST* be updated first
        s.p[:x] *= 1.0 - s.p[:u]
        s.p[:u] += s.p[:U] * (1.0 - s.p[:u])
    end

    s.p[:u] += ((s.p[:U] - s.p[:u]) / s.p[:tauf]) * dt
    s.p[:x] += ((1.0 - s.p[:x]) / s.p[:taud])  * dt
end
# =========================================================================== #
function update!(s::TsodyksSynapse, dt::Time, spike::Bool, tnow::Time)

    s.p[:u] = s.p[:U] + (s.p[:u] - s.p[:U]) * exp(-dt / s.p[:tauf])
    s.p[:x] = 1.0 + (s.p[:x] - 1.0) * exp(-dt / s.p[:taud])

    if spike
        s.weight = s.p.weight * s.p[:u] * s.p[:x]
        tmp = 1.0 - s.p[:u]
        s.p[:x] *= tmp
        s.p[:u] += s.p[:U] * tmp
    end
end
# =========================================================================== #
# function update!(s::TsodyksSynapse1, dt::Time, spike::Bool, tnow::Time)
# # https://github.com/nest/nest-simulator/blob/master/models/tsodyks_connection.h
# #   tauf: [ms] time constant for fascilitation
# #   taud: [ms] time constant for recovery
# #   U:    asymptotic value of probability of release
# #   x:    amount of resources in recovered state
# #   y:    amount of resources in active state
# #   u:    actual probability of release
# #
# #  U     - maximum probability of release [0,1]
# #  tau   - time constant of synaptic current in ms
# #  tauf  - time constant for facilitation in ms
# #  taud  - time constant for depression in ms
# #  x     - initial fraction of synaptic vesicles in the readily
# #          releasable pool [0,1]
# #  y     - initial fraction of synaptic vesicles in the synaptic
# #          cleft [0,1]
# # SEE ALSO:
# #https://github.com/nest/nest-simulator/blob/master/models/tsodyks2_connection.h
#
#     # NOTE: exact solution requires tnow.... which we don't actually have
#     h = -(tnow - s.onset)
#
#     pu = s.p[:tauf] == 0.0 ? 0.0 : exp(h / s.p[:tauf])
#     py = exp(h / s.p.tau)
#     pz = s.p[:taud] == 0.0 ? 0.0 : exp(h / s.p[:taud])
#
#     pxy = ((pz - 1.0) * s.p[:taud] - (py - 1.0) * s.p.tau) / (s.p.tau - s.p[:taud])
#
#     pxz = 1.0 - pz
#
#     z = 1.0 - s.p[:x] - s.p[:y]
#
#     @pcheck(pu); @pcheck(py); @pcheck(pz); @pcheck(pxy); @pcheck(pxz); @pcheck(z)
#
#     # NOTE: order is critical here
#     s.p[:u] *= pu
#     s.p[:x] += pxy * s.p[:y] + pxz * z
#     s.p[:y] *= py
#
#     s.p[:u] += s.p[:U] * (1.0 - s.p[:u])
#
#     dw = s.p[:u] * s.p[:x]
#
#     s.p[:x] -= dw
#     s.p[:y] += dw
#
#     s.weight = s.p.weight * dw
#
# end
# =========================================================================== #
# DCSynapse methods
# =========================================================================== #
function send!(pre::N, s::DCSynapse, spk::Bool, tnow::Time) where {N<:BaseNeuron}
    push!(s.queue, pre.vm)
    s.onset = tnow #NOTE: this is not needed...
end
# --------------------------------------------------------------------------- #
function receive!(post::N, s::DCSynapse, dt::Time, tnow::Time) where {N<:BaseNeuron}

    if length(s.queue) < floor(s.p.delay/dt) || tnow < dt #floor(1.0/dt)
        out = 0.0
    else
        out = (shift!(s.queue) - post.vm)*s.weight
    end

    return out
end
# =========================================================================== #
function test2(dt=0.05)
    tnow = 0.0
    dur = 1200.0

    pre = LIFNeuron()
    post = LIFNeuron()
    Lif.set!(pre, :xi, 0.0)
    Lif.set!(post, :xi, 0.0)
    Lif.set!(post, :tau, 10.0)

    #                    wi   delay tau   U   taud  tauf
    s = Synapse(Tsodyks, 0.05, 3.0, 3.0, 0.4, 100.0, 10.0)

    @show(s.p)

    times = tnow:dt:dur
    inp = zeros(length(times),)
    x = zeros(length(times),)
    w = zeros(length(times),)
    post_vm = zeros(length(times),)

    spikes = 33.0:33.0:1000.0
    nspike = length(spikes)
    inc = 2
    next_spike = spikes[1]

    for k = 1:length(times)

        if times[k] >= next_spike
            send!(pre, s, true, times[k])
            if inc <= nspike
                next_spike = spikes[inc]
                inc += 1
            else
                next_spike = NaN
            end
        else
            send!(pre, s, false, times[k])
        end

        inp[k] = receive!(post, s, dt, times[k])

        step!(pre, 0.0, dt, times[k])
        step!(post, inp[k], dt, times[k])
        post_vm[k] = post.vm
        x[k] = s.p[:x]
        w[k] = s.weight
    end

    return times, inp, x, w, post_vm
end
# =========================================================================== #
function test(dt=0.043)
    tnow = 0.0
    dur = 150.0

    #2 noise-less neurons
    pre = LIFNeuron()
    Lif.set!(pre, :xi, 0.0)
    # pre.p.tau = 1.0
    post = LIFNeuron()
    Lif.set!(post, :xi, 0.0)

    # s = DCSynapse(1.0, 0.1)
    # s = Synapse(.38, 3.0, 2.0, 0.4, 60)

    #                        wi  delay tau  dx   taux
    # s = Synapse(Facilitating, 0.05, 3.0, 2.0, 0.4, 40.0)

    # NOTE: taud (and tauf) *MUST* be > dt
    #   to approximate a static synapse, set U = 1, taud = 0.05, (tauf doesn't
    #   matter)
    #                    wi   delay tau   U   taud  tauf
    s = Synapse(Tsodyks, 0.05, 3.0, 2.0, 0.4, 100.0, 20.0)

    @show(s.p)

    times = tnow:dt:dur
    inp = zeros(length(times),)
    x = zeros(length(times),)
    w = zeros(length(times),)
    # u = zeros(length(times),)
    post_vm = zeros(length(times),)

    isi = 5.0
    nspike = 5
    spikes = [25.0:isi:(25.0+isi*(nspike-1)); 100.0]
    nspike += 1
    inc = 2
    next_spike = spikes[1]

    for k = 1:length(times)

        if times[k] >= next_spike
            send!(pre, s, true, times[k])
            if inc <= nspike
                next_spike = spikes[inc]
                inc += 1
            else
                next_spike = NaN
            end
        else
            send!(pre, s, false, times[k])
        end

        inp[k] = receive!(post, s, dt, times[k])

        step!(pre, 0.0, dt, times[k])
        step!(post, inp[k], dt, times[k])
        post_vm[k] = post.vm
        x[k] = s.p[:x]
        w[k] = s.weight
        # u[k] = s.p[:u]
    end

    return times, inp, x, w, post_vm
end
# =========================================================================== #

end #END MODULE

# delegation:
# https://github.com/JuliaLang/DataStructures.jl/blob/master/src/delegate.jl

# NOTE: synaptic plasticity
#   http://www.scholarpedia.org/article/Short-term_plasticity
#   http://www.briansimulator.org/docs/stp.html
#   Markram et al (1998). Differential signaling via the same axon of
#   neocortical pyramidal neurons, PNAS 95(9):5323-8
#
#   weight *= u*x # "Synaptic weights are modulated by the product u*x (in 0..1)
#                 #  which is taken before updating the variables"
#   dx/dt = (1-x)/taux
#   du/dt = (U-u)/tauu
#   x *= (1-u)    #NOTE: order of update matters, x is first
#   u += U*(1-u)
#
# where:
#   x = fraction of resources available following release
#   u = fraciton of available resources ready for useful
#   U = increment of u produced by a spike
#
# gamma distribution
# (x.^(k-1).*exp(-x/theta))./(theta^k*gamma(k))

################################################################################
# dt = 0.05
# tau1 = 10.0
# tau2 = 20.0
# t = 0:dt:1999*dt
# d1 = 1.2
# d2 = 1.0
# x = zeros(length(t))
# for k in eachindex(x)
#    d1 += dt * ((-d1 / tau1))
#    d2 += dt * ((-d2) / tau2)
#    x[k] = d1 - d2
# end

# dt = 0.05
# tau = 10.0
# t = 0:dt:1999*dt
# d = 1.2
# x = zeros(length(t))
# for k in eachindex(x)
#     d += dt * ((-d + 1.0) / tau)
#     x[k] = d
# end

ext: jl


text:
using Pkg
Pkg.activate(joinpath(@__DIR__,".."))
using LinearAlgebra
using BenchmarkTools, Test
using Cthulhu
using Plots ; theme(:dark)
using SmoothedSpectralAbscissa ; const SSA = SmoothedSpectralAbscissa

function rand_gauss(n::Integer)
  return randn(n,n) ./ sqrt(n)
end
function rand_nonnormal(n::Integer,(ud::Real)=1.01)
  mat = rand_gauss(n)
  sh = schur(mat)
  upd = diagm(0=>fill(1.0,n),1=>fill(ud,n-1))
  return sh.vectors*upd*sh.Schur*inv(upd)*sh.vectors'
end

##

sizes = [25,50,75,100,125,150]

ssa_times = similar(sizes,Float64)
ssa_allocs = similar(ssa_times)

for (k,n) in enumerate(sizes)
  alloc = SSA.SSAAlloc(n)
  A=rand_nonnormal(n,1.0)
  be = @benchmark SSA.ssa!($A,nothing,$alloc)
  ssa_times[k] = median(be.times)
  ssa_allocs[k] = be.allocs
end

##
bar(sizes,ssa_times;leg=false,yscale=:log10)

@show ssa_times ./ ssa_times_old


text:
# Uses the same optimiser as the base instance. Only works if that solver understands SOCP.
struct ESCB2Exact <: ESCB2OptimisationAlgorithm; end

function optimise_linear_sqrtlinear(instance::CombinatorialInstance{T}, ::ESCB2Exact,
                                    linear::Dict{T, Float64}, sqrtlinear::Dict{T, Float64},
                                    sqrtlinear_weight::Float64, ::Int;
                                    with_trace::Bool=false) where T
  # Get a LP formulation
  if ! has_lp_formulation(instance)
    error("The exact formulation for ESCB-2 relies on a LP formulation, which the solver associated to this instance cannot provide.")
  end

  m, obj, vars = get_lp_formulation(instance, linear)

  # Encode the square root as a geometric mean, i.e. a SOCP.
  #     t = sqrt(b^T x) = geomean(b^T x, 1)
  #     (2 t)² + (1 - b^T x)² <= (1 + b^T x)²
  #     (1 + b^T x, 1 - b^T x, 2 * t) in SOCP, following JuMP's conventions
  # Add this reformulation only once in the model to avoid JuMP complaining.
  if :CombinatorialESCB2 in keys(m.ext)
    dict = m.ext[:CombinatorialESCB2]
    confidence_bonus = dict[:confidence_bonus]
    confidence_bonus_linear = dict[:confidence_bonus_linear]
    eq = dict[:eq]
    cone = dict[:cone]

    for i in keys(sqrtlinear) # TODO: what about coefficients that become zero? May that happen? In that case, they won't necessarily be in sqrtlinear, and therefore not updated.
      # Due to normalisation, all variables are considered on the left-hand side,
      # hence the minus sign.
      set_normalized_coefficient(eq, vars[i], - sqrtlinear[i])
    end
  else
    @variable(m, confidence_bonus_linear >= 0) # b^T x
    @variable(m, confidence_bonus >= 0) # t

    @constraint(m, eq, confidence_bonus_linear == sqrtlinear_weight * sum(sqrtlinear[i] * vars[i] for i in keys(sqrtlinear)))
    @constraint(m, cone, [1 + confidence_bonus_linear, 1 - confidence_bonus_linear, 2 * confidence_bonus] in SecondOrderCone())

    m.ext[:CombinatorialESCB2] = Dict(:confidence_bonus_linear => confidence_bonus_linear, :confidence_bonus => confidence_bonus, :eq => eq, :cone => cone)
  end

  @objective(m, Max, obj + confidence_bonus)

  # Solve this program.
  set_silent(m)
  t0 = time_ns()
  optimize!(m)
  t1 = time_ns()

  if termination_status(m) != MOI.OPTIMAL
    error("The model was not solved correctly.")
  end

  # Retrieve the optimum value and, if needed, the trace of the solving time.
  arm_indices = all_arm_indices(instance)
  sol = T[i for i in arm_indices if value(vars[i]) > 0.5]

  if with_trace
    run_details = ESCB2Details()
    run_details.n_iterations = 1
    run_details.best_objective = value(obj + confidence_bonus)
    run_details.solver_time = (t1 - t0) / 1_000_000_000

    return sol, run_details
  else
    return sol
  end
end


text:
struct AdjointLayeredArray{layer, T, A <: AbstractLayeredArray{layer, T}} <: AbstractLayeredMatrix{layer, T}
    parent::A
    function AdjointLayeredArray{layer, T, A}(parent::A) where {layer, T, A}
        layer > 0 || throw(ArgumentError("layer must be `layer > 0`, got $layer"))
        new{layer, T, A}(parent)
    end
end

Base.parent(x::AdjointLayeredArray) = x.parent

Base.adjoint(x::A) where {layer, T, A <: AbstractLayeredArray{layer, T}} = AdjointLayeredArray{layer, T, A}(x)
Base.adjoint(x::AdjointLayeredArray) = parent(x)

Base.size(x::AdjointLayeredArray) = (size(parent(x), 2), size(parent(x), 1))
Base.axes(x::AdjointLayeredArray) = (axes(parent(x), 2), axes(parent(x), 1))

Base.IndexStyle(::Type{<: AdjointLayeredArray{<: Any, <: Any, A}}) where {A} = IndexStyle(A)

# IndexLinear
@inline function Base.getindex(x::AdjointLayeredArray, i::Integer)
    @boundscheck checkbounds(x, i)
    @inbounds parent(x)[i]
end
@inline function Base.setindex!(x::AdjointLayeredArray, v, i::Integer)
    @boundscheck checkbounds(x, i)
    @inbounds parent(x)[i] = v
    x
end

# IndexCartesian
@inline function Base.getindex(x::AdjointLayeredArray, i::Integer, j::Integer)
    @boundscheck checkbounds(x, i, j)
    @inbounds parent(x)[j, i]
end
@inline function Base.setindex!(x::AdjointLayeredArray, v, i::Integer, j::Integer)
    @boundscheck checkbounds(x, i, j)
    @inbounds parent(x)[j, i] = v
    x
end


text:
struct TCG <: TRSPSolver
end
summary(::TCG) = "Steihaug-Toint Truncated CG"

function (ms::TCG)(∇f, H, Δ::T, s, scheme, λ0=0; abstol=1e-10, maxiter=min(5,length(s)), κeasy=T(1)/10, κhard=T(2)/10) where T
	# We only know that TCG has its properties if we start with s = 0
	s .= T(0)
	# g is the gradient of the quadratic model with the step α*p
	g = copy(∇f)
	v = apply_inverse_preconditioner(InPlace(), P, copy(g), g)
	p = -v
	for i = 1:maxiter
		# Check for negative curvature
		κ = dot(p, H*p)
		if κ ⩽ T(0)
			# positive_root_of
			σ = positive_root_of()
			s .= s .+ σ.*p
		end
		α = dot(g, v)/κ
		step_length = M_norm(s .+ α.*p)
		# If next step is outside, step to the boundary
		if step_length ⫺ Δ
			σ = positive_root_of()
			s .= s .+ σ.*p
		end
		s .= s + α*p
		den = dot(g, v)
		g .= g .+ α.*(H*p)
		v = apply_inverse_preconditioner(InPlace(), P, v, g)
		β = dot(g, v)/den
		p .= -v .+ β.*p
	end
end


text:
# This file illustrates the simulation of MFSK signalling. The simulation is performed using
# waveform AWGN channel.

using Plots
using FFTW
using LinearAlgebra
using Statistics
using DigitalCommunications

# Settings
k = 2                                       # Bits per symbol
M = 2^3                                     # Constellation size
T = 1.                                      # Symbol duration
f = 1 / T                                   # Symbols per second
fs = 10                                     # Sampling frequency
ts = 1 / fs                                 # Sampling period
t = 0 : ts : T - ts                         # Time vector for symbol duration
l = length(t)                               # Number of samples per symbol
g(t) = 0 ≤ t ≤ T ? sqrt(2 * Eg / T) : 0     # Modulation pulse: rectangular pulse
Eg = 1                                      # Modulation pulse energy: rectangular pulse
nsymbols = Int(1e6)                         # Number of symbols
Δf = 1 / T                                  # 2 times the minimum frequency seperation for orthogonality

# Define basis
basis = map(m -> (t -> g(t) * exp(1im * 2π * m * Δf * t)), 1 : M)
basisvals = map(base -> base.(t), basis)
alphabet = map(m -> setindex!(zeros(M), 1, m), 1 : M)
references = map(1 : M) do mi
    Ai = alphabet[mi] * sqrt(2 * Eg)
    sum(Ai .* basisvals)
end

# Simulate the system
ebno = collect(0 : 1 : 10)
esno = ebno .+ 10 * log10(k)
m = rand(1 : M, nsymbols)
tx = map(m) do mi
    Ai = alphabet[mi] * sqrt(Eg)
    sum(Ai .* basisvals)
end
Es = mean(txi -> sum(abs.(txi).^2) * ts, tx)
symerr = zeros(length(esno))
for i = 1 : length(symerr)
    # Generate noise
    σ = sqrt(Es * fs / 2 / (10^(esno[i] / 10)))
    n = [σ * (randn(l) + 1im * randn(l)) for i = 1 : length(tx)]
    # Corrupt transmitted signal
    rx = tx + n
    # Detect symbols
    mr = map(rx) do rxi
        argmax(map(reference -> real(rxi ⋅ reference), references))
    end
    # Calculate symbol
    ber = sum(m .!= mr) / length(m)
    symerr[i] = ber
end

# Plots
plt = plot(title="$M-FSK", xlabel="esno [dB]", ylabel="Pe")
plot!(esno, berfsk.(esno, M), marker=:circle, yscale=:log10, label="theoretical")
plot!(esno, symerr, marker=:circle, yscale=:log10, label="montecarlo")


text:
"""
    r = RationalTransferFunction(num::AbstractPolynomial, den::AbstractPolynomial, Ts:Real)

Construct a rational transfer function model `r` from its numerator and denominator polynomials
`num` and `den`, respectively, and sampling time `Ts`.

If `r::RationalTransferFunction{T,λ,P <: Polynomial(T,λ),Ts}` is a rational transfer function system model
object defined as `r(λ) = num(λ)/den(λ)`, where  `num(λ)` and `den(λ)` are polynomials
with coefficients in `T` and with the indeterminate `λ`, and `Ts` is the sampling time, then:

`r.num` is the numerator polynomial `num(λ)`;

`r.den` is the denominator polynomial `den(λ)`.

The sampling time `Ts` can have the following values:

  - `Ts = 0` for a continuous-time system and
             `λ = s` is the complex variable in the Laplace transform;

  - `Ts > 0` or `Ts = -1` for a discrete-time system and
             `λ = z` is the complex variable in the `Z`-transform;
             `Ts = -1` indicates a discrete-time system with an unspecified sampling time.

The sampling time can be obtained as `r.Ts`.
The symbol (or _variable_) used for the indeterminate `λ` is the common symbol used for the
indeterminates of the polynomials `num(λ)` and `den(λ)` and can be obtained as `r.var`.
The roots of the numerator polynomial `num(λ)` (also called _zeros_ of `r(λ)`)
can be obtained as `r.zeros`, while the roots of the denominator polynomial `den(λ)`
(also called _poles_ of `r(λ)`) can be obtained as `r.poles`.
The ratio of the leading polynomial coefficients of `num(λ)` and `den(λ)`
(also called _gain_ of `r(λ)`) can be obtained as `r.gain`.
"""
struct RationalTransferFunction{T,X,P<:AbstractPolynomial{T,X},Ts} <: AbstractRationalFunction{T,X,P}
    ## Define RationalTransferFunction as a subtype of AbstractRationalFunction
    num::P
    den::P
    function RationalTransferFunction{T,X,P,Ts}(num::P, den::P) where{T, X, P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        new{T,X,P,Ts}(num,den)
    end
    function RationalTransferFunction{T,X,P,Ts}(num::P, den::P, ts::Union{Real,Nothing}) where {T, X, P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        check_Ts(Ts,ts)
        new{T,X,P,Ts}(num,den)
    end
    # can promote constants to polynomials if S <: T
    function  RationalTransferFunction{T,X,P,Ts}(num::S, den::P, ts::Union{Real,Nothing}) where {S, T, X, P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        check_Ts(Ts,ts)
        new{T,X,P,Ts}(convert(P,num),den)
    end
    function  RationalTransferFunction{T,X,P,Ts}(num::P, den::S, ts::Union{Real,Nothing}) where {S, T, X, P<:AbstractPolynomial{T,X}, Ts}
        check_den(den)
        check_Ts(Ts,ts)
        new{T,X,P,Ts}(num, convert(P,den))
    end
    function RationalTransferFunction{T,X,P}(num::P, den::P, Ts::Union{Real,Nothing}) where {T,X,P<:AbstractPolynomial{T,X}}
        check_den(den)
        Ts′ = standardize_Ts(Ts)
        #new{T,X,P,Val(Ts′)}(num,den)
        new{T,X,P,Ts′}(num,den)
    end
end

#Base.eltype(sys::RationalTransferFunction) = eltype(sys.num)
_eltype(R::RationalTransferFunction) = eltype(R.num)
_eltype(R::VecOrMat{<:RationalTransferFunction}) = length(R) == 0 ? eltype(eltype(Polynomials.coeffs.(numpoly.(R)))) : eltype(eltype(R[1]))
_eltype(R::VecOrMat{<:Polynomial}) = length(R) == 0 ? eltype(eltype(Polynomials.coeffs.(R))) : eltype(eltype(R[1]))

#(pq::RationalTransferFunction)(x) = Polynomials.eval_rationalfunction(x, pq)
(pq::RationalTransferFunction)(x) = pq.num(x)/pq.den(x)

function Base.getproperty(F::RationalTransferFunction{T,X,P,Ts}, d::Symbol)  where {T,X,P,Ts}
    if d === :Ts
        return Ts
    elseif d === :var
        return X
    elseif d === :zeros
        return roots(getfield(F, :num))
    elseif d === :poles
        return roots(getfield(F, :den))
    elseif d === :gain
        return last(getfield(F, :num)) / last(getfield(F, :den))
    else
        getfield(F, d)
    end
end
Base.propertynames(F::RationalTransferFunction) =
    (:Ts, :var, :zeros, :poles, :gain, fieldnames(typeof(F))...)
poles(F::RationalTransferFunction) = F.poles
poles(F::Polynomial{T}) where T = zeros(T,0)
gain(F::RationalTransferFunction) = F.gain
gain(F::Polynomial) = last(F.coeffs)
gpole(F::RationalTransferFunction) = [F.poles; Inf*ones(Int,max(0,degree(F.num)-degree(F.den)))]
gzero(F::RationalTransferFunction) = [F.zeros; Inf*ones(Int,max(0,degree(F.den)-degree(F.num)))]

function Base.convert(::Type{PQ}, pq::RationalTransferFunction) where {PQ <:RationalFunction}
    p,q = pq
    p//q
end

# alternate constructors
function RationalTransferFunction(p′::P, q′::Q, Ts::Union{Real,Nothing}) where {T, X, P<:AbstractPolynomial{T,X},
                                                                                S,  Q<:AbstractPolynomial{S,X}}
    p, q = promote(p′, q′)
    RationalTransferFunction{eltype(p),X,typeof(p)}(p,q,Ts)
end
function RationalTransferFunction(p′::P, q′::S, Ts::Union{Real,Nothing}) where {S, T, X, P<:AbstractPolynomial{T,X}}
    p, q = promote(p′, q′)
    RationalTransferFunction{eltype(p),X,typeof(p)}(p,q,Ts)
 end
 function RationalTransferFunction(p′::S, q′::P, Ts::Union{Real,Nothing}) where {S, T, X, P<:AbstractPolynomial{T,X}}
    p, q = promote(p′, q′)
    RationalTransferFunction{eltype(q),X,typeof(q)}(p,q,Ts)
 end

function Polynomials.rational_function(::Type{PQ}, p::P, q::Q) where {PQ <:RationalTransferFunction,
                                                          T, X, P<:AbstractPolynomial{T,X},
                                                          S,    Q<:AbstractPolynomial{S,X}}
    RationalTransferFunction(promote(p,q)..., sampling_time(PQ))
end

## helpers for constructors
# standardize Ts or throw error
function standardize_Ts(Ts)
    isnothing(Ts) || Ts >= 0 || Ts == -1 ||
        throw(ArgumentError("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)"))
    Ts′ = isnothing(Ts) ? Ts : Float64(Ts)
end
function check_Ts(Ts, ts)
    # ValT(Ts) == promote_Ts(ValT(Ts), ts) || throw(ArgumentError("sampling times have mismatch"))
    ValT(Val(Ts)) == promote_Ts(ValT(Val(Ts)), ts) || throw(ArgumentError("sampling times have mismatch"))
end
function check_den(den)
    iszero(den) && throw(ArgumentError("Cannot create a rational function with zero denominator"))
end

ValT(::Val{T}) where {T} = T
# sampling_time(pq::RationalTransferFunction{T,X,P,Ts}) where {T,X,P,Ts} = ValT(Ts)
# sampling_time(::Type{𝑷}) where {T,X,P,Ts, 𝑷<:RationalTransferFunction{T,X,P,Ts}} = ValT(Ts)
sampling_time(pq::RationalTransferFunction{T,X,P,Ts}) where {T,X,P,Ts} = ValT(Val(Ts))
sampling_time(::Type{𝑷}) where {T,X,P,Ts, 𝑷<:RationalTransferFunction{T,X,P,Ts}} = ValT(Val(Ts))
sampling_time(::Type{<:Number}) = nothing
sampling_time(p::Number) = nothing

"""
    r = rtf(num, den; Ts = rts, var = rvar )

Create the rational transfer function `r(λ) = num(λ)/den(λ)` with the polynomials `num(λ)` and `den(λ)`,
sampling time `rts` and variable name `rvar`, representing
the transfer function of a single-input single-output system of the form

    Y(λ) = r(λ) U(λ),

where `U(λ)` and `Y(λ)` are the Laplace or `Z` transformed input `u(t)` and output `y(t)`, respectively,
and `λ = s`, the complex variable in the Laplace transform, if `rts = 0`, or  `λ = z`,
the complex variable in the `Z` transform, if `rts ≠ 0`.
Both `num` and `den` can be real or complex numbers as well.

The resulting `r` is such that `r.Ts = rts` (default `rts = 0`) and `r.var = rvar`.
The default value of `rvar` is `rvar = Polynomials.indeterminate(num)` if `num` is a polynomial,
`rvar = Polynomials.indeterminate(den)` if `num` is a number and `den` is a polynomial, and
`rvar = :s` if both `num` and `den` are numbers.
"""
function rtf(num::Polynomial, den::Polynomial; Ts::Real = 0,
             var::Symbol = promote_var(num,den))
    T = promote_type(eltype(num),eltype(den))
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(num.coeffs), Polynomial{T,var}(den.coeffs))
end
function rtf(num::Polynomial, den::Number; Ts::Real = 0,  var::Symbol = Polynomials.indeterminate(num))
    T = promote_type(eltype(num),eltype(den))
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(num.coeffs), Polynomial{T,var}(den))
end
function rtf(num::Number, den::Polynomial; Ts::Real = 0,  var::Symbol = Polynomials.indeterminate(den))
    T = promote_type(eltype(den),eltype(num))
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(num),Polynomial{T,var}(den.coeffs))
end
function rtf(num::Number, den::Number; Ts::Real = 0, var::Symbol = :s)
    T = promote_type(eltype(num),eltype(den))
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(num),Polynomial{T,var}(den))
end

# display sys
Base.print(io::IO, sys::RationalTransferFunction) = show(io, sys)
Base.show(io::IO, sys::RationalTransferFunction) = show(io, MIME("text/plain"), sys)

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sys::RationalTransferFunction)
    summary(io, sys); println(io)
    p, q = Polynomials.pqs(sys)
    if order(sys) == 0
        printpoly(io, p)
        println(io, "\n\nStatic gain.")
        return
    end
    if q == one(q)
        printpoly(io, p, descending_powers=true)
    else
        count(p.coeffs .!= 0) <= 1 ? printpoly(io, p, descending_powers=true) :
             (print(io,"("); printpoly(io, p, descending_powers=true); print(io, ")") )
        print(io, " // ")
        count(q.coeffs .!= 0) <= 1 ? printpoly(io, q, descending_powers=true) :
             (print(io,"("); printpoly(io, q, descending_powers=true); print(io, ")") )
    end
    sys.Ts < 0 && println(io, "\n\nSample time: unspecified.")
    sys.Ts > 0 && println(io, "\n\nSample time: $(sys.Ts) second(s).")
    iszero(sys.Ts) ? println(io, "\n\nContinuous-time rational transfer function model.") :
                     println(io, "Discrete-time rational transfer function model.")
end

function show_rp(io::IO, sys::RationalTransferFunction)
    p,q = Polynomials.pqs(sys)
    if order(sys) == 0
        printpoly(io, p, descending_powers=true)
        return
    end
    if q == one(q)
        printpoly(io, p, descending_powers=true)
    else
        count(p.coeffs .!= 0) <= 1 ? printpoly(io, p, descending_powers=true) :
             (print(io,"("); printpoly(io, p, descending_powers=true); print(io, ")") )
        print(io, " // ")
        count(q.coeffs .!= 0) <= 1 ? printpoly(io, q, descending_powers=true) :
             (print(io,"("); printpoly(io, q, descending_powers=true); print(io, ")") )
    end
end

Base.print(io::IO, R::VecOrMat{<:RationalTransferFunction}) = show(io, R)
Base.show(io::IO, R::VecOrMat{<:RationalTransferFunction}) = show(io, MIME("text/plain"), R)

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, R::VecOrMat{<:RationalTransferFunction})
    summary(io, R); println(io)
    p, m = size(R,1), size(R,2)
    for j=1:m
         for i=1:p
            println(io, "From input $j to output $i")
            show_rp(io, R[i,j])
            print(io, "\n\n")
        end
    end
    Ts = sampling_time(R[1,1])
    Ts < 0 && println(io, "Sample time: unspecified.")
    Ts > 0 && println(io, "Sample time: $Ts second(s).")
    iszero(Ts) ? println(io, "Continuous-time rational transfer function matrix model.") :
                 println(io, "Discrete-time rational transfer function matrix model.")
end

function promote_var(num::Polynomial, den::Polynomial)
    m = length(num)
    n = length(den)
    if m > 1 &&  n > 1
       var = Polynomials.indeterminate(num)
       var != Polynomials.indeterminate(den) && error("The numerator and denominator polynomials must have the same variable")
    elseif m > 1
       var = Polynomials.indeterminate(num)
    elseif  n > 1
       var = Polynomials.indeterminate(den)
    else
       var = :s   # use the default variable for continuous-time systems
    end
    return var
end

"""
    r = rtf(f; Ts = rts, var = rvar)

Create the rational transfer function `r(λ) = f(λ)` with sampling time `rts` and variable name `rvar`
such that:

(1) if `f(λ)` is a rational transfer function, then `r.Ts = rts` (default `rts = f.Ts`) and `r.var = rvar` (default `rvar = f.var`);

(2) if `f(λ)` is a rational function, then `r.Ts = rts` (default `rts = 0`) and `r.var = rvar` (default `rvar = Polynomials.indeterminate(f.num)`);

(3) if `f(λ)` is a polynomial, then `r.Ts = rts` (default `rts = 0`) and `r.var = rvar` (default `rvar = Polynomials.indeterminate(f)`);

(4) if `f(λ)` is a ral or complex number, then `r.Ts = rts` (default `rts = 0`) and `r.var = rvar` (default `rvar = :s`);
"""
function rtf(f::RationalTransferFunction; Ts::Real = f.Ts, var::Symbol = f.var)
    T = eltype(eltype(f))
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(f.num.coeffs), Polynomial{T,var}(f.den.coeffs))
end
function rtf(f::AbstractRationalFunction; Ts::Real = 0, var::Symbol = Polynomials.indeterminate(f.num))
    T = eltype(eltype(f))
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(f.num.coeffs), Polynomial{T,var}(f.den.coeffs))
end
function rtf(p::P; Ts::Real = 0, var::Symbol = Polynomials.indeterminate(p)) where {T,X,P<:AbstractPolynomial{T,X}}
    #T = eltype(p)
    #return RationalTransferFunction{T}(Polynomial{T}(p.coeffs,var), Polynomial{T}(one(T),var), Ts)
    #return RationalTransferFunction{T,var}(Polynomial{T,var}(p.coeffs), Polynomial{T,var}(one(T)), Ts)
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(p.coeffs), Polynomial{T,var}(one(T)))
end
function rtf(num::Number; Ts::Real = 0, var::Symbol = :s)
    T = eltype(num)
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(num), Polynomial{T,var}(one(T)))
end

"""
    r = rtf(var; Ts = rts)
    r = rtf('s') or r = rtf('z'; Ts = rts)
    r = rtf("s") or r = rtf("z"; Ts = rts)
    r = rtf(:s) or r = rtf(:z; Ts = rts)

Create the rational transfer function `r(λ) = λ`, such that `r.var` and `r.Ts` are set as follows:

    (1) `r.var = :s` and `r.Ts = 0` if `var = 's'`, or `var = "s"` or `var = :s` ;

    (2) `r.var = :z` and `r.Ts = rts` if `var = 'z'`, or `var = "z"` or `var = :z`;

    (3) `r.var = var` and `r.Ts = rts` (default `rts = 0.`) otherwise.
"""
function rtf(::Type{T}, var::Union{AbstractString,Char,Symbol}; Ts::Union{Real,Missing} = missing) where T
    if var == "s" || var == 's' || var === :s
       ismissing(Ts) && (Ts = 0)
       Ts == 0. ? (return rtf(Polynomial{T}([0,1],:s),Ts=Ts)) : error("Ts must be zero for a continuous-time system")
    elseif var == "z" || var == 'z' || var === :z
       ismissing(Ts) && (Ts = -1)
       (Ts == -1 || Ts > 0) ? (return rtf(Polynomial{T}([0,1],:z),Ts = Ts)) : error("Ts must be either a positive number or -1 (unspecified)")
    else
       ismissing(Ts) && (Ts = 0.)
       Ts >= 0 || Ts == -1 || error("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)")
       return rtf(Polynomial{T}([0,1],var),Ts = Ts)
    end
end
rtf(var::Union{AbstractString,Char,Symbol}; Ts::Union{Real,Missing} = missing) = rtf(Int, var; Ts = Ts)
function Base.:/(p::AbstractPolynomial,q::AbstractPolynomial)
    RationalFunction(p,q)
end

"""
    r = rtf(z, p, k; Ts = rts, var = rvar)

Create from the roots (zeros) `z`, poles `p`, gain `k`, sampling time `rts` and variable name `rvar`
the rational transfer function `r(λ) = k*num(λ)/den(λ)`, where `num(λ)` and `den(λ)` are monic polynomials
with roots equal `z` and `p`, respectively, and such that `r.Ts = rts` (default `rts = 0`)
 and `r.var = rvar` (default `rvar = :s` if `Ts = 0` or `rvar = :z` if `Ts ≠ 0`).
"""
function rtf(zer::Vector{<:Number}, pol::Vector{<:Number}, k::Number; Ts::Real = 0,
             var::Symbol = (Ts == 0 ? :s : :z) )
    Ts >= 0. || Ts == -1. || error("Ts must be either a positive number, 0 (continuous system), or -1 (unspecified)")
    T = promote_type(eltype(zer),eltype(pol),eltype(k))
    reald = false   # check data correspond to real numerator and denominator
    if eltype(k) <: Real
       tempc = zer[imag.(zer) .> 0]
       if isempty(tempc)
          realz = true
       else
          tempc1 = conj(zer[imag.(zer) .< 0])
          realz = isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))])
       end
       if realz
          tempc = pol[imag.(pol) .> 0]
          if isempty(tempc)
             realp = true
          else
             tempc1 = conj(pol[imag.(pol) .< 0])
             realp = isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))])
          end
          reald = realp
       end
    end
    num = reald ? real(k*fromroots(zer)) :  k*fromroots(zer)
    den = reald ? real(fromroots(pol)) :  k*fromroots(pol)
    return RationalTransferFunction{T,var,Polynomial{T,var},Float64(Ts)}(Polynomial{T,var}(num.coeffs), Polynomial{T,var}(den.coeffs))
end

## ----
function Base.convert(::Type{PQ}, p::R) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts}, R<:AbstractRationalFunction}
    PQ(Polynomial(T.(p.num.coeffs),X), Polynomial(T.(p.den.coeffs),X), sampling_time(PQ))
end

function Base.convert(::Type{PQ}, p::P1) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts}, P1<:AbstractPolynomial}
    #PQ(p, one(p), sampling_time(PQ))
    PQ(Polynomial(T.(p.coeffs),X), Polynomial(one(T),X), sampling_time(PQ))
end
function Base.convert(::Type{PQ}, p::Number) where {PQ <: RationalTransferFunction}
    PQ(p, one(eltype(PQ)), sampling_time(PQ))
end

function promote_Ts(p,q)
    Ts1,Ts2 = sampling_time.((p,q))
    promote_Ts(Ts1, Ts2)
end

function promote_Ts(Ts1::Union{Real,Nothing}, Ts2::Union{Real,Nothing})
    isnothing(Ts1) && (return Ts2)
    isnothing(Ts2) && (return Ts1)
    Ts1 == Ts2 && (return Ts1)
    Ts1 == -1 && (Ts2 > 0 ? (return Ts2) : error("Sampling time mismatch"))
    Ts2 == -1 && (Ts1 > 0 ? (return Ts1) : error("Sampling time mismatch"))
    error("Sampling time mismatch")
end

function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts},
                                                           T′,X′,P′,Ts′,PQ′ <: RationalTransferFunction{T′,X′,P′,Ts′}}
    S = promote_type(T,T′)
    Polynomials.assert_same_variable(X,X′)
    Y = X
    Q = promote_type(P, P′)
    ts = promote_Ts(PQ, PQ′)
    #RationalTransferFunction{S,Y,Q,Val(ts)}
    RationalTransferFunction{S,Y,Q,ts}
end
function Base.promote_rule(::Type{PQ}, ::Type{PQ′}) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts},
                                                           T′,X′,P′,PQ′ <: RationalFunction{T′,X′,P′}}
    S = promote_type(T,T′)
    Polynomials.assert_same_variable(X,X′)
    Y = X
    Q = promote_type(P, P′)
    ts = sampling_time(PQ)
     #RationalTransferFunction{S,Y,Q,Val(ts)}
    RationalTransferFunction{S,Y,Q,ts}
end

function Base.promote_rule(::Type{PQ}, ::Type{P′}) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts},
                                                          T′,X′,P′ <: Polynomial{T′,X′}}
    S = promote_type(T,T′)
    Polynomials.assert_same_variable(X,X′)
    Y = X
    Q = promote_type(P, P′)
    ts = sampling_time(PQ)
    #RationalTransferFunction{S,Y,Q,Val(ts)}
    RationalTransferFunction{S,Y,Q,ts}
end

#Base.promote_rule(::Type{PQ}, ::Type{P}) where {PQ <: RationalTransferFunction, P<:AbstractPolynomial} = PQ
Base.promote_rule(::Type{PQ}, ::Type{P}) where {PQ <: RationalTransferFunction, P<:Number} = PQ

Base.zero(::Type{PQ}) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts}}  = RationalTransferFunction{T,X,P,Ts}(zero(P), one(P))
Base.zero(f::RationalTransferFunction) = zero(typeof(f))
Base.one(::Type{PQ}) where {T,X,P,Ts,PQ <: RationalTransferFunction{T,X,P,Ts}}  = RationalTransferFunction{T,X,P,Ts}(one(P), one(P))
Base.one(f::RationalTransferFunction) = one(typeof(f))

"""
     zpk(r) -> (z, p, k)

Compute the roots (zeros) `z`, poles `p` and gain `k` of the rational transfer function `r(λ)`.
"""
function zpk(F::RationalTransferFunction)
    return F.zeros, F.poles, F.gain
end
function zpk(F::Polynomial{T}) where T
    return roots(F), zeros(T,0), gain(F)
end

"""
     rt = adjoint(r)

Compute the adjoint `rt(λ)` of the rational transfer function `r(λ)` such that for
`r(λ) = num(λ)/den(λ)` we have:

    (1) `rt(λ) = conj(num(-λ))/conj(num(-λ))`, if `r.Ts = 0`;

    (2) `rt(λ) = conj(num(1/λ))/conj(num(1/λ))`, if `r.Ts = -1` or `r.Ts > 0`.
"""
function Base.adjoint(f::RationalTransferFunction)
    if f.Ts == 0
       p1 = copy(conj(f.num.coeffs))
       i1 = 2:2:length(p1)
       p1[i1] = -p1[i1]
       p2 = copy(conj(f.den.coeffs))
       i2 = 2:2:length(p2)
       p2[i2] = -p2[i2]
       return RationalTransferFunction(Polynomial(p1,Polynomials.indeterminate(f.num)), Polynomial(p2,Polynomials.indeterminate(f.den)), f.Ts)
    else
       m = degree(f.num)+1
       n = degree(f.den)+1
       p1 = reverse(conj(f.num.coeffs[1:m]))
       p2 = reverse(conj(f.den.coeffs[1:n]))
       if m >= n
          return RationalTransferFunction(Polynomial(p1,Polynomials.indeterminate(f.num)), Polynomial([zeros(eltype(p2),m-n); p2],Polynomials.indeterminate(f.den)), f.Ts)
       else
          return RationalTransferFunction(Polynomial([zeros(eltype(p1),n-m); p1],Polynomials.indeterminate(f.num)), Polynomial(p2,Polynomials.indeterminate(f.den)), f.Ts)
       end
    end
end
transpose(f::RationalTransferFunction) = f

denpoly(f::RationalTransferFunction) = f.den
numpoly(f::RationalTransferFunction) = f.num
denpoly(f::Number) = one(Polynomial{eltype(f)})
numpoly(f::Number) = Polynomial{eltype(f)}(f)

"""
    n = order(r)

Determine the order `n` of a rational transfer function `r` as the maximum of degrees of its
numerator and denominator polynomials (`n` is also known as the _McMillan degree_ of `r`).
"""
function order(r::RationalTransferFunction)
    max(degree(r.num),degree(r.den))
end

# operations
function Base.:+(p::R, q::R) where {T,X,P,Ts,R <: RationalTransferFunction{T,X,P,Ts}}
    p0,p1 = pqs(p)
    q0,q1 = pqs(q)
    Polynomials.rational_function(R, p0*q1 + p1*q0, p1*q1)
end
Base.:+(p::Number, q::RationalTransferFunction) = q + p
function Base.:+(p::R,  q::Number) where {T,X,Q,Ts,R <: RationalTransferFunction{T,X,Q,Ts}}
    Tx = promote_type(T,eltype(q))
    p0,p1 = pqs(p)
    #p + q*one(p)
    Polynomials.rational_function(RationalTransferFunction{Tx,X,Polynomial{Tx,X},Ts}, p0 + p1*q, p1)
end
Base.:+(p::Polynomial, q::RationalTransferFunction) = q + p
#Base.:+(p::RationalTransferFunction, q::RationalTransferFunction) = +(promote(p,q)...)
function Base.:+(p::R, q::P) where {T,X,Q,Ts,R <: RationalTransferFunction{T,X,Q,Ts}, T1, P <: Polynomial{T1,X}}
    Tx = promote_type(T,T1)
    p0,p1 = pqs(p)
    Polynomials.rational_function(RationalTransferFunction{Tx,X,Polynomial{Tx,X},Ts}, p0 + p1*q, p1)
end

function Base.:-(p::R, q::R) where {T,X,P,Ts,R <: RationalTransferFunction{T,X,P,Ts}}
    p0,p1 = pqs(p)
    q0,q1 = pqs(q)
    Polynomials.rational_function(R, p0*q1 - p1*q0, p1*q1)
end

function Base.:*(p::R, q::R) where {T,X,P,Ts,R <: RationalTransferFunction{T,X,P,Ts}}
    p0,p1 = pqs(p)
    q0,q1 = pqs(q)
    Polynomials.rational_function(R, p0*q0, p1*q1)
end

# equality
import Base: ==
function ==(p::RationalTransferFunction{T,X,P,T1}, q::RationalTransferFunction{S,Y,Q,T2}) where {T,X,P,T1,S,Y,Q,T2}
    isconstant(p.num) && isconstant(p.den) && isconstant(q.num) && isconstant(q.den) && p(0) == q(0) && return true
    (X == Y && T1 == T2) || return false
    p₀, p₁ = pqs(p)
    q₀, q₁ = pqs(q)
    p₀ * q₁ == q₀ * p₁ || return false
end
function ==(p::RationalTransferFunction{T,X,P,T1}, q::Polynomial{S,Y}) where {T,X,P,T1,S,Y}
    isconstant(p.num) && isconstant(p.den) && isconstant(q) && p(0) == q(0) && return true
    X == Y || return false
    p₀, p₁ = pqs(p)
    p₀ == p₁ * q || return false
end
==(p::Polynomial{S,Y}, q::RationalTransferFunction{T,X,P,T1}) where {S,Y,T,X,P,T1} = ==(q,p)

function Base.isapprox(pq₁::PQ₁, pq₂::PQ₂,
                       rtol::Real = sqrt(eps(float(real(promote_type(T,S))))),
                       atol::Real = zero(real(promote_type(T,S)))) where
                                {T,X,P,T1,PQ₁<:RationalTransferFunction{T,X,P,T1},
                                 S,Y,Q,T2,PQ₂<:RationalTransferFunction{S,Y,Q,T2}}

    p₁,q₁ = pqs(pq₁)
    p₂,q₂ = pqs(pq₂)

    return X == Y && T1 == T2 && isapprox(p₁*q₂, q₁*p₂; rtol=rtol, atol=atol)
end
function Base.isapprox(pq₁::PQ₁, q₂::Q₂,
                       rtol::Real = sqrt(eps(float(real(promote_type(T,S))))),
                       atol::Real = zero(real(promote_type(T,S)))) where
                        {T,X,P,T1,PQ₁<:RationalTransferFunction{T,X,P,T1}, S,Y,Q₂<:Polynomial{S,Y}}

    p₁,q₁ = pqs(pq₁)
    return X == Y && isapprox(p₁, q₁*q₂; rtol=rtol, atol=atol)
end
Base.isapprox(q₁::Q₁, pq₂::PQ₂,
              rtol::Real = sqrt(eps(float(real(promote_type(T,S))))),
              atol::Real = zero(real(promote_type(T,S)))) where
              {T,X,P,T2,PQ₂<:RationalTransferFunction{T,X,P,T2},S,Y,Q₁<:Polynomial{S,Y}} = isapprox(pq₂,q₁)
Base.isapprox(q₁::Number, pq₂::PQ₂,
              rtol::Real = sqrt(eps(float(real(promote_type(T,eltype(q₁)))))),
              atol::Real = zero(real(promote_type(T,eltype(q₁))))) where
              {T,X,P,T2,PQ₂<:RationalTransferFunction{T,X,P,T2}} = isapprox(pq₂,q₁*one(pq₂))
Base.isapprox(pq₂::PQ₂, q₁::Number,
              rtol::Real = sqrt(eps(float(real(promote_type(T,eltype(q₁)))))),
              atol::Real = zero(real(promote_type(T,eltype(q₁))))) where
              {T,X,P,T2,PQ₂<:RationalTransferFunction{T,X,P,T2}} = isapprox(pq₂,q₁*one(pq₂))

## various functions
"""
     normalize(r)

Normalize the rational transfer function `r(λ)` to have a monic denominator polynomial.
"""
function normalize(f::RationalTransferFunction)
    k = last(f.den.coeffs)
    rtf(f.num/k,f.den/k,Ts = f.Ts)
end
"""
     simplify(r; atol = 0, rtol = atol)

Simplify the rational transfer function `r(λ)` by cancellation of common divisors of numerator and denominator.
The keyword arguments `atol` and `rtol` are the absolute and relative tolerances for the nonzero
numerator and denominator coefficients.
"""
function simplify(r::RationalTransferFunction; atol::Real = 0, rtol::Real=10*eps(float(real(_eltype(r)))))
    m = degree(r.num)
    (m == 0 || degree(r.den) == 0) && (return r)
    pnum = r.num.coeffs
    pden = r.den.coeffs
    d, u, w,  = MatrixPencils.polgcdvw(pnum, pden, atol = atol, rtol = rtol, maxnit = 3)
    length(u) > m && (return r) # no cancelation, keep original r
    knum = pnum[end]/d[end]/u[end]
    kden = pden[end]/d[end]/w[end]
    return rtf(Polynomial(u*knum,Polynomials.indeterminate(r.num)), Polynomial(w*kden,Polynomials.indeterminate(r.den)), Ts = r.Ts)
end
"""
    Rval = evalfr(R,val)

Evaluate the rational transfer function matrix `R(λ)` for `λ = val`.
"""
function evalfr(R::VecOrMat{<:RationalTransferFunction}, val::Number)
    return pmeval(numpoly.(R),val) ./ pmeval(denpoly.(R),val)
end
"""
    Rval = evalfr(R; fval = 0)

Evaluate the rational transfer function matrix `R(λ)` for `λ = val`, where `val = im*fval`
for a continuous-time system or `val = exp(im*fval*abs(Ts))` for a discrete-time system,
with `Ts` the system sampling time.
"""
function evalfr(R::VecOrMat{<:RationalTransferFunction}; fval::Number = 0)
    length(R) > 0 ? Ts = R[1].Ts : Ts = 0
    val = Ts == 0 ? im*abs(fval) : exp(im*abs(fval*Ts))
    return pmeval(numpoly.(R),val) ./ pmeval(denpoly.(R),val)
end
dcgain(R::VecOrMat{<:RationalTransferFunction}) = _eltype(R) <: Complex ? evalfr(R) : real(evalfr(R))
"""
    rval = evalfr(r,val)

Evaluate the rational transfer function  `r(λ)` for `λ = val`.
"""
function evalfr(r::RationalTransferFunction,val::Number)
    return r.num(val) ./ r.den(val)
end
"""
    rval = evalfr(r; fval = 0)

Evaluate the rational transfer function  `r(λ)` for `λ = val`, where `val = im*fval`
for a continuous-time system or `val = exp(im*fval*Ts)` for a discrete-time system,
with `Ts` the system sampling time.
"""
function evalfr(r::RationalTransferFunction; fval::Number = 0)
    Ts = r.Ts
    val = Ts == 0 ? im*abs(fval) : exp(im*abs(fval*Ts))
    return r.num(val) ./ r.den(val)
end
dcgain(r::RationalTransferFunction) = _eltype(r) <: Complex ? evalfr(r) : real(evalfr(r))