text:
@testset "4.4.0 (a trg)^m (b cot)^n" begin
    (a, b, c, d, e, f, m, n, p, x, ) = @variables a b c d e f m n p x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(b*cot(c+d*x))^n=#

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*cot(c+d*x)^n=#

    @test_int [cot(a + b*x)^1, x, 1, log(sin(a + b*x))/b]
    @test_int [cot(a + b*x)^2, x, 2, -x - cot(a + b*x)/b]
    @test_int [cot(a + b*x)^3, x, 2, -(cot(a + b*x)^2/(2*b)) - log(sin(a + b*x))/b]
    @test_int [cot(a + b*x)^4, x, 3, x + cot(a + b*x)/b - cot(a + b*x)^3/(3*b)]
    @test_int [cot(a + b*x)^5, x, 3, cot(a + b*x)^2/(2*b) - cot(a + b*x)^4/(4*b) + log(sin(a + b*x))/b]
    @test_int [cot(a + b*x)^6, x, 4, -x - cot(a + b*x)/b + cot(a + b*x)^3/(3*b) - cot(a + b*x)^5/(5*b)]
    @test_int [cot(a + b*x)^7, x, 4, -(cot(a + b*x)^2/(2*b)) + cot(a + b*x)^4/(4*b) - cot(a + b*x)^6/(6*b) - log(sin(a + b*x))/b]
    @test_int [cot(a + b*x)^8, x, 5, x + cot(a + b*x)/b - cot(a + b*x)^3/(3*b) + cot(a + b*x)^5/(5*b) - cot(a + b*x)^7/(7*b)]

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

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*cot(c+d*x))^(n/3)=#

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

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

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

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

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

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

    @test_int [(a*cot(x)^2)^(3/2),x, 3, (-(1/2))*a*cot(x)*sqrt(a*cot(x)^2) - a*sqrt(a*cot(x)^2)*log(sin(x))*tan(x)]
    @test_int [sqrt(a*cot(x)^2), x, 2, sqrt(a*cot(x)^2)*log(sin(x))*tan(x)]
    @test_int [1/sqrt(a*cot(x)^2), x, 2, -((cot(x)*log(cos(x)))/sqrt(a*cot(x)^2))]
    @test_int [1/(a*cot(x)^2)^(3/2),x, 3, (cot(x)*log(cos(x)))/(a*sqrt(a*cot(x)^2)) + tan(x)/(2*a*sqrt(a*cot(x)^2))]

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

    @test_int [(a*cot(x)^3)^(3/2), x, 14, (2/3)*a*sqrt(a*cot(x)^3) + (a*atan(1 - sqrt(2)*sqrt(cot(x)))*sqrt(a*cot(x)^3))/(sqrt(2)*cot(x)^(3/2)) - (a*atan(1 + sqrt(2)*sqrt(cot(x)))*sqrt(a*cot(x)^3))/(sqrt(2)*cot(x)^(3/2)) - (2/7)*a*cot(x)^2*sqrt(a*cot(x)^3) - (a*sqrt(a*cot(x)^3)*log(1 - sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*cot(x)^(3/2)) + (a*sqrt(a*cot(x)^3)*log(1 + sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*cot(x)^(3/2))]
    @test_int [sqrt(a*cot(x)^3), x, 13, -((atan(1 - sqrt(2)*sqrt(cot(x)))*sqrt(a*cot(x)^3))/(sqrt(2)*cot(x)^(3/2))) + (atan(1 + sqrt(2)*sqrt(cot(x)))*sqrt(a*cot(x)^3))/(sqrt(2)*cot(x)^(3/2)) - (sqrt(a*cot(x)^3)*log(1 - sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*cot(x)^(3/2)) + (sqrt(a*cot(x)^3)*log(1 + sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*cot(x)^(3/2)) - 2*sqrt(a*cot(x)^3)*tan(x)]
    @test_int [1/sqrt(a*cot(x)^3), x, 13, (2*cot(x))/sqrt(a*cot(x)^3) - (atan(1 - sqrt(2)*sqrt(cot(x)))*cot(x)^(3/2))/(sqrt(2)*sqrt(a*cot(x)^3)) + (atan(1 + sqrt(2)*sqrt(cot(x)))*cot(x)^(3/2))/(sqrt(2)*sqrt(a*cot(x)^3)) + (cot(x)^(3/2)*log(1 - sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*sqrt(a*cot(x)^3)) - (cot(x)^(3/2)*log(1 + sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*sqrt(a*cot(x)^3))]
    @test_int [1/(a*cot(x)^3)^(3/2),x, 14, -(2/(3*a*sqrt(a*cot(x)^3))) + (atan(1 - sqrt(2)*sqrt(cot(x)))*cot(x)^(3/2))/(sqrt(2)*a*sqrt(a*cot(x)^3)) - (atan(1 + sqrt(2)*sqrt(cot(x)))*cot(x)^(3/2))/(sqrt(2)*a*sqrt(a*cot(x)^3)) + (cot(x)^(3/2)*log(1 - sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*a*sqrt(a*cot(x)^3)) - (cot(x)^(3/2)*log(1 + sqrt(2)*sqrt(cot(x)) + cot(x)))/(2*sqrt(2)*a*sqrt(a*cot(x)^3)) + (2*tan(x)^2)/(7*a*sqrt(a*cot(x)^3))]

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

    @test_int [(a*cot(x)^4)^(3/2),x, 5, (1/3)*a*cot(x)*sqrt(a*cot(x)^4) - (1/5)*a*cot(x)^3*sqrt(a*cot(x)^4) - a*sqrt(a*cot(x)^4)*tan(x) - a*x*sqrt(a*cot(x)^4)*tan(x)^2]
    @test_int [sqrt(a*cot(x)^4), x, 3, (-sqrt(a*cot(x)^4))*tan(x) - x*sqrt(a*cot(x)^4)*tan(x)^2]
    @test_int [1/sqrt(a*cot(x)^4), x, 3, cot(x)/sqrt(a*cot(x)^4) - (x*cot(x)^2)/sqrt(a*cot(x)^4)]
    @test_int [1/(a*cot(x)^4)^(3/2),x, 5, cot(x)/(a*sqrt(a*cot(x)^4)) - (x*cot(x)^2)/(a*sqrt(a*cot(x)^4)) - tan(x)/(3*a*sqrt(a*cot(x)^4)) + tan(x)^3/(5*a*sqrt(a*cot(x)^4))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(b*cot(c+d*x)^p)^n*with*n*symbolic=#

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

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

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

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

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

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

    @test_int [(a*sin(e + f*x))^m*(b*cot(e + f*x))^n, x, 2, -(((b*cot(e + f*x))^(1 + n)*HypergeometricFunctions._₂F₁((1 + n)/2, (1/2)*(1 - m + n), (3 + n)/2, cos(e + f*x)^2)*(a*sin(e + f*x))^m*(sin(e + f*x)^2)^((1/2)*(1 - m + n)))/(b*f*(1 + n)))]

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

    @test_int [(a*cos(e + f*x))^m*(b*cot(e + f*x))^n, x, 2, -(((a*cos(e + f*x))^m*(b*cot(e + f*x))^(1 + n)*HypergeometricFunctions._₂F₁((1 + n)/2, (1/2)*(1 + m + n), (1/2)*(3 + m + n), cos(e + f*x)^2)*(sin(e + f*x)^2)^((1 + n)/2))/(b*f*(1 + m + n)))]

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

    @test_int [(a*cot(e + f*x))^m*(b*cot(e + f*x))^n, x, 3, -(((a*cot(e + f*x))^(1 + m)*(b*cot(e + f*x))^n*HypergeometricFunctions._₂F₁(1, (1/2)*(1 + m + n), (1/2)*(3 + m + n), -cot(e + f*x)^2))/(a*f*(1 + m + n)))]

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

    @test_int [(a*sec(e + f*x))^m*(b*cot(e + f*x))^n, x, 3, -(((b*cot(e + f*x))^(1 + n)*HypergeometricFunctions._₂F₁((1 + n)/2, (1/2)*(1 - m + n), (1/2)*(3 - m + n), cos(e + f*x)^2)*(a*sec(e + f*x))^m*(sin(e + f*x)^2)^((1 + n)/2))/(b*f*(1 - m + n)))]

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

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

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

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

    @test_int [csc(e + f*x)^6*(d*cot(e + f*x))^n, x, 3, -((d*cot(e + f*x))^(1 + n)/(d*f*(1 + n))) - (2*(d*cot(e + f*x))^(3 + n))/(d^3*f*(3 + n)) - (d*cot(e + f*x))^(5 + n)/(d^5*f*(5 + n))]
    @test_int [csc(e + f*x)^4*(d*cot(e + f*x))^n, x, 3, -((d*cot(e + f*x))^(1 + n)/(d*f*(1 + n))) - (d*cot(e + f*x))^(3 + n)/(d^3*f*(3 + n))]
    @test_int [csc(e + f*x)^2*(d*cot(e + f*x))^n, x, 2, -((d*cot(e + f*x))^(1 + n)/(d*f*(1 + n)))]
    @test_int [sin(e + f*x)^2*(d*cot(e + f*x))^n, x, 2, -(((d*cot(e + f*x))^(1 + n)*HypergeometricFunctions._₂F₁(2, (1 + n)/2, (3 + n)/2, -cot(e + f*x)^2))/(d*f*(1 + n)))]
    @test_int [sin(e + f*x)^4*(d*cot(e + f*x))^n, x, 2, -(((d*cot(e + f*x))^(1 + n)*HypergeometricFunctions._₂F₁(3, (1 + n)/2, (3 + n)/2, -cot(e + f*x)^2))/(d*f*(1 + n)))]

    @test_int [csc(e + f*x)^3*(d*cot(e + f*x))^n, x, 1, -(((d*cot(e + f*x))^(1 + n)*csc(e + f*x)^3*HypergeometricFunctions._₂F₁((1 + n)/2, (4 + n)/2, (3 + n)/2, cos(e + f*x)^2)*(sin(e + f*x)^2)^((4 + n)/2))/(d*f*(1 + n)))]
    @test_int [csc(e + f*x)^1*(d*cot(e + f*x))^n, x, 1, -(((d*cot(e + f*x))^(1 + n)*csc(e + f*x)*HypergeometricFunctions._₂F₁((1 + n)/2, (2 + n)/2, (3 + n)/2, cos(e + f*x)^2)*(sin(e + f*x)^2)^((2 + n)/2))/(d*f*(1 + n)))]
    @test_int [sin(e + f*x)^1*(d*cot(e + f*x))^n, x, 1, -(((d*cot(e + f*x))^(1 + n)*HypergeometricFunctions._₂F₁(n/2, (1 + n)/2, (3 + n)/2, cos(e + f*x)^2)*sin(e + f*x)*(sin(e + f*x)^2)^(n/2))/(d*f*(1 + n)))]
    @test_int [sin(e + f*x)^3*(d*cot(e + f*x))^n, x, 1, -(((d*cot(e + f*x))^(1 + n)*HypergeometricFunctions._₂F₁((1/2)*(-2 + n), (1 + n)/2, (3 + n)/2, cos(e + f*x)^2)*sin(e + f*x)^3*(sin(e + f*x)^2)^((1/2)*(-2 + n)))/(d*f*(1 + n)))]

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

    @test_int [(a*csc(e + f*x))^m*(b*cot(e + f*x))^n, x, 1, -(((b*cot(e + f*x))^(1 + n)*(a*csc(e + f*x))^m*HypergeometricFunctions._₂F₁((1 + n)/2, (1/2)*(1 + m + n), (3 + n)/2, cos(e + f*x)^2)*(sin(e + f*x)^2)^((1/2)*(1 + m + n)))/(b*f*(1 + n)))]
    end

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.0 (a trg)^m (b cot)^n.jl


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 9f4d1108-1fed-4314-ae4e-560ec7a9be9c
using PlutoTest

# ╔═╡ 76d27f60-b0e8-11eb-3fd3-a7b73d0d4cfe
md"""
# [Day 14: Reindeer Olympics](https://adventofcode.com/2015/day/14)

This year is the Reindeer Olympics! Reindeer can fly at high speeds, but must rest occasionally to recover their energy. Santa would like to know which of his reindeer is fastest, and so he has them race.

Reindeer can only either be flying (always at their top speed) or resting (not moving at all), and always spend whole seconds in either state.

For example, suppose you have the following Reindeer:

-   Comet can fly 14 km/s for 10 seconds, but then must rest for 127 seconds.
-   Dancer can fly 16 km/s for 11 seconds, but then must rest for 162 seconds.

After one second, Comet has gone 14 km, while Dancer has gone 16 km. After ten seconds, Comet has gone 140 km, while Dancer has gone 160 km. On the eleventh second, Comet begins resting (staying at 140 km), and Dancer continues on for a total distance of 176 km. On the 12th second, both reindeer are resting. They continue to rest until the 138th second, when Comet flies for another ten seconds. On the 174th second, Dancer flies for another 11 seconds.

In this example, after the 1000th second, both reindeer are resting, and Comet is in the lead at `1120` km (poor Dancer has only gotten `1056` km by that point). So, in this situation, Comet would win (if the race ended at 1000 seconds).

Given the descriptions of each reindeer (in your puzzle input), after exactly `2503` seconds, **what distance has the winning reindeer traveled?**
"""

# ╔═╡ 7b3cffcf-5c04-4d13-9714-8dad9e4956f3
puzzle_input = "Vixen can fly 19 km/s for 7 seconds, but then must rest for 124 seconds.
Rudolph can fly 3 km/s for 15 seconds, but then must rest for 28 seconds.
Donner can fly 19 km/s for 9 seconds, but then must rest for 164 seconds.
Blitzen can fly 19 km/s for 9 seconds, but then must rest for 158 seconds.
Comet can fly 13 km/s for 7 seconds, but then must rest for 82 seconds.
Cupid can fly 25 km/s for 6 seconds, but then must rest for 145 seconds.
Dasher can fly 14 km/s for 3 seconds, but then must rest for 38 seconds.
Dancer can fly 3 km/s for 16 seconds, but then must rest for 37 seconds.
Prancer can fly 25 km/s for 6 seconds, but then must rest for 143 seconds."

# ╔═╡ 4c30acd3-078f-40a6-8e22-c90704abef38
function parseline(line)
	m = match(r"(\w+) can fly (\d+) km/s for (\d+) seconds, but then must rest for (\d+) seconds.", line)
	(km=parse(Int, m[2]), fly_seconds=parse(Int, m[3]), rest_seconds=parse(Int, m[4]))
end

# ╔═╡ 6259b41c-e73c-46b6-bf76-2f42c6428003
example_players = map(parseline, ["Comet can fly 14 km/s for 10 seconds, but then must rest for 127 seconds.", "Dancer can fly 16 km/s for 11 seconds, but then must rest for 162 seconds."])

# ╔═╡ 21d4c9cf-2502-4522-90fc-5136e5bd13bd
function distance(player, time)
	d = 0
	#while time > 0
	#	stepTime = min(fly_seconds, time)
	#	d += stepTime * km
	#	time -= stepTime + rest_seconds
	#end

	d = time ÷ (player.fly_seconds + player.rest_seconds) * (player.fly_seconds * player.km)
	d += min(time % (player.fly_seconds + player.rest_seconds), player.fly_seconds) * player.km
	d
end

# ╔═╡ bec1d683-76e6-4b80-968c-4179944073cc
maxdistance(players, time) = maximum(map(player->distance(player, time), players))

# ╔═╡ cac38821-4cef-4d15-a86b-7878a37392da
@assert maxdistance(example_players, 1000) == 1120

# ╔═╡ 315d6498-c7d6-44f6-876f-3cbc599d7393
puzzle_players = map(parseline, split(puzzle_input, "\n"))

# ╔═╡ 598ca1b1-acfc-449d-be51-6ae63efaeaa7
part1 = maxdistance(puzzle_players, 2503)

# ╔═╡ 40ae0576-75ae-4d63-b163-4fbefe0e0060
@test part1 == 2660

# ╔═╡ 17210bdf-21cd-4f48-8e86-b7c568ee3e2d
md"Your puzzle answer was `2660`."

# ╔═╡ 0975d5a4-675e-42b8-8ab1-4b31dc2de102
md"""
# Part Two

Seeing how reindeer move in bursts, Santa decides he's not pleased with the old scoring system.

Instead, at the end of each second, he awards one point to the reindeer currently in the lead. (If there are multiple reindeer tied for the lead, they each get one point.) He keeps the traditional 2503 second time limit, of course, as doing otherwise would be entirely ridiculous.

Given the example reindeer from above, after the first second, Dancer is in the lead and gets one point. He stays in the lead until several seconds into Comet's second burst: after the 140th second, Comet pulls into the lead and gets his first point. Of course, since Dancer had been in the lead for the 139 seconds before that, he has accumulated 139 points by the 140th second.

After the 1000th second, Dancer has accumulated `689` points, while poor Comet, our old champion, only has `312`. So, with the new scoring system, Dancer would win (if the race ended at 1000 seconds).

Again given the descriptions of each reindeer (in your puzzle input), after exactly `2503` seconds, **how many points does the winning reindeer have?**
"""

# ╔═╡ 9cd8509d-38c3-49d0-bc6d-84c862bd37e0
function maxscore(players, time)
	scores = Dict()

	for idx in 1:length(players)
		scores[idx] = 0
	end

	for t in 1:time
		distances = Dict()
		for idx in 1:length(players)
			distances[idx] = distance(players[idx], t)
		end
		maxdistance = maximum(values(distances))
		winners = filter(idx-> distances[idx] == maxdistance, keys(distances))
		for idx in winners
			scores[idx] += 1
		end
	end
	maximum(values(scores))

end

# ╔═╡ 14d5831a-f848-426d-a014-2a1e24ee1d68
@assert maxscore(example_players, 1000) == 689

# ╔═╡ 0dec0d65-0cc9-4d7c-bdf0-d2bbc3f6c956
part2 = maxscore(puzzle_players, 2503)

# ╔═╡ 09982bbe-252b-40e2-b818-c1f2f575f458
@test part2 == 1256

# ╔═╡ fa92058e-35cf-43b8-8552-71dd16e1c52e
md"Your puzzle answer was `1256`."

# ╔═╡ Cell order:
# ╠═9f4d1108-1fed-4314-ae4e-560ec7a9be9c
# ╟─76d27f60-b0e8-11eb-3fd3-a7b73d0d4cfe
# ╟─7b3cffcf-5c04-4d13-9714-8dad9e4956f3
# ╠═4c30acd3-078f-40a6-8e22-c90704abef38
# ╠═6259b41c-e73c-46b6-bf76-2f42c6428003
# ╠═21d4c9cf-2502-4522-90fc-5136e5bd13bd
# ╠═bec1d683-76e6-4b80-968c-4179944073cc
# ╠═cac38821-4cef-4d15-a86b-7878a37392da
# ╠═315d6498-c7d6-44f6-876f-3cbc599d7393
# ╠═598ca1b1-acfc-449d-be51-6ae63efaeaa7
# ╠═40ae0576-75ae-4d63-b163-4fbefe0e0060
# ╟─17210bdf-21cd-4f48-8e86-b7c568ee3e2d
# ╟─0975d5a4-675e-42b8-8ab1-4b31dc2de102
# ╠═9cd8509d-38c3-49d0-bc6d-84c862bd37e0
# ╠═14d5831a-f848-426d-a014-2a1e24ee1d68
# ╠═0dec0d65-0cc9-4d7c-bdf0-d2bbc3f6c956
# ╠═09982bbe-252b-40e2-b818-c1f2f575f458
# ╟─fa92058e-35cf-43b8-8552-71dd16e1c52e

ext: jl


text:
module RealRoots
using Random
using LinearAlgebra
import Roots: find_zero, AlefeldPotraShi
using Printf

## Find real roots of a polynomial using a DesCartes Method
##
## State of the art for this approach is summarized here:
## Cf. "Computing Real Roots of Real Polynomials ... and now For Real!" by Alexander Kobel, Fabrice Rouillier, Michael Rouillier
## (https://arxiv.org/abs/1605.00410)
##
## Earlier work of a Descartes-like method  [Vincent's theorem](http://en.wikipedia.org/wiki/Vincent's_theorem)
## are here "Efficient isolation of polynomial’s real roots" by
## Fabrice Rouillier; Paul Zimmermann
##
## This implementation doesn't take nearly enough care with the details, but takes some ideas
## to implement a means to find real roots of non-pathological polynomials (lowish degree, roots separated)
##
## XXX Needs work XXX

## Polynomial transformations
##
## The Taylor shift here is the most expensive operation
## https://arxiv.org/pdf/1605.00410.pdf has a better strategy
## of partial Taylor shifts, using just the nearby roots
##

using Polynomials
import Roots: fzero
using ..AGCD

# for interface
struct PolyType{T}
    p::Vector{T}
end
(p::PolyType)(x) = poly_eval(p.p, x)

Polynomials.degree(p::Vector{T}) where {T} = findlast(!iszero,p) - 1

## Interval with [a,b], N
## N is for Newton-Test
mutable struct Intervalab{T,S}
    a::T
    b::S
    N::Int
    bnd::Int
end
width(I::Intervalab) = I.b-I.a
midpoint(I::Intervalab) = I.a + (0.5) * width(I)

## how we store the state of our algorithm to find zero in [0,1]
struct State{T}
    Internal::Vector{Intervalab{T,T}}                    # DesBound > 1
    Isol::Vector{Intervalab{T,T}}                        # DesBound == 1
    Unresolved::Vector{Intervalab{T,T}}
    p::Vector{T}
end

State(p::Vector{T}) where {T} = State(Intervalab{T,T}[], Intervalab{T,T}[], Intervalab{T,T}[], p)
(st::State)(x) = poly_eval(st.p, x)

# polys over vectors mutating first

# poly is [p_0, p_1, ..., p_n] where there may be zeros in last terms

# p + q
function poly_add!(p::Vector{T}, q::Vector{S}) where {T, S}
    n,m = length(p), length(q)
    l = min(n, m)
    for i in 1:l
        p[i] += q[i]
    end
    for i in l:m
        push!(p, q[i])
    end
end

# p-q
poly_sub!(p::Vector{T}, q::Vector{S}) where {T, S} = poly_add(p, -q)

# may be more than one way
function poly_mul!(p1::Vector{T}, p2::Vector{S}) where {T, S}
    R = promote_type(T,S)

    n = length(p1)
    m = length(p2)

    a = zeros(R,m+n+1)

    for i = 1:n
        for j = 1:m
            a[i+j+1] += p1[i] * p2[j]
        end
    end
    p1[:] = a

end

function poly_deriv!(p::Vector{T}) where {T}
    for i = 1:length(p)-1
        p[i] = i * p[i+1]
    end
    p[end] = zero(T)
    p
end

function poly_eval(p::Vector{T}, x::S) where {T, S}
    R = promote_type(T,S)

    y = convert(R, p[end])

    for i = (lastindex(p)-1):-1:1
        y = p[i] + x*y
    end
    return y
end

" `p(x + λ)`: Translate polynomial left by λ "
function poly_translate!(p::Vector{T}, lambda=1) where {T}
    p1 = copy(p)
    p[1] = poly_eval(p1, lambda)
    m = one(T)
    for k in 2:length(p1)
        p[k] = poly_eval(poly_deriv!(p1), lambda) / m
        m *= k
    end
    p
end
Tλ(p, lambda=1)   = poly_translate!(copy(p), lambda)

" `R(p)` finds  `x^n p(1/x)` which is a reversal of coefficients "
function poly_reverse!(p)
    reverse!(p)
end
R(p) = poly_reverse!(copy(p))

" `p(λ x)`: scale x axis by λ "
function poly_scale!(p, lambda)
    for i in 2:length(p)
        p[i] *= lambda^(i-1)
    end
end
Hλ(p::Vector{T}, λ=one(T)) where {T} = poly_scale!(copy(p), λ)

"q = p(-x)"
function poly_flip!(p::Vector{T}) where {T}
    for i in 2:2:length(p)
        p[i] = -p[i]
    end
end

## Upper bound on size of real roots that is tighter than cauchy
## titan.princeton.edu/papers/claire/hertz-etal-99.ps
function upperbound(p::Vector{T}) where {T}
    p = p[findfirst(!iszero, p):end]
    descartes_bound(p) == 0 && return zero(T)

    p = p[findfirst(!iszero, p):end]

    q, d = p/p[end], length(p)-1

    d == 0 && error("degree 0 is a constant")
    d == 1 && abs(q[1])

    a1 = abs(q[d])
    B = maximum([abs(q[i]) for i in 1:(d-1)])

    a,b,c = 1, -(1+a1), a1-B
    (-b + sqrt(b^2 - 4a*c))/2
end

function lowerbound(p::Vector{T}) where {T <: Real}
    p = p[findfirst(!iszero, p):end]

    poly_flip!(p)
    ret = -upperbound(p)
    poly_flip!(p)
    ret
end

## Descartes Bounds

" Descartes bound on (0,oo). Just count sign changes"
function descartes_bound(p::Vector{T}) where {T}
    length(p) == 0 && return -1
    cnt, sgn = 0, sign(p[1])
    for i in 2:length(p)
        nsgn = sign(p[i])
        if nsgn * sgn < 0
            sgn = nsgn
            cnt += 1
        end
    end
    cnt
end

# shift polynomial so (a,b) -> (0, oo)
function translate_ab(p::Vector{T}, a, b) where {T}
    #    p0 = Tλ(p, a)
    #    p1 = Hλ(p0, b-a)
    #    p2 = Tλ(R(p1),1)

    p1 = copy(p)
    poly_translate!(p1, a)
    poly_scale!(p1, b-a)
    poly_reverse!(p1)
    poly_translate!(p1, one(T))
    p1
end

" Descartes bound on (a, b)"
descartes_bound_ab(p::Vector{T}, a, b) where {T} = descartes_bound(translate_ab(p, a, b))
DescartesBound_ab(st, node) = descartes_bound_ab(st.p, node.a, node.b)

## Tests

## Return true or false
zero_test(st::State, node)  =   DescartesBound_ab(st, node) == 0
one_test(st::State, node)  =   DescartesBound_ab(st, node) == 1

## return count -1 (can't tell), 0, 1, or more
zero_one_test(st::State, node) = DescartesBound_ab(st, node)

# find admissible point
# XXX improve me
function find_admissible_point(st::State{T},  I::Intervalab, m=midpoint(I), Ni::T=one(T), c::T=one(T)) where {T}
    N = ceil(Int, c * Polynomials.degree(st.p)/2)
    ep = min(m-I.a, I.b - m) / (4*Ni)
    mis = [m + i/N * ep for i in -N:N]
    curmin = min(norm(st(I.a)), norm(st(I.b)))/100
    for m in shuffle(mis)
        (m < I.b || m > I.a) || continue
        descartes_bound_ab(st.p, I.a, m) == -1 && continue
        descartes_bound_ab(st.p, m, I.b) == -1 && continue
        norm(st(m)) > 0 && return m
    end

    error("No admissible point found")
#    mx, i = findmax(norm.(st.p.(mis)))
#    mis[i]
end

# find splitting point
# find an admissible point that can be used to split interval. Needs to have all
# coefficients not straddle 0
# return (logical, two intervals)
function split_interval(st::State{T},I::Intervalab,  m=midpoint(I), Ni=one(T), c=one(T)) where {T}
    N = ceil(Int, c * Polynomials.degree(st.p)/2)
    ep = min(1, width(I)) / (16*Ni)
    mis = T[m + i/N * ep for i in -N:N]
    mis = filter(m -> m > I.a && m < I.b, mis)
    mx, i = findmax(norm.(st.(mis)))

    ## Now, we need a point that is bigger than max and leaves conclusive
    for i in eachindex(mis)
        mi = mis[i]
        abs(st(mi)) >= min(mx/4, one(T)) || continue
        ileft = Intervalab(I.a, mi, I.N, -2)
        nl = DescartesBound_ab(st, ileft)
        #nl = descartes_bound_ab(fatten(st.p), ileft.a, ileft.b)
        nl == -1 && continue
        iright = Intervalab(mi, I.b, I.N, -2)
        nr = DescartesBound_ab(st, iright)
        #nr = descartes_bound_ab(fatten(st.p), iright.a, iright.b)
        nr == -1 && continue

        # XXX improve this XXX
        # identify degenerate cases here. This is not good.
        # if nl + nr < I.bnd
        #     if nl == 0 || nr == 0
        #         if (nl == 0 && rem(I.bnd - nr,2) == 1) || ( nl == 0 && rem(I.bnd - nr,2) == 1)
        #             println("nl=$nl, nr=$nr, bnd=$(I.bnd) -- is this an error")
        #             return(false, I, I)
        #         end
        #     elseif nl + nr < I.bnd
        #         println("nl=$nl, nr=$nr, bnd=$(I.bnd) -- is this an error")
        #         return(false, I, I)
        #     end
        # end

        ileft.bnd = nl; iright.bnd=nr

##        println("Split $(I.a) < $mi < $(I.b): $(nl) & $(nr)")

        return (true, ileft, iright)
    end
    #    println("DEBUG: $I is a bad interval for splitting?")
    return (false, I, I)
end

## split interval
## adds to intervals if successful
## adds node to Unresolved if not
function linear_step(st::State{T}, node) where {T}
    succ, I1, I2 = split_interval(st, node)

    if succ
        push!(st.Internal, I1)
        push!(st.Internal, I2)
    else
        push!(st.Unresolved, node)
    end
    return true

end

## return (true, I), (false, node)
## I will be smaller interval containing all roots in node
function newton_test(st::State{T}, node) where {T}
    NMAX = 1024  # 2147483648 = 2^31
    (node.N > NMAX) && return (false, node)
    (zero_one_test(st, node) in (0,1)) && return(false, node)

    a, b, m, w, N, bnd  = node.a, node.b, midpoint(node), width(node), node.N, node.bnd

    pprime = poly_deriv!(copy(st.p))

    a1 = a - st(a) / poly_eval(pprime, a)
    b1 = b - st(b) / poly_eval(pprime, b)

    if a < a1 && zero_test(st, Intervalab(a, a1, N, 0))
        if b1 < b && zero_test(st, Intervalab(b1, b, N, 0))
            return (true, Intervalab(a1, b1, N*N, -2))
        else
            return(true, Intervalab(b1, b, N*N, -2))
        end
    elseif b1 < b && zero_test(st, Intervalab(b1, b, N, 0))
        return (true, Intervalab(a, b, N*N, -2))
    end

    ## boundary test?

    mlstar::T = find_admissible_point(st, node, a + w/(2N))
    if mlstar > a && zero_test(st, Intervalab(mlstar, b, N,-2))
        return (true, Intervalab(a, mlstar, N,-2))
    end

    mrstar::T = find_admissible_point(st, node, b - w/(2N))
    if mrstar < b && zero_test(st, Intervalab(a, mrstar, N,-2))
        return (true, Intervalab(mrstar, b, N,-2))
    end
    return (false, node)
end

## Add successors to I
## We have
function addSucc(st::State{T}, node) where {T}

    val, I = newton_test(st, node)

    if val
        ## a,b = node.a, node.b
        ## println("Newton test: a=$a <= $(I.a) <= $(I.b) <= $b=b")
        push!(st.Internal, I)
    else
        ##        println("linear step")
        succ = linear_step(st, node)
        if !succ
            warn("node $node was a failure")
        end
    end
    true
end

## m, M should bound the roots
## essentially algorithm 4
function ANewDsc(p::Vector{T}, m = lowerbound(p), M=upperbound(p)) where {T <: Real}

    st = State(p)

    base_node = Intervalab(m, M, 4, -2)
    base_node.bnd = DescartesBound_ab(st, base_node)

    if base_node.bnd == -1
        append!(st.Internal, break_up_interval(st, base_node, 4))
    else
        push!(st.Internal, base_node)
    end

    while length(st.Internal) > 0
        node = pop!(st.Internal)

        bnd = node.bnd
        if bnd == -2
            bnd = DescartesBound_ab(st, node)
            node.bnd = bnd
        end

        if bnd < 0
            # this is a bad node!
            warn("Bad node, how did it get here: $node")
        elseif bnd == 0
            continue
        elseif bnd == 1
            push!(st.Isol, node)
            continue
        else
            addSucc(st, node)
        end
    end
    st
end

# populate `Isol`
# p must not have any roots with even degree. (e.g. no (x-c)^(2k) exists as a factor for any c,k
# assumed square free (at least no roots of even multiplicity)
function isolate_roots(p::Vector{T}, m, M) where {T <: Real}

#    try
    st = ANewDsc(p, m, M)
    return st
    # catch err
    #     if  !(T <: BigFloat)
    #         try
    #             st = ANewDsc(convert(Poly{BigFloat}, p), m, M)
    #             return st
    #         catch err
    #             rethrow(err)
    #         end
    #     end
    # end

end

"""

     real_roots(p, [m], [M]; square_free=true)

Returns real roots of a polynomial presented via its coefficients `[p_0, p_1, ..., p_n]`.

* `p`: polynomial coefficients, `Vector{T<:Real}`
* `m`: lower bound on real roots. Defaults to `lowerbound(p)`
* `M`: upper bound on real roots. Defaults to `upperbound(p)`
* `square_free`::Bool. If false, the polynomial `agcd(p, polyder(p))` is used. This polynomial---in theory--- would have the
same real roots as `p`, however in practice the approximate `gcd` can be off.

"""
real_roots(p::Poly{T}, args...; kwargs...) where {T <: Real} = real_roots(p.a, args...; kwargs...)
function real_roots(p::Vector{T}, m = lowerbound(p), M=upperbound(p); square_free::Bool=false) where {T <: Real}

    # deflate zero
    nzroots = 0
    while iszero(p[1])
        popfirst!(p)
        nzroots += 1
    end

    if !square_free
        u,v,w,err = AGCD.agcd(Poly(p), polyder(Poly(p)))
        p = v.a
    end

    st = isolate_roots(p, m, M)

    if length(st.Unresolved) > 0
        println("Some intervals are unresolved:")
        println("------------------------------")
        for node in st.Unresolved
            @printf "* There may be up to %d roots in (%0.16f, %0.16f).\n" node.bnd node.a node.b
        end
        println("------------------------------")
    end

    rts = zeros(T, length(st.Isol))
    for i in eachindex(st.Isol)
        node = st.Isol[i]
        rt = find_zero(PolyType(p), (node.a, node.b), AlefeldPotraShi())
        #find_zero(PolyType(p), big.((node.a, node.b)), AlefeldPotraShi())
        rts[i] = rt
    end

    nzroots > 0 && push!(rts, zero(T))
    rts
end

end


text:
abstract type AbstractSolution end
abstract type AbstractAlgorithm end
abstract type AbstractParameters end
abstract type AbstractProblem end


text:
################# Parser functions ######################

function open_model_file(path::Q) where {Q <: AbstractString}

    #= Opens the model file and reads the contents, which are stored
       in a vector of strings. =#

    model_file = open(path)
    model_array = readlines(model_file)
    close(model_file)

    # Remove lines or parts of lines that have been commented out.

    for i in eachindex(model_array)
        if occursin("#",model_array[i]) == true
            components = split(model_array[i],"#")
            keep = prod(components[isodd.(eachindex(components))])
            model_array[i] = keep
        end
    end

    # Remove blank lines

    model_array = model_array[model_array .!= ""]

    return model_array

end

function find_term(model_array::Array{Q,1},term::Q) where {Q <: AbstractString}

    #= Finds the position in the String-vector where model-terms (states, jumps,
       parameters, equations, etc) are located. =#

    locations = findall(y->contains(y,term), model_array)
    if length(locations) == 1
        return locations[1]
    elseif length(locations) == 0
        error("The $term-designation does not appear in the model file.")
    else
        error("The $term-designation appears multiple times in the model file.")
    end

end

function find_end(model_array::Array{Q,1},startfrom::S) where {Q <: AbstractString,S <: Integer}

    #= Each of the key model terms must be followed by an 'end', whose location
       in the vector this function finds.  White space at the start of a line is
       stripped out. =#

    for end_location = startfrom:length(model_array)
        if startswith(strip(model_array[end_location]), "end") == true
            return end_location
        end
    end

    error("Model file is missing an 'end' term.")

end

function get_variables(model_array::Array{Q,1},term::Q) where {Q <: AbstractString}

    #= This function extracts the variable names and ensures that no names
       are repeated =#

    term_begin = find_term(model_array, term) + 1
    term_end   = find_end(model_array, term_begin) - 1

    if term_begin > term_end
        if term in ["shocks:","states:"]
            return String[]
        else
            error("The model file contains no $(term[1:end-1]).")
        end
    end

    term_block = model_array[term_begin:term_end]

    # Extract the variable names

    # Remove any trailing variable separators: "," or ";".

    for i = 1:length(term_block)
        if endswith(term_block[i],union(",",";"))
            term_block[i] = term_block[i][1:end-1]
        end
    end

    # Extract the names and place them in a vector

    variables = String.(setdiff(strip.(split(term_block[1], union(",",";"))),[""]))
    for i = 2:length(term_block)
        variables = [variables;String.(setdiff(strip.(split(term_block[i], union(",",";"))),[""]))]
    end

    # Check whether names are repeated

    if length(variables) != length(unique(variables))
        error("Some $(term[1:end-1]) are repreated.")
    end

    # Check to ensure that variables contains non-empty string elements

    if length(variables) == 1 && variables[1] == ""
        if term in ["shocks:","states:"]
            return String[]
        else
            error("The model file contains no $(term[1:end-1]).")
        end
    else
        return variables
    end

end

function combine_variables(x::Array{Q,1},y::Array{Q,1}) where {Q <: AbstractString}

    #= This function combines the states, "x", with the jump variables, "y", to
       generate the model's variables. =#

    if length(x) == 0 # There are no states
        return y
    elseif length(intersect(x,y)) > 0
        error("Some states and jumps have the same name.")
    else
        return [x;y]
    end

end

function get_parameters_and_values(model_array::Array{Q,1},term::Q) where {Q <: AbstractString}

    #= This function extracts the names and associated values for each of the
       model's parameters.  The parameter names are sorted so that larger names
       come first. =#

    parametersbegin = find_term(model_array, term) + 1
    parametersend   = find_end(model_array, parametersbegin) - 1
    if parametersbegin > parametersend
        error("The model file contains no $(term[1:end-1])")
    end

    parameterblock = model_array[parametersbegin:parametersend]

    # Extract the parameter names and values

    # Remove any trailing separators: "," or ";".

    for i = 1:length(parameterblock)
        if endswith(parameterblock[i],union(",",";"))
            parameterblock[i] = parameterblock[i][1:end-1]
        end
    end

    revised_parameterblock = String.(strip.(split(parameterblock[1], union(",",";"))))
    for i = 2:length(parameterblock)
        revised_parameterblock = [revised_parameterblock;String.(strip.(split(parameterblock[i], union(",",";"))))]
    end

    # Extract the parameter names and values

    unassigned_parameter_index = 1
    unassigned_parameters = Array{Q}(undef,0)

    parameters = Array{Q}(undef,length(revised_parameterblock))
    values     = Array{Q}(undef,length(revised_parameterblock))
    for i = 1:length(revised_parameterblock)
        if occursin("=",revised_parameterblock[i]) == false
            parameters[i] = revised_parameterblock[i]
            values[i] = "p[$unassigned_parameter_index]" # p is a reserved name
            push!(unassigned_parameters,revised_parameterblock[i])
            unassigned_parameter_index += 1
        else
            pair = strip.(split(revised_parameterblock[i],"="))
            parameters[i] = pair[1]
            values[i] = pair[2]
        end
    end

    parameter_order   = sortperm(length.(parameters),rev = true)
    sorted_parameters = parameters[parameter_order]
    sorted_values     = values[parameter_order]

    return sorted_parameters, sorted_values, unassigned_parameters

end

function get_equations(model_array::Array{Q,1},term::Q) where {Q <: AbstractString}

    # Extract the model's equations.

    equationsbegin  = find_term(model_array, term) + 1
    equationsend    = find_end(model_array, equationsbegin) - 1
    if equationsbegin > equationsend
        error("The model file contains no $(term[1:end-1])")
    end

    equation_block  = model_array[equationsbegin:equationsend]

    # Extract the equations

    # Remove any trailing separators: "," or ";".

    for i = 1:length(equation_block)
        if endswith(equation_block[i],union(",",";"))
            equation_block[i] = equation_block[i][1:end-1]
        end
    end

    # Extract the equations and place them in a vector

    equations = String.(strip.(split(equation_block[1], union(",",";"))))
    for i = 2:length(equation_block)
        equations = [equations;String.(strip.(split(equation_block[i], union(",",";"))))]
    end

    # For every model equation...

    for i = 1:length(equation_block)
        if occursin("[",equations[i]) == true # Replace open square bracket with open round parenthesis
            equations[i] = replace(equations[i],"["=>"(")
        elseif occursin("]",equations[i]) == true # Replace close square bracket with close round parenthesis
            equations[i] = replace(equations[i],"]"=>")")
        elseif occursin("{",equations[i]) == true # Replace open curly brace with open round parenthesis
            equations[i] = replace(equations[i],"{"=>"(")
        elseif occursin("}",equations[i]) == true # Replace close curly brace with close round parenthesis
            equations[i] = replace(equations[i],"}"=>")")
        end
        if occursin("=",equations[i]) == false # Check that each equation contains an equals sign
            error("Equation line $i does not contain an '=' sign.")
        elseif length(findall("(",equations[i])) != length(findall(")",equations[i]))
            error("Equation line $i has unbalanced parentheses")
        end
    end

    return equations

end

function reorder_equations(equations::Array{Q,1},shocks::Array{Q,1},states::Array{Q,1},jumps::Array{Q,1},parameters::Array{Q,1}) where {Q <: AbstractString}

    if isempty(shocks) == false # Case for stochastic models
        reordered_equations, reordered_states, reordered_shocks = reorder_equations_stochastic(equations,shocks,states,jumps,parameters)
        return reordered_equations, reordered_states, reordered_shocks
    else # Case for deterministic models
        reordered_equations, reordered_states = reorder_equations_deterministic(equations,states,jumps,parameters)
        return reordered_equations, reordered_states, shocks
    end

end

function reorder_equations_deterministic(equations::Array{Q,1},states::Array{Q,1},jumps::Array{Q,1},parameters::Array{Q,1}) where {Q <: AbstractString}

    #= This function reorders the model's equations so that the equations
       containing the shock processes appear first. =#

    reordered_equations = copy(equations)
    reordered_states    = copy(states)

    # Construct summary information about each equation

    trash_equations = copy(equations)

    variables = combine_variables(combine_variables(combine_variables(states,jumps),parameters),["exp","log"])
    variables = variables[sortperm(length.(variables),rev = true)] # sort variables so that longer names are first

    states_number = zeros(Int64,length(equations))
    jumps_number  = zeros(Int64,length(equations))
    for i = 1:length(equations)
      for j in variables
        if j in states
            states_number[i] += occursin(j,trash_equations[i])
            trash_equations[i] = replace(trash_equations[i],j => ":")
        elseif j in jumps
            jumps_number[i] += occursin(j,trash_equations[i])
            trash_equations[i] = replace(trash_equations[i],j => ":")
        else
            trash_equations[i] = replace(trash_equations[i],j => ":")
        end
      end
    end
    number_eqns_with_no_jumps = sum(jumps_number .== 0)

    # Put the equations with no jumps in them at the top

    ind = sortperm(jumps_number)
    reordered_equations .= reordered_equations[ind]
    states_number .= states_number[ind]
    jumps_number .= jumps_number[ind]

    # Now sort out the order of the states in the system

    states_that_have_been_ordered = Int64[]
    for k = 1:length(states)
        for j = 1:number_eqns_with_no_jumps
            if states_number[j] == k
                for i = 1:length(states)
                    if occursin(reordered_states[i],reordered_equations[j]) == true  && j != i && (i in states_that_have_been_ordered) == false
                        reordered_states[i], reordered_states[j] = reordered_states[j], reordered_states[i]
                        push!(states_that_have_been_ordered,j)
                        break
                    end
                end
            end
        end
    end

    return reordered_equations, reordered_states

end

function reorder_equations_stochastic(equations::Array{Q,1},shocks::Array{Q,1},states::Array{Q,1},jumps::Array{Q,1},parameters::Array{Q,1}) where {Q <: AbstractString}

    #= This function reorders the model's equations so that the equations
       containing the shock processes appear first. =#

    reordered_equations = copy(equations)
    reordered_states = copy(states)
    reordered_shocks = copy(shocks)

    # Construct summary information about each equation

    trash_equations = copy(equations)

    variables = combine_variables(combine_variables(combine_variables(combine_variables(states,jumps),shocks),parameters),["exp","log"])

    variables = variables[sortperm(length.(variables),rev = true)] # sort variables so that longer names are first

    shocks_number = zeros(Int64,length(equations))
    states_number = zeros(Int64,length(equations))
    jumps_number  = zeros(Int64,length(equations))
    for i = 1:length(equations)
      for j in variables
        if j in shocks
          shocks_number[i] += occursin(j,trash_equations[i])
          trash_equations[i] = replace(trash_equations[i],j => ":")
        elseif j in states
          states_number[i] += occursin(j,trash_equations[i])
          trash_equations[i] = replace(trash_equations[i],j => ":")
        elseif j in jumps
          jumps_number[i] += occursin(j,equations[i])
          trash_equations[i] = replace(trash_equations[i],j => ":")
        else
          trash_equations[i] = replace(trash_equations[i],j => ":")
        end
      end
    end
    number_eqns_with_shocks = sum(shocks_number .!= 0)

    # Put the equations with no jumps in them at the top

    ind = sortperm(jumps_number)
    reordered_equations .= reordered_equations[ind]
    states_number .= states_number[ind]
    jumps_number .= jumps_number[ind]

    # Put the shock equations with the fewest shocks at the top

    ind = sortperm(shocks_number[1:number_eqns_with_shocks])
    reordered_equations[1:number_eqns_with_shocks] .= reordered_equations[ind]
    shocks_number[1:number_eqns_with_shocks] .= shocks_number[ind]
    states_number[1:number_eqns_with_shocks] .= states_number[ind]
    jumps_number[1:number_eqns_with_shocks]  .= jumps_number[ind]

    # Order the shock processes to try to make shock's variance-covariance matrix have non-zero diagonals.

    shocks_that_have_been_ordered = Int64[]
    for k = 1:length(shocks)
        for j = 1:number_eqns_with_shocks
            if shocks_number[j] == k
                for i = 1:length(shocks)
                    if occursin(reordered_shocks[i],reordered_equations[j]) == true  && j != i && (i in shocks_that_have_been_ordered) == false
                        reordered_shocks[i], reordered_shocks[j] = reordered_shocks[j], reordered_shocks[i]
                        push!(shocks_that_have_been_ordered,j)
                        break
                    end
                end
            end
        end
    end

    # Now sort out the order of the states in the system

    states_that_have_been_ordered = Int64[]
    for k = 1:length(states)
        for j = 1:number_eqns_with_shocks
            if states_number[j] == k
                for i = 1:length(states)
                    if occursin(reordered_states[i],reordered_equations[j]) == true  && j != i && (i in states_that_have_been_ordered) == false
                        reordered_states[i], reordered_states[j] = reordered_states[j], reordered_states[i]
                        push!(states_that_have_been_ordered,j)
                        break
                    end
                end
            end
        end
    end

    return reordered_equations, reordered_states, reordered_shocks

end

function reorganize_equations(equations::Array{Q,1},states::Array{Q,1},jumps::Array{Q,1},variables::Array{Q,1},lag_variables::Array{Q,1}) where {Q <: AbstractString}

    #= This function replaces lagged variables with pseudo current variables,
       augmenting the state vector and the model's equations accordingly. =#

    reorganized_equations = copy(equations)
    reorganized_states    = copy(states)

    # First we determine if any equation contains a lagged variable.

    model_has_lags = false
    for i = 1:length(equations)
        for j in lag_variables
            if occursin(j,equations[i]) == true
                model_has_lags = true
                break
            end
        end
    end

    #= If a model contains lagged variables, then we introduce a pseudo variable
       in its place and augment the list of state variables and the set of model
       equations. =#

    if model_has_lags == true
        for j = 1:length(lag_variables)
            if sum(occursin.(lag_variables[j],equations)) != 0
                for i = 1:length(equations)
                    reorganized_equations[i] = replace(reorganized_equations[i], lag_variables[j] => string(variables[j],"lag"))
                end
                reorganized_states    = [reorganized_states;string(variables[j],"lag")]
                reorganized_equations = [reorganized_equations;string(variables[j],"lag(+1) = ",variables[j])]
            end
        end
    end

    reorganized_variables = [reorganized_states;jumps]

    return reorganized_equations, reorganized_states, reorganized_variables

end

function get_re_model_primatives(model_array::Array{Q,1}) where {Q <: AbstractString}

    #= This function takes the model-array read from a model file, extracts the
       critical model information, does some basic error checking, and returns
       it in a structure. =#

    states                        = get_variables(model_array,"states:")
    jumps                         = get_variables(model_array,"jumps:")
    shocks                        = get_variables(model_array,"shocks:")
    variables                     = combine_variables(states,jumps)
    equations                     = get_equations(model_array,"equations:")
    (parameters, parametervalues, unassigned_parameters) = get_parameters_and_values(model_array,"parameters:")

    for i in [variables;parameters]
        if sum(occursin.(i,equations)) == false
            println("Warning: $i is not in any equation.")
        end
    end

    combined_names = [parameters;variables;shocks]
    if length(unique(combined_names)) != length(combined_names)
        error("Some parameters, variables, or shocks have the same name.")
    end

    reserved_names = ("exp","log","x","p",":",";")
    for name in reserved_names
        if name in combined_names
            error("$name cannot be the name for a variable, a shock, or a parameter.")
        end
    end

    lag_variables  = string.(variables,"(-1)")

    reordered_equations, states, shocks = reorder_equations(equations,shocks,states,jumps,parameters)
    reorganized_equations, states, variables = reorganize_equations(reordered_equations,states,jumps,variables,lag_variables)

    re_model_primatives = REModelPrimatives(states,jumps,shocks,variables,parameters,parametervalues,reorganized_equations,unassigned_parameters)

    return re_model_primatives

end

function repackage_equations(model::ModelPrimatives)

    #= This function is critical for repackaging the model's equations, replacing
       parameter names with values, and numbering variables. =#

    equations       = model.equations
    shocks          = model.shocks
    variables       = model.variables
    parameters      = model.parameters
    parametervalues = model.parametervalues

    repackaged_equations = copy(equations)

    if length(shocks) != 0
        combined_names = [variables;parameters;shocks;["exp","log"]]
    else
        combined_names = [variables;parameters;["exp","log"]]
    end

    sorted_combined_names = combined_names[sortperm(length.(combined_names),rev = true)]

    #= Now we go through every equation and replace future variables, variables, and
      shocks with a numbered element of a vector, "x".  We also replace parameter
      names with parameter values. =#

    for j in sorted_combined_names
        if j in variables
            variable_index = findfirst(isequal(j),variables)
            for i = 1:length(repackaged_equations)
                repackaged_equations[i] = replace(repackaged_equations[i],"$j(+1)" => "x[$(length(variables) + variable_index)]")
                repackaged_equations[i] = replace(repackaged_equations[i],j => "x[$(variable_index)]")
            end
        elseif j in parameters
            parameter_index = findfirst(isequal(j),parameters)
            for i = 1:length(repackaged_equations)
                repackaged_equations[i] = replace(repackaged_equations[i],j => parametervalues[parameter_index])
            end
        elseif j in shocks # Okay even if there are no shocks
            shock_index = findfirst(isequal(j),shocks)
            for i = 1:length(repackaged_equations)
                repackaged_equations[i] = replace(repackaged_equations[i],j => "x[$(2*length(variables) + shock_index)]")
            end
        elseif j == "exp"
            for i = 1:length(repackaged_equations)
                repackaged_equations[i] = replace(repackaged_equations[i],j => ":")
            end
        elseif j == "log"
            for i = 1:length(repackaged_equations)
                repackaged_equations[i] = replace(repackaged_equations[i],j => ";")
            end
        end
    end

    #= Now go back through every equation and restore exp and log where necessary =#

    for i = 1:length(repackaged_equations)
        if occursin(":",repackaged_equations[i]) == true
            repackaged_equations[i] = replace(repackaged_equations[i],":" => "exp")
        elseif occursin(";",repackaged_equations[i]) == true
            repackaged_equations[i] = replace(repackaged_equations[i],";" => "log")
        end
    end

    return repackaged_equations

end

function create_steady_state_equations(model::ModelPrimatives)

    # Make the model static by replacing leads and lags with current variables
    # and setting shocks equal to zero.  Also, replace parameter names with
    # their associated value.

    equations       = model.equations
    variables       = model.variables
    shocks          = model.shocks
    parameters      = model.parameters
    parametervalues = model.parametervalues

    steady_state_equations = copy(equations)

    if length(shocks) != 0
        combined_names = [variables;parameters;shocks]
    else
        combined_names = [variables;parameters]
    end

    sorted_combined_names = combined_names[sortperm(length.(combined_names),rev = true)]

    #= First we go through every equation and replace exp with : and log with ;.
       This is to guard them during variables and parameter substitution. =#

    for i = 1:length(steady_state_equations)
      if occursin("exp",steady_state_equations[i]) == true
        steady_state_equations[i] = replace(steady_state_equations[i],"exp" => ":")
      elseif occursin("log",steady_state_equations[i]) == true
        steady_state_equations[i] = replace(steady_state_equations[i],"log" => ";")
      end
    end

    # Now we go through every equation and replace future variables, variables, and
    # shocks with a numbered element of a vector, "x".  We also replace parameter
    # names with parameter values

    for j in sorted_combined_names
        if j in variables
            variable_index = findfirst(isequal(j),variables)
            for i = 1:length(equations)
                steady_state_equations[i] = replace(steady_state_equations[i],"$j(-1)" => "x[$(variable_index)]")
                steady_state_equations[i] = replace(steady_state_equations[i],"$j(+1)" => "x[$(variable_index)]")
                steady_state_equations[i] = replace(steady_state_equations[i],j => "x[$(variable_index)]")
            end
        elseif j in parameters
            parameter_index = findfirst(isequal(j),parameters)
            for i = 1:length(equations)
                steady_state_equations[i] = replace(steady_state_equations[i],j => parametervalues[parameter_index])
            end
        elseif j in shocks
            shock_index = findfirst(isequal(j),shocks)
            for i = 1:length(equations)
                steady_state_equations[i] = replace(steady_state_equations[i],j => 0.0)
            end
        end
    end

    #= Finally, go back through every equation and restore exp and log where necessary =#

    for i = 1:length(steady_state_equations)
        if occursin(":",steady_state_equations[i]) == true
            steady_state_equations[i] = replace(steady_state_equations[i],":" => "exp")
        elseif occursin(";",steady_state_equations[i]) == true
            steady_state_equations[i] = replace(steady_state_equations[i],";" => "log")
        end
    end

    return steady_state_equations

end

function make_equations_equal_zero(equations::Array{Q,1}) where {Q <: AbstractString}

    # Reexpress all of the model's equations such that they equal zero.

    zeroed_equations = similar(equations)

    for i = 1:length(equations)
        pair = strip.(split(equations[i], "="))
        zeroed_equations[i] = string(pair[1], " - (", pair[2], ")")
    end

    return zeroed_equations

end

function create_projection_equations(equations::Array{Q,1},model::ModelPrimatives) where {Q <: AbstractString}

    projection_equations = copy(equations)

    nx = length(model.states)
    ny = length(model.jumps)
    ns = length(model.shocks)
    ne = length(equations)
    nv = nx + ny

    for j = 1:nx
        for i = 1:ne
            projection_equations[i] = replace(projection_equations[i], "x[$j]" => "state[$j]")
        end
    end

    for j = 1:nv
        for i = 1:ne
            projection_equations[i] = replace(projection_equations[i], "x[$(nx+j)]" => "x[$j]")
        end
    end

    for j = 1:ny
        for i = 1:ne
            projection_equations[i] = replace(projection_equations[i], "x[$(nx+nv+j)]" => "approx$j")
        end
    end

    for j = 1:ns
        for i = 1:ne
            projection_equations[i] = replace(projection_equations[i], "x[$(2*nv+j)]" => 0.0)
        end
    end

    jumps_to_be_approximated = Int64[]
    eqns_to_be_approximated  = Int64[]
    for i = 1:ny
        for j = 1:ne
            if occursin("approx$i",projection_equations[j]) == true
                push!(jumps_to_be_approximated,i)
                push!(eqns_to_be_approximated,j)
            end
        end
    end

    jumps_to_be_approximated = unique(jumps_to_be_approximated)
    eqns_to_be_approximated  = sort(unique(eqns_to_be_approximated))

    return projection_equations, jumps_to_be_approximated, eqns_to_be_approximated

end

function create_processed_model_file(model::ModelPrimatives,path::Q) where {Q <: AbstractString}

    # Takes the model's primatives and turns these into a processed-model file.
    # This file is saved as a text file in the same folder as the model file.

    # First, get or construct all the information needed for the processed-model file

    repackaged_equations = repackage_equations(model)

    nonlinear_equations, jumps_to_be_approximated, eqns_to_be_approximated = create_projection_equations(repackaged_equations,model)
    projection_equations = make_equations_equal_zero(nonlinear_equations)

    steady_state_equations = create_steady_state_equations(model)
    static_equations       = make_equations_equal_zero(steady_state_equations)
    dynamic_equations      = make_equations_equal_zero(repackaged_equations)

    number_states    = length(model.states)
    number_jumps     = length(model.jumps)
    number_shocks    = length(model.shocks)
    number_variables = length(model.variables)
    number_equations = length(model.equations)

    variables = model.variables

    # Build up the string containing the processed model information that gets saved

    # First, add the model's summary information

    model_string = "nx = $number_states \n \n"
    model_string = string(model_string,"ny = $number_jumps \n \n")
    model_string = string(model_string,"ns = $number_shocks \n \n")
    model_string = string(model_string,"nv = $number_variables \n \n")
    model_string = string(model_string,"ne = $number_equations \n \n")

    model_string = string(model_string,"jumps_to_approximate = $jumps_to_be_approximated \n \n")
    model_string = string(model_string,"eqns_to_approximate = $eqns_to_be_approximated \n \n")
    model_string = string(model_string,"variables = $variables \n \n")

    # Second, add the model's static information

    if length(model.unassigned_parameters) != 0
        nlsolve_static_string = "function nlsolve_static_equations(f::Array{T,1},x::Array{T,1},p::Array{T1,1}) where {T<:Number, T1<:Real} \n \n"
        static_string         = "function static_equations(x::Array{T,1},p::Array{T1,1}) where {T<:Number, T1<:Real} \n \n"
    else
        nlsolve_static_string = "function nlsolve_static_equations(f::Array{T,1},x::Array{T,1}) where {T<:Number} \n \n"
        static_string         = "function static_equations(x::Array{T,1}) where {T<:Number} \n \n"
    end
    static_string         = string(static_string,"  f = Array{T,1}(undef,length(x)) \n \n")
    for i = 1:length(static_equations)
        static_string         = string(static_string,"  f[$i] = ",static_equations[i],"\n")
        nlsolve_static_string = string(nlsolve_static_string,"  f[$i] = ",static_equations[i],"\n")
    end

    static_string         = string(static_string,"\n  return f \n \n","end")
    nlsolve_static_string = string(nlsolve_static_string,"\n","end")

    model_string = string(model_string,nlsolve_static_string," \n \n",static_string," \n \n")

    # Third, add the model's dynamic information for perturbation solvers

    if length(model.unassigned_parameters) != 0
        dynamic_string = "function dynamic_equations(x::Array{T,1},p::Array{T1,1}) where {T<:Number, T1<:Real} \n \n"
    else
        dynamic_string = "function dynamic_equations(x::Array{T,1}) where {T<:Number} \n \n"
    end
    dynamic_string = string(dynamic_string,"  f = Array{T,1}(undef,$number_equations) \n \n")
    for i = 1:number_equations
        dynamic_string = string(dynamic_string,"  f[$i] = ",dynamic_equations[i],"\n")
    end

    dynamic_string = string(dynamic_string,"\n  return f \n \n","end \n")

    each_equation_string = Array{String}(undef,number_equations)

    for i = 1:number_equations
        if length(model.unassigned_parameters) != 0
            each_equation_string[i] = "function dynamic_eqn_$i(x::Array{T,1},p::Array{T1,1}) where {T<:Number, T1<:Real} \n \n"
        else
            each_equation_string[i] = "function dynamic_eqn_$i(x::Array{T,1}) where {T<:Number} \n \n"
        end
        each_equation_string[i] = string(each_equation_string[i],"  f = ", dynamic_equations[i],"\n","\n  return f \n \n","end \n")
    end

    for i = 1:number_equations
        dynamic_string = string(dynamic_string, "\n", each_equation_string[i])
    end

    individual_equations_string = "individual_equations = Array{Function}(undef,$number_equations) \n"
    for i = 1:number_equations
        individual_equations_string = string(individual_equations_string,"individual_equations[$i] = dynamic_eqn_$i","\n")
    end

    dynamic_string = string(dynamic_string,"\n",individual_equations_string)
    model_string = string(model_string,dynamic_string)

    # Fourth, add the model's dynamic information for projection solvers

    # For Chebyshev and Smolyak

    if length(model.unassigned_parameters) != 0
        closure_string = "function closure_projection_equations(state,scaled_weights,order,domain,approximate,p) \n \n"
    else
        closure_string = "function closure_projection_equations(state,scaled_weights,order,domain,approximate) \n \n"
    end
    closure_string = string(closure_string, "  function projection_equations(f::Array{T,1},x::Array{T,1}) where {T<:Number} \n \n")
    weight_number = 1
    for i in jumps_to_be_approximated
        closure_string = string(closure_string, "    approx$i = approximate(scaled_weights[$weight_number],x[$number_jumps+1:end],order,domain)","\n")
        weight_number += 1
    end
    closure_string = string(closure_string,"\n","    #f = Array{T,1}(undef,$number_equations) \n \n")
    for i = 1:length(projection_equations)
        closure_string = string(closure_string,"    f[$i] = ", projection_equations[i], "\n")
    end
    closure_string = string(closure_string,"\n    #return f \n \n  end \n \n  return projection_equations \n \n","end \n")

    # For piecewise linear

    if length(model.unassigned_parameters) != 0
        if number_shocks == 0  # We need to separate the function generated for the stochastic and deterministic cases
            closure_pl_string = "function closure_projection_equations_pl(variables,grid,state,approximate,p) \n \n"
        else
            closure_pl_string = "function closure_projection_equations_pl(variables,grid,state,integrals,approximate,p) \n \n"
        end
    else
        if number_shocks == 0  # We need to separate the function generated for the stochastic and deterministic cases
            closure_pl_string = "function closure_projection_equations_pl(variables,grid,state,approximate) \n \n"
        else
            closure_pl_string = "function closure_projection_equations_pl(variables,grid,state,integrals,approximate) \n \n"
        end
    end

    closure_pl_string = string(closure_pl_string,"  function projection_equations_pl(f::Array{T,1},x::Array{T,1}) where {T<:Number} \n \n")
    for i in jumps_to_be_approximated
        if number_shocks == 0
            closure_pl_string = string(closure_pl_string,"    approx$i = approximate(variables[$i],grid,x[$number_jumps+1:end])","\n")
        else
            closure_pl_string = string(closure_pl_string,"    approx$i = approximate(variables[$i],grid,x[$number_jumps+1:end],integrals)","\n")
        end
    end

    closure_pl_string = string(closure_pl_string,"\n","    #f = Array{T,1}(undef,$number_equations) \n \n")
    for i = 1:length(projection_equations)
        closure_pl_string = string(closure_pl_string,"    f[$i] = ",projection_equations[i], "\n")
    end

    closure_pl_string = string(closure_pl_string,"\n    #return f \n \n  end \n \n  return projection_equations_pl \n \n","end \n")

    model_string = string(model_string,"\n",closure_string)
    model_string = string(model_string,"\n",closure_pl_string)
    model_string = string(model_string,"\n","unassigned_parameters = $(model.unassigned_parameters)" )

    model_path = replace(path,".txt" => "_processed.txt")
    open(model_path, "w") do io write(io, model_string) end

end

function process_re_model(model_array::Array{Q,1},path::Q) where {Q <: AbstractString}

    # Creates the processed model structure for rational expectations models
    # (anticipating that other types of models may come later).

    re_model_primatives = get_re_model_primatives(model_array)

    create_processed_model_file(re_model_primatives,path)

    println("The model's variables are now in this order: ",re_model_primatives.variables)
    if length(re_model_primatives.unassigned_parameters) != 0
        println("The following parameters do not have values assigned: $(re_model_primatives.unassigned_parameters)")
    end

end

function process_model(path::Q) where {Q <: AbstractString}

    # Main function used to open, read, and process a model file.  The processed model
    # in written to a file that contains all the information needed for the model
    # solvers.

    model_array = open_model_file(path)

    process_re_model(model_array, path)

end

function retrieve_processed_model(path::Q) where {Q <: AbstractString}

    if !occursin("_processed",path)
        path = replace(path,".txt" => "_processed.txt")
    end

    include(path) # The information included is placed in the global scope

    if length(unassigned_parameters) != 0
        dsge_model = REModelPartial(nx,ny,ns,nv,ne,jumps_to_approximate,eqns_to_approximate,variables,nlsolve_static_equations,static_equations,dynamic_equations,individual_equations,closure_projection_equations,closure_projection_equations_pl,unassigned_parameters)
    else
        dsge_model = REModel(nx,ny,ns,nv,ne,jumps_to_approximate,eqns_to_approximate,variables,nlsolve_static_equations,static_equations,dynamic_equations,individual_equations,closure_projection_equations,closure_projection_equations_pl)
    end

    return dsge_model

end

function assign_parameters(model,param::Array{T,1}) where {T <: Number}

    nx = model.number_states
    ny = model.number_jumps
    ns = model.number_shocks
    nv = model.number_variables
    ne = model.number_equations
    jumps_approx = model.jumps_approximated
    eqns_approx = model.eqns_approximated
    vars = model.variables

    nlsse(f,x) = model.nlsolve_static_function(f,x,param)
    sf(x)      = model.static_function(x,param)
    df(x)      = model.dynamic_function(x,param)

    ief        = Array{Function}(undef,ne)
    for i = 1:ne
        ffie(x) = model.each_eqn_function[i](x,param)
        ief[i] = ffie
    end

    cf(state,scaled_weights,order,domain,approximate) = model.closure_function(state,scaled_weights,order,domain,approximate,param)
    cfpl_stoch(variables,grid,state,integrals,approximate)  = model.closure_function_piecewise(variables,grid,state,integrals,approximate,param)
    cfpl_det(variables,grid,state,approximate)  = model.closure_function_piecewise(variables,grid,state,approximate,param)

    if ns != 0
        newmod = REModel(nx,ny,ns,nv,ne,jumps_approx,eqns_approx,vars,nlsse,sf,df,ief,cf,cfpl_stoch)
        return newmod
    else
        newmod = REModel(nx,ny,ns,nv,ne,jumps_approx,eqns_approx,vars,nlsse,sf,df,ief,cf,cfpl_det)
        return newmod
    end

end