text:
@testset "4.7.6 f^(a+b x+c x^2) trig(d+e x+f x^2)^n" begin
    (F, a, b, c, d, e, f, m, n, x, ) = @variables F a b c d e f m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*involving*products*of*exponentials*and*trig*functions=#

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

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

    @test_int [F^(c*(a + b*x))*sin(d + e*x)^n, x, 2, -((F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(-n, -((e*n + im*b*c*log(F))/(2*e)), (1/2)*(2 - n - (im*b*c*log(F))/e), ℯ^(2*im*(d + e*x)))*sin(d + e*x)^n)/((1 - ℯ^(2*im*(d + e*x)))^n*(im*e*n - b*c*log(F))))]

    @test_int [F^(c*(a + b*x))*sin(d + e*x)^3, x, 2, -((6*e^3*F^(c*(a + b*x))*cos(d + e*x))/(9*e^4 + 10*b^2*c^2*e^2*log(F)^2 + b^4*c^4*log(F)^4)) + (6*b*c*e^2*F^(c*(a + b*x))*log(F)*sin(d + e*x))/(9*e^4 + 10*b^2*c^2*e^2*log(F)^2 + b^4*c^4*log(F)^4) - (3*e*F^(c*(a + b*x))*cos(d + e*x)*sin(d + e*x)^2)/(9*e^2 + b^2*c^2*log(F)^2) + (b*c*F^(c*(a + b*x))*log(F)*sin(d + e*x)^3)/(9*e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*sin(d + e*x)^2, x, 2, (2*e^2*F^(c*(a + b*x)))/(b*c*log(F)*(4*e^2 + b^2*c^2*log(F)^2)) - (2*e*F^(c*(a + b*x))*cos(d + e*x)*sin(d + e*x))/(4*e^2 + b^2*c^2*log(F)^2) + (b*c*F^(c*(a + b*x))*log(F)*sin(d + e*x)^2)/(4*e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*sin(d + e*x)^1, x, 1, -((e*F^(c*(a + b*x))*cos(d + e*x))/(e^2 + b^2*c^2*log(F)^2)) + (b*c*F^(c*(a + b*x))*log(F)*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*csc(d + e*x)^1, x, 1, -((2*ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, (e - im*b*c*log(F))/(2*e), (1/2)*(3 - (im*b*c*log(F))/e), ℯ^(2*im*(d + e*x))))/(e - im*b*c*log(F)))]
    @test_int [F^(c*(a + b*x))*csc(d + e*x)^2, x, 1, -((4*ℯ^(2*im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/(2*e), 2 - (im*b*c*log(F))/(2*e), ℯ^(2*im*(d + e*x))))/(2*im*e + b*c*log(F)))]
    @test_int [F^(c*(a + b*x))*csc(d + e*x)^3, x, 2, -((F^(c*(a + b*x))*cot(d + e*x)*csc(d + e*x))/(2*e)) - (b*c*F^(c*(a + b*x))*csc(d + e*x)*log(F))/(2*e^2) - (ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, (e - im*b*c*log(F))/(2*e), (1/2)*(3 - (im*b*c*log(F))/e), ℯ^(2*im*(d + e*x)))*(e + im*b*c*log(F)))/e^2]
    @test_int [F^(c*(a + b*x))*csc(d + e*x)^4, x, 2, -((F^(c*(a + b*x))*cot(d + e*x)*csc(d + e*x)^2)/(3*e)) - (b*c*F^(c*(a + b*x))*csc(d + e*x)^2*log(F))/(6*e^2) + (2*ℯ^(2*im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/(2*e), 2 - (im*b*c*log(F))/(2*e), ℯ^(2*im*(d + e*x)))*(2*im*e - b*c*log(F)))/(3*e^2)]

    @test_int [ℯ^x*sin(x)^4, x, 3, (24*ℯ^x)/85 - (24/85)*ℯ^x*cos(x)*sin(x) + (12/85)*ℯ^x*sin(x)^2 - (4/17)*ℯ^x*cos(x)*sin(x)^3 + (1/17)*ℯ^x*sin(x)^4]

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

    @test_int [F^(c*(a + b*x))*cos(d + e*x)^n, x, 2, -((F^(c*(a + b*x))*cos(d + e*x)^n*HypergeometricFunctions._₂F₁(-n, -((e*n + im*b*c*log(F))/(2*e)), (1/2)*(2 - n - (im*b*c*log(F))/e), -ℯ^(2*im*(d + e*x))))/((1 + ℯ^(2*im*(d + e*x)))^n*(im*e*n - b*c*log(F))))]

    @test_int [F^(c*(a + b*x))*cos(d + e*x)^3, x, 2, (b*c*F^(c*(a + b*x))*cos(d + e*x)^3*log(F))/(9*e^2 + b^2*c^2*log(F)^2) + (6*b*c*e^2*F^(c*(a + b*x))*cos(d + e*x)*log(F))/(9*e^4 + 10*b^2*c^2*e^2*log(F)^2 + b^4*c^4*log(F)^4) + (3*e*F^(c*(a + b*x))*cos(d + e*x)^2*sin(d + e*x))/(9*e^2 + b^2*c^2*log(F)^2) + (6*e^3*F^(c*(a + b*x))*sin(d + e*x))/(9*e^4 + 10*b^2*c^2*e^2*log(F)^2 + b^4*c^4*log(F)^4)]
    @test_int [F^(c*(a + b*x))*cos(d + e*x)^2, x, 2, (2*e^2*F^(c*(a + b*x)))/(b*c*log(F)*(4*e^2 + b^2*c^2*log(F)^2)) + (b*c*F^(c*(a + b*x))*cos(d + e*x)^2*log(F))/(4*e^2 + b^2*c^2*log(F)^2) + (2*e*F^(c*(a + b*x))*cos(d + e*x)*sin(d + e*x))/(4*e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*cos(d + e*x)^1, x, 1, (b*c*F^(c*(a + b*x))*cos(d + e*x)*log(F))/(e^2 + b^2*c^2*log(F)^2) + (e*F^(c*(a + b*x))*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*sec(d + e*x)^1, x, 1, (2*ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, (e - im*b*c*log(F))/(2*e), (1/2)*(3 - (im*b*c*log(F))/e), -ℯ^(2*im*(d + e*x))))/(im*e + b*c*log(F))]
    @test_int [F^(c*(a + b*x))*sec(d + e*x)^2, x, 1, (4*ℯ^(2*im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/(2*e), 2 - (im*b*c*log(F))/(2*e), -ℯ^(2*im*(d + e*x))))/(2*im*e + b*c*log(F))]
    @test_int [F^(c*(a + b*x))*sec(d + e*x)^3, x, 2, -((ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, (e - im*b*c*log(F))/(2*e), (1/2)*(3 - (im*b*c*log(F))/e), -ℯ^(2*im*(d + e*x)))*(im*e - b*c*log(F)))/e^2) - (b*c*F^(c*(a + b*x))*log(F)*sec(d + e*x))/(2*e^2) + (F^(c*(a + b*x))*sec(d + e*x)*tan(d + e*x))/(2*e)]
    @test_int [F^(c*(a + b*x))*sec(d + e*x)^4, x, 2, -((2*ℯ^(2*im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/(2*e), 2 - (im*b*c*log(F))/(2*e), -ℯ^(2*im*(d + e*x)))*(2*im*e - b*c*log(F)))/(3*e^2)) - (b*c*F^(c*(a + b*x))*log(F)*sec(d + e*x)^2)/(6*e^2) + (F^(c*(a + b*x))*sec(d + e*x)^2*tan(d + e*x))/(3*e)]

    @test_int [ℯ^x*cos(x)^4, x, 3, (24*ℯ^x)/85 + (12/85)*ℯ^x*cos(x)^2 + (1/17)*ℯ^x*cos(x)^4 + (24/85)*ℯ^x*cos(x)*sin(x) + (4/17)*ℯ^x*cos(x)^3*sin(x)]

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

    @test_int [ℯ^(c*(a + b*x))*tan(d + e*x)^3, x, 6, (im*ℯ^(c*(a + b*x)))/(b*c) - (6*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), -ℯ^(2*im*(d + e*x))))/(b*c) + (12*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), -ℯ^(2*im*(d + e*x))))/(b*c) - (8*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(3, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), -ℯ^(2*im*(d + e*x))))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))*tan(d + e*x)^2, x, 5, -(ℯ^(c*(a + b*x))/(b*c)) + (4*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), -ℯ^(2*im*(d + e*x))))/(b*c) - (4*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), -ℯ^(2*im*(d + e*x))))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))*tan(d + e*x)^1, x, 4, -((im*ℯ^(c*(a + b*x)))/(b*c)) + (2*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), -ℯ^(2*im*(d + e*x))))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))*cot(d + e*x)^1, x, 4, (im*ℯ^(c*(a + b*x)))/(b*c) - (2*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), ℯ^(2*im*(d + e*x))))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))*cot(d + e*x)^2, x, 5, -(ℯ^(c*(a + b*x))/(b*c)) + (4*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), ℯ^(2*im*(d + e*x))))/(b*c) - (4*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), ℯ^(2*im*(d + e*x))))/(b*c)]
    @test_int [ℯ^(c*(a + b*x))*cot(d + e*x)^3, x, 6, -((im*ℯ^(c*(a + b*x)))/(b*c)) + (6*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(1, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), ℯ^(2*im*(d + e*x))))/(b*c) - (12*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), ℯ^(2*im*(d + e*x))))/(b*c) + (8*im*ℯ^(c*(a + b*x))*HypergeometricFunctions._₂F₁(3, -((im*b*c)/(2*e)), 1 - (im*b*c)/(2*e), ℯ^(2*im*(d + e*x))))/(b*c)]

    @test_int [F^(a + b*x)*tan(pi/4 - (c + d*x)/2), x, 5, (im*F^(a + b*x))/(b*log(F)) - (2*im*F^(a + b*x)*HypergeometricFunctions._₂F₁(1, -((im*b*log(F))/d), 1 - (im*b*log(F))/d, im*ℯ^(im*(c + d*x))))/(b*log(F))]

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

    @test_int [F^(c*(a + b*x))*sec(d + e*x)^n, x, 2, ((1 + ℯ^(2*im*(d + e*x)))^n*F^(a*c + b*c*x)*HypergeometricFunctions._₂F₁(n, (e*n - im*b*c*log(F))/(2*e), (1/2)*(2 + n - (im*b*c*log(F))/e), -ℯ^(2*im*(d + e*x)))*sec(d + e*x)^n)/(im*e*n + b*c*log(F))]

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

    @test_int [F^(c*(a + b*x))*csc(d + e*x)^n, x, 2, -(((1 - ℯ^(-2*im*(d + e*x)))^n*F^(a*c + b*c*x)*csc(d + e*x)^n*HypergeometricFunctions._₂F₁(n, (e*n + im*b*c*log(F))/(2*e), (1/2)*(2 + n + (im*b*c*log(F))/e), ℯ^(-2*im*(d + e*x))))/(im*e*n - b*c*log(F)))]

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

    #= [(f*x)^m*F^(c*(a + b*x))*sin(d + e*x)^3, x, 10, -((3*F^(a*c)*(f*x)^m*Gamma(1 + m, x*(im*e - b*c*log(F))))/(ℯ^(im*d)*(x*(im*e - b*c*log(F)))^m*(8*(e + im*b*c*log(F))))) + (F^(a*c)*(f*x)^m*Gamma(1 + m, x*(3*im*e - b*c*log(F))))/(ℯ^(3*im*d)*(x*(3*im*e - b*c*log(F)))^m*(8*(3*e + im*b*c*log(F)))) - (3*ℯ^(im*d)*F^(a*c)*(f*x)^m*Gamma(1 + m, (-x)*(im*e + b*c*log(F))))/(((-x)*(im*e + b*c*log(F)))^m*(8*(e - im*b*c*log(F)))) + (ℯ^(3*im*d)*F^(a*c)*(f*x)^m*Gamma(1 + m, (-x)*(3*im*e + b*c*log(F))))/(((-x)*(3*im*e + b*c*log(F)))^m*(8*(3*e - im*b*c*log(F))))]
    @test_int [(f*x)^m*F^(c*(a + b*x))*sin(d + e*x)^2, x, 7, (F^(a*c)*(f*x)^m*Gamma(1 + m, (-b)*c*x*log(F)))/(((-b)*c*x*log(F))^m*(2*b*c*log(F))) + (F^(a*c)*(f*x)^m*Gamma(1 + m, x*(2*im*e - b*c*log(F))))/(ℯ^(2*im*d)*(x*(2*im*e - b*c*log(F)))^m*(4*(2*im*e - b*c*log(F)))) - (ℯ^(2*im*d)*F^(a*c)*(f*x)^m*Gamma(1 + m, (-x)*(2*im*e + b*c*log(F))))/(((-x)*(2*im*e + b*c*log(F)))^m*(4*(2*im*e + b*c*log(F))))] =#
    @test_int [(f*x)^m*F^(c*(a + b*x))*sin(d + e*x)^1, x, -1, -((F^(a*c)*(f*x)^m*Gamma(1 + m, x*(im*e - b*c*log(F))))/(ℯ^(im*d)*(x*(im*e - b*c*log(F)))^m*(2*(e + im*b*c*log(F))))) - (ℯ^(im*d)*F^(a*c)*(f*x)^m*Gamma(1 + m, (-x)*(im*e + b*c*log(F))))/(((-x)*(im*e + b*c*log(F)))^m*(2*(e - im*b*c*log(F))))]
    @test_int [(f*x)^m*F^(c*(a + b*x))/sin(d + e*x)^1, x, 1, CannotIntegrate(F^(a*c + b*c*x)*(f*x)^m*csc(d + e*x), x)]
    @test_int [(f*x)^m*F^(c*(a + b*x))/sin(d + e*x)^2, x, 1, CannotIntegrate(F^(a*c + b*c*x)*(f*x)^m*csc(d + e*x)^2, x)]

    @test_int [f*F^(c*(a + b*x))*(f*x)^(-2 + m)*(e*x*cos(d + e*x) + (-1 + m + b*c*x*log(F))*sin(d + e*x)), x, 11, F^(a*c + b*c*x)*(f*x)^(-1 + m)*sin(d + e*x)]
    @test_int [f*F^(c*(a + b*x))*(f*x)^m*(e*x*cos(d + e*x) + (1 + m + b*c*x*log(F))*sin(d + e*x)), x, -6, f*F^(c*(a + b*x))*x*(f*x)^m*sin(d + e*x)]
    @test_int [(F^(c*(a + b*x))*(f*x)^m*(e*x*cos(d + e*x) + (m + b*c*x*log(F))*sin(d + e*x)))/x, x, 7, F^(a*c + b*c*x)*(f*x)^m*sin(d + e*x)]

    #= [F^(c*(a + b*x))*x*(e*x*cos(d + e*x) + (2 + b*c*x*log(F))*sin(d + e*x)), x, 0, F^(c*(a + b*x))*x^2*sin(d + e*x)] =#
    @test_int [F^(c*(a + b*x))*(e*x*cos(d + e*x) + (1 + b*c*x*log(F))*sin(d + e*x)), x, 14, F^(c*(a + b*x))*x*sin(d + e*x), (e^3*F^(a*c + b*c*x)*cos(d + e*x))/(e^2 + b^2*c^2*log(F)^2)^2 + (b^2*c^2*e*F^(a*c + b*c*x)*cos(d + e*x)*log(F)^2)/(e^2 + b^2*c^2*log(F)^2)^2 - (e*F^(a*c + b*c*x)*cos(d + e*x))/(e^2 + b^2*c^2*log(F)^2) - (b*c*e^2*F^(a*c + b*c*x)*log(F)*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2)^2 - (b^3*c^3*F^(a*c + b*c*x)*log(F)^3*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2)^2 + (e^2*F^(a*c + b*c*x)*x*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2) + (b*c*F^(a*c + b*c*x)*log(F)*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2) + (b^2*c^2*F^(a*c + b*c*x)*x*log(F)^2*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*(e*cos(d + e*x) + b*c*log(F)*sin(d + e*x)), x, 1, F^(c*(a + b*x))*sin(d + e*x)]
    @test_int [(F^(c*(a + b*x))*(e*x*cos(d + e*x) + (-1 + b*c*x*log(F))*sin(d + e*x)))/x^2, x, 6, (F^(a*c + b*c*x)*sin(d + e*x))/x]
    @test_int [(F^(c*(a + b*x))*(e*x*cos(d + e*x) + (-2 + b*c*x*log(F))*sin(d + e*x)))/x^3, x, 10, (F^(a*c + b*c*x)*sin(d + e*x))/x^2]

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

    @test_int [ℯ^(a + b*x)*cos(c + d*x)^1*sin(c + d*x)^1, x, 3, -((d*ℯ^(a + b*x)*cos(2*c + 2*d*x))/(b^2 + 4*d^2)) + (b*ℯ^(a + b*x)*sin(2*c + 2*d*x))/(2*(b^2 + 4*d^2))]
    @test_int [ℯ^(a + b*x)*cos(c + d*x)^1*sin(c + d*x)^2, x, 4, (b*ℯ^(a + b*x)*cos(c + d*x))/(4*(b^2 + d^2)) - (b*ℯ^(a + b*x)*cos(3*c + 3*d*x))/(4*(b^2 + 9*d^2)) + (d*ℯ^(a + b*x)*sin(c + d*x))/(4*(b^2 + d^2)) - (3*d*ℯ^(a + b*x)*sin(3*c + 3*d*x))/(4*(b^2 + 9*d^2))]
    @test_int [ℯ^(a + b*x)*cos(c + d*x)^1*sin(c + d*x)^3, x, 4, -((d*ℯ^(a + b*x)*cos(2*c + 2*d*x))/(2*(b^2 + 4*d^2))) + (d*ℯ^(a + b*x)*cos(4*c + 4*d*x))/(2*(b^2 + 16*d^2)) + (b*ℯ^(a + b*x)*sin(2*c + 2*d*x))/(4*(b^2 + 4*d^2)) - (b*ℯ^(a + b*x)*sin(4*c + 4*d*x))/(8*(b^2 + 16*d^2))]

    @test_int [ℯ^(a + b*x)*cos(c + d*x)^2*sin(c + d*x)^1, x, 4, -((d*ℯ^(a + b*x)*cos(c + d*x))/(4*(b^2 + d^2))) - (3*d*ℯ^(a + b*x)*cos(3*c + 3*d*x))/(4*(b^2 + 9*d^2)) + (b*ℯ^(a + b*x)*sin(c + d*x))/(4*(b^2 + d^2)) + (b*ℯ^(a + b*x)*sin(3*c + 3*d*x))/(4*(b^2 + 9*d^2))]
    @test_int [ℯ^(a + b*x)*cos(c + d*x)^2*sin(c + d*x)^2, x, 4, ℯ^(a + b*x)/(8*b) - (b*ℯ^(a + b*x)*cos(4*c + 4*d*x))/(8*(b^2 + 16*d^2)) - (d*ℯ^(a + b*x)*sin(4*c + 4*d*x))/(2*(b^2 + 16*d^2))]
    @test_int [ℯ^(a + b*x)*cos(c + d*x)^2*sin(c + d*x)^3, x, 5, -((d*ℯ^(a + b*x)*cos(c + d*x))/(8*(b^2 + d^2))) - (3*d*ℯ^(a + b*x)*cos(3*c + 3*d*x))/(16*(b^2 + 9*d^2)) + (5*d*ℯ^(a + b*x)*cos(5*c + 5*d*x))/(16*(b^2 + 25*d^2)) + (b*ℯ^(a + b*x)*sin(c + d*x))/(8*(b^2 + d^2)) + (b*ℯ^(a + b*x)*sin(3*c + 3*d*x))/(16*(b^2 + 9*d^2)) - (b*ℯ^(a + b*x)*sin(5*c + 5*d*x))/(16*(b^2 + 25*d^2))]

    @test_int [ℯ^(a + b*x)*cos(c + d*x)^3*sin(c + d*x)^1, x, 4, -((d*ℯ^(a + b*x)*cos(2*c + 2*d*x))/(2*(b^2 + 4*d^2))) - (d*ℯ^(a + b*x)*cos(4*c + 4*d*x))/(2*(b^2 + 16*d^2)) + (b*ℯ^(a + b*x)*sin(2*c + 2*d*x))/(4*(b^2 + 4*d^2)) + (b*ℯ^(a + b*x)*sin(4*c + 4*d*x))/(8*(b^2 + 16*d^2))]
    @test_int [ℯ^(a + b*x)*cos(c + d*x)^3*sin(c + d*x)^2, x, 5, (b*ℯ^(a + b*x)*cos(c + d*x))/(8*(b^2 + d^2)) - (b*ℯ^(a + b*x)*cos(3*c + 3*d*x))/(16*(b^2 + 9*d^2)) - (b*ℯ^(a + b*x)*cos(5*c + 5*d*x))/(16*(b^2 + 25*d^2)) + (d*ℯ^(a + b*x)*sin(c + d*x))/(8*(b^2 + d^2)) - (3*d*ℯ^(a + b*x)*sin(3*c + 3*d*x))/(16*(b^2 + 9*d^2)) - (5*d*ℯ^(a + b*x)*sin(5*c + 5*d*x))/(16*(b^2 + 25*d^2))]
    @test_int [ℯ^(a + b*x)*cos(c + d*x)^3*sin(c + d*x)^3, x, 4, -((3*d*ℯ^(a + b*x)*cos(2*c + 2*d*x))/(16*(b^2 + 4*d^2))) + (3*d*ℯ^(a + b*x)*cos(6*c + 6*d*x))/(16*(b^2 + 36*d^2)) + (3*b*ℯ^(a + b*x)*sin(2*c + 2*d*x))/(32*(b^2 + 4*d^2)) - (b*ℯ^(a + b*x)*sin(6*c + 6*d*x))/(32*(b^2 + 36*d^2))]

    @test_int [ℯ^x*x*sin(x), x, 4, (1/2)*ℯ^x*cos(x) - (1/2)*ℯ^x*x*cos(x) + (1/2)*ℯ^x*x*sin(x)]
    @test_int [ℯ^x*x^2*sin(x), x, 11, (-(1/2))*ℯ^x*cos(x) + ℯ^x*x*cos(x) - (1/2)*ℯ^x*x^2*cos(x) - (1/2)*ℯ^x*sin(x) + (1/2)*ℯ^x*x^2*sin(x)]

    @test_int [ℯ^x*x*cos(x), x, 4, (1/2)*ℯ^x*x*cos(x) - (1/2)*ℯ^x*sin(x) + (1/2)*ℯ^x*x*sin(x)]
    @test_int [ℯ^x*x^2*cos(x), x, 11, (-(1/2))*ℯ^x*cos(x) + (1/2)*ℯ^x*x^2*cos(x) + (1/2)*ℯ^x*sin(x) - ℯ^x*x*sin(x) + (1/2)*ℯ^x*x^2*sin(x)]

    @test_int [ℯ^(3*x)*(-5*cos(4*x) + 2*sin(4*x)), x, 4, (-(23/25))*ℯ^(3*x)*cos(4*x) - (14/25)*ℯ^(3*x)*sin(4*x)]

    @test_int [sin(x)/ℯ^x + ℯ^x*sin(x), x, 3, ((-(1/2))*cos(x))/ℯ^x - (1/2)*ℯ^x*cos(x) - ((1/2)*sin(x))/ℯ^x + (1/2)*ℯ^x*sin(x)]

    @test_int [F^(a + b*x)*(cos(c + d*x)/(e + e*sin(c + d*x))), x, 5, (im*F^(a + b*x))/(b*e*log(F)) - (2*im*F^(a + b*x)*HypergeometricFunctions._₂F₁(1, -((im*b*log(F))/d), 1 - (im*b*log(F))/d, im*ℯ^(im*(c + d*x))))/(b*e*log(F))]
    @test_int [F^(a + b*x)*(cos(c + d*x)/(e - e*sin(c + d*x))), x, 5, -((im*F^(a + b*x))/(b*e*log(F))) + (2*im*F^(a + b*x)*HypergeometricFunctions._₂F₁(1, -((im*b*log(F))/d), 1 - (im*b*log(F))/d, (-im)*ℯ^(im*(c + d*x))))/(b*e*log(F))]
    @test_int [F^(a + b*x)*(sin(c + d*x)/(e + e*cos(c + d*x))), x, 5, -((im*F^(a + b*x))/(b*e*log(F))) + (2*im*F^(a + b*x)*HypergeometricFunctions._₂F₁(1, -((im*b*log(F))/d), 1 - (im*b*log(F))/d, -ℯ^(im*(c + d*x))))/(b*e*log(F))]
    @test_int [F^(a + b*x)*(sin(c + d*x)/(e - e*cos(c + d*x))), x, 5, (im*F^(a + b*x))/(b*e*log(F)) - (2*im*F^(a + b*x)*HypergeometricFunctions._₂F₁(1, -((im*b*log(F))/d), 1 - (im*b*log(F))/d, ℯ^(im*(c + d*x))))/(b*e*log(F))]

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

    @test_int [ℯ^(x^2)*sin(b*x), x, 6, (1/4)*im*ℯ^(b^2/4)*sqrt(pi)*erfi((1/2)*((-im)*b + 2*x)) - (1/4)*im*ℯ^(b^2/4)*sqrt(pi)*erfi((1/2)*(im*b + 2*x))]
    @test_int [ℯ^(x^2)*cos(b*x), x, 6, (1/4)*ℯ^(b^2/4)*sqrt(pi)*erfi((1/2)*((-im)*b + 2*x)) + (1/4)*ℯ^(b^2/4)*sqrt(pi)*erfi((1/2)*(im*b + 2*x))]
    @test_int [ℯ^(x^2)*sin(a + b*x), x, 6, (1/4)*im*ℯ^((-im)*a + b^2/4)*sqrt(pi)*erfi((1/2)*((-im)*b + 2*x)) - (1/4)*im*ℯ^(im*a + b^2/4)*sqrt(pi)*erfi((1/2)*(im*b + 2*x))]
    @test_int [ℯ^(x^2)*cos(a + b*x), x, 6, (1/4)*ℯ^((-im)*a + b^2/4)*sqrt(pi)*erfi((1/2)*((-im)*b + 2*x)) + (1/4)*ℯ^(im*a + b^2/4)*sqrt(pi)*erfi((1/2)*(im*b + 2*x))]

    @test_int [ℯ^(2*x^2)*x*cos(2*x^2), x, 2, (1/8)*ℯ^(2*x^2)*cos(2*x^2) + (1/8)*ℯ^(2*x^2)*sin(2*x^2)]

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

    @test_int [ℯ^x*sin(ℯ^x), x, 2, -cos(ℯ^x)]

    @test_int [ℯ^x*csc(ℯ^x)*sec(ℯ^x), x, 3, log(tan(ℯ^x))]
    @test_int [ℯ^x*cos(ℯ^x), x, 2, sin(ℯ^x)]
    @test_int [ℯ^(2*x)*cos(ℯ^(2*x)), x, 2, sin(ℯ^(2*x))/2]
    @test_int [cos(ℯ^(-2*x))/ℯ^(2*x), x, 2, -sin(ℯ^(-2*x))/2]
    @test_int [ℯ^(2*x)*cos(ℯ^x), x, 3, cos(ℯ^x) + ℯ^x*sin(ℯ^x)]
    @test_int [ℯ^(-1 + 3*x)*cos(ℯ^(-1 + 3*x))*sin(1 + ℯ^(-1 + 3*x)), x, 4, -cos(1 + 2*ℯ^(-1 + 3*x))/12 + (ℯ^(-1 + 3*x)*sin(1))/6]

    @test_int [ℯ^x*tan(ℯ^x), x, 2, -log(cos(ℯ^x))]

    @test_int [ℯ^x*sec(ℯ^x), x, 2, atanh(sin(ℯ^x))]
    @test_int [ℯ^x*sec(ℯ^x)*tan(ℯ^x), x, 3, sec(ℯ^x)]

    @test_int [ℯ^x*csc(ℯ^x)^2, x, 3, -cot(ℯ^x)]

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

    @test_int [ℯ^x*sin(a + b*x), x, 1, -((b*ℯ^x*cos(a + b*x))/(1 + b^2)) + (ℯ^x*sin(a + b*x))/(1 + b^2)]
    @test_int [ℯ^x*sin(a + c*x^2), x, 6, ((-1)^(3/4)*ℯ^((im/4)*(4*a + c^(-1)))*sqrt(pi)*erf(((-1)^(1/4)*(1 + (2*im)*c*x))/(2*sqrt(c))))/(4*sqrt(c)) + ((-1)^(3/4)*sqrt(pi)*erfi(((-1)^(1/4)*(1 - (2*im)*c*x))/(2*sqrt(c))))/(4*sqrt(c)*ℯ^((im/4)*(4*a + c^(-1))))]
    @test_int [ℯ^x*sin(a + b*x + c*x^2), x, 6, ((-1)^(3/4)*ℯ^((1/4)*im*(4*a + (1 + im*b)^2/c))*sqrt(pi)*erf(((-1)^(1/4)*(1 + im*b + 2*im*c*x))/(2*sqrt(c))))/(4*sqrt(c)) + ((-1)^(3/4)*ℯ^((-im)*a + (im*(im + b)^2)/(4*c))*sqrt(pi)*erfi(((-1)^(1/4)*(1 - im*b - 2*im*c*x))/(2*sqrt(c))))/(4*sqrt(c))]

    @test_int [ℯ^x^2*sin(a + b*x), x, 6, (im/4)*ℯ^((-im)*a + b^2/4)*sqrt(pi)*erfi(((-im)*b + 2*x)/2) - (im/4)*ℯ^(im*a + b^2/4)*sqrt(pi)*erfi((im*b + 2*x)/2)]
    @test_int [ℯ^x^2*sin(a + c*x^2), x, 4, ((im/4)*sqrt(pi)*erfi(sqrt(1 - im*c)*x))/(sqrt(1 - im*c)*ℯ^(im*a)) - ((im/4)*ℯ^(im*a)*sqrt(pi)*erfi(sqrt(1 + im*c)*x))/sqrt(1 + im*c)]
    @test_int [ℯ^x^2*sin(a + b*x + c*x^2), x, 6, -((im*sqrt(pi)*erfi((im*b - 2*(1 - im*c)*x)/(2*sqrt(1 - im*c))))/(ℯ^(im*(a - b^2/(4*im + 4*c)))*(4*sqrt(1 - im*c)))) - (im*ℯ^(im*a + b^2/(4*(1 + im*c)))*sqrt(pi)*erfi((im*b + 2*(1 + im*c)*x)/(2*sqrt(1 + im*c))))/(4*sqrt(1 + im*c))]

    @test_int [f^(a + b*x)*sin(d + f*x^2), x, 8, ((-1)^(3/4)*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*((2*im)*f*x + b*log(f)))/(2*sqrt(f))))/4 - ((-1)^(3/4)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*((2*im)*f*x - b*log(f)))/(2*sqrt(f))))/(4*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f)))]
    @test_int [f^(a + b*x)*sin(d + f*x^2)^2, x, 9, (1/16 + im/16)*ℯ^((2*im)*d + ((im/8)*b^2*log(f)^2)/f)*f^(-1/2 + a)*sqrt(pi)*erf(((1/4 + im/4)*((4*im)*f*x + b*log(f)))/sqrt(f)) + ((1/16 + im/16)*f^(-1/2 + a)*sqrt(pi)*erfi(((1/4 + im/4)*((4*im)*f*x - b*log(f)))/sqrt(f)))/ℯ^((im/8)*(16*d + (b^2*log(f)^2)/f)) + f^(a + b*x)/(2*b*log(f))]
    @test_int [f^(a + b*x)*sin(d + f*x^2)^3, x, 14, (3*(-1)^(3/4)*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*((2*im)*f*x + b*log(f)))/(2*sqrt(f))))/16 + (1/16 - im/16)*ℯ^((3*im)*d + ((im/12)*b^2*log(f)^2)/f)*f^(-1/2 + a)*sqrt(pi/6)*erf(((1/2 + im/2)*((6*im)*f*x + b*log(f)))/(sqrt(6)*sqrt(f))) - (3*(-1)^(3/4)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*((2*im)*f*x - b*log(f)))/(2*sqrt(f))))/(16*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f))) - ((1/16 - im/16)*f^(-1/2 + a)*sqrt(pi/6)*erfi(((1/2 + im/2)*((6*im)*f*x - b*log(f)))/(sqrt(6)*sqrt(f))))/ℯ^((im/12)*(36*d + (b^2*log(f)^2)/f))]

    @test_int [f^(a + b*x)*sin(d + e*x + f*x^2), x, 8, ((-1)^(3/4)*ℯ^((im/4)*(4*d + (im*e + b*log(f))^2/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*(im*e + (2*im)*f*x + b*log(f)))/(2*sqrt(f))))/4 - ((-1)^(3/4)*ℯ^((-im)*d + ((im/4)*(e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*(im*e + (2*im)*f*x - b*log(f)))/(2*sqrt(f))))/4]
    @test_int [f^(a + b*x)*sin(d + e*x + f*x^2)^2, x, 9, (1/16 + im/16)*ℯ^((2*im)*d + ((im/8)*((2*im)*e + b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erf(((1/4 + im/4)*((2*im)*e + (4*im)*f*x + b*log(f)))/sqrt(f)) + (1/16 + im/16)*ℯ^((-2*im)*d + ((im/8)*(2*e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erfi(((1/4 + im/4)*((2*im)*e + (4*im)*f*x - b*log(f)))/sqrt(f)) + f^(a + b*x)/(2*b*log(f))]
    @test_int [f^(a + b*x)*sin(d + e*x + f*x^2)^3, x, 14, (3*(-1)^(3/4)*ℯ^((im/4)*(4*d + (im*e + b*log(f))^2/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*(im*e + (2*im)*f*x + b*log(f)))/(2*sqrt(f))))/16 + (1/16 - im/16)*ℯ^((3*im)*d + ((im/12)*((3*im)*e + b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi/6)*erf(((1/2 + im/2)*((3*im)*e + (6*im)*f*x + b*log(f)))/(sqrt(6)*sqrt(f))) - (3*(-1)^(3/4)*ℯ^((-im)*d + ((im/4)*(e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*(im*e + (2*im)*f*x - b*log(f)))/(2*sqrt(f))))/16 - (1/16 - im/16)*ℯ^((-3*im)*d + ((im/12)*(3*e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi/6)*erfi(((1/2 + im/2)*((3*im)*e + (6*im)*f*x - b*log(f)))/(sqrt(6)*sqrt(f)))]

    @test_int [f^(a + c*x^2)*sin(d + e*x), x, 8, If(13>=8, ((-im/4)*ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) - ((im/4)*ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))), (im*ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(4*sqrt(c)*sqrt(log(f))) - (im*ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + c*x^2)*sin(d + e*x)^2, x, 9, If(13>=8, (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^((-2*im)*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi((im*e - c*x*log(f))/(sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))) - (ℯ^((2*im)*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi((im*e + c*x*log(f))/(sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))), (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) - (ℯ^(-2*im*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - c*x*log(f))/(sqrt(c)*sqrt(log(f))))))/(8*sqrt(c)*sqrt(log(f))) - (ℯ^(2*im*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi((im*e + c*x*log(f))/(sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + c*x^2)*sin(d + e*x)^3, x, 14, If(13>=8, (((-3*im)/16)*ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) + ((im/16)*ℯ^((-3*im)*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) - (((3*im)/16)*ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) + ((im/16)*ℯ^((3*im)*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))), (3*im*ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) - (im*ℯ^(-3*im*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((3*im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) - (3*im*ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (im*ℯ^(3*im*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi((3*im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))))]

    @test_int [f^(a + c*x^2)*sin(d + f*x^2), x, 6, ((im/4)*f^a*sqrt(pi)*erf(x*sqrt(im*f - c*log(f))))/(ℯ^(im*d)*sqrt(im*f - c*log(f))) - ((im/4)*ℯ^(im*d)*f^a*sqrt(pi)*erfi(x*sqrt(im*f + c*log(f))))/sqrt(im*f + c*log(f))]
    @test_int [f^(a + c*x^2)*sin(d + f*x^2)^2, x, 7, (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) - (f^a*sqrt(pi)*erf(x*sqrt((2*im)*f - c*log(f))))/(8*ℯ^((2*im)*d)*sqrt((2*im)*f - c*log(f))) - (ℯ^((2*im)*d)*f^a*sqrt(pi)*erfi(x*sqrt((2*im)*f + c*log(f))))/(8*sqrt((2*im)*f + c*log(f)))]
    @test_int [f^(a + c*x^2)*sin(d + f*x^2)^3, x, 10, (((3*im)/16)*f^a*sqrt(pi)*erf(x*sqrt(im*f - c*log(f))))/(ℯ^(im*d)*sqrt(im*f - c*log(f))) - ((im/16)*f^a*sqrt(pi)*erf(x*sqrt((3*im)*f - c*log(f))))/(ℯ^((3*im)*d)*sqrt((3*im)*f - c*log(f))) - (((3*im)/16)*ℯ^(im*d)*f^a*sqrt(pi)*erfi(x*sqrt(im*f + c*log(f))))/sqrt(im*f + c*log(f)) + ((im/16)*ℯ^((3*im)*d)*f^a*sqrt(pi)*erfi(x*sqrt((3*im)*f + c*log(f))))/sqrt((3*im)*f + c*log(f))]

    @test_int [f^(a + c*x^2)*sin(d + e*x + f*x^2), x, 8, (im*ℯ^((-im)*d - e^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(4*sqrt(im*f - c*log(f))) - (im*ℯ^(im*d + e^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(4*sqrt(im*f + c*log(f)))]
    @test_int [f^(a + c*x^2)*sin(d + e*x + f*x^2)^2, x, 9, (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) - (ℯ^(-2*im*d - e^2/(2*im*f - c*log(f)))*f^a*sqrt(pi)*erf((im*e + x*(2*im*f - c*log(f)))/sqrt(2*im*f - c*log(f))))/(8*sqrt(2*im*f - c*log(f))) - (ℯ^(2*im*d + e^2/(2*im*f + c*log(f)))*f^a*sqrt(pi)*erfi((im*e + x*(2*im*f + c*log(f)))/sqrt(2*im*f + c*log(f))))/(8*sqrt(2*im*f + c*log(f)))]
    @test_int [f^(a + c*x^2)*sin(d + e*x + f*x^2)^3, x, 14, (3*im*ℯ^((-im)*d - e^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(16*sqrt(im*f - c*log(f))) - (im*ℯ^(-3*im*d - (9*e^2)/(4*(3*im*f - c*log(f))))*f^a*sqrt(pi)*erf((3*im*e + 2*x*(3*im*f - c*log(f)))/(2*sqrt(3*im*f - c*log(f)))))/(16*sqrt(3*im*f - c*log(f))) - (3*im*ℯ^(im*d + e^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(16*sqrt(im*f + c*log(f))) + (im*ℯ^(3*im*d + (9*e^2)/(4*(3*im*f + c*log(f))))*f^a*sqrt(pi)*erfi((3*im*e + 2*x*(3*im*f + c*log(f)))/(2*sqrt(3*im*f + c*log(f)))))/(16*sqrt(3*im*f + c*log(f)))]

    @test_int [f^(a + b*x + c*x^2)*sin(d + e*x), x, 8, If(13>=8, ((-im/4)*ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) - ((im/4)*ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))), (im*ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(4*sqrt(c)*sqrt(log(f))) - (im*ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + b*x + c*x^2)*sin(d + e*x)^2, x, 10, If(13>=8, (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^((-2*im)*d + (2*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((2*im)*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))) - (ℯ^((2*im)*d - ((2*im)*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((2*im)*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))), (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) - (ℯ^(-2*im*d + (2*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((2*im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(8*sqrt(c)*sqrt(log(f))) - (ℯ^(2*im*d - (2*im*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((2*im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + b*x + c*x^2)*sin(d + e*x)^3, x, 14, If(13>=8, (((-3*im)/16)*ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) + ((im/16)*ℯ^((-3*im)*d + (3*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) - (((3*im)/16)*ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))) + ((im/16)*ℯ^((3*im)*d - ((3*im)*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(sqrt(c)*sqrt(log(f))), (3*im*ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) - (im*ℯ^(-3*im*d + (3*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((3*im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) - (3*im*ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (im*ℯ^(3*im*d - (3*im*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((3*im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))))]

    @test_int [f^(a + b*x + c*x^2)*sin(d + f*x^2), x, 8, -((im*ℯ^((-im)*d + (b^2*log(f)^2)/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(4*sqrt(im*f - c*log(f)))) - (im*ℯ^(im*d - (b^2*log(f)^2)/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(4*sqrt(im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*sin(d + f*x^2)^2, x, 10, (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(-2*im*d + (b^2*log(f)^2)/(8*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(2*im*f - c*log(f)))/(2*sqrt(2*im*f - c*log(f)))))/(8*sqrt(2*im*f - c*log(f))) - (ℯ^(2*im*d - (b^2*log(f)^2)/(8*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(2*im*f + c*log(f)))/(2*sqrt(2*im*f + c*log(f)))))/(8*sqrt(2*im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*sin(d + f*x^2)^3, x, 14, -((3*im*ℯ^((-im)*d + (b^2*log(f)^2)/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(16*sqrt(im*f - c*log(f)))) + (im*ℯ^(-3*im*d + (b^2*log(f)^2)/(12*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(3*im*f - c*log(f)))/(2*sqrt(3*im*f - c*log(f)))))/(16*sqrt(3*im*f - c*log(f))) - (3*im*ℯ^(im*d - (b^2*log(f)^2)/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(16*sqrt(im*f + c*log(f))) + (im*ℯ^(3*im*d - (b^2*log(f)^2)/(4*(3*im*f + c*log(f))))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(3*im*f + c*log(f)))/(2*sqrt(3*im*f + c*log(f)))))/(16*sqrt(3*im*f + c*log(f)))]

    @test_int [f^(a + b*x + c*x^2)*sin(d + e*x + f*x^2), x, 8, (im*ℯ^((-im)*d - (e + im*b*log(f))^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e - b*log(f) + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(4*sqrt(im*f - c*log(f))) - (im*ℯ^(im*d + (e - im*b*log(f))^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(4*sqrt(im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*sin(d + e*x + f*x^2)^2, x, 10, (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) - (ℯ^(-2*im*d - (2*e + im*b*log(f))^2/(8*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((2*im*e - b*log(f) + 2*x*(2*im*f - c*log(f)))/(2*sqrt(2*im*f - c*log(f)))))/(8*sqrt(2*im*f - c*log(f))) - (ℯ^(2*im*d + (2*e - im*b*log(f))^2/(8*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((2*im*e + b*log(f) + 2*x*(2*im*f + c*log(f)))/(2*sqrt(2*im*f + c*log(f)))))/(8*sqrt(2*im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*sin(d + e*x + f*x^2)^3, x, 14, (3*im*ℯ^((-im)*d - (e + im*b*log(f))^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e - b*log(f) + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(16*sqrt(im*f - c*log(f))) - (im*ℯ^(-3*im*d - (3*e + im*b*log(f))^2/(4*(3*im*f - c*log(f))))*f^a*sqrt(pi)*erf((3*im*e - b*log(f) + 2*x*(3*im*f - c*log(f)))/(2*sqrt(3*im*f - c*log(f)))))/(16*sqrt(3*im*f - c*log(f))) - (3*im*ℯ^(im*d + (e - im*b*log(f))^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(16*sqrt(im*f + c*log(f))) + (im*ℯ^(3*im*d - (3*im*e + b*log(f))^2/(4*(3*im*f + c*log(f))))*f^a*sqrt(pi)*erfi((3*im*e + b*log(f) + 2*x*(3*im*f + c*log(f)))/(2*sqrt(3*im*f + c*log(f)))))/(16*sqrt(3*im*f + c*log(f)))]

    @test_int [f^(a + b*x + c*x^2)*sin(a + b*x + e*x^2), x, 8, If(13>=8, (im*sqrt(pi)*erf((b*(im - log(f)) + 2*x*(im*e - c*log(f)))/(2*sqrt(im*e - c*log(f)))))/(ℯ^((im - log(f))*(a - (b^2*(im - log(f)))/(4*im*e - 4*c*log(f))))*(4*sqrt(im*e - c*log(f)))) - (im*ℯ^((im + log(f))*(a - (b^2*(im + log(f)))/(4*im*e + 4*c*log(f))))*sqrt(pi)*erfi((b*(im + log(f)) + 2*x*(im*e + c*log(f)))/(2*sqrt(im*e + c*log(f)))))/(4*sqrt(im*e + c*log(f))), -((im*sqrt(pi)*erf(-((b*(im - log(f)) + 2*x*(im*e - c*log(f)))/(2*sqrt(im*e - c*log(f))))))/(ℯ^((im - log(f))*(a - (b^2*(im - log(f)))/(4*im*e - 4*c*log(f))))*(4*sqrt(im*e - c*log(f))))) - (im*ℯ^((im + log(f))*(a - (b^2*(im + log(f)))/(4*im*e + 4*c*log(f))))*sqrt(pi)*erfi((b*(im + log(f)) + 2*x*(im*e + c*log(f)))/(2*sqrt(im*e + c*log(f)))))/(4*sqrt(im*e + c*log(f))))]

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

    @test_int [ℯ^x*cos(a + b*x), x, 1, (ℯ^x*cos(a + b*x))/(1 + b^2) + (b*ℯ^x*sin(a + b*x))/(1 + b^2)]
    @test_int [ℯ^x*cos(a + c*x^2), x, 6, -((-1)^(1/4)*ℯ^((im/4)*(4*a + c^(-1)))*sqrt(pi)*erf(((-1)^(1/4)*(1 + (2*im)*c*x))/(2*sqrt(c))))/(4*sqrt(c)) + ((-1)^(1/4)*sqrt(pi)*erfi(((-1)^(1/4)*(1 - (2*im)*c*x))/(2*sqrt(c))))/(4*sqrt(c)*ℯ^((im/4)*(4*a + c^(-1))))]
    @test_int [ℯ^x*cos(a + b*x + c*x^2), x, 6, -(((-1)^(1/4)*ℯ^((1/4)*im*(4*a + (1 + im*b)^2/c))*sqrt(pi)*erf(((-1)^(1/4)*(1 + im*b + 2*im*c*x))/(2*sqrt(c))))/(4*sqrt(c))) + ((-1)^(1/4)*ℯ^((-im)*a + (im*(im + b)^2)/(4*c))*sqrt(pi)*erfi(((-1)^(1/4)*(1 - im*b - 2*im*c*x))/(2*sqrt(c))))/(4*sqrt(c))]

    @test_int [ℯ^x^2*cos(a + b*x), x, 6, (ℯ^((-im)*a + b^2/4)*sqrt(pi)*erfi(((-im)*b + 2*x)/2))/4 + (ℯ^(im*a + b^2/4)*sqrt(pi)*erfi((im*b + 2*x)/2))/4]
    @test_int [ℯ^x^2*cos(a + c*x^2), x, 4, (sqrt(pi)*erfi(sqrt(1 - im*c)*x))/(4*sqrt(1 - im*c)*ℯ^(im*a)) + (ℯ^(im*a)*sqrt(pi)*erfi(sqrt(1 + im*c)*x))/(4*sqrt(1 + im*c))]
    @test_int [ℯ^x^2*cos(a + b*x + c*x^2), x, 6, -((sqrt(pi)*erfi((im*b - 2*(1 - im*c)*x)/(2*sqrt(1 - im*c))))/(ℯ^(im*(a - b^2/(4*im + 4*c)))*(4*sqrt(1 - im*c)))) + (ℯ^(im*a + b^2/(4*(1 + im*c)))*sqrt(pi)*erfi((im*b + 2*(1 + im*c)*x)/(2*sqrt(1 + im*c))))/(4*sqrt(1 + im*c))]

    @test_int [f^(a + b*x)*cos(d + f*x^2), x, 8, -((-1)^(1/4)*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*((2*im)*f*x + b*log(f)))/(2*sqrt(f))))/4 - ((-1)^(1/4)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*((2*im)*f*x - b*log(f)))/(2*sqrt(f))))/(4*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f)))]
    @test_int [f^(a + b*x)*cos(d + f*x^2)^2, x, 9, (-1/16 - im/16)*ℯ^((2*im)*d + ((im/8)*b^2*log(f)^2)/f)*f^(-1/2 + a)*sqrt(pi)*erf(((1/4 + im/4)*((4*im)*f*x + b*log(f)))/sqrt(f)) - ((1/16 + im/16)*f^(-1/2 + a)*sqrt(pi)*erfi(((1/4 + im/4)*((4*im)*f*x - b*log(f)))/sqrt(f)))/ℯ^((im/8)*(16*d + (b^2*log(f)^2)/f)) + f^(a + b*x)/(2*b*log(f))]
    @test_int [f^(a + b*x)*cos(d + f*x^2)^3, x, 14, (-3*(-1)^(1/4)*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*((2*im)*f*x + b*log(f)))/(2*sqrt(f))))/16 - (1/16 + im/16)*ℯ^((3*im)*d + ((im/12)*b^2*log(f)^2)/f)*f^(-1/2 + a)*sqrt(pi/6)*erf(((1/2 + im/2)*((6*im)*f*x + b*log(f)))/(sqrt(6)*sqrt(f))) - (3*(-1)^(1/4)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*((2*im)*f*x - b*log(f)))/(2*sqrt(f))))/(16*ℯ^((im/4)*(4*d + (b^2*log(f)^2)/f))) - ((1/16 + im/16)*f^(-1/2 + a)*sqrt(pi/6)*erfi(((1/2 + im/2)*((6*im)*f*x - b*log(f)))/(sqrt(6)*sqrt(f))))/ℯ^((im/12)*(36*d + (b^2*log(f)^2)/f))]

    @test_int [f^(a + b*x)*cos(d + e*x + f*x^2), x, 8, -((-1)^(1/4)*ℯ^((im/4)*(4*d + (im*e + b*log(f))^2/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*(im*e + (2*im)*f*x + b*log(f)))/(2*sqrt(f))))/4 - ((-1)^(1/4)*ℯ^((-im)*d + ((im/4)*(e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*(im*e + (2*im)*f*x - b*log(f)))/(2*sqrt(f))))/4]
    @test_int [f^(a + b*x)*cos(d + e*x + f*x^2)^2, x, 9, (-1/16 - im/16)*ℯ^((2*im)*d + ((im/8)*((2*im)*e + b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erf(((1/4 + im/4)*((2*im)*e + (4*im)*f*x + b*log(f)))/sqrt(f)) - (1/16 + im/16)*ℯ^((-2*im)*d + ((im/8)*(2*e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erfi(((1/4 + im/4)*((2*im)*e + (4*im)*f*x - b*log(f)))/sqrt(f)) + f^(a + b*x)/(2*b*log(f))]
    @test_int [f^(a + b*x)*cos(d + e*x + f*x^2)^3, x, 14, (-3*(-1)^(1/4)*ℯ^((im/4)*(4*d + (im*e + b*log(f))^2/f))*f^(-1/2 + a)*sqrt(pi)*erf(((-1)^(1/4)*(im*e + (2*im)*f*x + b*log(f)))/(2*sqrt(f))))/16 - (1/16 + im/16)*ℯ^((3*im)*d + ((im/12)*((3*im)*e + b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi/6)*erf(((1/2 + im/2)*((3*im)*e + (6*im)*f*x + b*log(f)))/(sqrt(6)*sqrt(f))) - (3*(-1)^(1/4)*ℯ^((-im)*d + ((im/4)*(e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi)*erfi(((-1)^(1/4)*(im*e + (2*im)*f*x - b*log(f)))/(2*sqrt(f))))/16 - (1/16 + im/16)*ℯ^((-3*im)*d + ((im/12)*(3*e + im*b*log(f))^2)/f)*f^(-1/2 + a)*sqrt(pi/6)*erfi(((1/2 + im/2)*((3*im)*e + (6*im)*f*x - b*log(f)))/(sqrt(6)*sqrt(f)))]

    @test_int [f^(a + c*x^2)*cos(d + e*x), x, 8, If(13>=8, -(ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))), (ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + c*x^2)*cos(d + e*x)^2, x, 9, If(13>=8, (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) - (ℯ^((-2*im)*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi((im*e - c*x*log(f))/(sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))) + (ℯ^((2*im)*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi((im*e + c*x*log(f))/(sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))), (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(-2*im*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - c*x*log(f))/(sqrt(c)*sqrt(log(f))))))/(8*sqrt(c)*sqrt(log(f))) + (ℯ^(2*im*d + e^2/(c*log(f)))*f^a*sqrt(pi)*erfi((im*e + c*x*log(f))/(sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + c*x^2)*cos(d + e*x)^3, x, 14, If(13>=8, (-3*ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) - (ℯ^((-3*im)*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (3*ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (ℯ^((3*im)*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))), (3*ℯ^((-im)*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) + (ℯ^(-3*im*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((3*im*e - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) + (3*ℯ^(im*d + e^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (ℯ^(3*im*d + (9*e^2)/(4*c*log(f)))*f^a*sqrt(pi)*erfi((3*im*e + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))))]

    @test_int [f^(a + c*x^2)*cos(d + f*x^2), x, 6, (f^a*sqrt(pi)*erf(x*sqrt(im*f - c*log(f))))/(4*ℯ^(im*d)*sqrt(im*f - c*log(f))) + (ℯ^(im*d)*f^a*sqrt(pi)*erfi(x*sqrt(im*f + c*log(f))))/(4*sqrt(im*f + c*log(f)))]
    @test_int [f^(a + c*x^2)*cos(d + f*x^2)^2, x, 7, (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) + (f^a*sqrt(pi)*erf(x*sqrt((2*im)*f - c*log(f))))/(8*ℯ^((2*im)*d)*sqrt((2*im)*f - c*log(f))) + (ℯ^((2*im)*d)*f^a*sqrt(pi)*erfi(x*sqrt((2*im)*f + c*log(f))))/(8*sqrt((2*im)*f + c*log(f)))]
    @test_int [f^(a + c*x^2)*cos(d + f*x^2)^3, x, 10, (3*f^a*sqrt(pi)*erf(x*sqrt(im*f - c*log(f))))/(16*ℯ^(im*d)*sqrt(im*f - c*log(f))) + (f^a*sqrt(pi)*erf(x*sqrt((3*im)*f - c*log(f))))/(16*ℯ^((3*im)*d)*sqrt((3*im)*f - c*log(f))) + (3*ℯ^(im*d)*f^a*sqrt(pi)*erfi(x*sqrt(im*f + c*log(f))))/(16*sqrt(im*f + c*log(f))) + (ℯ^((3*im)*d)*f^a*sqrt(pi)*erfi(x*sqrt((3*im)*f + c*log(f))))/(16*sqrt((3*im)*f + c*log(f)))]

    @test_int [f^(a + c*x^2)*cos(d + e*x + f*x^2), x, 8, (ℯ^((-im)*d - e^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(4*sqrt(im*f - c*log(f))) + (ℯ^(im*d + e^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(4*sqrt(im*f + c*log(f)))]
    @test_int [f^(a + c*x^2)*cos(d + e*x + f*x^2)^2, x, 9, (f^a*sqrt(pi)*erfi(sqrt(c)*x*sqrt(log(f))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(-2*im*d - e^2/(2*im*f - c*log(f)))*f^a*sqrt(pi)*erf((im*e + x*(2*im*f - c*log(f)))/sqrt(2*im*f - c*log(f))))/(8*sqrt(2*im*f - c*log(f))) + (ℯ^(2*im*d + e^2/(2*im*f + c*log(f)))*f^a*sqrt(pi)*erfi((im*e + x*(2*im*f + c*log(f)))/sqrt(2*im*f + c*log(f))))/(8*sqrt(2*im*f + c*log(f)))]
    @test_int [f^(a + c*x^2)*cos(d + e*x + f*x^2)^3, x, 14, (3*ℯ^((-im)*d - e^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(16*sqrt(im*f - c*log(f))) + (ℯ^(-3*im*d - (9*e^2)/(4*(3*im*f - c*log(f))))*f^a*sqrt(pi)*erf((3*im*e + 2*x*(3*im*f - c*log(f)))/(2*sqrt(3*im*f - c*log(f)))))/(16*sqrt(3*im*f - c*log(f))) + (3*ℯ^(im*d + e^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(16*sqrt(im*f + c*log(f))) + (ℯ^(3*im*d + (9*e^2)/(4*(3*im*f + c*log(f))))*f^a*sqrt(pi)*erfi((3*im*e + 2*x*(3*im*f + c*log(f)))/(2*sqrt(3*im*f + c*log(f)))))/(16*sqrt(3*im*f + c*log(f)))]

    @test_int [f^(a + b*x + c*x^2)*cos(d + e*x), x, 8, If(13>=8, -(ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))), (ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(4*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + b*x + c*x^2)*cos(d + e*x)^2, x, 10, If(13>=8, (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) - (ℯ^((-2*im)*d + (2*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((2*im)*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))) + (ℯ^((2*im)*d - ((2*im)*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((2*im)*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))), (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(-2*im*d + (2*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((2*im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(8*sqrt(c)*sqrt(log(f))) + (ℯ^(2*im*d - (2*im*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((2*im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(8*sqrt(c)*sqrt(log(f))))]
    @test_int [f^(a + b*x + c*x^2)*cos(d + e*x)^3, x, 14, If(13>=8, (-3*ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) - (ℯ^((-3*im)*d + (3*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (3*ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (ℯ^((3*im)*d - ((3*im)*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(((3*im)*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))), (3*ℯ^((-im)*d + (e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) + (ℯ^(-3*im*d + (3*e + im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi(-((3*im*e - b*log(f) - 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f))))))/(16*sqrt(c)*sqrt(log(f))) + (3*ℯ^(im*d + (e - im*b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))) + (ℯ^(3*im*d - (3*im*e + b*log(f))^2/(4*c*log(f)))*f^a*sqrt(pi)*erfi((3*im*e + b*log(f) + 2*c*x*log(f))/(2*sqrt(c)*sqrt(log(f)))))/(16*sqrt(c)*sqrt(log(f))))]

    @test_int [f^(a + b*x + c*x^2)*cos(d + f*x^2), x, 8, -((ℯ^((-im)*d + (b^2*log(f)^2)/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(4*sqrt(im*f - c*log(f)))) + (ℯ^(im*d - (b^2*log(f)^2)/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(4*sqrt(im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*cos(d + f*x^2)^2, x, 10, (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) - (ℯ^(-2*im*d + (b^2*log(f)^2)/(8*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(2*im*f - c*log(f)))/(2*sqrt(2*im*f - c*log(f)))))/(8*sqrt(2*im*f - c*log(f))) + (ℯ^(2*im*d - (b^2*log(f)^2)/(8*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(2*im*f + c*log(f)))/(2*sqrt(2*im*f + c*log(f)))))/(8*sqrt(2*im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*cos(d + f*x^2)^3, x, 14, -((3*ℯ^((-im)*d + (b^2*log(f)^2)/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(16*sqrt(im*f - c*log(f)))) - (ℯ^(-3*im*d + (b^2*log(f)^2)/(12*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((b*log(f) - 2*x*(3*im*f - c*log(f)))/(2*sqrt(3*im*f - c*log(f)))))/(16*sqrt(3*im*f - c*log(f))) + (3*ℯ^(im*d - (b^2*log(f)^2)/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(16*sqrt(im*f + c*log(f))) + (ℯ^(3*im*d - (b^2*log(f)^2)/(4*(3*im*f + c*log(f))))*f^a*sqrt(pi)*erfi((b*log(f) + 2*x*(3*im*f + c*log(f)))/(2*sqrt(3*im*f + c*log(f)))))/(16*sqrt(3*im*f + c*log(f)))]

    @test_int [f^(a + b*x + c*x^2)*cos(d + e*x + f*x^2), x, 8, (ℯ^((-im)*d - (e + im*b*log(f))^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e - b*log(f) + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(4*sqrt(im*f - c*log(f))) + (ℯ^(im*d + (e - im*b*log(f))^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(4*sqrt(im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*cos(d + e*x + f*x^2)^2, x, 10, (f^(a - b^2/(4*c))*sqrt(pi)*erfi(((b + 2*c*x)*sqrt(log(f)))/(2*sqrt(c))))/(4*sqrt(c)*sqrt(log(f))) + (ℯ^(-2*im*d - (2*e + im*b*log(f))^2/(8*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((2*im*e - b*log(f) + 2*x*(2*im*f - c*log(f)))/(2*sqrt(2*im*f - c*log(f)))))/(8*sqrt(2*im*f - c*log(f))) + (ℯ^(2*im*d + (2*e - im*b*log(f))^2/(8*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((2*im*e + b*log(f) + 2*x*(2*im*f + c*log(f)))/(2*sqrt(2*im*f + c*log(f)))))/(8*sqrt(2*im*f + c*log(f)))]
    @test_int [f^(a + b*x + c*x^2)*cos(d + e*x + f*x^2)^3, x, 14, (3*ℯ^((-im)*d - (e + im*b*log(f))^2/(4*im*f - 4*c*log(f)))*f^a*sqrt(pi)*erf((im*e - b*log(f) + 2*x*(im*f - c*log(f)))/(2*sqrt(im*f - c*log(f)))))/(16*sqrt(im*f - c*log(f))) + (ℯ^(-3*im*d - (3*e + im*b*log(f))^2/(4*(3*im*f - c*log(f))))*f^a*sqrt(pi)*erf((3*im*e - b*log(f) + 2*x*(3*im*f - c*log(f)))/(2*sqrt(3*im*f - c*log(f)))))/(16*sqrt(3*im*f - c*log(f))) + (3*ℯ^(im*d + (e - im*b*log(f))^2/(4*im*f + 4*c*log(f)))*f^a*sqrt(pi)*erfi((im*e + b*log(f) + 2*x*(im*f + c*log(f)))/(2*sqrt(im*f + c*log(f)))))/(16*sqrt(im*f + c*log(f))) + (ℯ^(3*im*d - (3*im*e + b*log(f))^2/(4*(3*im*f + c*log(f))))*f^a*sqrt(pi)*erfi((3*im*e + b*log(f) + 2*x*(3*im*f + c*log(f)))/(2*sqrt(3*im*f + c*log(f)))))/(16*sqrt(3*im*f + c*log(f)))]

    @test_int [f^(a + b*x + c*x^2)*cos(a + b*x + e*x^2), x, 8, If(13>=8, (sqrt(pi)*erf((b*(im - log(f)) + 2*x*(im*e - c*log(f)))/(2*sqrt(im*e - c*log(f)))))/(ℯ^((im - log(f))*(a - (b^2*(im - log(f)))/(4*im*e - 4*c*log(f))))*(4*sqrt(im*e - c*log(f)))) + (ℯ^((im + log(f))*(a - (b^2*(im + log(f)))/(4*im*e + 4*c*log(f))))*sqrt(pi)*erfi((b*(im + log(f)) + 2*x*(im*e + c*log(f)))/(2*sqrt(im*e + c*log(f)))))/(4*sqrt(im*e + c*log(f))), -((sqrt(pi)*erf(-((b*(im - log(f)) + 2*x*(im*e - c*log(f)))/(2*sqrt(im*e - c*log(f))))))/(ℯ^((im - log(f))*(a - (b^2*(im - log(f)))/(4*im*e - 4*c*log(f))))*(4*sqrt(im*e - c*log(f))))) + (ℯ^((im + log(f))*(a - (b^2*(im + log(f)))/(4*im*e + 4*c*log(f))))*sqrt(pi)*erfi((b*(im + log(f)) + 2*x*(im*e + c*log(f)))/(2*sqrt(im*e + c*log(f)))))/(4*sqrt(im*e + c*log(f))))]

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

    @test_int [F^(c*(a + b*x))*(f + f*sin(d + e*x))^2, x, 8, (f^2*F^(a*c + b*c*x))/(b*c*log(F)) - (2*e*f^2*F^(a*c + b*c*x)*cos(d + e*x))/(e^2 + b^2*c^2*log(F)^2) + (2*e^2*f^2*F^(a*c + b*c*x))/(b*c*log(F)*(4*e^2 + b^2*c^2*log(F)^2)) + (2*b*c*f^2*F^(a*c + b*c*x)*log(F)*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2) - (2*e*f^2*F^(a*c + b*c*x)*cos(d + e*x)*sin(d + e*x))/(4*e^2 + b^2*c^2*log(F)^2) + (b*c*f^2*F^(a*c + b*c*x)*log(F)*sin(d + e*x)^2)/(4*e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*(f + f*sin(d + e*x))^1, x, 6, (f*F^(a*c + b*c*x))/(b*c*log(F)) - (e*f*F^(a*c + b*c*x)*cos(d + e*x))/(e^2 + b^2*c^2*log(F)^2) + (b*c*f*F^(a*c + b*c*x)*log(F)*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))/(f + f*sin(d + e*x))^1, x, 2, -((2*ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/e, 2 - (im*b*c*log(F))/e, im*ℯ^(im*(d + e*x))))/(f*(e - im*b*c*log(F))))]
    @test_int [F^(c*(a + b*x))/(f + f*sin(d + e*x))^2, x, 3, -((F^(c*(a + b*x))*cot(d/2 + pi/4 + (e*x)/2)*csc(d/2 + pi/4 + (e*x)/2)^2)/(6*e*f^2)) - (b*c*F^(c*(a + b*x))*csc(d/2 + pi/4 + (e*x)/2)^2*log(F))/(6*e^2*f^2) - (2*ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/e, 2 - (im*b*c*log(F))/e, im*ℯ^(im*(d + e*x)))*(e + im*b*c*log(F)))/(3*e^2*f^2)]

    @test_int [F^(c*(a + b*x))*(f + f*cos(d + e*x))^2, x, 8, (f^2*F^(a*c + b*c*x))/(b*c*log(F)) + (2*b*c*f^2*F^(a*c + b*c*x)*cos(d + e*x)*log(F))/(e^2 + b^2*c^2*log(F)^2) + (2*e^2*f^2*F^(a*c + b*c*x))/(b*c*log(F)*(4*e^2 + b^2*c^2*log(F)^2)) + (b*c*f^2*F^(a*c + b*c*x)*cos(d + e*x)^2*log(F))/(4*e^2 + b^2*c^2*log(F)^2) + (2*e*f^2*F^(a*c + b*c*x)*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2) + (2*e*f^2*F^(a*c + b*c*x)*cos(d + e*x)*sin(d + e*x))/(4*e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))*(f + f*cos(d + e*x))^1, x, 6, (f*F^(a*c + b*c*x))/(b*c*log(F)) + (b*c*f*F^(a*c + b*c*x)*cos(d + e*x)*log(F))/(e^2 + b^2*c^2*log(F)^2) + (e*f*F^(a*c + b*c*x)*sin(d + e*x))/(e^2 + b^2*c^2*log(F)^2)]
    @test_int [F^(c*(a + b*x))/(f + f*cos(d + e*x))^1, x, 2, (2*ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/e, 2 - (im*b*c*log(F))/e, -ℯ^(im*(d + e*x))))/(f*(im*e + b*c*log(F)))]
    @test_int [F^(c*(a + b*x))/(f + f*cos(d + e*x))^2, x, 3, -((2*ℯ^(im*(d + e*x))*F^(c*(a + b*x))*HypergeometricFunctions._₂F₁(2, 1 - (im*b*c*log(F))/e, 2 - (im*b*c*log(F))/e, -ℯ^(im*(d + e*x)))*(im*e - b*c*log(F)))/(3*e^2*f^2)) - (b*c*F^(c*(a + b*x))*log(F)*sec(d/2 + (e*x)/2)^2)/(6*e^2*f^2) + (F^(c*(a + b*x))*sec(d/2 + (e*x)/2)^2*tan(d/2 + (e*x)/2))/(6*e*f^2)]
    end

max_stars_repo_path: test/4 Trig functions/4.7 Miscellaneous/4.7.6 f^(a+b x+c x^2) trig(d+e x+f x^2)^n.jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Hartree-Fock Self-Consistent Field Theory

# +
"""Tutorial implementing a basic Hartree--Fock SCF program."""

__authors__ = ["D. Menendez", "D. A. Sirianni"]
__credits__ = ["D. Menendez", "D. G. A. Smith"]
__email__   = "danielmail7@gmail.com"

__copyright__ = "(c) 2014-2020, The Psi4Julia Developers"
__license__   = "BSD-3-Clause"
__date__      = "07/28/2020"
# -

# ## I. Theoretical Overview
# In this tutorial, we will seek to introduce the theory and implementation of the quantum chemical method known as Hartree-Fock Self-Consistent Field Theory (HF-SCF) with restricted orbitals and closed-shell systems (RHF).  This theory seeks to solve the pseudo-eigenvalue matrix equation
#
# $$\sum_{\nu} F_{\mu\nu}C_{\nu i} = \epsilon_i\sum_{\nu}S_{\mu\nu}C_{\nu i}$$
# $${\bf FC} = {\bf SC\epsilon},$$
#
# called the Roothan equations, which can be solved self-consistently for the orbital coefficient matrix **C** and the orbital energy eigenvalues $\epsilon_i$.  The Fock matrix, **F**, has elements $F_{\mu\nu}$ given (in the atomic orbital basis) as
#
# $$F_{\mu\nu} = H_{\mu\nu} + 2(\mu\,\nu\left|\,\lambda\,\sigma)D_{\lambda\sigma} - (\mu\,\lambda\,\right|\nu\,\sigma)D_{\lambda\sigma},$$
#
# where $D_{\lambda\sigma}$ is an element of the one-particle density matrix **D**, constructed from the orbital coefficient matrix **C**:
#
# $$D_{\lambda\sigma} = C_{\sigma i}C_{\lambda i}$$
#
# Formally, the orbital coefficient matrix **C** is a $N\times M$ matrix, where $N$ is the number of atomic basis functions, and $M$ is the total number of molecular orbitals.  Physically, this matrix describes the contribution of every atomic basis function (columns) to a particular molecular orbital (e.g., the $i^{\rm th}$ row).  The density matrix **D** is a square matrix describing the electron density contained in each orbital.  In the molecular orbital basis, the density matrix has elements
#
# $$D_{pq} = \left\{
# \begin{array}{ll}
# 2\delta_{pq} & p\; {\rm occupied} \\
# 0 & p\; {\rm virtual} \\
# \end{array}\right .$$
#
# The total RHF energy is given by
#
# $$E^{\rm RHF}_{\rm total} = E^{\rm RHF}_{\rm elec} + E^{\rm BO}_{\rm nuc},$$
#
# where $E^{\rm RHF}_{\rm elec}$ is the final electronic RHF energy, and $E^{\rm BO}_{\rm nuc}$ is the total nuclear repulsion energy within the Born-Oppenheimer approximation.  To compute the electronic energy, we may use the density matrix in the AO basis:
#
# $$E^{\rm RHF}_{\rm elec} = (F_{\mu\nu} + H_{\mu\nu})D_{\mu\nu},$$
#
# and the nuclear repulsion energy is simply
#
# $$E^{\rm BO}_{\rm nuc} = \sum_{A>B}\frac{Z_AZ_B}{r_{AB}}$$
#
# where $Z_A$ is the nuclear charge of atom $A$, and the sum runs over all unique nuclear pairs.

# ## II. Implementation
#
# Using the above overview, let's write a RHF program using  Psi4 , NumPy, TensorOperations, and Julia's LinearAlgebra standard library.  First, we need to import these modules:

# ==> Import Psi4 & NumPy <==
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
using TensorOperations: @tensor
using LinearAlgebra: Diagonal, Hermitian, eigen, tr
using Printf: @printf

# Next, using what you learned in the previous tutorial module, set the following  Psi4  and molecule options.
#
# Memory & Output specifications:
# - Give 500 Mb of memory to Psi4
# - Set Psi4 output file to "output.dat"
# - Set a variable `numpy_memory` to an acceptable amount of available memory for the working computer to use for storing tensors
#
# Molecule definition:
# - Define the "physicist's water molecule" (O-H bond length = 1.1 Angstroms, HOH bond angle = 104 degrees)
# - Molecular symmetry C1
#
# Computation options:
# - basis set cc-pVDZ
# - SCF type PK
# - Energy convergence criterion to 0.00000001
#

# +
# ==> Set Basic Psi4 Options <==
# Memory specification
psi4.set_memory(Int(5e8))
numpy_memory = 2

# Set output file
psi4.core.set_output_file("output.dat", false)

# Define Physicist's water -- don't forget C1 symmetry!
mol = psi4.geometry("""
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")

# Set computation options
psi4.set_options(Dict("basis" => "cc-pvdz",
                      "scf_type" => "pk",
                      "e_convergence" => 1e-8))
# -

# Since we will be writing our own, iterative RHF procedure, we will need to define options that we can use to tweak our convergence behavior.  For example, if something goes wrong and our SCF doesn't converge, we don't want to spiral into an infinite loop.  Instead, we can specify the maximum number of iterations allowed, and store this value in a variable called `maxiter`.  Here are some good default options for our program:
# ~~~python
# MAXITER = 40
# E_conv = 1.0e-6
# ~~~
# These are by no means the only possible values for these options, and it's encouraged to try different values and see for yourself how different choices affect the performance of our program.  For now, let's use the above as our default.

# ==> Set default program options <==
# Maximum SCF iterations
MAXITER = 40
# Energy convergence criterion
E_conv = 1.0e-6;

# Before we can build our Fock matrix, we'll need to compute the following static one- and two-electron quantities:
#
# - Electron repulsion integrals (ERIs) **I** between our AOs
# - Overlap matrix **S**
# - Core Hamiltonian matrix **H**
#
# Fortunately for us, we can do this using the machinery in  Psi4.  In the first module, you learned about `psi4.core.Wavefunction` and `psi4.core.MintsHelper` classes.  In the cell below, use these classes to perform the following:
#
# 1. Create Class Instances
#
#     a. Build a wavefunction for our molecule and basis set
#
#     b. Create an instance of the `MintsHelper` class with the basis set for the wavefunction
#
# 2. Build overlap matrix, **S**
#
#     a. Get the AO overlap matrix from `MintsHelper`, and cast it into a Julia array
#
#     b. Get the number of AO basis functions and number of doubly occupied orbitals from S and the wavefunciton
#
# 3. Compute ERI Tensor, **I**
#
#     a. Get ERI tensor from `MintsHelper`, and cast it into a Julia array
#
# 4. Build core Hamiltonian, **H**
#
#     a. Get AO kinetic energy matrix from `MintsHelper`, and cast it into a Julia array
#
#     b. Get AO potential energy matrix from `MintsHelper`, and cast it into a Julia array
#
#     c. Build core Hamiltonian from kinetic & potential energy matrices

# +
# ==> Compute static 1e- and 2e- quantities with Psi4 <==
# Class instantiation
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option("basis"))
mints = psi4.core.MintsHelper(wfn.basisset())

# Overlap matrix
S = np.asarray(mints.ao_overlap()) # we only need a copy

# Number of basis Functions & doubly occupied orbitals
nbf = size(S)[1]
ndocc = wfn.nalpha()

println("Number of occupied orbitals: ", ndocc)
println("Number of basis functions: ", nbf)

# Memory check for ERI tensor
I_size = nbf^4 * 8.e-9
println("\nSize of the ERI tensor will be $I_size GB.")
memory_footprint = I_size * 1.5
if I_size > numpy_memory
    psi4.core.clean()
    throw(OutOfMemoryError("Estimated memory utilization ($memory_footprint GB) exceeds " *
                           "allotted memory limit of $numpy_memory GB."))
end

# Build ERI Tensor
I = np.asarray(mints.ao_eri()) # we only need a copy

# Build core Hamiltonian
T = np.asarray(mints.ao_kinetic()) # we only need a copy
V = np.asarray(mints.ao_potential()) # we only need a copy
H = T + V;
# -

# The Roothan equations
#
# $${\bf FC} = {\bf SC\epsilon}$$
#
# are only *pseudo*-eigenvalue equations due to the presence of the overlap matrix **S** on the right hand side of the equation.  Normally, the AO basis set will not be orthonormal, so the overlap matrix **S** will not be unity and therefore cannot be ignored.  Let's check to see whether our AO basis is orthonormal:

# ==> Inspecting S for AO orthonormality <==
hope = S ≈ Diagonal(ones(size(S)[1]))
println("\nDo we have any hope that our AO basis is orthonormal? ", hope)

# Just as we'd expected -- looks like we can't ignore the AO overlap matrix.  Therefore, the Fock matrix **F** cannot simply be diagonalized to solve for the orbital coefficient matrix **C**.  There is still hope, however!  We can overcome this issue by transforming the AO basis so that all of our basis functions are orthonormal.  In other words, we seek a matrix **A** such that the transformation
#
# $${\bf A}^{\dagger}{\bf SA} = {\bf 1}$$
#
# One method of doing this is called *symmetric orthogonalization*, which lets ${\bf A} = {\bf S}^{-1/2}$.  Then,
#
# $${\bf A}^{\dagger}{\bf SA} = {\bf S}^{-1/2}{\bf SS}^{-1/2} = {\bf S}^{-1/2}{\bf S}^{1/2} = {\bf S}^0 = {\bf 1},$$
#
# and we see that this choice for **A** does in fact yield an orthonormal AO basis.  In the cell below, construct this transformation matrix using  Psi4's built-in `Matrix` class member function `power()` just like the following:
# ~~~python
# A = mints.ao_overlap()
# A.power(-0.5, 1.e-16)
# A = np.asarray(A)
# ~~~

# +
# ==> Construct AO orthogonalization matrix A <==
A = mints.ao_overlap()
A.power(-0.5, 1.e-16) # ≈ Julia's A^(-0.5) after psi4view()
A = np.asarray(A) # we only need a copy

# Check orthonormality
S_p = A * S * A
new_hope = S ≈ Diagonal(ones(size(S)[1]))

if new_hope
    println("There is a new hope for diagonalization!")
else
    println("Whoops...something went wrong. Check that you've correctly built the transformation matrix.")
    @show sum(S_p) - tr(S_p)
end
# -

# The drawback of this scheme is that we would now have to either re-compute the ERI and core Hamiltonian tensors in the newly orthogonal AO basis, or transform them using our **A** matrix (both would be overly costly, especially transforming **I**).  On the other hand, substitute ${\bf C} = {\bf AC}'$ into the Roothan equations:
#
# \begin{align}
# {\bf FAC'} &= {\bf SAC}'{\bf \epsilon}\\
# {\bf A}^{\dagger}({\bf FAC}')&= {\bf A}^{\dagger}({\bf SAC}'){\bf \epsilon}\\
# ({\bf A}^{\dagger}{\bf FA}){\bf C}'&= ({\bf A}^{\dagger}{\bf SA}){\bf C}'{\bf \epsilon}\\
# {\bf F}'{\bf C}' &= {\bf 1C}'{\bf \epsilon}\\
# {\bf F}'{\bf C}' &= {\bf C}'{\bf \epsilon}\\
# \end{align}
#
# Clearly, we have arrived at a canonical eigenvalue equation.  This equation can be solved directly for the transformed orbital coefficient matrix ${\bf C}'$ by diagonalizing the transformed Fock matrix, ${\bf F}'$, before transforming ${\bf C}'$ back into the original AO basis with ${\bf C} = {\bf AC}'$.
#
# Before we can get down to the business of using the Fock matrix **F** to compute the RHF energy, we first need to compute the orbital coefficient **C** matrix.  But, before we compute the **C** matrix, we first need to build **F**.  Wait...hold on a second.  Which comes first, **C** or **F**?  Looking at the Roothan equations more closely, we see that that both sides depend on the **C** matrix, since **F** is a function of the orbitals:
#
#
# $${\bf F}({\bf C}){\bf C} = {\bf SC\epsilon}\,;\;\;F_{\mu\nu} = H_{\mu\nu} + 2(\mu\,\nu\mid\lambda\,\sigma)C_{\sigma i}C_{\lambda i} - (\mu\,\lambda\,\mid\nu\,\sigma)C_{\sigma i}C_{\lambda i}.$$
#
# Therefore technically, *neither* **F** nor **C** can come first!  In order to proceed, we instead begin with a *guess* for the Fock matrix, from which we obtain a guess at the **C** matrix.  Without orbital coefficients (and therefore without electron densities), the most logical starting point for obtaining a guess at the Fock matrix is to begin with the only component of **F** that does *not* involve densities: the core Hamiltonian, **H**.  Below, using the `eigen()` function, and forcing hermitianess with `Hermitian()`, obtain coefficient and density matrices using the core guess:
#
# 1. Obtain ${\bf F}'$ by transforming the core Hamiltonian with the ${\bf A}$ matrix
# 2. Diagonalize the transformed Fock matrix for $\epsilon$ and ${\bf C}'$
# 3. Use doubly-occupied slice of coefficient matrix to build density matrix

# +
# ==> Compute C & D matrices with CORE guess <==
# Transformed Fock matrix
F_p = A * H * A

# Diagonalize F_p for eigenvalues & eigenvectors with NumPy
e, C_p = eigen(Hermitian(F_p))

# Transform C_p back into AO basis
C = A * C_p

# Grab occupied orbitals
C_occ = C[:, 1:ndocc]

# Build density matrix from occupied orbitals
D = C_occ * C_occ' ;
# -

# The final quantity we need to compute before we can proceed with our implementation of the SCF procedure is the Born-Oppenheimer nuclear repulsion energy, $E^{\rm BO}_{\rm nuc}$.  We could use the expression given above in $\S$1, however we can also obtain this value directly from  Psi4's `Molecule` class.  In the cell below, compute the nuclear repulsion energy using either method.

# ==> Nuclear Repulsion Energy <==
E_nuc = mol.nuclear_repulsion_energy()

# Within each SCF iteration, we'll have to perform a number of tensor contractions when building the Fock matrix, computing the total RHF energy, and performing several transformations.  Since the computational expense of this process is related to the number of unique indices, the most intensive step of computing the total electronic energy will be performing the four-index contractions corresponding to building Coulomb and Exchange matrices **J** and **K**, with elements
#
# \begin{align}
# J[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\,\nu\mid\lambda\,\sigma)D_{\lambda\sigma}\\
# K[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\,\lambda\mid\nu\,\sigma)D_{\lambda\sigma},
# \end{align}
#
# when building the Fock matrix.  Fortunately, once **J** and **K** have been built, the Fock matrix may be computed as a simple matrix addition, instead of element-wise:
#
# $$ {\bf F} = {\bf H} + 2{\bf J} - {\bf K} = {\bf H} + {\bf G}.$$
#
# Formation of the **J** and **K** matrices will be the most expensive step of the RHF procedure, scaling with respect to the number of AOs as ${\cal O}(N^4)$.  Strategies for building these marices efficiently, as well as different methods for handling these tensor contractions, will be discussed in greater detail in tutorials 2c and 2d in this module, respectively.
#
# Let's now write our SCF iterations according to the following algorithm:
#
# #### Algorithm 1: SCF Iteration
# for `scf_iter` less than `MAXITER`, do:
# 1. Build Fock matrix
#     - Build the two-electron Coulomb & Exchange matrix **G**
#     - Form the Fock matrix
# 2. RHF Energy
#     - Compute total RHF energy
#     - If change in RHF energy less than `E_conv`, break
#     - Save latest RHF energy as `E_old`
# 3. Compute new orbital guess
#     - Transform Fock matrix to orthonormal AO basis
#     - Diagonalize ${\bf F}'$ for $\epsilon$ and ${\bf C}'$
#     - Back transform ${\bf C}'$ to AO basis
#     - Form **D** from occupied orbital slice of **C**
#

# +
# ==> SCF Iterations <==
# Output and pre-iteration energy declarations
SCF_E = let SCF_E = 0.0, E_old = 0.0, D = D

   print("==> Starting SCF Iterations <==\n")

   # Begin Iterations
   for scf_iter in 1:MAXITER
       # Build Fock matrix
       @tensor G[p,q] := (2I[p,q,r,s] - I[p,r,q,s]) * D[r,s]
       F = H + G

       # Compute RHF energy
       SCF_E = tr((H + F) * D) + E_nuc
       @printf("SCF Iteration %3d: Energy = %4.16f dE = %1.5e \n",scf_iter, SCF_E, SCF_E - E_old)

       # SCF Converged?
       if abs(SCF_E - E_old) < E_conv
           break
       end
       E_old = SCF_E

       # Compute new orbital guess
       F_p =  A * F * A
       e, C_p = eigen(Hermitian(F_p))
       C = A * C_p
       C_occ = C[:, 1:ndocc]
       D = C_occ * C_occ'

       # MAXITER exceeded?
       if scf_iter == MAXITER
           psi4.core.clean()
           throw(MethodError("Maximum number of SCF iterations exceeded."))
       end
   end
   SCF_E # return RHF SCF energy
end

# Post iterations
println("\nSCF converged.")
println("Final RHF Energy: $SCF_E [Eh]")
# -

# Congratulations! You've written your very own Restricted Hartree-Fock program!  Finally, let's check your final RHF energy against  Psi4:

# Compare to Psi4
SCF_E_psi = psi4.energy("SCF")
psi4.compare_values(SCF_E_psi, SCF_E, 6, "SCF Energy")

# ## References
# 1. [[Szabo:1996](http://store.doverpublications.com/0486691861.html)] A. Szabo and N. S. Ostlund, *Modern Quantum Chemistry*, Introduction to Advanced Electronic Structure Theory. Courier Corporation, 1996.
# 2. [[Levine:2000](https://books.google.com/books?id=80RpQgAACAAJ&dq=levine%20quantum%20chemistry%205th%20edition&source=gbs_book_other_versions)] I. N. Levine, *Quantum Chemistry*. Prentice-Hall, New Jersey, 5th edition, 2000.
# 3. [[Helgaker:2000](https://books.google.com/books?id=lNVLBAAAQBAJ&pg=PT1067&dq=helgaker+molecular+electronic+structure+theory&hl=en&sa=X&ved=0ahUKEwj37I7MkofUAhWG5SYKHaoPAAkQ6AEIKDAA#v=onepage&q=helgaker%20molecular%20electronic%20structure%20theory&f=false)] T. Helgaker, P. Jorgensen, and J. Olsen, *Molecular Electronic Structure Theory*, John Wiley & Sons Inc, 2000.

ext: jl


text:
# algorithm

"""
    SinkhornGibbs()

Vanilla Sinkhorn algorithm.
"""
struct SinkhornGibbs <: Sinkhorn end

Base.show(io::IO, ::SinkhornGibbs) = print(io, "Sinkhorn algorithm")

# cache

struct SinkhornGibbsCache{U,V,KT}
    u::U
    v::V
    K::KT
    Kv::U
end

function build_cache(
    ::Type{T},
    ::SinkhornGibbs,
    size2::Tuple,
    μ::AbstractVecOrMat,
    ν::AbstractVecOrMat,
    C::AbstractMatrix,
    ε::Real,
) where {T}
    # compute Gibbs kernel (has to be mutable for ε-scaling algorithm)
    K = similar(C, T)
    @. K = exp(-C / ε)

    # create and initialize dual potentials
    u = similar(μ, T, size(μ, 1), size2...)
    v = similar(ν, T, size(ν, 1), size2...)
    fill!(u, one(T))
    fill!(v, one(T))

    # cache for next iterate of `u`
    Kv = similar(u)

    return SinkhornGibbsCache(u, v, K, Kv)
end

# use `SinkhornGibbs` as default algorithm
# TODO: remove deprecations
function sinkhorn(
    μ,
    ν,
    C,
    ε;
    tol=nothing,
    atol=tol,
    check_marginal_step=nothing,
    check_convergence=check_marginal_step,
    kwargs...,
)
    if tol !== nothing
        Base.depwarn(
            "keyword argument `tol` is deprecated, please use `atol` and `rtol`", :sinkhorn
        )
    end
    if check_marginal_step !== nothing
        Base.depwarn(
            "keyword argument `check_marginal_step` is deprecated, please use `check_convergence`",
            :sinkhorn,
        )
    end

    _check_convergence = check_convergence === nothing ? 10 : check_convergence

    return sinkhorn(
        μ,
        ν,
        C,
        ε,
        SinkhornGibbs();
        atol=atol,
        check_convergence=_check_convergence,
        kwargs...,
    )
end

function sinkhorn2(
    μ,
    ν,
    C,
    ε;
    tol=nothing,
    atol=tol,
    check_marginal_step=nothing,
    check_convergence=check_marginal_step,
    kwargs...,
)
    if tol !== nothing
        Base.depwarn(
            "keyword argument `tol` is deprecated, please use `atol` and `rtol`", :sinkhorn2
        )
    end
    if check_marginal_step !== nothing
        Base.depwarn(
            "keyword argument `check_marginal_step` is deprecated, please use `check_convergence`",
            :sinkhorn2,
        )
    end

    _check_convergence = check_convergence === nothing ? 10 : check_convergence

    return sinkhorn2(
        μ,
        ν,
        C,
        ε,
        SinkhornGibbs();
        atol=atol,
        check_convergence=_check_convergence,
        kwargs...,
    )
end

# interface

prestep!(::SinkhornSolver{SinkhornGibbs}, ::Int) = nothing

function plan(solver::SinkhornSolver{SinkhornGibbs})
    cache = solver.cache
    γ = cache.K .* add_singleton(cache.u, Val(2)) .* add_singleton(cache.v, Val(1))
    return γ
end


text:
using FastGaussQuadrature
using Serialization
#using Compat: @info

const GL_NODES, GL_WEIGHTS = gausslegendre(100)
const LGF_DIR  = joinpath(@__DIR__, "cache")
const LGF_FILE = joinpath(LGF_DIR, "lgftable.dat")

function load_lgf(N)
    if isfile(LGF_FILE)
        G = deserialize(open(LGF_FILE, "r"))
        if size(G,1) ≥ N
            return G
        end
    end
    build_lgf(N)
end

function build_lgf(N)
    @info "Building and caching LGF table"

    g = zeros(N, N)
    for y in 0:N-1, x in 0:y
        g[x+1,y+1] = lgf(x,y)
    end

    G = Symmetric(g)
    mkpath(LGF_DIR)
    #open(LGF_FILE, "w") do f
    #    serialize(f, G)
    #end
    serialize(open(LGF_FILE,"w"),G)
    G
end

quadgauss(f::Function) = dot(GL_WEIGHTS, f(GL_NODES))

function lgf(i, j)
    if i == j ==0
        return 0.0
    elseif i ≥ j
        v = quadgauss() do x
            if x == -1
                return sqrt(2)abs(i)
            else
                t = (x .+ 1)./2
                return 0.5real((1 .-
                                ( (t.-sqrt(1im))./(t.+sqrt(1im)) ).^(j.+abs(i)).*
                                ( (t.+sqrt(-1im))./(t.-sqrt(-1im)) ).^(j.-abs(i)) ))./t
            end

        end
        return 0.5v/pi
    else
        return lgf(j,i)
    end

end

const LGF_TABLE = load_lgf(1500)


text:
"""
```
impulse_responses(m, system)

impulse_responses(system, horizon)
```

Compute impulse responses for a single draw.

### Inputs

- `m::AbstractDSGEModel`: model object
- `system::System{S}`: state-space system matrices
- `horizon::Int`: number of periods ahead to forecast
- `flip_shocks::Bool`: Whether to compute IRFs in response to a positive shock (by default the shock magnitude is a negative 1 std. shock)

where `S<:AbstractFloat`

### Outputs

- `states::Array{S, 3}`: matrix of size `nstates` x `horizon` x `nshocks` of
  state impulse response functions
- `obs::Array{S, 3}`: matrix of size `nobs` x `horizon` x `nshocks` of
  observable impulse response functions
- `pseudo::Array{S, 3}`: matrix of size `npseudo` x `horizon` x `nshocks` of
  pseudo-observable impulse response functions
"""
function impulse_responses(m::AbstractRepModel, system::System{S};
                           flip_shocks::Bool = false) where {S<:AbstractFloat}
    horizon = impulse_response_horizons(m)
    states, obs, pseudo = impulse_responses(system, horizon, flip_shocks = flip_shocks)
    return states, obs, pseudo
end

function impulse_responses_augmented(m::Union{AbstractHetModel,RepDSGEGovDebt}, system::System{S};
                                     flip_shocks::Bool = false,
                                     use_alternate_consumption = false) where {S<:AbstractFloat}
    horizon = impulse_response_horizons(m)

    TTT_jump, TTT_state, eu  = klein(m)
    if eu == -1
        throw(KleinError())
    end

    TTT, RRR = klein_transition_matrices(m, TTT_state, TTT_jump)
    CCC      = zeros(n_model_states(m))

    TTT, RRR, CCC = if typeof(m) != RepDSGEGovDebt{Float64}
        augment_states(m, TTT, TTT_jump, RRR, CCC)
    else
        TTT, RRR, CCC
    end

    measurement_equation = measurement(m, TTT, RRR, CCC)
    transition_equation = Transition(TTT, RRR, CCC)
    system = System(transition_equation, measurement_equation)

    states, obs, pseudo = impulse_responses(system, horizon, flip_shocks = flip_shocks)

    state_indices_orig = stack_indices(m.endogenous_states_original, get_setting(m, :states))
    jump_indices_orig  = stack_indices(m.endogenous_states_original, get_setting(m, :jumps))

    if typeof(m) == RepDSGEGovDebt{Float64}
        return states, obs, pseudo
    end

    Qx = get_setting(m, :Qx)
    Qy = get_setting(m, :Qy)

    # In this case, the length of state_indices and jump_indices seems to give the number of
    # UNNORMALIZED states/jumps however I'm not sure if this will always be the case/if this is
    # really what we want saved into these settings...so might need to modify in future.
    states_unnormalized = Array{Float64}(undef, length(state_indices_orig),
                                         horizon, size(states, 3))
    jumps_unnormalized = Array{Float64}(undef, length(jump_indices_orig), horizon, size(states, 3))

    model_states_unnormalized = Array{Float64}(undef, length(state_indices_orig) +
                                               length(jump_indices_orig) +
                                               length(m.endogenous_states_augmented), horizon,
                                               size(states, 3))
    for i in 1:size(states, 3)
        state_inds = 1:n_backward_looking_states(m)
        states_unnormalized[:, :, i] = Matrix(Qx')*states[state_inds, :, i]
        jump_inds = n_backward_looking_states(m)+1:n_backward_looking_states(m) + n_jumps(m)
        jumps_unnormalized[:, :, i] = Matrix(Qy')*states[jump_inds, :, i]

        augmented_states = states[n_backward_looking_states(m) + n_jumps(m) + 1 : end, :, i]
        model_states_unnormalized[:, :, i] = vcat(states_unnormalized[:, :, i],
                                                  jumps_unnormalized[:, :, i],
                                                  augmented_states)
    end

    if use_alternate_consumption
        endo = m.endogenous_states_original
        c_implied = (m[:ystar]/m[:g]) - m[:xstar]
        IRFC_implied = m[:ystar]/(c_implied*m[:g])*(model_states_unnormalized[endo[:y′_t], :, :] -
                                                    model_states_unnormalized[endo[:g′_t], :, :])-
                                                    (m[:xstar]/c_implied)*
        model_states_unnormalized[endo[:I′_t], :, :]

        model_states_unnormalized[1229, :, :] = IRFC_implied

        z_consumption = model_states_unnormalized[endo[:z′_t], :, :]
        g_lag = cat(0, model_states_unnormalized[endo[:g′_t], :, :], dims = 2)[:, 1:40, :]
        IRFC_implied_lag = m[:ystar]/(c_implied*m[:g])*(model_states_unnormalized[endo[:y′_t1], :, :]  - g_lag) - (m[:xstar]/c_implied)*model_states_unnormalized[endo[:I′_t1], :, :]
        IRF_observable_implied = IRFC_implied - IRFC_implied_lag .+ z_consumption
        obs[m.observables[:obs_consumption], :, :] = IRF_observable_implied
    end

    return model_states_unnormalized, obs, pseudo
end

function impulse_responses(m::Union{AbstractHetModel,RepDSGEGovDebt}, system::System{S};
                           flip_shocks::Bool = false,
                           use_augmented_states::Bool = true,
                           use_alternate_consumption::Bool = false) where {S<:AbstractFloat}
    if use_augmented_states
        return impulse_responses_augmented(m, system, flip_shocks = flip_shocks,
                                           use_alternate_consumption = use_alternate_consumption)
    end

    horizon = impulse_response_horizons(m)
    states, obs, pseudo = impulse_responses(system, horizon, flip_shocks = flip_shocks)

    if typeof(m) == RepDSGEGovDebt
        return states, obs, pseudo
    end

    Qx = get_setting(m, :Qx)
    Qy = get_setting(m, :Qy)

    state_indices_orig = stack_indices(m.endogenous_states_original, get_setting(m, :states))
    jump_indices_orig  = stack_indices(m.endogenous_states_original, get_setting(m, :jumps))

    # In this case, the length of state_indices and jump_indices seems to give the number of
    # UNNORMALIZED states/jumps however I'm not sure if this will always be the case/if this is
    # really what we want saved into these settings...so might need to modify in future.
    states_unnormalized = Array{Float64}(undef, length(state_indices_orig),
                                         horizon, size(states, 3))
    jumps_unnormalized = Array{Float64}(undef, length(jump_indices_orig), horizon, size(states, 3))
    model_states_unnormalized = Array{Float64}(undef, length(state_indices_orig) +
                                               length(jump_indices_orig), horizon, size(states, 3))
    for i in 1:size(states, 3)
        state_inds = 1:n_backward_looking_states(m)
        states_unnormalized[:, :, i] = Matrix(Qx')*states[state_inds, :, i]
        jump_inds = n_backward_looking_states(m)+1:n_backward_looking_states(m) + n_jumps(m)
        jumps_unnormalized[:, :, i] = Matrix(Qy')*states[jump_inds, :, i]

        model_states_unnormalized[:, :, i] = vcat(states_unnormalized[:, :, i],
                                                  jumps_unnormalized[:, :, i])
    end

    if use_alternate_consumption
        endo_orig = m.endogenous_states_original
        c_implied = (m[:ystar]/m[:g]) - m[:xstar]
        IRFC_implied = m[:ystar]/(c_implied*m[:g])*(model_states_unnormalized[endo_orig[:y′_t], :, :]  - model_states_unnormalized[endo_orig[:g′_t], :, :]) - (m[:xstar]/c_implied)*model_states_unnormalized[endo_orig[:I′_t], :, :]
        model_states_unnormalized = vcat(model_states_unnormalized, IRFC_implied)

        z_consumption = model_states_unnormalized[m.endogenous_states_original[:z′_t], :, :]
        IRFC_implied_lag = cat(0, IRFC_implied, dims = 2)[:, 1:40, :]
        IRF_observable_implied = IRFC_implied - IRFC_implied_lag .+ z_consumption
        obs[m.observables[:obs_consumption], :, :] = IRF_observable_implied
    end

    return model_states_unnormalized, obs, pseudo
end

function impulse_responses(system::System{S}, horizon::Int;
                           flip_shocks::Bool = false) where {S<:AbstractFloat}
    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    for i = 1:nshocks
        # Isolate single shock
        shocks = zeros(S, nshocks, horizon)
        if flip_shocks
            shocks[i, 1] = sqrt(system[:QQ][i, i]) # a positive 1 s.d. shock
        else
            shocks[i, 1] = -sqrt(system[:QQ][i, i]) # a negative 1 s.d. shock
        end
        # Iterate state space forward
        states[:, :, i], obs[:, :, i], pseudo[:, :, i], _ = forecast(system, s_0, shocks)
    end

    return states, obs, pseudo
end

# Method for specifying the subset of shocks, and the size of each shock
function impulse_responses(m::AbstractDSGEModel, system::System{S},
                           horizon::Int, shock_names::Vector{Symbol},
                           shock_values::Vector{Float64}) where S<:AbstractFloat

    # Must provide a name and value for each shock
    @assert length(shock_names) == length(shock_values)

    # Setup
    exo          = m.exogenous_shocks
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    for (i, shock) in enumerate(shock_names)
        # Isolate single shock
        shocks = zeros(S, nshocks, horizon)
        shocks[exo[shock], 1] = shock_values[i]

        # Iterate state space forward
        states[:, :, exo[shock]], obs[:, :, exo[shock]], pseudo[:, :, exo[shock]], _ = forecast(system, s_0, shocks)
    end

    return states, obs, pseudo
end

# Method for picking a specific shock and the size of the desired shift in the state or
# observed variable using that shock. (Omitting the feature to do pseudo-observables
# now, for simplicity)
function impulse_responses(m::AbstractDSGEModel, system::System{S},
                           horizon::Int, shock_name::Symbol,
                           var_name::Symbol, var_value::Float64) where S<:AbstractFloat
    # Setup
    var_names, var_class =
    if var_name in keys(m.endogenous_states)
        m.endogenous_states, :states
    elseif var_name in keys(m.observables)
        m.observables, :obs
    end
    exo          = m.exogenous_shocks
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    # Isolate single shock
    shocks = zeros(S, nshocks, horizon)
    if var_class == :states
        shocks[exo[shock_name], 1] = obtain_shock_from_desired_state_value(var_value,
                                                                      var_names[var_name],
                                                                      exo[shock_name],
                                                                      system[:RRR])
    else # == :obs
        shocks[exo[shock_name], 1] = obtain_shock_from_desired_obs_value(var_value,
                                                                    var_names[var_name],
                                                                    exo[shock_name],
                                                                    system[:ZZ],
                                                                    system[:RRR])
    end

    # Iterate state space forward
    states[:, :, exo[shock_name]], obs[:, :, exo[shock_name]], pseudo[:, :, exo[shock_name]], _ = forecast(system, s_0, shocks)

    return states, obs, pseudo
end

function impulse_responses_peg(m::AbstractRepModel, system::System{S}, horizon::Int;
                           flip_shocks::Bool = false, H::Int = 0, peg::Symbol = :all_periods, real_rate = false) where {S<:AbstractFloat}
    ## H is peg horizon

    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    shocks = m.exogenous_shocks

    states = zeros(S, nstates, horizon) #, nshocks)
    obs    = zeros(S, nobs,    horizon) #, nshocks)
    pseudo = zeros(S, npseudo, horizon) #, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    ### IRFs to FFR peg for H periods
    # Set anticipated shocks to peg FFR from t+1 to t+H
    FFRpeg = -0.25/4 #/400

    PsiR1 = 0 #constant
    # ZZ matrix
    PsiR2 = zeros(nstates)
    PsiR2[m.endogenous_states[:R_t]] = 1
    if real_rate
        PsiR2[m.endogenous_states[:Eπ_t]] = -1
    end
    if H == 0
        Rht = system[:RRR][:,shocks[:rm_sh]]
    else
        Rht = system[:RRR][:,vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")])] #% Columns of Rh referring to anticipated monetary shocks
    end
    bb = zeros(H+1,1)
    MH = zeros(H+1,H+1)
    for hh = 1:H+1
        if peg == :all_periods
           # @show size(FFRpeg), size(PsiR1), size(PsiR2), size(system[:TTT]), hh, size(s_0)
            #@show size(FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_0)
            bb[hh,1] = (FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_0)
            #@show bb[hh, 1]
            MH[hh,:] = PsiR2'*(system[:TTT])^(hh-1)*Rht
        elseif peg == :some_periods
            if hh < H+1
                FFRpeg_ = 0
            else
                FFRpeg_ = FFRpeg
            end
            bb[hh,1] = (FFRpeg_ - PsiR1 - PsiR2'*(system[:TTT])^hh*s_0)[1]
            MH[hh,:] = PsiR2'*(system[:TTT])^(hh-1)*Rht
        end

    end
    monshocks = MH\bb
    #@show size(MH), bb, monshocks
    #% Simulate the model
    etpeg = zeros(nshocks,horizon)
    if H == 0
        etpeg[shocks[:rm_sh],1] = monshocks[1]
    else
        etpeg[vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")]),1] = monshocks
    end
    states[:, :], obs[:, :], pseudo[:, :] = forecast(system, s_0, etpeg)
    return states, obs, pseudo
end

# Given a desired value, x, for some state s^i_1, back out the necessary value of
# ϵ^j_1 needed s.t. s^i_1 = x.
# E.g. if I want the impulse response of output (s^i_1) to be 50 basis points
# based on a shock to productivity (ϵ^j_1), how large does that shock
# to productivity need to be?
# The answer is ϵ^j_1 = x/R_{i,j}
function obtain_shock_from_desired_state_value(state_value::Float64, state_ind::Int,
                                               shock_ind::Int, RRR::Matrix{Float64})
    return state_value/RRR[state_ind, shock_ind]
end

# Similar principle to the function for states
# The answer is ϵ^j_1 = x/{Σ_{k=1}^n Z_{i, k} R_{k, j}},
# where n is the total number of states, i is the index of the observable of interest
# y^i_1, and j is the index of the shock of interest, ϵ^j_1
function obtain_shock_from_desired_obs_value(obs_value::Float64, obs_ind::Int,
                                             shock_ind::Int, ZZ::Matrix{Float64},
                                             RRR::Matrix{Float64})
    return obs_value/dot(ZZ[obs_ind, :], RRR[:, shock_ind])
end


text:
import Base:
  # types
  Array,
  CartesianIndices,
  Vector,
  # symbols
  +,
  -,
  *,
  /,
  ==,
  <,
  >,
  # functions
  adjoint,
  axes,
  complex,
  convert,
  copy,
  copyto!,
  deepcopy,
  deleteat!,
  eachindex,
  eltype,
  fill!,
  getindex,
  hash,
  isapprox,
  isempty,
  isless,
  iterate,
  keys,
  lastindex,
  length,
  map,
  map!,
  ndims,
  push!,
  resize!,
  setindex!,
  show,
  similar,
  size,
  summary,
  truncate,
  zero

import Base.Broadcast:
  # types
  AbstractArrayStyle,
  Broadcasted,
  BroadcastStyle,
  DefaultArrayStyle,
  Style,
  # functions
  _broadcast_getindex,
  broadcasted,
  broadcastable,
  instantiate

import ITensors.ContractionSequenceOptimization:
  optimized_contraction_sequence

import HDF5:
  read,
  write

import LinearAlgebra:
  axpby!,
  axpy!,
  dot,
  eigen,
  exp,
  factorize,
  lmul!,
  mul!,
  norm,
  normalize!,
  qr,
  rmul!,
  svd,
  tr

import NDTensors:
  # Modules
  Strided, # to control threading
  # Methods
  array,
  blockdim,
  blockoffsets,
  contract,
  dense,
  dim,
  dims,
  disable_tblis,
  eachnzblock,
  enable_tblis,
  ind,
  inds,
  insertblock!,
  matrix,
  #maxdim,
  mindim,
  nblocks,
  nnz,
  nnzblocks,
  nzblock,
  nzblocks,
  outer,
  permuteblocks,
  polar,
  scale!,
  setblockdim!,
  sim,
  store,
  sum,
  tensor,
  truncate!,
  using_tblis,
  vector,
  # Deprecated
  addblock!

import Random:
  randn!


text:
# override Documenter.jl/src/Writers/HTMLWriter.jl
using Dates
import Documenter: Documents, Documenter, Writers, Utilities
import Documenter.Utilities.DOM: DOM, Tag, @tags
import Documenter.Writers.HTMLWriter:
    getpage,
    render_head,
    render_article,
    open_output,
    mdconvert,
    pagetitle,
    navhref,
    domify,
    render_settings,
    THEMES

const zh_CN_numbers = ["一", "二", "三", "四", "五", "六",
    "七", "八", "九", "十", "十一", "十二"]
const zh_CN_months = zh_CN_numbers[1:12] .* "月"

const zh_CN_months_abbrev = zh_CN_months
const zh_CN_days = [
    "周" .* zh_CN_numbers[1:6]
    "周日"
]

Dates.LOCALES["zh_CN"] = Dates.DateLocale(zh_CN_months, zh_CN_months_abbrev, zh_CN_days, [""])

#=
  1. 修改编辑翻译提示语
    - 改 icon 为 FontAwesome: fa-globe(f0ac)
    - 改 "Edit on Github" 为 “完善 Transifex 上的翻译"
    - 改为对应的 Transifex URL
  2. 修改前后翻页提示语
    - "Previous" => t_Previous
    - "Next" => t_Next
=#

const t_Edit_on_xx = " 完善 Transifex 上的翻译" # 开头加空格，避免与 icon 靠太近
const t_Icon = "\uf0ac" # fa-globe

# 生成对应文件在 Transifex 上的翻译地址
function transifex_url(rel_path)
    # 首页源文件在 GitHub 上
    if splitdir(rel_path) == splitdir("src/index.md")
        return "https://github.com/JuliaCN/JuliaZH.jl/blob/master/doc/src/index.md"
    end

    # rel_path => "src\\manual\\getting-started.md"

    paths, file_name = splitdir(rel_path)
    # paths     => "src\\manual"
    # file_name => "getting-started.md"

    _, parent_fold = splitdir(paths)
    # parent_fold => "manual"

    page_name = replace(file_name, '.' => "")
    # page_name => "getting-startedmd"

    # final URL =>
    #   "https://www.transifex.com/juliacn/manual-zh_cn/translate/#zh_CN/getting-startedmd"
    "https://www.transifex.com/juliacn/$parent_fold-zh_cn/translate/#zh_CN/$page_name"
end

# workaround on Documenter/#977
function Documenter.Writers.HTMLWriter.render_navbar(ctx, navnode, edit_page_link::Bool)
    @tags div header nav ul li a span

    # The breadcrumb (navigation links on top)
    navpath = Documents.navpath(navnode)
    header_links = map(navpath) do nn
        title = mdconvert(pagetitle(ctx, nn); droplinks=true)
        nn.page === nothing ? li(a[".is-disabled"](title)) : li(a[:href => navhref(ctx, nn, navnode)](title))
    end
    header_links[end] = header_links[end][".is-active"]
    breadcrumb = nav[".breadcrumb"](
        ul[".is-hidden-mobile"](header_links),
        ul[".is-hidden-tablet"](header_links[end]) # when on mobile, we only show the page title, basically
    )

    # The "Edit on GitHub" links and the hamburger to open the sidebar (on mobile) float right
    navbar_right = div[".docs-right"]

    # Set the logo and name for the "Edit on.." button.
    if edit_page_link && (ctx.settings.edit_link !== nothing) && !ctx.settings.disable_git
        host_type = Utilities.repo_host_from_url(ctx.doc.user.repo)
        logo = t_Icon

        pageurl = get(getpage(ctx, navnode).globals.meta, :EditURL, getpage(ctx, navnode).source)
        edit_branch = isa(ctx.settings.edit_link, String) ? ctx.settings.edit_link : nothing

        # 1.3: 改 URL
        url = transifex_url(getpage(ctx, navnode).source)
        if url !== nothing
            title = t_Edit_on_xx
            push!(navbar_right.nodes,
                a[".docs-edit-link", :href => url, :title => title](
                    span[".docs-icon.fab"](logo),
                    span[".docs-label.is-hidden-touch"](title)
                )
            )
        end
    end

    # Settings cog
    push!(navbar_right.nodes, a[
        "#documenter-settings-button.docs-settings-button.fas.fa-cog",
        :href => "#", :title => "设置",
    ])

    # Hamburger on mobile
    push!(navbar_right.nodes, a[
        "#documenter-sidebar-button.docs-sidebar-button.fa.fa-bars.is-hidden-desktop",
        :href => "#"
    ])

    # Construct the main <header> node that should be the first element in div.docs-main
    header[".docs-navbar"](breadcrumb, navbar_right)
end

function render_settings(ctx)
    @tags div header section footer p button hr span select option label a

    theme_selector = p(
        label[".label"]("选择主题"),
        div[".select"](
            select["#documenter-themepicker"](option[:value=>theme](theme) for theme in THEMES)
        )
    )

    now_full = Dates.format(now(), "Y U d E HH:MM"; locale="zh_CN")
    now_short = Dates.format(now(), "Y U d E"; locale="zh_CN")
    buildinfo = p(
        "本文档在",
        span[".colophon-date", :title => now_full](now_short),
        "由",
        a[:href => "https://github.com/JuliaDocs/Documenter.jl"]("Documenter.jl"),
        "使用$(Base.VERSION)版本的Julia生成。"
    )

    div["#documenter-settings.modal"](
        div[".modal-background"],
        div[".modal-card"](
            header[".modal-card-head"](
                p[".modal-card-title"]("设置"),
                button[".delete"]()
            ),
            section[".modal-card-body"](
                theme_selector, hr(), buildinfo
            ),
            footer[".modal-card-foot"]()
        )
    )
end


text:
@extend struct ParametersV2016 <: Parameters
    pbacktime::Array{Float64,1} # Backstop price
    cpricebase::Array{Float64,1} # Carbon price in base case
    rr::Array{Float64,1} # Average utility social discount rate
    optlrsav::Float64 # Optimal savings rate
    ψ₂::JuMP.NonlinearParameter
    cumtree::Array{Float64,1} # Cumulative from land
end

function Base.show(io::IO, opt::ParametersV2016)
    ctx = IOContext(io, :limit => true, :compact => true);
    println(io, "Calculated Parameters for 2016 versions.");
    println(io, "Optimal savings rate: $(opt.optlrsav)");
    println(io, "Carbon cycle transition matrix coefficients");
    println(io, "ϕ₁₁: $(opt.ϕ₁₁), ϕ₂₁: $(opt.ϕ₂₁), ϕ₂₂: $(opt.ϕ₂₂), ϕ₃₂: $(opt.ϕ₃₂), ϕ₃₃: $(opt.ϕ₃₃)");
    println(io, "2015 Carbon intensity: $(opt.σ₀)");
    println(io, "Climate model parameter: $(opt.λ)");
    println(io, "Backstop price: ");
    show(ctx, opt.pbacktime);
    println(io, "");
    print(io, "Growth rate of productivity: ");
    show(ctx, opt.gₐ);
    println(io, "");
    print(io, "Emissions from deforestation: ");
    show(ctx, opt.Etree);
    println(io, "");
    print(io, "Avg utility social discout rate: ");
    show(ctx, opt.rr);
    println(io, "");
    print(io, "Base case carbon price: ");
    show(ctx, opt.cpricebase);
    println(io, "");
    print(io, "Population and labour: ");
    show(ctx, opt.L);
    println(io, "");
    print(io, "Total factor productivity: ");
    show(ctx, opt.A);
    println(io, "");
    print(io, "Δσ: ");
    show(ctx, opt.gσ);
    println(io, "");
    print(io, "σ: ");
    show(ctx, opt.σ);
    println(io, "");
    print(io, "cumtree: ");
    show(ctx, opt.cumtree);
    println(io, "");
    print(io, "θ₁: ");
    show(ctx, opt.θ₁);
    println(io, "");
    print(io, "Exogenious forcing: ");
    show(ctx, opt.fₑₓ);
end

@extend struct OptionsV2016R <: Options
    e₀::Float64 #Industrial emissions 2015 (GtCO2 per year)
    μ₀::Float64 #Initial emissions control rate for base case 2015
    tnopol::Float64 #Period before which no emissions controls base
    cprice₀::Float64 #Initial base carbon price (2010$ per tCO2)
    gcprice::Float64 #Growth rate of base carbon price per year
    ψ₁₀::Float64 #Initial damage intercept
end

@extend struct OptionsV2016R2 <: Options
    e₀::Float64 #Industrial emissions 2015 (GtCO2 per year)
    μ₀::Float64 #Initial emissions control rate for base case 2015
    tnopol::Float64 #Period before which no emissions controls base
    cprice₀::Float64 #Initial base carbon price (2010$ per tCO2)
    gcprice::Float64 #Growth rate of base carbon price per year
    ψ₁₀::Float64 #Initial damage intercept
end

@extend struct VariablesV2016 <: Variables
    Eind::Array{JuMP.VariableRef,1} # Industrial emissions (GtCO2 per year)
    Ω::Array{JuMP.VariableRef,1} # Damages as fraction of gross output
    Λ::Array{JuMP.VariableRef,1} # Cost of emissions reductions  (trillions 2005 USD per year)
    CPRICE::Array{JuMP.VariableRef,1} # Carbon price (2005$ per ton of CO2)
    CEMUTOTPER::Array{JuMP.VariableRef,1} # Period utility
    CCATOT::Array{JuMP.VariableRef,1} # Total carbon emissions (GTC)
end

@extend struct EquationsV2016 <: Equations
    cc::Array{JuMP.ConstraintRef,1} # Output Consumption
end

include("Scenarios2016.jl")
include("Results2016.jl")

include("2016R.jl")
include("2016R2.jl")