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

@test integrate(x^4*arccos(a*x), x) == :(-1//5*a^-5*(1+-1*a^2*x^2)^(1/2)+-1//25*a^-5*(1+-1*a^2*x^2)^5//2+1//5*x^5*arccos(a*x)+2//15*a^-5*(1+-1*a^2*x^2)^3//2)
@test integrate(x^3*arccos(a*x), x) == :(1//4*x^4*arccos(a*x)+3//32*a^-4*arcsin(a*x)+-3//32*x*a^-3*(1+-1*a^2*x^2)^(1/2)+-1//16*a^-1*x^3*(1+-1*a^2*x^2)^(1/2))
@test integrate(x^2*arccos(a*x), x) == :(-1//3*a^-3*(1+-1*a^2*x^2)^(1/2)+1//3*x^3*arccos(a*x)+1//9*a^-3*(1+-1*a^2*x^2)^3//2)
@test integrate(x*arccos(a*x), x) == :((1/2)*x^2*arccos(a*x)+1//4*a^-2*arcsin(a*x)+-1//4*x*a^-1*(1+-1*a^2*x^2)^(1/2))
@test integrate(arccos(a*x), x) == :(x*arccos(a*x)+-1*a^-1*(1+-1*a^2*x^2)^(1/2))
@test integrate(x^-1*arccos(a*x), x) == :(arccos(a*x)*log(1+exp(2*im*arccos(a*x)))+-1//2*im*arccos(a*x)^2+-1//2*im*Polylogarithms.polylog(2,-1*exp(2*im*arccos(a*x))))
@test integrate(x^-2*arccos(a*x), x) == :(a*arctanh((1+-1*a^2*x^2)^(1/2))+-1*x^-1*arccos(a*x))
@test integrate(x^-3*arccos(a*x), x) == :(-1//2*x^-2*arccos(a*x)+(1/2)*a*x^-1*(1+-1*a^2*x^2)^(1/2))
@test integrate(x^-4*arccos(a*x), x) == :(-1//3*x^-3*arccos(a*x)+1//6*a^3*arctanh((1+-1*a^2*x^2)^(1/2))+1//6*a*x^-2*(1+-1*a^2*x^2)^(1/2))
@test integrate(x^-5*arccos(a*x), x) == :(-1//4*x^-4*arccos(a*x)+1//6*a^3*x^-1*(1+-1*a^2*x^2)^(1/2)+1//12*a*x^-3*(1+-1*a^2*x^2)^(1/2))
@test integrate(x^-6*arccos(a*x), x) == :(-1//5*x^-5*arccos(a*x)+3//40*a^5*arctanh((1+-1*a^2*x^2)^(1/2))+1//20*a*x^-4*(1+-1*a^2*x^2)^(1/2)+3//40*a^3*x^-2*(1+-1*a^2*x^2)^(1/2))
@test integrate(x^4*arccos(a*x)^2, x) == :(-2//125*x^5+-16//75*x*a^-4+-8//225*a^-2*x^3+1//5*x^5*arccos(a*x)^2+-16//75*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+-8//75*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+-2//25*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^3*arccos(a*x)^2, x) == :(-1//32*x^4+-3//32*a^-4*arccos(a*x)^2+-3//32*a^-2*x^2+1//4*x^4*arccos(a*x)^2+-3//16*x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+-1//8*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^2*arccos(a*x)^2, x) == :(-2//27*x^3+-4//9*x*a^-2+1//3*x^3*arccos(a*x)^2+-4//9*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+-2//9*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x*arccos(a*x)^2, x) == :(-1//4*x^2+(1/2)*x^2*arccos(a*x)^2+-1//4*a^-2*arccos(a*x)^2+-1//2*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(arccos(a*x)^2, x) == :(-2x+x*arccos(a*x)^2+-2*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^-1*arccos(a*x)^2, x) == :((1/2)*Polylogarithms.polylog(3,-1*exp(2*im*arccos(a*x)))+arccos(a*x)^2*log(1+exp(2*im*arccos(a*x)))+-1//3*im*arccos(a*x)^3+-1*im*Polylogarithms.polylog(2,-1*exp(2*im*arccos(a*x)))*arccos(a*x))
@test integrate(x^-2*arccos(a*x)^2, x) == :(-1*x^-1*arccos(a*x)^2+-2*im*a*Polylogarithms.polylog(2,im*exp(im*arccos(a*x)))+2*im*a*Polylogarithms.polylog(2,-1*im*exp(im*arccos(a*x)))+-4*im*a*arccos(a*x)*arctan(exp(im*arccos(a*x))))
@test integrate(x^-3*arccos(a*x)^2, x) == :(a^2*log(x)+-1//2*x^-2*arccos(a*x)^2+a*x^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^-4*arccos(a*x)^2, x) == :(-1//3*a^2*x^-1+-1//3*x^-3*arccos(a*x)^2+-1//3*im*a^3*Polylogarithms.polylog(2,im*exp(im*arccos(a*x)))+1//3*im*a^3*Polylogarithms.polylog(2,-1*im*exp(im*arccos(a*x)))+-2//3*im*a^3*arccos(a*x)*arctan(exp(im*arccos(a*x)))+1//3*a*x^-2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^-5*arccos(a*x)^2, x) == :(-1//4*x^-4*arccos(a*x)^2+-1//12*a^2*x^-2+1//3*a^4*log(x)+1//3*a^3*x^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+1//6*a*x^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^4*arccos(a*x)^3, x) == :(-76//1125*a^-5*(1+-1*a^2*x^2)^3//2+-6//125*x^5*arccos(a*x)+1//5*x^5*arccos(a*x)^3+6//625*a^-5*(1+-1*a^2*x^2)^5//2+298//375*a^-5*(1+-1*a^2*x^2)^(1/2)+-16//25*x*a^-4*arccos(a*x)+-8//25*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2+-8//75*a^-2*x^3*arccos(a*x)+-4//25*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2+-3//25*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2)
@test integrate(x^3*arccos(a*x)^3, x) == :(-45//256*a^-4*arcsin(a*x)+-3//32*a^-4*arccos(a*x)^3+-3//32*x^4*arccos(a*x)+1//4*x^4*arccos(a*x)^3+-9//32*a^-2*x^2*arccos(a*x)+3//128*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)+45//256*x*a^-3*(1+-1*a^2*x^2)^(1/2)+-9//32*x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2+-3//16*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2)
@test integrate(x^2*arccos(a*x)^3, x) == :(-2//9*x^3*arccos(a*x)+-2//27*a^-3*(1+-1*a^2*x^2)^3//2+1//3*x^3*arccos(a*x)^3+14//9*a^-3*(1+-1*a^2*x^2)^(1/2)+-4//3*x*a^-2*arccos(a*x)+-2//3*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2+-1//3*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2)
@test integrate(x*arccos(a*x)^3, x) == :((1/2)*x^2*arccos(a*x)^3+-3//4*x^2*arccos(a*x)+-3//8*a^-2*arcsin(a*x)+-1//4*a^-2*arccos(a*x)^3+3//8*x*a^-1*(1+-1*a^2*x^2)^(1/2)+-3//4*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2)
@test integrate(arccos(a*x)^3, x) == :(x*arccos(a*x)^3+-6*x*arccos(a*x)+6*a^-1*(1+-1*a^2*x^2)^(1/2)+-3*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2)
@test integrate(x^-1*arccos(a*x)^3, x) == :(arccos(a*x)^3*log(1+exp(2*im*arccos(a*x)))+-1//4*im*arccos(a*x)^4+3//2*Polylogarithms.polylog(3,-1*exp(2*im*arccos(a*x)))*arccos(a*x)+3//4*im*Polylogarithms.polylog(4,-1*exp(2*im*arccos(a*x)))+-3//2*im*arccos(a*x)^2*Polylogarithms.polylog(2,-1*exp(2*im*arccos(a*x))))
@test integrate(x^-2*arccos(a*x)^3, x) == :(-1*x^-1*arccos(a*x)^3+-6*a*Polylogarithms.polylog(3,-1*im*exp(im*arccos(a*x)))+6*a*Polylogarithms.polylog(3,im*exp(im*arccos(a*x)))+-6*im*a*arccos(a*x)^2*arctan(exp(im*arccos(a*x)))+-6*im*a*Polylogarithms.polylog(2,im*exp(im*arccos(a*x)))*arccos(a*x)+6*im*a*Polylogarithms.polylog(2,-1*im*exp(im*arccos(a*x)))*arccos(a*x))
@test integrate(x^-3*arccos(a*x)^3, x) == :(-1//2*x^-2*arccos(a*x)^3+3*a^2*arccos(a*x)*log(1+exp(2*im*arccos(a*x)))+-3//2*im*a^2*arccos(a*x)^2+-3//2*im*a^2*Polylogarithms.polylog(2,-1*exp(2*im*arccos(a*x)))+3//2*a*x^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2)
@test integrate(x^-4*arccos(a*x)^3, x) == :(a^3*Polylogarithms.polylog(3,im*exp(im*arccos(a*x)))+a^3*arctanh((1+-1*a^2*x^2)^(1/2))+-1*a^3*Polylogarithms.polylog(3,-1*im*exp(im*arccos(a*x)))+-1//3*x^-3*arccos(a*x)^3+-1*a^2*x^-1*arccos(a*x)+im*a^3*Polylogarithms.polylog(2,-1*im*exp(im*arccos(a*x)))*arccos(a*x)+(1/2)*a*x^-2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2+-1*im*a^3*arccos(a*x)^2*arctan(exp(im*arccos(a*x)))+-1*im*a^3*Polylogarithms.polylog(2,im*exp(im*arccos(a*x)))*arccos(a*x))
@test integrate(x^-5*arccos(a*x)^3, x) == :(-1//4*x^-4*arccos(a*x)^3+a^4*arccos(a*x)*log(1+exp(2*im*arccos(a*x)))+-1//2*im*a^4*arccos(a*x)^2+-1//2*im*a^4*Polylogarithms.polylog(2,-1*exp(2*im*arccos(a*x)))+-1//4*a^2*x^-2*arccos(a*x)+1//4*a^3*x^-1*(1+-1*a^2*x^2)^(1/2)+(1/2)*a^3*x^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2+1//4*a*x^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^2)
@test integrate(x^5*arccos(a*x)^4, x) == :(1//324*x^6+-5//96*a^-6*arccos(a*x)^4+-1//18*x^6*arccos(a*x)^2+1//6*x^6*arccos(a*x)^4+65//3456*a^-2*x^4+245//1152*a^-6*arccos(a*x)^2+245//1152*a^-4*x^2+-5//16*a^-4*x^2*arccos(a*x)^2+-5//48*a^-2*x^4*arccos(a*x)^2+-5//24*x*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+-5//36*a^-3*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+-1//9*a^-1*x^5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+1//54*a^-1*x^5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+65//864*a^-3*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+245//576*x*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^4*arccos(a*x)^4, x) == :(24//3125*x^5+-12//125*x^5*arccos(a*x)^2+1//5*x^5*arccos(a*x)^4+1088//16875*a^-2*x^3+16576//5625*x*a^-4+-32//25*x*a^-4*arccos(a*x)^2+-32//75*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+-16//75*a^-2*x^3*arccos(a*x)^2+16576//5625*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+-16//75*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+-4//25*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+24//625*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+1088//5625*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^3*arccos(a*x)^4, x) == :(3//128*x^4+-3//16*x^4*arccos(a*x)^2+-3//32*a^-4*arccos(a*x)^4+1//4*x^4*arccos(a*x)^4+45//128*a^-4*arccos(a*x)^2+45//128*a^-2*x^2+-9//16*a^-2*x^2*arccos(a*x)^2+-3//8*x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+-1//4*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+3//32*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+45//64*x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^2*arccos(a*x)^4, x) == :(8//81*x^3+-4//9*x^3*arccos(a*x)^2+1//3*x^3*arccos(a*x)^4+160//27*x*a^-2+-8//3*x*a^-2*arccos(a*x)^2+-8//9*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+160//27*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)+-4//9*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+8//27*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x*arccos(a*x)^4, x) == :(3//4*x^2+(1/2)*x^2*arccos(a*x)^4+-3//2*x^2*arccos(a*x)^2+-1//4*a^-2*arccos(a*x)^4+3//4*a^-2*arccos(a*x)^2+-1*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+3//2*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(arccos(a*x)^4, x) == :(24x+x*arccos(a*x)^4+-12*x*arccos(a*x)^2+-4*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3+24*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x))
@test integrate(x^-1*arccos(a*x)^4, x) == :(-3//2*Polylogarithms.polylog(5,-1*exp(2*im*arccos(a*x)))+arccos(a*x)^4*log(1+exp(2*im*arccos(a*x)))+3*arccos(a*x)^2*Polylogarithms.polylog(3,-1*exp(2*im*arccos(a*x)))+-1//5*im*arccos(a*x)^5+-2*im*arccos(a*x)^3*Polylogarithms.polylog(2,-1*exp(2*im*arccos(a*x)))+3*im*Polylogarithms.polylog(4,-1*exp(2*im*arccos(a*x)))*arccos(a*x))
@test integrate(x^-2*arccos(a*x)^4, x) == :(-1*x^-1*arccos(a*x)^4+-24*im*a*Polylogarithms.polylog(4,-1*im*exp(im*arccos(a*x)))+-24*a*Polylogarithms.polylog(3,-1*im*exp(im*arccos(a*x)))*arccos(a*x)+24*im*a*Polylogarithms.polylog(4,im*exp(im*arccos(a*x)))+24*a*Polylogarithms.polylog(3,im*exp(im*arccos(a*x)))*arccos(a*x)+-12*im*a*arccos(a*x)^2*Polylogarithms.polylog(2,im*exp(im*arccos(a*x)))+-8*im*a*arccos(a*x)^3*arctan(exp(im*arccos(a*x)))+12*im*a*arccos(a*x)^2*Polylogarithms.polylog(2,-1*im*exp(im*arccos(a*x))))
@test integrate(x^-3*arccos(a*x)^4, x) == :(3*a^2*Polylogarithms.polylog(3,-1*exp(2*im*arccos(a*x)))+-1//2*x^-2*arccos(a*x)^4+-2*im*a^2*arccos(a*x)^3+6*a^2*arccos(a*x)^2*log(1+exp(2*im*arccos(a*x)))+-6*im*a^2*Polylogarithms.polylog(2,-1*exp(2*im*arccos(a*x)))*arccos(a*x)+2*a*x^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3)
@test integrate(x^-4*arccos(a*x)^4, x) == :(-1//3*x^-3*arccos(a*x)^4+-4*im*a^3*Polylogarithms.polylog(2,im*exp(im*arccos(a*x)))+-4*im*a^3*Polylogarithms.polylog(4,-1*im*exp(im*arccos(a*x)))+-4*a^3*Polylogarithms.polylog(3,-1*im*exp(im*arccos(a*x)))*arccos(a*x)+-2*a^2*x^-1*arccos(a*x)^2+4*im*a^3*Polylogarithms.polylog(2,-1*im*exp(im*arccos(a*x)))+4*im*a^3*Polylogarithms.polylog(4,im*exp(im*arccos(a*x)))+4*a^3*Polylogarithms.polylog(3,im*exp(im*arccos(a*x)))*arccos(a*x)+-8*im*a^3*arccos(a*x)*arctan(exp(im*arccos(a*x)))+-2*im*a^3*arccos(a*x)^2*Polylogarithms.polylog(2,im*exp(im*arccos(a*x)))+2*im*a^3*arccos(a*x)^2*Polylogarithms.polylog(2,-1*im*exp(im*arccos(a*x)))+-4//3*im*a^3*arccos(a*x)^3*arctan(exp(im*arccos(a*x)))+2//3*a*x^-2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3)
@test integrate(x^6*arccos(a*x)^-1, x) == :(-9//64*a^-7*sinintegral(3*arccos(a*x))+-5//64*a^-7*sinintegral(5*arccos(a*x))+-5//64*a^-7*sinintegral(arccos(a*x))+-1//64*a^-7*sinintegral(7*arccos(a*x)))
@test integrate(x^5*arccos(a*x)^-1, x) == :(-5//32*a^-6*sinintegral(2*arccos(a*x))+-1//8*a^-6*sinintegral(4*arccos(a*x))+-1//32*a^-6*sinintegral(6*arccos(a*x)))
@test integrate(x^4*arccos(a*x)^-1, x) == :(-3//16*a^-5*sinintegral(3*arccos(a*x))+-1//8*a^-5*sinintegral(arccos(a*x))+-1//16*a^-5*sinintegral(5*arccos(a*x)))
@test integrate(x^3*arccos(a*x)^-1, x) == :(-1//4*a^-4*sinintegral(2*arccos(a*x))+-1//8*a^-4*sinintegral(4*arccos(a*x)))
@test integrate(x^2*arccos(a*x)^-1, x) == :(-1//4*a^-3*sinintegral(3*arccos(a*x))+-1//4*a^-3*sinintegral(arccos(a*x)))
@test integrate(x*arccos(a*x)^-1, x) == :(-1//2*a^-2*sinintegral(2*arccos(a*x)))
@test integrate(arccos(a*x)^-1, x) == :(-1*a^-1*sinintegral(arccos(a*x)))
@test integrate(x^6*arccos(a*x)^-2, x) == :(-27//64*a^-7*cosintegral(3*arccos(a*x))+-25//64*a^-7*cosintegral(5*arccos(a*x))+-7//64*a^-7*cosintegral(7*arccos(a*x))+-5//64*a^-7*cosintegral(arccos(a*x))+a^-1*x^6*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1)
@test integrate(x^5*arccos(a*x)^-2, x) == :(-5//16*a^-6*cosintegral(2*arccos(a*x))+-3//16*a^-6*cosintegral(6*arccos(a*x))+-1//2*a^-6*cosintegral(4*arccos(a*x))+a^-1*x^5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1)
@test integrate(x^4*arccos(a*x)^-2, x) == :(-9//16*a^-5*cosintegral(3*arccos(a*x))+-5//16*a^-5*cosintegral(5*arccos(a*x))+-1//8*a^-5*cosintegral(arccos(a*x))+a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1)
@test integrate(x^3*arccos(a*x)^-2, x) == :(-1//2*a^-4*cosintegral(2*arccos(a*x))+-1//2*a^-4*cosintegral(4*arccos(a*x))+a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1)
@test integrate(x^2*arccos(a*x)^-2, x) == :(-3//4*a^-3*cosintegral(3*arccos(a*x))+-1//4*a^-3*cosintegral(arccos(a*x))+a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1)
@test integrate(x*arccos(a*x)^-2, x) == :(-1*a^-2*cosintegral(2*arccos(a*x))+x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1)
@test integrate(arccos(a*x)^-2, x) == :(-1*a^-1*cosintegral(arccos(a*x))+a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1)
@test integrate(x^4*arccos(a*x)^-3, x) == :(1//16*a^-5*sinintegral(arccos(a*x))+5//2*x^5*arccos(a*x)^-1+25//32*a^-5*sinintegral(5*arccos(a*x))+27//32*a^-5*sinintegral(3*arccos(a*x))+-2*a^-2*x^3*arccos(a*x)^-1+(1/2)*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-2)
@test integrate(x^3*arccos(a*x)^-3, x) == :(a^-4*sinintegral(4*arccos(a*x))+(1/2)*a^-4*sinintegral(2*arccos(a*x))+2*x^4*arccos(a*x)^-1+-3//2*a^-2*x^2*arccos(a*x)^-1+(1/2)*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-2)
@test integrate(x^2*arccos(a*x)^-3, x) == :(1//8*a^-3*sinintegral(arccos(a*x))+3//2*x^3*arccos(a*x)^-1+9//8*a^-3*sinintegral(3*arccos(a*x))+-1*x*a^-2*arccos(a*x)^-1+(1/2)*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-2)
@test integrate(x*arccos(a*x)^-3, x) == :(a^-2*sinintegral(2*arccos(a*x))+x^2*arccos(a*x)^-1+-1//2*a^-2*arccos(a*x)^-1+(1/2)*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-2)
@test integrate(arccos(a*x)^-3, x) == :((1/2)*x*arccos(a*x)^-1+(1/2)*a^-1*sinintegral(arccos(a*x))+(1/2)*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-2)
@test integrate(x^4*arccos(a*x)^-4, x) == :(1//48*a^-5*cosintegral(arccos(a*x))+5//6*x^5*arccos(a*x)^-2+27//32*a^-5*cosintegral(3*arccos(a*x))+125//96*a^-5*cosintegral(5*arccos(a*x))+-2//3*a^-2*x^3*arccos(a*x)^-2+2*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+-25//6*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+1//3*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3)
@test integrate(x^3*arccos(a*x)^-4, x) == :(1//3*a^-4*cosintegral(2*arccos(a*x))+2//3*x^4*arccos(a*x)^-2+4//3*a^-4*cosintegral(4*arccos(a*x))+-1//2*a^-2*x^2*arccos(a*x)^-2+x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+-8//3*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+1//3*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3)
@test integrate(x^2*arccos(a*x)^-4, x) == :((1/2)*x^3*arccos(a*x)^-2+1//24*a^-3*cosintegral(arccos(a*x))+9//8*a^-3*cosintegral(3*arccos(a*x))+-1//3*x*a^-2*arccos(a*x)^-2+1//3*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+-3//2*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+1//3*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3)
@test integrate(x*arccos(a*x)^-4, x) == :(-1//6*a^-2*arccos(a*x)^-2+1//3*x^2*arccos(a*x)^-2+2//3*a^-2*cosintegral(2*arccos(a*x))+-2//3*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+1//3*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3)
@test integrate(arccos(a*x)^-4, x) == :(1//6*x*arccos(a*x)^-2+1//6*a^-1*cosintegral(arccos(a*x))+-1//6*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1+1//3*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3)
@test integrate(x^4*arccos(a*x)^(1/2), x) == :(1//5*x^5*arccos(a*x)^(1/2)+-1//16*2^(1/2)*pi^(1/2)*a^-5*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//96*6^(1/2)*pi^(1/2)*a^-5*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//800*10^(1/2)*pi^(1/2)*a^-5*fresnelc(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^3*arccos(a*x)^(1/2), x) == :(-3//32*a^-4*arccos(a*x)^(1/2)+1//4*x^4*arccos(a*x)^(1/2)+-1//16*pi^(1/2)*a^-4*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+-1//128*2^(1/2)*pi^(1/2)*a^-4*fresnelc(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^2*arccos(a*x)^(1/2), x) == :(1//3*x^3*arccos(a*x)^(1/2)+-1//8*2^(1/2)*pi^(1/2)*a^-3*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//72*6^(1/2)*pi^(1/2)*a^-3*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x*arccos(a*x)^(1/2), x) == :((1/2)*x^2*arccos(a*x)^(1/2)+-1//4*a^-2*arccos(a*x)^(1/2)+-1//8*pi^(1/2)*a^-2*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(arccos(a*x)^(1/2), x) == :(x*arccos(a*x)^(1/2)+-1//2*2^(1/2)*pi^(1/2)*a^-1*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^4*arccos(a*x)^3//2, x) == :(1//5*x^5*arccos(a*x)^3//2+-4//25*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+-3//50*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+-2//25*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+1//192*6^(1/2)*pi^(1/2)*a^-5*fresnels(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+3//32*2^(1/2)*pi^(1/2)*a^-5*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+3//8000*10^(1/2)*pi^(1/2)*a^-5*fresnels(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^3*arccos(a*x)^3//2, x) == :(-3//32*a^-4*arccos(a*x)^3//2+1//4*x^4*arccos(a*x)^3//2+3//64*pi^(1/2)*a^-4*fresnels(2*pi^-1//2*arccos(a*x)^(1/2))+-9//64*x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+-3//32*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+3//1024*2^(1/2)*pi^(1/2)*a^-4*fresnels(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^2*arccos(a*x)^3//2, x) == :(1//3*x^3*arccos(a*x)^3//2+-1//3*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+-1//6*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+1//144*6^(1/2)*pi^(1/2)*a^-3*fresnels(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+3//16*2^(1/2)*pi^(1/2)*a^-3*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x*arccos(a*x)^3//2, x) == :((1/2)*x^2*arccos(a*x)^3//2+-1//4*a^-2*arccos(a*x)^3//2+3//32*pi^(1/2)*a^-2*fresnels(2*pi^-1//2*arccos(a*x)^(1/2))+-3//8*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2))
@test integrate(arccos(a*x)^3//2, x) == :(x*arccos(a*x)^3//2+-3//2*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^(1/2)+3//4*2^(1/2)*pi^(1/2)*a^-1*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^4*arccos(a*x)^5//2, x) == :(-3//100*x^5*arccos(a*x)^(1/2)+1//5*x^5*arccos(a*x)^5//2+-4//15*a^-5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+-2//5*x*a^-4*arccos(a*x)^(1/2)+-1//15*a^-2*x^3*arccos(a*x)^(1/2)+-2//15*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+-1//10*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+3//16000*10^(1/2)*pi^(1/2)*a^-5*fresnelc(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+5//1152*6^(1/2)*pi^(1/2)*a^-5*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+15//64*2^(1/2)*pi^(1/2)*a^-5*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^3*arccos(a*x)^5//2, x) == :(-15//256*x^4*arccos(a*x)^(1/2)+-3//32*a^-4*arccos(a*x)^5//2+1//4*x^4*arccos(a*x)^5//2+225//2048*a^-4*arccos(a*x)^(1/2)+-45//256*a^-2*x^2*arccos(a*x)^(1/2)+15//256*pi^(1/2)*a^-4*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+-15//64*x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+-5//32*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+15//8192*2^(1/2)*pi^(1/2)*a^-4*fresnelc(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^2*arccos(a*x)^5//2, x) == :(-5//36*x^3*arccos(a*x)^(1/2)+1//3*x^3*arccos(a*x)^5//2+-5//6*x*a^-2*arccos(a*x)^(1/2)+-5//9*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+-5//18*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+5//864*6^(1/2)*pi^(1/2)*a^-3*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+15//32*2^(1/2)*pi^(1/2)*a^-3*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x*arccos(a*x)^5//2, x) == :((1/2)*x^2*arccos(a*x)^5//2+-15//32*x^2*arccos(a*x)^(1/2)+-1//4*a^-2*arccos(a*x)^5//2+15//64*a^-2*arccos(a*x)^(1/2)+15//128*pi^(1/2)*a^-2*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+-5//8*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2)
@test integrate(arccos(a*x)^5//2, x) == :(x*arccos(a*x)^5//2+-15//4*x*arccos(a*x)^(1/2)+-5//2*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^3//2+15//8*2^(1/2)*pi^(1/2)*a^-1*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^4*arccos(a*x)^-1//2, x) == :(-1//8*2^(1/2)*pi^(1/2)*a^-5*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//16*6^(1/2)*pi^(1/2)*a^-5*fresnels(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//80*10^(1/2)*pi^(1/2)*a^-5*fresnels(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^3*arccos(a*x)^-1//2, x) == :(-1//4*pi^(1/2)*a^-4*fresnels(2*pi^-1//2*arccos(a*x)^(1/2))+-1//16*2^(1/2)*pi^(1/2)*a^-4*fresnels(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^2*arccos(a*x)^-1//2, x) == :(-1//4*2^(1/2)*pi^(1/2)*a^-3*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//12*6^(1/2)*pi^(1/2)*a^-3*fresnels(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x*arccos(a*x)^-1//2, x) == :(-1//2*pi^(1/2)*a^-2*fresnels(2*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(arccos(a*x)^-1//2, x) == :(-1*2^(1/2)*pi^(1/2)*a^-1*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^6*arccos(a*x)^-3//2, x) == :(2*a^-1*x^6*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+-9//32*6^(1/2)*pi^(1/2)*a^-7*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-5//32*2^(1/2)*pi^(1/2)*a^-7*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-5//32*10^(1/2)*pi^(1/2)*a^-7*fresnelc(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//32*14^(1/2)*pi^(1/2)*a^-7*fresnelc(14^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^5*arccos(a*x)^-3//2, x) == :(-5//8*pi^(1/2)*a^-6*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+2*a^-1*x^5*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+-1//2*2^(1/2)*pi^(1/2)*a^-6*fresnelc(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//8*3^(1/2)*pi^(1/2)*a^-6*fresnelc(2*3^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^4*arccos(a*x)^-3//2, x) == :(2*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+-3//8*6^(1/2)*pi^(1/2)*a^-5*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//4*2^(1/2)*pi^(1/2)*a^-5*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//8*10^(1/2)*pi^(1/2)*a^-5*fresnelc(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^3*arccos(a*x)^-3//2, x) == :(-1*pi^(1/2)*a^-4*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+2*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+-1//2*2^(1/2)*pi^(1/2)*a^-4*fresnelc(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^2*arccos(a*x)^-3//2, x) == :(2*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+-1//2*2^(1/2)*pi^(1/2)*a^-3*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+-1//2*6^(1/2)*pi^(1/2)*a^-3*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x*arccos(a*x)^-3//2, x) == :(-2*pi^(1/2)*a^-2*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+2*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2)
@test integrate(arccos(a*x)^-3//2, x) == :(2*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+-2*2^(1/2)*pi^(1/2)*a^-1*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^4*arccos(a*x)^-5//2, x) == :(20//3*x^5*arccos(a*x)^-1//2+-16//3*a^-2*x^3*arccos(a*x)^-1//2+1//6*2^(1/2)*pi^(1/2)*a^-5*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+2//3*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3//2+3//4*6^(1/2)*pi^(1/2)*a^-5*fresnels(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+5//12*10^(1/2)*pi^(1/2)*a^-5*fresnels(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^3*arccos(a*x)^-5//2, x) == :(16//3*x^4*arccos(a*x)^-1//2+-4*a^-2*x^2*arccos(a*x)^-1//2+4//3*pi^(1/2)*a^-4*fresnels(2*pi^-1//2*arccos(a*x)^(1/2))+2//3*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3//2+4//3*2^(1/2)*pi^(1/2)*a^-4*fresnels(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^2*arccos(a*x)^-5//2, x) == :(4*x^3*arccos(a*x)^-1//2+-8//3*x*a^-2*arccos(a*x)^-1//2+6^(1/2)*pi^(1/2)*a^-3*fresnels(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+1//3*2^(1/2)*pi^(1/2)*a^-3*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+2//3*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3//2)
@test integrate(x*arccos(a*x)^-5//2, x) == :(-4//3*a^-2*arccos(a*x)^-1//2+8//3*x^2*arccos(a*x)^-1//2+8//3*pi^(1/2)*a^-2*fresnels(2*pi^-1//2*arccos(a*x)^(1/2))+2//3*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3//2)
@test integrate(arccos(a*x)^-5//2, x) == :(4//3*x*arccos(a*x)^-1//2+2//3*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-3//2+4//3*2^(1/2)*pi^(1/2)*a^-1*fresnels(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^4*arccos(a*x)^-7//2, x) == :(4//3*x^5*arccos(a*x)^-3//2+-16//15*a^-2*x^3*arccos(a*x)^-3//2+-40//3*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+1//15*2^(1/2)*pi^(1/2)*a^-5*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+2//5*a^-1*x^4*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-5//2+5//6*10^(1/2)*pi^(1/2)*a^-5*fresnelc(10^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+9//10*6^(1/2)*pi^(1/2)*a^-5*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+32//5*a^-3*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2)
@test integrate(x^3*arccos(a*x)^-7//2, x) == :(16//15*x^4*arccos(a*x)^-3//2+-4//5*a^-2*x^2*arccos(a*x)^-3//2+16//15*pi^(1/2)*a^-4*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+-128//15*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+2//5*a^-1*x^3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-5//2+16//5*x*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+32//15*2^(1/2)*pi^(1/2)*a^-4*fresnelc(2*2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^2*arccos(a*x)^-7//2, x) == :(4//5*x^3*arccos(a*x)^-3//2+-8//15*x*a^-2*arccos(a*x)^-3//2+16//15*a^-3*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+-24//5*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+2//5*a^-1*x^2*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-5//2+2//15*2^(1/2)*pi^(1/2)*a^-3*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2))+6//5*6^(1/2)*pi^(1/2)*a^-3*fresnelc(6^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x*arccos(a*x)^-7//2, x) == :(-4//15*a^-2*arccos(a*x)^-3//2+8//15*x^2*arccos(a*x)^-3//2+32//15*pi^(1/2)*a^-2*fresnelc(2*pi^-1//2*arccos(a*x)^(1/2))+-32//15*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+2//5*x*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-5//2)
@test integrate(arccos(a*x)^-7//2, x) == :(4//15*x*arccos(a*x)^-3//2+-8//15*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-1//2+2//5*a^-1*(1+-1*a^2*x^2)^(1/2)*arccos(a*x)^-5//2+8//15*2^(1/2)*pi^(1/2)*a^-1*fresnelc(2^(1/2)*pi^-1//2*arccos(a*x)^(1/2)))
@test integrate(x^3*arccos(a*x)^n, x) == :(2^(-6+-2n)*a^-4*(im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,4*im*arccos(a*x))+2^(-6+-2n)*a^-4*(-1*im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,-4*im*arccos(a*x))+2^(-4+-1n)*a^-4*(im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,2*im*arccos(a*x))+2^(-4+-1n)*a^-4*(-1*im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,-2*im*arccos(a*x)))
@test integrate(x^2*arccos(a*x)^n, x) == :(1//8*a^-3*(im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,im*arccos(a*x))+1//8*a^-3*(-1*im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,-1*im*arccos(a*x))+1//8*3^(-1+-1n)*a^-3*(im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,3*im*arccos(a*x))+1//8*3^(-1+-1n)*a^-3*(-1*im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,-3*im*arccos(a*x)))
@test integrate(x*arccos(a*x)^n, x) == :(2^(-3+-1n)*a^-2*(im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,2*im*arccos(a*x))+2^(-3+-1n)*a^-2*(-1*im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,-2*im*arccos(a*x)))
@test integrate(arccos(a*x)^n, x) == :((1/2)*a^-1*(im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,im*arccos(a*x))+(1/2)*a^-1*(-1*im*arccos(a*x))^(-1n)*arccos(a*x)^n*SpecialFunctions.gamma(1+n,-1*im*arccos(a*x)))
@test integrate(x^3*(a+b*arccos(c*x)), x) == :(1//4*x^4*(a+b*arccos(c*x))+3//32*b*c^-4*arcsin(c*x)+-3//32*b*x*c^-3*(1+-1*c^2*x^2)^(1/2)+-1//16*b*c^-1*x^3*(1+-1*c^2*x^2)^(1/2))
@test integrate(x^2*(a+b*arccos(c*x)), x) == :(1//3*x^3*(a+b*arccos(c*x))+-1//3*b*c^-3*(1+-1*c^2*x^2)^(1/2)+1//9*b*c^-3*(1+-1*c^2*x^2)^3//2)
@test integrate(x*(a+b*arccos(c*x)), x) == :((1/2)*x^2*(a+b*arccos(c*x))+1//4*b*c^-2*arcsin(c*x)+-1//4*b*x*c^-1*(1+-1*c^2*x^2)^(1/2))
@test integrate(a+b*arccos(c*x), x) == :(a*x+b*x*arccos(c*x)+-1*b*c^-1*(1+-1*c^2*x^2)^(1/2))
@test integrate(x^-1*(a+b*arccos(c*x)), x) == :((a+b*arccos(c*x))*log(1+exp(2*im*arccos(c*x)))+-1//2*im*b*Polylogarithms.polylog(2,-1*exp(2*im*arccos(c*x)))+-1//2*im*b^-1*(a+b*arccos(c*x))^2)
@test integrate(x^-2*(a+b*arccos(c*x)), x) == :(-1*x^-1*(a+b*arccos(c*x))+b*c*arctanh((1+-1*c^2*x^2)^(1/2)))
@test integrate(x^-3*(a+b*arccos(c*x)), x) == :(-1//2*x^-2*(a+b*arccos(c*x))+(1/2)*b*c*x^-1*(1+-1*c^2*x^2)^(1/2))
@test integrate(x^-4*(a+b*arccos(c*x)), x) == :(-1//3*x^-3*(a+b*arccos(c*x))+1//6*b*c^3*arctanh((1+-1*c^2*x^2)^(1/2))+1//6*b*c*x^-2*(1+-1*c^2*x^2)^(1/2))
@test integrate(x^2*(a+b*arccos(c*x))^2, x) == :(-2//27*b^2*x^3+1//3*x^3*(a+b*arccos(c*x))^2+-4//9*x*b^2*c^-2+-4//9*b*c^-3*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))+-2//9*b*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x)))
@test integrate(x*(a+b*arccos(c*x))^2, x) == :((1/2)*x^2*(a+b*arccos(c*x))^2+-1//4*b^2*x^2+-1//4*c^-2*(a+b*arccos(c*x))^2+-1//2*b*x*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x)))
@test integrate((a+b*arccos(c*x))^2, x) == :(x*(a+b*arccos(c*x))^2+-2*x*b^2+-2*b*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x)))
@test integrate(x^-1*(a+b*arccos(c*x))^2, x) == :((a+b*arccos(c*x))^2*log(1+exp(2*im*arccos(c*x)))+(1/2)*b^2*Polylogarithms.polylog(3,-1*exp(2*im*arccos(c*x)))+-1//3*im*b^-1*(a+b*arccos(c*x))^3+-1*im*b*(a+b*arccos(c*x))*Polylogarithms.polylog(2,-1*exp(2*im*arccos(c*x))))
@test integrate(x^-2*(a+b*arccos(c*x))^2, x) == :(-1*x^-1*(a+b*arccos(c*x))^2+-2*im*c*b^2*Polylogarithms.polylog(2,im*exp(im*arccos(c*x)))+2*im*c*b^2*Polylogarithms.polylog(2,-1*im*exp(im*arccos(c*x)))+-4*im*b*c*(a+b*arccos(c*x))*arctan(exp(im*arccos(c*x))))
@test integrate(x^2*(a+b*arccos(c*x))^3, x) == :(1//3*x^3*(a+b*arccos(c*x))^3+-2//9*b^2*x^3*(a+b*arccos(c*x))+-2//27*b^3*c^-3*(1+-1*c^2*x^2)^3//2+14//9*b^3*c^-3*(1+-1*c^2*x^2)^(1/2)+-4//3*a*x*b^2*c^-2+-4//3*x*b^3*c^-2*arccos(c*x)+-2//3*b*c^-3*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^2+-1//3*b*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^2)
@test integrate(x*(a+b*arccos(c*x))^3, x) == :((1/2)*x^2*(a+b*arccos(c*x))^3+-1//4*c^-2*(a+b*arccos(c*x))^3+-3//4*b^2*x^2*(a+b*arccos(c*x))+-3//8*b^3*c^-2*arcsin(c*x)+3//8*x*b^3*c^-1*(1+-1*c^2*x^2)^(1/2)+-3//4*b*x*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^2)
@test integrate((a+b*arccos(c*x))^3, x) == :(x*(a+b*arccos(c*x))^3+-6*a*x*b^2+-6*x*b^3*arccos(c*x)+6*b^3*c^-1*(1+-1*c^2*x^2)^(1/2)+-3*b*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^2)
@test integrate(x^-1*(a+b*arccos(c*x))^3, x) == :((a+b*arccos(c*x))^3*log(1+exp(2*im*arccos(c*x)))+-1//4*im*b^-1*(a+b*arccos(c*x))^4+3//2*b^2*(a+b*arccos(c*x))*Polylogarithms.polylog(3,-1*exp(2*im*arccos(c*x)))+3//4*im*b^3*Polylogarithms.polylog(4,-1*exp(2*im*arccos(c*x)))+-3//2*im*b*(a+b*arccos(c*x))^2*Polylogarithms.polylog(2,-1*exp(2*im*arccos(c*x))))
@test integrate(x^-2*(a+b*arccos(c*x))^3, x) == :(-1*x^-1*(a+b*arccos(c*x))^3+-6*c*b^3*Polylogarithms.polylog(3,-1*im*exp(im*arccos(c*x)))+6*c*b^3*Polylogarithms.polylog(3,im*exp(im*arccos(c*x)))+-6*im*b*c*(a+b*arccos(c*x))^2*arctan(exp(im*arccos(c*x)))+-6*im*c*b^2*(a+b*arccos(c*x))*Polylogarithms.polylog(2,im*exp(im*arccos(c*x)))+6*im*c*b^2*(a+b*arccos(c*x))*Polylogarithms.polylog(2,-1*im*exp(im*arccos(c*x))))
@test integrate(x^2*(a+b*arccos(c*x))^-1, x) == :(-1//4*b^-1*c^-3*cos(a*b^-1)*sinintegral(a*b^-1+arccos(c*x))+-1//4*b^-1*c^-3*cos(3*a*b^-1)*sinintegral(3*arccos(c*x)+3*a*b^-1)+1//4*b^-1*c^-3*cosintegral(3*arccos(c*x)+3*a*b^-1)*sin(3*a*b^-1)+1//4*b^-1*c^-3*cosintegral(a*b^-1+arccos(c*x))*sin(a*b^-1))
@test integrate(x*(a+b*arccos(c*x))^-1, x) == :((1/2)*b^-1*c^-2*cosintegral(2*arccos(c*x)+2*a*b^-1)*sin(2*a*b^-1)+-1//2*b^-1*c^-2*cos(2*a*b^-1)*sinintegral(2*arccos(c*x)+2*a*b^-1))
@test integrate((a+b*arccos(c*x))^-1, x) == :(b^-1*c^-1*cosintegral(b^-1*(a+b*arccos(c*x)))*sin(a*b^-1)+-1*b^-1*c^-1*cos(a*b^-1)*sinintegral(b^-1*(a+b*arccos(c*x))))
@test integrate(x^2*(a+b*arccos(c*x))^-2, x) == :(-3//4*b^-2*c^-3*cos(3*a*b^-1)*cosintegral(3*arccos(c*x)+3*a*b^-1)+-3//4*b^-2*c^-3*sin(3*a*b^-1)*sinintegral(3*arccos(c*x)+3*a*b^-1)+-1//4*b^-2*c^-3*cos(a*b^-1)*cosintegral(a*b^-1+arccos(c*x))+-1//4*b^-2*c^-3*sin(a*b^-1)*sinintegral(a*b^-1+arccos(c*x))+b^-1*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-1)
@test integrate(x*(a+b*arccos(c*x))^-2, x) == :(-1*b^-2*c^-2*cos(2*a*b^-1)*cosintegral(2*arccos(c*x)+2*a*b^-1)+-1*b^-2*c^-2*sin(2*a*b^-1)*sinintegral(2*arccos(c*x)+2*a*b^-1)+x*b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-1)
@test integrate((a+b*arccos(c*x))^-2, x) == :(b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-1+-1*b^-2*c^-1*cos(a*b^-1)*cosintegral(a*b^-1+arccos(c*x))+-1*b^-2*c^-1*sin(a*b^-1)*sinintegral(a*b^-1+arccos(c*x)))
@test integrate(x^2*(a+b*arccos(c*x))^-3, x) == :(3//2*b^-2*x^3*(a+b*arccos(c*x))^-1+b^-3*c^-3*cosintegral(b^-1*(a+b*arccos(c*x)))*sin(a*b^-1)+-1*x*b^-2*c^-2*(a+b*arccos(c*x))^-1+-1*b^-3*c^-3*cos(a*b^-1)*sinintegral(b^-1*(a+b*arccos(c*x)))+-9//8*b^-3*c^-3*cosintegral(3*arccos(c*x)+3*a*b^-1)*sin(3*a*b^-1)+-9//8*b^-3*c^-3*cosintegral(a*b^-1+arccos(c*x))*sin(a*b^-1)+9//8*b^-3*c^-3*cos(a*b^-1)*sinintegral(a*b^-1+arccos(c*x))+9//8*b^-3*c^-3*cos(3*a*b^-1)*sinintegral(3*arccos(c*x)+3*a*b^-1)+(1/2)*b^-1*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-2)
@test integrate(x*(a+b*arccos(c*x))^-3, x) == :(b^-2*x^2*(a+b*arccos(c*x))^-1+-1//2*b^-2*c^-2*(a+b*arccos(c*x))^-1+b^-3*c^-2*cos(2*a*b^-1)*sinintegral(2*arccos(c*x)+2*a*b^-1)+-1*b^-3*c^-2*cosintegral(2*arccos(c*x)+2*a*b^-1)*sin(2*a*b^-1)+(1/2)*x*b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-2)
@test integrate((a+b*arccos(c*x))^-3, x) == :((1/2)*x*b^-2*(a+b*arccos(c*x))^-1+(1/2)*b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-2+(1/2)*b^-3*c^-1*cos(a*b^-1)*sinintegral(b^-1*(a+b*arccos(c*x)))+-1//2*b^-3*c^-1*cosintegral(b^-1*(a+b*arccos(c*x)))*sin(a*b^-1))
@test integrate(x^2*(a+b*arccos(c*x))^(1/2), x) == :(1//3*x^3*(a+b*arccos(c*x))^(1/2)+-1//8*2^(1/2)*pi^(1/2)*b^(1/2)*c^-3*cos(a*b^-1)*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1//8*2^(1/2)*pi^(1/2)*b^(1/2)*c^-3*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+-1//72*6^(1/2)*pi^(1/2)*b^(1/2)*c^-3*cos(3*a*b^-1)*fresnelc(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1//72*6^(1/2)*pi^(1/2)*b^(1/2)*c^-3*fresnels(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(3*a*b^-1))
@test integrate(x*(a+b*arccos(c*x))^(1/2), x) == :((1/2)*x^2*(a+b*arccos(c*x))^(1/2)+-1//4*c^-2*(a+b*arccos(c*x))^(1/2)+-1//8*pi^(1/2)*b^(1/2)*c^-2*cos(2*a*b^-1)*fresnelc(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1//8*pi^(1/2)*b^(1/2)*c^-2*fresnels(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(2*a*b^-1))
@test integrate((a+b*arccos(c*x))^(1/2), x) == :(x*(a+b*arccos(c*x))^(1/2)+-1//2*2^(1/2)*pi^(1/2)*b^(1/2)*c^-1*cos(a*b^-1)*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1//2*2^(1/2)*pi^(1/2)*b^(1/2)*c^-1*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1))
@test integrate(x^2*(a+b*arccos(c*x))^3//2, x) == :(1//3*x^3*(a+b*arccos(c*x))^3//2+-1//3*b*c^-3*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^(1/2)+-1//6*b*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^(1/2)+-3//16*2^(1/2)*pi^(1/2)*b^3//2*c^-3*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+-1//144*6^(1/2)*pi^(1/2)*b^3//2*c^-3*fresnelc(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(3*a*b^-1)+1//144*6^(1/2)*pi^(1/2)*b^3//2*c^-3*cos(3*a*b^-1)*fresnels(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+3//16*2^(1/2)*pi^(1/2)*b^3//2*c^-3*cos(a*b^-1)*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate(x*(a+b*arccos(c*x))^3//2, x) == :((1/2)*x^2*(a+b*arccos(c*x))^3//2+-1//4*c^-2*(a+b*arccos(c*x))^3//2+-3//8*b*x*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^(1/2)+-3//32*pi^(1/2)*b^3//2*c^-2*fresnelc(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(2*a*b^-1)+3//32*pi^(1/2)*b^3//2*c^-2*cos(2*a*b^-1)*fresnels(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate((a+b*arccos(c*x))^3//2, x) == :(x*(a+b*arccos(c*x))^3//2+-3//2*b*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^(1/2)+-3//4*2^(1/2)*pi^(1/2)*b^3//2*c^-1*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+3//4*2^(1/2)*pi^(1/2)*b^3//2*c^-1*cos(a*b^-1)*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate(x^2*(a+b*arccos(c*x))^5//2, x) == :(1//3*x^3*(a+b*arccos(c*x))^5//2+-5//36*b^2*x^3*(a+b*arccos(c*x))^(1/2)+-5//6*x*b^2*c^-2*(a+b*arccos(c*x))^(1/2)+-5//9*b*c^-3*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^3//2+-5//18*b*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^3//2+5//864*6^(1/2)*pi^(1/2)*b^5//2*c^-3*cos(3*a*b^-1)*fresnelc(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+5//864*6^(1/2)*pi^(1/2)*b^5//2*c^-3*fresnels(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(3*a*b^-1)+15//32*2^(1/2)*pi^(1/2)*b^5//2*c^-3*cos(a*b^-1)*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+15//32*2^(1/2)*pi^(1/2)*b^5//2*c^-3*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1))
@test integrate(x*(a+b*arccos(c*x))^5//2, x) == :((1/2)*x^2*(a+b*arccos(c*x))^5//2+-1//4*c^-2*(a+b*arccos(c*x))^5//2+-15//32*b^2*x^2*(a+b*arccos(c*x))^(1/2)+15//64*b^2*c^-2*(a+b*arccos(c*x))^(1/2)+-5//8*b*x*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^3//2+15//128*pi^(1/2)*b^5//2*c^-2*cos(2*a*b^-1)*fresnelc(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+15//128*pi^(1/2)*b^5//2*c^-2*fresnels(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(2*a*b^-1))
@test integrate((a+b*arccos(c*x))^5//2, x) == :(x*(a+b*arccos(c*x))^5//2+-15//4*x*b^2*(a+b*arccos(c*x))^(1/2)+-5//2*b*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^3//2+15//8*2^(1/2)*pi^(1/2)*b^5//2*c^-1*cos(a*b^-1)*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+15//8*2^(1/2)*pi^(1/2)*b^5//2*c^-1*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1))
@test integrate(x^2*(a+b*arccos(c*x))^-1//2, x) == :(-1//4*2^(1/2)*pi^(1/2)*b^-1//2*c^-3*cos(a*b^-1)*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1//12*6^(1/2)*pi^(1/2)*b^-1//2*c^-3*cos(3*a*b^-1)*fresnels(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+1//4*2^(1/2)*pi^(1/2)*b^-1//2*c^-3*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+1//12*6^(1/2)*pi^(1/2)*b^-1//2*c^-3*fresnelc(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(3*a*b^-1))
@test integrate(x*(a+b*arccos(c*x))^-1//2, x) == :((1/2)*pi^(1/2)*b^-1//2*c^-2*fresnelc(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(2*a*b^-1)+-1//2*pi^(1/2)*b^-1//2*c^-2*cos(2*a*b^-1)*fresnels(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate((a+b*arccos(c*x))^-1//2, x) == :(2^(1/2)*pi^(1/2)*b^-1//2*c^-1*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+-1*2^(1/2)*pi^(1/2)*b^-1//2*c^-1*cos(a*b^-1)*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate(x^2*(a+b*arccos(c*x))^-3//2, x) == :(2*b^-1*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-1//2+-1//2*2^(1/2)*pi^(1/2)*b^-3//2*c^-3*cos(a*b^-1)*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1//2*2^(1/2)*pi^(1/2)*b^-3//2*c^-3*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+-1//2*6^(1/2)*pi^(1/2)*b^-3//2*c^-3*cos(3*a*b^-1)*fresnelc(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1//2*6^(1/2)*pi^(1/2)*b^-3//2*c^-3*fresnels(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(3*a*b^-1))
@test integrate(x*(a+b*arccos(c*x))^-3//2, x) == :(-2*pi^(1/2)*b^-3//2*c^-2*cos(2*a*b^-1)*fresnelc(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-2*pi^(1/2)*b^-3//2*c^-2*fresnels(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(2*a*b^-1)+2*x*b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-1//2)
@test integrate((a+b*arccos(c*x))^-3//2, x) == :(2*b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-1//2+-2*2^(1/2)*pi^(1/2)*b^-3//2*c^-1*cos(a*b^-1)*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-2*2^(1/2)*pi^(1/2)*b^-3//2*c^-1*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1))
@test integrate(x^2*(a+b*arccos(c*x))^-5//2, x) == :(4*b^-2*x^3*(a+b*arccos(c*x))^-1//2+-8//3*x*b^-2*c^-2*(a+b*arccos(c*x))^-1//2+2//3*b^-1*c^-1*x^2*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-3//2+6^(1/2)*pi^(1/2)*b^-5//2*c^-3*cos(3*a*b^-1)*fresnels(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))+-1*6^(1/2)*pi^(1/2)*b^-5//2*c^-3*fresnelc(6^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(3*a*b^-1)+-1//3*2^(1/2)*pi^(1/2)*b^-5//2*c^-3*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+1//3*2^(1/2)*pi^(1/2)*b^-5//2*c^-3*cos(a*b^-1)*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate(x*(a+b*arccos(c*x))^-5//2, x) == :(-4//3*b^-2*c^-2*(a+b*arccos(c*x))^-1//2+8//3*b^-2*x^2*(a+b*arccos(c*x))^-1//2+-8//3*pi^(1/2)*b^-5//2*c^-2*fresnelc(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(2*a*b^-1)+2//3*x*b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-3//2+8//3*pi^(1/2)*b^-5//2*c^-2*cos(2*a*b^-1)*fresnels(2*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate((a+b*arccos(c*x))^-5//2, x) == :(4//3*x*b^-2*(a+b*arccos(c*x))^-1//2+2//3*b^-1*c^-1*(1+-1*c^2*x^2)^(1/2)*(a+b*arccos(c*x))^-3//2+-4//3*2^(1/2)*pi^(1/2)*b^-5//2*c^-1*fresnelc(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2))*sin(a*b^-1)+4//3*2^(1/2)*pi^(1/2)*b^-5//2*c^-1*cos(a*b^-1)*fresnels(2^(1/2)*pi^-1//2*b^-1//2*(a+b*arccos(c*x))^(1/2)))
@test integrate((d*x)^5//2*(a+b*arccos(c*x)), x) == :(2//7*d^-1*(d*x)^7//2*(a+b*arccos(c*x))+-4//49*b*c^-1*(d*x)^5//2*(1+-1*c^2*x^2)^(1/2)+20//147*b*c^-7//2*d^5//2*Elliptic.F(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1)+-20//147*b*c^-3*d^2*(d*x)^(1/2)*(1+-1*c^2*x^2)^(1/2))
@test integrate((d*x)^3//2*(a+b*arccos(c*x)), x) == :(2//5*d^-1*(d*x)^5//2*(a+b*arccos(c*x))+-12//25*b*c^-5//2*d^3//2*Elliptic.F(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1)+-4//25*b*c^-1*(d*x)^3//2*(1+-1*c^2*x^2)^(1/2)+12//25*b*c^-5//2*d^3//2*Elliptic.E(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1))
@test integrate((d*x)^(1/2)*(a+b*arccos(c*x)), x) == :(2//3*d^-1*(d*x)^3//2*(a+b*arccos(c*x))+-4//9*b*c^-1*(d*x)^(1/2)*(1+-1*c^2*x^2)^(1/2)+4//9*b*c^-3//2*d^(1/2)*Elliptic.F(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1))
@test integrate((d*x)^-1//2*(a+b*arccos(c*x)), x) == :(2*d^-1*(d*x)^(1/2)*(a+b*arccos(c*x))+-4*b*c^-1//2*d^-1//2*Elliptic.F(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1)+4*b*c^-1//2*d^-1//2*Elliptic.E(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1))
@test integrate((d*x)^-3//2*(a+b*arccos(c*x)), x) == :(d^-1*(d*x)^-1//2*(-2a+-2*b*arccos(c*x))+-4*b*c^(1/2)*d^-3//2*Elliptic.F(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1))
@test integrate((d*x)^-5//2*(a+b*arccos(c*x)), x) == :(-1//3*d^-1*(d*x)^-3//2*(2a+2*b*arccos(c*x))+-4//3*b*c^3//2*d^-5//2*Elliptic.F(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1)+4//3*b*c^3//2*d^-5//2*Elliptic.E(arcsin(c^(1/2)*d^-1//2*(d*x)^(1/2)),-1)+4//3*b*c*d^-2*(d*x)^-1//2*(1+-1*c^2*x^2)^(1/2))

max_stars_repo_path: test/5 Inverse trig functions/5.2 Inverse cosine/5.2.2 (d x)^m (a+b arccos(c x))^n.jl


text:
module FastAI

using Base: NamedTuple
using Colors: colormaps_sequential
using Reexport
@reexport using DLPipelines
@reexport using FluxTraining
@reexport using DataLoaders
@reexport using Flux

using Animations
using Makie
using Colors
using DataAugmentation
using DataAugmentation: getbounds, Bounds
import DLPipelines: methoddataset, methodmodel, methodlossfn, methoddataloaders,
    mockmodel, mocksample, predict, predictbatch, mockmodel, encode, encodeinput,
    encodetarget, decodeŷ, decodey
using IndirectArrays: IndirectArray
using LearnBase: getobs, nobs
using FilePathsBase
using FixedPointNumbers
using Flux
using Flux.Optimise
import Flux.Optimise: apply!, Optimiser, WeightDecay
using FluxTraining: Learner, handle
using FluxTraining.Events
using JLD2: jldsave, jldopen
using Markdown
using MLDataPattern
using Parameters
using PrettyTables
using StaticArrays
using Setfield
using ShowCases
using Tables
using Statistics
using Test: @testset, @test, @test_nowarn

include("plotting.jl")
include("learner.jl")

# Data block API
include("datablock/block.jl")
include("datablock/encoding.jl")
include("datablock/method.jl")
include("datablock/describe.jl")
include("datablock/checks.jl")
include("datablock/wrappers.jl")

# submodules
include("datasets/Datasets.jl")
@reexport using .Datasets

include("models/Models.jl")
using .Models

# Encodings
include("encodings/tabularpreprocessing.jl")
include("encodings/onehot.jl")
include("encodings/imagepreprocessing.jl")
include("encodings/projective.jl")
include("encodings/keypointpreprocessing.jl")

# Training interface
include("datablock/models.jl")
include("datablock/loss.jl")
include("datablock/plot.jl")

# training
include("training/paramgroups.jl")
include("training/discriminativelrs.jl")
include("training/utils.jl")
include("training/onecycle.jl")
include("training/finetune.jl")
include("training/lrfind.jl")
include("training/metrics.jl")

include("serialization.jl")

include("fasterai/methodregistry.jl")
include("fasterai/learningmethods.jl")
include("fasterai/defaults.jl")

export
    # submodules
    Datasets,
    Models,
    datasetpath,
    loadtaskdata,
    mapobs,
    groupobs,
    filterobs,
    shuffleobs,
    datasubset,

    # method API
    methodmodel,
    methoddataset,
    methoddataloaders,
    methodlossfn,
    getobs,
    nobs,
    predict,
    predictbatch,

    # plotting API
    plotbatch,
    plotsample,
    plotsamples,
    plotpredictions,
    makebatch,

    # blocks
    Image,
    Mask,
    Label,
    LabelMulti,
    Keypoints,
    Many,
    TableRow,
    Continuous,

    # encodings
    encode,
    decode,
    setup,
    ProjectiveTransforms,
    ImagePreprocessing,
    OneHot,
    KeypointPreprocessing,
    Only,
    Named,
    augs_projection, augs_lighting,
    TabularPreprocessing,

    BlockMethod,
    describemethod,
    checkblock,

    # learning methods
    findlearningmethods,
    ImageClassificationSingle,
    ImageClassificationMulti,
    ImageSegmentation,
    ImageKeypointRegression,
    TabularClassificationSingle,
    TabularRegression,

    # training
    methodlearner,
    Learner,
    fit!,
    fitonecycle!,
    finetune!,
    lrfind,
    savemethodmodel,
    loadmethodmodel,
    accuracy_thresh,

    gpu,
    plot

end  # module

ext: jl


text:
"""
    mutable struct PropertyStore
        gmaps::Dict{String, Any}
        emaps::Dict{String,AEdgeMap}
        vmaps::Dict{String,AVertexMap}
    end

A type storing properties associated to networks.
"""
mutable struct PropertyStore
    gmaps::Dict{String, Any}
    emaps::Dict{String,AEdgeMap}
    vmaps::Dict{String,AVertexMap}
end

PropertyStore() = PropertyStore(Dict{String, Any}(), Dict{String,AEdgeMap}(), Dict{String,AVertexMap}())

### GRAPH
function set_graph_property!(p::PropertyStore, name::AbstractString, x)
    p.gmaps[name] = x
end

function rem_graph_property!(p::PropertyStore, name::AbstractString)
    !haskey(p.gmaps, name) && error("Property $name not present.")
    delete!(p.gmaps, name)
    p.gmaps
end

graph_property(p::PropertyStore, name::AbstractString) = p.gmaps[name]
graph_property(p::PropertyStore) = p.gmaps

### EDGE
function add_edge_property!(p::PropertyStore, name::AbstractString, emap::AEdgeMap)
    haskey(p.emaps, name) && error("Property $name already present.")
    p.emaps[name] = emap
end

function rem_edge_property!(p::PropertyStore, name::AbstractString)
    !haskey(p.emaps, name) && error("Property $name not present.")
    delete!(p.emaps, name)
    p.emaps
end

edge_property(p::PropertyStore, name::AbstractString) = p.emaps[name]

edge_property(p::PropertyStore) = p.emaps

function edge_property(p::PropertyStore, e::AEdge)
    pe = Dict{String, Any}()
    for (name, em) in edge_property(p)
        if haskey(em, e)
            pe[name] = em[e]
        end
    end
    return pe
end

## VERTEX
function add_vertex_property!(p::PropertyStore, name::AbstractString, vmap::AVertexMap)
    haskey(p.vmaps, name) && error("Property $name already present.")
    p.vmaps[name] = vmap
end
function rem_vertex_property!(p::PropertyStore, name::AbstractString)
    !haskey(p.vmaps, name) && error("Property $name not present.")
    delete!(p.vmaps, name)
    p.vmaps
end

vertex_property(p::PropertyStore, name::AbstractString) = p.vmaps[name]
vertex_property(p::PropertyStore) = p.vmaps

function vertex_property(p::PropertyStore, i::Integer)
    pi = Dict{String, Any}()
    for (name, vm) in vertex_property(p)
        if haskey(vm, i)
            pi[name] = vm[i]
        end
    end
    return pi
end

###

function swap_vertices!(props::PropertyStore, u::Integer, v::Integer)
    for vmap in values(props.vmaps)
        hasu, hasv = haskey(vmap, u), haskey(vmap, v)
        if hasu && hasv
            vmap[u], vmap[v] = vmap[v], vmap[u]
        else
            if hasu
                vmap[v] = vmap[u]
            end
            if hasv
                vmap[u] = vmap[v]
            end
        end
    end
    #TODO should swap edges for non indexed graphs
end

has_graph_property(props::PropertyStore, name::AbstractString) = haskey(props.gmaps, name)
has_vertex_property(props::PropertyStore, name::AbstractString) = haskey(props.vmaps, name)
has_vertex_property(props::PropertyStore, name::AbstractString, v::Integer) =
            haskey(props.vmaps, name) && haskey(props.vmaps[name], v)

has_edge_property(props::PropertyStore, name::AbstractString) = haskey(props.emaps, name)
has_edge_property(props::PropertyStore, name::AbstractString, e::AEdge) =
            haskey(props.emaps, name) && haskey(props.emaps[name], e)

function ==(p1::PropertyStore, p2::PropertyStore)
    oke = true
    for name in keys(p1.emaps)
        if !haskey(p2.emaps, name)
            oke = false
            break
        end
        oke &= p1.emaps[name] == p2.emaps[name]
    end
    return (p1.gmaps == p2.gmaps) && (p1.emaps == p2.emaps) && oke
end


text:
md"""
# Advanced Tutorial

The previous examples have been focused on giving an overview of the purpose of this library using the simple high level API. For more advanced usage, there are
many of ways to customize the behavior and in other ways alter or hook in to the functionality. Some of the more important concepts are described below.

To make it more convenient to extend `NaiveNASlib`, the two submodules `NaiveNASlib.Advanced` and `NaiveNASlib.Extend` export most of the useful stuff, such as
abstract types and composable strategies. For now they are also part of the public API, but in future releases they might be moved to separate subpackages so
that they can be versioned separately (e.g NaiveNASlibCore).
"""
using NaiveNASlib.Advanced, NaiveNASlib.Extend

md"""
## Strategies

For more or less all operations which mutate the graph, it is possible achieve fine grained control of the operation through selecting a strategy.

Here is an example of strategies for changing the size.
"""

@testset "Strategies" begin #src
# First we make a simple graph where one vertex has a constraint for changing the size.
invertex = inputvertex("in", 3)
layer1 = linearvertex(invertex, 4)
joined = conc(scalarmult(layer1, 2), scalarmult(layer1, 3), dims=1)
# Note that `joined` can only change in steps of 2.
#
# Strategy to try to change it by one and throw an error when not successful.
exact_or_fail = ΔNoutExact(joined => 1; fallback=ThrowΔSizeFailError("Size change failed!!"))

# Note that we now call [`Δsize!`](@ref) instead of [`Δnout!`](@ref) as the wanted action is given by the strategy.
@test_throws NaiveNASlib.ΔSizeFailError Δsize!(exact_or_fail)
@test nout(joined) == 2*nout(layer1) == 8 # No change was made.

# Try to change by one and fail silently when not successful.
exact_or_noop = ΔNoutExact(joined=>1;fallback=ΔSizeFailNoOp())

@test !Δsize!(exact_or_noop)
@test nout(joined) == 2*nout(layer1) == 8 # No change was made.

# In many cases it is ok to not get the exact change which was requested.
relaxed_or_fail = ΔNoutRelaxed(joined=>1;fallback=ThrowΔSizeFailError("This should not happen!!"))

@test Δsize!(relaxed_or_fail)
## Changed by two as this was the smallest possible change.
@test nout(joined) == 2*nout(layer1) == 10

# Logging when fallback is applied is also possible.
using Logging: Info
exact_or_log_then_relax = ΔNoutExact(joined=>1;
                                        fallback=LogΔSizeExec(
                                                        "Exact failed, relaxing",
                                                        Info,
                                                        relaxed_or_fail))

@test_logs (:info, "Exact failed, relaxing") Δsize!(exact_or_log_then_relax)
@test nout(joined) == 2*nout(layer1) == 12

# If one wants to see every size change we can set up an [`AfterΔSizeCallback`](@ref) strategy to log it for us like this:
exact_or_log_then_relax_verbose = logafterΔsize(v -> "some vertex";base=exact_or_log_then_relax)

@test_logs(
    (:info, "Exact failed, relaxing"),
    (:info, r"Change nin of some vertex"),
    (:info, r"Change nout of some vertex"),
    match_mode=:any,
    Δsize!(exact_or_log_then_relax_verbose))
end #src

md"""
A similar pattern is used for most other mutating operations. See the advanced reference documentation for the complete set.
"""
md"""
## Traits

A variant (bastardization?) of the [holy trait](https://docs.julialang.org/en/v1/manual/methods/#Trait-based-dispatch-1) pattern is used to
annotate the type of a vertex. The core idea is discussed a bit in the [Terminology](@ref) section, but it is
also possible to attach other information and behaviors by freeriding on this mechanism.

This is done by adding the argument `traitdecoration` when creating a vertex and supplying a function which takes a trait and return a new trait
(which typically wraps the input).
"""

@testset "Traits" begin #src
# Naming vertices is so useful for logging and debugging I almost made it mandatory.
#
# If a vertex does not have the named trait, [`name`](@ref) will return a generic string. Compare
noname = linearvertex(inputvertex("in", 2), 2)
@test name(noname) == "MutationVertex::SizeAbsorb"
# with
hasname = absorbvertex(LinearLayer(2, 3), inputvertex("in", 2), traitdecoration = t -> NamedTrait("named layer", t))
@test name(hasname) == "named layer"

# which is basically what the convenience methods do under the hood.
hasnametoo = absorbvertex("also named", LinearLayer(2, 3), inputvertex("in", 2))
@test name(hasnametoo) == "also named"

# [`AfterΔSizeTrait`](@ref) can be used to attach an [`AbstractAfterΔSizeStrategy`](@ref) to an individual vertex.
# In this case we use [`logafterΔsize`](@ref) from the example above.
verbose_vertex_info(v) = string(name(v),
                                 " with inputs=[", join(name.(inputs(v)),  ", "),
                                "] and outputs=[", join(name.(outputs(v)), ", "), ']')
named_verbose_logging(t) = AfterΔSizeTrait(
                                        logafterΔsize(verbose_vertex_info),
                                        NamedTrait("layer1", t))
layer1 = absorbvertex(  LinearLayer(2, 3),
                        inputvertex("in", 2),
                        traitdecoration = named_verbose_logging)
# The above is a mouthful, but `NaiveNASlib.Advanced` exports the [`named`](@ref) and [`logged`](@ref) functions for convenience
layer2 = absorbvertex(  LinearLayer(nout(layer1), 4),
                        layer1;
                        traitdecoration = logged(name) ∘ named("layer2"))

# Now logs for `layer2` are less verbose than logs for `layer1` due to `name` being used to print the
# vertex instead of `verbose_vertex_info`.
@test_logs(
(:info, "Change nout of layer1 with inputs=[in] and outputs=[layer2] by [1, 2, 3, -1]"),
(:info, "Change nin of layer2 by [1, 2, 3, -1]"),
Δnout!(layer1, 1))

# For more elaborate traits with elementwise operations one can use traitconf and `>>`.
add = traitconf(logged(verbose_vertex_info) ∘ named("layer1+layer2")) >> layer1 + layer2
@test name(add) == "layer1+layer2"

@test_logs(
(:info, "Change nout of layer1 with inputs=[in] and outputs=[layer2, layer1+layer2] by [1, 2, 3, 4, -1]"),
(:info, "Change nin of layer2 by [1, 2, 3, 4, -1]"),
(:info, "Change nout of layer2 by [1, 2, 3, 4, -1]"),
(:info, "Change nin of layer1+layer2 with inputs=[layer1, layer2] and outputs=[] by [1, 2, 3, 4, -1] and [1, 2, 3, 4, -1]"),
(:info, "Change nout of layer1+layer2 with inputs=[layer1, layer2] and outputs=[] by [1, 2, 3, 4, -1]"),
Δnout!(layer1, 1))

# When creating own trait wrappers, remember to subtype [`DecoratingTrait`](@ref) or else there will be pain.

struct PainfulTrait{T<:MutationTrait} <: MutationTrait
    base::T
end
painlayer = absorbvertex(   LinearLayer(2, 3),
                            inputvertex("in", 2);
                            traitdecoration = PainfulTrait)

# Wrong! Not a subtype of [`DecoratingTrait`](@ref).
@test_throws MethodError Δnout!(painlayer, 1)

# Now one must implement a lot of methods for `PainfulTrait`...
#
# Lets try that again:
struct SmoothSailingTrait{T<:MutationTrait} <: DecoratingTrait
    base::T
end
# Right! Is a subtype of [`DecoratingTrait`](@ref).
# Just implement [`base`](@ref base(t::DecoratingTrait)) and all will be fine.
NaiveNASlib.base(t::SmoothSailingTrait) = t.base

smoothlayer = absorbvertex( LinearLayer(2, 3),
                            inputvertex("in", 2);
                            traitdecoration = SmoothSailingTrait)

@test Δnout!(smoothlayer, 1)
@test nout(smoothlayer) == 4
end #src

md"""
## Graph instrumentation and modification

In many cases it is desirable to change things like traits of an existing graph. This can be achieved through [Functors.jl](https://github.com/FluxML/Functors.jl), often through clever usage of the `walk` function.

Depending on what one wants to achieve, it can be more or less messy. Here is a pretty messy example:
"""

@testset "Graph instrumentation and modification" begin #src
invertex = inputvertex("in", 2)
layer1 = linearvertex(invertex, 3)
layer2 = linearvertex(layer1, 4)

graph = CompGraph(invertex, layer2)

@test name.(vertices(graph)) == ["in", "MutationVertex::SizeAbsorb", "MutationVertex::SizeAbsorb"]

# Ok, lets add names to `layer1` and `layer2` and change the name of `invertex` by using `Functors`.
import Functors

# Add a name to `layer1` and `layer2`
function walk(f, v::MutationVertex)
    ## This is probably not practical to do in a real graph, so make sure you have names when first creating it...
    name = v == layer1 ? "layer1" : "layer2"
    addname(x) = x
    ## SizeAbsorb has no fields, otherwise we would have had to use walk to wrap it
    addname(t::SizeAbsorb) = NamedTrait(name, t)
    Functors._default_walk(v) do x
        Functors.fmap(addname, x; walk)
    end
end

# Change name of `invertex` once we get there.
# We could also just have made a string version of `addname` above
# since there are no other strings in the graph, but this is safer.
function walk(f, v::InputVertex)
    rename(x) = x
    rename(s::String) = "in changed"
    Functors._default_walk(v) do x
        Functors.fmap(rename, x; walk)
    end
end

# Everything else just gets functored as normal.
walk(f, x) = Functors._default_walk(f, x)

# I must admit that thinking about what this does makes me a bit dizzy...
namedgraph = Functors.fmap(identity, graph; walk)

@test name.(vertices(namedgraph)) == ["in changed", "layer1", "layer2"]
end #src


text:
# Protein Data Bank (PDB) format
# ==============================

import Base: readline, start, done, next

type State
    io::IO
    record::Symbol
    linen::Int
    # the number of ATOM/HETATM records
    n_atom_records::Int
    max_atom_serial::Int
    line::ASCIIString
    State(io::IO) = new(io, :Init, 0, 0, 0)
end

macro pdb_assert(ex, msg)
    quote
        @assert ex error("$($msg) (line:$($s.linen))")
    end
end

function readline(s::State)
    s.line = readline(s.io)
    s.linen += 1
    s.record = (
        beginswith(s.line, "HEADER") ? :HEADER :
        beginswith(s.line, "OBSLTE") ? :OBSLTE :
        beginswith(s.line, "TITLE ") ? :TITLE  :
        beginswith(s.line, "SPLIT ") ? :SPLIT  :
        beginswith(s.line, "CAVEAT") ? :CAVEAT :
        beginswith(s.line, "COMPND") ? :COMPND :
        beginswith(s.line, "SOURCE") ? :SOURCE :
        beginswith(s.line, "KEYWDS") ? :KEYWDS :
        beginswith(s.line, "EXPDTA") ? :EXPDTA :
        beginswith(s.line, "NUMMDL") ? :NUMMDL :
        beginswith(s.line, "MDLTYP") ? :MDLTYP :
        beginswith(s.line, "AUTHOR") ? :AUTHOR :
        beginswith(s.line, "REVDAT") ? :REVDAT :
        beginswith(s.line, "SPRSDE") ? :SPRSDE :
        beginswith(s.line, "JRNL  ") ? :JRNL   :
        beginswith(s.line, "REMARK") ? :REMARK :
        beginswith(s.line, "DBREF ") ? :DBREF  :
        beginswith(s.line, "DBREF1") ? :DBREF1 :
        beginswith(s.line, "DBREF2") ? :DBREF2 :
        beginswith(s.line, "SEQADV") ? :SEQADV :
        beginswith(s.line, "SEQRES") ? :SEQRES :
        beginswith(s.line, "MODRES") ? :MODRES :
        beginswith(s.line, "HET   ") ? :HET    :
        beginswith(s.line, "HETNAM") ? :HETNAM :
        beginswith(s.line, "HETSYN") ? :HETSYN :
        beginswith(s.line, "FORMUL") ? :FORMUL :
        beginswith(s.line, "HELIX ") ? :HELIX  :
        beginswith(s.line, "SHEET ") ? :SHEET  :
        beginswith(s.line, "SSBOND") ? :SSBOND :
        beginswith(s.line, "LINK  ") ? :LINK   :
        beginswith(s.line, "CISPEP") ? :CISPEP :
        beginswith(s.line, "SITE  ") ? :SITE   :
        beginswith(s.line, "CRYST1") ? :CRYST1 :
        beginswith(s.line, "ORIGX1") ? :ORIGX1 :
        beginswith(s.line, "ORIGX2") ? :ORIGX2 :
        beginswith(s.line, "ORIGX3") ? :ORIGX3 :
        beginswith(s.line, "SCALE1") ? :SCALE1 :
        beginswith(s.line, "SCALE2") ? :SCALE2 :
        beginswith(s.line, "SCALE3") ? :SCALE3 :
        beginswith(s.line, "MTRIX1") ? :MTRIX1 :
        beginswith(s.line, "MTRIX2") ? :MTRIX2 :
        beginswith(s.line, "MTRIX3") ? :MTRIX3 :
        beginswith(s.line, "MODEL ") ? :MODEL  :
        beginswith(s.line, "ATOM  ") ? :ATOM   :
        beginswith(s.line, "ANISOU") ? :ANISOU :
        beginswith(s.line, "TER   ") ? :TER    :
        beginswith(s.line, "HETATM") ? :HETATM :
        beginswith(s.line, "ENDMDL") ? :ENDMDL :
        beginswith(s.line, "CONECT") ? :CONECT :
        beginswith(s.line, "MASTER") ? :MASTER :
        beginswith(s.line, "END   ") ? :END    : :Unknown
    )
    s.line
end

type Title
    # HEADER
    classification::ASCIIString
    dep_date::ASCIIString
    idcode::ASCIIString
    # TITLE
    title::ASCIIString
    # OBSLTE
    rep_date::ASCIIString
    r_idcodes::Vector{ASCIIString}
    function Title(classification, dep_date, idcode, titles)
        new(classification, dep_date, idcode, titles)
    end
    function Title(classification, dep_date, idcode, titles, rep_date, r_idcodes)
        new(classification, dep_date, idcode, titles, rep_date, r_idcodes)
    end
end

isobsolete(t::Title) = isdefined(t, :rep_date)

type Remark
    remarks::Dict{Int,Vector{ASCIIString}}
end

type DBReference
    idcode::ASCIIString
    chainID::Char
    seqrange::UnitRange{Int}
    ins::(Char,Char)
    database::ASCIIString
    dbaccession::ASCIIString
    dbidcode::ASCIIString
    dbseqrange::UnitRange{Int}
    dbins::(Char,Char)
end

type PrimaryStructure
    dbref::Nullable{DBReference}
    seqres::Dict{Char,Vector{ASCIIString}}
end

function show(io::IO, primstr::PrimaryStructure)
    chains = sort(collect(keys(primstr.seqres)))
    n_chains = length(chains)
    print(io, "Primary Structure ($n_chains chains)")
    for chain in chains
        println(io)
        residues = primstr.seqres[chain]
        n_residues = length(residues)
        if n_residues > 4
            seq = join(vcat(residues[1:2], "...", residues[end-1:end]), "-")
        else
            seq = join(residues, "-")
        end
        print(io, " * $chain ($n_residues residues): $seq")
    end
end

type Het
    hetID::ASCIIString
    chainID::Char
    seqnum::Int
    icode::Char
    numhetatoms::Int
    text::ASCIIString
end

type Hetnam
    hetID::ASCIIString
    text::ASCIIString
end

type Hetsyn
    hetID::ASCIIString
    hetsyn::ASCIIString
end

type Formula
    compnum::Int
    hetID::ASCIIString
    asterisk::Bool
    text::ASCIIString
end

type Heterogen
    hets::Vector{Het}
    hetnams::Vector{Hetnam}
    hetsyns::Vector{Hetsyn}
    formulae::Vector{Formula}
end

function show(io::IO, hetgen::Heterogen)
    print(io, "Heterogeneous Components ($(length(hetgen.hets)) molecules/residues)")
    for het in hetgen.hets
        println(io)
        i = findfirst(x -> x.hetID == het.hetID, hetgen.hetnams)
        hettxt = i > 0 ? hetgen.hetnams[i].text : ""
        print(io, " * $(het.hetID) ($(het.numhetatoms) atoms): $(hettxt)")
    end
end

type Helix
   sernum::Int
   helixID::ASCIIString
   initresname::ASCIIString
   initchainID::Char
   initseqnum::Int
   initicode::Char
   endresname::ASCIIString
   endchainid::Char
end

type Atom
    hetero::Bool
    serial::Int
    name::ASCIIString
    altloc::Char
    resname::ASCIIString
    chainID::Char
    resseq::Int
    icode::Char
    x::Float32
    y::Float32
    z::Float32
    occupancy::Float32
    tempfactor::Float32
    element::Element
    charge::Float32
    anisou::Matrix{Int}  # optional

    # constructors
    Atom(hetero, serial, name, altloc, resname, chainID, resseq, icode, x, y, z, occupancy, tempfactor, element, charge) =
        new(hetero, serial, name, altloc, resname, chainID, resseq, icode, x, y, z, occupancy, tempfactor, element, charge)
    Atom(hetero, serial, name, altloc, resname, chainID, resseq, icode, x, y, z, occupancy, tempfactor, element, charge, anisou) =
        new(hetero, serial, name, altloc, resname, chainID, resseq, icode, x, y, z, occupancy, tempfactor, element, charge, anisou)
end

ishetero(a::Atom) = a.hetero

type Molecule
    name::ASCIIString
    seq::Int
    range::UnitRange{Int}
end

typealias Residue Molecule

type Chain
    chainID::Char
    residues::Vector{Residue}
    range::UnitRange{Int}
end

# Model - the core object to store coordinate data
# if a PDB entry has only one model, all ATOMs and HETATMs are stored in this object;
# otherwise a Model object stores coordinates between MODEL and ENDMDL records
type Model
    # model serial found in the MODEL record
    # if an entry has only one model and there is no MODEL record, this serial is set to 0
    serial::Int

    # chains terminated by a TER record
    chains::Vector{Chain}

    # molecules which aren't a monomer within a chain (e.g. metal ions, ligand, water, etc.)
    molecules::Vector{Molecule}

    # atoms found in ATOM/HETATM records
    # if there are multiple conformations, only the first occurring locations are
    # stored in this field: alternative locations are stored in altloc field (see below)
    atoms::Vector{Atom}

    # alternatives locations (key: serial of atoms, value: atoms found in alternative locations)
    altlocs::Dict{Int,Vector{Atom}}
end

function Model(serial::Int, chains::Vector{(Char,Vector{Atom})}, molecules::Vector{Atom})
    chains′ = Chain[]
    atoms′ = Atom[]
    altlocs′ = Dict{Int,Vector{Atom}}()
    for (chainID, atoms) in chains
        chain_start = endof(atoms′) + 1
        residues = read_mols!(atoms′, altlocs′, atoms, Residue)
        chain_stop = endof(atoms′)
        push!(chains′, Chain(chainID, residues, chain_start:chain_stop))
    end
    molecules′ = read_mols!(atoms′, altlocs′, molecules, Molecule)
    Model(serial, chains′, molecules′, atoms′, altlocs′)
end

function read_mols!{T<:Union(Molecule,Residue)}(dst, altlocs, src, ::Type{T})
    first_altloc_i = findfirst(a -> a.altloc != ' ', src)
    first_altloc = first_altloc_i > 0 ? src[first_altloc_i].altloc : '\0'
    mol_start = endof(dst) + 1
    mol_stop = endof(dst)
    mols = T[]
    for (i, atom) in enumerate(src)
        if atom.altloc == ' ' || atom.altloc == first_altloc
            if i > 1 && atom.resseq != src[i-1].resseq
                mol_stop = endof(dst)
                push!(mols, T(src[i-1].resname, src[i-1].resseq, mol_start:mol_stop))
                mol_start = mol_stop + 1
            end
            push!(dst, atom)
        else
            # this atom is in an alternative location
            if !haskey(altlocs, atom.serial)
                altlocs[atom.serial] = Atom[]
            end
            push!(altlocs[atom.serial], atom)
        end
    end
    if mol_stop < endof(dst)
        mol_stop = endof(dst)
        push!(mols, T(src[end].resname, src[end].resseq, mol_start:mol_stop))
    end
    mols
end

name(a::Atom) = a.name
name(r::Residue) = r.name
name(c::Chain) = c.chainID

type Coordinate
    models::Vector{Model}
end

type Connectivity
    connection::SparseMatrixCSC{Bool,Int}
end

isconnected(a::Atom, b::Atom) = connection[a.serial,b.serial] || connection[b.serial,a.serial]

type Entry
    title::Title
    remark::Remark
    primarystructure::PrimaryStructure
    heterogen::Heterogen
    coordinate::Coordinate
    connectivity::Connectivity
end

name(entry::Entry) = entry.title.idcode
models(entry::Entry) = entry.coordinate.models

function show(io::IO, ent::Entry)
    println(io, "Protein: $(ent.title.classification) [$(ent.title.idcode)]")
    models = ent.coordinate.models
    print(io, "  $(length(models)) models")
    if length(models) == 1
        println(io)
        model = first(models)
        n_chains = length(model.chains)
        print(io, "  $n_chains chains - ")
        println(join(map(c -> c.chainID, model.chains), " "))
        print(io, "  $(length(model.atoms)) atoms")
    end
end

# Iterators
# ---------

type AtomIterator{T}
    model::Model
    view::T
end

start(::AtomIterator) = 1

eachatom(model::Model, atom::Atom) = AtomIterator(model, atom)
done(iter::AtomIterator{Atom}, i) = i > 1
next(iter::AtomIterator{Atom}, i) = iter.view, i + 1

eachatom{T<:Union(Molecule,Chain)}(model::Model, view::T) = AtomIterator(model, view)
done{T<:Union(Molecule,Chain)}(iter::AtomIterator{T}, i) = i > endof(iter.view.range)
next{T<:Union(Molecule,Chain)}(iter::AtomIterator{T}, i) = iter.model.atoms[iter.view.range[i]], i + 1

eachatom(model::Model) = AtomIterator(model, model)
done(iter::AtomIterator{Model}, i) = i > endof(iter.model.atoms)
next(iter::AtomIterator{Model}, i) = iter.model.atoms[i], i + 1

type MoleculeIterator{T}
    model::Model
    view::T
end

eachmol(model::Model, mol::Molecule) = MoleculeIterator(model, mol)
start(::MoleculeIterator{Molecule}) = 1
done(iter::MoleculeIterator{Molecule}, i) = i > 1
next(iter::MoleculeIterator{Molecule}, i) = iter.view, i + 1

eachmol(model::Model, chain::Chain) = MoleculeIterator(model, chain)
start(::MoleculeIterator{Chain}) = 1
done(iter::MoleculeIterator{Chain}, i) = i > endof(iter.view.residues)
next(iter::MoleculeIterator{Chain}, i) = iter.view.residue[i], i + 1

eachmol(model::Model) = MoleculeIterator(model, model)
start(::MoleculeIterator{Model}) = (1, 1)
done(iter::MoleculeIterator{Model}, ij) = ij[1] > endof(iter.view.chains) && ij[2] > endof(iter.view.molecules)
function next(iter::MoleculeIterator{Model}, ij)
    i, j = ij
    if i <= endof(iter.view.chains)
        # first, iterate through molecules (residues) within chains
        chain = iter.view.chains[i]
        mol = chain.residues[j]
        j += 1
        if j > endof(chain.residues)
            i += 1
            j = 1
        end
    else
        # second, iterate through molecules out of chains
        mol = iter.view.molecules[j]
        j += 1
    end
    return mol, (i, j)
end

type ChainIterator{T}
    model::Model
    view::T
end

start(::ChainIterator) = 1

eachchain(model::Model, chain::Chain) = ChainIterator(model, chain)
done(iter::ChainIterator{Chain}, i) = i > 1
next(iter::ChainIterator{Chain}, i) = iter.view, i + 1

eachchain(model::Model) = ChainIterator(model, model)
done(iter::ChainIterator{Model}, i) = i > endof(iter.model.chains)
next(iter::ChainIterator{Model}, i) = iter.model.chains[i], i + 1

# Parsers
# -------
#
# Every parser takes a current state and returns parsed object.
# Since a PDB entry is divided into several sections in its format specs, an entry parser
# is also separated into sub-parsers corresponding to each section:
#   * parse_title_section
#   * parse_remark_section
#   * parse_primary_structure_section
#   * parse_coodinate_section
# And each section parser calls record parsers contained in its section.
# The specs of the PDB file format can be found in:
#     "Atomic Coordinate Entry Format Description"  Version 3.3: July, 2011
#     http://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html
# And other official documents about annotations:
#     "wwPDB Processing Procedures and Policies Document"
#     Section A: http://www.wwpdb.org/documentation/procedure  January 2014 Version 2.7
#     Section B: http://www.wwpdb.org/documentation/policy  April 2015 version 3.6
#
# PDB-101 would also help you understand the format:
#     http://www.rcsb.org/pdb/101/structural_view_of_biology.do

function parse_entry(io::IO)
    s = State(io)
    readline(s)
    title = parse_title_section(s)
    # TODO
    while s.record !== :REMARK; readline(s); end
    remark = parse_remark_section(s)
    primarystructure = parse_primary_structure_section(s)
    heterogen = parse_heterogen_section(s)
    # TODO
    while s.record !== :ATOM && s.record !== :MODEL; readline(s); end
    coordinate = parse_coodinate_section(s)
    connectivity = parse_connectivity_section(s)
    parse_bookkeeping_section(s)
    @assert eof(s.io)
    Entry(title, remark, primarystructure, heterogen, coordinate, connectivity)
end

# Title Section
# -------------

function parse_title_section(s)
    @assert s.record === :HEADER "the first record should be HEADER"
    classification, dep_date, idcode = parse_header(s)
    readline(s)
    obsolete = false
    if s.record === :OBSLTE
        rep_date, idcode′, r_idcodes = parse_obslte(s)
        @assert idcode′ == idcode "obsoleted idCode should match the header idCode (header: $idcode, obsoleted: $idcode′)"
        obsolete = true
    end
    titles = parse_title(s)
    obsolete ? Title(classification, dep_date, idcode, titles, rep_date, r_idcodes) :
               Title(classification, dep_date, idcode, titles)
end

function parse_header(s)
    # classification, depDate, idCode
    strip(s.line[11:50]), s.line[51:59], s.line[63:66]
end

function parse_obslte(s)
    rep_date = s.line[12:20]
    id_code = s.line[22:25]
    r_id_codes = ASCIIString[]
    while s.record === :OBSLTE
        for i in 32:5:72
            r_id_codes = s.line[i:i+3]
            r_id_codes == "    " && break
            push(r_id_codes, r_id_codes)
        end
        readline(s)
    end
    rep_date, id_code, r_id_codes
end

function parse_title(s)
    titles = ASCIIString[]
    while s.record === :TITLE
        push!(titles, rstrip(s.line[11:80]))
        readline(s)
    end
    join(titles)
end

# Remark Section
# --------------

parse_remark_section(s) = parse_remark(s)

function parse_remark(s)
    remarks = Dict{Int,Vector{ASCIIString}}()
    while s.record === :REMARK
        remarknum = int(s.line[8:10])
        if !haskey(remarks, remarknum)
            remarks[remarknum] = ASCIIString[]
        end
        push!(remarks[remarknum], s.line[12:79])
        readline(s)
    end
    Remark(remarks)
end

# Primary Structure Section
# -------------------------

function parse_primary_structure_section(s)
    if s.record === :DBREF
        while s.record === :DBREF
            dbref = Nullable(parse_dbref(s))
            readline(s)
        end
    else
        dbref = Nullable{DBReference}()
    end
    while s.record === :DBREF1 || s.record === :DBREF2
        # TODO
        readline(s)
    end
    while s.record === :SEQADV
        parse_seqadv(s)
        readline(s)
    end
    seqres = parse_seqres(s)
    while s.record === :MODRES
        parse_modres(s)
        readline(s)
    end
    PrimaryStructure(dbref, seqres)
end

function parse_dbref(s)
    idcode = s.line[8:11]
    chainID = s.line[13]
    seqbegin = int(s.line[15:18])
    insertbegin = s.line[19]
    seqend = int(s.line[21:24])
    insertend = s.line[25]
    database = strip(s.line[27:32])
    dbaccession = strip(s.line[34:41])
    dbidcode = strip(s.line[43:54])
    dbseqbegin = int(s.line[56:60])
    idbnsbeg = s.line[61]
    dbseqend = int(s.line[63:67])
    dbinsend = s.line[68]
    DBReference(
        idcode, chainID, seqbegin:seqend, (insertbegin,insertend),
        database, dbaccession, dbidcode, dbseqbegin:dbseqend, (idbnsbeg,dbinsend)
    )
end

function parse_seqadv(s)
    idcode = s.line[8:11]
    resname = strip(s.line[13:15])
    chainID = s.line[17]
    seqnum = int(s.line[19:22])
    icode = s.line[23]
    database = strip(s.line[25:28])
    dbaccession = strip(s.line[30:38])
    dbres = strip(s.line[40:42])
    dbseq = int(s.line[44:48])
    conflict = strip(s.line[50:70])
end

function parse_seqres(s)
    seqres = Dict{Char,Vector{ASCIIString}}()
    while s.record === :SEQRES
        current_chain = s.line[12]
        chain = s.line[12]
        numres = int(s.line[14:17])
        residues = ASCIIString[]
        while chain == current_chain
            sernum = int(s.line[8:10])
            for i in 20:4:68
                residue = s.line[i:i+2]
                residue == "   " && break
                push!(residues, residue)
            end
            readline(s)
            chain = s.line[12]
        end
        @assert length(residues) == numres
        seqres[current_chain] = residues
    end
    seqres
end

function parse_modres(s)
    idcode = s.line[8:11]
    resname = strip(s.line[13:15])
    chainID = s.line[17]
    seqnum = int(s.line[19:22])
    icode = s.line[23]
    stdres = strip(s.line[25:27])
    comment = strip(s.line[30:70])
end

# Heterogen Section
# -----------------

function parse_heterogen_section(s)
    hets = Het[]
    while s.record === :HET
        push!(hets, parse_het(s))
        readline(s)
    end
    hetnams = Hetnam[]
    while s.record === :HETNAM
        push!(hetnams, parse_hetnam(s))
        readline(s)
    end
    hetsyns = Hetsyn[]
    while s.record === :HETSYN
        push!(hetsyns, parse_hetsyn(s))
        readline(s)
    end
    formulae = Formula[]
    while s.record === :FORMUL
        push!(formulae, parse_formul(s))
        readline(s)
    end
    Heterogen(hets, hetnams, hetsyns, formulae)
end

function parse_het(s)
    hetID = strip(s.line[8:10])
    chainID = s.line[13]
    seqnum = int(s.line[14:17])
    icode = s.line[18]
    numhetatoms = int(s.line[21:25])
    text = strip(s.line[31:70])
    Het(hetID, chainID, seqnum, icode, numhetatoms, text)
end

function parse_hetnam(s)
    # TODO: continuation
    hetID = strip(s.line[12:14])
    text = rstrip(s.line[16:70])
    Hetnam(hetID, text)
end

function parse_hetsyn(s)
    # TODO: continuation
    hetID = strip(s.line[12:14])
    hetsyn = rstrip(s.line[16:70])
    Hetsyn(hetID, hetsyn)
end

function parse_formul(s)
    compnum = int(s.line[9:10])
    hetID = strip(s.line[13:15])
    # TODO: continuation
    asterisk = s.line[19] == '*'
    text = rstrip(s.line[20:70])
    Formula(compnum, hetID, asterisk, text)
end

# Coordinate Section
# ------------------

function parse_coodinate_section(s)
    models = Model[]
    @assert s.n_atom_records == 0
    while s.record === :MODEL || s.record === :ATOM
        multimodel = s.record === :MODEL
        if multimodel
            mserial = parse_model(s)
            readline(s)
        else
            mserial = 0
        end
        atoms = Atom[]
        chains = (Char,Vector{Atom})[]
        local atom::Atom
        while (multimodel && s.record !== :ENDMDL) || (!multimodel && !(s.record === :CONECT || s.record === :MASTER))
            if s.record === :ATOM || s.record === :HETATM
                atom = parse_atom(s)
                s.n_atom_records += 1
                s.max_atom_serial = max(atom.serial, s.max_atom_serial)
                push!(atoms, atom)
            elseif s.record === :ANISOU
                serial, name, altloc, resname, chainID, resseq, icode, mtx, element, charge = parse_anisou(s)
                # assuming an ATOM record is immediately followed by an associated ANISOU record
                @assert serial == atom.serial
                @assert name == atom.name
                @assert altloc == atom.altloc
                @assert resname == atom.resname
                @assert chainID == atom.chainID
                @assert resseq == atom.resseq
                @assert icode == atom.icode
                @assert element == atom.element
                @assert charge === atom.charge
                atom.anisou = mtx
            elseif s.record === :TER
                # end of a chain
                serial, resname, chainID, resseq, icode = parse_ter(s)
                @assert !isempty(atoms) "found TER record with no ATOM/HETATM records"
                @assert serial == atom.serial + 1
                @assert resname == atom.resname
                @assert chainID == atom.chainID
                @assert resseq == atom.resseq
                @assert icode == atom.icode
                push!(chains, (chainID, atoms))
                atoms = Atom[]
            else
                if multimodel
                    error("invalid record MODEL and ENDMDL: $(string(s.record))")
                else
                    error("invalid record before CONECT: $(string(s.record))")
                end
            end
            readline(s)
        end
        push!(models, Model(mserial, chains, atoms))
        if multimodel
            @assert s.record === :ENDMDL
            readline(s)
        end
    end
    Coordinate(models)
end

function parse_model(s)
    # model serial number
    int(s.line[11:14])
end

function parse_atom(s)
    hetero = s.record === :HETATM
    serial = int(s.line[7:11])
    name = strip(s.line[13:16])
    altloc = s.line[17]
    resname = strip(s.line[18:20])
    chainID = s.line[22]
    resseq = int(s.line[23:26])
    icode = s.line[27]
    x = float32(s.line[31:38])
    y = float32(s.line[39:46])
    z = float32(s.line[47:54])
    occupancy = float32(s.line[55:60])
    tempfactor = float32(s.line[61:66])
    element = parse_element(s.line[77:78])
    # e.g. 2+, 1-
    charge = s.line[79:80] == "  " ? NaN32 : float32(s.line[80:-1:79])
    Atom(hetero, serial, name, altloc, resname, chainID, resseq, icode,
         x, y, z, occupancy, tempfactor, element, charge)
end

parse_hetatm(s) = parse_atom(s)

function parse_anisou(s)
    serial = int(s.line[7:11])
    name = strip(s.line[13:16])
    altloc = s.line[17]
    resname = strip(s.line[18:20])
    chainID = s.line[22]
    resseq = int(s.line[23:26])
    icode = s.line[27]
    u00 = int(s.line[29:35])
    u11 = int(s.line[36:42])
    u22 = int(s.line[43:49])
    u01 = int(s.line[50:56])
    u02 = int(s.line[57:63])
    u12 = int(s.line[64:70])
    mtx = [u00 u01 u02;
           u01 u11 u12;
           u02 u12 u22]
    element = parse_element(s.line[77:78])
    charge = s.line[79:80] == "  " ? NaN32 : float32(s.line[80:-1:79])
    return serial, name, altloc, resname, chainID, resseq, icode, mtx, element, charge
end

function parse_ter(s)
    # serial, resName, chainID, resSeq, iCode
    int(s.line[7:11]), s.line[18:20], s.line[22], int(s.line[23:26]), s.line[27]
end

# Connectivity Section
# --------------------

function parse_connectivity_section(s)
    Connectivity(parse_conect(s))
end

function parse_conect(s)
    connection = spzeros(Bool, s.max_atom_serial, s.max_atom_serial)
    while s.record === :CONECT
        a = int(s.line[7:11])
        for i in 12:5:27
            num = s.line[i:i+4]
            num == "     " && break
            b = int(num)
            connection[a,b] = true
        end
        readline(s)
    end
    connection
end

# Bookkeeping Section
# -------------------

function parse_bookkeeping_section(s)
    parse_master(s)
    readline(s)
    parse_end(s)
end

function parse_master(s)
    # TODO: check
    n_remark = int(s.line[11:15])
    n_het = int(s.line[21:25])
    n_helix = int(s.line[26:30])
    n_sheet = int(s.line[31:35])
    # deprecated: n_turn = int(s.line[36:40])
    n_site = int(s.line[41:45])
    n_xform = int(s.line[46:50])
    n_coord = int(s.line[51:55])
    n_ter = int(s.line[56:60])
    n_conect = int(s.line[61:65])
    n_seq = int(s.line[66:70])
    # 148L doesn't match???
    #@assert s.n_atoms == n_coord
end

function parse_end(s)
    @assert s.record === :END
end


text:
# %% Comparison with analytically resolved models
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

datasets = Vector(undef, 6)
labels = Vector{String}(undef, 6)

N = Int(1e5)

datasets[1] = :roessler_chaotic
labels[1] = "Rössler"

datasets[2] = :henonheiles_chaotic
labels[2] = "Hénon-Heiles"

datasets[3] = :lorenz96_chaotic
labels[3] = "Lorenz96 (\$D\$=8)"

datasets[4] = :henon_chaotic
labels[4] = "Hénon map"
datasets[5] = :towel_chaotic
labels[5] = "Towel map"
datasets[6] = :coupled_logistics
labels[6] = "CLM (\$D\$=8)"

qH = 1
qC = 2
Cmethod = "standard" # bueno or standard. Decides εmax for correlation sum.

eHs, eCs, Hs, Cs = [Vector{Float64}[] for i in 1:4]
for i in 1:length(datasets)
    # These are the parameters that change between the different elements of the plot
    data = datasets[i]
    # Here we simply pack all parameters into a dictionary
    # (other parameters are (probably) globals)
    params = @strdict N qH qC data
    if (data == :lorenz96_chaotic) || (data == :henonheiles_chaotic)
        params["D"] = 8
    end
    if Cmethod ≠ "standard"
        params["Cmethod"] = Cmethod
    end

    # This is the main call that calculates everything
    output, s = produce_or_load(
        datadir("main"), params, make_C_H;
        prefix = string(data), suffix = "jld2", force = false,
        ignores = ["data"], storepatch = false,
    )
    @unpack eH, eC, H, C = output
    push!(eHs, eH); push!(Hs, H); push!(eCs, eC); push!(Cs, C)
end

legendtitle = "known dynamical rule"

fig, axs = mainplot(
    Hs, Cs, eHs, eCs, labels, legendtitle;
    qH, qC, tol = 0.25,

    # This chooses the kind of fit we do to the linear region.
    # The option `FractalDimension.logarithmic_corrected_fit_lsqfit`
    # is the method suggested by Sprott.
    # The alternative is `FractalDimension.linear_regression_fit_glm`,
    # which is standard linear regression.
    # By default, the correction is only done for the correlation sum method,
    # because this is what the Sprott paper says, but also because we found
    # (by experimentation) that this correction makes the entropy method
    # less accurate when used. For the correlation sum method the
    # accuracy of the corrected version is indeed better.

    # Feel free to experiment by altering any of these two:
    dimension_fit_C = FractalDimension.logarithmic_corrected_fit_lsqfit,
    dimension_fit_H = FractalDimension.linear_regression_fit_glm,
)

wsave(plotsdir("paper", "different_systems_$N"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
# An implicitly nested structure - represents columns, tables, multiple levels of grouped-by tables...
struct Tabular{I, T}
    index::I
    data::T
end

@inline index(t::Tabular) = t.index

@inline indices(t::Tabular) = _indices((index(t),), first(t))

@inline _indices(out::Tuple, t::Tabular) = _indices((out..., index(t)), first(t))
@inline _indices(out::Tuple, x) = out


text:
using FastL2LiR
using Test

@testset "FastL2LiR.jl" begin

    # Array size tests

    X = rand(10, 200)
    Y = rand(10, 100)
    Xp = rand(8, 200)

    model = fit(X, Y, 100.0)
    pred = predict(model, Xp)

    @test size(model.W) == (200, 100)
    @test size(model.b) == (1, 100)
    @test size(pred) == (8, 100)

    model = fit(X, Y, 100.0, 100)

    @test size(model.W) == (200, 100)
    @test size(model.b) == (1, 100)

    X = rand(10, 200)
    Y = rand(10, 5, 4, 3)
    Xp = rand(8, 200)

    model = fit(X, Y, 100.0)
    pred = predict(model, Xp)

    @test size(model.W) == (200, 5, 4, 3)
    @test size(model.b) == (1, 5, 4, 3)
    @test size(pred) == (8, 5, 4, 3)

    model = fit(X, Y, 100.0, 100)

    @test size(model.W) == (200, 5, 4, 3)
    @test size(model.b) == (1, 5, 4, 3)
end


text:
"""
    A tree in which every node has 0 or 1 children
"""
struct OneTree
    nodes::Vector{Int}
end
AbstractTrees.treekind(::Type{OneTree}) = AbstractTrees.IndexedTree()
AbstractTrees.siblinglinks(::Type{OneTree}) = AbstractTrees.StoredSiblings()
AbstractTrees.relative_state(t::OneTree, _, __::Int) = 1
Base.getindex(t::OneTree, idx) = t.nodes[idx]
AbstractTrees.childindices(tree::OneTree, node::Int) =
    (ret = (node == 0 || tree[node] == 0) ? () : (tree[node],))
AbstractTrees.children(tree::OneTree) = AbstractTrees.children(tree, tree)
AbstractTrees.rootstate(tree::OneTree) = 1
AbstractTrees.printnode(io::IO, t::OneTree) =
    AbstractTrees.printnode(io::IO, t[AbstractTrees.rootstate(t)])
Base.eltype(::Type{<:TreeIterator{OneTree}}) = Int
Base.IteratorEltype(::Type{<:TreeIterator{OneTree}}) = Base.HasEltype()

@testset "OneTree" begin
    ot = OneTree([2,3,4,0])
    @test repr_tree(ot) == "2\n└─ 3\n   └─ 4\n      └─ 0\n"
    @test @inferred(collect(Leaves(ot))) == [0]
    @test eltype(collect(Leaves(ot))) === Int
    @test collect(PreOrderDFS(ot)) == [2,3,4,0]
    @test collect(PostOrderDFS(ot)) == [0,4,3,2]
end

"""
    Stores an explicit parent for some other kind of tree
"""
struct ParentTree{T}
    tree::T
    parents::Vector{Int}
end
AbstractTrees.treekind(::Type{ParentTree{T}}) where {T} = AbstractTrees.treekind(T)
AbstractTrees.parentlinks(::Type{ParentTree{T}}) where {T} = AbstractTrees.StoredParents()
AbstractTrees.siblinglinks(::Type{ParentTree{T}}) where {T} = AbstractTrees.siblinglinks(T)
AbstractTrees.relative_state(t::ParentTree, x, r::Int) =
    AbstractTrees.relative_state(t.tree, x, r)
Base.getindex(t::ParentTree, idx) = t.tree[idx]
AbstractTrees.childindices(tree::ParentTree, node::Int) = AbstractTrees.childindices(tree.tree, node)
AbstractTrees.children(tree::ParentTree) = AbstractTrees.children(tree, tree)
AbstractTrees.rootstate(tree::ParentTree) = AbstractTrees.rootstate(tree.tree)
AbstractTrees.parentind(tree::ParentTree, node::Int) = tree.parents[node]
AbstractTrees.printnode(io::IO, t::ParentTree) =
    AbstractTrees.printnode(io::IO, t[AbstractTrees.rootstate(t)])

@testset "ParentTree" begin
    ot = OneTree([2,3,4,0])
    pt = ParentTree(ot,[0,1,2,3])

    @test repr_tree(pt) == "2\n└─ 3\n   └─ 4\n      └─ 0\n"
    @test collect(Leaves(pt)) == [0]
    @test collect(PreOrderDFS(pt)) == [2,3,4,0]
    @test collect(PostOrderDFS(pt)) == [0,4,3,2]

    # Test modification while iterating over PreOrderDFS
    a = [1,[2,[3]]]
    b = treemap!(PreOrderDFS(a)) do node
        !isa(node, Vector) && return node
        ret = pushfirst!(copy(node),0)
        # And just for good measure stomp over the old node to make sure nothing
        # is cached.
        empty!(node)
        ret
    end
    @test b == Any[0,1,Any[0,2,[0,3]]]
end

struct IntTree
    num::Int
    children::Vector{IntTree}
end
==(x::IntTree,y::IntTree) = x.num == y.num && x.children == y.children
AbstractTrees.children(itree::IntTree) = itree.children
Base.eltype(::Type{<:TreeIterator{IntTree}}) = IntTree
Base.IteratorEltype(::Type{<:TreeIterator{IntTree}}) = Base.HasEltype()
AbstractTrees.nodetype(::IntTree) = IntTree

@testset "IntTree" begin
    itree = IntTree(1, [IntTree(2, IntTree[])])
    iter = Leaves(itree)
    @test @inferred(first(iter)) == IntTree(2, IntTree[])
    val, state = iterate(iter)
    @test Base.return_types(iterate, Tuple{typeof(iter), typeof(state)}) ==
        [Union{Nothing, Tuple{IntTree,typeof(state)}}]
end

#=
@test treemap(PostOrderDFS(tree)) do ind, x, children
    IntTree(isa(x,Int) ? x : mapreduce(x->x.num,+,0,children),
        isempty(children) ? IntTree[] : children)
end == IntTree(6,[IntTree(1,IntTree[]),IntTree(5,[IntTree(2,IntTree[]),IntTree(3,IntTree[])])])
=#

@test collect(PostOrderDFS([])) == Any[[]]


text:
include("../examples/rijke.jl")
include("../src/lss.jl")
using OrdinaryDiffEq
using JLD
using SharedArrays
function Rijke_tangent_sensitivity(beta)
	s = [beta, 0.2]
	n = 2000
	d = N

	du_trj = zeros(d, d, n)
	dJ_trj = zeros(2,d,n)
	f_trj = zeros(d,n)
	X_trj = zeros(d,n)

	dJds = zeros(2)
	vsh = zeros(d, n)
	nRunup = 1000000
	u = Rijke_ODE(rand(d),s,nRunup)

	for i = 1:n
    vel_p_i = view(u, 1:2*Ng)
    dJ_trj[1,1:2*Ng,i] = 0.5*vel_p_i
    dJ_trj[2,1:Ng,i] = cjpixf
    X_trj[:,i] = perturbation(u, s, 1.e-6)
    f!(view(f_trj,:,i), u, s, 1.)
    du_trj[:,:, i] = dRijke(u, s, 1.e-6)
    u = Rijke_ODE(u, s, 1)
	end
	y, xi = lss(du_trj, X_trj, f_trj, s, 2)
	dJds = compute_sens(y, xi, dJ_trj, f_trj)
	return dJds[1]
end
function average_sensitivity(beta)
	n_samples = 100
	dJds_avg = @distributed (+) for i = 1:n_samples
        Rijke_tangent_sensitivity(beta)
	end
	dJds_avg /= n_samples
end
function optimize()
    nNewton = 20
    beta_path = zeros(nNewton)
    dJds = zeros(nNewton)
    gamma = 0.1
    beta_path[1] = 6.5
    beta = beta_path[1]
    for nN = 1:nNewton
	dJds[nN] = average_sensitivity(beta)
       	beta = beta - gamma*dJds[nN]
	beta_path[nN] = beta
    end
    save("../data/param_optim/beta_dJds_1.jld",
	 "beta_path", beta_path,
	 "dJds_path", dJds)
end
function compute_Eac(beta)
	nRunup = 1000000
	s = [beta, 0.2]
	d = N
	u = Rijke_ODE(rand(d),s,nRunup)
	nSteps = 20000
	Eac = 0.
	for i = 1:nSteps
		if(i==1000)
			println("executing step")
		end
		u = Rijke_ODE(u, s, 1)
		Eac += 0.25*norm(u[1:2*Ng])^2.0/nSteps
	end
	return Eac
end
function compute_Eac_along_path()
	X = load(string("../data/param_optim/",
			"beta_dJds_1.jld"))
	Y = load(string("../data/param_optim/",
			"beta_dJds.jld"))
	beta1 = X["beta_path"]
	n1, =  size(beta1)
	beta = Y["beta_path"]
	beta2 = [beta1; beta]
	n, = size(beta2)
	Eac = SharedArray{Float64}(n)
	t = @distributed for i=1:n
		Eac[i] = compute_Eac(beta2[i])
	end
	wait(t)
	save("../data/param_optim/Eac_along_path.jld",
		"Eac1", Eac[1:n1],
		"Eac", Eac[n1+1:end],
		"beta", beta,
		"beta1", beta1)
end

#=
    n_samples = 1000
    dJacds_ens, dJrayds_ens = zeros(2,n_samples),
                    zeros(2,n_samples)
    up = zeros(d)
    um = zeros(d)
    eps = 1.e-3
    n = 100 # produces a growth of ≈ e
    for i = 1:n_samples
        @time dJacds_ens[:,i] = Zygote.gradient(s ->Jac(u, s, n), s)[1]
    dJrayds_ens[:,i] = Zygote.gradient(s ->Jray(u, s, n),
                        	   s)[1]
    u .= Rijke(u, [7.0, 0.2], n)
    end
=#

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
function sample_field(N)
    uhat = zeros(ComplexF64,2*N); # preallocate space

    ξ = randn(N); # generate the random variables

    # construct the eigenvalues
    k = 1:N;
    λ = @. 1/(π*k)^2;

    # fill in the nonzero entries
    # NOTE we need to multiply by 2 *N for FFT scaling
    @. uhat[2:N+1] = 2 * N * sqrt(λ) * sqrt(2) * ξ;

    # invert and get the relevant imaginary part
    u = imag.(ifft(uhat))[N+2:end];
    return u
end


text:
abstract type IntegratorRK{dType, tType} <: DeterministicIntegrator{dType, tType} end
abstract type IntegratorPRK{dType, tType} <: IntegratorRK{dType, tType} end

@inline equation(integrator::IntegratorRK) = integrator.params.equ
@inline timestep(integrator::IntegratorRK) = integrator.params.Δt
@inline tableau(integrator::IntegratorRK)  = integrator.params.tab

@inline nstages(integrator::IntegratorRK)  = nstages(tableau(integrator))
@inline eachstage(integrator::IntegratorRK) = 1:nstages(integrator)