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

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*Problems*Involving*Inverse*Hyperbolic*Cosecants=#

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

    @test_int [x^3*acsch(a + b*x), x, 8, -(((2 - 17*a^2)*(a + b*x)*sqrt(1 + 1/(a + b*x)^2))/(12*b^4)) + (x^2*(a + b*x)*sqrt(1 + 1/(a + b*x)^2))/(12*b^2) - (a*(a + b*x)^2*sqrt(1 + 1/(a + b*x)^2))/(3*b^4) - (a^4*acsch(a + b*x))/(4*b^4) + (1/4)*x^4*acsch(a + b*x) + (a*(1 - 2*a^2)*atanh(sqrt(1 + 1/(a + b*x)^2)))/(2*b^4)]
    @test_int [x^2*acsch(a + b*x), x, 7, -((5*a*(a + b*x)*sqrt(1 + 1/(a + b*x)^2))/(6*b^3)) + (x*(a + b*x)*sqrt(1 + 1/(a + b*x)^2))/(6*b^2) + (a^3*acsch(a + b*x))/(3*b^3) + (1/3)*x^3*acsch(a + b*x) - ((1 - 6*a^2)*atanh(sqrt(1 + 1/(a + b*x)^2)))/(6*b^3)]
    @test_int [x^1*acsch(a + b*x), x, 6, ((a + b*x)*sqrt(1 + 1/(a + b*x)^2))/(2*b^2) - (a^2*acsch(a + b*x))/(2*b^2) + (1/2)*x^2*acsch(a + b*x) - (a*atanh(sqrt(1 + 1/(a + b*x)^2)))/b^2]
    @test_int [acsch(a + b*x)/x^1, x, 14, acsch(a + b*x)*log(1 - (a*ℯ^acsch(a + b*x))/(1 - sqrt(1 + a^2))) + acsch(a + b*x)*log(1 - (a*ℯ^acsch(a + b*x))/(1 + sqrt(1 + a^2))) - acsch(a + b*x)*log(1 - ℯ^(2*acsch(a + b*x))) + PolyLog(2, (a*ℯ^acsch(a + b*x))/(1 - sqrt(1 + a^2))) + PolyLog(2, (a*ℯ^acsch(a + b*x))/(1 + sqrt(1 + a^2))) - (1/2)*PolyLog(2, ℯ^(2*acsch(a + b*x)))]
    @test_int [acsch(a + b*x)/x^2, x, 6, -((b*acsch(a + b*x))/a) - acsch(a + b*x)/x + (2*b*atanh((a + tanh((1/2)*acsch(a + b*x)))/sqrt(1 + a^2)))/(a*sqrt(1 + a^2))]
    @test_int [acsch(a + b*x)/x^3, x, 8, (b*(a + b*x)*sqrt(1 + 1/(a + b*x)^2))/(2*a*(1 + a^2)*x) + (b^2*acsch(a + b*x))/(2*a^2) - acsch(a + b*x)/(2*x^2) - ((1 + 2*a^2)*b^2*atanh((a + tanh((1/2)*acsch(a + b*x)))/sqrt(1 + a^2)))/(a^2*(1 + a^2)^(3/2))]

    @test_int [(e + f*x)^3*(a + b*acsch(c + d*x))^2, x, 20, (b^2*f^2*(d*e - c*f)*x)/d^3 + (b^2*f^3*(c + d*x)^2)/(12*d^4) - (b*f^3*(c + d*x)*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/(3*d^4) + (3*b*f*(d*e - c*f)^2*(c + d*x)*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/d^4 + (b*f^2*(d*e - c*f)*(c + d*x)^2*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/d^4 + (b*f^3*(c + d*x)^3*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/(6*d^4) - ((d*e - c*f)^4*(a + b*acsch(c + d*x))^2)/(4*d^4*f) + ((e + f*x)^4*(a + b*acsch(c + d*x))^2)/(4*f) - (2*b*f^2*(d*e - c*f)*(a + b*acsch(c + d*x))*atanh(ℯ^acsch(c + d*x)))/d^4 + (4*b*(d*e - c*f)^3*(a + b*acsch(c + d*x))*atanh(ℯ^acsch(c + d*x)))/d^4 - (b^2*f^3*log(c + d*x))/(3*d^4) + (3*b^2*f*(d*e - c*f)^2*log(c + d*x))/d^4 - (b^2*f^2*(d*e - c*f)*PolyLog(2, -ℯ^acsch(c + d*x)))/d^4 + (2*b^2*(d*e - c*f)^3*PolyLog(2, -ℯ^acsch(c + d*x)))/d^4 + (b^2*f^2*(d*e - c*f)*PolyLog(2, ℯ^acsch(c + d*x)))/d^4 - (2*b^2*(d*e - c*f)^3*PolyLog(2, ℯ^acsch(c + d*x)))/d^4]
    @test_int [(e + f*x)^2*(a + b*acsch(c + d*x))^2, x, 17, (b^2*f^2*x)/(3*d^2) + (2*b*f*(d*e - c*f)*(c + d*x)*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/d^3 + (b*f^2*(c + d*x)^2*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/(3*d^3) - ((d*e - c*f)^3*(a + b*acsch(c + d*x))^2)/(3*d^3*f) + ((e + f*x)^3*(a + b*acsch(c + d*x))^2)/(3*f) - (2*b*f^2*(a + b*acsch(c + d*x))*atanh(ℯ^acsch(c + d*x)))/(3*d^3) + (4*b*(d*e - c*f)^2*(a + b*acsch(c + d*x))*atanh(ℯ^acsch(c + d*x)))/d^3 + (2*b^2*f*(d*e - c*f)*log(c + d*x))/d^3 - (b^2*f^2*PolyLog(2, -ℯ^acsch(c + d*x)))/(3*d^3) + (2*b^2*(d*e - c*f)^2*PolyLog(2, -ℯ^acsch(c + d*x)))/d^3 + (b^2*f^2*PolyLog(2, ℯ^acsch(c + d*x)))/(3*d^3) - (2*b^2*(d*e - c*f)^2*PolyLog(2, ℯ^acsch(c + d*x)))/d^3]
    @test_int [(e + f*x)^1*(a + b*acsch(c + d*x))^2, x, 11, (b*f*(c + d*x)*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/d^2 - ((d*e - c*f)^2*(a + b*acsch(c + d*x))^2)/(2*d^2*f) + ((e + f*x)^2*(a + b*acsch(c + d*x))^2)/(2*f) + (4*b*(d*e - c*f)*(a + b*acsch(c + d*x))*atanh(ℯ^acsch(c + d*x)))/d^2 + (b^2*f*log(c + d*x))/d^2 + (2*b^2*(d*e - c*f)*PolyLog(2, -ℯ^acsch(c + d*x)))/d^2 - (2*b^2*(d*e - c*f)*PolyLog(2, ℯ^acsch(c + d*x)))/d^2]
    @test_int [(e + f*x)^0*(a + b*acsch(c + d*x))^2, x, 8, ((c + d*x)*(a + b*acsch(c + d*x))^2)/d + (4*b*(a + b*acsch(c + d*x))*atanh(ℯ^acsch(c + d*x)))/d + (2*b^2*PolyLog(2, -ℯ^acsch(c + d*x)))/d - (2*b^2*PolyLog(2, ℯ^acsch(c + d*x)))/d]
    @test_int [(a + b*acsch(c + d*x))^2/(e + f*x)^1, x, 17, -(((a + b*acsch(c + d*x))^2*log(1 - ℯ^(2*acsch(c + d*x))))/f) + ((a + b*acsch(c + d*x))^2*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/f + ((a + b*acsch(c + d*x))^2*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/f - (b*(a + b*acsch(c + d*x))*PolyLog(2, ℯ^(2*acsch(c + d*x))))/f + (2*b*(a + b*acsch(c + d*x))*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/f + (2*b*(a + b*acsch(c + d*x))*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/f + (b^2*PolyLog(3, ℯ^(2*acsch(c + d*x))))/(2*f) - (2*b^2*PolyLog(3, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/f - (2*b^2*PolyLog(3, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/f]
    @test_int [(a + b*acsch(c + d*x))^2/(e + f*x)^2, x, 12, (d*(a + b*acsch(c + d*x))^2)/(f*(d*e - c*f)) - (a + b*acsch(c + d*x))^2/(f*(e + f*x)) - (2*b*d*(a + b*acsch(c + d*x))*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/((d*e - c*f)*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)) + (2*b*d*(a + b*acsch(c + d*x))*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/((d*e - c*f)*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)) - (2*b^2*d*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/((d*e - c*f)*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)) + (2*b^2*d*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/((d*e - c*f)*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))]
    @test_int [(a + b*acsch(c + d*x))^2/(e + f*x)^3, x, 23, -((b*d^2*f*sqrt(1 + 1/(c + d*x)^2)*(a + b*acsch(c + d*x)))/((d*e - c*f)*(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)*(f + (d*e - c*f)/(c + d*x)))) + (d^2*(a + b*acsch(c + d*x))^2)/(2*f*(d*e - c*f)^2) - (a + b*acsch(c + d*x))^2/(2*f*(e + f*x)^2) + (b*d^2*f^2*(a + b*acsch(c + d*x))*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/((d*e - c*f)^2*(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)^(3/2)) - (2*b*d^2*(a + b*acsch(c + d*x))*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/((d*e - c*f)^2*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)) - (b*d^2*f^2*(a + b*acsch(c + d*x))*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/((d*e - c*f)^2*(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)^(3/2)) + (2*b*d^2*(a + b*acsch(c + d*x))*log(1 + (ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))))/((d*e - c*f)^2*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)) + (b^2*d^2*f*log(f + (d*e - c*f)/(c + d*x)))/((d*e - c*f)^2*(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)) + (b^2*d^2*f^2*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/((d*e - c*f)^2*(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)^(3/2)) - (2*b^2*d^2*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f - sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/((d*e - c*f)^2*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)) - (b^2*d^2*f^2*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/((d*e - c*f)^2*(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)^(3/2)) + (2*b^2*d^2*PolyLog(2, -((ℯ^acsch(c + d*x)*(d*e - c*f))/(f + sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2)))))/((d*e - c*f)^2*sqrt(d^2*e^2 - 2*c*d*e*f + (1 + c^2)*f^2))]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*acsch(a*x^n)=#

    @test_int [x^3*acsch(sqrt(x)), x, 4, -((sqrt(-1 - x)*sqrt(x))/(4*sqrt(-x))) - ((-1 - x)^(3/2)*sqrt(x))/(4*sqrt(-x)) - (3*(-1 - x)^(5/2)*sqrt(x))/(20*sqrt(-x)) - ((-1 - x)^(7/2)*sqrt(x))/(28*sqrt(-x)) + (1/4)*x^4*acsch(sqrt(x))]
    @test_int [x^2*acsch(sqrt(x)), x, 4, (sqrt(-1 - x)*sqrt(x))/(3*sqrt(-x)) + (2*(-1 - x)^(3/2)*sqrt(x))/(9*sqrt(-x)) + ((-1 - x)^(5/2)*sqrt(x))/(15*sqrt(-x)) + (1/3)*x^3*acsch(sqrt(x))]
    @test_int [x^1*acsch(sqrt(x)), x, 4, -((sqrt(-1 - x)*sqrt(x))/(2*sqrt(-x))) - ((-1 - x)^(3/2)*sqrt(x))/(6*sqrt(-x)) + (1/2)*x^2*acsch(sqrt(x))]
    @test_int [x^0*acsch(sqrt(x)), x, 3, (sqrt(-1 - x)*sqrt(x))/sqrt(-x) + x*acsch(sqrt(x))]
    @test_int [acsch(sqrt(x))/x^1, x, 7, acsch(sqrt(x))^2 - 2*acsch(sqrt(x))*log(1 - ℯ^(2*acsch(sqrt(x)))) - PolyLog(2, ℯ^(2*acsch(sqrt(x))))]
    @test_int [acsch(sqrt(x))/x^2, x, 5, sqrt(-1 - x)/(2*sqrt(-x)*sqrt(x)) - acsch(sqrt(x))/x - (sqrt(x)*atan(sqrt(-1 - x)))/(2*sqrt(-x))]
    @test_int [acsch(sqrt(x))/x^3, x, 6, sqrt(-1 - x)/(8*sqrt(-x)*x^(3/2)) - (3*sqrt(-1 - x))/(16*sqrt(-x)*sqrt(x)) - acsch(sqrt(x))/(2*x^2) + (3*sqrt(x)*atan(sqrt(-1 - x)))/(16*sqrt(-x))]
    @test_int [acsch(sqrt(x))/x^4, x, 7, sqrt(-1 - x)/(18*sqrt(-x)*x^(5/2)) - (5*sqrt(-1 - x))/(72*sqrt(-x)*x^(3/2)) + (5*sqrt(-1 - x))/(48*sqrt(-x)*sqrt(x)) - acsch(sqrt(x))/(3*x^3) - (5*sqrt(x)*atan(sqrt(-1 - x)))/(48*sqrt(-x))]

    @test_int [acsch(1/x), x, 3, -sqrt(1 + x^2) + x*asinh(x)]

    @test_int [acsch(a*x^n)/x, x, 7, acsch(a*x^n)^2/(2*n) - (acsch(a*x^n)*log(1 - ℯ^(2*acsch(a*x^n))))/n - PolyLog(2, ℯ^(2*acsch(a*x^n)))/(2*n)]
    @test_int [acsch(a*x^5)/x, x, 7, (1/10)*acsch(a*x^5)^2 - (1/5)*acsch(a*x^5)*log(1 - ℯ^(2*acsch(a*x^5))) - (1/10)*PolyLog(2, ℯ^(2*acsch(a*x^5)))]

    #= ::Section::Closed:: =#
    #=Integrands*involving*inverse*hyperbolic*cosecants*of*exponentials=#

    @test_int [acsch(c*ℯ^(a + b*x)), x, 7, acsch(c*ℯ^(a + b*x))^2/(2*b) - (acsch(c*ℯ^(a + b*x))*log(1 - ℯ^(2*acsch(c*ℯ^(a + b*x)))))/b - PolyLog(2, ℯ^(2*acsch(c*ℯ^(a + b*x))))/(2*b)]

    #= ::Section::Closed:: =#
    #=Integrands*involving*exponentials*of*inverse*hyperbolic*cosecants=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*ℯ^acsch(a*x^p)=#

    @test_int [x^m*ℯ^acsch(a*x), x, 4, x^m/(a*m) + (x^(1 + m)*HypergeometricFunctions._₂F₁(-(1/2), (1/2)*(-1 - m), (1 - m)/2, -(1/(a^2*x^2))))/(1 + m)]

    @test_int [x^4*ℯ^acsch(a*x), x, 4, -((2*(1 + 1/(a^2*x^2))^(3/2)*x^3)/(15*a^2)) + x^4/(4*a) + (1/5)*(1 + 1/(a^2*x^2))^(3/2)*x^5]
    @test_int [x^3*ℯ^acsch(a*x), x, 7, (sqrt(1 + 1/(a^2*x^2))*x^2)/(8*a^2) + x^3/(3*a) + (1/4)*sqrt(1 + 1/(a^2*x^2))*x^4 - atanh(sqrt(1 + 1/(a^2*x^2)))/(8*a^4)]
    @test_int [x^2*ℯ^acsch(a*x), x, 3, x^2/(2*a) + (1/3)*(1 + 1/(a^2*x^2))^(3/2)*x^3]
    @test_int [x^1*ℯ^acsch(a*x), x, 6, x/a + (1/2)*sqrt(1 + 1/(a^2*x^2))*x^2 + atanh(sqrt(1 + 1/(a^2*x^2)))/(2*a^2)]
    @test_int [x^0*ℯ^acsch(a*x), x, 5, ℯ^acsch(a*x)*x - acsch(a*x)/a + log(x)/a, sqrt(1 + 1/(a^2*x^2))*x - acsch(a*x)/a + log(x)/a]
    @test_int [ℯ^acsch(a*x)/x^1, x, 6, -sqrt(1 + 1/(a^2*x^2)) - 1/(a*x) + atanh(sqrt(1 + 1/(a^2*x^2)))]
    @test_int [ℯ^acsch(a*x)/x^2, x, 5, -(1/(2*a*x^2)) - sqrt(1 + 1/(a^2*x^2))/(2*x) - (1/2)*a*acsch(a*x)]
    @test_int [ℯ^acsch(a*x)/x^3, x, 3, (-(1/3))*a^2*(1 + 1/(a^2*x^2))^(3/2) - 1/(3*a*x^3)]
    @test_int [ℯ^acsch(a*x)/x^4, x, 6, -(1/(4*a*x^4)) - sqrt(1 + 1/(a^2*x^2))/(4*x^3) - (a^2*sqrt(1 + 1/(a^2*x^2)))/(8*x) + (1/8)*a^3*acsch(a*x)]
    @test_int [ℯ^acsch(a*x)/x^5, x, 5, (1/3)*a^4*(1 + 1/(a^2*x^2))^(3/2) - (1/5)*a^4*(1 + 1/(a^2*x^2))^(5/2) - 1/(5*a*x^5)]

    @test_int [x^m*ℯ^acsch(a*x^2), x, 4, -(x^(-1 + m)/(a*(1 - m))) + (x^(1 + m)*HypergeometricFunctions._₂F₁(-(1/2), (1/4)*(-1 - m), (3 - m)/4, -(1/(a^2*x^4))))/(1 + m)]

    @test_int [x^4*ℯ^acsch(a*x^2), x, 8, -((2*sqrt(1 + 1/(a^2*x^4)))/(5*a^2*(a + 1/x^2)*x)) + (2*sqrt(1 + 1/(a^2*x^4))*x)/(5*a^2) + x^3/(3*a) + (1/5)*sqrt(1 + 1/(a^2*x^4))*x^5 + (2*sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.E(2*acot(sqrt(a)*x), 1/2))/(5*a^(7/2)*sqrt(1 + 1/(a^2*x^4))) - (sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.F(2*acot(sqrt(a)*x), 1/2))/(5*a^(7/2)*sqrt(1 + 1/(a^2*x^4)))]
    @test_int [x^3*ℯ^acsch(a*x^2), x, 6, x^2/(2*a) + (1/4)*sqrt(1 + 1/(a^2*x^4))*x^4 + atanh(sqrt(1 + 1/(a^2*x^4)))/(4*a^2)]
    @test_int [x^2*ℯ^acsch(a*x^2), x, 5, x/a + (1/3)*sqrt(1 + 1/(a^2*x^4))*x^3 - (sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.F(2*acot(sqrt(a)*x), 1/2))/(3*a^(5/2)*sqrt(1 + 1/(a^2*x^4)))]
    @test_int [x^1*ℯ^acsch(a*x^2), x, 6, (1/2)*sqrt(1 + 1/(a^2*x^4))*x^2 - acsch(a*x^2)/(2*a) + log(x)/a]
    @test_int [x^0*ℯ^acsch(a*x^2), x, 7, -(1/(a*x)) - (2*sqrt(1 + 1/(a^2*x^4)))/((a + 1/x^2)*x) + sqrt(1 + 1/(a^2*x^4))*x + (2*sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.E(2*acot(sqrt(a)*x), 1/2))/(a^(3/2)*sqrt(1 + 1/(a^2*x^4))) - (sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.F(2*acot(sqrt(a)*x), 1/2))/(a^(3/2)*sqrt(1 + 1/(a^2*x^4)))]
    @test_int [ℯ^acsch(a*x^2)/x^1, x, 6, (-(1/2))*sqrt(1 + 1/(a^2*x^4)) - 1/(2*a*x^2) + (1/2)*atanh(sqrt(1 + 1/(a^2*x^4)))]
    @test_int [ℯ^acsch(a*x^2)/x^2, x, 5, -(1/(3*a*x^3)) - sqrt(1 + 1/(a^2*x^4))/(3*x) - (sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.F(2*acot(sqrt(a)*x), 1/2))/(3*sqrt(a)*sqrt(1 + 1/(a^2*x^4)))]
    @test_int [ℯ^acsch(a*x^2)/x^3, x, 6, -(1/(4*a*x^4)) - sqrt(1 + 1/(a^2*x^4))/(4*x^2) - (1/4)*a*acsch(a*x^2)]
    @test_int [ℯ^acsch(a*x^2)/x^4, x, 7, -(1/(5*a*x^5)) - sqrt(1 + 1/(a^2*x^4))/(5*x^3) - (2*a^2*sqrt(1 + 1/(a^2*x^4)))/(5*(a + 1/x^2)*x) + (2*sqrt(a)*sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.E(2*acot(sqrt(a)*x), 1/2))/(5*sqrt(1 + 1/(a^2*x^4))) - (sqrt(a)*sqrt((a^2 + 1/x^4)/(a + 1/x^2)^2)*(a + 1/x^2)*Elliptic.F(2*acot(sqrt(a)*x), 1/2))/(5*sqrt(1 + 1/(a^2*x^4)))]
    @test_int [ℯ^acsch(a*x^2)/x^5, x, 3, (-(1/6))*a^2*(1 + 1/(a^2*x^4))^(3/2) - 1/(6*a*x^6)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*ℯ^(n*acsch(a*x))=#

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

    @test_int [x^m*ℯ^(2*acsch(a*x)), x, 5, -((2*x^(-1 + m))/(a^2*(1 - m))) + x^(1 + m)/(1 + m) + (2*x^m*HypergeometricFunctions._₂F₁(-(1/2), -(m/2), 1 - m/2, -(1/(a^2*x^2))))/(a*m)]

    @test_int [x^4*ℯ^(2*acsch(a*x)), x, 8, (sqrt(1 + 1/(a^2*x^2))*x^2)/(4*a^3) + (2*x^3)/(3*a^2) + (sqrt(1 + 1/(a^2*x^2))*x^4)/(2*a) + x^5/5 - atanh(sqrt(1 + 1/(a^2*x^2)))/(4*a^5)]
    @test_int [x^3*ℯ^(2*acsch(a*x)), x, 4, x^2/a^2 + (2*(1 + 1/(a^2*x^2))^(3/2)*x^3)/(3*a) + x^4/4]
    @test_int [x^2*ℯ^(2*acsch(a*x)), x, 7, (2*x)/a^2 + (sqrt(1 + 1/(a^2*x^2))*x^2)/a + x^3/3 + atanh(sqrt(1 + 1/(a^2*x^2)))/a^3]
    @test_int [x^1*ℯ^(2*acsch(a*x)), x, 6, (2*sqrt(1 + 1/(a^2*x^2))*x)/a + x^2/2 - (2*acsch(a*x))/a^2 + (2*log(x))/a^2]
    @test_int [x^0*ℯ^(2*acsch(a*x)), x, 7, -((2*sqrt(1 + 1/(a^2*x^2)))/a) - 2/(a^2*x) + x + (2*atanh(sqrt(1 + 1/(a^2*x^2))))/a]
    @test_int [ℯ^(2*acsch(a*x))/x^1, x, 6, -(1/(a^2*x^2)) - sqrt(1 + 1/(a^2*x^2))/(a*x) - acsch(a*x) + log(x)]
    @test_int [ℯ^(2*acsch(a*x))/x^2, x, 4, (-(2/3))*a*(1 + 1/(a^2*x^2))^(3/2) - 2/(3*a^2*x^3) - 1/x, (-(1/2))*a*sqrt(1 + 1/(a^2*x^2)) - (1/6)*a*(sqrt(1 + 1/(a^2*x^2)) + 1/(a*x))^3 - 1/(2*x)]
    @test_int [ℯ^(2*acsch(a*x))/x^3, x, 7, -(1/(2*a^2*x^4)) - sqrt(1 + 1/(a^2*x^2))/(2*a*x^3) - 1/(2*x^2) - (a*sqrt(1 + 1/(a^2*x^2)))/(4*x) + (1/4)*a^2*acsch(a*x)]
    @test_int [ℯ^(2*acsch(a*x))/x^4, x, 6, (2/3)*a^3*(1 + 1/(a^2*x^2))^(3/2) - (2/5)*a^3*(1 + 1/(a^2*x^2))^(5/2) - 2/(5*a^2*x^5) - 1/(3*x^3)]
    @test_int [ℯ^(2*acsch(a*x))/x^5, x, 8, -(1/(3*a^2*x^6)) - sqrt(1 + 1/(a^2*x^2))/(3*a*x^5) - 1/(4*x^4) - (a*sqrt(1 + 1/(a^2*x^2)))/(12*x^3) + (a^3*sqrt(1 + 1/(a^2*x^2)))/(8*x) - (1/8)*a^4*acsch(a*x)]

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

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*ℯ^(n*acsch(a*x)) / (1-a^2*x^2)=#

    @test_int [(d*x)^m*ℯ^acsch(c*x)/(1 + c^2*x^2), x, 4, -((d*(d*x)^(-1 + m)*HypergeometricFunctions._₂F₁(1/2, (1 - m)/2, (3 - m)/2, -(1/(c^2*x^2))))/(c^2*(1 - m))) + ((d*x)^m*HypergeometricFunctions._₂F₁(1, m/2, (2 + m)/2, (-c^2)*x^2))/(c*m)]

    @test_int [x^5*ℯ^acsch(c*x)/(1 + c^2*x^2), x, 9, -(x/c^5) - (3*sqrt(1 + 1/(c^2*x^2))*x^2)/(8*c^4) + x^3/(3*c^3) + (sqrt(1 + 1/(c^2*x^2))*x^4)/(4*c^2) + atan(c*x)/c^6 + (3*atanh(sqrt(1 + 1/(c^2*x^2))))/(8*c^6)]
    @test_int [x^4*ℯ^acsch(c*x)/(1 + c^2*x^2), x, 6, -((2*sqrt(1 + 1/(c^2*x^2))*x)/(3*c^4)) + x^2/(2*c^3) + (sqrt(1 + 1/(c^2*x^2))*x^3)/(3*c^2) - log(1 + c^2*x^2)/(2*c^5)]
    @test_int [x^3*ℯ^acsch(c*x)/(1 + c^2*x^2), x, 7, x/c^3 + (sqrt(1 + 1/(c^2*x^2))*x^2)/(2*c^2) - atan(c*x)/c^4 - atanh(sqrt(1 + 1/(c^2*x^2)))/(2*c^4)]
    @test_int [x^2*ℯ^acsch(c*x)/(1 + c^2*x^2), x, 3, (sqrt(1 + 1/(c^2*x^2))*x)/c^2 + log(1 + c^2*x^2)/(2*c^3)]
    @test_int [x^1*ℯ^acsch(c*x)/(1 + c^2*x^2), x, 5, atan(c*x)/c^2 + atanh(sqrt(1 + 1/(c^2*x^2)))/c^2]
    @test_int [x^0*ℯ^acsch(c*x)/(1 + c^2*x^2), x, 7, -(acsch(c*x)/c) + log(x)/c - log(1 + c^2*x^2)/(2*c)]
    @test_int [ℯ^acsch(c*x)/(x^1*(1 + c^2*x^2)), x, 4, -sqrt(1 + 1/(c^2*x^2)) - 1/(c*x) - atan(c*x)]
    @test_int [ℯ^acsch(c*x)/(x^2*(1 + c^2*x^2)), x, 7, -(1/(2*c*x^2)) - sqrt(1 + 1/(c^2*x^2))/(2*x) + (1/2)*c*acsch(c*x) - c*log(x) + (1/2)*c*log(1 + c^2*x^2)]
    @test_int [ℯ^acsch(c*x)/(x^3*(1 + c^2*x^2)), x, 7, c^2*sqrt(1 + 1/(c^2*x^2)) - (1/3)*c^2*(1 + 1/(c^2*x^2))^(3/2) - 1/(3*c*x^3) + c/x + c^2*atan(c*x)]

    #= ::Section::Closed:: =#
    #=Miscellaneous*integrands*involving*inverse*hyperbolic*cosecants=#

    @test_int [acsch(a + b*x)/((a*d)/b + d*x), x, 8, acsch(a + b*x)^2/(2*d) - (acsch(a + b*x)*log(1 - ℯ^(2*acsch(a + b*x))))/d - PolyLog(2, ℯ^(2*acsch(a + b*x)))/(2*d)]

    @test_int [x^3*acsch(a + b*x^4), x, 6, ((a + b*x^4)*acsch(a + b*x^4))/(4*b) + atanh(sqrt(1 + 1/(a + b*x^4)^2))/(4*b)]

    @test_int [x^(n-1)*acsch(a + b*x^n), x, 6, ((a + b*x^n)*acsch(a + b*x^n))/(b*n) + atanh(sqrt(1 + 1/(a + b*x^n)^2))/(b*n)]
    end

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.6 Inverse hyperbolic cosecant/7.6.2 Inverse hyperbolic cosecant functions.jl


text:
function sigmoid(z)
  1 ./ (1 .+ e .^ (-z))
end

function sigmoid_grad(z)
  sigmoid(z) .* (1 - sigmoid(z))
end

function init_theta(num_input, num_output, nums_hidden)
  num_units = [num_input, nums_hidden..., num_output]
  layers = length(num_units)
  Theta = [rand(num_units[i], num_units[i-1]+1) - 0.5 for i = 2:layers]
end

function cost(Y, Y_pred)
  len = size(Y, 1)
 J = (-1.0 / len ) * sum((Y .* log(Y_pred)) + ((1-Y) .* log(1-Y_pred)))
end

function forward_prop(X, Theta)
  A = {}
  H = {}
  push!(H, X)
  push!(A, zeros(1, 1))

  for t=1:length(Theta)
    # add bias
    H[t] = [ones(size(H[t], 1), 1) H[t]]
    # pre-activation
    push!(A, H[t] * Theta[t]')
    # activation
    push!(H, sigmoid(A[t+1]))
  end
  A, H, H[end]
end

function back_prop(X, Y, Theta)
  len = size(X, 1)
  T = length(Theta)

  Theta_grad = {}
  for i=1:T
    push!(Theta_grad, zeros(size(Theta[i])))
  end

  A, H, Y_pred = forward_prop(X, Theta)

  delta_N = {}

  for t=1:T
    push!(delta_N, {})
  end

  delta = Y_pred - Y

  delta_N[T] = delta

  for t=0:T-2
    delta = (delta*Theta[T-t][:, 2:end]) .* sigmoid_grad(A[T-t])
    delta_N[T-t-1] = delta
  end

  for t=1:T
    Theta_grad[t] = Theta_grad[t] + delta_N[t]' * H[t]
  end

  Y_pred, Theta_grad
end

function fit(X, Y, Theta, learning_rate, epochs)
  J_list = zeros(epochs)
  for i=1:epochs
    Y_pred, Theta_grad = back_prop(X, Y, Theta)
    J_list[i] = cost(Y_train, Y_pred)
    for t=1:length(Theta)
      Theta[t] = Theta[t] - learning_rate*Theta_grad[t]
    end
  end
  Theta, J_list
end

function predict(X, Theta)
  _, _, y = forward_prop(X, Theta)
  y
end

ext: jl


text:
function check_valid_sat_pure(model,P_sat,V_l,V_v,T)
    _,dpdvl = p∂p∂V(model,V_l,T,SA[1.0])
    _,dpdvv = p∂p∂V(model,V_v,T,SA[1.0])
    (dpdvl > 0) | (dpdvv > 0) && return false
    ε = abs(V_l-V_v)/(eps(typeof(V_l-V_v)))
    #if ΔV > ε then Vl and Vv are different values
    return ε > 5e9
end

function try_sat_pure(model,V0,f!,T,result,error_val,method = LineSearch(Newton()))
    if !isfinite(V0[1]) | !isfinite(V0[2])
        return false
    end
    try
        res = sat_pure(model,V0,f!,T,method)
        result[] = res
    catch e #normally, failures occur near the critical point
        error_val[] = e
        return false
    end

    (P_sat,V_l,V_v) = result[]
    return check_valid_sat_pure(model,P_sat,V_l,V_v,T)
end

"""
    saturation_pressure(model::EoSModel, T, V0 = x0_sat_pure(model,T))

Performs a single component saturation equilibrium calculation, at the specified temperature `T`, of one mol of pure sustance specified by `model`

Returns `(p₀, Vₗ, Vᵥ)` where `p₀` is the saturation pressure (in Pa), `Vₗ` is the liquid saturation volume (in m³) and `Vᵥ` is the vapour saturation volume (in m³).

If the calculation fails, returns  `(NaN, NaN, NaN)`

`V0` is `[log10(Vₗ₀),log10(Vᵥ₀)]` , where `Vₗ₀`  and `Vᵥ₀` are initial guesses for the liquid and vapour volumes.
"""
function saturation_pressure(model::EoSModel, T, V0 = x0_sat_pure(model,T))
    T = T*T/T
    V_lb = lb_volume(model,SA[1.0])
    TYPE = promote_type(typeof(T),typeof(V_lb))
    nan = zero(TYPE)/zero(TYPE)
    scales = scale_sat_pure(model)
    f! = (F,x) -> Obj_Sat(model, F, T, exp10(x[1]), exp10(x[2]),scales)
    res0 = (nan,nan,nan)
    result = Ref(res0)
    error_val = Ref{Any}(nothing)
    converged = try_sat_pure(model,V0,f!,T,result,error_val)
    #did not converge, but didnt error.
    if converged
        return result[]
    end
    (T_c, p_c, V_c) = crit_pure(model)
    if abs(T_c-T) < eps(typeof(T))
        return (p_c,V_c,V_c)
    end
    if T_c < T
        @error "initial temperature $T greater than critical temperature $T_c. returning NaN"
    else
        V0 = x0_sat_pure_crit(model,T,T_c,p_c,V_c)
        converged = try_sat_pure(model,V0,f!,T,result,error_val)
        if converged
            return result[]
        end
    end
    #not converged, even trying with better critical aprox.
    return res0
end

#with the critical point, we can perform a
#corresponding states approximation with the
#propane reference equation of state
function x0_sat_pure_crit(model,T,T_c,P_c,V_c)
    h = V_c*5000
    T0 = 369.89*T/T_c
    Vl0 = (1.0/_propaneref_rholsat(T0))*h
    Vv0 = (1.0/_propaneref_rhovsat(T0))*h
    _1 = SA[1.0]
    #μ_l = only(VT_chemical_potential(model,Vl0,T,_1))
    #μ_v = only(VT_chemical_potential(model,Vv0,T,_1))
    #@show (μ_l < μ_v,T/T_c)
    #if μ_l < μ_v
      #@show μ_l,μ_v
    #end
    # _,dpdvv = p∂p∂V(model,Vv0,T,SA[1.0])
    # @show dpdvv*Vv0
    # _,dpdvv = p∂p∂V(model,2*Vv0,T,SA[1.0])
    # @show dpdvv*Vv0
    return [log10(Vl0),log10(Vv0)]
end
function sat_pure(model::EoSModel,V0,f!,T,method =LineSearch(Newton()))
    r = Solvers.nlsolve(f!, V0 ,method )
    Vsol = Solvers.x_sol(r)
    V_l = exp10(Vsol[1])
    V_v = exp10(Vsol[2])
    P_sat = pressure(model,V_v,T)
    return (P_sat,V_l,V_v)
end

function Obj_Sat(model::EoSModel, F, T, V_l, V_v,scales)
    fun(x) = eos(model, x[2], T,SA[x[1]])
    df(x)  = ForwardDiff.gradient(fun,x)
    df_l = df(SA[one(V_l*T),V_l*one(T)])
    df_v = df(SA[one(V_v),V_v*one(T)])
    (p_scale,μ_scale) = scales
    #T̄ = T/T_scale(model)
    F[1] = (df_l[2]-df_v[2])*p_scale
    F[2] = (df_l[1]-df_v[1])*μ_scale
    return F
end

function saturation_temperature(model::EoSModel,p)
    f(z) = Obj_sat_pure_T(model,z,p)
    Tc,pc,vc = crit_pure(model)
    if p>pc
        return (NaN,NaN,NaN)
    else
        T = Roots.find_zero(f,(0.3*Tc,0.99*Tc))
        p,v_l,v_v = saturation_pressure(model,T)
        return T,v_l,v_v
    end
end

function Obj_sat_pure_T(model,T,p)
    p̃,v_l,v_v = saturation_pressure(model,T)
    return p̃-p
end

"""
    enthalpy_vap(model::EoSModel, T)

Calculates `ΔH`, the difference between saturated vapour and liquid enthalpies at temperature `T`, in J
"""
function enthalpy_vap(model::EoSModel, T)
    (P_sat,V_l,V_v) = saturation_pressure(model,T)
    H_v = VT_enthalpy(model,V_v,T)
    H_l = VT_enthalpy(model,V_l,T)
    H_vap=H_v -H_l
    return H_vap
end

function choose_p(p1,p2)
    if p1 > 0 && p2 > 0
        p = (p1+p2)/2
    elseif p1 < 0
        p = p2
    elseif p2 < 0
        p = p1
    elseif isnan(p1) || isnan(p2)
        p = zero(p1)/zero(p2)
    end
    return p
end
#=
    an alternative algorithm for saturation pressure.

"""
    FugSucessiveSatPure(;max_iters=100,ptol= 1e-10,vtol=1e-15,V0=nothing)

Performs a saturation pressure calculation, by sucessive substitution:
```julia
    vl = pressure(model,p,T,phase=:l)
    vv = pressure(model,p,T,phase=:v)
    p = (eos(model,vl,T) - eos(model,vv,T))/(vv-vl)
```
It's more stable than the default volume based method, but it is slower.
if no initial `V0` is passed, it will be calculated via `x0_sat_pure(model,T)`.if a number is passed, it will be used as an initial pressure.
"""
Base.@kwdef struct FugSucessiveSatPure{T}
    max_iters::Int = 100
    ptol::Float64 = 1e-10
    vtol::Float64 = 1e-15
    V0::T = nothing
end

export FugSucessiveSatPure

function saturation_pressure(model::EoSModel,T,method::FugSucessiveSatPure)
    nan = zero(T)/zero(T)
    if isnothing(method.V0) || length(method.V0) == 2
        if isnothing(method.V0)
            V1,V2 = x0_sat_pure(model,T)
        else
            V1,V2 = method.V0
        end
        V1,V2 = exp10(V1),exp10(V2)
        p1,p2 = pressure(model,V1,T), pressure(model,V2,T)
    elseif (method.V0) isa Real
        p = method.V0
        p1,p2 = p
        V1 = volume(model,p,T;phase=:l)
        V2 = volume(model,p,T;phase=:v)
    else
        throw(error("invalid initial point for saturation pressure"))
    end
    T = T*one(nan)
    A(V, t) = eos(model, V, t,SA[1.0])
    p = choose_p(p1,p2)
    if isnan(p)
        T_c,P_c,V_c = crit_pure(model)
        if T_c < T
            @error "initial temperature $T greater than critical temperature $T_c. returning NaN"
        end
        V1,V2 = x0_sat_pure_crit(model,T,T_c,P_c,V_c)
        V1,V2 = exp10(V1),exp10(V2)
        p1,p2 = pressure(model,V1,T), pressure(model,V2,T)
        p = choose_p(p1,p2)
    end
    V1old = nan
    V2old = nan
    if isnan(p)
        @error("cannot find initial pressure")
        return (nan,nan,nan)
    end
    pold = one(T)/zero(T)
    iters = method.max_iters
    ptol = method.ptol
    vtol = method.vtol

    function Fug(model,V,T,p)
        μres = ForwardDiff.derivative(z->eos_res(model,V,T,SA[z]),1.0)
        Z = p*V/(R̄*T)
        exp(μres/R̄/T)/Z
    end
    for i = 1:iters
        V1 = volume(model,p,T;phase=:l)
        V2 = volume(model,p,T;phase=:v)
        vvtol = max(abs(V1 - V1old)/V1, abs(V2 - V2old)/V2)
        pptol = abs(pold - p)/p
        if vvtol < vtol && i > 1
            @show i
            if check_valid_sat_pure(model,p,V1,V2,T)
                return (p,V1,V2)
            else
                return (nan,nan,nan)
            end
        end
        pold = p
        ϕl =  Fug(model,V1,T,p)
        ϕv =  Fug(model,V2,T,p)

        p = p*ϕl/ϕv

        if pptol < ptol && i > 1
            if check_valid_sat_pure(model,p,V1,V2,T)
                return (p,V1,V2)
            else
                return (nan,nan,nan)
            end
        end
        V1old = V1
        V2old = V2
    end
    return (nan,nan,nan)
end
=#


text:
using StaticArrays
using Test, StaticRanges, Dates, Documenter, IntervalSets

using StaticRanges: grow_beg, grow_beg!, grow_end, grow_end!,
    shrink_beg, shrink_beg!, shrink_end, shrink_end!

using ChainedFixes
using OffsetArrays
using OffsetArrays: IdOffsetRange

using ArrayInterface
using ArrayInterface: can_change_size, known_first, known_last, known_step, axes_types
using ArrayInterface.Static

# Uniqueness methods
using StaticRanges: ArrayInterface.ismutable

using Base: OneTo, step_hp

catch_nothing(x) = x
catch_nothing(x::Nothing) = 0

@test isempty(setdiff(detect_ambiguities(StaticRanges,Base,Core), detect_ambiguities(StaticArrays,Base,Core)))

#=
When using ranges find_all will produce either an AbstractUnitRange or GapRange.
Therefore, in order test type stability at this level we just ensure it only ever
comes down to one of these two types and convert them to Vector.
=#
to_vec(x::AbstractUnitRange{Int}) = collect(x)
to_vec(x::GapRange) = collect(x)
to_vec(x) = x

@test StaticRanges.ArrayInterface.ismutable(DynamicAxis(10))

@testset "checkindex" begin
    r = 1:5
    all_inbounds = 1:5
    top_outbounds = 1:6
    bottom_outbounds = 0:5
    @testset "checkindexlo" begin
        @test StaticRanges.checkindexlo(r, all_inbounds)
        @test StaticRanges.checkindexlo(r, top_outbounds)
        @test !StaticRanges.checkindexlo(r, bottom_outbounds)
        @test StaticRanges.checkindexlo(r, CartesianIndex((1,)))
    end

    @testset "checkindexhi" begin
        @test StaticRanges.checkindexhi(r, all_inbounds)
        @test !StaticRanges.checkindexhi(r, top_outbounds)
        @test StaticRanges.checkindexhi(r, bottom_outbounds)
        @test StaticRanges.checkindexhi(r, CartesianIndex((1,)))
    end
end

include("length_tests.jl")
include("find.jl")
include("gaprange_tests.jl")

include("promotion_tests.jl")
# FIXME include("nitty_gritty_promotion.jl")

@test first(GapRange(2:5, 7:10)) == 2
@test last(GapRange(2:5, 7:10)) == 10

@testset "Range Interface" begin
    for r in (DynamicAxis(10),
              mutable(1:10),
              MutableRange(StepRangeLen(1, 1, 10)))
        @testset "$(typeof(r))" begin
            @testset "first" begin
                @test @inferred(first(r)) == 1
            end
            @testset "firstindex" begin
                @test @inferred(firstindex(r)) == 1
            end
            @testset "step" begin
                @test @inferred(step(r)) == 1
            end
            @testset "last" begin
                @test @inferred(last(r)) == 10
            end
            @testset "lastindex" begin
                @test @inferred(lastindex(r)) == 10
            end
            @testset "length" begin
                @test @inferred(length(r)) == 10
            end
        end
    end
end
include("broadcast.jl")
include("dynamic_axis.jl")
include("intersect_tests.jl")

@testset "empty" begin
    for r in (static(OneTo(10)),
              DynamicAxis(10),
              srange(1.0, step=1, stop=10.0),
              mrange(1.0, step=1, stop=10.0),
              static(StepRange(1, 1, 10)),
              static(UnitRange(1, 10)),
              static(LinRange(1.5, 5.5, 9)))
        @test isempty(empty(r))
        if can_change_size(r)
            empty!(r)
            @test isempty(r)
        end
    end
end

#=
@testset "UnitRange" begin
    r = MutableRange(1:10)
    b = UnitRange(1, 10)
    rfloat = AbstractUnitRange{Float64}(r)
    @test eltype(rfloat) == Float64
    @test isa(rfloat, R)
    @test R{Int}(r) == r  # we don't use `===` because dynamic ranges should construct different ones
    @test R{Float64}(r) == R(1., 10.)
    @test eltype(R{Int}(UnitRange(UInt(1), UInt(10)))) == Int
    @test R(UnitRange(UInt(1), UInt(10))) == R(UInt(1), UInt(10))
    @test first(r) == r.start
    @test last(r) == r.stop

    @test intersect(r, r[2]) == intersect(b, b[2])
    @test intersect(r, r[2]) == intersect(b, b[2])

    @test @inferred(getindex(UnitMRange(1.0, 10.0), 2)) == 2.0

    @test_throws ErrorException r.notfield
    #=
    setproperty!(r, :start, 2)
    @test r == UnitMRange(2, 10)

    setproperty!(r, :stop, 8)
    @test r == UnitMRange(2, 8)

    @test_throws ErrorException setproperty!(r, :anything, 3)
    =#
end
=#

for frange in (mrange, srange)
    @testset "$frange" begin
        @test_throws ArgumentError frange(1)
        @test_throws ArgumentError frange(1, step=1)
        @test_throws ArgumentError frange(nothing)
        @test_throws ArgumentError frange(nothing, length=1)
        @test_throws ArgumentError frange(1, step=1, stop=1, length=1)
        if frange == mrange
            # cannot infer a static parameter from construction
            @testset "colon" begin
                @inferred(frange(10, step=1, stop=0))
                @inferred(frange(1, step=.2, stop=2))
                @inferred(frange(1., step=.2, stop=2.))
                @inferred(frange(2, step=-.2, stop=1))
                @inferred(frange(1, 0))
                @inferred(frange(0.0, -0.5))
            end
        end
        @testset "indexing" begin
            L32 = frange(Int32(1), stop=Int32(4), length=4)
            L64 = frange(Int64(1), stop=Int64(4), length=4)
            @test @inferred(L32[1]) === 1.0 && @inferred(L64[1]) === 1.0
            @test L32[2] == 2 && L64[2] == 2
            @test L32[3] == 3 && L64[3] == 3
            @test L32[4] == 4 && L64[4] == 4
            @test frange(1.0, stop=2.0, length=2)[1] === 1.0
            @test frange(1.0f0, stop=2.0f0, length=2)[1] === 1.0f0
            @test frange(Float16(1.0), stop=Float16(2.0), length=2)[1] === Float16(1.0)

            let r = frange(5, step=-1, stop=1)
                @test r[1]==5
                @test r[2]==4
                @test r[3]==3
                @test r[4]==2
                @test r[5]==1
            end
            @test @inferred(frange(0.1, step=0.1, stop=0.3)[2]) === 0.2
            @test @inferred(frange(0.1f0, step=0.1f0, stop=0.3f0)[2]) === 0.2f0

            @test frange(1, 5)[1:4] == 1:4
            @test frange(1.0, 5)[1:4] == 1.0:4
            @test frange(2, 6)[1:4] == 2:5
            @test frange(1, 6)[2:5] == 2:5
            @test frange(1, 6)[2:2:5] == 2:2:4
            @test frange(1, step=2, stop=13)[2:6] == 3:2:11
            @test frange(1, step=2, stop=13)[2:3:7] == 3:6:13

            @test isempty(frange(1, 4)[5:4])
            #@test_throws BoundsError frange(1:10)[8:-1:-2]
        end
       @testset "keys/values" begin
            keytype_is_correct(r) = keytype(r) == eltype(keys(r))
            valtype_is_correct(r) = valtype(r) == eltype(values(r))
            @test keytype_is_correct(frange(1, 3))
            @test keytype_is_correct(frange(1, step=.3, stop=4))
            @test keytype_is_correct(frange(.11, step=.1, stop=.3))
            @test keytype_is_correct(frange(Int8(1), Int8(5)))
            @test keytype_is_correct(frange(Int16(1), Int8(5)))
            @test keytype_is_correct(frange(Int16(1), step=Int8(3), stop=Int8(5)))
            @test keytype_is_correct(frange(Int8(1), step=Int16(3), stop=Int8(5)))
            @test keytype_is_correct(frange(Int8(1), step=Int8(3), stop=Int16(5)))
            @test keytype_is_correct(frange(Int64(1), Int64(5)))
            @test keytype_is_correct(frange(Int64(1), Int64(5)))
            @test keytype_is_correct(frange(Int128(1), Int128(5)))
            @test valtype_is_correct(1:3)
            @test valtype_is_correct(1:.3:4)
            @test valtype_is_correct(.1:.1:.3)
            @test valtype_is_correct(frange(Int8(1), Int8(5)))
            @test valtype_is_correct(frange(Int16(1), Int8(5)))
            @test valtype_is_correct(frange(Int16(1), step=Int8(3), stop=Int8(5)))
            @test valtype_is_correct(frange(Int8(1), step=Int16(3), stop=Int8(5)))
            @test valtype_is_correct(frange(Int8(1), step=Int8(3), stop=Int16(5)))
            @test valtype_is_correct(frange(Int64(1), Int64(5)))
            @test valtype_is_correct(frange(Int64(1), Int64(5)))
            @test valtype_is_correct(frange(Int128(1), Int128(5)))

            if frange isa typeof(mrange)
                @test keytype_is_correct(DynamicAxis(4))
                @test keytype_is_correct(DynamicAxis(Int32(4)))
                @test valtype_is_correct(DynamicAxis(4))
                @test valtype_is_correct(DynamicAxis(Int32(4)))
            else
                @test keytype_is_correct(static(OneTo((4))))
                @test keytype_is_correct(static(OneTo((Int32(4)))))
                @test valtype_is_correct(static(OneTo((4))))
                @test valtype_is_correct(static(OneTo((Int32(4)))))
            end
        end

        @testset "findall(::Base.Fix2{typeof(in)}, ::Array)" begin
            @test findall(in(3:20), [5.2, 3.3]) == findall(in(Vector(3:20)), [5.2, 3.3])

            # FIXME shouldn't both be reversed?
            let span = frange(5, 20),
                r = frange(-7, step=3, stop=42)
                @test findall(in(span), r) == 5:10

                r = frange(15, step=-2, stop=-38)
                @test findall(in(span), r) == 6:-1:1
            end
        end

        if VERSION > v"1.2"
            @testset "issubset" begin
                @test issubset(frange(1, 3), 1:typemax(Int)) #32461
                @test issubset(frange(1, 3), 1:3)
                @test issubset(frange(1, 3), 1:4)
                @test issubset(frange(1, 3), 0:3)
                @test issubset(frange(1, 3), 0:4)
                @test !issubset(frange(1, 5), 2:5)
                @test !issubset(frange(1, 5), 1:4)
                @test !issubset(frange(1, 5), 2:4)
                @test issubset(frange(1, step=3, stop=10), 1:10)
                @test !issubset(frange(1, 10), 1:3:10)

                @test issubset(1:3, frange(1, typemax(Int))) #32461
                @test issubset(1:3, frange(1, 3))
                @test issubset(1:3, frange(1, 4))
                @test issubset(1:3, frange(0, 3))
                @test issubset(1:3, frange(0, 4))
                @test !issubset(1:5, frange(2, 5))
                @test !issubset(1:5, frange(1, 4))
                @test !issubset(1:5, frange(2, 4))
                @test issubset(1:3:10, frange(1, 10))
                @test !issubset(1:10, frange(1, step=3, stop=10))
                if frange isa typeof(mrange)
                    @test issubset(DynamicAxis(5), DynamicAxis(10))
                    @test !issubset(DynamicAxis(10), DynamicAxis(5))
                    @test issubset(DynamicAxis(5), DynamicAxis(10))
                    @test !issubset(DynamicAxis(10), DynamicAxis(5))
                else
                    @test issubset(static(OneTo(5)), static(OneTo((10))))
                    @test !issubset(static(OneTo(10)), static(OneTo((5))))
                    @test issubset(static(OneTo(5)), static(OneTo((10))))
                    @test !issubset(static(OneTo(10)), static(OneTo((5))))
                end

            end
        end

        @testset "in" begin
            @test 0 in frange(UInt(0), step=100, stop=typemax(UInt))
            @test last(frange(UInt(0), step=100, stop=typemax(UInt))) in frange(UInt(0), step=100, stop=typemax(UInt))
            @test -9223372036854775790 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test -9223372036854775690 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test -90 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 10 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 110 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 9223372036854775610 in frange(-9223372036854775790, step=100, stop=9223372036854775710)
            @test 9223372036854775710 in frange(-9223372036854775790, step=100, stop=9223372036854775710)

            @test !(3.5 in frange(1, 5))
            @test (3 in frange(1, 5))
            @test (3 in frange(5, step=-1, stop=1))
            #@test (3 in 3+0*(1:5))
            #@test !(4 in 3+0*(1:5))

            let r = frange(0.0, step=0.01, stop=1.0)
                @test (r[30] in r)
            end
            let r = frange((-4*Int64(maxintfloat(Int === Int32 ? Float32 : Float64))), 5)
                @test (3 in r)
                @test (3.0 in r)
            end

            @test !(1 in frange(1, 0))
            @test !(1.0 in frange(1.0, 0.0))
        end

        @testset "in() works across types, including non-numeric types (#21728)" begin
            @test 1//1 in frange(1, 3)
            @test 1//1 in frange(1.0, 3.0)
            @test !(5//1 in frange(1, 3))
            @test !(5//1 in frange(1.0, 3.0))
            @test Complex(1, 0) in frange(1, 3)
            @test Complex(1, 0) in frange(1.0, 3.0)
            @test Complex(1.0, 0.0) in frange(1, 3)
            @test Complex(1.0, 0.0) in frange(1.0, 3.0)
            @test !(Complex(1, 1) in frange(1, 3))
            @test !(Complex(1, 1) in frange(1.0, 3.0))
            @test !(Complex(1.0, 1.0) in frange(1, 3))
            @test !(Complex(1.0, 1.0) in frange(1.0, 3.0))
            @test !(π in frange(1, 3))
            @test !(π in frange(1.0, 3.0))
            @test !("a" in frange(1, 3))
            @test !("a" in frange(1.0, 3.0))
            @test !(1 in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
            @test !(Complex(1, 0) in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
            @test !(π in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
            @test !("a" in frange(Date(2017, 01, 01), step=Dates.Day(1), stop=Date(2017, 01, 05)))
        end
        @testset "sums of ranges" begin
            @test sum(frange(1, 100)) == 5050
            @test sum(frange(0, 100)) == 5050
            @test sum(frange(-100, 100)) == 0
            @test sum(frange(0, step=2, stop=100)) == 2550
        end
        @testset "Tricky sums of StepRangeLen #8272" begin
            @test sum(frange(10000., step=-0.0001, stop=0)) == 5.00000005e11
            @test sum(frange(0, step=0.001, stop=1)) == 500.5
            @test sum(frange(0, step=0.000001, stop=1)) == 500000.5
            @test sum(frange(0, step=0.1, stop=10)) == 505.
        end
        @testset "operations between ranges and arrays" begin
            @test all(([frange(1, 5);] + (frange(5, step=-1, stop=1))) .== 6)
            @test all(((frange(5, step=-1, stop=1)) + [frange(1, 5);]) .== 6)
            @test all(([frange(1, 5);] - (frange(1, 5))) .== 0)
            @test all((frange(1, 5) - [frange(1, 5);]) .== 0)
        end

        @testset "broadcasted operations with scalars" begin
            @test broadcast(-, frange(1, 3)) == -1:-1:-3
            @test broadcast(-, frange(1, 3), 2) == -1:1
            @test broadcast(-, frange(1, 3), 0.25) == 1-0.25:3-0.25
            @test broadcast(+, frange(1, 3)) == 1:3
            @test broadcast(+, frange(1, 3), 2) == 3:5
            @test broadcast(+, frange(1, 3), 0.25) == 1+0.25:3+0.25
            @test broadcast(+, frange(1, step=2, stop=6), 1) == 2:2:6
            @test broadcast(+, frange(1, step=2, stop=6), 0.3) == 1+0.3:2:5+0.3
            @test broadcast(-, frange(1, step=2, stop=6), 1) == 0:2:4
            @test broadcast(-, frange(1, step=2, stop=6), 0.3) == 1-0.3:2:5-0.3
            @test broadcast(-, 2, frange(1, 3)) == 1:-1:-1
        end

        @testset "loops involving typemin/typemax" begin
            n = 0
            s = 0
            # loops ending at typemax(Int)
            for i = frange((typemax(Int)-1), typemax(Int))
                s += 1
                @test s <= 2
            end
            @test s == 2

            s = 0
            for i = frange((typemax(Int)-2), (typemax(Int)-1))
                s += 1
                @test s <= 2
            end
            @test s == 2

            s = 0
            for i = frange(typemin(Int), (typemin(Int)+1))
                s += 1
                @test s <= 2
            end
            @test s == 2

            # loops covering the full range of integers
            s = 0
            for i = frange(typemin(UInt8), typemax(UInt8))
                s += 1
            end
            @test s == 256

            s = 0
            for i = frange(typemin(UInt), typemax(UInt))
                i == 10 && break
                s += 1
            end
            @test s == 10

            s = 0
            for i = frange(typemin(UInt8), step=one(UInt8), stop=typemax(UInt8))
                s += 1
            end
            @test s == 256

            s = 0
            for i = frange(typemin(UInt), step=1, stop=typemax(UInt))
                i == 10 && break
                s += 1
            end
            @test s == 10

            # loops past typemax(Int)
            n = 0
            s = Int128(0)
            for i = frange(typemax(UInt64)-2, typemax(UInt64))
                n += 1
                s += i
            end
            @test n == 3
            @test s == 3*Int128(typemax(UInt64)) - 3

            # loops over empty ranges
            s = 0
            for i = 0xff:0x00
                s += 1
            end
            @test s == 0

            s = 0
            for i = frange(Int128(typemax(Int128)), Int128(typemin(Int128)))
                s += 1
            end
            @test s == 0
        end

        function range_fuzztests(::Type{T}, niter, nrange) where {T}
            for i = 1:niter, n in nrange
                strt, Δ = randn(T), randn(T)
                Δ == 0 && continue
                stop = strt + (n-1)*Δ
                # `n` is not necessarily unique s.t. `strt + (n-1)*Δ == stop`
                # so test that `length(strt:Δ:stop)` satisfies this identity
                # and is the closest value to `(stop-strt)/Δ` to do so
                lo = hi = n
                while strt + (lo-1)*Δ == stop; lo -= 1; end
                while strt + (hi-1)*Δ == stop; hi += 1; end
                m = clamp(round(Int, (stop-strt)/Δ) + 1, lo+1, hi-1)
                r = strt:Δ:stop
                @test m == length(r)
                @test strt == first(r)
                @test Δ == step(r)
                @test_skip stop == last(r)
                l = range(strt, stop=stop, length=n)
                @test n == length(l)
                @test strt == first(l)
                @test stop  == last(l)
            end
        end

        @testset "range fuzztests for $T" for T = (Float32, Float64,)
            range_fuzztests(T, 2^15, frange(1, 5))
        end

        @testset "range with very large endpoints for type $T" for T = (Float32, Float64)
            largeint = Int(min(maxintfloat(T), typemax(Int)))
            a = floatmax()
            for i = 1:5
                @test [frange(a, stop=a, length=1);] == [a]
                @test [frange(-a, stop=-a, length=1);] == [-a]
                b = floatmax()
                for j = 1:5
                    @test [frange(-a, stop=b, length=0);] == []
                    @test [frange(-a, stop=b, length=2);] == [-a,b]
                    @test [frange(-a, stop=b, length=3);] == [-a,(b-a)/2,b]
                    @test [frange(a, stop=-b, length=0);] == []
                    @test [frange(a, stop=-b, length=2);] == [a,-b]
                    @test [frange(a, stop=-b, length=3);] == [a,(a-b)/2,-b]
                    for c = largeint-3:largeint
                        s = range(-a, stop=b, length=c)
                        @test first(s) == -a
                        @test last(s) == b
                        @test length(s) == c
                        s = range(a, stop=-b, length=c)
                        @test first(s) == a
                        @test last(s) == -b
                        @test length(s) == c
                    end
                    b = prevfloat(b)
                end
                a = prevfloat(a)
            end
        end

        @testset "ranges with very small endpoints for type $T" for T = (Float32, Float64)
            z = zero(T)
            u = eps(z)
            @test first(frange(u, stop=u, length=0)) == u
            @test last(frange(u, stop=u, length=0)) == u
            @test first(frange(-u, stop=u, length=0)) == -u
            @test last(frange(-u, stop=u, length=0)) == u
            @test [frange(-u, stop=u, length=0);] == []
            @test [frange(-u, stop=-u, length=1);] == [-u]
            @test [frange(-u, stop=u, length=2);] == [-u,u]
            @test [frange(-u, stop=u, length=3);] == [-u,0,u]
            @test first(frange(-u, stop=-u, length=0)) == -u
            @test last(frange(-u, stop=-u, length=0)) == -u
            @test first(frange(u, stop=-u, length=0)) == u
            @test last(frange(u, stop=-u, length=0)) == -u
            @test [frange(u, stop=-u, length=0);] == []
            @test [frange(u, stop=u, length=1);] == [u]
            @test [frange(u, stop=-u, length=2);] == [u,-u]
            @test [frange(u, stop=-u, length=3);] == [u,0,-u]
            v = frange(-u, stop=u, length=12)
            @test length(v) == 12
            @test [-3u:u:3u;] == [frange(-3u, stop=3u, length=7);] == [-3:3;].*u
            @test [3u:-u:-3u;] == [frange(3u, stop=-3u, length=7);] == [3:-1:-3;].*u
        end

        @testset "tricky floating-point ranges" begin
            for (start, step, stop, len) in ((1, 1, 3, 3),
                                             (0, 1, 3, 4),
                                             (3, -1, -1, 5),
                                             (1, -1, -3, 5),
                                             (0, 1, 10, 11),
                                             (0, 7, 21, 4),
                                             (0, 11, 33, 4),
                                             (1, 11, 34, 4),
                                             (0, 13, 39, 4),
                                             (1, 13, 40, 4),
                                             (11, 11, 33, 3),
                                             (3, 1, 11, 9),
                                             (0, 10, 55, 0),
                                             (0, -1, 5, 0),
                                             (0, 10, 5, 0),
                                             (0, 1, 5, 0),
                                             (0, -10, 5, 0),
                                             (0, -10, 0, 1),
                                             (0, -1, 1, 0),
                                             (0, 1, -1, 0),
                                             (0, -1, -10, 11))
                r = frange(start/10, step=step/10, stop=stop/10)
                a = Vector(frange(start, step=step, stop=stop))./10
                ra = Vector(r)

                @test r == a
                @test isequal(r, a)

                @test r == ra
                @test isequal(r, ra)

                @test hash(r) == hash(a)
                @test hash(r) == hash(ra)

                if len > 0
                    l = frange(start/10, stop=stop/10, length=len)
                    la = Vector(l)

                    @test a == l
                    @test r == l
                    @test isequal(a, l)
                    @test isequal(r, l)

                    @test l == la
                    @test isequal(l, la)

                    @test hash(l) == hash(a)
                    @test hash(l) == hash(la)
                end
            end

            @test 1.0:1/49:27.0 == range(1.0, stop=27.0, length=1275) == [49:1323;]./49
            @test isequal(1.0:1/49:27.0, range(1.0, stop=27.0, length=1275))
            @test isequal(frange(1.0, step=1/49, stop=27.0), Vector(49:1323)./49)
            @test hash(1.0:1/49:27.0) == hash(frange(1.0, stop=27.0, length=1275)) == hash(Vector(frange(49, 1323))./49)

            @test [frange(prevfloat(0.1), step=0.1, stop=0.3);] == [prevfloat(0.1), 0.2, 0.3]
            @test [frange(nextfloat(0.1), step=0.1, stop=0.3);] == [nextfloat(0.1), 0.2]
            @test [frange(prevfloat(0.0), step=0.1, stop=0.3);] == [prevfloat(0.0), 0.1, 0.2]
            @test [frange(nextfloat(0.0), step=0.1, stop=0.3);] == [nextfloat(0.0), 0.1, 0.2]
            @test [frange(0.1, step=0.1, stop=prevfloat(0.3));] == [0.1, 0.2]
            @test [frange(0.1, step=0.1, stop=nextfloat(0.3));] == [0.1, 0.2, nextfloat(0.3)]
            @test [frange(0.0, step=0.1, stop=prevfloat(0.3));] == [0.0, 0.1, 0.2]
            @test [frange(0.0, step=0.1, stop=nextfloat(0.3));] == [0.0, 0.1, 0.2, nextfloat(0.3)]
            @test [frange(0.1, step=prevfloat(0.1), stop=0.3);] == [0.1, 0.2, 0.3]
            @test [frange(0.1, step=nextfloat(0.1), stop=0.3);] == [0.1, 0.2]
            @test [frange(0.0, step=prevfloat(0.1), stop=0.3);] == [0.0, prevfloat(0.1), prevfloat(0.2), 0.3]
            @test [frange(0.0, step=nextfloat(0.1), stop=0.3);] == [0.0, nextfloat(0.1), nextfloat(0.2)]
        end

        @testset "overflowing sums (see #5798)" begin
            if Sys.WORD_SIZE == 64
                @test sum(frange(Int128(1), 10^18)) == div(10^18 * (Int128(10^18)+1), 2)
                @test sum(frange(Int128(1), 10^18-1)) == div(10^18 * (Int128(10^18)-1), 2)
            else
                @test sum(frange(Int64(1), 10^9)) == div(10^9 * (Int64(10^9)+1), 2)
                @test sum(frange(Int64(1), 10^9-1)) == div(10^9 * (Int64(10^9)-1), 2)
            end
        end

        @testset "issue #20373 (unliftable ranges with exact end points)" begin
            @test [3*0.05:0.05:0.2;]    == [frange(3*0.05, stop=0.2, length=2);]   == [3*0.05,0.2]
            @test [0.2:-0.05:3*0.05;]   == [frange(0.2, stop=3*0.05, length=2);]   == [0.2,3*0.05]
            @test [-3*0.05:-0.05:-0.2;] == [frange(-3*0.05, stop=-0.2, length=2);] == [-3*0.05,-0.2]
            @test [-0.2:0.05:-3*0.05;]  == [frange(-0.2, stop=-3*0.05, length=2);] == [-0.2,-3*0.05]
        end

        @testset "issue #6973" begin
            r1 = frange(1.0, step=0.1, stop=2.0)
            r2 = frange(1.0f0, step=0.2f0, stop=3.0f0)
            r3 = frange(1, step=2, stop=21)
            @test r1 + r1 == 2*r1
            @test r1 + r2 == 2.0:0.3:5.0
            @test (r1 + r2) - r2 == r1
            @test r1 + r3 == convert(StepRangeLen{Float64}, r3) + r1
            @test r3 + r3 == 2 * r3
        end

       @testset "Inexact errors on 32 bit architectures. #22613" begin
            @test first(frange(log(0.2), stop=log(10.0), length=10)) == log(0.2)
            @test last(frange(log(0.2), stop=log(10.0), length=10)) == log(10.0)
            # not used internally for StaticRanges
            #@test length(Base.floatrange(-3e9, 1.0, 1, 1.0)) == 1
        end

        # issue #7426
        @test [frange(typemax(Int), step=1, stop=typemax(Int));] == [typemax(Int)]

        @testset "issue #7387" begin
            for r in (frange(0, 1), frange(0.0, 1.0))
                local r
                @test [r .+ im;] == [r;] .+ im
                @test [r .- im;] == [r;] .- im
                @test [r * im;] == [r;] * im
                @test [r / im;] == [r;] / im
            end
        end

        # near-equal ranges
        @test frange(0.0, step=0.1, stop=1.0) != 0.0f0:0.1f0:1.0f0

        #issue #7484
        let r7484 = frange(0.1, step=0.1, stop=1)
            @test [reverse(r7484);] == reverse([r7484;])
        end

        # issue #2959
        @test frange(1.0, 1.5) == 1.0:1.0:1.5 == 1.0:1.0

        @testset "comparing UnitRanges and OneTo" begin
            @test frange(1, step=2, stop=10) == 1:2:10 != 1:3:10 != 1:3:13 != frange(2, step=3, stop=13) == 2:3:11 != frange(2, 11)
            @test frange(1, step=1, stop=10) == 1:10 == 1:10 == DynamicAxis(10) == static(OneTo(10))
            @test 1:10 != frange(2, 10) != 2:11 != Base.OneTo(11)
            @test DynamicAxis(10) != static(OneTo(11)) != frange(1, 10)
        end

        @testset "issue #7114" begin
            let r = frange(-0.004532318104333742, step=1.2597349521122731e-5, stop=0.008065031416788989)
                @test length(r[1:end-1]) == length(r) - 1
                @test isa(r[1:2:end],AbstractRange) && length(r[1:2:end]) == div(length(r)+1, 2)
                @test r[3:5][2] ≈ r[4]
                @test r[5:-2:1][2] ≈ r[3]
                @test_throws BoundsError r[0:10]
                @test_throws BoundsError r[1:10000]
            end

            let r = frange(1/3, stop=5/7, length=6)
                @test length(r) == 6
                @test r[1] == 1/3
                @test abs(r[end] - 5/7) <= eps(5/7)
            end

            let r = frange(0.25, stop=0.25, length=1)
                @test length(r) == 1
                @test_throws ArgumentError frange(0.25, stop=0.5, length=1)
            end
        end

        # Preservation of high precision upon addition
        let r = range(-0.1, step=0.1, stop=0.3) + broadcast(+, -0.3:0.1:0.1, 1e-12)
            @test r[3] == 1e-12
        end

        @testset "range with 1 or 0 elements (whose step length is NaN)" begin
            @test issorted(frange(1, stop=1, length=0))
            @test issorted(frange(1, stop=1, length=1))
        end
        @testset "indexing range with empty range (#4309)" begin
            @test frange(3, 6)[5:4] == 7:6
            @test_throws BoundsError frange(3, 6)[5:5]
            @test_throws BoundsError frange(3, 6)[5]
            @test frange(0, step=2, stop=10)[7:6] == 12:2:10
            @test_throws BoundsError (0:2:10)[7:7]
        end
        # indexing with negative ranges (#8351)
        for a=AbstractRange[3:6, frange(0, step=2, stop=10)], b=AbstractRange[frange(0, 1), frange(2, step=-1, stop=0)]
            @test_throws BoundsError a[b]
        end
        @testset "sort/sort!/partialsort" begin
            @test sort(frange(1, 2)) == UnitRange(1,2)
            if frange == mrange
                @test sort!(frange(1, 2)) == UnitRange(1,2)
            end
            @test sort(frange(1, 10), rev=true) == 10:-1:1
            @test sort(frange(-3, 3), by=abs) == [0,-1,1,-2,2,-3,3]
            @test partialsort(frange(1, 10), 4) == 4
        end

        if frange == mrange
            # this works with `srange` but takes way too long to run.
            function loop_range_values(::Type{T}) where T
                for a = -5:25,
                    s = [-5:-1; 1:25; ],
                    d = 1:25,
                    n = -1:15

                    denom = convert(T, d)
                    strt = convert(T, a)/denom
                    Δ     = convert(T, s)/denom
                    stop  = convert(T, (a + (n - 1) * s)) / denom
                    vals  = T[a:s:(a + (n - 1) * s); ] ./ denom
                    r = strt:Δ:stop
                    @test [r;] == vals
                    @test [frange(strt, stop=stop, length=length(r));] == vals
                    n = length(r)
                    @test [r[1:n];] == [r;]
                    @test [r[2:n];] == [r;][2:end]
                    @test [r[1:3:n];] == [r;][1:3:n]
                    @test [r[2:2:n];] == [r;][2:2:n]
                    @test [r[n:-1:2];] == [r;][n:-1:2]
                    @test [r[n:-2:1];] == [r;][n:-2:1]
                end
            end

            @testset "issue #7420 for type $T" for T = (Float32, Float64,) # BigFloat),
                loop_range_values(T)
            end
        end

        # comparing and hashing ranges
        @testset "comparing and hashing ranges" begin
            Rs = AbstractRange[]
            for r in (frange(1, 1),
                      frange(1, step=1, stop=1),
                      frange(1, 2),
                      frange(1, step=1, stop=2),
                      map(Int32, frange(1, step=3, stop=17)),
                      map(Int64, frange(1, step=3, stop=17)),
                      frange(1, 0),
                      frange(1, step=-1, stop=0),
                      frange(17, 0, step=-3),
                      frange(0.0, step=0.1, stop=1.0),
                      map(Float32, frange(0.0, step=0.1, stop=1.0)),
                      frange(1.0, step=eps(), stop=1.0) .+ 10eps(),
                      frange(9007199254740990., step=1.0, stop=9007199254740994),
                      frange(0, stop=1, length=20),
                      map(Float32, frange(0, stop=1, length=20)))
                local r
                ar = Vector(r)
                @test r == ar
                @test isequal(r,ar)
                @test hash(r) == hash(ar)
                for s in Rs
                    as = Vector(s)
                    @test isequal(r,s) == (hash(r)==hash(s))
                    @test (r==s) == (ar==as)
                end
            end
        end
    end
end

#=  TODO mod with ranges
@testset "mod with ranges" begin
    for n in -10:10
        @test mod(n, 0:4) == mod(n, 5)
        @test mod(n, 1:5) == mod1(n, 5)
        @test mod(n, 2:6) == 2 + mod(n-2, 5)
        @test mod(n, Base.OneTo(5)) == mod1(n, 5)
    end
    @test mod(Int32(3), 1:5) == 3
    @test mod(big(typemax(Int))+99, 0:4) == mod(big(typemax(Int))+99, 5)
    @test_throws MethodError mod(3.141, 1:5)
    @test_throws MethodError mod(3, UnitRange(1.0,5.0))
    @test_throws MethodError mod(3, 1:2:7)
    @test_throws DivideError mod(3, 1:0)
end
=#
for frange in (mrange,srange)
    @testset "reverse-$frange" begin
        @test reverse(reverse(frange(1, 10))) == 1:10
        @test reverse(reverse(frange(typemin(Int), typemax(Int)))) == typemin(Int):typemax(Int)
        @test reverse(reverse(frange(typemin(Int), step=2, stop=typemax(Int)))) == typemin(Int):2:typemax(Int)
    end
end

r = MutableRange(1:1:10)
@test reverse!(r) == reverse(1:1:10)
@test reverse!(r) == 1:1:10

@testset "resize tests" begin
    include("resize_tests.jl")
end

if VERSION > v"1.6" && sizeof(Int) === 8
    @testset "docs" begin
        doctest(StaticRanges)
    end
end
include("count_tests.jl")

#=
#@test 1.0:(.3-.1)/.1 == 1.0:2.0
=#

#=
detect_ambiguities(StaticRanges,Base,Core)

=#


text:
"""
    `evaluate(solver::AbstractCFRSolver)`

Evaluate full tree traversed by CFR solver. \n
Returns tuple corresponding to utilities for both players.
"""
function evaluate(solver::AbstractCFRSolver)
    finalize_strategies!(solver)

    ih = initialhist(solver.game)

    p1_eval = evaluate(solver, ih, 1, 0, 1.0, 1.0)
    p2_eval = evaluate(solver, ih, 2, 0, 1.0, 1.0)

    return (p1_eval, p2_eval)
end

function evaluate(sol::AbstractCFRSolver, p::Int)
    ih = initialhist(sol.game)
    return evaluate(sol, ih, p, 0, 1.0, 1.0)
end

function evaluate(solver::AbstractCFRSolver, h, i, t, π_1, π_2)
    game = solver.game
    if isterminal(game, h)
        return utility(game, i, h)
    elseif iszero(player(game, h)) # chance player
        A = chance_actions(game, h)
        s = 0.0
        for a in A
            s += evaluate(solver, next_hist(game, h, a), i, t, π_1, π_2)
        end
        return s / length(A)
    end

    I = infokey(game, h)
    A = actions(game, h)

    v_σ = 0.0

    σ = strategy(solver, I)
    for (k,a) in enumerate(A)
        v_σ_Ia = 0.0
        h′ = next_hist(game, h, a)
        if player(game, h) === 1
            v_σ_Ia = evaluate(solver, h′, i, t, σ[k]*π_1, π_2)
        else
            v_σ_Ia = evaluate(solver, h′, i, t, π_1, σ[k]*π_2)
        end
        v_σ += σ[k]*v_σ_Ia
    end

    return v_σ
end

"""
    `MonteCarloEvaluate(solver::AbstractCFRSolver, N::Int)`

Monte Carlo evaluation sampling chance player actions. \n
Returns tuple corresponding to utilities for both players.
"""
function MonteCarloEvaluate(solver::AbstractCFRSolver, N::Int)
    finalize_strategies!(solver)

    p1_eval = 0.0
    p2_eval = 0.0

    ih = initialhist(solver.game)
    for _ in 1:N
        p1_eval += MonteCarloEvaluate(solver, ih, 1, 0, 1.0, 1.0)
        p2_eval += MonteCarloEvaluate(solver, ih, 2, 0, 1.0, 1.0)
    end

    p1_eval /= N
    p2_eval /= N

    return (p1_eval, p2_eval)
end

function MonteCarloEvaluate(solver::AbstractCFRSolver, h, i, t, π_1, π_2)
    game = solver.game
    if isterminal(game, h)
        return utility(game, i, h)
    elseif iszero(player(game, h)) # chance player
        a = chance_action(game, h)
        h′ = next_hist(game,h,a)
        return MonteCarloEvaluate(solver, h′, i, t, π_1, π_2)
    end

    I = infokey(game, h)
    A = actions(game, h)

    v_σ = 0.0
    σ = strategy(solver, I)
    for (k,a) in enumerate(A)
        v_σ_Ia = 0.0
        h′ = next_hist(game, h, a)
        if player(game, h) === 1
            v_σ_Ia = MonteCarloEvaluate(solver, h′, i, t, σ[k]*π_1, π_2)
        else
            v_σ_Ia = MonteCarloEvaluate(solver, h′, i, t, π_1, σ[k]*π_2)
        end
        v_σ += σ[k]*v_σ_Ia
    end

    return v_σ
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
@testset "Interval matrix exponential" begin

    @test quadratic_expansion(-3..3, 1.0, 2.0) == Interval(-0.125, 21)

    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])

    for i in 0:4
        _truncated_exponential_series(m, 1.0, i)
    end

    overapp1 = exp_overapproximation(m, 1.0, 4)
    overapp2 = horner(m, 10)
    overapp3 = scale_and_square(m, 5, 1.0, 4)
    underapp = exp_underapproximation(m, 1.0, 4)

    @test underapp isa IntervalMatrix
    for overapp in [overapp1, overapp2, overapp3]
        @test overapp isa IntervalMatrix
    end
end

@testset "Interval matrix correction terms" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    f = correction_hull(m, 1e-3, 5)
    f2 = input_correction(m, 1e-3, 5)
    f = correction_hull(mid(m), 1e-3, 5)
end

@testset "Interval matrix square" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])

    a = m * m
    b = square(m)
    @test b ⊆ a
end

@testset "Interval matrix power" begin
    m = IntervalMatrix([2.0..2.0 2.0..3.0; 0.0..0.0 -1.0..1.0])
    pow = IntervalMatrixPower(m)
    increment!(pow)  # power of two
    increment!(pow, algorithm="multiply")
    increment!(pow)  # power of two
    increment!(pow, algorithm="power")
    increment!(pow, algorithm="decompose_binary")
    increment!(pow, algorithm="intersect")
end

@testset "Transmission line model" begin
    A = transmission_line() |> Matrix
    expA = exp(A)
    @test opnorm(expA, Inf)< 1e-6
    Aint = IntervalMatrix(interval.(A))
    expA_int = exp_overapproximation(Aint, 1.0, 10)
    @test opnorm(expA_int, Inf) > 1e122
end


text:
module TestPartial
using Test
using Configurations

@option struct MyOption
    x::Int
    y::Int = x + 1
end

@test from_dict(MyOption, Dict("x"=>1)) == MyOption(1, 2)

end


text:
module StylishArrays

abstract type StylishArray{T, N} <: AbstractArray{T, N} end

using Base.Broadcast: Broadcasted, AbstractArrayStyle, preprocess, BroadcastStyle, DefaultArrayStyle
using Base.Broadcast: instantiate
using LinearAlgebra
using Swizzles
using Swizzles.NullArrays
using Swizzles.ScalarArrays

export StylishArray, myidentity, Styled

struct Styled{Style, Arg}
    arg::Arg
end

@inline Styled{Style}(arg::Arg) where {Style, Arg} = Styled{Style, Arg}(arg)
@inline Styled(arg) = Styled{typeof(BroadcastStyle(typeof(arg)))}(arg)

myidentity(x) = x
@inline Base.similar(src::Styled{Style}, args...) where {Style} = similar(Broadcasted{Style}(identity, (src.arg,), axes(src.arg)), args...)
Base.@propagate_inbounds Base.copy(src::Styled) = copyto!(similar(src), src)
Base.@propagate_inbounds function Base.copyto!(dst::AbstractArray, src::Styled{Style}) where {Style}
    @boundscheck axes(dst) == axes(src.arg) || error("TODO")
    copyto!(dst, Broadcasted{Style}(myidentity, (src.arg,), axes(dst)))
end

"""
    StylishArray <: AbstractArray

A convenience type for defining array types which are implemented using
`similar(::Broadcasted)`, `copyto!(::Any, ::Broadcasted)` and
`copy(::Broadcasted)`. Consult the [Interfaces chapter](@ref
man-interfaces-broadcasting) on broadcasting for more info about broadcast. Many
Base functions are implemented, ` including `map`, `getindex`, and `foreach`.
Power users who wish to avoid broadcasting semantics can also intercept
`copyto!(::Any, ::Styled)`, `copy(::Styled)`, and `similar(::Styled)`

See also: [`similar`](@ref), [`copy`](@ref), [`copyto!`](@ref)
"""
abstract type StylishArray{T, N} <: AbstractArray{T, N} end

#These two are more for intercept purposes than StylishArray purposes.
Base.@propagate_inbounds Base.Broadcast.materialize(A::StylishArray) = identity.(A)
Base.@propagate_inbounds Base.Broadcast.materialize!(dst, A::StylishArray) = dst .= A

#Beware infinite recursion!

Base.@propagate_inbounds Base.similar(src::StylishArray) = similar(Styled(src))

Base.@propagate_inbounds Base.copy(src::StylishArray) = copy(Styled(src))

Base.@propagate_inbounds Base.copyto!(dst::StylishArray, src::AbstractArray) = _copyto!(dst, src)
Base.@propagate_inbounds Base.copyto!(dst::AbstractArray, src::StylishArray) = _copyto!(dst, src)
Base.@propagate_inbounds Base.copyto!(dst::StylishArray, src::StylishArray) = _copyto!(dst, src)
Base.@propagate_inbounds function _copyto!(dst::AbstractArray, src)
    @boundscheck axes(dst) == axes(src) || error("I don't like it when copyto! has mismatched axes lololol")
    #dst .= myidentity.(reshape(arrayify(src), axes(dst)))
    copyto!(dst, Styled(src))
    return dst
end

#The following nonsense means that generated arrays can override getindex or they can override copyto!(view)
Base.@propagate_inbounds Base.getindex(arr::StylishArray, I::Integer)::eltype(arr) = _getindex(arr, I)
Base.@propagate_inbounds Base.getindex(arr::StylishArray, I...) = _getindex(arr, I...)

Base.@propagate_inbounds function _getindex(arr, I...)
    src = Styled(view(arr, I...))
    dst = copyto!(similar(src), src)
    if ndims(dst) == 0
        return dst[]
    else
        return dst
    end
end

#The following nonsense means that stylish arrays can override getindex or they can override copyto!(view)
Base.@propagate_inbounds Base.setindex!(arr::StylishArray, v, I::Integer) = _setindex!(arr, v, I)
Base.@propagate_inbounds Base.setindex!(arr::StylishArray, v, I...) = _setindex!(arr, v, I...)

Base.@propagate_inbounds function _setindex!(arr, v, I...)
    dst = view(arr, I...)
    if ndims(dst) == 0
        copyto!(dst, Styled(ScalarArray(v)))[]
    else
        copyto!(dst, Styled(v))
    end
end

Base.@propagate_inbounds function Base.map(f::F, arr::StylishArray, tail...) where {F}
    src = @_ f.(arr, tail...)
    copyto!(similar(src), src)
end

Base.@propagate_inbounds function Base.foreach(f::F, arr::StylishArray, tail...) where {F}
    src = @_ f(arr, tail...)
    copyto!(NullArray(axes(src)), src)
    return nothing
end

Base.@propagate_inbounds function Base.reduce(op::Op, arr::StylishArray; dims=:, kwargs...) where {Op}
    return _reduce(op, arr, dims, kwargs.data)
end

Base.@propagate_inbounds function _reduce(op::Op, arr::StylishArray, dims, nt::NamedTuple{()}) where {Op}
    res = copy(Reduce(op, dims)(arr))
    return dims == Colon() ? res[] : res
end
Base.@propagate_inbounds function _reduce(op::Op, arr::StylishArray, dims, nt::NamedTuple{(:init,)}) where {Op}
    res = copy(Reduce(op, dims)(nt.init, arr))
    return dims == Colon() ? res[] : res
end

Base.@propagate_inbounds function Base.mapreduce(f::F, op::Op, arr::StylishArray; dims=:, kwargs...) where {F, Op}
    return _mapreduce(f, op, arr, dims, kwargs.data)
end

Base.@propagate_inbounds function _mapreduce(f::F, op::Op, arr::StylishArray, dims, nt::NamedTuple{()}) where {F, Op}
    res = copy(Reduce(op, dims)(@_ f.(arr)))
    return dims == Colon() ? res[] : res
end
Base.@propagate_inbounds function _mapreduce(f::F, op::Op, arr::StylishArray, dims, nt::NamedTuple{(:init,)}) where {F, Op}
    res = copy(Reduce(op, dims)(nt.init, @_ f.(arr)))
    return dims == Colon() ? res[] : res
end

Base.@propagate_inbounds function Base.sum(arr::StylishArray; dims=:, kwargs...)
    return _sum(arr, dims, kwargs.data)
end

Base.@propagate_inbounds function _sum(arr::StylishArray, dims, nt::NamedTuple{()})
    res = copy(Sum(dims)(arr))
    return dims == Colon() ? res[] : res
end
Base.@propagate_inbounds function _sum(arr::StylishArray, dims, nt::NamedTuple{(:init,)})
    res = copy(Sum(dims)(nt.init, arr))
    return dims == Colon() ? res[] : res
end

Base.@propagate_inbounds Base.maximum(arr::StylishArray; kwargs...) = reduce(max, arr; kwargs...)

Base.@propagate_inbounds Base.minimum(arr::StylishArray; kwargs...) = reduce(min, arr; kwargs...)

Base.@propagate_inbounds function LinearAlgebra.dot(x::StylishArray, y::AbstractArray)
    res = copy(Sum()(@_ x .* y))
    return ndims(res) == 0 ? res[] : res
end

Base.@propagate_inbounds function LinearAlgebra.dot(x::AbstractArray, y::StylishArray)
    res = copy(Sum()(@_ x .* y))
    return ndims(res) == 0 ? res[] : res
end

Base.@propagate_inbounds function LinearAlgebra.dot(x::StylishArray, y::StylishArray)
    res = copy(Sum()(@_ x .* y))
    return ndims(res) == 0 ? res[] : res
end

Base.@propagate_inbounds function LinearAlgebra.adjoint(x::StylishArray)
    return Beam(2, 1)(@_ adjoint.(x))
end

Base.@propagate_inbounds function LinearAlgebra.transpose(x::StylishArray)
    return Beam(2, 1)(@_ transpose.(x))
end

Base.@propagate_inbounds function Base.permutedims(x::StylishArray)
    return Beam(2, 1)(x)
end

Base.@propagate_inbounds function Base.permutedims(x::StylishArray, perm)
    return Focus(perm...)(x)
end

Ts = (:(StylishArray{<:Any, 1}), :(StylishArray{<:Any, 2}), :(AbstractArray{<:Any, 1}), :(AbstractArray{<:Any, 2}))
for (A, B) in Iterators.product(Ts, Ts)
    if :(StylishArray{<:Any, 1}) in (A, B) || :(StylishArray{<:Any, 2}) in (A, B)
        @eval begin
            Base.@propagate_inbounds function Base.:*(A::$A, B::$B)
                res = copy(SumOut(2)(@_ A .* Beam(2, 3)(B)))
                return ndims(res) == 0 ? res[] : res
            end
        end
    end
end

for (Y, A, B) in Iterators.product(Ts, Ts, Ts)
    if :(StylishArray{<:Any, 1}) in (Y, A, B) || :(StylishArray{<:Any, 2}) in (Y, A, B)
        @eval begin
            Base.@propagate_inbounds function LinearAlgebra.mul!(Y::$Y, A::$A, B::$B)
                res = copyto!(Y, SumOut(2)(@_ A .* Beam(2, 3)(B)))
                return ndims(res) == 0 ? res[] : res
            end
        end
    end
end

struct Square{T, S} <: Number
    arg::T
    scale::S
end

@inline square(x) = Square(sign(x)^2, norm(x))

@inline root(x::Square) = sqrt(x.arg) * x.scale

@inline Base.zero(::Type{Square{T, S}}) where {T, S} = Square{T, S}(zero(T), zero(S))

function Base.promote_rule(::Type{Square{T1, S1}}, ::Type{Square{T2, S2}}) where {T1, S1, T2, S2}
    return Square{promote_type(T1, T2), promote_type(S1, S2)}
end

function Base.convert(::Type{Square{T, S}}, x::Square) where {T, S}
    return Square(convert(T, x.arg), convert(S, x.scale))
end

@inline function Base.:+(x::T, y::T) where {T <: Square}
    if x.scale < y.scale
        (x, y) = (y, x)
    end
    if x.scale > y.scale
        if iszero(y.scale)
            return Square(x.arg + zero(y.arg) * (one(y.scale)/one(x.scale))^1, x.scale)
        else
            return Square(x.arg + y.arg * (y.scale/x.scale)^2, x.scale)
        end
    else
        return Square(x.arg + y.arg * (one(y.scale)/one(x.scale))^1, x.scale)
    end
end

struct Power{T, S, E} <: Number
    arg::T
    scale::S
    exponent::E
end

@inline power(x, p) = Power(sign(x)^p, norm(x), p)

@inline root(x::Power) = x.arg ^ inv(x.exponent) * x.scale

@inline Base.zero(::Type{Power{T, S, E}}) where {T, S, E} = Power{T, S, E}(zero(T), zero(S), one(E))
@inline Base.zero(x::Power) = Power(zero(x.arg), zero(x.scale), x.exponent)

function Base.promote_rule(::Type{Power{T1, S1, E1}}, ::Type{Power{T2, S2, E2}}) where {T1, S1, E1, T2, S2, E2}
    return Power{promote_type(T1, T2), promote_type(S1, S2), promote_type(E1, E2)}
end

function Base.convert(::Type{Power{T, S, E}}, x::Power) where {T, S, E}
    return Power(convert(T, x.arg), convert(S, x.scale), convert(E, x.exponent))
end

@inline function Base.:+(x::T, y::T) where {T <: Power}
    if x.exponent != y.exponent
        if iszero(x.arg) && iszero(x.scale)
            (x, y) = (y, x)
        end
        if iszero(y.arg) && iszero(y.scale)
            y = Power(y.arg, y.scale, x.exponent)
        else
            ArgumentError("Cannot accurately add Powers with different exponents")
        end
    end
    #TODO handle negative exponent
    if x.scale < y.scale
        (x, y) = (y, x)
    end
    if x.scale > y.scale
        if iszero(y.scale)
            return Power(x.arg + zero(y.arg) * (one(y.scale)/one(x.scale))^one(y.exponent), x.scale, x.exponent)
        else
            return Power(x.arg + y.arg * (y.scale/x.scale)^y.exponent, x.scale, x.exponent)
        end
    else
        return Power(x.arg + y.arg * (one(y.scale)/one(x.scale))^one(y.exponent), x.scale, x.exponent)
    end
end

Base.@propagate_inbounds function LinearAlgebra.norm(arr::StylishArray; kwargs...)
    norm(arr, 2; kwargs...)
end

Base.@propagate_inbounds function LinearAlgebra.norm(arr::StylishArray, p::Real)
    if p == 2
        return root(copy(Sum()(@_ square.(arr)))[])
    elseif p == 1
        return copy(Sum()(@_ norm.(arr)))[]
    elseif p == Inf
        return copy(Reduce(max)(@_ norm.(arr)))[]
    elseif p == 0
        return copy(Sum()(@_ norm.(norm.(arr), 0)))[]
    elseif p == -Inf
        return copy(Reduce(min)(@_ norm.(arr)))[]
    else
        return root(copy(Sum()(@_ power.(arr, p)))[])
    end
end

Base.@propagate_inbounds function distance(x, y)
    return distance(x, y, 2)
end
Base.@propagate_inbounds function distance(x, y, p)
    return norm(arrayify(@_ x .- y), p)
end

end


text:
using Oceananigans
using Oceananigans.Units
using Statistics
using Printf
using Oceananigans.MultiRegion
using LinearAlgebra, SparseArrays
using Oceananigans.Solvers: constructors, unpack_constructors
using Oceananigans.Grids: architecture
using Oceananigans.Utils

function geostrophic_adjustment_test(free_surface, grid; regions = 1)

    if architecture(grid) isa GPU
        devices = (0, 0)
    else
        devices = nothing
    end
    mrg = MultiRegionGrid(grid, partition = XPartition(regions), devices = devices)

    coriolis = FPlane(f = 1e-4)

    model = HydrostaticFreeSurfaceModel(grid = mrg,
                                        coriolis = coriolis,
                                        free_surface = free_surface)

    gaussian(x, L) = exp(-x^2 / 2L^2)

    U = 0.1 # geostrophic velocity
    L  = grid.Lx / 40 # gaussian width
    x₀ = grid.Lx / 4 # gaussian center

    vᵍ(x, y, z) = -U * (x - x₀) / L * gaussian(x - x₀, L)

    g = model.free_surface.gravitational_acceleration
    η = model.free_surface.η

    η₀ = coriolis.f * U * L / g # geostrohpic free surface amplitude

    ηᵍ(x) = η₀ * gaussian(x - x₀, L)

    ηⁱ(x, y) = 2 * ηᵍ(x)

    set!(model, v = vᵍ)
    @apply_regionally set!(η, ηⁱ)

    gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed
    Δt = 2 * model.grid.Δxᶜᵃᵃ / gravity_wave_speed

    for step in 1:10
        time_step!(model, Δt)
    end

    return η
end

Lh = 100kilometers
Lz = 400meters

for arch in archs

    free_surface   = ImplicitFreeSurface(solver_method = :HeptadiagonalIterativeSolver, maximum_iterations = 64 * 3)
    topology_types = (Bounded, Periodic, Bounded), (Periodic, Periodic, Bounded)

    @testset "Testing multi region implicit free surface" begin
        for topology_type in topology_types
            grid = RectilinearGrid(arch,
                        size = (64, 3, 1),
                        x = (0, Lh), y = (0, Lh), z = (-Lz, 0),
                        topology = topology_type)

            ηs = geostrophic_adjustment_test(free_surface, grid);
            ηs = Array(interior(ηs));

            for regions in [2, 4]
                @info "  Testing $regions partitions on $(topology_type) on the $arch"
                η = geostrophic_adjustment_test(free_surface, grid, regions = regions)
                η = Array(interior(reconstruct_global_field(η)))

                @test all(η .≈ ηs)
            end
        end
    end
end


text:
export CuVec
module CuVec
using ..Mocha

const THREADS_PER_BLOCK = 128
function cuda_geometry(len::Int)
  x_block = round(Int, ceil(convert(Float64, len)/THREADS_PER_BLOCK))
  return (x_block, THREADS_PER_BLOCK)
end

for (ctype, dtype) in [(:float, Float32), (:double, Float64)]
  # define add!, sub!, mul!, div!, div2!
  for name in [:add, :sub, :mul, :div, :div2]
    @eval begin
      function $(Symbol("$(name)!"))(backend::GPUBackend, ::Type{$dtype}, X, Y, len::Int)
        X = convert(Ptr{Void},X)
        Y = convert(Ptr{Void},Y)
        cuda_dim = cuda_geometry(len)
        kernel = backend.mocha.$(Symbol("elem_$(name)_$ctype"))
        CUDA.launch(kernel, cuda_dim..., (X, Y, len))
      end
    end
  end

  # define add_scal!
  @eval begin
    function add_scal!(backend::GPUBackend, ::Type{$dtype}, X, Y, len::Int)
      X = convert(Ptr{Void}, X)
      Y = convert($dtype, Y)
      cuda_dim = cuda_geometry(len)
      kernel = backend.mocha.$(Symbol("add_scal_$ctype"))
      CUDA.launch(kernel, cuda_dim..., (X,Y,len))
    end
  end

  # define log!
  @eval begin
    function log!(backend::GPUBackend, ::Type{$dtype}, X, len::Int)
      X = convert(Ptr{Void}, X)
      cuda_dim = cuda_geometry(len)
      kernel = backend.mocha.$(Symbol("elem_log_$ctype"))
      CUDA.launch(kernel, cuda_dim..., (X,len))
    end
  end

  # define exp!
  @eval begin
    function exp!(backend::GPUBackend, ::Type{$dtype}, X, len::Int)
      X = convert(Ptr{Void}, X)
      cuda_dim = cuda_geometry(len)
      kernel = backend.mocha.$(Symbol("elem_exp_$ctype"))
      CUDA.launch(kernel, cuda_dim..., (X,len))
    end
  end

  # define mul_scal!
  @eval begin
    function mul_scal!(backend::GPUBackend, ::Type{$dtype}, X, Y, len::Int)
      X = convert(Ptr{Void}, X)
      Y = convert($dtype, Y)
      cuda_dim = cuda_geometry(len)
      kernel = backend.mocha.$(Symbol("mul_scal_$ctype"))
      CUDA.launch(kernel, cuda_dim..., (X,Y,len))
    end
  end
end

# define add!, sub!, mul!, div!, div2! for blobs
for name in [:add, :sub, :mul, :div, :div2]
  @eval begin
    function $(Symbol("$(name)!")){T}(backend::GPUBackend, X::CuTensorBlob{T}, Y::CuTensorBlob{T})
      len = length(X)
      $(Symbol("$(name)!"))(backend, T, X.ptr.p, Y.ptr.p, len)
    end
  end
end
function add_scal!{T}(backend::GPUBackend, X::CuTensorBlob{T}, Y)
  Y = convert(T, Y)
  len = length(X)
  add_scal!(backend, T, X.ptr.p, Y, len)
end
function mul_scal!{T}(backend::GPUBackend, X::CuTensorBlob{T}, Y)
  Y = convert(T, Y)
  len = length(X)
  mul_scal!(backend, T, X.ptr.p, Y, len)
end
function log!{T}(backend::GPUBackend, X::CuTensorBlob{T})
  log!(backend, T, X.ptr.p, length(X))
end
function exp!{T}(backend::GPUBackend, X::CuTensorBlob{T})
  exp!(backend, T, X.ptr.p, length(X))
end

for (postfix, dt1, dt2) in [(:fi, Float32, Int), (:di, Float64, Int),
                            (:ff, Float32, Float32), (:dd, Float64, Float64)]
  @eval begin
    function pow!(backend::GPUBackend, ::Type{$dt1}, X, Y::$dt2, len::Int)
      X = convert(Ptr{Void}, X)
      cuda_dim = cuda_geometry(len)
      kernel = backend.mocha.$(Symbol("elem_pow_$postfix"))
      CUDA.launch(kernel, cuda_dim..., (X,Y,len))
    end
    function pow!{T}(backend::GPUBackend, X::CuTensorBlob{T}, Y::$dt2)
      pow!(backend, T, X.ptr.p, Y, length(X))
    end
  end
end

end # module cuVec