text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(tanh(a+b*x)^6, x) == :(x+-1*b^-1*tanh(a+b*x)+-1//3*b^-1*tanh(a+b*x)^3+-1//5*b^-1*tanh(a+b*x)^5)
@test integrate(tanh(a+b*x)^5, x) == :(b^-1*log(cosh(a+b*x))+-1//2*b^-1*tanh(a+b*x)^2+-1//4*b^-1*tanh(a+b*x)^4)
@test integrate(tanh(a+b*x)^4, x) == :(x+-1*b^-1*tanh(a+b*x)+-1//3*b^-1*tanh(a+b*x)^3)
@test integrate(tanh(a+b*x)^3, x) == :(b^-1*log(cosh(a+b*x))+-1//2*b^-1*tanh(a+b*x)^2)
@test integrate(tanh(a+b*x)^2, x) == :(x+-1*b^-1*tanh(a+b*x))
@test integrate(tanh(a+b*x), x) == :(b^-1*log(cosh(a+b*x)))
@test integrate(coth(a+b*x), x) == :(b^-1*log(sinh(a+b*x)))
@test integrate(coth(a+b*x)^2, x) == :(x+-1*b^-1*coth(a+b*x))
@test integrate(coth(a+b*x)^3, x) == :(b^-1*log(sinh(a+b*x))+-1//2*b^-1*coth(a+b*x)^2)
@test integrate(coth(a+b*x)^4, x) == :(x+-1*b^-1*coth(a+b*x)+-1//3*b^-1*coth(a+b*x)^3)
@test integrate(coth(a+b*x)^5, x) == :(b^-1*log(sinh(a+b*x))+-1//2*b^-1*coth(a+b*x)^2+-1//4*b^-1*coth(a+b*x)^4)
@test integrate(coth(a+b*x)^6, x) == :(x+-1*b^-1*coth(a+b*x)+-1//3*b^-1*coth(a+b*x)^3+-1//5*b^-1*coth(a+b*x)^5)
@test integrate((b*tanh(c+d*x))^7//2, x) == :(b^7//2*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2))+b^7//2*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2))+-2*b^3*d^-1*(b*tanh(c+d*x))^(1/2)+-2//5*b*d^-1*(b*tanh(c+d*x))^5//2)
@test integrate((b*tanh(c+d*x))^5//2, x) == :(b^5//2*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2))+-1*b^5//2*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2))+-2//3*b*d^-1*(b*tanh(c+d*x))^3//2)
@test integrate((b*tanh(c+d*x))^3//2, x) == :(b^3//2*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2))+b^3//2*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2))+-2*b*d^-1*(b*tanh(c+d*x))^(1/2))
@test integrate((b*tanh(c+d*x))^(1/2), x) == :(b^(1/2)*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2))+-1*b^(1/2)*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2)))
@test integrate((b*tanh(c+d*x))^-1//2, x) == :(b^-1//2*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2))+b^-1//2*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2)))
@test integrate((b*tanh(c+d*x))^-3//2, x) == :(b^-3//2*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2))+-1*b^-3//2*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2))+-2*b^-1*d^-1*(b*tanh(c+d*x))^-1//2)
@test integrate((b*tanh(c+d*x))^-5//2, x) == :(b^-5//2*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2))+b^-5//2*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2))+-2//3*b^-1*d^-1*(b*tanh(c+d*x))^-3//2)
@test integrate((b*tanh(c+d*x))^-7//2, x) == :(b^-7//2*d^-1*arctanh(b^-1//2*(b*tanh(c+d*x))^(1/2))+-1*b^-7//2*d^-1*arctan(b^-1//2*(b*tanh(c+d*x))^(1/2))+-2*b^-3*d^-1*(b*tanh(c+d*x))^-1//2+-2//5*b^-1*d^-1*(b*tanh(c+d*x))^-5//2)
@test integrate(tanh(8x)^1//3, x) == :(-1//16*log(1+-1*tanh(8x)^2//3)+1//32*log(1+tanh(8x)^2//3+tanh(8x)^4//3)+-1//16*3^(1/2)*arctan(1//3*3^(1/2)*(1+2*tanh(8x)^2//3)))
@test integrate((a*tanh(x)^2)^3//2, x) == :(-1//2*a*(a*tanh(x)^2)^(1/2)*tanh(x)+a*(a*tanh(x)^2)^(1/2)*coth(x)*log(cosh(x)))
@test integrate((a*tanh(x)^2)^(1/2), x) == :((a*tanh(x)^2)^(1/2)*coth(x)*log(cosh(x)))
@test integrate((a*tanh(x)^2)^-1//2, x) == :((a*tanh(x)^2)^-1//2*log(sinh(x))*tanh(x))
@test integrate((-1*tanh(c+d*x)^2)^5//2, x) == :(-1//2*d^-1*(-1*tanh(c+d*x)^2)^(1/2)*tanh(c+d*x)+-1//4*d^-1*(-1*tanh(c+d*x)^2)^(1/2)*tanh(c+d*x)^3+d^-1*(-1*tanh(c+d*x)^2)^(1/2)*coth(c+d*x)*log(cosh(c+d*x)))
@test integrate((-1*tanh(c+d*x)^2)^3//2, x) == :((1/2)*d^-1*(-1*tanh(c+d*x)^2)^(1/2)*tanh(c+d*x)+-1*d^-1*(-1*tanh(c+d*x)^2)^(1/2)*coth(c+d*x)*log(cosh(c+d*x)))
@test integrate((-1*tanh(c+d*x)^2)^(1/2), x) == :(d^-1*(-1*tanh(c+d*x)^2)^(1/2)*coth(c+d*x)*log(cosh(c+d*x)))
@test integrate((-1*tanh(c+d*x)^2)^-1//2, x) == :(d^-1*(-1*tanh(c+d*x)^2)^-1//2*log(sinh(c+d*x))*tanh(c+d*x))
@test integrate((-1*tanh(c+d*x)^2)^-3//2, x) == :((1/2)*d^-1*(-1*tanh(c+d*x)^2)^-1//2*coth(c+d*x)+-1*d^-1*(-1*tanh(c+d*x)^2)^-1//2*log(sinh(c+d*x))*tanh(c+d*x))
@test integrate((-1*tanh(c+d*x)^2)^-5//2, x) == :(-1//2*d^-1*(-1*tanh(c+d*x)^2)^-1//2*coth(c+d*x)+-1//4*d^-1*(-1*tanh(c+d*x)^2)^-1//2*coth(c+d*x)^3+d^-1*(-1*tanh(c+d*x)^2)^-1//2*log(sinh(c+d*x))*tanh(c+d*x))
@test integrate((tanh(x)^3)^(1/2), x) == :(-2*(tanh(x)^3)^(1/2)*coth(x)+(tanh(x)^3)^(1/2)*tanh(x)^-3//2*arctan(tanh(x)^(1/2))+(tanh(x)^3)^(1/2)*tanh(x)^-3//2*arctanh(tanh(x)^(1/2)))
@test integrate((a*tanh(x)^3)^3//2, x) == :(-2//3*a*(a*tanh(x)^3)^(1/2)+-2//7*a*(a*tanh(x)^3)^(1/2)*tanh(x)^2+a*(a*tanh(x)^3)^(1/2)*tanh(x)^-3//2*arctanh(tanh(x)^(1/2))+-1*a*(a*tanh(x)^3)^(1/2)*tanh(x)^-3//2*arctan(tanh(x)^(1/2)))
@test integrate((a*tanh(x)^3)^(1/2), x) == :(-2*(a*tanh(x)^3)^(1/2)*coth(x)+(a*tanh(x)^3)^(1/2)*tanh(x)^-3//2*arctan(tanh(x)^(1/2))+(a*tanh(x)^3)^(1/2)*tanh(x)^-3//2*arctanh(tanh(x)^(1/2)))
@test integrate((a*tanh(x)^3)^-1//2, x) == :(-2*(a*tanh(x)^3)^-1//2*tanh(x)+(a*tanh(x)^3)^-1//2*tanh(x)^3//2*arctanh(tanh(x)^(1/2))+-1*(a*tanh(x)^3)^-1//2*tanh(x)^3//2*arctan(tanh(x)^(1/2)))
@test integrate((a*tanh(x)^4)^3//2, x) == :(-1*a*(a*tanh(x)^4)^(1/2)*coth(x)+-1//3*a*(a*tanh(x)^4)^(1/2)*tanh(x)+-1//5*a*(a*tanh(x)^4)^(1/2)*tanh(x)^3+a*x*(a*tanh(x)^4)^(1/2)*coth(x)^2)
@test integrate((a*tanh(x)^4)^(1/2), x) == :(-1*(a*tanh(x)^4)^(1/2)*coth(x)+x*(a*tanh(x)^4)^(1/2)*coth(x)^2)
@test integrate((a*tanh(x)^4)^-1//2, x) == :(-1*(a*tanh(x)^4)^-1//2*tanh(x)+x*(a*tanh(x)^4)^-1//2*tanh(x)^2)
@test integrate((a+a*tanh(c+d*x))^5, x) == :(16*x*a^5+-8*a^5*d^-1*tanh(c+d*x)+-2*a*d^-1*(a^2+a^2*tanh(c+d*x))^2+16*a^5*d^-1*log(cosh(c+d*x))+-2//3*a^2*d^-1*(a+a*tanh(c+d*x))^3+-1//4*a*d^-1*(a+a*tanh(c+d*x))^4)
@test integrate((a+a*tanh(c+d*x))^4, x) == :(-1*d^-1*(a^2+a^2*tanh(c+d*x))^2+8*x*a^4+-4*a^4*d^-1*tanh(c+d*x)+8*a^4*d^-1*log(cosh(c+d*x))+-1//3*a*d^-1*(a+a*tanh(c+d*x))^3)
@test integrate((a+a*tanh(c+d*x))^3, x) == :(4*x*a^3+-2*a^3*d^-1*tanh(c+d*x)+4*a^3*d^-1*log(cosh(c+d*x))+-1//2*a*d^-1*(a+a*tanh(c+d*x))^2)
@test integrate((a+a*tanh(c+d*x))^2, x) == :(2*x*a^2+-1*a^2*d^-1*tanh(c+d*x)+2*a^2*d^-1*log(cosh(c+d*x)))
@test integrate((a+a*tanh(c+d*x))^-1, x) == :((1/2)*x*a^-1+-1//2*d^-1*(a+a*tanh(c+d*x))^-1)
@test integrate((a+a*tanh(c+d*x))^-2, x) == :(-1//4*d^-1*(a+a*tanh(c+d*x))^-2+-1//4*d^-1*(a^2+a^2*tanh(c+d*x))^-1+1//4*x*a^-2)
@test integrate((a+a*tanh(c+d*x))^-3, x) == :(-1//6*d^-1*(a+a*tanh(c+d*x))^-3+-1//8*d^-1*(a^3+a^3*tanh(c+d*x))^-1+1//8*x*a^-3+-1//8*a^-1*d^-1*(a+a*tanh(c+d*x))^-2)
@test integrate((a+a*tanh(c+d*x))^-4, x) == :(-1//8*d^-1*(a+a*tanh(c+d*x))^-4+-1//16*d^-1*(a^2+a^2*tanh(c+d*x))^-2+-1//16*d^-1*(a^4+a^4*tanh(c+d*x))^-1+1//16*x*a^-4+-1//12*a^-1*d^-1*(a+a*tanh(c+d*x))^-3)
@test integrate((a+a*tanh(c+d*x))^-5, x) == :(-1//10*d^-1*(a+a*tanh(c+d*x))^-5+-1//32*d^-1*(a^5+a^5*tanh(c+d*x))^-1+1//32*x*a^-5+-1//16*a^-1*d^-1*(a+a*tanh(c+d*x))^-4+-1//24*a^-2*d^-1*(a+a*tanh(c+d*x))^-3+-1//32*a^-1*d^-1*(a^2+a^2*tanh(c+d*x))^-2)
@test integrate((1+tanh(x))^7//2, x) == :(-8*(1+tanh(x))^(1/2)+-4//3*(1+tanh(x))^3//2+-2//5*(1+tanh(x))^5//2+8*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^5//2, x) == :(-4*(1+tanh(x))^(1/2)+-2//3*(1+tanh(x))^3//2+4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^3//2, x) == :(-2*(1+tanh(x))^(1/2)+2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^(1/2), x) == :(2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^-1//2, x) == :(-1*(1+tanh(x))^-1//2+(1/2)*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^-3//2, x) == :(-1//2*(1+tanh(x))^-1//2+-1//3*(1+tanh(x))^-3//2+1//4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^-5//2, x) == :(-1//4*(1+tanh(x))^-1//2+-1//5*(1+tanh(x))^-5//2+-1//6*(1+tanh(x))^-3//2+1//8*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((a+b*tanh(c+d*x))^5, x) == :(a*x*(a^4+5*b^4+10*a^2*b^2)+-1//4*b*d^-1*(a+b*tanh(c+d*x))^4+b*d^-1*(b^4+5*a^4+10*a^2*b^2)*log(cosh(c+d*x))+-2//3*a*b*d^-1*(a+b*tanh(c+d*x))^3+-1//2*b*d^-1*(a+b*tanh(c+d*x))^2*(b^2+3*a^2)+-4*a*b^2*d^-1*(a^2+b^2)*tanh(c+d*x))
@test integrate((a+b*tanh(c+d*x))^4, x) == :(x*(a^4+b^4+6*a^2*b^2)+-1//3*b*d^-1*(a+b*tanh(c+d*x))^3+-1*a*b*d^-1*(a+b*tanh(c+d*x))^2+-1*b^2*d^-1*(b^2+3*a^2)*tanh(c+d*x)+4*a*b*d^-1*(a^2+b^2)*log(cosh(c+d*x)))
@test integrate((a+b*tanh(c+d*x))^3, x) == :(a*x*(a^2+3*b^2)+-1//2*b*d^-1*(a+b*tanh(c+d*x))^2+b*d^-1*(b^2+3*a^2)*log(cosh(c+d*x))+-2*a*b^2*d^-1*tanh(c+d*x))
@test integrate((a+b*tanh(c+d*x))^2, x) == :(x*(a^2+b^2)+-1*b^2*d^-1*tanh(c+d*x)+2*a*b*d^-1*log(cosh(c+d*x)))
@test integrate((a+b*tanh(c+d*x))^-1, x) == :(a*x*(a^2+-1*b^2)^-1+-1*b*d^-1*(a^2+-1*b^2)^-1*log(a*cosh(c+d*x)+b*sinh(c+d*x)))
@test integrate((a+b*tanh(c+d*x))^-2, x) == :(x*(a^2+-1*b^2)^-2*(a^2+b^2)+b*d^-1*(a+b*tanh(c+d*x))^-1*(a^2+-1*b^2)^-1+-2*a*b*d^-1*(a^2+-1*b^2)^-2*log(a*cosh(c+d*x)+b*sinh(c+d*x)))
@test integrate((a+b*tanh(c+d*x))^-3, x) == :(a*x*(a^2+-1*b^2)^-3*(a^2+3*b^2)+b*d^-1*(a+b*tanh(c+d*x))^-2*(-2*b^2+2*a^2)^-1+-1*b*d^-1*(a^2+-1*b^2)^-3*(b^2+3*a^2)*log(a*cosh(c+d*x)+b*sinh(c+d*x))+2*a*b*d^-1*(a+b*tanh(c+d*x))^-1*(a^2+-1*b^2)^-2)
@test integrate((a+b*tanh(c+d*x))^-4, x) == :(x*(a^2+-1*b^2)^-4*(a^4+b^4+6*a^2*b^2)+b*d^-1*(a+b*tanh(c+d*x))^-3*(-3*b^2+3*a^2)^-1+a*b*d^-1*(a+b*tanh(c+d*x))^-2*(a^2+-1*b^2)^-2+b*d^-1*(a+b*tanh(c+d*x))^-1*(a^2+-1*b^2)^-3*(b^2+3*a^2)+-4*a*b*d^-1*(a^2+-1*b^2)^-4*(a^2+b^2)*log(a*cosh(c+d*x)+b*sinh(c+d*x)))
@test integrate((4+6*tanh(c+d*x))^-1, x) == :(-1//5*x+3//10*d^-1*log(2*cosh(c+d*x)+3*sinh(c+d*x)))
@test integrate((4+-6*tanh(c+d*x))^-1, x) == :(-1//5*x+-3//10*d^-1*log(-3*sinh(c+d*x)+2*cosh(c+d*x)))
@test integrate((a+b*tanh(c+d*x))^(1/2), x) == :(d^-1*(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*tanh(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^(1/2)*arctanh((a+-1b)^-1//2*(a+b*tanh(c+d*x))^(1/2)))
@test integrate((a+b*tanh(c+d*x))^-1//2, x) == :(d^-1*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*tanh(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*tanh(c+d*x))^(1/2)))
@test integrate((1+tanh(x))^-1*sinh(x)^4, x) == :(-1*(8+-8*tanh(x))^-1+-3*(16+16*tanh(x))^-1+-1//24*(1+tanh(x))^-3+1//16*x+1//32*(1+-1*tanh(x))^-2+5//32*(1+tanh(x))^-2)
@test integrate((1+tanh(x))^-1*sinh(x)^3, x) == :(-1//3*cosh(x)^3+-1//5*sinh(x)^5+1//5*cosh(x)^5)
@test integrate((1+tanh(x))^-1*sinh(x)^2, x) == :((4+4*tanh(x))^-1+(8+-8*tanh(x))^-1+-1//8*x+-1//8*(1+tanh(x))^-2)
@test integrate((1+tanh(x))^-1*sinh(x), x) == :(-1//3*sinh(x)^3+1//3*cosh(x)^3)
@test integrate((1+tanh(x))^-1*csch(x), x) == :(-1*arctanh(cosh(x))+-1*sinh(x)+cosh(x))
@test integrate((1+tanh(x))^-1*csch(x)^2, x) == :(-1*coth(x)+-1*log(tanh(x))+log(1+tanh(x)))
@test integrate((1+tanh(x))^-1*csch(x)^3, x) == :(-1//2*arctanh(cosh(x))+-1//2*coth(x)*csch(x)+csch(x))
@test integrate((1+tanh(x))^-1*csch(x)^4, x) == :((1/2)*coth(x)^2+-1//3*coth(x)^3)
@test integrate((1+tanh(x))^-1*csch(x)^5, x) == :(1//3*csch(x)^3+1//8*arctanh(cosh(x))+-1//4*csch(x)^3*coth(x)+-1//8*coth(x)*csch(x))
@test integrate((1+tanh(x))^-1*csch(x)^6, x) == :(-1//2*coth(x)^2+-1//5*coth(x)^5+1//3*coth(x)^3+1//4*coth(x)^4)
@test integrate((1+tanh(x))^-1*csch(x)^7, x) == :(-1//16*arctanh(cosh(x))+1//5*csch(x)^5+-1//6*csch(x)^5*coth(x)+-1//24*csch(x)^3*coth(x)+1//16*coth(x)*csch(x))
@test integrate((a+b*tanh(x))^-1*sinh(x)^4, x) == :(-1*(-4*b^2+4*a^2)^-1*cosh(x)^4*(b+-1*a*tanh(x))+1//8*(a^2+-1*b^2)^-2*cosh(x)^2*(4*b*(-1*b^2+2*a^2)+-1*a*(-1*b^2+5*a^2)*tanh(x))+-1*b*a^4*(a^2+-1*b^2)^-3*log(a+b*tanh(x))+-1//16*a*(a+b)^-3*(b+3a)*log(1+-1*tanh(x))+1//16*a*(a+-1b)^-3*(-1b+3a)*log(1+tanh(x)))
@test integrate((a+b*tanh(x))^-1*sinh(x)^3, x) == :(a*(-3*b^2+3*a^2)^-1*cosh(x)^3+-1*a*(a^2+-1*b^2)^-1*cosh(x)+-1*b*(-3*b^2+3*a^2)^-1*sinh(x)^3+b*a^2*(a^2+-1*b^2)^-2*sinh(x)+-1*a*b^2*(a^2+-1*b^2)^-2*cosh(x)+-1*b*a^3*(a^2+-1*b^2)^-5//2*arctan((a^2+-1*b^2)^-1//2*(a*sinh(x)+b*cosh(x))))
@test integrate((a+b*tanh(x))^-1*sinh(x)^2, x) == :(-1*(-2*b^2+2*a^2)^-1*cosh(x)^2*(b+-1*a*tanh(x))+-1//4*a*(a+-1b)^-2*log(1+tanh(x))+1//4*a*(a+b)^-2*log(1+-1*tanh(x))+b*a^2*(a^2+-1*b^2)^-2*log(a+b*tanh(x)))
@test integrate((a+b*tanh(x))^-1*sinh(x), x) == :(a*(a^2+-1*b^2)^-1*cosh(x)+-1*b*(a^2+-1*b^2)^-1*sinh(x)+a*b*(a^2+-1*b^2)^-3//2*arctan((a^2+-1*b^2)^-1//2*(a*sinh(x)+b*cosh(x))))
@test integrate((a+b*tanh(x))^-1*csch(x), x) == :(-1*a^-1*arctanh(cosh(x))+-1*b*a^-1*(a^2+-1*b^2)^-1//2*arctan((a^2+-1*b^2)^-1//2*(a*sinh(x)+b*cosh(x))))
@test integrate((a+b*tanh(x))^-1*csch(x)^2, x) == :(-1*a^-1*coth(x)+b*a^-2*log(a+b*tanh(x))+-1*b*a^-2*log(tanh(x)))
@test integrate((a+b*tanh(x))^-1*csch(x)^3, x) == :((1/2)*a^-1*arctanh(cosh(x))+b*a^-2*csch(x)+-1*a^-3*b^2*arctanh(cosh(x))+-1//2*a^-1*coth(x)*csch(x)+b*a^-3*(a^2+-1*b^2)^(1/2)*arctan((a^2+-1*b^2)^-1//2*(a*sinh(x)+b*cosh(x))))
@test integrate((a+b*tanh(x))^-1*csch(x)^4, x) == :(-1//3*a^-1*coth(x)^3+a^-3*(a^2+-1*b^2)*coth(x)+(1/2)*b*a^-2*coth(x)^2+b*a^-4*(a^2+-1*b^2)*log(tanh(x))+-1*b*a^-4*(a^2+-1*b^2)*log(a+b*tanh(x)))
@test integrate((a+b*tanh(x))^-1*csch(x)^5, x) == :(-3//8*a^-1*arctanh(cosh(x))+a^-5*b^4*sech(x)+a^-4*b^3*arctan(sinh(x))+-1*b*a^-2*arctan(sinh(x))+-1*b*a^-2*csch(x)+-1*a^-5*b^4*arctanh(cosh(x))+-3//2*a^-3*b^2*sech(x)+-1//4*a^-1*csch(x)^3*coth(x)+1//3*b*a^-2*csch(x)^3+3//2*a^-4*b^3*csch(x)+3//2*a^-3*b^2*arctanh(cosh(x))+3//8*a^-1*coth(x)*csch(x)+b*a^-4*(a^2+-1*b^2)*arctan(sinh(x))+a^-5*b^2*(a^2+-1*b^2)*sech(x)+-1*b*a^-5*(a^2+-1*b^2)^3//2*arctan((a^2+-1*b^2)^-1//2*(a*sinh(x)+b*cosh(x)))+-1//2*a^-4*b^3*sech(x)^2*csch(x)+-1//2*a^-4*b^3*sech(x)*tanh(x)+-1//2*a^-3*b^2*csch(x)^2*sech(x))
@test integrate((a+b*tanh(x))^-1*csch(x)^6, x) == :(-1//5*a^-1*coth(x)^5+-1*a^-5*(a^2+-1*b^2)^2*coth(x)+1//3*a^-3*coth(x)^3*(-1*b^2+2*a^2)+1//4*b*a^-2*coth(x)^4+b*a^-6*(a^2+-1*b^2)^2*log(a+b*tanh(x))+-1*b*a^-6*(a^2+-1*b^2)^2*log(tanh(x))+-1//2*b*a^-4*coth(x)^2*(-1*b^2+2*a^2))
@test integrate((im+tanh(x))^-1*csch(x), x) == :(im*arctanh(cosh(x))+-1//2*im*2^(1/2)*arctanh((1/2)*2^(1/2)*(im*sinh(x)+cosh(x))))
@test integrate((1+tanh(x))^-1*cosh(x)^4, x) == :((8+-8*tanh(x))^-1+-3*(16+16*tanh(x))^-1+-3//32*(1+tanh(x))^-2+-1//24*(1+tanh(x))^-3+1//32*(1+-1*tanh(x))^-2+5//16*x)
@test integrate((1+tanh(x))^-1*cosh(x)^3, x) == :(4//5*sinh(x)+4//15*sinh(x)^3+-1*(5+5*tanh(x))^-1*cosh(x)^3)
@test integrate((1+tanh(x))^-1*cosh(x)^2, x) == :((8+-8*tanh(x))^-1+-1*(4+4*tanh(x))^-1+-1//8*(1+tanh(x))^-2+3//8*x)
@test integrate((1+tanh(x))^-1*cosh(x), x) == :(2//3*sinh(x)+-1*(3+3*tanh(x))^-1*cosh(x))
@test integrate((1+tanh(x))^-1*sech(x), x) == :(-1*(1+tanh(x))^-1*sech(x))
@test integrate((1+tanh(x))^-1*sech(x)^2, x) == :(log(1+tanh(x)))
@test integrate((1+tanh(x))^-1*sech(x)^3, x) == :(arctan(sinh(x))+sech(x))
@test integrate((1+tanh(x))^-1*sech(x)^4, x) == :(-1//2*tanh(x)^2+tanh(x))
@test integrate((1+tanh(x))^-1*sech(x)^5, x) == :((1/2)*arctan(sinh(x))+1//3*sech(x)^3+(1/2)*sech(x)*tanh(x))
@test integrate((1+tanh(x))^-1*sech(x)^6, x) == :(-2//3*(1+-1*tanh(x))^3+1//4*(1+-1*tanh(x))^4)
@test integrate((1+tanh(x))^-1*sech(x)^7, x) == :(1//5*sech(x)^5+3//8*arctan(sinh(x))+1//4*sech(x)^3*tanh(x)+3//8*sech(x)*tanh(x))
@test integrate((a+b*tanh(x))^-1*sech(x)^8, x) == :(-1//6*b^-1*tanh(x)^6+-1*b^-7*(a^2+-1*b^2)^3*log(a+b*tanh(x))+-1//2*b^-5*tanh(x)^2*(a^4+3*b^4+-3*a^2*b^2)+-1//4*b^-3*tanh(x)^4*(a^2+-3*b^2)+1//5*a*b^-2*tanh(x)^5+a*b^-6*(a^4+3*b^4+-3*a^2*b^2)*tanh(x)+1//3*a*b^-4*tanh(x)^3*(a^2+-3*b^2))
@test integrate((a+b*tanh(x))^-1*sech(x)^6, x) == :(1//4*b^-1*tanh(x)^4+b^-5*(a^2+-1*b^2)^2*log(a+b*tanh(x))+(1/2)*b^-3*tanh(x)^2*(a^2+-2*b^2)+-1//3*a*b^-2*tanh(x)^3+-1*a*b^-4*(a^2+-2*b^2)*tanh(x))
@test integrate((a+b*tanh(x))^-1*sech(x)^4, x) == :(-1//2*b^-1*tanh(x)^2+a*b^-2*tanh(x)+-1*b^-3*(a^2+-1*b^2)*log(a+b*tanh(x)))
@test integrate((a+b*tanh(x))^-1*sech(x)^2, x) == :(b^-1*log(a+b*tanh(x)))
@test integrate((a+b*tanh(x))^-1, x) == :(a*x*(a^2+-1*b^2)^-1+-1*b*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x)))
@test integrate((a+b*tanh(x))^-1*cosh(x)^2, x) == :(b^3*(a^2+-1*b^2)^-2*log(a+b*tanh(x))+-1*(-2*b^2+2*a^2)^-1*cosh(x)^2*(b+-1*a*tanh(x))+-1//4*(a+b)^-2*(a+2b)*log(1+-1*tanh(x))+1//4*(a+-1b)^-2*(a+-2b)*log(1+tanh(x)))
@test integrate((a+b*tanh(x))^-1*cosh(x)^4, x) == :(-1*b^5*(a^2+-1*b^2)^-3*log(a+b*tanh(x))+-1*(-4*b^2+4*a^2)^-1*cosh(x)^4*(b+-1*a*tanh(x))+-1//16*(a+b)^-3*(3*a^2+8*b^2+9*a*b)*log(1+-1*tanh(x))+1//8*(a^2+-1*b^2)^-2*cosh(x)^2*(4*b^3+-1*a*b^2*(7+-3*a^2*b^-2)*tanh(x))+1//16*(a+-1b)^-3*(3*a^2+8*b^2+-9*a*b)*log(1+tanh(x)))
@test integrate((a+b*tanh(x))^-1*sech(x)^3, x) == :(b^-1*sech(x)+a*b^-2*arctan(sinh(x))+-1*b^-2*(a^2+-1*b^2)^(1/2)*arctan((a^2+-1*b^2)^-1//2*(b+a*tanh(x))*cosh(x)))
@test integrate((a+b*tanh(x))^-1*sech(x), x) == :((a^2+-1*b^2)^-1//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tanh(x))*cosh(x)))
@test integrate((a+b*tanh(x))^-1*cosh(x), x) == :(a*(a^2+-1*b^2)^-1*sinh(x)+-1*b*(a^2+-1*b^2)^-1*cosh(x)+-1*b^2*(a^2+-1*b^2)^-3//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tanh(x))*cosh(x)))
@test integrate((a+b*tanh(x))^-1*cosh(x)^3, x) == :(a*(a^2+-1*b^2)^-1*sinh(x)+a*(-3*b^2+3*a^2)^-1*sinh(x)^3+b^3*(a^2+-1*b^2)^-2*cosh(x)+b^4*(a^2+-1*b^2)^-5//2*arctan((a^2+-1*b^2)^-1//2*(b+a*tanh(x))*cosh(x))+-1*b*(-3*b^2+3*a^2)^-1*cosh(x)^3+-1*a*b^2*(a^2+-1*b^2)^-2*sinh(x))
@test integrate((1+tanh(x))^-1*tanh(x)^5, x) == :(tanh(x)^2+-2*log(cosh(x))+-5//2*tanh(x)+-5//6*tanh(x)^3+5//2*x+(2+2*tanh(x))^-1*tanh(x)^4)
@test integrate((1+tanh(x))^-1*tanh(x)^4, x) == :(-1*tanh(x)^2+2*log(cosh(x))+-3//2*x+3//2*tanh(x)+(2+2*tanh(x))^-1*tanh(x)^3)
@test integrate((1+tanh(x))^-1*tanh(x)^3, x) == :(-1*log(cosh(x))+-3//2*tanh(x)+3//2*x+(2+2*tanh(x))^-1*tanh(x)^2)
@test integrate((1+tanh(x))^-1*tanh(x)^2, x) == :(-1*(2+2*tanh(x))^-1+-1//2*x+log(cosh(x)))
@test integrate((1+tanh(x))^-1*tanh(x), x) == :((2+2*tanh(x))^-1+(1/2)*x)
@test integrate((1+tanh(x))^-1, x) == :((1/2)*x+-1*(2+2*tanh(x))^-1)
@test integrate((1+tanh(x))^-1*coth(x), x) == :((2+2*tanh(x))^-1+-1//2*x+log(sinh(x)))
@test integrate((1+tanh(x))^-1*coth(x)^2, x) == :(-1*log(sinh(x))+-3//2*coth(x)+3//2*x+(2+2*tanh(x))^-1*coth(x))
@test integrate((1+tanh(x))^-1*coth(x)^3, x) == :(-1*coth(x)^2+2*log(sinh(x))+-3//2*x+3//2*coth(x)+(2+2*tanh(x))^-1*coth(x)^2)
@test integrate((1+tanh(x))^-1*coth(x)^4, x) == :(coth(x)^2+-2*log(sinh(x))+-5//2*coth(x)+-5//6*coth(x)^3+5//2*x+(2+2*tanh(x))^-1*coth(x)^3)
@test integrate((1+tanh(x))^3//2*tanh(x), x) == :(-2*(1+tanh(x))^(1/2)+-2//3*(1+tanh(x))^3//2+2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^(1/2)*tanh(x), x) == :(-2*(1+tanh(x))^(1/2)+2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^-1//2*tanh(x), x) == :((1+tanh(x))^-1//2+(1/2)*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^-3//2*tanh(x), x) == :(-1//2*(1+tanh(x))^-1//2+1//3*(1+tanh(x))^-3//2+1//4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^3//2*tanh(x)^2, x) == :(-2*(1+tanh(x))^(1/2)+-2//5*(1+tanh(x))^5//2+2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^(1/2)*tanh(x)^2, x) == :(-2//3*(1+tanh(x))^3//2+2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^-1//2*tanh(x)^2, x) == :(-1*(1+tanh(x))^-1//2+-2*(1+tanh(x))^(1/2)+(1/2)*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((1+tanh(x))^-3//2*tanh(x)^2, x) == :(-1//3*(1+tanh(x))^-3//2+3//2*(1+tanh(x))^-1//2+1//4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+tanh(x))^(1/2)))
@test integrate((a+b*tanh(x))^-1*tanh(x)^5, x) == :(-1//3*b^-1*tanh(x)^3+a*(a^2+-1*b^2)^-1*log(cosh(x))+(1/2)*a*b^-2*tanh(x)^2+-1*b*x*(a^2+-1*b^2)^-1+-1*b^-3*(a^2+b^2)*tanh(x)+a^5*b^-4*(a^2+-1*b^2)^-1*log(a+b*tanh(x)))
@test integrate((a+b*tanh(x))^-1*tanh(x)^4, x) == :(-1//2*b^-1*tanh(x)^2+a*x*(a^2+-1*b^2)^-1+a*b^-2*tanh(x)+-1*b*(a^2+-1*b^2)^-1*log(cosh(x))+-1*a^4*b^-3*(a^2+-1*b^2)^-1*log(a+b*tanh(x)))
@test integrate((a+b*tanh(x))^-1*tanh(x)^3, x) == :(-1*b^-1*tanh(x)+a*(a^2+-1*b^2)^-1*log(cosh(x))+-1*b*x*(a^2+-1*b^2)^-1+a^3*b^-2*(a^2+-1*b^2)^-1*log(a+b*tanh(x)))
@test integrate((a+b*tanh(x))^-1*tanh(x)^2, x) == :(b^-1*log(cosh(x))+-1*a*x*b^-2+x*a^3*b^-2*(a^2+-1*b^2)^-1+-1*a^2*b^-1*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x)))
@test integrate((a+b*tanh(x))^-1*tanh(x), x) == :(a*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x))+-1*b*x*(a^2+-1*b^2)^-1)
@test integrate((a+b*tanh(x))^-1, x) == :(a*x*(a^2+-1*b^2)^-1+-1*b*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x)))
@test integrate((a+b*tanh(x))^-1*coth(x), x) == :(a^-1*log(sinh(x))+-1*b*x*(a^2+-1*b^2)^-1+a^-1*b^2*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x)))
@test integrate((a+b*tanh(x))^-1*coth(x)^2, x) == :(-1*a^-1*coth(x)+a*x*(a^2+-1*b^2)^-1+-1*b*a^-2*log(sinh(x))+-1*a^-2*b^3*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x)))
@test integrate((a+b*tanh(x))^-1*coth(x)^3, x) == :(-1//2*a^-1*coth(x)^2+b*a^-2*coth(x)+a^-3*(a^2+b^2)*log(sinh(x))+-1*b*x*(a^2+-1*b^2)^-1+a^-3*b^4*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x)))
@test integrate((a+b*tanh(x))^-1*coth(x)^4, x) == :(-1//3*a^-1*coth(x)^3+a*x*(a^2+-1*b^2)^-1+(1/2)*b*a^-2*coth(x)^2+-1*a^-3*(a^2+b^2)*coth(x)+-1*b*a^-4*(a^2+b^2)*log(sinh(x))+-1*a^-4*b^5*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x)))
@test integrate(x*(a+b*tanh(x))^-2*sech(x)^2, x) == :(-1*(a^2+-1*b^2)^-1*log(a*cosh(x)+b*sinh(x))+-1*x*b^-1*(a+b*tanh(x))^-1+a*x*b^-1*(a^2+-1*b^2)^-1)
@test integrate(x*(a+b*tanh(c+d*x)^2)^-1*sech(c+d*x)^2, x) == :(-1//4*b^-1//2*d^-2*(-1a)^-1//2*Polylogarithms.polylog(2,-1*(a+-1b+2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+1//4*b^-1//2*d^-2*(-1a)^-1//2*Polylogarithms.polylog(2,-1*(a+-1b+-2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+(1/2)*x*b^-1//2*d^-1*(-1a)^-1//2*log(1+(a+-1b+-2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+-1//2*x*b^-1//2*d^-1*(-1a)^-1//2*log(1+(a+-1b+2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x)))
@test integrate(x^2*(a+b*tanh(c+d*x)^2)^-1*sech(c+d*x)^2, x) == :(-1//4*b^-1//2*d^-3*(-1a)^-1//2*Polylogarithms.polylog(3,-1*(a+-1b+-2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+1//4*b^-1//2*d^-3*(-1a)^-1//2*Polylogarithms.polylog(3,-1*(a+-1b+2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+(1/2)*x*b^-1//2*d^-2*(-1a)^-1//2*Polylogarithms.polylog(2,-1*(a+-1b+-2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+(1/2)*b^-1//2*d^-1*x^2*(-1a)^-1//2*log(1+(a+-1b+-2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+-1//2*x*b^-1//2*d^-2*(-1a)^-1//2*Polylogarithms.polylog(2,-1*(a+-1b+2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x))+-1//2*b^-1//2*d^-1*x^2*(-1a)^-1//2*log(1+(a+-1b+2*b^(1/2)*(-1a)^(1/2))^-1*(a+b)*exp(2c+2*d*x)))
@test integrate(x^-1*tanh(a+b*log(c*x^n)), x) == :(b^-1*n^-1*log(cosh(a+b*log(c*x^n))))
@test integrate(x^-1*tanh(a+b*log(c*x^n))^2, x) == :(-1*b^-1*n^-1*tanh(a+b*log(c*x^n))+log(x))
@test integrate(x^-1*tanh(a+b*log(c*x^n))^3, x) == :(b^-1*n^-1*log(cosh(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*tanh(a+b*log(c*x^n))^2)
@test integrate(x^-1*tanh(a+b*log(c*x^n))^4, x) == :(-1*b^-1*n^-1*tanh(a+b*log(c*x^n))+-1//3*b^-1*n^-1*tanh(a+b*log(c*x^n))^3+log(x))
@test integrate(x^-1*tanh(a+b*log(c*x^n))^5, x) == :(b^-1*n^-1*log(cosh(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*tanh(a+b*log(c*x^n))^2+-1//4*b^-1*n^-1*tanh(a+b*log(c*x^n))^4)
@test integrate(x^-1*tanh(a+b*log(c*x^n))^5//2, x) == :(b^-1*n^-1*arctanh(tanh(a+b*log(c*x^n))^(1/2))+-1*b^-1*n^-1*arctan(tanh(a+b*log(c*x^n))^(1/2))+-2//3*b^-1*n^-1*tanh(a+b*log(c*x^n))^3//2)
@test integrate(x^-1*tanh(a+b*log(c*x^n))^3//2, x) == :(b^-1*n^-1*arctan(tanh(a+b*log(c*x^n))^(1/2))+b^-1*n^-1*arctanh(tanh(a+b*log(c*x^n))^(1/2))+-2*b^-1*n^-1*tanh(a+b*log(c*x^n))^(1/2))
@test integrate(x^-1*tanh(a+b*log(c*x^n))^(1/2), x) == :(b^-1*n^-1*arctanh(tanh(a+b*log(c*x^n))^(1/2))+-1*b^-1*n^-1*arctan(tanh(a+b*log(c*x^n))^(1/2)))
@test integrate(x^-1*tanh(a+b*log(c*x^n))^-1//2, x) == :(b^-1*n^-1*arctan(tanh(a+b*log(c*x^n))^(1/2))+b^-1*n^-1*arctanh(tanh(a+b*log(c*x^n))^(1/2)))
@test integrate(x^-1*tanh(a+b*log(c*x^n))^-3//2, x) == :(b^-1*n^-1*arctanh(tanh(a+b*log(c*x^n))^(1/2))+-1*b^-1*n^-1*arctan(tanh(a+b*log(c*x^n))^(1/2))+-2*b^-1*n^-1*tanh(a+b*log(c*x^n))^-1//2)
@test integrate(x^-1*tanh(a+b*log(c*x^n))^-5//2, x) == :(b^-1*n^-1*arctan(tanh(a+b*log(c*x^n))^(1/2))+b^-1*n^-1*arctanh(tanh(a+b*log(c*x^n))^(1/2))+-2//3*b^-1*n^-1*tanh(a+b*log(c*x^n))^-3//2)
@test integrate((a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*tanh(x)^5, x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2a+tanh(x)^2*(b+2c)))+-1//2*c^-1*(a+b*tanh(x)^2+c*tanh(x)^4)^(1/2)+1//4*c^-3//2*(b+-2c)*arctanh((1/2)*c^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2*c*tanh(x)^2)))
@test integrate((a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*tanh(x)^3, x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2a+tanh(x)^2*(b+2c)))+-1//2*c^-1//2*arctanh((1/2)*c^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2*c*tanh(x)^2)))
@test integrate((a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*tanh(x), x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2a+tanh(x)^2*(b+2c))))
@test integrate((a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*coth(x), x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2a+tanh(x)^2*(b+2c)))+-1//2*a^-1//2*arctanh((1/2)*a^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(2a+b*tanh(x)^2)))
@test integrate((a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*coth(x)^3, x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2a+tanh(x)^2*(b+2c)))+-1//2*a^-1//2*arctanh((1/2)*a^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(2a+b*tanh(x)^2))+-1//2*a^-1*(a+b*tanh(x)^2+c*tanh(x)^4)^(1/2)*coth(x)^2+1//4*b*a^-3//2*arctanh((1/2)*a^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(2a+b*tanh(x)^2)))
@test integrate((a+b*tanh(x)^2+c*tanh(x)^4)^(1/2)*tanh(x), x) == :(-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^(1/2)+(1/2)*(a+b+c)^(1/2)*arctanh((1/2)*(a+b+c)^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2a+tanh(x)^2*(b+2c)))+-1//4*c^-1//2*(b+2c)*arctanh((1/2)*c^-1//2*(a+b*tanh(x)^2+c*tanh(x)^4)^-1//2*(b+2*c*tanh(x)^2)))
@test integrate(tanh(a+b*x)^4*exp(a+b*x), x) == :(b^-1*exp(a+b*x)+-3*b^-1*arctan(exp(a+b*x))+5*b^-1*(1+exp(2a+2*b*x))^-1*exp(a+b*x)+-14//3*b^-1*(1+exp(2a+2*b*x))^-2*exp(a+b*x)+8//3*b^-1*(1+exp(2a+2*b*x))^-3*exp(a+b*x))
@test integrate(tanh(a+b*x)^3*exp(a+b*x), x) == :(b^-1*exp(a+b*x)+-3*b^-1*arctan(exp(a+b*x))+-2*b^-1*(1+exp(2a+2*b*x))^-2*exp(a+b*x)+3*b^-1*(1+exp(2a+2*b*x))^-1*exp(a+b*x))
@test integrate(tanh(a+b*x)^2*exp(a+b*x), x) == :(b^-1*exp(a+b*x)+-2*b^-1*arctan(exp(a+b*x))+2*b^-1*(1+exp(2a+2*b*x))^-1*exp(a+b*x))
@test integrate(exp(a+b*x)*tanh(a+b*x), x) == :(b^-1*exp(a+b*x)+-2*b^-1*arctan(exp(a+b*x)))
@test integrate(coth(a+b*x)*exp(a+b*x), x) == :(b^-1*exp(a+b*x)+-2*b^-1*arctanh(exp(a+b*x)))
@test integrate(coth(a+b*x)^2*exp(a+b*x), x) == :(b^-1*exp(a+b*x)+-2*b^-1*arctanh(exp(a+b*x))+2*b^-1*(1+-1*exp(2a+2*b*x))^-1*exp(a+b*x))
@test integrate(coth(a+b*x)^3*exp(a+b*x), x) == :(b^-1*exp(a+b*x)+-3*b^-1*arctanh(exp(a+b*x))+-2*b^-1*(1+-1*exp(2a+2*b*x))^-2*exp(a+b*x)+3*b^-1*(1+-1*exp(2a+2*b*x))^-1*exp(a+b*x))
@test integrate(coth(a+b*x)^4*exp(a+b*x), x) == :(b^-1*exp(a+b*x)+-3*b^-1*arctanh(exp(a+b*x))+5*b^-1*(1+-1*exp(2a+2*b*x))^-1*exp(a+b*x)+-14//3*b^-1*(1+-1*exp(2a+2*b*x))^-2*exp(a+b*x)+8//3*b^-1*(1+-1*exp(2a+2*b*x))^-3*exp(a+b*x))
@test integrate(tanh(2x)^2*exp(x), x) == :((1+exp(4x))^-1*exp(x)+-1//4*2^(1/2)*arctan(1+2^(1/2)*exp(x))+-1//8*2^(1/2)*log(1+2^(1/2)*exp(x)+exp(2x))+1//4*2^(1/2)*arctan(1+-1*2^(1/2)*exp(x))+1//8*2^(1/2)*log(1+-1*2^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate(exp(x)*tanh(2x), x) == :((1/2)*2^(1/2)*arctan(1+-1*2^(1/2)*exp(x))+-1//2*2^(1/2)*arctan(1+2^(1/2)*exp(x))+-1//4*2^(1/2)*log(1+2^(1/2)*exp(x)+exp(2x))+1//4*2^(1/2)*log(1+-1*2^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate(coth(2x)*exp(x), x) == :(-1*arctan(exp(x))+-1*arctanh(exp(x))+exp(x))
@test integrate(coth(2x)^2*exp(x), x) == :(-1//2*arctan(exp(x))+-1//2*arctanh(exp(x))+(1+-1*exp(4x))^-1*exp(x)+exp(x))
@test integrate(tanh(3x)^2*exp(x), x) == :(-2//9*arctan(exp(x))+-1//9*arctan(3^(1/2)+2*exp(x))+1//9*arctan(3^(1/2)+-2*exp(x))+2*(3+3*exp(6x))^-1*exp(x)+-1//18*3^(1/2)*log(1+3^(1/2)*exp(x)+exp(2x))+1//18*3^(1/2)*log(1+-1*3^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate(exp(x)*tanh(3x), x) == :(-2//3*arctan(exp(x))+-1//3*arctan(3^(1/2)+2*exp(x))+1//3*arctan(3^(1/2)+-2*exp(x))+-1//6*3^(1/2)*log(1+3^(1/2)*exp(x)+exp(2x))+1//6*3^(1/2)*log(1+-1*3^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate(coth(3x)*exp(x), x) == :(-2//3*arctanh(exp(x))+-1//6*log(1+exp(x)+exp(2x))+1//6*log(1+-1*exp(x)+exp(2x))+-1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+2*exp(x)))+1//3*3^(1/2)*arctan(1//3*3^(1/2)*(1+-2*exp(x)))+exp(x))
@test integrate(coth(3x)^2*exp(x), x) == :(-2//9*arctanh(exp(x))+-1//18*log(1+exp(x)+exp(2x))+1//18*log(1+-1*exp(x)+exp(2x))+2*(3+-3*exp(6x))^-1*exp(x)+-1//9*3^(1/2)*arctan(1//3*3^(1/2)*(1+2*exp(x)))+1//9*3^(1/2)*arctan(1//3*3^(1/2)*(1+-2*exp(x)))+exp(x))
@test integrate(tanh(4x)^2*exp(x), x) == :((2+2*exp(8x))^-1*exp(x)+-1//8*(4+-2*2^(1/2))^-1//2*arctan((2+2^(1/2))^-1//2*((2+-1*2^(1/2))^(1/2)+2*exp(x)))+-1//8*(4+2*2^(1/2))^-1//2*arctan((2+-1*2^(1/2))^-1//2*((2+2^(1/2))^(1/2)+2*exp(x)))+-1//32*(2+2^(1/2))^(1/2)*log(1+(2+2^(1/2))^(1/2)*exp(x)+exp(2x))+-1//32*(2+-1*2^(1/2))^(1/2)*log(1+(2+-1*2^(1/2))^(1/2)*exp(x)+exp(2x))+1//8*(4+-2*2^(1/2))^-1//2*arctan((2+2^(1/2))^-1//2*((2+-1*2^(1/2))^(1/2)+-2*exp(x)))+1//8*(4+2*2^(1/2))^-1//2*arctan((2+-1*2^(1/2))^-1//2*((2+2^(1/2))^(1/2)+-2*exp(x)))+1//32*(2+2^(1/2))^(1/2)*log(1+-1*(2+2^(1/2))^(1/2)*exp(x)+exp(2x))+1//32*(2+-1*2^(1/2))^(1/2)*log(1+-1*(2+-1*2^(1/2))^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate(exp(x)*tanh(4x), x) == :((1/2)*(4+-2*2^(1/2))^-1//2*arctan((2+2^(1/2))^-1//2*((2+-1*2^(1/2))^(1/2)+-2*exp(x)))+(1/2)*(4+2*2^(1/2))^-1//2*arctan((2+-1*2^(1/2))^-1//2*((2+2^(1/2))^(1/2)+-2*exp(x)))+-1//2*(4+-2*2^(1/2))^-1//2*arctan((2+2^(1/2))^-1//2*((2+-1*2^(1/2))^(1/2)+2*exp(x)))+-1//2*(4+2*2^(1/2))^-1//2*arctan((2+-1*2^(1/2))^-1//2*((2+2^(1/2))^(1/2)+2*exp(x)))+-1//8*(2+2^(1/2))^(1/2)*log(1+(2+2^(1/2))^(1/2)*exp(x)+exp(2x))+-1//8*(2+-1*2^(1/2))^(1/2)*log(1+(2+-1*2^(1/2))^(1/2)*exp(x)+exp(2x))+1//8*(2+2^(1/2))^(1/2)*log(1+-1*(2+2^(1/2))^(1/2)*exp(x)+exp(2x))+1//8*(2+-1*2^(1/2))^(1/2)*log(1+-1*(2+-1*2^(1/2))^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate(coth(4x)*exp(x), x) == :(-1//2*arctan(exp(x))+-1//2*arctanh(exp(x))+-1//4*2^(1/2)*arctan(1+2^(1/2)*exp(x))+-1//8*2^(1/2)*log(1+2^(1/2)*exp(x)+exp(2x))+1//4*2^(1/2)*arctan(1+-1*2^(1/2)*exp(x))+1//8*2^(1/2)*log(1+-1*2^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate(coth(4x)^2*exp(x), x) == :(-1//8*arctan(exp(x))+-1//8*arctanh(exp(x))+(2+-2*exp(8x))^-1*exp(x)+-1//16*2^(1/2)*arctan(1+2^(1/2)*exp(x))+-1//32*2^(1/2)*log(1+2^(1/2)*exp(x)+exp(2x))+1//16*2^(1/2)*arctan(1+-1*2^(1/2)*exp(x))+1//32*2^(1/2)*log(1+-1*2^(1/2)*exp(x)+exp(2x))+exp(x))
@test integrate((a+-1*tanh(2x))^-1*exp(x), x) == :(-1*(1+-1a)^-1*exp(x)+(1+a)^-1//2*(1+-1a)^-1*(1+-1*a^2)^-1//4*arctan((1+a)^-1//4*(1+-1a)^1//4*exp(x))+(1+a)^-1//2*(1+-1a)^-1*(1+-1*a^2)^-1//4*arctanh((1+a)^-1//4*(1+-1a)^1//4*exp(x)))
@test integrate((a+-1*tanh(2x))^-2*exp(x), x) == :((1+-1a)^-2*exp(x)+(1+a)^-1*(1+-1a)^-2*(1+a+(-1+a)*exp(4x))^-1*exp(x)+-1//2*(1+a)^-3//2*(1+-1a)^-2*(1+-1*a^2)^-1//4*(1+4a)*arctan((1+a)^-1//4*(1+-1a)^1//4*exp(x))+-1//2*(1+a)^-3//2*(1+-1a)^-2*(1+-1*a^2)^-1//4*(1+4a)*arctanh((1+a)^-1//4*(1+-1a)^1//4*exp(x)))
@test integrate((tanh(a*c+b*c*x)^2)^5//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*coth(a*c+b*c*x)*exp(c*(a+b*x))+-15//4*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*arctan(exp(c*(a+b*x)))*coth(a*c+b*c*x)+-4*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-4*coth(a*c+b*c*x)*exp(c*(a+b*x))+-55//6*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-2*coth(a*c+b*c*x)*exp(c*(a+b*x))+25//4*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-1*coth(a*c+b*c*x)*exp(c*(a+b*x))+26//3*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-3*coth(a*c+b*c*x)*exp(c*(a+b*x)))
@test integrate((tanh(a*c+b*c*x)^2)^3//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*coth(a*c+b*c*x)*exp(c*(a+b*x))+-3*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*arctan(exp(c*(a+b*x)))*coth(a*c+b*c*x)+-2*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-2*coth(a*c+b*c*x)*exp(c*(a+b*x))+3*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*(1+exp(2*c*(a+b*x)))^-1*coth(a*c+b*c*x)*exp(c*(a+b*x)))
@test integrate((tanh(a*c+b*c*x)^2)^(1/2)*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*coth(a*c+b*c*x)*exp(c*(a+b*x))+-2*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^(1/2)*arctan(exp(c*(a+b*x)))*coth(a*c+b*c*x))
@test integrate((tanh(a*c+b*c*x)^2)^-1//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-2*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*arctanh(exp(c*(a+b*x)))*tanh(a*c+b*c*x))
@test integrate((tanh(a*c+b*c*x)^2)^-3//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-3*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*arctanh(exp(c*(a+b*x)))*tanh(a*c+b*c*x)+-2*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*(1+-1*exp(2*c*(a+b*x)))^-2*exp(c*(a+b*x))*tanh(a*c+b*c*x)+3*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*(1+-1*exp(2*c*(a+b*x)))^-1*exp(c*(a+b*x))*tanh(a*c+b*c*x))
@test integrate((tanh(a*c+b*c*x)^2)^-5//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-15//4*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*arctanh(exp(c*(a+b*x)))*tanh(a*c+b*c*x)+-4*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*(1+-1*exp(2*c*(a+b*x)))^-4*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-55//6*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*(1+-1*exp(2*c*(a+b*x)))^-2*exp(c*(a+b*x))*tanh(a*c+b*c*x)+25//4*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*(1+-1*exp(2*c*(a+b*x)))^-1*exp(c*(a+b*x))*tanh(a*c+b*c*x)+26//3*b^-1*c^-1*(tanh(a*c+b*c*x)^2)^-1//2*(1+-1*exp(2*c*(a+b*x)))^-3*exp(c*(a+b*x))*tanh(a*c+b*c*x))
@test integrate(sin(tanh(a+b*x))^3, x) == :(-3//8*b^-1*cosintegral(1+-1*tanh(a+b*x))*sin(1)+-3//8*b^-1*cosintegral(1+tanh(a+b*x))*sin(1)+-1//8*b^-1*cos(3)*sinintegral(3+-3*tanh(a+b*x))+-1//8*b^-1*cos(3)*sinintegral(3+3*tanh(a+b*x))+1//8*b^-1*cosintegral(3+-3*tanh(a+b*x))*sin(3)+1//8*b^-1*cosintegral(3+3*tanh(a+b*x))*sin(3)+3//8*b^-1*cos(1)*sinintegral(1+-1*tanh(a+b*x))+3//8*b^-1*cos(1)*sinintegral(1+tanh(a+b*x)))
@test integrate(sin(tanh(a+b*x))^2, x) == :(-1//4*b^-1*log(1+-1*tanh(a+b*x))+1//4*b^-1*log(1+tanh(a+b*x))+-1//4*b^-1*cos(2)*cosintegral(2+2*tanh(a+b*x))+-1//4*b^-1*sin(2)*sinintegral(2+2*tanh(a+b*x))+1//4*b^-1*cos(2)*cosintegral(2+-2*tanh(a+b*x))+1//4*b^-1*sin(2)*sinintegral(2+-2*tanh(a+b*x)))
@test integrate(sin(tanh(a+b*x)), x) == :((1/2)*b^-1*cos(1)*sinintegral(1+-1*tanh(a+b*x))+(1/2)*b^-1*cos(1)*sinintegral(1+tanh(a+b*x))+-1//2*b^-1*cosintegral(1+-1*tanh(a+b*x))*sin(1)+-1//2*b^-1*cosintegral(1+tanh(a+b*x))*sin(1))
@test integrate(cos(tanh(a+b*x))^3, x) == :(-3//8*b^-1*cos(1)*cosintegral(1+-1*tanh(a+b*x))+-3//8*b^-1*sin(1)*sinintegral(1+-1*tanh(a+b*x))+-1//8*b^-1*cos(3)*cosintegral(3+-3*tanh(a+b*x))+-1//8*b^-1*sin(3)*sinintegral(3+-3*tanh(a+b*x))+1//8*b^-1*cos(3)*cosintegral(3+3*tanh(a+b*x))+1//8*b^-1*sin(3)*sinintegral(3+3*tanh(a+b*x))+3//8*b^-1*cos(1)*cosintegral(1+tanh(a+b*x))+3//8*b^-1*sin(1)*sinintegral(1+tanh(a+b*x)))
@test integrate(cos(tanh(a+b*x))^2, x) == :(-1//4*b^-1*log(1+-1*tanh(a+b*x))+1//4*b^-1*log(1+tanh(a+b*x))+-1//4*b^-1*cos(2)*cosintegral(2+-2*tanh(a+b*x))+-1//4*b^-1*sin(2)*sinintegral(2+-2*tanh(a+b*x))+1//4*b^-1*cos(2)*cosintegral(2+2*tanh(a+b*x))+1//4*b^-1*sin(2)*sinintegral(2+2*tanh(a+b*x)))
@test integrate(cos(tanh(a+b*x)), x) == :((1/2)*b^-1*cos(1)*cosintegral(1+tanh(a+b*x))+(1/2)*b^-1*sin(1)*sinintegral(1+tanh(a+b*x))+-1//2*b^-1*cos(1)*cosintegral(1+-1*tanh(a+b*x))+-1//2*b^-1*sin(1)*sinintegral(1+-1*tanh(a+b*x)))

max_stars_repo_path: test/6 Hyperbolic functions/6.3 Hyperbolic tangent/6.3.2 Hyperbolic tangent functions.jl


text:
struct IndexedCounts{M<:AbstractMer,C<:CountMode}
    mer_index::Vector{M}            # Ordered list of kmer that contain counts
    count_names::Vector{Symbol}     # Names of the count vectors
    counts::Vector{Vector{UInt16}}  # Count vectors, each contains an entry per mer in the index
    name::Symbol                    # The name of this MerCounts
end

# Construction of empty MerCounts with a name
@inline function IndexedCounts{M,C}(name::Symbol) where {M<:AbstractMer,C<:CountMode}
    return IndexedCounts{M,C}(Vector{M}(), Vector{Symbol}(), Vector{Vector{UInt16}}(), name)
end

@inline function IndexedCounts(::Type{M}, name::Symbol, mode::C = CANONICAL) where {M<:AbstractMer,C<:CountMode}
    return IndexedCounts{M,C}(name)
end

# Construction from a sequence distance graph
function IndexedCounts(name::Symbol, sequences, mode::C = CANONICAL) where {M<:AbstractMer,C<:CountMode}
    @info string("Creating an empty kmer count datastore called ", name)
    mc = IndexedCounts{M,C}(name, mode)
    index!(mc, sequences)
    return mc
end

function _count_and_collapse!(mers::Vector{M}, counts::Vector{UInt16}) where {M<:AbstractMer}
    wi = 1
    ri = 1
    stop = lastindex(mers) + 1
    @inbounds while ri < stop
        mers[wi] = mers[ri]
        ci = one(UInt16)
        while (ri += 1) < stop && mers[ri] == mers[wi]
            ci = ci + one(UInt16)
        end
        push!(counts, ci)
        wi = wi + 1
    end
    resize!(mers, ri)
    return nothing
end

function index!(mc::IndexedCounts{M,C}, ref_sequences::Vector{<:LongSequence}) where {M<:AbstractMer,C<:CountMode}
    @info "Indexing kmer counts of reference sequences"
    @info "Populating index with reference sequence mers"
    # Add all Kmers from reference sequences.
    # t = expected number of kmers.
    t = sum(length(x) for x in ref_sequences) - BioSequences.ksize(M) + 1
    mcindex = mc.mer_index
    resize!(mcindex, t)
    vec_i = 1
    f = C()
    for seq in ref_sequences
        for mer in each(M, seq)
            mcindex[vec_i] = f(mer)
            vec_i = vec_i + 1
        end
    end
    resize!(mcindex, vec_i)
    @info "Sorting the index"
    # Sort the index
    sort!(mcindex)
    @info "Counting and collapsing mers in index"
    # Create a first count;
    # Collapse the kmer index, saving the coverage to the first count.
    first_count = Vector{UInt16}()
    sizehint!(first_count, length(mcindex))
    push!(mc.counts, first_count)
    push!(mc.count_names, :sdg)
    _count_and_collapse!(mcindex, first_count)
    return mc
end

function Base.summary(io::IO, mc::IndexedCounts{M}) where {M}
    print(io, BioSequences.ksize(M), "-mer Indexed Counts Datastore '",  mc.name, "': ", length(mc.counts), " stored ", BioSequences.ksize(M), "-mer counts")
end

function Base.show(io::IO, mc::IndexedCounts{M}) where {M}
    summary(io, mc)
end

ext: jl


text:
## Dimensions and descriptions

desc = doc"
# Describing `RLEVector` Objects
For an RLEVector `x = RLEVector([4,5,6],[3,6,9])`
 * `length(x)` The full length of the vector, uncompressed
 * `size(x)` Same as `length`, as for any other vector
 * `size(x,dim)` Returns `(length(x),1) for dim == 1`
 * `rstart(x)` The index of the beginning of each run
 * `rwidth(x)` The width of each run
 * `rstop(x)` The index of the end of each run
 * `rlast(x)` The data value for each run
 * `isempty(x)` Returns boolean, as for any other vector
 * `nrun(x)` Returns the number of runs represented in the array
 * `endtype(x)` Returns the element type of the run ends

## See also
length, size, rstart, rwidth, rstop, rvalue, isempty, nrun, endtype
"

@doc desc ->
function ndims(x::RLEVector)
    return(1)
end

@doc desc ->
function nrun(x::RLEVector)
  length(x.runends)
end

@doc desc ->
function length{T1,T2<:Integer}(x::RLEVector{T1,T2})
  re = x.runends
  ind = endof(re)
  if (ind > 0)
    @inbounds len = re[ind]
  else
    len = zero(T2)
  end
  return(len)
end

@doc desc ->
function size(x::RLEVector)
  (length(x),)
end

@doc desc ->
function size(x::RLEVector, dim::Integer)
  len = length(x)
  if dim == 1
    return(len)
  else
    return( (len,1) )
  end
end

@doc desc ->
function isempty(x::RLEVector)
  isempty(x.runends)
end

### Getters
@doc desc ->
function rfirst(x::RLEVector)
  re = x.runends
  rval = similar(re)
  prev = zero(eltype(re))
    @inbounds for i in eachindex(re)
        rval[i] = prev + 1
        prev = re[i]
  end
  return(rval)
end

@doc desc ->
function rfirst(x::RLEVector, run::Integer)
    num_one = one(eltype(x.runends))
    run == 1 ? num_one : x.runends[run-1] + num_one
end

@doc desc ->
function rwidth(x::RLEVector)
  re = x.runends
  rval = similar(re)
  prev = zero(eltype(re))
  @inbounds for i in eachindex(re)
      rei = re[i]
      rval[i] = rei - prev
      prev = rei
  end
  return(rval)
end

function rwidth(x::RLEVector, run::Integer)
  run == 1 ? x.runends[1] : x.runends[run] - x.runends[run-1]
end

@doc desc ->
function rlast(x::RLEVector)
  x.runends
end

@doc desc ->
function rvalue(x::RLEVector)
  x.runvalues
end

@doc desc ->
endtype(x::RLEVector) = eltype(rlast(x))


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

# test continuum nonlinear elasticity with surface load

X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [1.0, 1.0, 0.0],
    4 => [0.0, 1.0, 0.0],
    5 => [0.0, 0.0, 1.0],
    6 => [1.0, 0.0, 1.0],
    7 => [1.0, 1.0, 1.0],
    8 => [0.0, 1.0, 1.0])

element1 = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
element2 = Element(Quad4, (5, 6, 7, 8))
update!((element1, element2), "geometry", X)
update!(element1, "youngs modulus", 900.0)
update!(element1, "poissons ratio", 0.25)
update!(element2, "displacement traction force", [0.0, 0.0, -100.0])

problem = Problem(Elasticity, "solve continuum block", 3)
problem.properties.finite_strain = true
add_elements!(problem, element1, element2)

symxy = Element(Quad4, (1, 2, 3, 4))
symxz = Element(Quad4, (1, 2, 6, 5))
symyz = Element(Quad4, (1, 4, 8, 5))
update!((symxy, symxz, symyz), "geometry", X)
update!(symxy, "displacement 3", 0.0)
update!(symxz, "displacement 2", 0.0)
update!(symyz, "displacement 1", 0.0)

bc = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
add_elements!(bc, symxy, symxz, symyz)

analysis = Analysis(Nonlinear)
add_problems!(analysis, problem, bc)
run!(analysis)

u = element1("displacement", [1.0, 1.0, 1.0], 0.0)
u_expected = [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]
@debug("displacement at tip", u, u_expected)

# verified using Code Aster:
# 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm
@test isapprox(u, u_expected; rtol=1.0e-4)


text:
using OptimalTransport

using Distances
using ForwardDiff
using LogExpFunctions
using PythonOT: PythonOT

using LinearAlgebra
using Random
using Test

const POT = PythonOT

Random.seed!(100)

@testset "sinkhorn_epsscaling.jl" begin
    # size of source and target
    M = 250
    N = 200

    # create two random histograms
    μ = normalize!(rand(M), 1)
    ν = normalize!(rand(N), 1)

    # create random cost matrix
    C = pairwise(SqEuclidean(), rand(1, M), rand(1, N); dims=2)

    # regularization parameter
    ε = 0.01

    @testset "example" begin
        # compute optimal transport plan and cost with POT
        γ_pot = POT.sinkhorn(μ, ν, C, ε; method="sinkhorn_stabilized", numItermax=5_000)
        c_pot = POT.sinkhorn2(μ, ν, C, ε; method="sinkhorn_stabilized", numItermax=5_000)[1]

        for alg in (SinkhornGibbs(), SinkhornStabilized())
            # compute optimal transport plan and cost
            γ = sinkhorn(μ, ν, C, ε, SinkhornEpsilonScaling(alg); maxiter=5_000)
            c = sinkhorn2(μ, ν, C, ε, SinkhornEpsilonScaling(alg); maxiter=5_000)

            # check that plan and cost are consistent
            @test c ≈ dot(γ, C)

            # compare with Sinkhorn algorithm without ε-scaling
            γ_wo_epsscaling = sinkhorn(μ, ν, C, ε, alg; maxiter=5_000)
            c_wo_epsscaling = sinkhorn2(μ, ν, C, ε, alg; maxiter=5_000)
            @test γ ≈ γ_wo_epsscaling rtol = 1e-6
            @test c ≈ c_wo_epsscaling

            # compare with POT
            @test γ ≈ γ_pot rtol = 1e-6
            @test c ≈ c_pot rtol = 1e-6
        end
    end

    @testset "AD" begin
        # compute gradients with respect to source and target marginals separately and
        # together
        for alg in (SinkhornGibbs(), SinkhornStabilized())
            epsscaling_alg = SinkhornEpsilonScaling(alg)
            ForwardDiff.gradient(zeros(N)) do xs
                return sinkhorn2(μ, softmax(xs), C, ε, epsscaling_alg)
            end
            ForwardDiff.gradient(zeros(M)) do xs
                return sinkhorn2(softmax(xs), ν, C, ε, epsscaling_alg)
            end
            ForwardDiff.gradient(zeros(M + N)) do xs
                return sinkhorn2(
                    softmax(xs[1:M]), softmax(xs[(M + 1):end]), C, ε, epsscaling_alg
                )
            end
        end
    end

    @testset "deprecations" begin
        alg = SinkhornEpsilonScaling(
            SinkhornStabilized(; absorb_tol=10); factor=0.7, steps=2
        )
        γ = sinkhorn(μ, ν, C, ε, alg)
        @test (@test_deprecated sinkhorn_stabilized_epsscaling(
            μ, ν, C, ε; absorb_tol=10, lambda=0.7, k=2
        )) == γ
        @test (@test_deprecated sinkhorn_stabilized_epsscaling(
            μ, ν, C, ε; absorb_tol=10, scaling_factor=0.7, scaling_steps=2
        )) == γ
    end
end


text:
# X.F. Yang's Double Cosine kernel (2013)
function doublecosine(r, h)
    # Kernel constant
    k = 2
    if abs(r) < 2h
        # Normalizing factor
        σ = 1/(6k * h)

        # Checking sign
        sign = r == 0 ? 1 : r/abs(r)

        # Kernel computation
        W = (4cos(π/(2h) * r) + cos(π/h * r) + 3) * σ

        # Derivatives computation
        ∂x_W = (4π/(2h) * sin(π/(2h) * r) + π/h * sin(π/h * r)) * -sign * σ

        return (W,∂x_W)
    else
        return (0,0)
    end
    res
end

# J.J. Monaghan's Cubic Spline kernel (1992)
function cubicspline(r, h)
    if abs(r) < 2h
        # Reduced radius
        q = r/h

        # Normalizing factor
        σ = 2/(3h)

        # Checking sign
        sign = r == 0 ? 1 : r/abs(r)

        if abs(r) < h
            # Kernel computation
            W = (1 - 1.5q^2 + 0.75q^3) * σ

            # Derivatives computation
            ∂x_W = (2q - 1.5q^2) * sign/h^2

        else #if h <abs(r) < 2h
            # Kernel computation
            W = (2-q)^3/4 * σ

            # Derivatives computation
            ∂x_W = (q-2)^2/2 * sign/h^2
        end
        return (W,∂x_W)
    else
        return (0,0)
    end
    res
end


text:
module Display

# The methodology of these routines is based on an initial implementation by
# @roelmatthysen for the display of complicated structures in the package
# BasisFunctions.jl

export composite_show,
    composite_show_compact,
    composite_displaystencil

using ..CompositeTypes

typename(object) = string(typeof(object))

"What is the standard display symbol of the object?"
displaysymbol(object) = typename(object)[1]
displaysymbol(object::Function) = 'F'

# some common types
displaysymbol(A::AbstractArray) = 'A'
displaysymbol(A::AbstractVector) = 'v'
displaysymbol(x::Number) = 'c'

object_parentheses(a::Complex) = true

"Representation of a symbol."
struct Symbol{S}
    sym ::  S
end

displaysymbol(s::Symbol) = s.sym

Comma() = Symbol(',')
Times() = Symbol('×')
Plus() = Symbol('+')

"""
For composite objects, what is their combination symbol?

By default objects are comma separated, of the form `F(a, b, c)` where `F` is
the constructor function. But they could also have the form `a + b + c`. In the
former case the combination symbol is a comma, in the latter case it is `+`.
"""
combinationsymbol(object) = Comma()

is_comma_separated(object) = combinationsymbol(object) == Comma()

# For the combination of composite objects, what is the infix string?
# It is ", " for comma-separated objects, but " + " with an extra space for other symbols.
infix_string(object) = _infix_string(displaysymbol(combinationsymbol(object)))
_infix_string(s) = s == ',' ? ", " : " " * s * " "

maximum_depth(object) = 3
maximum_replength(object) = 30
maximum_textlength(object) = 60
maximum_children(object) = 5

istext(el) = false
istext(el::Char) = true
istext(el::String) = true

objectref(a) = istext(a) ? TextObject(a) : a

function concatenate_components(object; reversecomponents = false)
    list = components(object)
    if reversecomponents
        list = reverse(list)
    end
    concatenate(list, infix_string(object), maximum_children(object))
end

function concatenate(list, sep, maxchildren = 5)
    A = Any[]
    push!(A, objectref(list[1]))
    if length(list) > 1
        if length(list) <= maxchildren
            for i in 2:length(list)
                push!(A, sep)
                push!(A, objectref(list[i]))
            end
        else
            for i in 2:3
                push!(A, sep)
                push!(A, objectref(list[i]))
            end
            push!(A, sep)
            push!(A, "...")
            push!(A, sep)
            push!(A, objectref(list[end]))
        end
    end
    A
end

"Does the object have a stencil?"
hasstencil(object) = length(displaystencil(object)) > 0

"""
Return the stencil of the object as an array.

The stencil of an object determines how it is displayed. The array may
contain characters or strings, and references to other objects. The concatenation
of all elements, with the objects replaced by their string representation, forms
the representation of the given object.

Example: say an object `obj` consists of two parts, `obj[1]` and `obj[2]`,
combined by a function commonly noted as `I`. Then the stencil may be:
`displaystencil(obj) = ["I(", obj[1], ", ", obj[2], ')']`.

In the string representation of the object, both `obj[1]` and `obj[2]` are
replaced by a string, or by a symbol that has its own string representation
separately.
"""
displaystencil(object) = []

constructorname(object) = typename(object)

"""
Default display stencil for composite objects.

The default leads to the representation `typename(component1, components2, ...)`,
unless the `combinationsymbol` of the object is defined differently. In the
latter case, if the combination symbol is '+' the stencil leads to
`component1 + component2 + ...`.
"""
function composite_displaystencil(object; kwargs...)
    if is_comma_separated(object)
        A1 = [constructorname(object), '(']
        A2 = concatenate_components(object; kwargs...)
        A3 = [')']
        vcat(A1, A2, A3)
    else
        concatenate_components(object; kwargs...)
    end
end

"Objects in a display stencil can be wrapped in certain circumstances."
abstract type WrappedObject end

hasstencil(object::WrappedObject) = hasstencil(object.object)
displaystencil(object::WrappedObject) = displaystencil(object.object)
stencil_parentheses(object::WrappedObject) = stencil_parentheses(object.object)
compact_repr(object::WrappedObject) = compact_repr(object.object)

"""
Display stencils treat text differently from other objects. If one of the
components of an object is text, it can be wrapped in a `TextObject` in order
to treat it as an object, rather than as text, in a stencil.
"""
struct TextObject <: WrappedObject
    object
end

"""
Wrap an object along with a symbol for it.

This can be useful in displaystencils of the form `F(a,b)` when `F` itself has
a stencil. It is usually best to represent `F` by a symbol, and to list the
string representation of `F` separately. In such a case the stencil of the object
`O` represented by `F(a,b)` would be `[SymbolObject(F), '(', a, b, ')']`.

In other words, `SymbolObject` forces an object to be represented by a symbol in
a stencil.
"""
struct SymbolObject{S} <: WrappedObject
    object
    sym     ::  S
end
SymbolObject(object) = SymbolObject(object, displaysymbol(object))
displaysymbol(object::SymbolObject) = object.sym

"Return all objects appearing in the display stencil."
stencil_objects(object) = objects(displaystencil(object))
objects(A::Vector) = unique([el for el in A if !istext(el)])

"""
The conventional compact representation of an object.

This representation is based on `show(io, d)`. If an object has overriden
that method, it can also override this one to avoid a loop.
"""
compact_repr(object) = repr(object)

"""
Recursively compute all stencils of the objects linked to by the given object,
up to a maximum depth of recursion. There can be self-referential objects
in the stencils.

The result is returned as a dictionary that maps objects to stencils.
"""
function recursive_stencils(object, depth = 1, maxdepth = maximum_depth(object),
                                stencils = Dict{Any,Any}())
    if depth > maxdepth
        if object ∉ keys(stencils)
            stencils[object] = [compact_repr(object)]
        end
        return stencils
    end
    if object ∉ keys(stencils)
        if hasstencil(object)
            A = displaystencil(object)
            stencils[object] = A
            for child in objects(A)
                recursive_stencils(child, depth+1, maxdepth, stencils)
            end
        else
            stencils[object] = [compact_repr(object)]
        end
    end
    stencils
end

"If an object in the stencil is complicated, should we put parentheses around it?"
function stencil_parentheses(object)
    if hasstencil(object)
        iscomposite(object) && !is_comma_separated(object)
    else
        true
    end
end

"If this object appears in a more complicated expression, does it require parentheses?"
function object_parentheses(object)
    if iscomposite(object)
        !is_comma_separated(object)
    else
        false
    end
end

"Compute representations of all stencils linked to by the object."
function recursive_stencil_reps(object, stencils, maxreplength, reps = Dict{Any,Any}(), subs = Dict{Any,Any}())
    if object ∈ keys(reps)
        return reps, subs
    end
    if object ∉ keys(stencils)
        # the object may have been skipped due to recursion depth
        reps[object] = compact_repr(object)
        return reps, subs
    end
    A = stencils[object]
    rep = ""
    for element in A
        if istext(element)
            rep *= element
        else
            recursive_stencil_reps(element, stencils, maxreplength, reps, subs)
            rep *= reps[element]
        end
    end
    if (length(rep) > maxreplength) || object isa SymbolObject
        sym = displaysymbol(object)
        reps[object] = string(sym)
        subs[object] = sym
    else
        reps[object] = rep
    end
    reps, subs
end

# Different operators with the same symbol get added subscripts
function subscript(i::Integer)
    @assert i >= 0
    if i < 10
        join('₀'+d for d in reverse(digits(i)))
    else
        subscript(div(i,10)) * subscript(mod(i,10))
    end
end

"Assign symbols to all objects in the list, adding subscripts if necessary."
function allocate_symbols(subs)
    invsubs = Dict{Any,Any}()
    for k in keys(subs)
        if subs[k] in keys(invsubs)
            push!(invsubs[subs[k]],k)
        else
            invsubs[subs[k]] = Any[k]
        end
    end
    syms = Dict{Any,Any}()
    for (sym,objects) in invsubs
        if length(objects) == 1
            syms[objects[1]] = sym
        else
            for i in 1:length(objects)
                syms[objects[i]] = sym * subscript(i)
            end
        end
    end
    syms
end

function substitute_symbols(object, stencils, syms, subs, reps = Dict{Any,Any}())
    if object in keys(reps)
        return reps
    end
    if object ∉ keys(stencils)
        reps[object] = compact_repr(object)
        return reps
    end
    A = stencils[object]
    rep = ""
    for element in A
        if istext(element)
            rep *= element
        elseif element in keys(subs)
            substitute_symbols(element, stencils, syms, subs, reps)
            rep *= syms[element]
        else
            substitute_symbols(element, stencils, syms, subs, reps)
            if object_parentheses(element) && stencil_parentheses(object)
                rep = rep * '(' * reps[element] * ')'
            else
                rep *= reps[element]
            end
        end
    end
    reps[object] = rep
    reps
end

"""
Display multi-line structured information about a composite object,
using the stencil of the object and, recursively, the stencils of any objects
therein.
"""
function composite_show(io::IO, ::MIME"text/plain", object)
    # get(io, :compact) -> this is true when a compact representation is requested,
    #  typically within an array or some other container object. This is not true
    #  for vectors (see https://github.com/JuliaLang/julia/pull/22981)
    # haskey(io, :typeinfo) -> this is true when type information has already
    #  been shown for an array. This is also true for vectors.
    if get(io, :compact, false) || haskey(io, :typeinfo)
        show(io, object)
    else
        composite_show(io, object)
    end
end

function composite_show(io::IO, object)
    stencils = recursive_stencils(object)
    reps, subs = recursive_stencil_reps(object, stencils, maximum_replength(object))
    if object in keys(subs)
        pop!(subs, object)
    end
    syms = allocate_symbols(subs)
    finalreps = substitute_symbols(object, stencils, syms, subs)

    print(io, finalreps[object])
    n = length(keys(syms))
    if n > 0
        print(io, '\n')
        print(io, '\n')
        i = 0
        for (object,sym) in syms
            print(io, sym, " = ", finalreps[object])
            i += 1
            if i < n
                print(io, '\n')
            end
        end
    end
    nothing
end

"""
Display structured information about a composite object using its stencil,
but not recursively.
"""
function composite_show_compact(io::IO, object)
    @assert hasstencil(object)
    A = displaystencil(object)
    for element in A
        if istext(element)
            print(io, element)
        else
            if object_parentheses(element) && stencil_parentheses(object)
                print(io, '(', compact_repr(element), ')')
            else
                print(io, compact_repr(element))
            end
        end
    end
end

end # module


text:
using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear, prob_ode_linear, prob_ode_2Dlinear

probArr = [prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear]
testTol = 0.25
dts = 1 .//(2 .^(10:-1:5))

@testset "Nordsieck Convergence Tests" begin
  for i in eachindex(probArr)
    sim = test_convergence(dts,probArr[i],AN5())
    @test sim.𝒪est[:final] ≈ 5 atol=testTol
    @test sim.𝒪est[:l2] ≈ 5 atol=testTol
    @test sim.𝒪est[:l∞] ≈ 5 atol=testTol
  end
end

probArr = [prob_ode_linear,
           prob_ode_2Dlinear]
@testset "Nordsieck Adaptivity Tests: AN5" begin
  for i in eachindex(probArr)
    prob = probArr[i]
    sol = solve(prob, AN5(), reltol=1e-6)
    @test length(sol.t) < 11
    exact = prob.f(Val{:analytic}, prob.u0, prob.p, prob.tspan[end])
    @test exact ≈ sol[end] atol=1e-5
  end
end

@testset "Nordsieck Adaptivity Tests: JVODE" begin
  for i in eachindex(probArr), sol = [JVODE_Adams(), JVODE_BDF()]
    prob = probArr[i]
    sol = solve(prob, sol, reltol=1e-4, abstol=1e-7)
    @test length(sol.t) < 22
    exact = prob.f(Val{:analytic}, prob.u0, prob.p, prob.tspan[end])
    @test norm(exact - sol[end], Inf) < 3e-3
  end
end


text:
function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.EqualTo)
    set_lower_bound(info, set.value)
    set_upper_bound(info, set.value)
end

function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.GreaterThan)
    set_lower_bound(info, set.lower)
end

function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.LessThan)
    set_upper_bound(info, set.upper)
end

function set_bounds(info::Union{VariableInfo, ConstraintInfo}, set::MOI.Interval)
    set_lower_bound(info, set.lower)
    set_upper_bound(info, set.upper)
end

MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{<:Bounds}) = true

function MOI.add_constraint(model::Optimizer, f::SV, set::S) where {S <: Bounds}
    variable_info = find_variable_info(model, f.variable)
    set_bounds(variable_info, set)
    return CI{SV, S}(f.variable.value)
end

MOI.supports_constraint(::Optimizer, ::Type{<:Union{SAF, SQF}}, ::Type{<:Bounds}) = true

function MOI.add_constraint(model::Optimizer, f::F, set::S) where {F <: Union{SAF, SQF}, S <: Bounds}
    ci = ConstraintInfo()
    ci.expression = to_expr(f)
    set_bounds(ci, set)
    push!(model.inner.constraint_info, ci)
    return CI{F, S}(length(model.inner.constraint_info))
end

MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{CI}) = true

function MOI.set(model::Optimizer, attr::MOI.ConstraintName, ci::CI{SV}, value)
    error("No support for naming constraints imposed on variables.")
end
function MOI.set(model::Optimizer, attr::MOI.ConstraintName, ci::CI, value)
    check_constraint_indices(model, ci)
    model.inner.constraint_info[ci.value].name = value
end
function MOI.get(model::Optimizer, ::MOI.ConstraintName, ci::CI)
     return model.inner.constraint_info[ci.value].name
end

function MOI.get(model::Optimizer, ::Type{MathOptInterface.ConstraintIndex}, name::String)
    for (i,c) in enumerate(model.inner.constraint_info)
        if name == c.name
            return CI(i)
        end
    end
    error("Unrecognized constraint name $name.")
end

MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{MOI.ZeroOne}) = true
MOI.supports_constraint(::Optimizer, ::Type{SV}, ::Type{MOI.Integer}) = true

function MOI.add_constraint(model::Optimizer, f::SV, set::S) where {S <: Union{MOI.ZeroOne, MOI.Integer}}
    variable_info = find_variable_info(model, f.variable)
    if set isa MOI.ZeroOne
        variable_info.category = :Bin
    elseif set isa MOI.Integer
        variable_info.category = :Int
    else
        error("Unsupported variable type $set.")
    end
    return CI{SV, S}(f.variable.value)
end

# MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

# function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_block_data::MOI.NLPBlockData)
#     model.nlp_block_data = nlp_block_data
#     return
# end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct NetworkInterfaceReferenceProperties <: SwaggerModel
    primary::Any # spec type: Union{ Nothing, Bool } # spec name: primary

    function NetworkInterfaceReferenceProperties(;primary=nothing)
        o = new()
        validate_property(NetworkInterfaceReferenceProperties, Symbol("primary"), primary)
        setfield!(o, Symbol("primary"), primary)
        o
    end
end # type NetworkInterfaceReferenceProperties

const _property_map_NetworkInterfaceReferenceProperties = Dict{Symbol,Symbol}(Symbol("primary")=>Symbol("primary"))
const _property_types_NetworkInterfaceReferenceProperties = Dict{Symbol,String}(Symbol("primary")=>"Bool")
Base.propertynames(::Type{ NetworkInterfaceReferenceProperties }) = collect(keys(_property_map_NetworkInterfaceReferenceProperties))
Swagger.property_type(::Type{ NetworkInterfaceReferenceProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkInterfaceReferenceProperties[name]))}
Swagger.field_name(::Type{ NetworkInterfaceReferenceProperties }, property_name::Symbol) =  _property_map_NetworkInterfaceReferenceProperties[property_name]

function check_required(o::NetworkInterfaceReferenceProperties)
    true
end

function validate_property(::Type{ NetworkInterfaceReferenceProperties }, name::Symbol, val)
end


text:
using Test
using FlightMechanics

# TEST rigid_body_velocity
# Test without angular velocity
vel_P = [10.0, 0.0, 0.0]
ω = [0.0, 0.0, 0.0]
r_PQ = [5.0, 5.0, 5.0]

exp_vel_Q = [10.0, 0.0, 0.0]
vel_Q = rigid_body_velocity(vel_P, ω, r_PQ)
@test isapprox(vel_Q, exp_vel_Q)

# Only angular velocity in Z axis
vel_P = [0.0, 0.0, 0.0]
ω = [0.0, 0.0, 1.0]
r_PQ = [5.0, 5.0, 5.0]

exp_vel_Q = [-sqrt(50) * sin(pi/4), sqrt(50) * cos(pi/4), 0.0]
vel_Q = rigid_body_velocity(vel_P, ω, r_PQ)
@test isapprox(vel_Q, exp_vel_Q)

# X linear velocity and Z angular velocity
vel_P = [10.0, 0.0, 0.0]
ω = [0.0, 0.0, 1.0]
r_PQ = [5.0, 5.0, 5.0]

exp_vel_Q = [-sqrt(50) * sin(pi/4) + 10.0, sqrt(50) * cos(pi/4), 0.0]
vel_Q = rigid_body_velocity(vel_P, ω, r_PQ)
@test isapprox(vel_Q, exp_vel_Q)

# TEST rigid_body_acceleration
# Only X linear acceleration
acc_P = [10.0, 0.0, 0.0]
ω = [0.0, 0.0, 0.0]
ω_dot = [0.0, 0.0, 0.0]
r_PQ = [5.0, 5.0, 5.0]

exp_accel_Q = [10.0, 0.0, 0.0]
accel_Q = rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)
@test isapprox(accel_Q, exp_accel_Q)

# Only Z angular acceleration
acc_P = [0.0, 0.0, 0.0]
ω = [0.0, 0.0, 0.0]
ω_dot = [0.0, 0.0, 1.0]
r_PQ = [5.0, 5.0, 5.0]

exp_accel_Q = [-sqrt(50) * sin(pi/4), sqrt(50) * cos(pi/4), 0.0]
accel_Q = rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)
@test isapprox(accel_Q, exp_accel_Q)

# Only Z angular velocity
acc_P = [0.0, 0.0, 0.0]
ω = [0.0, 0.0, 2.0]
ω_dot = [0.0, 0.0, 0.0]
r_PQ = [5.0, 5.0, 5.0]

exp_accel_Q = [-4.0 * sqrt(50) * sin(pi/4), -4.0 * sqrt(50) * cos(pi/4), 0.0]
accel_Q = rigid_body_acceleration(acc_P, ω, ω_dot, r_PQ)
@test isapprox(accel_Q, exp_accel_Q)

# Steiner theorem
# No point translation
point0 = [0, 0, 0]
point1 = point0
inertia = [1 0 0;
           0 2 0;
           0 0 3]
mass = 10
@test isapprox(steiner_inertia(point0, inertia, mass, point1), inertia)

# Point with no inertia
point0 = [0, 0, 0]
point1 = [10, 0, 0]
inertia = zeros(3, 3)
mass = 10

exp_inertia = [0       0     0;
               0    1000     0;
               0       0  1000]
@test isapprox(steiner_inertia(point0, inertia, mass, point1), exp_inertia)

# Point with no inertia
point0 = [0, 0, 0]
point1 = [10, 0, 0]
inertia = [1 0 0;
           0 2 0;
           0 0 3]
mass = 10

exp_inertia = [0       0     0;
               0    1000     0;
               0       0  1000] + inertia
@test isapprox(steiner_inertia(point0, inertia, mass, point1), exp_inertia)


text:
"""
    mutable struct Structure_Basic  <: Str

This contains information for model constructing of a series of basic structures.

# Supported list:
 - structure_sc() (Simple Cubic)
 - structure_bcc() (Body Centered Cubic)
 - structure_fcc() (Face Centered Cubic)
 - structure_dc() (Diamond Cubic)
"""
mutable struct Structure_Basic <: Str
    atom_vec::Matrix{Float64}
    cell_vec::Matrix{Float64}
    atom_type
    atom_charge
    num_atoms
    num_atom_types
    res_name
end

function transform(atom_vec, atom_type, atom_charge)
    num_atoms = size(atom_vec)[1]
    base = convert.(Int64, ones(num_atoms))

    if length(atom_type) == 1
        atom_type = base .* atom_type
        atom_charge = base .* atom_charge
        num_atom_types = 1
    else
        num_atom_types = length(unique(atom_type))
        if length(atom_charge) == 1
            atom_charge = base .* atom_charge
        end
    end

    atom_type, atom_charge, num_atoms, num_atom_types
end

"""
    structure_sc(;atom_type=1, atom_charge=0, lattice_const=1)

Do this will return a Structure_Basic type which contains all information needed to build a SC structure.

# Example
```julia-replr
str = structure_sc(atom_type=1, atom_charge=0, lattice_const=3)
```
"""
function structure_sc(;atom_type=1, atom_charge=0, lattice_const=1)
    atom_vec = [
        0   0   0
    ]
    cell_vec = diag([
        1   1   1
    ] .* lattice_const)
    atom_type, atom_charge, num_atoms, num_atom_types = transform(atom_vec, atom_type, atom_charge)
    res_name = "SC"
    Structure_Basic(atom_vec, cell_vec, atom_type, atom_charge, num_atoms, num_atom_types, res_name)
end

"""
    structure_bcc(;atom_type=1, atom_charge=0, lattice_const=1)

Do this will return a Structure_Basic type which contains all information needed to build a BCC structure.

# Arguments
- `atom_type`: type of each atoms, could be a interger or a vector represent type of all atoms or each atoms repectively.
- `atom_charge`: charge of each atoms, could be a float or a vector, same as `atom_type` Data_Unit: e.
- `lattice_const`: lattice constance of structure

# Example
```julia-repl
str = structure_bcc(atom_type=1, atom_charge=-1, lattice_const=3) # All atom in cell will have same type and charge
```
Or
```julia-replr
str = structure_bcc(atom_type=[1 2], atom_charge=[0 2], lattice_const=3)
```
"""
function structure_bcc(;atom_type=1, atom_charge=0, lattice_const=1)
    atom_vec = [
        0   0   0
        0.5 0.5 0.5
    ]
    cell_vec = diag([
        1   1   1
    ] .* lattice_const)
    atom_type, atom_charge, num_atoms, num_atom_types = transform(atom_vec, atom_type, atom_charge)
    res_name = "BCC"
    Structure_Basic(atom_vec, cell_vec, atom_type, atom_charge, num_atoms, num_atom_types, res_name)
end

"""
    structure_fcc(;atom_type=1, atom_charge=0, lattice_const=1)

Do this will return a Structure_Basic type which contains all information needed to build a FCC structure.

# Arguments
- `atom_type`: type of each atoms, could be a interger or a vector represent type of all atoms or each atoms repectively.
- `atom_charge`: charge of each atoms, could be a float or a vector, same as `atom_type` Data_Unit: e.
- `lattice_const`: lattice constance of structure

# Example
```julia-repl
str = structure_fcc(atom_type=1, atom_charge=-1, lattice_const=3) # All atom in cell will have same type and charge
```
Or
```julia-repl
str = structure_fcc(atom_type=[1 2 3 4], atom_charge=[0 0 0 0], lattice_const=3)
```
"""
function structure_fcc(;atom_type=1, atom_charge=0, lattice_const=1)
    atom_vec = [
        0   0   0
        0   0.5 0.5
        0.5 0   0.5
        0.5 0.5 0
    ]
    cell_vec = diag([
        1   1   1
    ] .* lattice_const)
    atom_type, atom_charge, num_atoms, num_atom_types = transform(atom_vec, atom_type, atom_charge)
    res_name = "FCC"
    Structure_Basic(atom_vec, cell_vec, atom_type, atom_charge, num_atoms, num_atom_types, res_name)
end

"""
    structure_dc(;atom_type=1, atom_charge=0, lattice_const=1)

Do this will return a Structure_Basic type which contains all information needed to build a DC structure.

# Arguments
- `atom_type`: type of each atoms, could be a interger or a vector represent type of all atoms or each atoms repectively.
- `atom_charge`: charge of each atoms, could be a float or a vector, same as `atom_type` Data_Unit: e.
- `lattice_const`: lattice constance of structure

# Example
```julia-repl
str = structure_dc(atom_type=1, atom_charge=-1, lattice_const=3) # All atom in cell will have same type and charge
```
Or
```julia-repl
str = structure_dc(atom_type=[1 2 3 4 5 6 7 8], atom_charge=[0 0 0 0 0 0 0 0], lattice_const=3)
```
"""
function structure_dc(;atom_type=1, atom_charge=0, lattice_const=1)
    atom_vec = [
        0    0    0
        0    0.5  0.5
        0.5  0    0.5
        0.5  0.5  0
        0.25 0.25 0.25
        0.25 0.75 0.75
        0.75 0.25 0.75
        0.75 0.75 0.25
    ]
    cell_vec = diag([
        1   1   1
    ] .* lattice_const)
    atom_type, atom_charge, num_atoms, num_atom_types = transform(atom_vec, atom_type, atom_charge)
    res_name = "DC"
    Structure_Basic(atom_vec, cell_vec, atom_type, atom_charge, num_atoms, num_atom_types, res_name)
end


text:
# Wrapper for level-1 BLIS matrix routines.
#

# Level1v API of common forms are put together.
#
macro blis_group_level1m_form1(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          diagoffa, BliDoff,
                          diaga,    BliDiag,
                          uploa,    BliUpLo,
                          transa,   BliTrans,
                          m,        BliDim,
                          n,        BliDim,
                          a,        Ptr{xType}, rsa, BliInc, csa, BliInc,
                          b,        Ptr{xType}, rsb, BliInc, csb, BliInc)
    end
end

@blis_group_level1m_form1 addm
@blis_group_level1m_form1 copym
@blis_group_level1m_form1 subm

macro blis_group_level1m_form2(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          diagoffa, BliDoff,
                          diaga,    BliDiag,
                          uploa,    BliUpLo,
                          transa,   BliTrans,
                          m,        BliDim,
                          n,        BliDim,
                          α,        Ptr{xType},
                          a,        Ptr{xType}, rsa, BliInc, csa, BliInc,
                          b,        Ptr{xType}, rsb, BliInc, csb, BliInc)
    end
end

@blis_group_level1m_form2 axpym
@blis_group_level1m_form2 scal2m

macro blis_group_level1m_form3(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          conjα,    BliConj,
                          diagoffa, BliDoff,
                          diaga,    BliDiag,
                          uploa,    BliUpLo,
                          m,        BliDim,
                          n,        BliDim,
                          α,        Ptr{xType},
                          a,        Ptr{xType}, rsa, BliInc, csa, BliInc)
    end
end

@blis_group_level1m_form3 scalm
@blis_group_level1m_form3 setm


text:
function setinterval(lo::ArbFloat{P}, hi::ArbFloat{P}) where {P}
    lo > hi && return setinterval(hi, lo)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_interval_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}, Ref{ArbFloat}, Clong), z, lo, hi, P)
    return z
end

function setinterval(lo::ArbReal{P}, hi::ArbReal{P}) where {P}
    lo > hi && return setinterval(hi, lo)
    setinterval(lowerbound(lo, ArbFloat), upperbound(hi, ArbFloat))
end

function interval(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    ArbFloat{P}(lowerbound(x)), ArbFloat{P}(upperbound(x))
end

function interval(x::ArbReal{P}) where {P}
   lowerbound(x), upperbound(x)
end

function setball(mid::ArbFloat{P}, rad::ArbFloat{P}) where {P}
    signbit(rad) && throw(ErrorException("nonnegative radius required ($rad)"))
    lbound = mid - rad
    ubound = mid + rad
    setinterval(lbound, ubound)
end

function setball(mid::ArbReal{P}, rad::ArbReal{P}) where {P}
    signbit(rad) && throw(ErrorException("nonnegative radius required ($rad)"))
    setball(ArbFloat{P}(mid), ArbFloat{P}(rad))
end

function ball(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    ArbFloat{P}(midpoint(x)), ArbFloat{P}(radius(x))
end

function ball(x::ArbReal{P}) where {P}
    midpoint(x), radius(x)
end

function ball(x::ArbFloat{P}) where {P}
    x, zero(typeof(x))
end

@inline function upperbound(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_ubound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function upperbound(x::ArbReal{P}) where {P}
    w = upperbound(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

@inline function lowerbound(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_lbound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function lowerbound(x::ArbReal{P}) where {P}
    w = lowerbound(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

@inline function upperbound_abs(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_abs_ubound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function upperbound_abs(x::ArbReal{P}) where {P}
    w = upperbound_abs(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

@inline function lowerbound_abs(x::ArbReal{P}, ::Type{ArbFloat}) where {P}
    w = ArbFloat{P}()
    ccall(@libarb(arb_get_abs_lbound_arf), Cvoid, (Ref{ArbFloat}, Ref{ArbReal}, Clong), w, x, P)
    return w
end

function lowerbound_abs(x::ArbReal{P}) where {P}
    w = lowerbound(x, ArbFloat)
    z = ArbReal{P}()
    ccall(@libarb(arb_set_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), z, w)
    return z
end

function increase_radius(x::ArbReal{P}, err::ArbFloat{P}) where {P}
    err >= 0 || throw(ErrorException("nonnegative err required ($err)"))
    ccall(@libarb(arb_add_error_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), x, err)
    return x
end

function increase_radius(x::ArbReal{P}, err::ArbReal{P}) where {P}
    err >= 0 || throw(ErrorException("nonnegative err required ($err)"))
    ccall(@libarb(arb_add_error_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, err)
    return x
end

increase_radius(x::ArbReal{P}) where {P} = increase_radius(x, ulp(x))

function decrease_radius(x::ArbReal{P}, err::ArbFloat{P}) where {P}
    err = -abs(err)
    ccall(@libarb(arb_add_error_arf), Cvoid, (Ref{ArbReal}, Ref{ArbFloat}), x, err)
    return x
end

function decrease_radius(x::ArbReal{P}, err::ArbReal{P}) where {P}
    err = -abs(err)
    ccall(@libarb(arb_add_error_arb), Cvoid, (Ref{ArbReal}, Ref{ArbReal}), x, err)
    return x
end

decrease_radius(x::ArbReal{P}) where {P} = decrease_radius(x, ulp(x))