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

@test integrate((b*x+c*x^2)^7//2, x) == :(1//16*c^-1*(b*x+c*x^2)^7//2*(b+2*c*x)+35//16384*b^8*c^-9//2*arctanh(x*c^(1/2)*(b*x+c*x^2)^-1//2)+-35//16384*b^6*c^-4*(b*x+c*x^2)^(1/2)*(b+2*c*x)+-7//384*b^2*c^-2*(b*x+c*x^2)^5//2*(b+2*c*x)+35//6144*b^4*c^-3*(b*x+c*x^2)^3//2*(b+2*c*x))
@test integrate((4*x^2+3*im*x)^7//2, x) == :((4*x^2+3*im*x)^7//2*(1//8*x+3//64*im)+1//2048*(4*x^2+3*im*x)^5//2*(63im+168x)+1//131072*(4*x^2+3*im*x)^3//2*(2835im+7560x)+1//4194304*(4*x^2+3*im*x)^(1/2)*(76545im+204120x)+229635//16777216*im*arcsin(1+-8//3*im*x))
@test integrate((4*x^2+3*im*x)^5//2, x) == :((4*x^2+3*im*x)^5//2*(1//6*x+1//16*im)+1//1024*(4*x^2+3*im*x)^3//2*(45im+120x)+1//32768*(4*x^2+3*im*x)^(1/2)*(1215im+3240x)+3645//131072*im*arcsin(1+-8//3*im*x))
@test integrate((4*x^2+3*im*x)^3//2, x) == :((4*x^2+3*im*x)^3//2*(1//4*x+3//32*im)+1//1024*(4*x^2+3*im*x)^(1/2)*(81im+216x)+243//4096*im*arcsin(1+-8//3*im*x))
@test integrate((4*x^2+3*im*x)^(1/2), x) == :((4*x^2+3*im*x)^(1/2)*((1/2)*x+3//16*im)+9//64*im*arcsin(1+-8//3*im*x))
@test integrate((-4*x^2+3x)^7//2, x) == :(-229635//16777216*arcsin(1+-8//3*x)+-1*(-4*x^2+3x)^7//2*(3//64+-1//8*x)+-1//2048*(-4*x^2+3x)^5//2*(63+-168x)+-1//131072*(-4*x^2+3x)^3//2*(2835+-7560x)+-1//4194304*(-4*x^2+3x)^(1/2)*(76545+-204120x))
@test integrate((-4*x^2+3x)^5//2, x) == :(-3645//131072*arcsin(1+-8//3*x)+-1*(-4*x^2+3x)^5//2*(1//16+-1//6*x)+-1//1024*(-4*x^2+3x)^3//2*(45+-120x)+-1//32768*(-4*x^2+3x)^(1/2)*(1215+-3240x))
@test integrate((-4*x^2+3x)^3//2, x) == :(-243//4096*arcsin(1+-8//3*x)+-1*(-4*x^2+3x)^3//2*(3//32+-1//4*x)+-1//1024*(-4*x^2+3x)^(1/2)*(81+-216x))
@test integrate((-4*x^2+3x)^(1/2), x) == :(-9//64*arcsin(1+-8//3*x)+(-4*x^2+3x)^(1/2)*(-3//16+(1/2)*x))
@test integrate((-1*x^2+6x)^(1/2), x) == :(-9//2*arcsin(1+-1//3*x)+(-1*x^2+6x)^(1/2)*(-3//2+(1/2)*x))
@test integrate((-9*x^2+5x)^(1/2), x) == :(-25//216*arcsin(1+-18//5*x)+(-9*x^2+5x)^(1/2)*(-5//36+(1/2)*x))
@test integrate((x+-1*x^2)^3//2, x) == :(-3//128*arcsin(1+-2x)+(x+-1*x^2)^(1/2)*(-3//64+3//32*x)+-1*(x+-1*x^2)^3//2*(1//8+-1//4*x))
@test integrate((x^2+4x)^(1/2), x) == :(-4*arctanh(x*(x^2+4x)^-1//2)+(x^2+4x)^(1/2)*(1+(1/2)*x))
@test integrate((x^2+-8x)^(1/2), x) == :(-16*arctanh(x*(x^2+-8x)^-1//2)+(x^2+-8x)^(1/2)*(-2+(1/2)*x))
@test integrate((x^2+-1x)^(1/2), x) == :(-1//4*arctanh(x*(x^2+-1x)^-1//2)+(x^2+-1x)^(1/2)*(-1//4+(1/2)*x))
@test integrate((b*x+c*x^2)^-7//2, x) == :(-1//5*b^-2*(b*x+c*x^2)^-5//2*(2b+4*c*x)+-256//15*b^-6*c^2*(b*x+c*x^2)^-1//2*(b+2*c*x)+32//15*c*b^-4*(b*x+c*x^2)^-3//2*(b+2*c*x))
@test integrate((4*x^2+3*im*x)^-1//2, x) == :((1/2)*im*arcsin(1+-8//3*im*x))
@test integrate((4*x^2+3*im*x)^-3//2, x) == :(1//9*(4*x^2+3*im*x)^-1//2*(6im+16x))
@test integrate((4*x^2+3*im*x)^-5//2, x) == :(1//27*(4*x^2+3*im*x)^-3//2*(6im+16x)+1//243*(4*x^2+3*im*x)^-1//2*(192im+512x))
@test integrate((4*x^2+3*im*x)^-7//2, x) == :(1//45*(4*x^2+3*im*x)^-5//2*(6im+16x)+1//1215*(4*x^2+3*im*x)^-3//2*(384im+1024x)+1//10935*(4*x^2+3*im*x)^-1//2*(12288im+32768x))
@test integrate((-4*x^2+3x)^-1//2, x) == :(-1//2*arcsin(1+-8//3*x))
@test integrate((-4*x^2+3x)^-3//2, x) == :(-1//9*(-4*x^2+3x)^-1//2*(6+-16x))
@test integrate((-4*x^2+3x)^-5//2, x) == :(-1//27*(-4*x^2+3x)^-3//2*(6+-16x)+-1//243*(-4*x^2+3x)^-1//2*(192+-512x))
@test integrate((-4*x^2+3x)^-7//2, x) == :(-1//45*(-4*x^2+3x)^-5//2*(6+-16x)+-1//1215*(-4*x^2+3x)^-3//2*(384+-1024x)+-1//10935*(-4*x^2+3x)^-1//2*(12288+-32768x))
@test integrate((b*x+-1*b^2*x^2)^-1//2, x) == :(-1*b^-1*arcsin(1+-2*b*x))
@test integrate((b*x+b^2*x^2)^-1//2, x) == :(2*b^-1*arctanh(b*x*(b*x+b^2*x^2)^-1//2))
@test integrate((-1*x^2+6x)^-1//2, x) == :(-1*arcsin(1+-1//3*x))
@test integrate((x^2+4x)^-1//2, x) == :(2*arctanh(x*(x^2+4x)^-1//2))
@test integrate((x^2+-2x)^-1//2, x) == :(2*arctanh(x*(x^2+-2x)^-1//2))
@test integrate((b*x+c*x^2)^4//3, x) == :(3//22*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-4//3*(-1*c*x*b^-2*(b+c*x))^4//3*(b*x+c*x^2)^4//3*(b+2*c*x)+3//55*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-4//3*(-1*c*x*b^-2*(b+c*x))^1//3*(b*x+c*x^2)^4//3*(b+2*c*x)+1//55*2^1//3*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^-4//3*(2+-1*3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^4//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^1//3, x) == :(3//10*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-1//3*(-1*c*x*b^-2*(b+c*x))^1//3*(b*x+c*x^2)^1//3*(b+2*c*x)+1//10*2^1//3*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^-1//3*(2+-1*3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^1//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^-2//3, x) == :(2^1//3*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^2//3*(2+-1*3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^-2//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^-5//3, x) == :((1/2)*c^-1*(-1*c*b^-2*(b*x+c*x^2))^5//3*(-1*c*x*b^-2*(b+c*x))^-2//3*(b*x+c*x^2)^-5//3*(3b+6*c*x)+2^1//3*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^5//3*(2+-1*3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^-5//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^-8//3, x) == :(1//5*c^-1*(-1*c*b^-2*(b*x+c*x^2))^8//3*(-1*c*x*b^-2*(b+c*x))^-5//3*(b*x+c*x^2)^-8//3*(3b+6*c*x)+1//5*c^-1*(-1*c*b^-2*(b*x+c*x^2))^8//3*(-1*c*x*b^-2*(b+c*x))^-2//3*(b*x+c*x^2)^-8//3*(21b+42*c*x)+14//5*2^1//3*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^8//3*(2+-1*3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^-8//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^5//3, x) == :(3//26*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-5//3*(-1*c*x*b^-2*(b+c*x))^5//3*(b*x+c*x^2)^5//3*(b+2*c*x)+15//364*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-5//3*(-1*c*x*b^-2*(b+c*x))^2//3*(b*x+c*x^2)^5//3*(b+2*c*x)+-1//364*2^2//3*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-5//3*(b*x+c*x^2)^5//3*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(15b+30*c*x)+5//182*2^1//6*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^-5//3*(b+2*c*x)^-1*(b*x+c*x^2)^5//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2))+-15//728*2^2//3*3^1//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^-5//3*(2+3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^5//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.E(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^2//3, x) == :(3//14*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-2//3*(-1*c*x*b^-2*(b+c*x))^2//3*(b*x+c*x^2)^2//3*(b+2*c*x)+-1//14*2^2//3*c^-1*(-1*c*b^-2*(b*x+c*x^2))^-2//3*(b*x+c*x^2)^2//3*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(3b+6*c*x)+1//7*2^1//6*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^-2//3*(b+2*c*x)^-1*(b*x+c*x^2)^2//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2))+-3//28*2^2//3*3^1//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^-2//3*(2+3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^2//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.E(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^-1//3, x) == :(-1//2*2^2//3*c^-1*(-1*c*b^-2*(b*x+c*x^2))^1//3*(b*x+c*x^2)^-1//3*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(3b+6*c*x)+2^1//6*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^1//3*(b+2*c*x)^-1*(b*x+c*x^2)^-1//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2))+-3//4*2^2//3*3^1//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^1//3*(2+3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^-1//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.E(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^-4//3, x) == :(c^-1*(-1*c*b^-2*(b*x+c*x^2))^4//3*(-1*c*x*b^-2*(b+c*x))^-1//3*(b*x+c*x^2)^-4//3*(3b+6*c*x)+3*2^2//3*c^-1*(-1*c*b^-2*(b*x+c*x^2))^4//3*(b*x+c*x^2)^-4//3*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(b+2*c*x)+-2*2^1//6*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^4//3*(b+2*c*x)^-1*(b*x+c*x^2)^-4//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2))+3//2*2^2//3*3^1//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^4//3*(2+3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^-4//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.E(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^-7//3, x) == :((1/2)*c^-1*(-1*c*b^-2*(b*x+c*x^2))^7//3*(-1*c*x*b^-2*(b+c*x))^-1//3*(b*x+c*x^2)^-7//3*(15b+30*c*x)+1//4*c^-1*(-1*c*b^-2*(b*x+c*x^2))^7//3*(-1*c*x*b^-2*(b+c*x))^-4//3*(b*x+c*x^2)^-7//3*(3b+6*c*x)+(1/2)*2^2//3*c^-1*(-1*c*b^-2*(b*x+c*x^2))^7//3*(b*x+c*x^2)^-7//3*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(15b+30*c*x)+-5*2^1//6*3^3//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^7//3*(b+2*c*x)^-1*(b*x+c*x^2)^-7//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.F(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2))+15//4*2^2//3*3^1//4*b^2*c^-1*((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3+2*2^1//3*(-1*c*x*b^-2*(b+c*x))^2//3))^(1/2)*(-1*(1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-2*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3))^-1//2*(-1*c*b^-2*(b*x+c*x^2))^7//3*(2+3^(1/2))^(1/2)*(b+2*c*x)^-1*(b*x+c*x^2)^-7//3*(1+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)*Elliptic.E(arcsin((1+-1*3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)^-1*(1+3^(1/2)+-1*2^2//3*(-1*c*x*b^-2*(b+c*x))^1//3)),-7+4*3^(1/2)))
@test integrate((b*x+c*x^2)^5//4, x) == :(1//7*c^-1*(b*x+c*x^2)^5//4*(b+2*c*x)+-5//84*b^2*c^-2*(b*x+c*x^2)^1//4*(b+2*c*x)+5//168*2^(1/2)*b^5*c^-3*(-1*c*b^-2*(b*x+c*x^2))^3//4*(b*x+c*x^2)^-3//4*Elliptic.F((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((b*x+c*x^2)^3//4, x) == :(1//5*c^-1*(b*x+c*x^2)^3//4*(b+2*c*x)+-3//20*2^(1/2)*b^3*c^-2*(-1*c*b^-2*(b*x+c*x^2))^1//4*(b*x+c*x^2)^-1//4*Elliptic.E((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((b*x+c*x^2)^1//4, x) == :(1//3*c^-1*(b*x+c*x^2)^1//4*(b+2*c*x)+-1//6*2^(1/2)*b^3*c^-2*(-1*c*b^-2*(b*x+c*x^2))^3//4*(b*x+c*x^2)^-3//4*Elliptic.F((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((b*x+c*x^2)^-1//4, x) == :(b*2^(1/2)*c^-1*(-1*c*b^-2*(b*x+c*x^2))^1//4*(b*x+c*x^2)^-1//4*Elliptic.E((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((b*x+c*x^2)^-3//4, x) == :(2*b*2^(1/2)*c^-1*(-1*c*b^-2*(b*x+c*x^2))^3//4*(b*x+c*x^2)^-3//4*Elliptic.F((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((b*x+c*x^2)^-5//4, x) == :(-1*b^-2*(b*x+c*x^2)^-1//4*(4b+8*c*x)+4*2^(1/2)*b^-1*(-1*c*b^-2*(b*x+c*x^2))^1//4*(b*x+c*x^2)^-1//4*Elliptic.E((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((b*x+c*x^2)^-9//4, x) == :(-1//5*b^-2*(b*x+c*x^2)^-5//4*(4b+8*c*x)+48//5*c*b^-4*(b*x+c*x^2)^-1//4*(b+2*c*x)+-48//5*c*2^(1/2)*b^-3*(-1*c*b^-2*(b*x+c*x^2))^1//4*(b*x+c*x^2)^-1//4*Elliptic.E((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((b*x+c*x^2)^-13//4, x) == :(-1//9*b^-2*(b*x+c*x^2)^-9//4*(4b+8*c*x)+-448//15*b^-6*c^2*(b*x+c*x^2)^-1//4*(b+2*c*x)+112//45*c*b^-4*(b*x+c*x^2)^-5//4*(b+2*c*x)+448//15*2^(1/2)*b^-5*c^2*(-1*c*b^-2*(b*x+c*x^2))^1//4*(b*x+c*x^2)^-1//4*Elliptic.E((1/2)*arcsin(1+2*c*x*b^-1),2))
@test integrate((a+c*x^2)^4, x) == :(x*a^4+1//9*c^4*x^9+4//3*c*a^3*x^3+4//7*a*c^3*x^7+6//5*a^2*c^2*x^5)
@test integrate((a+c*x^2)^3, x) == :(x*a^3+1//7*c^3*x^7+c*a^2*x^3+3//5*a*c^2*x^5)
@test integrate((a+c*x^2)^2, x) == :(x*a^2+1//5*c^2*x^5+2//3*a*c*x^3)
@test integrate(a+c*x^2, x) == :(a*x+1//3*c*x^3)
@test integrate((a+c*x^2)^-1, x) == :(a^-1//2*c^-1//2*arctan(x*a^-1//2*c^(1/2)))
@test integrate((a+c*x^2)^-2, x) == :((1/2)*x*a^-1*(a+c*x^2)^-1+(1/2)*a^-3//2*c^-1//2*arctan(x*a^-1//2*c^(1/2)))
@test integrate((a+c*x^2)^-3, x) == :(1//4*x*a^-1*(a+c*x^2)^-2+3//8*x*a^-2*(a+c*x^2)^-1+3//8*a^-5//2*c^-1//2*arctan(x*a^-1//2*c^(1/2)))
@test integrate((a+c*x^2)^5//2, x) == :(1//6*x*(a+c*x^2)^5//2+5//16*x*a^2*(a+c*x^2)^(1/2)+5//16*a^3*c^-1//2*arctanh(x*c^(1/2)*(a+c*x^2)^-1//2)+5//24*a*x*(a+c*x^2)^3//2)
@test integrate((a+c*x^2)^3//2, x) == :(1//4*x*(a+c*x^2)^3//2+3//8*a*x*(a+c*x^2)^(1/2)+3//8*a^2*c^-1//2*arctanh(x*c^(1/2)*(a+c*x^2)^-1//2))
@test integrate((a+c*x^2)^(1/2), x) == :((1/2)*x*(a+c*x^2)^(1/2)+(1/2)*a*c^-1//2*arctanh(x*c^(1/2)*(a+c*x^2)^-1//2))
@test integrate((a+c*x^2)^-1//2, x) == :(c^-1//2*arctanh(x*c^(1/2)*(a+c*x^2)^-1//2))
@test integrate((a+c*x^2)^-3//2, x) == :(x*a^-1*(a+c*x^2)^-1//2)
@test integrate((a+c*x^2)^-5//2, x) == :(1//3*x*a^-1*(a+c*x^2)^-3//2+2//3*x*a^-2*(a+c*x^2)^-1//2)
@test integrate((a+c*x^2)^-7//2, x) == :(1//5*x*a^-1*(a+c*x^2)^-5//2+4//15*x*a^-2*(a+c*x^2)^-3//2+8//15*x*a^-3*(a+c*x^2)^-1//2)
@test integrate((a+c*x^2)^-9//2, x) == :(1//7*x*a^-1*(a+c*x^2)^-7//2+6//35*x*a^-2*(a+c*x^2)^-5//2+8//35*x*a^-3*(a+c*x^2)^-3//2+16//35*x*a^-4*(a+c*x^2)^-1//2)
@test integrate((4+9*x^2+12x)^3//2, x) == :((4+9*x^2+12x)^3//2*(1//6+1//4*x))
@test integrate((4+9*x^2+12x)^(1/2), x) == :((4+9*x^2+12x)^(1/2)*(1//3+(1/2)*x))
@test integrate((4+9*x^2+12x)^-1//2, x) == :(1//3*(4+9*x^2+12x)^-1//2*(2+3x)*log(2+3x))
@test integrate((4+9*x^2+12x)^-3//2, x) == :(-1*(12+18x)^-1*(4+9*x^2+12x)^-1//2)
@test integrate((4+-12x+9*x^2)^(1/2), x) == :((4+-12x+9*x^2)^(1/2)*(-1//3+(1/2)*x))
@test integrate((4+-12x+9*x^2)^-1//2, x) == :(-1//3*(4+-12x+9*x^2)^-1//2*(2+-3x)*log(2+-3x))
@test integrate((-4+-9*x^2+12x)^(1/2), x) == :((-4+-9*x^2+12x)^(1/2)*(-1//3+(1/2)*x))
@test integrate((-4+-9*x^2+12x)^-1//2, x) == :(-1//3*(-4+-9*x^2+12x)^-1//2*(2+-3x)*log(2+-3x))
@test integrate((-4+-12x+-9*x^2)^(1/2), x) == :((-4+-12x+-9*x^2)^(1/2)*(1//3+(1/2)*x))
@test integrate((-4+-12x+-9*x^2)^-1//2, x) == :(1//3*(-4+-12x+-9*x^2)^-1//2*(2+3x)*log(2+3x))
@test integrate((b*x+c*x^2+1//4*c^-1*(-1+b^2))^5, x) == :(-5//896*c^-6*(1+-1b+-2*c*x)^7+-5//2304*c^-6*(1+-1b+-2*c*x)^9+-1//22528*c^-6*(1+-1b+-2*c*x)^11+1//384*c^-6*(1+-1b+-2*c*x)^6+1//2048*c^-6*(1+-1b+-2*c*x)^10+5//1024*c^-6*(1+-1b+-2*c*x)^8)
@test integrate((b*x+c*x^2+1//4*c^-1*(-4+b^2))^5, x) == :(-5//56*c^-6*(2+-1b+-2*c*x)^7+-5//576*c^-6*(2+-1b+-2*c*x)^9+-1//22528*c^-6*(2+-1b+-2*c*x)^11+1//12*c^-6*(2+-1b+-2*c*x)^6+1//1024*c^-6*(2+-1b+-2*c*x)^10+5//128*c^-6*(2+-1b+-2*c*x)^8)
@test integrate((b*x+c*x^2+1//4*c^-1*(-9+b^2))^5, x) == :(-405//896*c^-6*(3+-1b+-2*c*x)^7+-5//256*c^-6*(3+-1b+-2*c*x)^9+-1//22528*c^-6*(3+-1b+-2*c*x)^11+3//2048*c^-6*(3+-1b+-2*c*x)^10+81//128*c^-6*(3+-1b+-2*c*x)^6+135//1024*c^-6*(3+-1b+-2*c*x)^8)
@test integrate((b*x+c*x^2+1//4*c^-1*(-16+b^2))^5, x) == :(-10//7*c^-6*(4+-1b+-2*c*x)^7+-5//144*c^-6*(4+-1b+-2*c*x)^9+-1//22528*c^-6*(4+-1b+-2*c*x)^11+1//512*c^-6*(4+-1b+-2*c*x)^10+5//16*c^-6*(4+-1b+-2*c*x)^8+8//3*c^-6*(4+-1b+-2*c*x)^6)
@test integrate((2+3*x^2+4x)^-1, x) == :((1/2)*2^(1/2)*arctan((1/2)*2^(1/2)*(2+3x)))
@test integrate((4+x^2+-2*x*3^(1/2))^-1, x) == :(-1*arctan(3^(1/2)+-1x))
@test integrate((2+-3*x^2+4x)^-1, x) == :(-1//10*10^(1/2)*arctanh(1//10*10^(1/2)*(2+-3x)))
@test integrate((2+3*x^2+5x)^-1, x) == :(-1*log(1+x)+log(2+3x))
@test integrate((2+-3*x^2+5x)^-1, x) == :(-1//7*log(2+-1x)+1//7*log(1+3x))
@test integrate((1+2*x^2+pi*x)^-1, x) == :(-2*(-8+pi^2)^-1//2*arctanh((-8+pi^2)^-1//2*(pi+4x)))
@test integrate((1+-2*x^2+pi*x)^-1, x) == :(-2*(8+pi^2)^-1//2*arctanh((8+pi^2)^-1//2*(pi+-4x)))
@test integrate((1+3*x^2+pi*x)^-1, x) == :(2*(12+-1*pi^2)^-1//2*arctan((12+-1*pi^2)^-1//2*(pi+6x)))
@test integrate((1+-3*x^2+pi*x)^-1, x) == :(-2*(12+pi^2)^-1//2*arctanh((12+pi^2)^-1//2*(pi+-6x)))
@test integrate((a+b*x^2+c*x)^-1, x) == :(2*(-1*c^2+4*a*b)^-1//2*arctan((-1*c^2+4*a*b)^-1//2*(c+2*b*x)))
@test integrate((b+b*x^2+2*a*x)^-1, x) == :(-1*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*x)))
@test integrate((b+-1*b*x^2+2*a*x)^-1, x) == :(-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*x)))
@test integrate((2+3*x^2+4x)^-2, x) == :((8+12*x^2+16x)^-1*(2+3x)+3//8*2^(1/2)*arctan((1/2)*2^(1/2)*(2+3x)))
@test integrate((2+-3*x^2+4x)^-2, x) == :(-1*(40+-60*x^2+80x)^-1*(2+-3x)+-3//200*10^(1/2)*arctanh(1//10*10^(1/2)*(2+-3x)))
@test integrate((2+3*x^2+5x)^-2, x) == :(-6*log(2+3x)+6*log(1+x)+-1*(2+3*x^2+5x)^-1*(5+6x))
@test integrate((2+-3*x^2+5x)^-2, x) == :(-6//343*log(2+-1x)+6//343*log(1+3x)+-1*(98+-147*x^2+245x)^-1*(5+-6x))
@test integrate((a+b*x^2+c*x)^-2, x) == :((-1*c^2+4*a*b)^-1*(a+b*x^2+c*x)^-1*(c+2*b*x)+4*b*(-1*c^2+4*a*b)^-3//2*arctan((-1*c^2+4*a*b)^-1//2*(c+2*b*x)))
@test integrate((b+b*x^2+2*a*x)^-2, x) == :((1/2)*b*(a^2+-1*b^2)^-3//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*x))+-1*(-2*b^2+2*a^2)^-1*(b+b*x^2+2*a*x)^-1*(a+b*x))
@test integrate((b+-1*b*x^2+2*a*x)^-2, x) == :(-1*(2*a^2+2*b^2)^-1*(b+-1*b*x^2+2*a*x)^-1*(a+-1*b*x)+-1//2*b*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(a+-1*b*x)))
@test integrate((x^2+(a*b^-1)^(2*n^-1)+-2*x*(a*b^-1)^(n^-1)*cos(n^-1*(pi+-2*pi*k)))^-1, x) == :(-1*(a*b^-1)^(-1*n^-1)*arctan(-1*x*(a*b^-1)^(-1*n^-1)*csc(n^-1*(pi+-2*pi*k))+cot(n^-1*(pi+-2*pi*k)))*csc(n^-1*(pi+-2*pi*k)))
@test integrate((1+x^2+2*x*cos(1//7))^-1, x) == :(arctan((x+cos(1//7))*csc(1//7))*csc(1//7))
@test integrate((1+x^2+2*x*cos(1//7*pi))^-1, x) == :(arctan(x*csc(1//7*pi)+cot(1//7*pi))*csc(1//7*pi))
@test integrate((5+-6x+9*x^2)^(1/2), x) == :(2//3*arcsinh(-1//2+3//2*x)+(5+-6x+9*x^2)^(1/2)*(-1//6+(1/2)*x))
@test integrate((3+-4x+-4*x^2)^(1/2), x) == :((3+-4x+-4*x^2)^(1/2)*(1//4+(1/2)*x)+arcsin(1/2+x))
@test integrate((-8+6x+9*x^2)^(1/2), x) == :(-3//2*arctanh((-8+6x+9*x^2)^-1//2*(1+3x))+(-8+6x+9*x^2)^(1/2)*(1//6+(1/2)*x))
@test integrate((2+3*x^2+4x)^(1/2), x) == :((2+3*x^2+4x)^(1/2)*(1//3+(1/2)*x)+1//9*3^(1/2)*arcsinh((1/2)*2^(1/2)*(2+3x)))
@test integrate((2+-3*x^2+4x)^(1/2), x) == :((2+-3*x^2+4x)^(1/2)*(-1//3+(1/2)*x)+-5//9*3^(1/2)*arcsin(1//10*10^(1/2)*(2+-3x)))
@test integrate((2+3*x^2+5x)^(1/2), x) == :((2+3*x^2+5x)^(1/2)*(5//12+(1/2)*x)+-1//72*3^(1/2)*arctanh(1//6*3^(1/2)*(2+3*x^2+5x)^-1//2*(5+6x)))
@test integrate((2+-3*x^2+5x)^(1/2), x) == :((2+-3*x^2+5x)^(1/2)*(-5//12+(1/2)*x)+-49//72*3^(1/2)*arcsin(5//7+-6//7*x))
@test integrate((-2+3*x^2+4x)^(1/2), x) == :((-2+3*x^2+4x)^(1/2)*(1//3+(1/2)*x)+-5//9*3^(1/2)*arctanh(1//3*3^(1/2)*(-2+3*x^2+4x)^-1//2*(2+3x)))
@test integrate((-2+-3*x^2+4x)^(1/2), x) == :((-2+-3*x^2+4x)^(1/2)*(-1//3+(1/2)*x)+1//9*3^(1/2)*arctan(1//3*3^(1/2)*(-2+-3*x^2+4x)^-1//2*(2+-3x)))
@test integrate((-2+3*x^2+5x)^(1/2), x) == :((-2+3*x^2+5x)^(1/2)*(5//12+(1/2)*x)+-49//72*3^(1/2)*arctanh(1//6*3^(1/2)*(-2+3*x^2+5x)^-1//2*(5+6x)))
@test integrate((-2+-3*x^2+5x)^(1/2), x) == :((-2+-3*x^2+5x)^(1/2)*(-5//12+(1/2)*x)+-1//72*3^(1/2)*arcsin(5+-6x))
@test integrate((5+-6x+9*x^2)^-1//2, x) == :(1//3*arcsinh(-1//2+3//2*x))
@test integrate((3+-4x+-4*x^2)^-1//2, x) == :((1/2)*arcsin(1/2+x))
@test integrate((-8+6x+9*x^2)^-1//2, x) == :(1//3*arctanh((-8+6x+9*x^2)^-1//2*(1+3x)))
@test integrate((2+3*x^2+4x)^-1//2, x) == :(1//3*3^(1/2)*arcsinh((1/2)*2^(1/2)*(2+3x)))
@test integrate((2+-3*x^2+4x)^-1//2, x) == :(-1//3*3^(1/2)*arcsin(1//10*10^(1/2)*(2+-3x)))
@test integrate((2+3*x^2+5x)^-1//2, x) == :(1//3*3^(1/2)*arctanh(1//6*3^(1/2)*(2+3*x^2+5x)^-1//2*(5+6x)))
@test integrate((2+-3*x^2+5x)^-1//2, x) == :(-1//3*3^(1/2)*arcsin(5//7+-6//7*x))
@test integrate((-2+3*x^2+4x)^-1//2, x) == :(1//3*3^(1/2)*arctanh(1//3*3^(1/2)*(-2+3*x^2+4x)^-1//2*(2+3x)))
@test integrate((-2+-3*x^2+4x)^-1//2, x) == :(-1//3*3^(1/2)*arctan(1//3*3^(1/2)*(-2+-3*x^2+4x)^-1//2*(2+-3x)))
@test integrate((-2+3*x^2+5x)^-1//2, x) == :(1//3*3^(1/2)*arctanh(1//6*3^(1/2)*(-2+3*x^2+5x)^-1//2*(5+6x)))
@test integrate((-2+-3*x^2+5x)^-1//2, x) == :(-1//3*3^(1/2)*arcsin(5+-6x))
@test integrate((b*x+c*x^2+1//4*c^-1*(b^2+4c))^-1//2, x) == :(c^-1//2*arcsinh((1/2)*c^-1//2*(b+2*c*x)))
@test integrate((b*x+-1*c*x^2+1//4*c^-1*(-1*b^2+4c))^-1//2, x) == :(-1*c^-1//2*arcsin((1/2)*c^-1//2*(b+-2*c*x)))
@test integrate((b*x+-1*c*x^2+1//4*c^-1*(c+-1*b^2))^-1//2, x) == :(-1*c^-1//2*arcsin(c^-1//2*(b+-2*c*x)))
@test integrate((2+x^2+3x)^-3//2, x) == :((2+x^2+3x)^-1//2*(-6+-4x))
@test integrate((27+-24x+4*x^2)^-3//2, x) == :(1//9*(27+-24x+4*x^2)^-1//2*(3+-1x))
@test integrate(x*(5+-1*x^2+-4x)^-3//2, x) == :(1//9*(5+-1*x^2+-4x)^-1//2*(5+-2x))
@test integrate((5+-1*x^2+-4x)^-5//2, x) == :(1//27*(5+-1*x^2+-4x)^-3//2*(2+x)+1//243*(5+-1*x^2+-4x)^-1//2*(4+2x))
@test integrate((3+4x)^p, x) == :((3+4x)^(1+p)*(4+4p)^-1)

max_stars_repo_path: test/1 Algebraic functions/1.2 Trinomial products/1.2.1 Quadratic/1.2.1.1 (a+b x+c x^2)^p.jl


text:
module ProgressBarsModule

"""
Customisable progressbar decorator for iterators.
Copied from https://github.com/cloud-oak/ProgressBars.jl to allow for custom modifications.
Usage:
> using ProgressBars
> for i in ProgressBar(1:10)
> ....
> end
"""

import Printf: @sprintf

EIGHTS = Dict(
    0 => ' ', 1 => '▏', 2 => '▎', 3 => '▍', 4 => '▌', 5 => '▋', 6 => '▊', 7 => '▉', 8 => '█'
)

# Split this because UTF-8 indexing is horrible otherwise
# IDLE = collect("◢◤ ")
IDLE = collect("╱   ")

PRINTING_DELAY = 0.05 * 1e9

# export ProgressBar, tqdm, set_description, set_postfix, set_multiline_postfix
"""
Decorate an iterable object, returning an iterator which acts exactly
like the original iterable, but prints a dynamically updating
progressbar every time a value is requested.
"""
mutable struct ProgressBar
    wrapped::Any
    total::Int
    current::Int
    width::Int
    fixwidth::Bool
    leave::Bool
    start_time::UInt
    last_print::UInt
    description::AbstractString
    postfix::NamedTuple
    extra_lines::Int
    last_extra_lines::Int
    multilinepostfix::AbstractString
    mutex::Threads.SpinLock

    function ProgressBar(wrapped::Any; total::Int=-2, width=nothing, leave=true)
        this = new()
        this.wrapped = wrapped
        if width === nothing
            this.width = displaysize(stdout)[2]
            this.fixwidth = false
        else
            this.width = width
            this.fixwidth = true
        end
        this.leave = leave
        this.start_time = time_ns()
        this.last_print = this.start_time - 2 * PRINTING_DELAY
        this.description = ""
        this.postfix = NamedTuple()
        this.multilinepostfix = ""
        this.extra_lines = 0
        this.last_extra_lines = 0
        this.mutex = Threads.SpinLock()
        this.current = 0

        if total == -2  # No total given
            try
                this.total = length(wrapped)
            catch
                this.total = -1
            end
        else
            this.total = total
        end

        return this
    end
end

# Keep the old name as an alias
tqdm = ProgressBar

function format_time(seconds)
    if isfinite(seconds)
        mins, s = divrem(round(Int, seconds), 60)
        h, m = divrem(mins, 60)
    else
        h = 0
        m = Inf
        s = Inf
    end
    if h != 0
        return @sprintf("%02d:%02d:%02d", h, m, s)
    else
        return @sprintf("%02d:%02d", m, s)
    end
end

function display_progress(t::ProgressBar)
    seconds = (time_ns() - t.start_time) * 1e-9
    iteration = t.current - 1

    elapsed = format_time(seconds)
    speed = iteration / seconds
    if seconds == 0
        # Dummy value of 1 it/s if no time has elapsed
        speed = 1
    end
    iterations_per_second = @sprintf("%.1f it/s", speed)

    barwidth = t.width - 2 # minus two for the separators

    postfix_string = postfix_repr(t.postfix)

    # Reset Cursor to beginning of the line
    for line in 1:(t.extra_lines)
        move_up_1_line()
    end
    go_to_start_of_line()

    if t.description != ""
        barwidth -= length(t.description) + 1
        print(t.description * " ")
    end

    if (t.total <= 0)
        status_string = "$(t.current)it $elapsed [$iterations_per_second$postfix_string]"
        barwidth -= length(status_string) + 1
        if barwidth < 0
            barwidth = 0
        end

        print("┣")
        print(join(IDLE[1 + ((i + t.current) % length(IDLE))] for i in 1:barwidth))
        print("┫ ")
        print(status_string)
    else
        ETA = (t.total - t.current) / speed

        percentage_string = string(@sprintf("%.1f%%", t.current / t.total * 100))

        eta = format_time(ETA)
        status_string = "$(t.current)/$(t.total) [$elapsed<$eta, $iterations_per_second$postfix_string]"

        barwidth -= length(status_string) + length(percentage_string) + 1
        if barwidth < 0
            barwidth = 0
        end

        cellvalue = t.total / barwidth
        full_cells, remain = divrem(t.current, cellvalue)

        print(percentage_string)
        print("┣")
        print(repeat("█", Int(full_cells)))
        if (full_cells < barwidth)
            part = Int(floor(9 * remain / cellvalue))
            print(EIGHTS[part])
            print(repeat(" ", Int(barwidth - full_cells - 1)))
        end

        print("┫ ")
        print(status_string)
    end
    multiline_postfix_string = newline_to_spaces(t.multilinepostfix, t.width)
    t.last_extra_lines = t.extra_lines
    t.extra_lines = ceil(Int, length(multiline_postfix_string) / t.width) + 1
    print(multiline_postfix_string)
    return println() #Newline is required for Python to read in.
end

erase_to_end_of_line() = print("\033[K")
move_up_1_line() = print("\033[1A")
move_down_1_line() = print("\033[1B")
go_to_start_of_line() = print("\r")
erase_line() = begin
    go_to_start_of_line()
    erase_to_end_of_line()
end

# Clear the progress bar
function clear_progress(t::ProgressBar)
    # Reset cursor, fill width with empty spaces, and then reset again
    if t.last_extra_lines > t.extra_lines
        for line in 1:(t.last_extra_lines - t.extra_lines)
            move_down_1_line()
        end
    end
    for line in 1:max(t.extra_lines, t.last_extra_lines)
        erase_line()
        move_up_1_line()
    end
    return erase_line()
end

function set_multiline_postfix(t::ProgressBar, postfix::AbstractString)
    mistakenly_used_newline_at_start = postfix[1] == '\n' && length(postfix) > 1
    if mistakenly_used_newline_at_start
        postfix = postfix[2:end]
    end
    return t.multilinepostfix = postfix
end

function postfix_repr(postfix::NamedTuple)::AbstractString
    return join(map(tpl -> ", $(tpl[1]): $(tpl[2])", zip(keys(postfix), postfix)))
end

function Base.iterate(iter::ProgressBar)
    if displaysize(stdout)[2] != iter.width && !iter.fixwidth
        iter.width = displaysize(stdout)[2]
        print("\n"^(iter.extra_lines + 2))
    end
    iter.start_time = time_ns() - PRINTING_DELAY
    iter.current = 0
    display_progress(iter)
    return iterate(iter.wrapped)
end

function Base.iterate(iter::ProgressBar, s)
    if displaysize(stdout)[2] != iter.width && !iter.fixwidth
        iter.width = displaysize(stdout)[2]
        print("\n"^(iter.extra_lines + 2))
    end
    iter.current += 1
    if (time_ns() - iter.last_print > PRINTING_DELAY)
        display_progress(iter)
        iter.last_print = time_ns()
    end
    state = iterate(iter.wrapped, s)
    if state === nothing
        if iter.total > 0
            iter.current = iter.total
        end
        display_progress(iter)
        if iter.leave
            println()
        else
            clear_progress(iter)
        end
        return nothing
    end
    return state
end
Base.length(iter::ProgressBar) = length(iter.wrapped)
Base.eltype(iter::ProgressBar) = eltype(iter.wrapped)

function Base.unsafe_getindex(iter::ProgressBar, index::Int64)
    """
    Base.unsafe_getindex is used by the `Threads.@threads for ... in ...` macro
    in julia 1.3.
    This wrapper will do weird things when used directly.
    """
    item = Base.unsafe_getindex(iter.wrapped, index)
    lock(iter.mutex)
    iter.current += 1
    if time_ns() - iter.last_print > PRINTING_DELAY
        display_progress(iter)
        iter.last_print = time_ns()
    elseif iter.current == iter.total
        # Reached end of iteration
        display_progress(iter)
        if iter.leave
            println()
        else
            clear_progress(iter)
        end
    end
    unlock(iter.mutex)
    return item
end

function Base.firstindex(iter::ProgressBar)
    lock(iter.mutex)
    iter.start_time = time_ns() - PRINTING_DELAY
    iter.current = 0
    display_progress(iter)
    unlock(iter.mutex)
    return Base.firstindex(iter.wrapped)
end

function Base.getindex(iter::ProgressBar, index::Int64)
    """
    Base.getindex is used by the `Threads.@threads for ... in ...` macro
    from julia 1.4 on.
    This wrapper will do weird things when used directly.
    """
    item = Base.getindex(iter.wrapped, index)
    lock(iter.mutex)
    iter.current += 1
    if time_ns() - iter.last_print > PRINTING_DELAY
        display_progress(iter)
        iter.last_print = time_ns()
    elseif iter.current == iter.total
        # Reached end of iteration
        display_progress(iter)
        if iter.leave
            println()
        else
            clear_progress(iter)
        end
    end
    unlock(iter.mutex)
    return item
end

function newline_to_spaces(string, terminal_width)
    new_string = ""
    width_cumulator = 0
    for c in string
        if c == '\n'
            spaces_required = terminal_width - width_cumulator
            new_string *= " "^spaces_required
            width_cumulator = 0
        else
            new_string *= c
            width_cumulator += 1
        end
        if width_cumulator == terminal_width
            width_cumulator = 0
        end
    end
    return new_string
end

# end

end

ext: jl


text:
#=
A unival tree (which stands for "universal value") is a tree where all nodes under it have the same value.

Given the root to a binary tree, count the number of unival subtrees.

For example, the following tree has 5 unival subtrees:

   0
  / \
 1   0
    / \
   1   0
  / \
 1   1
=#

using Test
include("Solutions/problem8_univalsubtreescounter.jl")

# The example
#=
    0
   / \
  1   0
     / \
    1   0
   / \
  1   1
=#

node = Node(1, Node(1),Node(1))
node = Node(0, node, Node(0))
node = Node(0, Node(1), node)

@test count_unival_subtrees(node) == 5

# Other example
#=      0
      /  \
     0   1
    / \
   1  1
  /\  /\
 1 1 1 1
 =#

node = Node(1, Node(1), Node(1))
node = Node(0, node, node)
node = Node(0, node, Node(1))

@test count_unival_subtrees(node) == 8

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using LinearAlgebra         # 引入线性代数库
x=zeros(10)                 # 10 个元素的零矢量，列矢量
a=ones(3)                   # 3 个元素的矢量，全填充 1.0，列矢量
H=[1.0 2.0 3.0              # 直接构造矩阵
   4.0 5.0 6.0
   7.0 8.0 9.0]
H[2,1]                      # 4.0
ndims(H)                    # H 的维度，2
eltype(H)                   # H 中元素的类型，Float64
length(H)                   # H 中元素数目，9
size(H)                     # H 每一维度的大小所构成的元组，(3,3)
sum(H, dims=2)              # 对第2维求和，相似地，该用法可以用于以下函数： mean, std,
                            # prod, minimum, maximum, any, all
                            # 使用统计函数前需要 using Statistics
H[:,1]                      # 矩阵的第1列
H[1,:]                      # 矩阵的第1行（仍是列矢量）
Hc=H; Hc[1,1]=100.0; H[1,1] # 100.0，数组的赋值是按引用复制
Hc=copy(H)                  # 这样才是真的复制一新的矩阵
H*a                         # 矩阵乘法
exp(H)                      # 矩阵的函数
a*a                         # 出错的矩阵乘法
                            # 列矢量和列矢量的乘法行列不对应无法进行
a⋅a                         # 但是可以做点积
                            # ⋅ 用 \cdot<Tab> 可以打出来（在 REPL、Jupyter 中）
                            # 或者借助插件（见前文开发环境部分）
a×a                         # 还可以做叉积
sum(a)
H'                          # 共轭转置，相当于 adjoint(H)
typeof(H')                  # 注意共轭转置是惰性求值
                            # 它并不直接分配一块内存并写入计算的结果
                            # 这对计算效率有很大的帮助
G=rand(1000,1000); G=G+G'   # 构造一个实对称阵
eigvals(G)                  # 这样的小矩阵，我们可以直接求本征值
eigvecs(G)                  # 本征矢

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
export DiaSymSemiseparableChol, trinv

struct DiaSymSemiseparableChol <: SymSemiseparableCholesky
    n::Int64
    p::Int64
    U::AbstractArray
    W::AbstractArray
    ds::AbstractArray
end

# Constuctors
DiaSymSemiseparableChol(U::AbstractArray, W::AbstractArray, ds::AbstractArray) = DiaSymSemiseparableChol(size(U,1),size(U,2),U,W,ds);

# Mappings
mul!(y::AbstractArray, L::DiaSymSemiseparableChol, 		            b::AbstractArray) =   dss_tri_mul!(y, L.U,   L.W,   L.ds,   b);
mul!(y::AbstractArray, L::AdjointOperator{DiaSymSemiseparableChol}, b::AbstractArray) =  dssa_tri_mul!(y, L.A.U, L.A.W, L.A.ds, b);
inv!(y::AbstractArray, L::DiaSymSemiseparableChol, 		      		b::AbstractArray) =   dss_forward!(y, L.U,   L.W,   L.ds,   b);
inv!(y::AbstractArray, L::AdjointOperator{DiaSymSemiseparableChol}, b::AbstractArray) =  dssa_backward!(y, L.A.U, L.A.W, L.A.ds, b);
newlogdet(L::DiaSymSemiseparableChol) = dss_logdet(L.ds)
newlogdet(L::AdjointOperator{DiaSymSemiseparableChol}) = dss_logdet(L.A.ds)

#### Inverse of a EGRQSCholesky using the Cholesky factorization ####
function inv(L::DiaSymSemiseparableChol)
	return L'\(L\Diagonal(ones(L.n)))
end
function inv(L::DiaSymSemiseparableChol, b::AbstractArray)
	return L'\(L\b)
end

# Traces, norms and determinants
function fro_norm_L(L::DiaSymSemiseparableChol)
	return sum(squared_norm_cols(L.U, L.W, L.ds))
end

function trinv(L::DiaSymSemiseparableChol)
	dbar = L.ds;
	Y, Z = dss_create_yz(L.U, L.W, dbar)
	return sum(squared_norm_cols(Y,Z, dbar.^(-1)))
end
function tr(L::DiaSymSemiseparableChol)
	return sum(L.ds)
end
function tr(Ky::DiaSymSemiseparableChol, K::SymSemiseparable)
	n = Ky.n;
	p = Ky.p;
	c = Ky.ds;
	U = K.U;
	V = K.V;
	Y, Z = dss_create_yz(Ky.U, Ky.W, Ky.ds);
	b = 0;
	P = zeros(p,p);
	R = zeros(p,p);
	for k = 1:Ky.n
		yk = Y[k,:];
		zk = Z[k,:];
		uk = U[k,:];
		vk = V[k,:];
		cki = c[k]^(-1);
		b += yk'*P*yk + 2*yk'*R*uk*cki + uk'*vk*(cki^2);
		P += ((uk'*vk)*zk)*zk' + zk*(R*uk)' + (R*uk)*zk';
		R += zk*vk';
	end
	return b
end

########################################################################
#### Cholesky factorization of Higher-order quasiseparable matrices ####
########################################################################

#### Creating W and ds s.t. L = tril(UW',-1) + diag(ds) ####
function dss_create_wdbar(U::AbstractArray, V::AbstractArray, d::AbstractArray)
    n, m = size(U);
    P  = zeros(m, m);
    W  = zeros(n, m);
    dbar = zeros(n);
    for i = 1:n
        tmpU  = U[i,:]
        tmpW  = V[i,:] - P*tmpU;
        tmpds = sqrt(tmpU'*tmpW + d[i]);
        tmpW  = tmpW/tmpds
        W[i,:] = tmpW;
        dbar[i] = tmpds;
        P += tmpW*tmpW';
    end
    return W, dbar
end
#### Multiplying with Ld ####
function dss_tri_mul!(Y::AbstractArray,U::AbstractArray,W::AbstractArray,
                     ds::AbstractArray,X::AbstractArray)
     n, m = size(U)
     mx = size(X, 2)
     Wbar = zeros(m, mx)
     for i = 1:n
         tmpW = W[i,:]
         tmpU = U[i,:]
         tmpX = X[i:i,:];
         Y[i,:] = tmpU'*Wbar + ds[i]*tmpX;
         Wbar  +=  tmpW*tmpX;
     end
end
#### Adjoint of Ld ####
function dssa_tri_mul!(Y::AbstractArray,U::AbstractArray,W::AbstractArray,
                      ds::AbstractArray,X::AbstractArray)
     n, m = size(U)
     mx = size(X, 2)
     Ubar = zeros(m,mx)
     for i = n:-1:1
         tmpW = W[i,:]
         tmpU = U[i,:]
         tmpX = X[i:i,:]
         Y[i,:] = tmpW'*Ubar + ds[i]*tmpX;
         Ubar = Ubar + tmpU*tmpX;
     end
end
#### Forward substitution ####
function dss_forward!(X::AbstractArray,U::AbstractArray,W::AbstractArray,
                     ds::AbstractArray,B::AbstractArray)
    n, m = size(U)
    mx = size(B,2)
    Wbar = zeros(m,mx);
    for i = 1:n
        tmpU = U[i,:];
        tmpW = W[i,:];
        X[i:i,:] = (B[i:i,:] - tmpU'*Wbar)/ds[i];
        Wbar += tmpW .* X[i:i,:];
    end
end
#### Backward substitution ####
function dssa_backward!(X::AbstractArray,U::AbstractArray,W::AbstractArray,
                       ds::AbstractArray,B::AbstractArray)
    n, m = size(U)
    mx = size(B,2)
    Ubar = zeros(m,mx);
    for i = n:-1:1
        tmpU = U[i,:];
        tmpW = W[i,:];
        X[i:i,:] = (B[i:i,:] - tmpW'*Ubar)/ds[i];
        Ubar += tmpU .* X[i:i,:];
    end
end

#### Squared norm of columns of L = tril(UW',-1) + diag(dbar) ####
function squared_norm_cols(U::AbstractArray,W::AbstractArray,
                        dbar::AbstractArray)
    n, m = size(U)
    P = zeros(m, m)
    c = zeros(n)
    for i = n:-1:1
        tmpW = W[i,:]
        tmpU = U[i,:]
        c[i]  = dbar[i]^2 + tmpW'*P*tmpW
        P += tmpU*tmpU'
    end
    return c
end
#### Implicit inverse of  L = tril(UW',-1) + diag(dbar) ####
function dss_create_yz(U::AbstractArray, W::AbstractArray,
                    dbar::AbstractArray)
    n, m = size(U)
    Y = zeros(n,m)
    Z = zeros(n,m)
    dss_forward!(Y, U, W, dbar, U)
    dssa_backward!(Z, U, W, dbar, W)
    # Probably best not to use inv
    return Y, Z*inv(U'*Z - I)
end
#### Log-determinant ####
function dss_logdet(d::AbstractArray)
    return sum(log.(d))
end


text:
using Test
using Theta
using LinearAlgebra

@testset "Theta" begin
    include("theta_test.jl")
    include("lattice_test.jl")
    include("ellipsoid_test.jl")
    include("riemann_matrix_test.jl")
    include("characteristics_test.jl")
    include("accola_test.jl")
    include("fgsm_test.jl")
    include("schottky4_test.jl")
    include("schottky_null_test.jl")
end


text:
# misc functions
function epoch(;data,tbl,τ,sfreq,kwargs...)
    epoch(data,tbl,τ,sfreq;kwargs...)
end

function epoch(data::Array{T,1},tbl::DataFrame,τ::Tuple{Number,Number},sfreq;kwargs...) where T <:Union{Missing,Number}
    data_r = reshape(data,(1,:))
    epoch(data_r,tbl,τ,sfreq;kwargs... )
end

epoch(data::Matrix, tbl::DataFrame, τ::Vector,sfreq; eventtime) = epoch(data, tbl, Tuple(τ...), sfreq; eventtime=eventtime)
epoch(data::Matrix, tbl::DataFrame, τ::Tuple{Number,Number},sfreq; eventtime::String="latency") = epoch(data, tbl, τ, sfreq; eventtime=Symbol(eventtime))

function epoch(data::Array{T,2},tbl::DataFrame,τ::Tuple{Number,Number},sfreq;eventtime::Symbol=:latency) where T <:Union{Missing,Number}
# data: channels x times

# partial taken from EEG.jl

    numEpochs = size(tbl,1)
    τ = round_times(τ,sfreq)

    times = range(τ[1],stop=τ[2],step=1 ./sfreq)
    lenEpochs = length(times)
    numChans = size(data,1)
    epochs = Array{Union{Missing,Float64}}(missing,Int(numChans),Int(lenEpochs), Int(numEpochs))

    # User feedback
    @debug "Creating epochs: $numChans x $lenEpochs x $numEpochs"

    for si = 1:size(tbl,1)
        #eventonset = tbl[si,eventtime] # in samples
        #d_start = eventonset
        d_start = Int(round(tbl[si,eventtime] - sum(times.<0)))
        d_end =  Int(round(tbl[si,eventtime]+lenEpochs-1 - sum(times.<0)))
        e_start = 1
        e_end = lenEpochs
            #println("d: $(size(data)),e: $(size(epochs)) | $d_start,$d_end,$e_start,$e_end | $(tbl[si,eventtime])")
        if d_start<1
            e_start = e_start + (-d_start+1)
            d_start = 1
        end
        if d_end >size(data,2)
            e_end = e_end - (d_end-size(data,2))
            d_end = size(data,2)
        end
        #println("d: $(size(data)),e: $(size(epochs)) | $d_start,$d_end,$e_start,$e_end | $(tbl[si,eventtime])")
        epochs[:, e_start:e_end,si] = data[:,d_start:d_end]
    end
    return(epochs,times)
end

function round_times(τ,sfreq)
    # function to round τ to sfreq samples. This specifies the epoch length.
    # its a function to be the same for epoch & timeexpanded analyses
    return round.(τ .* sfreq) ./ sfreq

end
function dropMissingEpochs(X,y)
    missingIx = .!any(ismissing.(y),dims=(1,2))
    print(size(missingIx))
    goodIx = dropdims(missingIx,dims=(1,2))
    return X[goodIx,:],Array{Float64}(y[:,:,goodIx])
end

"""
$(SIGNATURES)
Flatten a 1D array from of a 2D/3D array. Also drops the empty dimension
"""
function linearize(x::AbstractArray)
    # used in condense_long to generate the long format
    return dropdims(reshape(x,:,1),dims=2)
end
function linearize(x::String)
    return x
end

"""
$(SIGNATURES)
Equates the length of data and designmatrix by cutting the shorter one

The reason we need this is because when generating the designmatrix, we do not know how long the data actually are. We only assume that event-latencies are synchronized with the data
"""
function zeropad(X,data::AbstractArray{T,2})where {T<:Union{Missing, <:Number}}
    @debug("2d zeropad")
    if size(X,1) > size(data,2)
        X = X[1:size(data,2),:]
    else
        data =data[:,1:size(X,1)]
    end
    return X,data
end
function zeropad(X,data::AbstractArray{T,3})where {T<:Union{Missing, <:Number}}
    @debug("3d zeropad")
    if size(X,1) > size(data,3)
        X = X[1:size(data,3),:]
    else
        data =data[:,:,1:size(X,1)]
    end
    return X,data
end

function clean_data(data::AbstractArray{T,2},winrej::AbstractArray{<:Number,2}) where T <:Union{Float64,Missing}
    data = Array{Union{Float64,Missing}}(data)
    for row in 1:size(winrej,1)
        data[:,Int.(winrej[row,1]:winrej[row,2])] .= missing
    end
    return data
end


text:
module TestLinqMacro

using Test
using DataFrames
using DataFramesMeta
using Statistics
using Random

Random.seed!(100)
n = 100
df = DataFrame(a = rand(1:3, n),
               b = ["a","b","c","d"][rand(1:4, n)],
               x = rand(n))

x = @where(df, :a .> 2, :b .!= "c")
x = @transform(x, y = 10 * :x)
x = @by(x, :b, meanX = mean(:x), meanY = mean(:y))
x = @orderby(x, -:meanX)
x = @select(x, var = :b, :meanX, :meanY)

x1 = @linq transform(where(df, :a .> 2, :b .!= "c"), y = 10 * :x)
x1 = @linq by(x1, :b, meanX = mean(:x), meanY = mean(:y))
x1 = @linq select(orderby(x1, -:meanX), var = :b, :meanX, :meanY)

## chaining
xlinq = @linq df  |>
    where(:a .> 2, :b .!= "c")  |>
    transform(y = 10 * :x)  |>
    by(:b, meanX = mean(:x), meanY = mean(:y))  |>
    orderby(-:meanX)  |>
    select(var = :b, :meanX, :meanY)

@test x == x1
@test x == xlinq

xlinq2 = @linq df  |>
    where(:a .> 2, :b .!= "c")  |>
    transform(y = 10 * :x)  |>
    groupby(:b) |>
    orderby(-mean(:x))  |>
    based_on(meanX = mean(:x), meanY = mean(:y))

@test xlinq2[!, [:meanX, :meanY]] == xlinq[!, [:meanX, :meanY]]

xlinq3 = @linq df  |>
    where(:a .> 2, :b .!= "c")  |>
    transform(y = 10 * :x)  |>
    DataFrames.groupby(:b) |>
    orderby(-mean(:x))  |>
    based_on(meanX = mean(:x), meanY = mean(:y))

@test xlinq3[!, [:meanX, :meanY]] == xlinq[!, [:meanX, :meanY]]

@test (@linq df |> with(:a)) == df.a

@testset "@linq with `cols`" begin
    df = DataFrame(
            a = [1, 2, 3, 4],
            b = ["a", "b", "c", "d"],
            x = [10, 20, 30, 40],
            y = [40, 50, 60, 70]
        )

    a_sym = :a
    b_str = "b"
    x_sym = :x
    y_str = "y"
    xlinq3 = @linq df  |>
        where(cols(a_sym) .> 2, :b .!= "c")  |>
        transform(cols(y_str) = 10 * cols(x_sym))  |>
        DataFrames.groupby(b_str) |>
        orderby(-mean(cols(x_sym)))  |>
        based_on(cols("meanX") = mean(:x), meanY = mean(:y))

    @test isequal(xlinq3, DataFrame(b = "d", meanX = 40.0, meanY = 400.0))
end

end # module


text:
using Flux
using Flux.Tracker
using Flux: @epochs

using JuliaRL
using JuliaRL.Environments
using Random
using LinearAlgebra
using ProgressMeter

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
#=
Utility files for sorting the position vectors and optimizing which
indices of the wavefunction and space to use when doing
the Husimi computation.
=#
export xspan, sparse_sampling

function findclosest(val, A)
    i = 1
    d = abs(val - A[i])
    @inbounds for j in eachindex(A)
        dd = abs(val - A[j])
        if dd < d
            i = j
            d = dd
        end
    end
    return i
end

"""
    xslices(r) → i
Given `r` (a vector of `SVector`s), *which is assumed to be sorted*,
find the indices `i` that denote the changes in `x` slices.
"""
function xslices(r)
    xcol = [1]
    @inbounds for i in 2:length(r)
        r[i][1] == r[i-1][1] && continue
        push!(xcol, i)
    end
    return xcol
end

"""
    xspan(r0, r, toskip) → span
Find the `span` of `r` (the position vectors) that are within ± `toskip`
from the `x` location of `r0`, or within ± `toskip` from the maximum and
minimum `x` of `r0`, if `r0` is a vector of `SVector`s.
"""
function xspan(r0::SVector, ris, toskip)
    xmin = xmax = r0[1]
    xspan(xmin, xmax, ris, toskip)
end
function xspan(r0s::Vector{<:SVector}, ris, toskip)
    xmin, xmax = extrema(r[1] for r in r0s)
    xspan(xmin, xmax, ris, toskip)
end

function xspan(xmin, xmax, ris, toskip)
    xcol = xslices(ris)

    k = 0; i = 1 # left index, for minimum x
    @inbounds for (z, n) in enumerate(xcol)
        if ris[n][1] ≥ xmin - toskip
            i = n
            k = z
            break
        end
    end

    it = 0;  j=xcol[end] # right index, for maximum x
    @inbounds for n in @view xcol[k:end]
        it += 1
        if ris[n][1] ≥ xmax + toskip
            j = n
            break
        end
    end

    N = length(ris)
    # This line ensures that you get all the y elements
    # of the right-most column:
    j = (j==xcol[end] ? N : xcol[it+k+1]-1)
    return i:j
end

# TODO:
# Now we also limit the span to ys (assuming approximately same
# y span for every column)

"""
    sparse_sampling(ris; e = 2, ex = e, ey = e)
Return indices `i` that sample `ris` (the positions) by sparsely
sampling every `ex` x and every `ey` y positions (in integers).
"""
function sparse_sampling(ris; e = 2, ex = e, ey = e)
    xcol = xslices(ris)
    is = Int[]
    sizehint!(is, round(Int, length(ris)/(ex*ey)))
    for j in 1:ex:length(xcol)-1
        append!(is, xcol[j]:ey:xcol[j+1]-1)
    end
    return is
end


text:
module MdvspCg

using GLPK

include("instance.jl")

export Instance,
    pullout_cost,
    pullin_cost,
    deadhead_cost,
    read_huisman_instance

include("main.jl")

export build_main_problem,
    add_column!

include("subproblem.jl")

export optimize!

function gencol()
    inst = read_huisman_instance("../mdvsp-colgen/instances/m4n500s0.inp")
    #inst = read_huisman_instance("../mdvsp-colgen/instances/m8n1500s4.inp")

    main = build_main_problem(inst)
    sp = map(1:inst.num_depots) do k
        Subproblem(inst, k)
    end

    timings = open("timings.csv", "w")
    write(timings, "event,iter,time,bytes,gctime\n")
    function tm_append(msg, iter, tm)
        write(timings, "$(msg),$(iter),$(tm.time),$(tm.bytes),$(tm.gctime)\n")
    end

    # Configures the optimizer
    set_optimizer(main, GLPK.Optimizer)

    for iter in 1:10000
        tm = @timed optimize!(main)
        tm_append("solve_main", iter, tm)

        println("Iter: $iter  Main obj: $(objective_value(main))")

        depot_duals = dual.(main[:depot_cap])
        task_duals = dual.(main[:task_asg])
        newcols = Vector{Any}(undef, inst.num_depots)

        tm = @timed begin
            Threads.@threads for k in 1:inst.num_depots
                newcols[k] = optimize!(sp[k], depot_duals[k], task_duals)
            end
        end
        tm_append("solve_pricing", iter, tm)

        converged = true
        tm = @timed begin
            for k in 1:inst.num_depots
                if !isempty(newcols[k].path)
                    converged = false
                    add_column!(main, k, newcols[k].cost, newcols[k].path)
                end
            end
        end
        tm_append("add_columns", iter, tm)

        if converged
            if !main[:asg_relax]
                println("** COLGEN CONVERGED! **")
                break
            else
                println("*** CHANGING RELAXATION STATUS ***")
                tm = @timed set_main_asg_equals!(main)
                tm_append("set_equals", iter, tm)
            end
        end

        # write_to_file(main, "main.$iter.lp")
    end

    write_to_file(main, "main.final.lp")
    close(timings)

    return main
end

end # module

max_stars_repo_name: afkummer/mdvsp-colgen-study

max_issues_repo_name: afkummer/mdvsp-colgen-study

max_forks_repo_name: afkummer/mdvsp-colgen-study


text:
using Logging

function check(p, i::Integrated, rng=Random.default_rng())
    genAs = [ generate(rng, freeze(i)) for _ in 1:numTests[] ]
    res = findCounterexample(p, genAs)
    res === nothing && return true
    @debug res
    entry = last(res)
    if errored(entry)
        return entry.example, exception(entry)
    else
        return entry.example
    end
end

const ExState = Tuple{<:Exception,Vector{Base.StackTraces.StackFrame}}
struct CheckEntry{T}
    example::T
    exception::Union{Nothing,ExState}
end

errored(ce::CheckEntry) = ce.exception !== nothing
exception(ce::CheckEntry) = errored(ce) && first(ce.exception)

function minimize!(log, f, t::Tree{T}, initEx) where {T}
    ex::Union{Nothing,ExState} = initEx
    while true
        r = root(t)
        @debug "Possible shrink value" r
        push!(log, CheckEntry(r, ex))
        subs = subtrees(t)
        filteredSubtrees = Iterators.filter(first, Iterators.map(f, subs))
        !any(first, filteredSubtrees) && break
        _, t, ex = first(filteredSubtrees)
    end
    infomsg = "$(length(log)) counterexamples found"
    errcount = count(errored, log)
    if !iszero(errcount)
        unique_errors = filter(!isnothing, unique(exception, log))
        distinct_errors = length(unique_errors)
        infomsg *= ", of which $errcount threw $distinct_errors distinct exception types"
        @info infomsg Errors=map(exception, unique_errors)
    else
        @info infomsg
    end
    log
end

function findCounterexample(f, trees::Vector{<:Tree})
    function _f(tree)
        try
            ((!f ∘ root)(tree), tree, nothing)
        catch ex
            trace = stacktrace(catch_backtrace())
            (true, tree, (ex, trace)) # if we threw, the property failed
        end
    end
    checkedTrees = map(_f, trees)
    filter!(first, checkedTrees)
    isempty(checkedTrees) && return nothing
    @info "Found counterexample for '$f', beginning shrinking..."
    _, initTree, initEx = first(checkedTrees)
    log = CheckEntry{eltype(initTree)}[]
    minimize!(log, _f, initTree, initEx)
end


text:
using ProjectAssigner: match
using Test
using DataFrames

@testset "ProjectAssigner.jl" begin
    @testset "Basic match" begin
        m = match(students="students.csv", projects="projects.csv", output="out.csv")

        d = Dict(r[:name]=>r[:project] for r in Tables.rows(m))
        @test d["Alice"] == "Airplane"
        @test d["Bob"] == "Airplane"
        @test d["Charlie"] == "Cubesat"
        @test d["Dale"] == "Cubesat"
    end

    @testset "Force" begin
        m = match(students="students.csv", projects="projects.csv", force=["Charlie"=>"Cubesat"])
        d = Dict(r[:name]=>r[:project] for r in Tables.rows(m))
        @test d["Alice"] == d["Bob"]
        @test d["Charlie"] == "Cubesat"

        m = match(students="students.csv", projects="projects.csv", force=["Charlie"=>"Airplane"])
        d = Dict(r[:name]=>r[:project] for r in Tables.rows(m))
        @test d["Alice"] == d["Bob"]
        @test d["Charlie"] == "Airplane"
    end

    @testset "Errors" begin
        students = DataFrame()
        @test_throws ArgumentError m = match(students=students, projects="projects.csv")

        projects = DataFrame()
        @test_throws ArgumentError m = match(students="students.csv", projects=projects)

        projects = DataFrame("name"=>[], "min"=>String[], "max"=>String[], "skill:test"=>String[])
        @test_throws Any m = match(students="students.csv", projects=projects)
    end
end


text:
# FrequencyDomainAnalysis.jl
# Tools for spectral density estimation and analysis of phase relationships
# between sets of signals.

# Copyright (C) 2013   Simon Kornblith

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.

# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

export multitaper, psd, xspec, coherence

# Perform tapered FFT of continuous signals
# A is samples x channels x trials
function multitaper(A::Union{AbstractVector{T}, AbstractMatrix{T}, AbstractArray{T,3}},
                             stats:: Vector{S}, fs::Real=1;
                             tapers::Union{Vector{T}, Matrix{T}}=dpss(size(A, 1), 4), nfft::Int=nextpow2(size(A, 1)),
                             freqrange::Range{Int}=0:-1) where {T<:Real, S<:TransformStatistic}

    # We can't guarantee order of kwargs, so this has to be here
    if freqrange == 0:-1
        freqrange = 1:(nfft >> 1 + 1)
    end

    n = size(A, 1)
    nout = nfft >> 1 + 1
    ntapers = size(tapers, 2)
    nchannels = size(A, 2)

    multiplier = convert(T, sqrt(2/fs))
    firstsamplemultiplier = freqrange[1] == 1 ? convert(T, sqrt(1/fs)) : multiplier
    lastsamplemultiplier = freqrange[end] == nout && iseven(nfft) ?
                           convert(T, sqrt(1/fs)) : multiplier

    for stat in stats
        init(stat, length(freqrange), nchannels, ntapers)
    end

    # dtype = outputtype(T)
    dtype = Float64
    fftin = zeros(Complex{dtype}, nfft, size(A, 2))
    fftout = Array{Complex{dtype}}(nout, size(A, 2))
    #fftview = sub(fftout, freqrange, 1:size(A,2))
    fftview = zeros(Complex{dtype}, length(freqrange), size(A, 2))

    # p = FFTW.plan_fft(fftin, 1, FFTW.ESTIMATE, FFTW.NO_TIMELIMIT)

    for j = 1:size(A, 3), i = 1:ntapers
        @inbounds for k = 1:nchannels, l = 1:n
            fftin[l, k] = A[l, k, j]*tapers[l, i]
        end

        # FFTW.execute(dtype, p.plan)
        fftout=fft(fftin, 1)
        # Copy FFT from one buffer to another, including only
        # frequencies in freqrange
        @inbounds for k = 1:size(fftout, 2)
            # Correct edges of FFT so that all power is in the same
            # units. This is necessary because the FFT divides power
            # associated with real signals between positive and
            # negative frequencies, except for power at 0 and FMAX. We
            # don't compute the negative frequencies, but we do need
            # to scale the FFT so that the power is correct.
            fftview[1, k] = fftout[freqrange[1], k]*firstsamplemultiplier
            for l = 2:size(fftview, 1)-1
                fftview[l, k] = fftout[freqrange[l], k]*multiplier
            end
            fftview[end, k] = fftout[freqrange[end], k]*lastsamplemultiplier
        end

        for stat in stats
            accum(stat, fftview, i)
        end
    end

    [finish(stat) for stat in stats]
end

# Calling with a single type or instance
multitaper(A::Union{AbstractVector{T}, AbstractMatrix{T}, AbstractArray{T,3}},
                    stat::S, args...; kw...) where {T<:Real, S <: TransformStatistic} =
    multitaper(A, [stat] , args...; kw...)[1]

# Convenience functions
psd{T<:Real}(A::AbstractArray{T}, args...; kw...) =
    multitaper(A, (PowerSpectrum{T}(),), args...; kw...)[1]
xspec{T<:Real}(A::Vector{T}, B::Vector{T}, args...; kw...) =
    multitaper(hcat(A, B), (CrossSpectrum{T}(),), args...; kw...)[1]
coherence{T<:Real}(A::Vector{T}, B::Vector{T}, args...; kw...) =
    multitaper(hcat(A, B), (Coherence{T}(),), args...; kw...)[1]

# using DSP
# using FFTW
# include("transform_stats.jl")
# include("multitaper.jl")
# x = 0:63
# signal = zeros(length(x), 1, 35)
# for i = 1:size(signal, 3)
# 	signal[:, i] = cos.(0.2pi*x+rand()*2pi)
# end
# c = multitaper([signal signal], Coherence())


text:
using FinancialToolbox, DualNumbers, FinancialFFT, FinancialMonteCarlo

A = 600.0;
N = 16;

method = CarrMadanMethod(A, N);
method_cml = CarrMadanLewisMethod(A, N);
method_l = LewisMethod(A, N);
# S0 = 100.0;
S0 = hyper(100.0, 1.0, 1.0, 0.0);
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.0;
sigma = 0.2;
# sigma = hyper(0.2, 1.0, 1.0, 0.0);
# sigma = dual(0.2, 1.0);
lam = 5.0;
mu1 = 0.03;
sigma1 = 0.02;
zero_rate = ZeroRate(r);

Model = MertonProcess(sigma, lam, mu1, sigma1, Underlying(S0, d));
# Model = NormalInverseGaussianProcess(sigma, 0.01, 0.01, Underlying(S0, d));
#Model = BlackScholesProcess(sigma, Underlying(S0, d));
options = Array{FinancialMonteCarlo.EuropeanOption}(undef, 10)
options[1:3] = [EuropeanOption(T, K) for i = 1:3];
options[8:10] = [EuropeanOption(T, K) for i = 1:3];
options[4:7] = [EuropeanOption(T, 70.0) for i = 4:7];
Nsim = 10000;
Nstep = 30;
mc = MonteCarloConfiguration(Nsim, Nstep);

@time @show pricer(Model, zero_rate, method, options);
@time @show pricer(Model, zero_rate, method_cml, options);
@time @show pricer(Model, zero_rate, method_l, options);
@time @show pricer(Model, zero_rate, mc, options);
typeof(Model) <: BlackScholesProcess ? @time(@show(blsprice(S0, K, r, T, sigma, d))) : nothing;


text:
using Weave
using Test

function weavestring(source; doctype = "pandoc", informat="markdown", mod=Main)
    p1 = Weave.parse_doc(source, informat)
    doc = Weave.WeaveDoc("dummy1.jmd", p1, Dict())
    return Weave.run(doc, doctype=doctype, mod=mod)
end

smod = """
```julia
module TestMod
    x = 3

    function printx()
        print("x")
    end

    export printx
end
```

```julia
using .TestMod
printx()
```
"""

wdoc = weavestring(smod)
@test wdoc.chunks[1].output == "Main.TestMod\n"
@test wdoc.chunks[2].output == "x"

sdoc = weavestring(smod, mod=:sandbox)
@test occursin(r"Main.WeaveSandBox[0-9]*.TestMod\n", sdoc.chunks[1].output)
@test sdoc.chunks[2].output == "x"


text:
using PhysicalOceanography
using Test, Documenter

@testset "PhysicalOceanography.jl" begin

    @test isapprox(f(45),1.03125867182e-4)
    @test isapprox(β(45),1.61687441686e-11)

end

@testset "doctests" begin
    doctest(PhysicalOceanography; manual = false)
end


text:
#--------------------
#--------------------
#### Symbolic
#--------------------
abstract type Symbolic{T} end

TermInterface.exprhead(x::Symbolic) = :call

TermInterface.symtype(x::Number) = typeof(x)
TermInterface.symtype(::Symbolic{T}) where {T} = T

TermInterface.metadata(s::Symbolic) = s.metadata
TermInterface.metadata(s::Symbolic, meta) = Setfield.@set! s.metadata = meta

function hasmetadata(s::Symbolic, ctx)
    metadata(s) isa AbstractDict && haskey(metadata(s), ctx)
end

function getmetadata(s::Symbolic, ctx)
    md = metadata(s)
    if md isa AbstractDict
        md[ctx]
    else
        throw(ArgumentError("$s does not have metadata for $ctx"))
    end
end

function getmetadata(s::Symbolic, ctx, default)
    md = metadata(s)
    md isa AbstractDict ? get(md, ctx, default) : default
end

# pirated for Setfield purposes:
Base.ImmutableDict(d::ImmutableDict{K,V}, x, y)  where {K, V} = ImmutableDict{K,V}(d, x, y)

assocmeta(d::Dict, ctx, val) = (d=copy(d); d[ctx] = val; d)
function assocmeta(d::Base.ImmutableDict, ctx, val)::ImmutableDict{DataType,Any}
    # optimizations
    # If using upto 3 contexts, things stay compact
    if isdefined(d, :parent)
        d.key === ctx && return @set d.value = val
        d1 = d.parent
        if isdefined(d1, :parent)
            d1.key === ctx && return @set d.parent.value = val
            d2 = d1.parent
            if isdefined(d2, :parent)
                d2.key === ctx && return @set d.parent.parent.value = val
            end
        end
    end
    Base.ImmutableDict{DataType, Any}(d, ctx, val)
end

function setmetadata(s::Symbolic, ctx::DataType, val)
    if s.metadata isa AbstractDict
        @set s.metadata = assocmeta(s.metadata, ctx, val)
    else
        # fresh Dict
        @set s.metadata = Base.ImmutableDict{DataType, Any}(ctx, val)
    end
end

Base.isequal(::Symbolic, x) = false
Base.isequal(x, ::Symbolic) = false
Base.isequal(::Symbolic, ::Symbolic) = false

### End of interface

function to_symbolic(x)
    Base.depwarn("`to_symbolic(x)` is deprecated, define the interface for your " *
                 "symbolic structure using `istree(x)`, `operation(x)`, `arguments(x)` " *
                 "and `similarterm(::YourType, f, args, symtype)`", :to_symbolic, force=true)

    x
end

Base.one( s::Symbolic) = one( symtype(s))

Base.zero(s::Symbolic) = zero(symtype(s))

"""
    promote_symtype(f, Ts...)

The result of applying `f` to arguments of [`symtype`](#symtype) `Ts...`

```julia
julia> promote_symtype(+, Real, Real)
Real

julia> promote_symtype(+, Complex, Real)
Number

julia> @syms f(x)::Complex
(f(::Number)::Complex,)

julia> promote_symtype(f, Number)
Complex
```

When constructing [`Term`](#Term)s without an explicit symtype,
`promote_symtype` is used to figure out the symtype of the Term.
"""
promote_symtype(f, Ts...) = Any

const NO_METADATA = nothing

#--------------------
#--------------------
#### Syms
#--------------------
"""
    Sym{T}(name::Symbol)

A named variable of type `T`. Type `T` can be `FnType{X,Y}` which
means the variable is a function with the type signature X -> Y where
`X` is a tuple type of arguments and `Y` is any type.
"""
struct Sym{T, M} <: Symbolic{T}
    name::Symbol
    metadata::M
end

TermInterface.issym(s::Sym) = true
Base.nameof(s::Sym) = s.name

ConstructionBase.constructorof(s::Type{<:Sym{T}}) where {T} = (n,m) -> Sym{T}(n, metadata=m)

function (::Type{Sym{T}})(name; metadata=NO_METADATA) where {T}
    Sym{T, typeof(metadata)}(name, metadata)
end

Base.hash(s::Sym{T}, u::UInt) where {T} = hash(T, hash(s.name, u))

function Base.isequal(a::Sym, b::Sym)
    symtype(a) !== symtype(b) && return false
    isequal(nameof(a), nameof(b))
end

Base.show(io::IO, v::Sym) = Base.show_unquoted(io, v.name)

#---------------------------
#---------------------------
#### Function-like variables
#---------------------------

# Maybe don't even need a new type, can just use Sym{FnType}
struct FnType{X<:Tuple,Y} end

(f::Symbolic{<:FnType})(args...) = Term{promote_symtype(f, symtype.(args)...)}(f, [args...])

function (f::Symbolic)(args...)
    error("Sym $f is not callable. " *
          "Use @syms $f(var1, var2,...) to create it as a callable.")
end

"""
    promote_symtype(f::Sym{FnType{X,Y}}, arg_symtypes...)

The output symtype of applying variable `f` to arugments of symtype `arg_symtypes...`.
if the arguments are of the wrong type then this function will error.
"""
function promote_symtype(f::Symbolic{FnType{X,Y}}, args...) where {X, Y}
    if X === Tuple
        return Y
    end

    # This is to handle `Tuple{T} where T`, so we cannot reliably query the type
    # parameters of the `Tuple` in `FnType`.
    t = Tuple{args...}
    if !(t <: X)
        error("$t is not a subtype of $X.")
    end
    return Y
end

"""
    @syms <lhs_expr>[::T1] <lhs_expr>[::T2]...

For instance:

    @syms foo::Real bar baz(x, y::Real)::Complex

Create one or more variables. `<lhs_expr>` can be just a symbol in which case
it will be the name of the variable, or a function call in which case a function-like
variable which has the same name as the function being called. The Sym type, or
in the case of a function-like Sym, the output type of calling the function
can be set using the `::T` syntax.

# Examples:

- `@syms foo bar::Real baz::Int` will create
variable `foo` of symtype `Number` (the default), `bar` of symtype `Real`
and `baz` of symtype `Int`
- `@syms f(x) g(y::Real, x)::Int h(a::Int, f(b))` creates 1-arg `f` 2-arg `g`
and 2 arg `h`. The second argument to `h` must be a one argument function-like
variable. So, `h(1, g)` will fail and `h(1, f)` will work.
"""
macro syms(xs...)
    defs = map(xs) do x
        n, t = _name_type(x)
        :($(esc(n)) = Sym{$(esc(t))}($(Expr(:quote, n))))
        nt = _name_type(x)
        n, t = nt.name, nt.type
        :($(esc(n)) = Sym{$(esc(t))}($(Expr(:quote, n))))
    end
    Expr(:block, defs...,
         :(tuple($(map(x->esc(_name_type(x).name), xs)...))))
end

function syms_syntax_error()
    error("Incorrect @syms syntax. Try `@syms x::Real y::Complex g(a) f(::Real)::Real` for instance.")
end

function _name_type(x)
    if x isa Symbol
        return (name=x, type=Number)
    elseif x isa Expr && x.head === :(::)
        if length(x.args) == 1
            return (name=nothing, type=x.args[1])
        end
        lhs, rhs = x.args[1:2]
        if lhs isa Expr && lhs.head === :call
            # e.g. f(::Real)::Unreal
            type = map(x->_name_type(x).type, lhs.args[2:end])
            return (name=lhs.args[1], type=:($FnType{Tuple{$(type...)}, $rhs}))
        else
            return (name=lhs, type=rhs)
        end
    elseif x isa Expr && x.head === :ref
        ntype = _name_type(x.args[1]) # a::Number
        N = length(x.args)-1
        return (name=ntype.name,
                type=:(Array{$(ntype.type), $N}),
                array_metadata=:(Base.Slice.(($(x.args[2:end]...),))))
    elseif x isa Expr && x.head === :call
        return _name_type(:($x::Number))
    else
        syms_syntax_error()
    end
end

function Base.show(io::IO, f::Symbolic{<:FnType{X,Y}}) where {X,Y}
    print(io, nameof(f))
    # Use `Base.unwrap_unionall` to handle `Tuple{T} where T`. This is not the
    # best printing, but it's better than erroring.
    argrepr = join(map(t->"::"*string(t), Base.unwrap_unionall(X).parameters), ", ")
    print(io, "(", argrepr, ")")
    print(io, "::", Y)
end

#--------------------
#--------------------
#### Terms
#--------------------
"""
    Term{T}(f, args::AbstractArray)

or
    Term(f, args::AbstractArray)

Symbolic expression representing the result of calling `f(args...)`.

- `operation(t::Term)` returns `f`
- `arguments(t::Term)` returns `args`
- `symtype(t::Term)` returns `T`

If `T` is not provided during construction, it is queried by calling
`SymbolicUtils.promote_symtype(f, map(symtype, args)...)`.

See [promote_symtype](#promote_symtype)
"""
struct Term{T, M} <: Symbolic{T}
    f::Any
    arguments::Any
    metadata::M
    hash::Ref{UInt} # hash cache
end

function ConstructionBase.constructorof(s::Type{<:Term{T}}) where {T}
    function (f, args, meta, hash)
        Term{T, typeof(meta)}(f, args, meta, hash)
    end
end

function (::Type{Term{T}})(f, args; metadata=NO_METADATA) where {T}
    Term{T, typeof(metadata)}(f, args, metadata, Ref{UInt}(0))
end

TermInterface.istree(t::Type{<:Term}) = true

function Term(f, args; metadata=NO_METADATA)
    Term{_promote_symtype(f, args)}(f, args, metadata=metadata)
end

TermInterface.operation(x::Term) = getfield(x, :f)

TermInterface.arguments(x::Term) = getfield(x, :arguments)

function Base.isequal(t1::Term, t2::Term)
    t1 === t2 && return true
    symtype(t1) !== symtype(t2) && return false

    a1 = arguments(t1)
    a2 = arguments(t2)

    isequal(operation(t1), operation(t2)) &&
        length(a1) == length(a2) &&
        all(isequal(l,r) for (l, r) in zip(a1,a2))
end

## This is much faster than hash of an array of Any
hashvec(xs, z) = foldr(hash, xs, init=z)

function Base.hash(t::Term{T}, salt::UInt) where {T}
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    h′ = hashvec(arguments(t), hash(operation(t), hash(T, salt)))
    t.hash[] = h′
    return h′
end

isassociative(::Any) = false
isassociative(::Union{typeof(+),typeof(*)}) = true

_promote_symtype(f::Sym, args) = promote_symtype(f, map(symtype, args)...)
function _promote_symtype(f, args)
    if length(args) == 0
        promote_symtype(f)
    elseif length(args) == 1
        promote_symtype(f, symtype(args[1]))
    elseif length(args) == 2
        promote_symtype(f, symtype(args[1]), symtype(args[2]))
    elseif isassociative(f)
        mapfoldl(symtype, (x,y) -> promote_symtype(f, x, y), args)
    else
        promote_symtype(f, map(symtype, args)...)
    end
end

function term(f, args...; type = nothing)
    if type === nothing
        T = _promote_symtype(f, args)
    else
        T = type
    end
    Term{T}(f, [args...])
end

"""
    similarterm(t, f, args, symtype; metadata=nothing)

Create a term that is similar in type to `t`. Extending this function allows packages
using their own expression types with SymbolicUtils to define how new terms should
be created. Note that `similarterm` may return an object that has a
different type than `t`, because `f` also influences the result.

## Arguments

- `t` the reference term to use to create similar terms
- `f` is the operation of the term
- `args` is the arguments
- The `symtype` of the resulting term. Best effort will be made to set the symtype of the
  resulting similar term to this type.
"""
TermInterface.similarterm(t::Type{<:Symbolic}, f, args; metadata=nothing, exprhead=:call) =
    similarterm(t, f, args, _promote_symtype(f, args); metadata=metadata, exprhead=exprhead)

TermInterface.similarterm(t::Type{<:Term}, f, args, symtype; metadata=nothing, exprhead=:call) =
    Term{_promote_symtype(f, args)}(f, args; metadata=metadata)

#--------------------
#--------------------
####  Pretty printing
#--------------------
const show_simplified = Ref(false)

Base.show(io::IO, t::Term) = show_term(io, t)

isnegative(t::Real) = t < 0
function isnegative(t)
    if istree(t) && operation(t) === (*)
        coeff = first(arguments(t))
        return isnegative(coeff)
    end
    return false
end

setargs(t, args) = Term{symtype(t)}(operation(t), args)
cdrargs(args) = setargs(t, cdr(args))

print_arg(io, x::Union{Complex, Rational}; paren=true) = print(io, "(", x, ")")
isbinop(f) = istree(f) && !istree(operation(f)) && Base.isbinaryoperator(nameof(operation(f)))
function print_arg(io, x; paren=false)
    if paren && isbinop(x)
        print(io, "(", x, ")")
    else
        print(io, x)
    end
end
print_arg(io, s::String; paren=true) = show(io, s)
function print_arg(io, f, x)
    f !== (*) && return print_arg(io, x)
    if Base.isbinaryoperator(nameof(f))
        print_arg(io, x, paren=true)
    else
        print_arg(io, x)
    end
end

function remove_minus(t)
    !istree(t) && return -t
    @assert operation(t) == (*)
    args = arguments(t)
    @assert args[1] < 0
    [-args[1], args[2:end]...]
end

function show_add(io, args)
    negs = filter(isnegative, args)
    nnegs = filter(!isnegative, args)
    for (i, t) in enumerate(nnegs)
        i != 1 && print(io, " + ")
        print_arg(io, +,  t)
    end

    for (i, t) in enumerate(negs)
        if i==1 && isempty(nnegs)
            print_arg(io, -, t)
        else
            print(io, " - ")
            show_mul(io, remove_minus(t))
        end
    end
end

function show_pow(io, args)
    base, ex = args

    if base isa Real && base < 0
        print(io, "(")
        print_arg(io, base)
        print(io, ")")
    else
        print_arg(io, base, paren=true)
    end
    print(io, "^")
    print_arg(io, ex, paren=true)
end

function show_mul(io, args)
    length(args) == 1 && return print_arg(io, *, args[1])

    minus = args[1] isa Number && args[1] == -1
    unit = args[1] isa Number && args[1] == 1

    paren_scalar = (args[1] isa Complex && !_iszero(imag(args[1]))) ||
                   args[1] isa Rational ||
                   (args[1] isa Number && !isfinite(args[1]))

    nostar = minus || unit ||
            (!paren_scalar && args[1] isa Number && !(args[2] isa Number))

    for (i, t) in enumerate(args)
        if i != 1
            if i==2 && nostar
            else
                print(io, "*")
            end
        end
        if i == 1 && minus
            print(io, "-")
        elseif i == 1 && unit
        else
            print_arg(io, *, t)
        end
    end
end

function show_ref(io, f, args)
    x = args[1]
    idx = args[2:end]

    istree(x) && print(io, "(")
    print(io, x)
    istree(x) && print(io, ")")
    print(io, "[")
    for i=1:length(idx)
        print_arg(io, idx[i])
        i != length(idx) && print(io, ", ")
    end
    print(io, "]")
end

function show_call(io, f, args)
    fname = istree(f) ? Symbol(repr(f)) : nameof(f)
    binary = Base.isbinaryoperator(fname)
    if binary
        for (i, t) in enumerate(args)
            i != 1 && print(io, " $fname ")
            print_arg(io, t, paren=true)
        end
    else
        if f isa Sym
            Base.show_unquoted(io, nameof(f))
        else
            Base.show(io, f)
        end
        print(io, "(")
        for i=1:length(args)
            print(io, args[i])
            i != length(args) && print(io, ", ")
        end
        print(io, ")")
    end
end

function show_term(io::IO, t)
    if get(io, :simplify, show_simplified[])
        return print(IOContext(io, :simplify=>false), simplify(t))
    end

    f = operation(t)
    args = arguments(t)

    if f === (+)
        show_add(io, args)
    elseif f === (*)
        show_mul(io, args)
    elseif f === (^)
        show_pow(io, args)
    elseif f === (getindex)
        show_ref(io, f, args)
    else
        show_call(io, f, args)
    end

    return nothing
end

showraw(io, t) = Base.show(IOContext(io, :simplify=>false), t)
showraw(t) = showraw(stdout, t)

######   Add Mul and Pow

sdict(kv...) = Dict{Any, Number}(kv...)

const SN = Symbolic{<:Number}
"""
    Add(T, coeff, dict::Dict)

Represents `coeff + (key1 * val1) + (key2 * val2) + ...`

where keys and values come from the dictionary (`dict`).
where `coeff` and the vals are `<:Number` and keys are symbolic.

- `operation(::Add)` -- returns `+`.
- `symtype(::Add)` -- returns `T`.
- `arguments(::Add)` -- returns a totally ordered vector of arguments. i.e.
  `[coeff, keyM*valM, keyN*valN...]`
"""
struct Add{X<:Number, T<:Number, D, M} <: Symbolic{X}
    coeff::T
    dict::D
    sorted_args_cache::Ref{Any}
    hash::Ref{UInt}
    metadata::M
end

function Add(T, coeff, dict; metadata=NO_METADATA)
    if isempty(dict)
        return coeff
    elseif _iszero(coeff) && length(dict) == 1
        k,v = first(dict)
        return _isone(v) ? k : Mul(T, makemul(v, k)...)
    end

    Add{T, typeof(coeff), typeof(dict), typeof(metadata)}(coeff, dict, Ref{Any}(nothing), Ref{UInt}(0), metadata)
end

TermInterface.symtype(a::Add{X}) where {X} = X

TermInterface.istree(a::Type{Add}) = true

TermInterface.operation(a::Add) = +

function TermInterface.unsorted_arguments(a::Add)
    args = [v*k for (k,v) in a.dict]
    iszero(a.coeff) ? args : vcat(a.coeff, args)
end

function TermInterface.arguments(a::Add)
    a.sorted_args_cache[] !== nothing && return a.sorted_args_cache[]
    args = sort!([v*k for (k,v) in a.dict], lt=<ₑ)
    a.sorted_args_cache[] = iszero(a.coeff) ? args : vcat(a.coeff, args)
end

Base.isequal(a::Add, b::Add) = a.coeff == b.coeff && isequal(a.dict, b.dict)

Base.show(io::IO, a::Add) = show_term(io, a)

"""
    makeadd(sign, coeff::Number, xs...)

Any Muls inside an Add should always have a coeff of 1
and the key (in Add) should instead be used to store the actual coefficient
"""
function makeadd(sign, coeff, xs...)
    d = sdict()
    for x in xs
        if x isa Add
            coeff += x.coeff
            _merge!(+, d, x.dict, filter=_iszero)
            continue
        end
        if x isa Number
            coeff += x
            continue
        end
        if x isa Mul
            k = Mul(symtype(x), 1, x.dict)
            v = sign * x.coeff + get(d, k, 0)
        else
            k = x
            v = sign + get(d, x, 0)
        end
        if iszero(v)
            delete!(d, k)
        else
            d[k] = v
        end
    end
    coeff, d
end

add_t(a,b) = promote_symtype(+, symtype(a), symtype(b))
sub_t(a,b) = promote_symtype(-, symtype(a), symtype(b))
sub_t(a) = promote_symtype(-, symtype(a))

function +(a::SN, b::SN)
    if a isa Add
        coeff, dict = makeadd(1, 0, b)
        T = promote_symtype(+, symtype(a), symtype(b))
        return Add(add_t(a,b), a.coeff + coeff, _merge(+, a.dict, dict, filter=_iszero))
    elseif b isa Add
        return b + a
    end
    Add(add_t(a,b), makeadd(1, 0, a, b)...)
end

+(a::Number, b::SN) = Add(add_t(a,b), makeadd(1, a, b)...)

+(a::SN, b::Number) = Add(add_t(a,b), makeadd(1, b, a)...)

+(a::SN) = a

+(a::Add, b::Add) = Add(add_t(a,b),
                        a.coeff + b.coeff,
                        _merge(+, a.dict, b.dict, filter=_iszero))

+(a::Number, b::Add) = iszero(a) ? b : Add(add_t(a,b), a + b.coeff, b.dict)

+(b::Add, a::Number) = iszero(a) ? b : Add(add_t(a,b), a + b.coeff, b.dict)

-(a::Add) = Add(sub_t(a), -a.coeff, mapvalues((_,v) -> -v, a.dict))

-(a::SN) = Add(sub_t(a), makeadd(-1, 0, a)...)

-(a::Add, b::Add) = Add(sub_t(a,b),
                        a.coeff - b.coeff,
                        _merge(-, a.dict, b.dict, filter=_iszero))

-(a::SN, b::SN) = a + (-b)

-(a::Number, b::SN) = a + (-b)

-(a::SN, b::Number) = a + (-b)

"""
    Mul(T, coeff, dict)

Represents coeff * (key1 ^ val1) * (key2 ^ val2) * ....

where coeff is a <:Number and keys and values come from the dictionary (`dict`).
where `coeff` and the vals are `<:Number` and keys are symbolic.

- `symtype(::Mul)` -- returns `T`.
- `operation(::Mul)` -- returns `*`.
- `arguments(::Mul)` -- returns a totally ordered vector of arguments. i.e.
  `[coeff, keyM^valM, keyN^valN...]`
"""
struct Mul{X<:Number, T<:Number, D, M} <: Symbolic{X}
    coeff::T
    dict::D
    sorted_args_cache::Ref{Any}
    hash::Ref{UInt}
    metadata::M
end

for S in [Add, Mul]
    @eval function ConstructionBase.constructorof(s::Type{<:$S{T}}) where {T}
        function (coeff, dict, argscache, hash, m)
            $S{T,
                typeof(coeff),
                typeof(dict),
                typeof(m)}(coeff,
            dict,
            argscache,
            hash,
            m)
        end
    end
end

function Mul(T, a,b; metadata=NO_METADATA)
    isempty(b) && return a
    if _isone(a) && length(b) == 1
        pair = first(b)
        if _isone(last(pair)) # first value
            return first(pair)
        else
            return unstable_pow(first(pair), last(pair))
        end
    else
        Mul{T, typeof(a), typeof(b), typeof(metadata)}(a,b, Ref{Any}(nothing), Ref{UInt}(0), metadata)
    end
end

TermInterface.symtype(a::Mul{X}) where {X} = X

TermInterface.istree(a::Type{Mul}) = true

TermInterface.operation(a::Mul) = *

unstable_pow(a, b) = a isa Integer && b isa Integer ? (a//1) ^ b : a ^ b

function TermInterface.unsorted_arguments(a::Mul)
    args = [unstable_pow(k, v) for (k,v) in a.dict]
    isone(a.coeff) ? args : vcat(a.coeff, args)
end

function TermInterface.arguments(a::Mul)
    a.sorted_args_cache[] !== nothing && return a.sorted_args_cache[]
    args = sort!([unstable_pow(k, v) for (k,v) in a.dict], lt=<ₑ)
    a.sorted_args_cache[] = isone(a.coeff) ? args : vcat(a.coeff, args)
end

Base.isequal(a::Mul, b::Mul) = a.coeff == b.coeff && isequal(a.dict, b.dict)

Base.show(io::IO, a::Mul) = show_term(io, a)

function makemul(coeff, xs...; d=sdict())
    for x in xs
        if x isa Pow && x.exp isa Number
            d[x.base] = x.exp + get(d, x.base, 0)
        elseif x isa Number
            coeff *= x
        elseif x isa Mul
            coeff *= x.coeff
            _merge!(+, d, x.dict, filter=_iszero)
        else
            v = 1 + get(d, x, 0)
            if _iszero(v)
                delete!(d, x)
            else
                d[x] = v
            end
        end
    end
    (coeff, d)
end

mul_t(a,b) = promote_symtype(*, symtype(a), symtype(b))
mul_t(a) = promote_symtype(*, symtype(a))

*(a::SN) = a

function *(a::SN, b::SN)
    # Always make sure Div wraps Mul
    if a isa Div && b isa Div
        Div(a.num * b.num, a.den * b.den)
    elseif a isa Div
        Div(a.num * b, a.den)
    elseif b isa Div
        Div(a * b.num, b.den)
    else
        Mul(mul_t(a,b), makemul(1, a, b)...)
    end
end

*(a::Mul, b::Mul) = Mul(mul_t(a, b),
                        a.coeff * b.coeff,
                        _merge(+, a.dict, b.dict, filter=_iszero))

function *(a::Number, b::SN)
    if iszero(a)
        a
    elseif isone(a)
        b
    elseif b isa Div
        Div(a*b.num, b.den)
    elseif b isa Add
        # 2(a+b) -> 2a + 2b
        T = promote_symtype(+, typeof(a), symtype(b))
        Add(T, b.coeff * a, Dict(k=>v*a for (k, v) in b.dict))
    else
        Mul(mul_t(a, b), makemul(a, b)...)
    end
end

*(a::SN, b::Number) = b * a

\(a::SN, b::Union{Number, SN}) = b / a

\(a::Number, b::SN) = b / a

/(a::SN, b::Number) = (b isa Integer ? 1//b : inv(b)) * a

//(a::Union{SN, Number}, b::SN) = a / b

//(a::SN, b::T) where {T <: Number} = (one(T) // b) * a

"""
    Div(numerator_factors, denominator_factors, simplified=false)

"""
struct Div{T,N,D, M} <: Symbolic{T}
    num::N
    den::D
    simplified::Bool
    metadata::M
end

Base.hash(x::Div, u::UInt64) = hash(x.num, hash(x.den, u))
Base.isequal(x::Div, y::Div) = isequal(x.num, y.num) && isequal(x.den, y.den)

const Rat = Union{Rational, Integer}

ratcoeff(x) = false, NaN
ratcoeff(x::Rat) = true, x
ratcoeff(x::Mul) = ratcoeff(x.coeff)
ratio(x::Integer,y::Integer) = iszero(rem(x,y)) ? div(x,y) : x//y
ratio(x::Rat,y::Rat) = x//y
function maybe_intcoeff(x::Mul)
    x.coeff isa Rational && isone(x.coeff.den) ? Setfield.@set!(x.coeff = x.coeff.num) : x
end
maybe_intcoeff(x::Rational) = isone(x.den) ? x.num : x
maybe_intcoeff(x) = x

function (::Type{Div{T}})(n, d, simplified=false; metadata=nothing) where {T}
    _iszero(n) && return zero(typeof(n))
    _isone(d) && return n

    if n isa Div && d isa Div
        return Div{T}(n.num * d.den, n.den * d.num)
    elseif n isa Div
        return Div{T}(n.num, n.den * d)
    elseif d isa Div
        return Div{T}(n * d.den, d.num)
    end

    d isa Number && _isone(-d) && return -1 * n
    n isa Rat && d isa Rat && return n // d # maybe called by oblivious code in simplify

    # GCD coefficient upon construction
    rat, nc = ratcoeff(n)
    if rat
        rat, dc = ratcoeff(d)
        if rat
            g = gcd(nc, dc) * sign(dc) # make denominator positive
            invdc = ratio(1, g)
            n = maybe_intcoeff(invdc * n)
            d = maybe_intcoeff(invdc * d)
        end
    end

    Div{T, typeof(n), typeof(d), typeof(metadata)}(n, d, simplified, metadata)
end

function Div(n,d, simplified=false; kw...)
    Div{promote_symtype((/), symtype(n), symtype(d))}(n,d, simplified; kw...)
end

numerators(x) = istree(x) && operation(x) == (*) ? arguments(x) : [x]
numerators(d::Div) = numerators(d.num)

denominators(x) = [1]
denominators(d::Div) = numerators(d.den)

TermInterface.istree(d::Type{Div}) = true

TermInterface.operation(d::Div) = (/)

function TermInterface.arguments(d::Div)
    [d.num, d.den]
end

Base.show(io::IO, d::Div) = show_term(io, d)

/(a::Union{SN,Number}, b::SN) = Div(a,b)

"""
    Pow(base, exp)

Represents `base^exp`, a lighter version of `Mul(1, Dict(base=>exp))`
"""
struct Pow{X, B, E, M} <: Symbolic{X}
    base::B
    exp::E
    metadata::M
end

function ConstructionBase.constructorof(::Type{<:Pow{X}}) where {X}
    (base, exp, m) ->
    Pow{promote_symtype(^, symtype(base), symtype(exp)), typeof(base), typeof(exp), typeof(m)}(base,exp,m)
end

function (::Type{<:Pow{T}})(a, b; metadata=NO_METADATA) where {T}
    _iszero(b) && return 1
    _isone(b) && return a
    Pow{T, typeof(a), typeof(b), typeof(metadata)}(a,b,metadata)
end
function Pow(a, b; metadata=NO_METADATA)
    Pow{promote_symtype(^, symtype(a), symtype(b))}(makepow(a, b)..., metadata=metadata)
end
TermInterface.symtype(a::Pow{X}) where {X} = X

TermInterface.istree(a::Type{Pow}) = true

TermInterface.operation(a::Pow) = ^

# Use `Union` to avoid promoting the base and exponent to the same type.
# For instance, if `a.base` is a multivariate polynomial and  `a.exp` is a number,
# we don't want to promote `a.exp` to a multivariate polynomial.
TermInterface.arguments(a::Pow) = Union{typeof(a.base), typeof(a.exp)}[a.base, a.exp]

Base.hash(p::Pow, u::UInt) = hash(p.exp, hash(p.base, u))

Base.isequal(p::Pow, b::Pow) = isequal(p.base, b.base) && isequal(p.exp, b.exp)

Base.show(io::IO, p::Pow) = show_term(io, p)

function makepow(a, b)
    base = a
    exp = b
    if a isa Pow
        base = a.base
        exp = a.exp * b
    end
    return (base, exp)
end

^(a::SN, b) = Pow(a, b)

^(a::SN, b::SN) = Pow(a, b)

^(a::Number, b::SN) = Pow(a, b)

function ^(a::Mul, b::Number)
    coeff = unstable_pow(a.coeff, b)
    Mul(promote_symtype(^, symtype(a), symtype(b)),
        coeff, mapvalues((k, v) -> b*v, a.dict))
end

function *(a::Mul, b::Pow)
    if b.exp isa Number
        Mul(mul_t(a, b),
            a.coeff, _merge(+, a.dict, Base.ImmutableDict(b.base=>b.exp), filter=_iszero))
    else
        Mul(mul_t(a, b),
            a.coeff, _merge(+, a.dict, Base.ImmutableDict(b=>1), filter=_iszero))
    end
end

*(a::Pow, b::Mul) = b * a

function copy_similar(d, others)
    K = promote_type(keytype(d), keytype.(others)...)
    V = promote_type(valtype(d), valtype.(others)...)
    Dict{K, V}(d)
end

_merge(f, d, others...; filter=x->false) = _merge!(f, copy_similar(d, others), others...; filter=filter)
function _merge!(f, d, others...; filter=x->false)
    acc = d
    for other in others
        for (k, v) in other
            v = f(v)
            if haskey(acc, k)
                v = acc[k] + v
            end
            if filter(v)
                delete!(acc, k)
            else
                acc[k] = v
            end
        end
    end
    acc
end

function mapvalues(f, d1::AbstractDict)
    d = copy(d1)
    for (k, v) in d
        d[k] = f(k, v)
    end
    d
end

const NumericTerm = Union{Term{<:Number}, Mul{<:Number},
                          Add{<:Number}, Pow{<:Number}, Div{<:Number}}

function TermInterface.similarterm(t::Type{P}, f, args, symtype; metadata=nothing, exprhead=:call) where P<:NumericTerm
    T = symtype
    if T === nothing
        T = _promote_symtype(f, args)
    end
    if f === (+)
        Add(T, makeadd(1, 0, args...)...; metadata=metadata)
    elseif f == (*)
        Mul(T, makemul(1, args...)...; metadata=metadata)
    elseif f == (/)
        @assert length(args) == 2
        Div{T}(args...; metadata=metadata)
    elseif f == (^) && length(args) == 2
        Pow{T}(makepow(args...)...; metadata=metadata)
    else
        Term{T}(f, args; metadata=metadata)
    end
end

function Base.hash(t::Union{Add,Mul}, u::UInt64)
    !iszero(u) && return hash(hash(t, zero(UInt64)), u)
    h = t.hash[]
    !iszero(h) && return h
    hashoffset = t isa Add ? 0xaddaddaddaddadda : 0xaaaaaaaaaaaaaaaa
    h′= hash(hashoffset, hash(t.coeff, hash(t.dict, u)))
    t.hash[] = h′
    return h′
end

import AbstractTrees

struct TreePrint
    op
    x
end
AbstractTrees.children(x::Term) = arguments(x)
function AbstractTrees.children(x::Union{Add, Mul})
    children = Any[x.coeff]
    for (key, coeff) in pairs(x.dict)
        if coeff == 1
            push!(children, key)
        else
            push!(children, TreePrint(x isa Add ? (:*) : (:^), (key, coeff)))
        end
    end
    return children
end
AbstractTrees.children(x::Union{Pow}) = [x.base, x.exp]
AbstractTrees.children(x::TreePrint) = [x.x[1], x.x[2]]

print_tree(x; show_type=false, maxdepth=Inf, kw...) = print_tree(stdout, x; show_type=show_type, maxdepth=maxdepth, kw...)
function print_tree(_io::IO, x::Union{Term, Add, Mul, Pow, Div}; show_type=false, kw...)
    AbstractTrees.print_tree(_io, x; withinds=true, kw...) do io, y, inds
        if istree(y)
            print(io, operation(y))
        elseif y isa TreePrint
            print(io, "(", y.op, ")")
        else
            print(io, y)
        end
        if !(y isa TreePrint) && show_type
            print(io, " [", typeof(y), "]")
        end
    end
end

TermInterface.istree(t::Type{<:Sym}) = false
TermInterface.istree(t::Type{<:Symbolic}) = true


text:
# # Splash on a sphere using Oceananigans.HydrostaticFreeSurfaceModel
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, GLMakie"
# ```

# ## A spherical domain
#
# We two-dimensional latitude-longitude grid on a sphere of unit radius,

using Oceananigans
using Oceananigans.Grids: RegularLatitudeLongitudeGrid

grid = RegularLatitudeLongitudeGrid(size = (360, 60, 1),
                                    radius = 1,
                                    latitude = (-60, 60),
                                    longitude = (-180, 180),
                                    z = (-1, 0))

using Oceananigans.Coriolis: HydrostaticSphericalCoriolis

coriolis = HydrostaticSphericalCoriolis(rotation_rate=10.0)

# ## Building a `HydrostaticFreeSurfaceModel`
#
# We use `grid` and `coriolis` to build a simple `HydrostaticFreeSurfaceModel`,

using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel
using Oceananigans.Models.HydrostaticFreeSurfaceModels: VectorInvariant, ExplicitFreeSurface

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    momentum_advection = VectorInvariant(),
                                    free_surface = ExplicitFreeSurface(gravitational_acceleration=1),
                                    tracers = (),
                                    buoyancy = nothing,
                                    coriolis = coriolis,
                                    closure = nothing)

# ## Making a splash
#
# We make a splash by imposing a Gaussian height field,

Gaussian(λ, ϕ, L) = exp(-(λ^2 + ϕ^2) / 2L^2)

a = 0.01 # splash amplitude
L = 10 # degree
ϕ₀ = 5 # degrees

splash(λ, ϕ, z) = a * Gaussian(λ, ϕ - ϕ₀, L)

set!(model, η=splash)

# ## Running a `Simulation`
#
# We pick a time-step that resolves the surface dynamics,

using Oceananigans.Utils: prettytime

g = model.free_surface.gravitational_acceleration

gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

# Time-scale for gravity wave propagation across the smallest grid cell
wave_propagation_time_scale = min(grid.radius * cosd(maximum(abs, grid.ϕᵃᶜᵃ)) * deg2rad(grid.Δλ),
                                  grid.radius * deg2rad(grid.Δϕ)) / gravity_wave_speed

simulation = Simulation(model,
                        Δt = 0.05wave_propagation_time_scale,
                        stop_iteration = 4000,
                        iteration_interval = 100,
                        progress = s -> @info "Iteration = $(s.model.clock.iteration) / $(s.stop_iteration)")

# ## Output
#
# We output the velocity field and free surface displacement,

output_fields = merge(model.velocities, (η=model.free_surface.η,))

using Oceananigans.OutputWriters: JLD2OutputWriter, TimeInterval, IterationInterval

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = IterationInterval(10),
                                                      prefix = "spherical_splash",
                                                      force = true)

run!(simulation)

# ## Visualizing the results

using JLD2, Printf, Oceananigans.Grids, GLMakie
using Oceananigans.Utils: hours

λ, ϕ, r = nodes(model.free_surface.η, reshape=true)

λ = λ .+ 180  # Convert to λ ∈ [0°, 360°]
ϕ = 90 .- ϕ   # Convert to ϕ ∈ [0°, 180°] (0° at north pole)

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

iter = Node(0)
plot_title = @lift @sprintf("Oceananigans.jl on the sphere! Spherical splash: u, v, η @ time = %s",
                            prettytime(file["timeseries/t/" * string($iter)]))

u = @lift file["timeseries/u/" * string($iter)][:, :, 1]
v = @lift file["timeseries/v/" * string($iter)][:, :, 1]
η = @lift file["timeseries/η/" * string($iter)][:, :, 1]

# Plot on the unit sphere to align with the spherical wireframe.
x3 = @. cosd(λ) * sind(ϕ)
y3 = @. sind(λ) * sind(ϕ)
z3 = @. cosd(ϕ) * λ ./ λ

x = @lift (1 .+ 20 .* file["timeseries/η/" * string($iter)][:, :, 1]) .* x3[:, :, 1]
y = @lift (1 .+ 20 .* file["timeseries/η/" * string($iter)][:, :, 1]) .* y3[:, :, 1]
z = @lift (1 .+ 20 .* file["timeseries/η/" * string($iter)][:, :, 1]) .* z3[:, :, 1]

fig = Figure(resolution = (1920, 1080))

clims = [(-0.003, 0.003), (-0.003, 0.003), (-0.01, 0.01)]

for (n, var) in enumerate((u, v, η))
    ax = fig[1, n] = LScene(fig, title="$n")
    wireframe!(ax, Sphere(Point3f0(0), 0.98f0), show_axis=false)
    surface!(ax, x, y, z, color=var, colormap=:balance, colorrange=clims[n])
    rotate_cam!(ax.scene, (2π/3, 0, 0))
    zoom!(ax.scene, (0, 0, 0), 2, false)
end

supertitle = fig[0, :] = Label(fig, plot_title, textsize=30)

record(fig, "spherical_splash.mp4", iterations, framerate=30) do i
    @info "Animating iteration $i/$(iterations[end])..."
    iter[] = i
end