text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## boolean conversions ##

convert(::Type{Bool}, x::Bool) = x
convert(::Type{Bool}, x::Float16) = x==0 ? false : x==1 ? true : throw(InexactError())
convert(::Type{Bool}, x::Real) = x==0 ? false : x==1 ? true : throw(InexactError())

# promote Bool to any other numeric type
promote_rule{T<:Number}(::Type{Bool}, ::Type{T}) = T

typemin(::Type{Bool}) = false
typemax(::Type{Bool}) = true

## boolean operations ##

function !(x::Bool)
    ## We need a better heuristic to detect this automatically
    @_pure_meta
    return box(Bool,not_int(unbox(Bool,x)))
end

(~)(x::Bool) = !x
(&)(x::Bool, y::Bool) = box(Bool,and_int(unbox(Bool,x),unbox(Bool,y)))
(|)(x::Bool, y::Bool) = box(Bool,or_int(unbox(Bool,x),unbox(Bool,y)))
xor(x::Bool, y::Bool) = (x!=y)

>>(x::Bool, c::Unsigned) = Int(x) >> c
<<(x::Bool, c::Unsigned) = Int(x) << c
>>>(x::Bool, c::Unsigned) = Int(x) >>> c

>>(x::Bool, c::Int) = Int(x) >> c
<<(x::Bool, c::Int) = Int(x) << c
>>>(x::Bool, c::Int) = Int(x) >>> c

>>(x::Bool, c::Integer) = Int(x) >> c
<<(x::Bool, c::Integer) = Int(x) << c
>>>(x::Bool, c::Integer) = Int(x) >>> c

signbit(x::Bool) = false
sign(x::Bool) = x
abs(x::Bool) = x
abs2(x::Bool) = x

<(x::Bool, y::Bool) = y&!x
<=(x::Bool, y::Bool) = y|!x

## do arithmetic as Int ##

+(x::Bool) =  Int(x)
-(x::Bool) = -Int(x)

+(x::Bool, y::Bool) = Int(x) + Int(y)
-(x::Bool, y::Bool) = Int(x) - Int(y)
*(x::Bool, y::Bool) = x & y
^(x::Bool, y::Bool) = x | !y
^(x::Integer, y::Bool) = ifelse(y, x, one(x))

function +{T<:AbstractFloat}(x::Bool, y::T)::promote_type(Bool,T)
    return ifelse(x, one(y) + y, y)
end
+(y::AbstractFloat, x::Bool) = x + y

function *{T<:Number}(x::Bool, y::T)::promote_type(Bool,T)
    return ifelse(x, y, copysign(zero(y), y))
end
function *{T<:Unsigned}(x::Bool, y::T)::promote_type(Bool,T)
    return ifelse(x, y, zero(y))
end
*(y::Number, x::Bool) = x * y

div(x::Bool, y::Bool) = y ? x : throw(DivideError())
fld(x::Bool, y::Bool) = div(x,y)
cld(x::Bool, y::Bool) = div(x,y)
rem(x::Bool, y::Bool) = y ? false : throw(DivideError())
mod(x::Bool, y::Bool) = rem(x,y)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# We have rotations along one, two or three axes (e.g. RotX, RotXY and RotXYZ).
# They compose together nicely, so the user can make Euler angles by:
#
#     RotX(θx) * RotY(θy) * RotZ(θz) -> RotXYZ(θx, θy, θz)
#
# and never get confused by the order of application, etc.

#########################
# Single axis rotations #
#########################
"""
    immutable RotX{T} <: Rotation{T}
    RotX(theta)

A 3×3 rotation matrix which represents a rotation by `theta` about the X axis.
"""
immutable RotX{T} <: Rotation{T}
    theta::T
end

@inline convert{R<:RotX}(::Type{R}, r::RotX) = R(r.theta)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{R}){R<:RotX}(t::NTuple{9}) = error("Cannot construct a cardinal axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotX{T}, i::Integer)
    if i == 1
        one(T)
    elseif i < 5
        zero(T)
    elseif i == 5
        cos(r.theta)
    elseif i == 6
        sin(r.theta)
    elseif i == 7
        zero(T)
    elseif i == 8
        -sin(r.theta)
    elseif i == 9
        cos(r.theta)
    else
        throw(BoundsError(r,i))
    end
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotX{T})
    (one(T),  zero(T),      zero(T),   # transposed representation
     zero(T), cos(r.theta), sin(r.theta),
     zero(T), -sin(r.theta), cos(r.theta))
end

@inline function Base.:*(r::RotX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    ct, st = cos(r.theta), sin(r.theta)
    T = promote_type(eltype(r), eltype(v))
    return similar_type(v,T)(v[1],
                             v[2] * ct - v[3] * st,
                             v[3] * ct + v[2] * st)
end

@inline Base.:*(r1::RotX, r2::RotX) = RotX(r1.theta + r2.theta)

@inline inv(r::RotX) = RotX(-r.theta)

# define null rotations for convenience
@inline eye(::Type{RotX}) = RotX(0.0)
@inline eye{T}(::Type{RotX{T}}) = RotX{T}(zero(T))

"""
    immutable RotY{T} <: Rotation{T}
    RotY(theta)

A 3×3 rotation matrix which represents a rotation by `theta` about the Y axis.
"""
immutable RotY{T} <: Rotation{T}
    theta::T
end

@inline convert{R<:RotY}(::Type{R}, r::RotY) = R(r.theta)

@inline (::Type{R}){R<:RotY}(t::NTuple{9}) = error("Cannot construct a cardinal axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotY{T}, i::Integer)
    if i == 1
        cos(r.theta)
    elseif i == 2
        zero(T)
    elseif i == 3
        -sin(r.theta)
    elseif i == 4
        zero(T)
    elseif i == 5
        one(T)
    elseif i == 6
        zero(T)
    elseif i == 7
        sin(r.theta)
    elseif i == 8
        zero(T)
    elseif i == 9
        cos(r.theta)
    else
        throw(BoundsError(r,i))
    end
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotY{T})
    (cos(r.theta), zero(T), -sin(r.theta),   # transposed representation
     zero(T),      one(T),   zero(T),
     sin(r.theta), zero(T),  cos(r.theta))
end

@inline function Base.:*(r::RotY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    ct, st = cos(r.theta), sin(r.theta)
    T = promote_type(eltype(r), eltype(v))
    return similar_type(v,T)(v[1] * ct + v[3] * st,
                             v[2],
                             v[3] * ct - v[1] * st)
end

@inline Base.:*(r1::RotY, r2::RotY) = RotY(r1.theta + r2.theta)

@inline inv(r::RotY) = RotY(-r.theta)

# define null rotations for convenience
@inline eye(::Type{RotY}) = RotY(0.0)
@inline eye{T}(::Type{RotY{T}}) = RotY{T}(zero(T))

"""
    immutable RotZ{T} <: Rotation{T}
    RotZ(theta)

A 3×3 rotation matrix which represents a rotation by `theta` about the Z axis.
"""
immutable RotZ{T} <: Rotation{T}
    theta::T
end

@inline convert{R<:RotZ}(::Type{R}, r::RotZ) = R(r.theta)

@inline (::Type{R}){R<:RotZ}(t::NTuple{9}) = error("Cannot construct a cardinal axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotZ{T}, i::Integer)
    if i == 1
        cos(r.theta)
    elseif i == 2
        sin(r.theta)
    elseif i == 3
        zero(T)
    elseif i == 4
        -sin(r.theta)
    elseif i == 5
        cos(r.theta)
    elseif i < 9
        zero(T)
    elseif i == 9
        one(T)
    else
        throw(BoundsError(r,i))
    end
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotZ{T})
    ( cos(r.theta), sin(r.theta), zero(T),   # transposed representation
     -sin(r.theta), cos(r.theta), zero(T),
      zero(T),      zero(T),      one(T))
end

@inline function Base.:*(r::RotZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    ct, st = cos(r.theta), sin(r.theta)
    T = promote_type(eltype(r), eltype(v))
    return similar_type(v,T)(v[1] * ct - v[2] * st,
                             v[2] * ct + v[1] * st,
                             v[3])
end

@inline Base.:*(r1::RotZ, r2::RotZ) = RotZ(r1.theta + r2.theta)

# define null rotations for convenience
@inline eye(::Type{RotZ}) = RotZ(0.0)
@inline eye{T}(::Type{RotZ{T}}) = RotZ{T}(zero(T))

@inline inv(r::RotZ) = RotZ(-r.theta)

function Base.rand{R <: Union{RotX,RotY,RotZ}}(::Type{R})
    T = eltype(R)
    if T == Any
        T = Float64
    end

    return R(2*pi*rand(T))
end

################################################################################
################################################################################

######################
# Two axis rotations #
######################

"""
    immutable RotXY{T} <: Rotation{T}
    RotXY(theta_x, theta_y)

A 3×3 rotation matrix which represents a rotation by `theta_y` about the Y axis,
followed by a rotation by `theta_x` about the X axis.
"""
immutable RotXY{T} <: Rotation{T}
    theta1::T
    theta2::T
end

@inline (::Type{R}){R<:RotXY}(r::RotXY) = R(r.theta1, r.theta2)
@inline convert{R<:RotXY}(::Type{R}, r::RotXY) = R(r.theta1, r.theta2)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 2 inputs (not 9).
@inline (::Type{RotXY}){X,Y}(x::X, y::Y) = RotXY{promote_type(X, Y)}(x, y)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{R}){R<:RotXY}(t::NTuple{9}) = error("Cannot construct a two-axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotXY{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotXY{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    # transposed representation
    (cosθ₂,  sinθ₁*sinθ₂,  cosθ₁*-sinθ₂,
     zero(T),      cosθ₁,          sinθ₁,
     sinθ₂,  -sinθ₁*cosθ₂,   cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotXY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    return similar_type(v,T)(cosθ₂*v[1] + sinθ₂*v[3],
                             sinθ₁*sinθ₂*v[1] + cosθ₁*v[2] + -sinθ₁*cosθ₂*v[3],
                             cosθ₁*-sinθ₂*v[1] + sinθ₁*v[2] + cosθ₁*cosθ₂*v[3])
end

@inline Base.:*(r1::RotX, r2::RotY) = RotXY(r1.theta, r2.theta)
@inline Base.:*(r1::RotXY, r2::RotY) = RotXY(r1.theta1, r1.theta2 + r2.theta)
@inline Base.:*(r1::RotX, r2::RotXY) = RotXY(r1.theta + r2.theta1, r2.theta2)

@inline inv(r::RotXY) = RotYX(-r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotXY}) = RotXY(0.0, 0.0)
@inline eye{T}(::Type{RotXY{T}}) = RotXY{T}(zero(T), zero(T))

"""
    immutable RotYX{T} <: Rotation{T}
    RotYX(theta_y, theta_x)

A 3×3 rotation matrix which represents a rotation by `theta_x` about the X axis,
followed by a rotation by `theta_y` about the Y axis.
"""
immutable RotYX{T} <: Rotation{T}
    theta1::T
    theta2::T
end

@inline (::Type{R}){R<:RotYX}(r::RotYX) = R(r.theta1, r.theta2)
@inline convert{R<:RotYX}(::Type{R}, r::RotYX) = R(r.theta1, r.theta2)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 2 inputs (not 9).
@inline (::Type{RotYX}){Y,X}(y::Y, x::X) = RotYX{promote_type(Y, X)}(y, x)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{R}){R<:RotYX}(t::NTuple{9}) = error("Cannot construct a two-axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotYX{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotYX{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    # transposed representation
    (cosθ₁,        zero(T),       -sinθ₁,
     sinθ₁*sinθ₂,  cosθ₂,   cosθ₁*sinθ₂,
     sinθ₁*cosθ₂,  -sinθ₂,  cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotYX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    return similar_type(v,T)(cosθ₁*v[1] + sinθ₁*sinθ₂*v[2] + sinθ₁*cosθ₂*v[3],
                             cosθ₂*v[2] + -sinθ₂*v[3],
                             -sinθ₁*v[1] + cosθ₁*sinθ₂*v[2] + cosθ₁*cosθ₂*v[3])
end

@inline Base.:*(r1::RotY, r2::RotX) = RotYX(r1.theta, r2.theta)
@inline Base.:*(r1::RotYX, r2::RotX) = RotYX(r1.theta1, r1.theta2 + r2.theta)
@inline Base.:*(r1::RotY, r2::RotYX) = RotYX(r1.theta + r2.theta1, r2.theta2)

@inline inv(r::RotYX) = RotXY(-r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotYX}) = RotYX(0.0, 0.0)
@inline eye{T}(::Type{RotYX{T}}) = RotYX{T}(zero(T), zero(T))

"""
    immutable RotXZ{T} <: Rotation{T}
    RotXZ(theta_x, theta_z)

A 3×3 rotation matrix which represents a rotation by `theta_z` about the Z axis,
followed by a rotation by `theta_x` about the X axis.
"""
immutable RotXZ{T} <: Rotation{T}
    theta1::T
    theta2::T
end

@inline (::Type{R}){R<:RotXZ}(r::RotXZ) = R(r.theta1, r.theta2)
@inline convert{R<:RotXZ}(::Type{R}, r::RotXZ) = R(r.theta1, r.theta2)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 2 inputs (not 9).
@inline (::Type{RotXZ}){X,Z}(x::X, z::Z) = RotXZ{promote_type(X, Z)}(x, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{R}){R<:RotXZ}(t::NTuple{9}) = error("Cannot construct a two-axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotXZ{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotXZ{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    # transposed representation
    (cosθ₂,   cosθ₁*sinθ₂,  sinθ₁*sinθ₂,
     -sinθ₂,  cosθ₁*cosθ₂,  sinθ₁*cosθ₂,
     zero(T),       -sinθ₁,       cosθ₁)
end

@inline function Base.:*(r::RotXZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    return similar_type(v,T)(cosθ₂*v[1] + -sinθ₂*v[2],
                             cosθ₁*sinθ₂*v[1] + cosθ₁*cosθ₂*v[2] + -sinθ₁*v[3],
                             sinθ₁*sinθ₂*v[1] + sinθ₁*cosθ₂*v[2] + cosθ₁*v[3])
end

@inline Base.:*(r1::RotX, r2::RotZ) = RotXZ(r1.theta, r2.theta)
@inline Base.:*(r1::RotXZ, r2::RotZ) = RotXZ(r1.theta1, r1.theta2 + r2.theta)
@inline Base.:*(r1::RotX, r2::RotXZ) = RotXZ(r1.theta + r2.theta1, r2.theta2)

@inline inv(r::RotXZ) = RotZX(-r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotXZ}) = RotXZ(0.0, 0.0)
@inline eye{T}(::Type{RotXZ{T}}) = RotXZ{T}(zero(T), zero(T))

"""
    immutable RotZX{T} <: Rotation{T}
    RotZX(theta_z, theta_x)

A 3×3 rotation matrix which represents a rotation by `theta_x` about the X axis,
followed by a rotation by `theta_z` about the Z axis.
"""
immutable RotZX{T} <: Rotation{T}
    theta1::T
    theta2::T
end

@inline (::Type{R}){R<:RotZX}(r::RotZX) = R(r.theta1, r.theta2)
@inline convert{R<:RotZX}(::Type{R}, r::RotZX) = R(r.theta1, r.theta2)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 2 inputs (not 9).
@inline (::Type{RotZX}){Z,X}(z::Z, x::X) = RotZX{promote_type(Z, X)}(z, x)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{R}){R<:RotZX}(t::NTuple{9}) = error("Cannot construct a two-axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotZX{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotZX{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    # transposed representation
    ( cosθ₁,         sinθ₁,         zero(T),
     -sinθ₁*cosθ₂,   cosθ₁*cosθ₂,   sinθ₂,
      sinθ₁*sinθ₂,   cosθ₁*-sinθ₂,  cosθ₂)
end

@inline function Base.:*(r::RotZX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    return similar_type(v,T)(cosθ₁*v[1] + -sinθ₁*cosθ₂*v[2] + sinθ₁*sinθ₂*v[3],
                             sinθ₁*v[1] + cosθ₁*cosθ₂*v[2] + cosθ₁*-sinθ₂*v[3],
                             sinθ₂*v[2] + cosθ₂*v[3])
end

@inline Base.:*(r1::RotZ, r2::RotX) = RotZX(r1.theta, r2.theta)
@inline Base.:*(r1::RotZX, r2::RotX) = RotZX(r1.theta1, r1.theta2 + r2.theta)
@inline Base.:*(r1::RotZ, r2::RotZX) = RotZX(r1.theta + r2.theta1, r2.theta2)

@inline inv(r::RotZX) = RotXZ(-r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotZX}) = RotZX(0.0, 0.0)
@inline eye{T}(::Type{RotZX{T}}) = RotZX{T}(zero(T), zero(T))

"""
    immutable RotZY{T} <: Rotation{T}
    RotZY(theta_z, theta_y)

A 3×3 rotation matrix which represents a rotation by `theta_y` about the Y axis,
followed by a rotation by `theta_z` about the Z axis.
"""
immutable RotZY{T} <: Rotation{T}
    theta1::T
    theta2::T
end

@inline (::Type{R}){R<:RotZY}(r::RotZY) = R(r.theta1, r.theta2)
@inline convert{R<:RotZY}(::Type{R}, r::RotZY) = R(r.theta1, r.theta2)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 2 inputs (not 9).
@inline (::Type{RotZY}){Z,Y}(z::Z, y::Y) = RotZY{promote_type(Z, Y)}(z, y)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{R}){R<:RotZY}(t::NTuple{9}) = error("Cannot construct a two-axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotZY{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotZY{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    # transposed representation
    ( cosθ₁*cosθ₂,  sinθ₁*cosθ₂, -sinθ₂,
     -sinθ₁,        cosθ₁,        zero(T),
      cosθ₁*sinθ₂,  sinθ₁*sinθ₂,  cosθ₂)
end

@inline function Base.:*(r::RotZY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    return similar_type(v,T)(cosθ₁*cosθ₂*v[1] + -sinθ₁*v[2] + cosθ₁*sinθ₂*v[3],
                             sinθ₁*cosθ₂*v[1] + cosθ₁*v[2] + sinθ₁*sinθ₂*v[3],
                             -sinθ₂*v[1] + cosθ₂*v[3])
end

@inline Base.:*(r1::RotZ, r2::RotY) = RotZY(r1.theta, r2.theta)
@inline Base.:*(r1::RotZY, r2::RotY) = RotZY(r1.theta1, r1.theta2 + r2.theta)
@inline Base.:*(r1::RotZ, r2::RotZY) = RotZY(r1.theta + r2.theta1, r2.theta2)

@inline inv(r::RotZY) = RotYZ(-r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotZY}) = RotZY(0.0, 0.0)
@inline eye{T}(::Type{RotZY{T}}) = RotZY{T}(zero(T), zero(T))

"""
    immutable RotYZ{T} <: Rotation{T}
    RotYZ(theta_y, theta_z)

A 3×3 rotation matrix which represents a rotation by `theta_z` about the Z axis,
followed by a rotation by `theta_y` about the Y axis.
"""
immutable RotYZ{T} <: Rotation{T}
    theta1::T
    theta2::T
end

@inline (::Type{R}){R<:RotYZ}(r::RotYZ) = R(r.theta1, r.theta2)
@inline convert{R<:RotYZ}(::Type{R}, r::RotYZ) = R(r.theta1, r.theta2)

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 2 inputs (not 9).
@inline (::Type{RotYZ}){Y,Z}(y::Y, z::Z) = RotYZ{promote_type(Y, Z)}(y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{R}){R<:RotYZ}(t::NTuple{9}) = error("Cannot construct a two-axis rotation from a matrix")
@inline function Base.getindex{T}(r::RotYZ{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotYZ{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    # transposed representation
    (cosθ₁*cosθ₂,   sinθ₂,    -sinθ₁*cosθ₂,
     cosθ₁*-sinθ₂,  cosθ₂,     sinθ₁*sinθ₂,
     sinθ₁,         zero(T),      cosθ₁)
end

@inline function Base.:*(r::RotYZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)

    return similar_type(v,T)(cosθ₁*cosθ₂*v[1] + cosθ₁*-sinθ₂*v[2] + sinθ₁*v[3],
                             sinθ₂*v[1] + cosθ₂*v[2],
                             -sinθ₁*cosθ₂*v[1] + sinθ₁*sinθ₂*v[2] + cosθ₁*v[3])
end

@inline Base.:*(r1::RotY, r2::RotZ) = RotYZ(r1.theta, r2.theta)
@inline Base.:*(r1::RotYZ, r2::RotZ) = RotYZ(r1.theta1, r1.theta2 + r2.theta)
@inline Base.:*(r1::RotY, r2::RotYZ) = RotYZ(r1.theta + r2.theta1, r2.theta2)

# define null rotations for convenience
@inline eye(::Type{RotYZ}) = RotYZ(0.0, 0.0)
@inline eye{T}(::Type{RotYZ{T}}) = RotYZ{T}(zero(T), zero(T))

@inline inv(r::RotYZ) = RotZY(-r.theta2, -r.theta1)

function Base.rand{R <: Union{RotXY,RotYZ,RotZX, RotXZ, RotYX, RotZY}}(::Type{R})
    T = eltype(R)
    if T == Any
        T = Float64
    end

    # Not really sure what this distribution is, but it's also not clear what
    # it should be! rand(RotXY) *is* invariant to pre-rotations by a RotX and
    # post-rotations by a RotY...
    return R(2*pi*rand(T), 2*pi*rand(T))
end

################################################################################
################################################################################

##########################
# Proper Euler Rotations #
##########################

"""
    immutable RotXYX{T} <: Rotation{T}
    RotXYX(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" XYX Euler angle convention,
consisting of first a rotation about the X axis by `theta3`, followed by a
rotation about the Y axis by `theta2`, and finally a rotation about the X axis
by `theta1`.
"""
immutable RotXYX{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotXYX}){X,Y,Z}(x::X, y::Y, z::Z) = RotXYX{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotXYX}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[2, 1], (-R[3, 1] + eps(t[1])) - eps(t[1]))  # TODO: handle denormal numbers better, as atan2(0,0) != atan2(0,-0)
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2((R[1, 2] * R[1, 2] + R[1, 3] * R[1, 3])^(1/2), R[1, 1]),
        atan2(- R[2, 3]*ct1 - R[3, 3]*st1, R[2, 2]*ct1 + R[3, 2]*st1))
end
@inline function Base.getindex{T}(r::RotXYX{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotXYX{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₂,        sinθ₁*sinθ₂,                        cosθ₁*-sinθ₂,
     sinθ₂*sinθ₃,  cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃,   sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃,
     sinθ₂*cosθ₃,  cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃,  sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotXYX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        cosθ₂*v[1] + sinθ₂*sinθ₃*v[2] + sinθ₂*cosθ₃*v[3],
        sinθ₁*sinθ₂*v[1] + (cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃)*v[2] + (cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃)*v[3],
        cosθ₁*-sinθ₂*v[1] + (sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃)*v[2] + (sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃)*v[3])
end

@inline Base.:*(r1::RotX, r2::RotYX) = RotXYX(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotXY, r2::RotX) = RotXYX(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotXYX, r2::RotX) = RotXYX(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotX, r2::RotXYX) = RotXYX(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotXYX) = RotXYX(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotXYX}) = RotXYX(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotXYX{T}}) = RotXYX{T}(zero(T), zero(T), zero(T))

"""
    immutable RotXZX{T} <: Rotation{T}
    RotXZX(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" XZX Euler angle convention,
consisting of first a rotation about the X axis by `theta3`, followed by a
rotation about the Z axis by `theta2`, and finally a rotation about the X axis
by `theta1`.
"""
immutable RotXZX{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotXZX}){X,Y,Z}(x::X, y::Y, z::Z) = RotXZX{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotXZX}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[3, 1], R[2, 1])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2((R[1, 2] * R[1, 2] + R[1, 3] * R[1, 3])^(1/2), R[1, 1]),
        atan2(R[3, 2]*ct1 - R[2, 2]*st1, R[3, 3]*ct1 - R[2, 3]*st1))
end
@inline function Base.getindex{T}(r::RotXZX{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotXZX{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₂,         cosθ₁*sinθ₂,                        sinθ₁*sinθ₂,
     -sinθ₂*cosθ₃,  cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃,   sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃,
     sinθ₂*sinθ₃,   cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃,  sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotXZX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        cosθ₂*v[1] + -sinθ₂*cosθ₃*v[2] + sinθ₂*sinθ₃*v[3],
        cosθ₁*sinθ₂*v[1] + (cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃)*v[2] + (cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃)*v[3],
        sinθ₁*sinθ₂*v[1] + (sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[2] + (sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

@inline Base.:*(r1::RotX, r2::RotZX) = RotXZX(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotXZ, r2::RotX) = RotXZX(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotXZX, r2::RotX) = RotXZX(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotX, r2::RotXZX) = RotXZX(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotXZX) = RotXZX(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotXZX}) = RotXZX(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotXZX{T}}) = RotXZX{T}(zero(T), zero(T), zero(T))

"""
    immutable RotYXY{T} <: Rotation{T}
    RotYXY(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" YXY Euler angle convention,
consisting of first a rotation about the Y axis by `theta3`, followed by a
rotation about the X axis by `theta2`, and finally a rotation about the Y axis
by `theta1`.
"""
immutable RotYXY{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotYXY}){X,Y,Z}(x::X, y::Y, z::Z) = RotYXY{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotYXY}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[1, 2], R[3, 2])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2((R[2, 1] * R[2, 1] + R[2, 3] * R[2, 3])^(1/2), R[2, 2]),
        atan2(R[1, 3]*ct1 - R[3, 3]*st1, R[1, 1]*ct1 - R[3, 1]*st1))
end
@inline function Base.getindex{T}(r::RotYXY{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotYXY{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₃ + sinθ₁*cosθ₂*-sinθ₃,  sinθ₂*sinθ₃,   -sinθ₁*cosθ₃ + cosθ₁*cosθ₂*-sinθ₃,
     sinθ₁*sinθ₂,                       cosθ₂,         cosθ₁*sinθ₂,
     cosθ₁*sinθ₃ + sinθ₁*cosθ₂*cosθ₃,   -sinθ₂*cosθ₃,  -sinθ₁*sinθ₃ + cosθ₁*cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotYXY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        (cosθ₁*cosθ₃ + sinθ₁*cosθ₂*-sinθ₃)*v[1] + sinθ₁*sinθ₂*v[2] + (cosθ₁*sinθ₃ + sinθ₁*cosθ₂*cosθ₃)*v[3],
        sinθ₂*sinθ₃*v[1] + cosθ₂*v[2] + -sinθ₂*cosθ₃*v[3],
        (-sinθ₁*cosθ₃ + cosθ₁*cosθ₂*-sinθ₃)*v[1] + cosθ₁*sinθ₂*v[2] + (-sinθ₁*sinθ₃ + cosθ₁*cosθ₂*cosθ₃)*v[3])
end

@inline Base.:*(r1::RotY, r2::RotXY) = RotYXY(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotYX, r2::RotY) = RotYXY(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotYXY, r2::RotY) = RotYXY(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotY, r2::RotYXY) = RotYXY(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotYXY) = RotYXY(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotYXY}) = RotYXY(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotYXY{T}}) = RotYXY{T}(zero(T), zero(T), zero(T))

"""
    immutable RotYZY{T} <: Rotation{T}
    RotYZY(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" YXY Euler angle convention,
consisting of first a rotation about the Y axis by `theta3`, followed by a
rotation about the Z axis by `theta2`, and finally a rotation about the Y axis
by `theta1`.
"""
immutable RotYZY{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotYZY}){X,Y,Z}(x::X, y::Y, z::Z) = RotYZY{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotYZY}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[3, 2], -R[1, 2])  # TODO: handle denormal numbers better, as atan2(0,0) != atan2(0,-0)
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2((R[2, 1] * R[2, 1] + R[2, 3] * R[2, 3])^(1/2), R[2, 2]),
        atan2(- R[3, 1]*ct1 - R[1, 1]*st1, R[3, 3]*ct1 + R[1, 3]*st1))
end
@inline function Base.getindex{T}(r::RotYZY{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotYZY{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₂*cosθ₃ + sinθ₁*-sinθ₃,  sinθ₂*cosθ₃,  -sinθ₁*cosθ₂*cosθ₃ + cosθ₁*-sinθ₃,
     cosθ₁*-sinθ₂,                      cosθ₂,        sinθ₁*sinθ₂,
     cosθ₁*cosθ₂*sinθ₃ + sinθ₁*cosθ₃,   sinθ₂*sinθ₃,  -sinθ₁*cosθ₂*sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotYZY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        (cosθ₁*cosθ₂*cosθ₃ + sinθ₁*-sinθ₃)*v[1] + cosθ₁*-sinθ₂*v[2] + (cosθ₁*cosθ₂*sinθ₃ + sinθ₁*cosθ₃)*v[3],
        sinθ₂*cosθ₃*v[1] + cosθ₂*v[2] + sinθ₂*sinθ₃*v[3],
        (-sinθ₁*cosθ₂*cosθ₃ + cosθ₁*-sinθ₃)*v[1] + sinθ₁*sinθ₂*v[2] + (-sinθ₁*cosθ₂*sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

@inline Base.:*(r1::RotY, r2::RotZY) = RotYZY(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotYZ, r2::RotY) = RotYZY(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotYZY, r2::RotY) = RotYZY(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotY, r2::RotYZY) = RotYZY(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotYZY) = RotYZY(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotYZY}) = RotYZY(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotYZY{T}}) = RotYZY{T}(zero(T), zero(T), zero(T))

"""
    immutable RotZXZ{T} <: Rotation{T}
    RotZXZ(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" ZXZ Euler angle convention,
consisting of first a rotation about the Z axis by `theta3`, followed by a
rotation about the X axis by `theta2`, and finally a rotation about the Z axis
by `theta1`.
"""
immutable RotZXZ{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotZXZ}){X,Y,Z}(x::X, y::Y, z::Z) = RotZXZ{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotZXZ}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[1, 3], (-R[2, 3] + eps()) - eps())  # TODO: handle denormal numbers better, as atan2(0,0) != atan2(0,-0)
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2((R[3, 1] * R[3, 1] + R[3, 2] * R[3, 2])^(1/2), R[3, 3]),
        atan2(- R[1, 2]*ct1 - R[2, 2]*st1, R[1, 1]*ct1 + R[2, 1]*st1))
end

@inline function Base.getindex{T}(r::RotZXZ{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotZXZ{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃,   sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃,   sinθ₂*sinθ₃,
     cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃,  sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃,  sinθ₂*cosθ₃,
     sinθ₁*sinθ₂,                        cosθ₁*-sinθ₂,                      cosθ₂)
end

@inline function Base.:*(r::RotZXZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
         (cosθ₁*cosθ₃ + -sinθ₁*cosθ₂*sinθ₃)*v[1] + (cosθ₁*-sinθ₃ + -sinθ₁*cosθ₂*cosθ₃)*v[2] + sinθ₁*sinθ₂*v[3],
         (sinθ₁*cosθ₃ + cosθ₁*cosθ₂*sinθ₃)*v[1] + (sinθ₁*-sinθ₃ + cosθ₁*cosθ₂*cosθ₃)*v[2] + cosθ₁*-sinθ₂*v[3],
         sinθ₂*sinθ₃*v[1] + sinθ₂*cosθ₃*v[2] + cosθ₂*v[3])
end

@inline Base.:*(r1::RotZ, r2::RotXZ) = RotZXZ(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotZX, r2::RotZ) = RotZXZ(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotZXZ, r2::RotZ) = RotZXZ(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotZ, r2::RotZXZ) = RotZXZ(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotZXZ) = RotZXZ(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotZXZ}) = RotZXZ(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotZXZ{T}}) = RotZXZ{T}(zero(T), zero(T), zero(T))

"""
    immutable RotZYZ{T} <: Rotation{T}
    RotZYZ(theta1, theta2, theta3)

A 3×3 rotation matrix parameterized by the "proper" ZXZ Euler angle convention,
consisting of first a rotation about the Z axis by `theta3`, followed by a
rotation about the Y axis by `theta2`, and finally a rotation about the Z axis
by `theta1`.
"""
immutable RotZYZ{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotZYZ}){X,Y,Z}(x::X, y::Y, z::Z) = RotZYZ{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotZYZ}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[2, 3], R[1, 3])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2((R[3, 1] * R[3, 1] + R[3, 2] * R[3, 2])^(1/2), R[3, 3]),
        atan2(R[2, 1]*ct1 - R[1, 1]*st1, R[2, 2]*ct1 - R[1, 2]*st1))
end

@inline function Base.getindex{T}(r::RotZYZ{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotZYZ{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃,   sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃,   -sinθ₂*cosθ₃,
     cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃,  sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃,  sinθ₂*sinθ₃,
     cosθ₁*sinθ₂,                        sinθ₁*sinθ₂,                       cosθ₂)
end

@inline function Base.:*(r::RotZYZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        (cosθ₁*cosθ₂*cosθ₃ + -sinθ₁*sinθ₃)*v[1] + (cosθ₁*cosθ₂*-sinθ₃ + -sinθ₁*cosθ₃)*v[2] + cosθ₁*sinθ₂*v[3],
        (sinθ₁*cosθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[1] + (sinθ₁*cosθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[2] + sinθ₁*sinθ₂*v[3],
        -sinθ₂*cosθ₃*v[1] + sinθ₂*sinθ₃*v[2] + cosθ₂*v[3])
end

@inline Base.:*(r1::RotZ, r2::RotYZ) = RotZYZ(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotZY, r2::RotZ) = RotZYZ(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotZYZ, r2::RotZ) = RotZYZ(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotZ, r2::RotZYZ) = RotZYZ(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotZYZ) = RotZYZ(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotZYZ}) = RotZYZ(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotZYZ{T}}) = RotZYZ{T}(zero(T), zero(T), zero(T))

###############################
# Tait-Bryant Euler Rotations #
###############################

"""
    immutable RotXYZ{T} <: Rotation{T}
    RotXYZ(theta1, theta2, theta3)
    RotXYZ(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" XYZ Euler angle
convention, consisting of first a rotation about the Z axis by `theta3`,
followed by a rotation about the Y axis by `theta2`, and finally a rotation
about the X axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in XYZ order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
immutable RotXYZ{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotXYZ}){X,Y,Z}(x::X, y::Y, z::Z) = RotXYZ{promote_type(promote_type(X, Y), Z)}(x, y, z)
@inline (::Type{Rot}){Rot<:RotXYZ}(; roll=0, pitch=0, yaw=0) = Rot(roll, pitch, yaw)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotXYZ}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(-R[2, 3], R[3, 3])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2(R[1, 3], (R[1, 1] * R[1, 1] + R[1, 2] * R[1, 2])^(1/2)),
        atan2(R[2, 1]*ct1 + R[3, 1]*st1, R[2, 2]*ct1 + R[3, 2]*st1))
end

@inline function Base.getindex{T}(r::RotXYZ{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotXYZ{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₂*cosθ₃,   sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃,   cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃,
     cosθ₂*-sinθ₃,  sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃,  cosθ₁*-sinθ₂*-sinθ₃ + sinθ₁*cosθ₃,
     sinθ₂,         -sinθ₁*cosθ₂,                      cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotXYZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        cosθ₂*cosθ₃*v[1] + cosθ₂*-sinθ₃*v[2] + sinθ₂*v[3],
        (sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[1] + (sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[2] + -sinθ₁*cosθ₂*v[3],
        (cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃)*v[1] + (cosθ₁*sinθ₂*sinθ₃ + sinθ₁*cosθ₃)*v[2] + cosθ₁*cosθ₂*v[3])
end

@inline Base.:*(r1::RotX, r2::RotYZ) = RotXYZ(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotXY, r2::RotZ) = RotXYZ(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotXYZ, r2::RotZ) = RotXYZ(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotX, r2::RotXYZ) = RotXYZ(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotXYZ) = RotZYX(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotXYZ}) = RotXYZ(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotXYZ{T}}) = RotXYZ{T}(zero(T), zero(T), zero(T))

"""
    immutable RotZYX{T} <: Rotation{T}
    RotZYX(theta1, theta2, theta3)
    RotZYX(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" ZYX Euler angle
convention, consisting of first a rotation about the X axis by `theta3`,
followed by a rotation about the Y axis by `theta2`, and finally a rotation
about the Z axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in ZYX order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
immutable RotZYX{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotZYX}){X,Y,Z}(x::X, y::Y, z::Z) = RotZYX{promote_type(promote_type(X, Y), Z)}(x, y, z)
@inline (::Type{Rot}){Rot<:RotZYX}(; roll=0, pitch=0, yaw=0) = Rot(yaw, pitch, roll)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotZYX}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[2, 1], R[1, 1])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2(-R[3, 1], (R[3, 2] * R[3, 2] + R[3, 3] * R[3, 3])^(1/2)),
        atan2(R[1, 3]*st1 - R[2, 3]*ct1, R[2, 2]*ct1 - R[1, 2]*st1))
end

@inline function Base.getindex{T}(r::RotZYX{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotZYX{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    ( cosθ₁*cosθ₂,                       sinθ₁*cosθ₂,                      -sinθ₂,
     -sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃,   cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃,   cosθ₂*sinθ₃,
      sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃,   cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃,  cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotZYX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        cosθ₁*cosθ₂*v[1] + (-sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃)*v[2] + (sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃)*v[3],
        sinθ₁*cosθ₂*v[1] + (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃)*v[2] + (cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃)*v[3],
        -sinθ₂*v[1] + cosθ₂*sinθ₃*v[2] + cosθ₂*cosθ₃*v[3])
end

@inline Base.:*(r1::RotZ, r2::RotYX) = RotZYX(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotZY, r2::RotX) = RotZYX(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotZYX, r2::RotX) = RotZYX(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotZ, r2::RotZYX) = RotZYX(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotZYX) = RotXYZ(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotZYX}) = RotZYX(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotZYX{T}}) = RotZYX{T}(zero(T), zero(T), zero(T))

"""
    immutable RotXZY{T} <: Rotation{T}
    RotXZY(theta1, theta2, theta3)
    RotXZY(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" XZY Euler angle
convention, consisting of first a rotation about the Y axis by `theta3`,
followed by a rotation about the Z axis by `theta2`, and finally a rotation
about the X axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in XZY order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
immutable RotXZY{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotXZY}){X,Y,Z}(x::X, y::Y, z::Z) = RotXZY{promote_type(promote_type(X, Y), Z)}(x, y, z)
@inline (::Type{Rot}){Rot<:RotXZY}(; roll=0, pitch=0, yaw=0) = Rot(roll, yaw, pitch)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotXZY}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[3, 2], R[2, 2])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2(-R[1, 2], (R[1, 1] * R[1, 1] + R[1, 3] * R[1, 3])^(1/2)),
        atan2(R[2, 1]*st1 - R[3, 1]*ct1, R[3, 3]*ct1 - R[2, 3]*st1))
end

@inline function Base.getindex{T}(r::RotXZY{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotXZY{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    ( cosθ₂*cosθ₃,  cosθ₁*sinθ₂*cosθ₃ + sinθ₁*sinθ₃,   sinθ₁*sinθ₂*cosθ₃ + cosθ₁*-sinθ₃,
     -sinθ₂,        cosθ₁*cosθ₂,                       sinθ₁*cosθ₂,
      cosθ₂*sinθ₃,  cosθ₁*sinθ₂*sinθ₃ + -sinθ₁*cosθ₃,  sinθ₁*sinθ₂*sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotXZY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        cosθ₂*cosθ₃*v[1] + -sinθ₂*v[2] + cosθ₂*sinθ₃*v[3],
        (cosθ₁*sinθ₂*cosθ₃ + sinθ₁*sinθ₃)*v[1] + cosθ₁*cosθ₂*v[2] + (cosθ₁*sinθ₂*sinθ₃ + -sinθ₁*cosθ₃)*v[3],
        (sinθ₁*sinθ₂*cosθ₃ + cosθ₁*-sinθ₃)*v[1] + sinθ₁*cosθ₂*v[2] + (sinθ₁*sinθ₂*sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

@inline Base.:*(r1::RotX, r2::RotZY) = RotXZY(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotXZ, r2::RotY) = RotXZY(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotXZY, r2::RotY) = RotXZY(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotX, r2::RotXZY) = RotXZY(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotXZY) = RotYZX(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotXZY}) = RotXZY(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotXZY{T}}) = RotXZY{T}(zero(T), zero(T), zero(T))

"""
    immutable RotYZX{T} <: Rotation{T}
    RotYZX(theta1, theta2, theta3)
    RotYZX(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" YZX Euler angle
convention, consisting of first a rotation about the X axis by `theta3`,
followed by a rotation about the Z axis by `theta2`, and finally a rotation
about the Y axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in YZX order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
immutable RotYZX{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotYZX}){X,Y,Z}(x::X, y::Y, z::Z) = RotYZX{promote_type(promote_type(X, Y), Z)}(x, y, z)
@inline (::Type{Rot}){Rot<:RotYZX}(; roll=0, pitch=0, yaw=0) = Rot(pitch, yaw, roll)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotYZX}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(-R[3, 1], R[1, 1])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2(R[2, 1], (R[2, 2] * R[2, 2] + R[2, 3] * R[2, 3])^(1/2)),
        atan2(R[3, 2]*ct1 + R[1, 2]*st1, R[3, 3]*ct1 + R[1, 3]*st1))
end

@inline function Base.getindex{T}(r::RotYZX{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotYZX{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₂,                        sinθ₂,         -sinθ₁*cosθ₂,
     cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃,   cosθ₂*cosθ₃,   sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃,
     cosθ₁*sinθ₂*sinθ₃ + sinθ₁*cosθ₃,    cosθ₂*-sinθ₃,  sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃)
end

@inline function Base.:*(r::RotYZX, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        cosθ₁*cosθ₂*v[1] + (cosθ₁*-sinθ₂*cosθ₃ + sinθ₁*sinθ₃)*v[2] + (cosθ₁*-sinθ₂*-sinθ₃ + sinθ₁*cosθ₃)*v[3],
        sinθ₂*v[1] + cosθ₂*cosθ₃*v[2] + cosθ₂*-sinθ₃*v[3],
        -sinθ₁*cosθ₂*v[1] + (sinθ₁*sinθ₂*cosθ₃ + cosθ₁*sinθ₃)*v[2] + (sinθ₁*sinθ₂*-sinθ₃ + cosθ₁*cosθ₃)*v[3])
end

@inline Base.:*(r1::RotY, r2::RotZX) = RotYZX(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotYZ, r2::RotX) = RotYZX(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotYZX, r2::RotX) = RotYZX(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotY, r2::RotYZX) = RotYZX(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotYZX) = RotXZY(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotYZX}) = RotYZX(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotYZX{T}}) = RotYZX{T}(zero(T), zero(T), zero(T))

"""
    immutable RotYXZ{T} <: Rotation{T}
    RotYXZ(theta1, theta2, theta3)
    RotYXZ(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" YXZ Euler angle
convention, consisting of first a rotation about the Z axis by `theta3`,
followed by a rotation about the X axis by `theta2`, and finally a rotation
about the Y axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in YXZ order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
immutable RotYXZ{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotYXZ}){X,Y,Z}(x::X, y::Y, z::Z) = RotYXZ{promote_type(promote_type(X, Y), Z)}(x, y, z)
@inline (::Type{Rot}){Rot<:RotYXZ}(; roll=0, pitch=0, yaw=0) = Rot(pitch, roll, yaw)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotYXZ}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(R[1, 3], R[3, 3])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2(-R[2, 3], (R[2, 1] * R[2, 1] + R[2, 2] * R[2, 2])^(1/2)),
        atan2(R[3, 2]*st1 - R[1, 2]*ct1, R[1, 1]*ct1 - R[3, 1]*st1))
end

@inline function Base.getindex{T}(r::RotYXZ{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotYXZ{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃,   cosθ₂*sinθ₃,  -sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃,
     cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃,  cosθ₂*cosθ₃,  sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃,
     sinθ₁*cosθ₂,                       -sinθ₂,       cosθ₁*cosθ₂)
end

@inline function Base.:*(r::RotYXZ, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*sinθ₃)*v[1] + (cosθ₁*-sinθ₃ + sinθ₁*sinθ₂*cosθ₃)*v[2] + sinθ₁*cosθ₂*v[3],
        cosθ₂*sinθ₃*v[1] + cosθ₂*cosθ₃*v[2] + -sinθ₂*v[3],
        (-sinθ₁*cosθ₃ + cosθ₁*sinθ₂*sinθ₃)*v[1] + (sinθ₁*sinθ₃ + cosθ₁*sinθ₂*cosθ₃)*v[2] + cosθ₁*cosθ₂*v[3])
end

@inline Base.:*(r1::RotY, r2::RotXZ) = RotYXZ(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotYX, r2::RotZ) = RotYXZ(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotYXZ, r2::RotZ) = RotYXZ(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotY, r2::RotYXZ) = RotYXZ(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotYXZ) = RotZXY(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotYXZ}) = RotYXZ(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotYXZ{T}}) = RotYXZ{T}(zero(T), zero(T), zero(T))

"""
    immutable RotZXY{T} <: Rotation{T}
    RotZXY(theta1, theta2, theta3)
    RotZXY(roll=r, pitch=p, yaw=y)

A 3×3 rotation matrix parameterized by the "Tait-Bryant" ZXY Euler angle
convention, consisting of first a rotation about the Y axis by `theta3`,
followed by a rotation about the X axis by `theta2`, and finally a rotation
about the Z axis by `theta1`.

The keyword argument form applies roll, pitch and yaw to the X, Y and Z axes
respectively, in ZXY order. (Because it is a right-handed coordinate system,
note that positive pitch is heading in the negative Z axis).
"""
immutable RotZXY{T} <: Rotation{T}
    theta1::T
    theta2::T
    theta3::T
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{RotZXY}){X,Y,Z}(x::X, y::Y, z::Z) = RotZXY{promote_type(promote_type(X, Y), Z)}(x, y, z)
@inline (::Type{Rot}){Rot<:RotZXY}(; roll=0, pitch=0, yaw=0) = Rot(yaw, roll, pitch)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline function (::Type{Rot}){Rot <: RotZXY}(t::NTuple{9})
    R = SMatrix{3,3}(t)

    t1 = atan2(-R[1, 2], R[2, 2])
    ct1, st1 = cos(t1), sin(t1)

    Rot(t1,
        atan2(R[3, 2], (R[3, 1] * R[3, 1] + R[3, 3] * R[3, 3])^(1/2)),
        atan2(R[1, 3]*ct1 + R[2, 3]*st1, R[1, 1]*ct1 + R[2, 1]*st1))
end

@inline function Base.getindex{T}(r::RotZXY{T}, i::Integer)
    Tuple(r)[i] # Slow...
end

@inline function Base.convert{T}(::Type{Tuple}, r::RotZXY{T})
    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    # transposed representation
    ( cosθ₁*cosθ₃ + sinθ₁*sinθ₂*-sinθ₃,  sinθ₁*cosθ₃ + cosθ₁*-sinθ₂*-sinθ₃,  cosθ₂*-sinθ₃,
     -sinθ₁*cosθ₂,                       cosθ₁*cosθ₂,                        sinθ₂,
      cosθ₁*sinθ₃ + sinθ₁*sinθ₂*cosθ₃,   sinθ₁*sinθ₃ + cosθ₁*-sinθ₂*cosθ₃,   cosθ₂*cosθ₃)
end

@inline function Base.:*(r::RotZXY, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    T = promote_type(eltype(r), eltype(v))

    sinθ₁ = sin(r.theta1)
    cosθ₁ = cos(r.theta1)
    sinθ₂ = sin(r.theta2)
    cosθ₂ = cos(r.theta2)
    sinθ₃ = sin(r.theta3)
    cosθ₃ = cos(r.theta3)

    return similar_type(v,T)(
        (cosθ₁*cosθ₃ + sinθ₁*sinθ₂*-sinθ₃)*v[1] + -sinθ₁*cosθ₂*v[2] + (cosθ₁*sinθ₃ + sinθ₁*sinθ₂*cosθ₃)*v[3],
        (sinθ₁*cosθ₃ + cosθ₁*-sinθ₂*-sinθ₃)*v[1] + cosθ₁*cosθ₂*v[2] + (sinθ₁*sinθ₃ + cosθ₁*-sinθ₂*cosθ₃)*v[3],
         cosθ₂*-sinθ₃*v[1] + sinθ₂*v[2] + cosθ₂*cosθ₃*v[3])
end

@inline Base.:*(r1::RotZ, r2::RotXY) = RotZXY(r1.theta, r2.theta1, r2.theta2)
@inline Base.:*(r1::RotZX, r2::RotY) = RotZXY(r1.theta1, r1.theta2, r2.theta)
@inline Base.:*(r1::RotZXY, r2::RotY) = RotZXY(r1.theta1, r1.theta2, r1.theta3 + r2.theta)
@inline Base.:*(r1::RotZ, r2::RotZXY) = RotZXY(r1.theta + r2.theta1, r2.theta2, r2.theta3)

@inline inv(r::RotZXY) = RotYXZ(-r.theta3, -r.theta2, -r.theta1)

# define null rotations for convenience
@inline eye(::Type{RotZXY}) = RotZXY(0.0, 0.0, 0.0)
@inline eye{T}(::Type{RotZXY{T}}) = RotZXY{T}(zero(T), zero(T), zero(T))

ext: jl


text:
__precompile__()

## TODO:
## * tidy up code

module SymPy

using Compat

"""

SymPy package to interface with Python's [SymPy package](http://www.sympy.org) through `PyCall`.

The basic idea is that a new type -- `Sym` -- is made to hold symbolic
objects.  For this type, the basic operators and appropriate functions
of `Julia` are overloaded for `Sym` objects so that the expressions
are treated symbolically and not evaluated immediately. Instances of
this type are created by the constructor `Sym` or `symbols` or the macro
`@vars`.

As well, many -- but not all -- of the SymPy functions are ported to
allow them to be called on `Sym` objects. Mostly these are implemented
through metaprogramming, so adding missing functions is not hard. They
are not generated automatically though, rather added by hand.

To find documentation on SymPy functions, one should refer to
SymPy's [website](http://docs.sympy.org/latest/index.html).

Plotting is provided through the `Plots` interface. For details, see the help page for `sympy_plotting`.

The package tutorial provides many examples. This can be read on
[GitHub](https://github.com/jverzani/SymPy.jl/blob/master/examples/tutorial.ipynb).

"""
SymPy

using Compat

using PyCall

import Base: show
import Base: convert, promote_rule
import Base: getindex
import Base: start, next, done
import Base: complex
import Base: sin, cos, tan, sinh, cosh, tanh, asin, acos,
       atan, asinh, acosh, atanh, sec, csc, cot, asec,
       acsc, acot, sech, csch, coth, asech, acsch, acoth,
       sinc, cosc, cosd, cotd, cscd, secd, sind, tand,
       acosd, acotd, acscd, asecd, asind, atand, atan2,
       sinpi, cospi,
       log, log2,
       log10, log1p, exponent, exp, exp2, expm1, cbrt, sqrt,
       erf, erfc, erfcx, erfi, erfinv, erfcinv, dawson, ceil, floor,
       trunc, round, significand,
       abs, max, min, maximum, minimum,
       sign, dot,
       besseli, besselj, besselk, bessely,
       airyai, airybi,
       zero, one
import Base: transpose
import Base: diff
import Base: factorial, gcd, lcm, isqrt
import Base: gamma, beta
import Base: length,  size
import Base: factor, expand, collect
import Base: !=, ==
import Base:  inv, conj,
              cross, eigvals, eigvecs, trace, norm
import Base: promote_rule
import Base: match, replace, round
import Base: +, -, *, /, //, \
import Base: ^, .^
import Base: &, |, !, >, >=, ==, <=, <
## poly.jl
import Base: div
import Base: trunc
import Base: isinf, isnan
import Base: real, imag
import Base: expm
import Base: nullspace

export sympy, sympy_meth, object_meth, call_matrix_meth
export Sym, @syms, @vars, @osyms, symbols
export pprint,  jprint
export SymFunction, @symfuns,
       SymMatrix,
       evalf, N,  subs,
       simplify, nsimplify,
       expand, factor, trunc,
       collect, separate,
       fraction,
       primitive, sqf, resultant, cancel,
       together, square,
       solve,
       limit,
       series, integrate,
       summation,
       dsolve,
       poly,  nroots, real_roots, polyroots,
       ∨, ∧, ¬,
       rhs, lhs, args,
       jacobian, hessian,
       Max, Min,
       rref
export PI, E, IM, oo
export relation, piecewise, Piecewise, piecewise_fold
export members, doc, _str

## Following PyPlot, we initialize our variables outside _init_
const sympy  = PyCall.PyNULL()
const mpmath = PyCall.PyNULL()

include("types.jl")
include("utils.jl")
include("mathops.jl")
include("core.jl")
include("logical.jl")
include("math.jl")
include("mpmath.jl")
include("specialfuns.jl")
include("solve.jl")
include("dsolve.jl")
include("subs.jl")
include("patternmatch.jl")
include("simplify.jl")
include("series.jl")
include("integrate.jl")
include("assumptions.jl")
include("poly.jl")
include("matrix.jl")
include("ntheory.jl")
include("sets.jl")
include("display.jl")
include("lambdify.jl")

## add call interface depends on version
VERSION >= v"0.5.0-" && include("call.jl")
v"0.4.0" <= VERSION < v"0.5.0-" && include("call-0.4.jl")

include("plot_recipes.jl") # hook into Plots

## create some methods

for meth in union(core_sympy_methods,
                  simplify_sympy_meths,
                  expand_sympy_meths,
                  functions_sympy_methods,
                  series_sympy_meths,
                  integrals_sympy_methods,
                  summations_sympy_methods,
                  logic_sympy_methods,
                  polynomial_sympy_methods,
                  ntheory_sympy_methods
                  )

    meth_name = string(meth)
    @eval begin
        @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth){T<:SymbolicObject}(ex::T, args...; kwargs...) = sympy_meth($meth_name, ex, args...; kwargs...)

    end
    eval(Expr(:export, meth))
end

for meth in union(core_object_methods,
                  integrals_instance_methods,
                  summations_instance_methods,
                  polynomial_instance_methods
                  )

    meth_name = string(meth)
    @eval begin
        @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(ex::SymbolicObject, args...; kwargs...) = object_meth(ex, $meth_name, args...; kwargs...)
    end
    eval(Expr(:export, meth))
end

for prop in union(core_object_properties,
                  summations_object_properties,
                  polynomial_predicates)

    prop_name = string(prop)
    @eval ($prop)(ex::Sym) = ex[@compat(Symbol($prop_name))]
    eval(Expr(:export, prop))
end

## For precompilation we must put PyCall instances in __init__:
function __init__()

    ## Define sympy, mpmath, ...
    copy!(sympy, PyCall.pyimport_conda("sympy", "sympy"))

    ## mappings from PyObjects to types.
    basictype = sympy[:basic]["Basic"]
    pytype_mapping(basictype, Sym)

    polytype = sympy[:polys]["polytools"]["Poly"]
    pytype_mapping(polytype, Sym)

    try
        matrixtype = sympy[:matrices]["MatrixBase"]
        pytype_mapping(matrixtype, SymMatrix)
        pytype_mapping(sympy[:Matrix], SymMatrix)
    catch e
    end

    ## Makes it possible to call in a sympy method, witout worrying about Sym objects

    global call_sympy_fun(fn::Function, args...; kwargs...) = fn(args...; kwargs...)
    global call_sympy_fun(fn::PyCall.PyObject, args...; kwargs...) = call_sympy_fun(convert(Function, fn), args...; kwargs...)

    ## Main interface to methods in sympy
    ## sympy_meth(:name, ars, kwars...)
    global sympy_meth(meth, args...; kwargs...) = begin
        ans = call_sympy_fun(convert(Function, sympy[@compat(Symbol(meth))]), args...; kwargs...)
        ## make nicer...
        try
            if isa(ans, Vector)
                ans = Sym[i for i in ans]
            end
        catch err
        end
        ans
    end
    global object_meth(object::SymbolicObject, meth, args...; kwargs...)  =  begin
        call_sympy_fun(project(object)[@compat(Symbol(meth))],  args...; kwargs...)
    end
    global call_matrix_meth(object::SymbolicObject, meth, args...; kwargs...) = begin
        out = object_meth(object, meth, args...; kwargs...)
        if isa(out, SymMatrix)
            convert(Array{Sym}, out)
        elseif  length(out) == 1
            out
        else
            map(u -> isa(u, SymMatrix) ? convert(Array{Sym}, u) : u, out)
        end
    end

    ##
    init_logical()
    init_math()
    init_mpmath()
    init_sets()
    init_lambdify()
end

end


text:
const LEFT_PARTICLES_COUNT = 320
const RIGHT_PARTICLES_COUNT = 80
const PARTICLES_COUNT = LEFT_PARTICLES_COUNT + RIGHT_PARTICLES_COUNT

const RIGHT_MIN = 0
const RIGHT_MAX = 0.6
const LEFT_MIN = -0.6
const LEFT_MAX = RIGHT_MIN - (RIGHT_MIN - LEFT_MIN) / (LEFT_PARTICLES_COUNT + 2)

const LEFT_DENSITY = 1.0
const LEFT_VELOCITY= 0.0
const LEFT_ENERGY = 2.5
const LEFT_PRESSURE = 1.0

const RIGHT_DENSITY = 0.25
const RIGHT_VELOCITY= 0.0
const RIGHT_ENERGY = 1.795
const RIGHT_PRESSURE = 0.1795

const PARTICLE_MASS = 0.001875
const GAMMA = 1.4

const DX = 0.6/80
const HSML = 2.0 * DX

const TIME_STEPS = 40
const TIME_DELTA = 0.005

const FIXED_SIDES_PARTICLES_COUNT = 20

function generateQuantities()
    # Position
    x = append!(
        [(i)*DX/4 + LEFT_MIN for i in 1:LEFT_PARTICLES_COUNT],
        [(i-LEFT_PARTICLES_COUNT-1)*DX + DX/2 for i in LEFT_PARTICLES_COUNT+1:PARTICLES_COUNT]
    )
    #[(i-LEFT_PARTICLES_COUNT)*DX for i in LEFT_PARTICLES_COUNT+1:PARTICLES_COUNT]
    # Density
    ρ = append!(LEFT_DENSITY .* ones(LEFT_PARTICLES_COUNT), RIGHT_DENSITY .* ones(RIGHT_PARTICLES_COUNT))
    # Energy
    u = append!(LEFT_ENERGY .* ones(LEFT_PARTICLES_COUNT), RIGHT_PRESSURE .* ones(RIGHT_PARTICLES_COUNT))
    # Velocity
    v = zeros(PARTICLES_COUNT)
    # Pressure
    p = append!(LEFT_PRESSURE .* ones(LEFT_PARTICLES_COUNT), RIGHT_PRESSURE .* ones(RIGHT_PARTICLES_COUNT))

    (x, ρ, u, v, p)
end


text:
module HomotopyContinuation

    import DynamicPolynomials
    import ElasticArrays
    import FixedPolynomials
    import LinearAlgebra
    import MixedSubdivisions
    import MultivariatePolynomials
    import PrettyTables
    import Printf
    import ProjectiveVectors
    import Random
    import StaticArrays
    import StaticPolynomials
    import TreeViews

    import LinearAlgebra: cond
    import Parameters: @pack!, @unpack
    import DynamicPolynomials: @polyvar, subs, differentiate
    import ProjectiveVectors: PVector
    import StaticArrays: SVector, @SVector
    import Test: @test
    import MixedSubdivisions: mixed_volume

    const FP = FixedPolynomials
    const MP = MultivariatePolynomials
    const SP = StaticPolynomials

    export @polyvar, subs, differentiate
    export mixed_volume
    export cond

    include("progress_meter.jl")
    import .ProgressMeter

    include("utilities.jl")
    include("affine_patches.jl")

    include("systems_and_homotopies.jl")
    include("input.jl")
    include("problems.jl")
    include("totaldegree.jl")

    include("newton.jl")
    include("predictors.jl")
    include("correctors.jl")

    include("core_tracker.jl")
    include("path_tracker.jl")
    include("polyhedral.jl")
    include("solve.jl")
    include("monodromy.jl")
end


text:
################################################################
# Code part of FARADS2D                                        #
# by dominik.bueschgens                                        #
# March 2021                                                   #
################################################################

abstract type Pair2D <: AbstractVector{Real} end

struct Point2D{T<:AbstractFloat} <: Pair2D
    x::T
    y::T
end

Point2D(x::Real, y::Real) = Point2D(promote(x,y)...)
Point2D(x::Integer, y::Integer) = Point2D(convert(Float64,x), convert(Float64,y))

struct Index2D{T<:Integer} <: Pair2D
    x::T
    y::T
end

Pair2D(x::T,y::T) where T <: AbstractFloat = Point2D(x,y)
Pair2D(x::Integer,y::T) where T <: AbstractFloat = Point2D(x,y)
Pair2D(x::T,y::Integer) where T <: AbstractFloat = Point2D(x,y)
Pair2D(x::T,y::T) where T <: Integer = Index2D(x,y)

import Base: size, getindex, setindex!, +, -, *, /

size(v::Pair2D) = (2, )
getindex(v::Pair2D, i::Int) = i < 3 ? (i == 1 ? v.x : v.y) : error("max index = 2")

function setindex!(v::Pair2D, a::Real, i::Integer)
    if i < 3
        if i == 1
            v.x = a
        else
            v.y = a
        end
    else
        error("max index = 2")
    end
    return v
end

+(a::T, b::T) where T <: Pair2D = T(a.x+b.x,a.y+b.y)
-(a::T, b::T) where T <: Pair2D = T(a.x-b.x,a.y-b.y)
*(a::T1, b::T2) where {T1 <: Pair2D, T2 <: Real} = T1(a.x*b,a.y*b)
/(a::T1, b::T2) where {T1 <: Pair2D, T2 <: Real} = T1(a.x/b,a.y/b)


text:
module PABASTO
### Mailboxes ###
import Base: put!, wait, isready, take!, fetch
using IntervalTrees
using Parameters

type IntervalMailbox <: AbstractChannel
	data::Dict
	lock::ReentrantLock
	function IntervalMailbox()
		new(Dict(),ReentrantLock())
	end
end

function put!{T}(m::IntervalMailbox,v::T,low::Real,high::Real)
	lock(m.lock)
	if !haskey(m.data,T)
		m.data[T]=IntervalMap{Real,T}()
	end
	push!(m.data[T],IntervalValue{Real,T}(low,high,v))
	unlock(m.lock)
	return m
end

function take!(m::IntervalMailbox,position)
	lock(m.lock)
	n=length(values(m.data))
	all_channels=shuffle(collect(m.data))
	for (T,channel) in all_channels
		i=IntervalTrees.firstintersection(channel,Interval{Real}(position-0.1,position+0.1))

		if i.index!=0
			leafnode=i.node
			interval=leafnode.entries.data[i.index]
			c=channel.n
			IntervalTrees.deletefirst!(channel,Interval{Real}(interval.first,interval.last))
			@assert channel.n==c-1
			unlock(m.lock)
			return interval.value
		end
	end
	unlock(m.lock)
	return nothing
end

type Mailbox <: AbstractChannel
	data::Dict
	function Mailbox()
		new(Dict())
	end
end

function put!{T}(m::Mailbox, v::T)
	if !haskey(m.data,T)
		m.data[T]=Channel(100)
	end
	put!(m.data[T],v)
	return m
end

function take!(m::Mailbox)
	n=length(values(m.data))
	all_channels=filter(x->isready(x[2]),collect(m.data))
	if length(all_channels)>0
		priorities=map(x->priority(x[1],length(x[2].data)),all_channels)
		mtype,c=all_channels[indmax(priorities)]
		return take!(c)
	else
		return nothing
	end
end

### Parameters ###
abstract Parameter
abstract Gradient
type ConcreteParameter <: Parameter
end

type ConcreteGradient <: Gradient
end

### Messages ###

type GradientUpdateMessage
	gradient::Future
end

function GradientUpdateMessage(gradient::Gradient)
	GradientUpdateMessage(@spawnat myid() gradient)
end

type ParameterUpdateRequestMessage
	worker_mailbox::RemoteChannel
end

type ParameterUpdateMessage
	parameters
end

type ParameterRequestMessage
end

type ExampleIndicesMessage
	indices::Array{Int}
end
type TestExampleIndicesMessage
	indices::Array{Int}
end
type TestLossMessage
	loss::Real
end

type ExamplesRequestMessage
	id::Int
	worker_mailbox::RemoteChannel
end

type CeaseOperationMessage
end

type CeasedOperationMessage
	id::Int
end

type FinishOperationMessage
end

type FinishedOperationMessage
	id::Int
end

type AdaptiveControlPolicyMessage
	tau::Float64
	num_workers::Int
	num_paramservers::Int
	example_batch_size::Int
	gossip_time::Float64
end

function update(p::Parameter, g::Gradient)
	#update p with gradient g
end

function add_procs(count)
	return fetch(@spawnat 1 Main.add_pabasto_procs(count))
end

##Priorities##
#The default priority of a message type is just the number
#of messages of that type

#To define priority of a new type, follow the example of
#ParameterUpdateMessage
function priority(message_type,queue_length)
	return queue_length
end
function priority(x::Type{TestExampleIndicesMessage},queue_length)
	return 5*queue_length
end

function priority(x::Type{ParameterUpdateMessage},queue_length)
	return 10*queue_length
end

function priority(x::Type{FinishOperationMessage},queue_length)
	return -1
end

# todo: split into separate modules
include("master.jl")
include("paramserver.jl")
include("worker.jl")
end


text:
module NSDEBase

############################################################################################
#                                         EXPORTS                                          #
############################################################################################

export AbstractNSDEProblem
export AbstractNSDESolver
export AbstractNSDESolution

export AbstractInitialValueProblem
export AbstractInitialValueSolver
export AbstractInitialValueSolution
export AbstractRightHandSideFunction

export InitialValueProblem, IVP
export RightHandSideFunction, RHS
export SplitRightHandSideFunction, SRHS

export Dahlquist
export Logistic
export SimplePendulum
export VanderPol
export Rössler
export Lorenz
export Lorenz96

export solve, solve!

############################################################################################
#                                           CORE                                           #
############################################################################################

using ArrowMacros
using LinearAlgebra
using ForwardDiff
using FiniteDifferences

include("abstract.jl")
include("utils.jl")
include("rhs.jl")
include("srhs.jl")
include("ivp.jl")
include("odes.jl")

end


text:
using DBInterface, Test

@test_throws MethodError DBInterface.connect(Int64)

# test @sql_str macro (does nothing)
@test sql"SELECT * FROM MyTable" == "SELECT * FROM MyTable"


text:
# # Chapter 1

# ## Using Types
using BenchmarkTools

# Declare type of function argument
iam(x::Integer) = "an integer"
iam(x::String) = "a string"

function addme(a, b)
  #Declare type of local variable x
  x::Int64 = 2
  #Type of variable y will be inferred
  y = (a+b) / x
  return y
end

iam(1)

iam("1")

iam(1.5)

# ## Multiple Dispatch
sumsqr(x, y) = x^2 + y^2

sumsqr(1, 2)

sumsqr(1.5, 2.5)

sumsqr(1 + 2im , 2 + 3im)

sumsqr(2 + 2im, 2.5)

#Composite Types

struct Pixel
    x::Int64
    y::Int64
    color::Int64
end

p = Pixel(5,5, 100)

p.x = 10;

p.x

# immutable types

mutable struct MPixel
    x::Int64
    y::Int64
    color::Int64
end

p = MPixel(5,5, 100)

p.x=10;

p.x

# TYpe Parameters

struct PPixel{T}
    x::Int64
    y::Int64
    color::T
end

#Type Inference

[x for x=1:5]

# Type Stability

function pos(x)
   if x < 0
      return 0
   else
      return x
   end
end

pos(-1)

pos(-2.5)

pos(2.5)

typeof(pos(2.5))

typeof(pos(-2.5))

function pos_fixed(x)
    if x < 0
        return zero(x)
    else
        return x
    end
end

pos_fixed(-2.4)

pos_fixed(-2)

typeof(pos_fixed(-2.4))

typeof(pos_fixed(-2))

@btime pos(2.5)
@btime pos_fixed(2.5)

@code_warntype pos(2.5)
@code_warntype pos_fixed(2.5)
@code_llvm pos(2.5)
@code_llvm pos_fixed(2.5)

@code_native pos(2.5)
@code_native pos_fixed(2.5)

# Loop variables

function sumsqrtn(n)
    r = 0
    for i = 1:n
        r = r + sqrt(i)
    end
    return r
end

@code_warntype sumsqrtn(5)

function sumsqrtn_fixed(n)
     r = 0.0
     for i = 1:n
         r = r + sqrt(i)
     end
     return r
end

@code_warntype sumsqrtn_fixed(5)

@btime sumsqrtn(1000_000)

@btime sumsqrtn_fixed(1000_000)

#Kernel Methods

function string_zeros(s::AbstractString)
    n=1000_000
    x = s=="Int64" ?
        Vector{Int64}(undef,n) :
        Vector{Float64}(undef, n)
    for i in 1:length(x)
        x[i] = 0
    end
    return x
end

@btime string_zeros("Int64");

function string_zeros_stable(s::AbstractString)
    n = 1000_000
    x = s=="Int64" ?
        Vector{Int64}(undef, n) :
        Vector{Float64}(undef, n)
    return fill_zeros(x)
end

function fill_zeros(x)
    for i in 1:length(x)
        x[i] = 0
    end
    return x
end

@btime string_zeros_stable("Int64");

#Types in storage
## Arrays

a = Int64[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
b = Number[1,2,3,4,5,6,7,8,9,10]

function arr_sumsqr(x::Array{T}) where T <: Number
    r = zero(T)
    for i = 1:length(x)
        r = r + x[i] ^ 2
    end
    return r
end

 @btime arr_sumsqr(a)
 @btime arr_sumsqr(b)

## Composite Types
struct Point
    x
    y
end

struct ConcretePoint
    x::Float64
    y::Float64
end

function sumsqr_points(a)
    s=0.0
    for x in a
        s = s + x.x^2 + x.y^2
    end
    return s
end

p_array = [Point(rand(), rand()) for i in 1:1000_000];
cp_array = [ConcretePoint(rand(), rand()) for i in 1:1000_000];

@btime sumsqr_points(p_array)
@btime sumsqr_points(cp_array)

## Parametric Composite Types

struct PointWithAbstract
    x::AbstractFloat
    y::AbstractFloat
end

struct ParametricPoint{T <: AbstractFloat}
    x::T
    y::T
end

pp_array = [ParametricPoint(rand(), rand()) for i in 1:1000_000];

@btime sumsqr_points(pp_array)


text:
using Test
using LIKWID

const MARKERFILE = "markerfile.txt"

Marker.init()
@test LIKWID.MarkerFile.read(MARKERFILE)
@test LIKWID.MarkerFile.numregions() == 2
@test LIKWID.MarkerFile.regiontag(0) == "matmul"
@test LIKWID.MarkerFile.regiontag(1) == "eigen"
@test LIKWID.MarkerFile.regiongroup(0) == 0
@test LIKWID.MarkerFile.regionevents(0) == 7
@test LIKWID.MarkerFile.regionmetrics(0) == 10
@test LIKWID.MarkerFile.regionthreads(0) == 1
@test LIKWID.MarkerFile.regiontime(0,0) ≈ 14.57989
@test LIKWID.MarkerFile.regiontime(1,0) ≈ 10.75343
@test LIKWID.MarkerFile.regioncount(0,0) == 1
@test LIKWID.MarkerFile.regionresult(0,0,0) ≈ 3.341458e9
@test LIKWID.MarkerFile.regionmetric(0,0,0) ≈ 14.57989
Marker.close()


text:
# AES tests from:
# https://git.lysator.liu.se/nettle/nettle/blob/master/testsuite/aes-test.c

# AES 128
for (key,text,encrypted) in [
    (
        "00010203050607080A0B0C0D0F101112",
        "506812A45F08C889B97F5980038B8359",
        "D8F532538289EF7D06B506A4FD5BE9C9"
    ),(
        "14151617191A1B1C1E1F202123242526",
        "5C6D71CA30DE8B8B00549984D2EC7D4B",
        "59AB30F4D4EE6E4FF9907EF65B1FB68C"
    ),(
        "28292A2B2D2E2F30323334353738393A",
        "53F3F4C64F8616E4E7C56199F48F21F6",
        "BF1ED2FCB2AF3FD41443B56D85025CB1",
    ),(
        "A0A1A2A3A5A6A7A8AAABACADAFB0B1B2",
        "F5F4F7F684878689A6A7A0A1D2CDCCCF",
        "CE52AF650D088CA559425223F4D32694",
    ),(
        "2b7e151628aed2a6abf7158809cf4f3c",
        "6bc1bee22e409f96e93d7e117393172aae2d8a571e03ac9c9eb76fac45af8e5130c81c46a35ce411e5fbc1191a0a52eff69f2445df4f9b17ad2b417be66c3710",
        "3ad77bb40d7a3660a89ecaf32466ef97f5d3d58503b9699de785895a96fdbaaf43b1cd7f598ece23881b00e3ed0306887b0c785e27e8ad3f8223207104725dd4",
    )
]
    @test encrypt("aes128", hex2bytes(key), hex2bytes(text)) == hex2bytes(encrypted)
    @test decrypt("aes128", hex2bytes(key), hex2bytes(encrypted)) == hex2bytes(text)
end

# AES192
for (key,text,encrypted) in [
    (
        "00010203050607080A0B0C0D0F10111214151617191A1B1C",
        "2D33EEF2C0430A8A9EBF45E809C40BB6",
        "DFF4945E0336DF4C1C56BC700EFF837F",
    ),(
        "8e73b0f7da0e6452c810f32b809079e562f8ead2522c6b7b",
        "6bc1bee22e409f96e93d7e117393172aae2d8a571e03ac9c9eb76fac45af8e5130c81c46a35ce411e5fbc1191a0a52eff69f2445df4f9b17ad2b417be66c3710",
        "bd334f1d6e45f25ff712a214571fa5cc974104846d0ad3ad7734ecb3ecee4eefef7afd2270e2e60adce0ba2face6444e9a4b41ba738d6c72fb16691603c18e0e",
    )
]
    @test encrypt("aes192", hex2bytes(key), hex2bytes(text)) == hex2bytes(encrypted)
    @test decrypt("aes192", hex2bytes(key), hex2bytes(encrypted)) == hex2bytes(text)
end

# AES256
for (key,text,encrypted) in [
    (
        "00010203050607080A0B0C0D0F10111214151617191A1B1C1E1F202123242526",
        "834EADFCCAC7E1B30664B1ABA44815AB",
        "1946DABF6A03A2A2C3D0B05080AED6FC",
    ),(
        "603deb1015ca71be2b73aef0857d77811f352c073b6108d72d9810a30914dff4",
        "6bc1bee22e409f96e93d7e117393172aae2d8a571e03ac9c9eb76fac45af8e5130c81c46a35ce411e5fbc1191a0a52eff69f2445df4f9b17ad2b417be66c3710",
        "f3eed1bdb5d2a03c064b5a7e3db181f8591ccb10d410ed26dc5ba74a31362870b6ed21b99ca6f4f9f153e7b1beafed1d23304b7a39f9f3ff067d8d8f9e24ecc7",
    )
]
    @test encrypt("aes256", hex2bytes(key), hex2bytes(text)) == hex2bytes(encrypted)
    @test decrypt("aes256", hex2bytes(key), hex2bytes(encrypted)) == hex2bytes(text)
end

# Test lower-level API
key = "this key's exactly 32 bytes long"
enc = Encryptor("AES256", key)
plaintext = "this is 16 chars"
ciphertext = encrypt(enc, plaintext)

dec = Decryptor("AES256", key)
deciphertext = decrypt(dec, ciphertext)
plaintext == bytestring(deciphertext)

# Test errors
@test_throws ArgumentError Encryptor("this is not a cipher", key)
@test_throws ArgumentError Decryptor("this is not a cipher either", key)

@test_throws ArgumentError Encryptor("AES256", "bad key length")
@test_throws ArgumentError Decryptor("AES256", "bad key length")

enc = Encryptor("AES256", key)
dec = Decryptor("AES256", key)
result = Array(UInt8, length(plaintext) - 1)
@test_throws ArgumentError encrypt!(enc, result, plaintext)
@test_throws ArgumentError decrypt!(dec, result, plaintext)

# Test show methods
println("Testing cipher show methods:")
println(get_cipher_types()["AES256"])
println(Encryptor("AES256", key))
println(Decryptor("AES256", key))


text:
"""
Toy program to demonstrate/practice parallel programming techniques.
"""
module GameOfLife

include("strategies.jl")
include("simulation.jl")

export step!, prepare!, life_rule!
export Serial, ThreadParallel, ProcParallel, LocalTasks, DistributedTasks

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Shaderc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Shaderc")
JLLWrappers.@generate_main_file("Shaderc", UUID("6bef4a77-805f-5404-b838-981e48d6e429"))
end  # module Shaderc_jll


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "VerizonEdgecastTokenBuilder"
version = v"0.1.0"

# Collection of sources required to build VerizonEdgecastTokenBuilder
sources = [
    "https://github.com/VerizonDigital/ectoken.git" =>
    "7b8812d476f5be5b290fe2832859b9b7636f43ae",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
tar -xzf *.tar.gz
export OPENSSL_INCLUDE="-I $WORKSPACE/srcdir/include"
export OPENSSL_LIBS="$WORKSPACE/srcdir/lib/libssl.a $WORKSPACE/srcdir/lib/libcrypto.a"
cd ectoken/c-ectoken/ecencrypt/
gcc -m64 -O2 -Wall -Werror -std=gnu99 ec_encrypt.c ectoken_v3.c base64.c -o 64/ectoken3 $OPENSSL_LIBS $OPENSSL_INCLUDE -lm
cp 64/ectoken3 $prefix/bin/
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:musl),
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "ectoken3", :ectoken3)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/OpenSSL-v1.1.1%2Bc%2B0/build_OpenSSL.v1.1.1+c.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: tanmaykm/VerizonEdgecastTokenBuilder

max_issues_repo_name: tanmaykm/VerizonEdgecastTokenBuilder

max_forks_repo_name: tanmaykm/VerizonEdgecastTokenBuilder


text:
using DeterministicPolicyGradient, ControlSystems
using Base.Test

const G = tf(0.1,[1,0.01,0.1])

const T = 20
const t = 1:T
const m = 1
const n = 2
const p = 5
x0 = [0,1.]

# Cost function             ==========================================
const Q1 = diagm([5,1])
const Q2 = 1eye(m)
# r(s,a) = -sum(abs(Q1*s)) - sum(abs(Q2*a))
r(s,a,t) = (-(s)'Q1*(s) - a'Q2*a)[1]
L = lqr(ss(G),Q1,Q2)
μlqr(x) = -L*x

# Initialize functions      ==========================================
cp = linspace(-5,5,p)
cv = linspace(-5,5,p)
grid1 = meshgrid(cp,cv)
const c1 = [grid1[1][:] grid1[2][:]]
P = size(c1,1)

function ϕ(s)
    a = exp(-1/(2*2)*sum((s'.-c1).^2,2))#exp(-1/2*(s.-c).^2)
    a ./= sum(a)
    a
    # [a; s[2]]
end

μ(s,Θ,t)          = Θ'ϕ(s)
∇μ(s)           = ϕ(s)
β(s,Θ,noise,i)  = Θ'ϕ(s) + noise[i]
ϕ(s,a,Θ)        = ∇μ(s)*(a-μ(s,Θ,t))
V(s,v)          = v'ϕ(s) # It's a good idea to bias V to some mean of the final landscape
Q(s,a,v,w,Θ,t,c)    = (ϕ(s,a,Θ)'w + V(s,v))[1]

function gradients(s1,s,a1,a,Θ,w,v,t,c)
    ∇μ = ϕ(s)
    ∇aQ = ∇μ'w
    ∇wQ = ∇μ*(a-Θ'∇μ)
    ∇vQ = ∇μ
    ∇aQ, ∇wQ, ∇vQ, ∇μ
end

simulate(Θ,x0, noise) = lsim(G, (t,s)->β(s,Θ,noise,t), t, x0)[3:4]
simulate(Θ,x0) = lsim(G, (t,s)->μ(s,Θ,t), t, x0)[3:4]
exploration(σβ) = filt(ones(5),[5],σβ*randn(T))
funs            = DPGfuns(μ,Q, gradients, simulate, exploration, r)

Θ               = zeros(P*m) # Weights
w               = 0.00001randn(P*m)
v               = 0.00001randn(P)
initial_state   = DPGstate(Θ,w,v)

opts = DPGopts(m,
σβ                  = 1,
αΘ                  = 0.001,
γ                   = 0.9999,
τ                   = 0.01,
iters               = 250,
critic_update       = :rls,
λrls                = 0.99999,
stepreduce_interval = 1,
stepreduce_factor   = 0.995,
hold_actor          = 10,
experience_replay   = 0,
experience_ratio    = 10,
momentum            = 0.8,
rmsprop             = true,
rmspropfactor       = 0.9,
critic_update_interval = 3,
eval_interval       = 10,
divergence_threshold  = 1.5)

cost, Θ, w, v, mem = dpg(opts, funs, initial_state, x0, c1)
@test minimum(cost) < 0.4cost[2]

opts = DPGopts(m,
σβ                  = 2,
αΘ                  = 0.001,
γ                   = 0.999,
τ                   = 0.01,
iters               = 250,
critic_update       = :kalman,
λrls                = 0.99999,
stepreduce_interval = 1,
stepreduce_factor   = 0.995,
hold_actor          = 10,
experience_replay   = 0,
experience_ratio    = 10,
momentum            = 0.8,
rmsprop             = true,
rmspropfactor       = 0.9,
critic_update_interval = 1,
eval_interval       = 10,
divergence_threshold  = 1.5)

cost, Θ, w, v, mem = dpg(opts, funs, initial_state, x0, c1)
@test minimum(cost) < 0.4cost[2]

opts = DPGopts(m,
σβ                  = 2,
αΘ                  = 0.001,
αw                  = 0.01,
αv                  = 0.01,
γ                   = 0.99,
τ                   = 0.01,
iters               = 3000,
critic_update       = :gradient,
stepreduce_interval = 10,
stepreduce_factor   = 0.995,
hold_actor          = 500,
experience_replay   = 0,
experience_ratio    = 1,
momentum            = 0.8,
rmsprop             = true,
rmspropfactor       = 0.9,
eval_interval       = 100,
divergence_threshold  = 1.5)
# Det var något vajsing med att δ aldrig blev uträknad vid gradient training
cost, Θ, w, v, mem = dpg(opts, funs, initial_state, x0, c1)
@test minimum(cost) < 0.5cost[2]

mem

## Test experience replay
n = 2
s = randn(n)
s1 = randn(n)
a = randn(n)
reward = randn()
δ1 = 1.
δ2 = 2.
δ3 = 3.

t1 = Transition(s,s1,a,reward,δ1,1,1)
t2 = Transition(s,s1,a,reward,δ2,2,2)
t3 = Transition(s,s1,a,reward,δ3,3,3)

@test t1 < t2

N = 3
mem = SortedReplayMemory(N)

push!(mem,t2)
push!(mem,t1)
push!(mem,t3)
@test length(mem.mem) == N
@test Collections.isheap(mem.mem, Base.Order.Reverse)
sort!(mem)
@test Collections.isheap(mem.mem, Base.Order.Reverse)

tt = sample_greedy!(mem)
@test tt.δ == 3.


text:
function test_user_knn()
    println("-- Testing UserKNN rating prediction")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    recommender = UserKNN(data, 1)
    build!(recommender)

    @test isapprox(recommender.sim[1, 2], 0.447, atol=1e-3)
    @test isapprox(recommender.sim[2, 3], -0.693, atol=1e-3)

    rec = recommend(recommender, 1, 3, [i for i in 1:8])
    @test rec[1] == (5 => 3.0)
    @test rec[2] == (8 => 3.0)
    @test rec[3] == (2 => 2.0)

    recommender = UserKNN(data, 1, true)
    build!(recommender)

    @test isapprox(recommender.sim[1, 2], 0.447, atol=1e-3)
    @test isapprox(recommender.sim[2, 3], -0.693, atol=1e-3)

    rec = recommend(recommender, 1, 3, [i for i in 1:8])
    @test rec[1] == (5 => 3.0)
    @test rec[2] == (8 => 3.0)
    @test rec[3] == (7 => 2.2)
end

test_user_knn()


text:
module Day12

export get_inputs, get_solution1, get_solution2

## Input getting
function get_inputs()
    test_input1 = test_input2 = read_inputs("test_input1.txt")
    test_output1 = 25
    test_output2 = 286
    data = read_inputs("input.txt")
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

function read_inputs(f_name)
    lines = readlines(joinpath(@__DIR__, f_name))
    return [Symbol(line[1])=>parse(Int, line[2:end]) for line in lines]
end

## Solution functions
# instruction => (trans, rot)
const move_dict = Dict(
    :N => ( im,   0),
    :S => (-im,   0),
    :E => (  1,   0),
    :W => ( -1,   0),
    :L => (  0,  im),
    :R => (  0, -im),
)

Base.@kwdef mutable struct Ship
    pos::Complex{Int} = 0
    dir::Complex{Int} = 1
    waypoint::Complex{Int} = 10+im
end

function get_solution(data, fun=move!)
    ship = Ship()
    for line in data
        fun(ship, line)
    end
    pos = ship.pos
    return abs(pos.re) + abs(pos.im)
end

# Part 1
function move!(ship, (instruction, amount))
    if instruction==:F
        ship.pos += amount * ship.dir
    else
        trans, rot = move_dict[instruction]
        ship.pos += amount * trans
        ship.dir *= rot ^ (amount ÷ 90)
    end
    return ship
end

get_solution1(data) = get_solution(data)

# Part 2
function move_with_waypoint!(ship, (instruction, amount))
    if instruction==:F
        ship.pos += ship.waypoint * amount
    else
        trans, rot = move_dict[instruction]
        ship.waypoint += amount * trans
        ship.waypoint *= rot ^ (amount ÷ 90)
    end
    return ship
end

get_solution2(data) = get_solution(data, move_with_waypoint!)

end


text:
module CellValuesGalleryTests

using Test
using Gridap
using Gridap.CellValuesGallery
using TensorValues
using StaticArrays

v = [ VectorValue(i,i) for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_number( cv, v )

v = [ SVector(i,i) for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_number( cv, v )

v = CartesianIndices((3,10))
cv = CellValueFromArray(v)
test_index_cell_value( cv, v )

v = [ [i,i,i] for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_array( cv, v )

v = [ [i,i,i] for i in 1:10 ]
cv = CellValueFromArray(v)
test_index_cell_array( cv, v )

data = [2,3,1,3,6,7,3,2,5,6,3,4]
ptrs = [1,4,4,7,13]
ca = CellVectorFromDataAndPtrs(data,ptrs)
a = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
test_index_cell_array( ca, a )

d,p = compress(ca)
@test d === ca.data
@test p === ca.ptrs

a = [ data[3*(i-1)+1:3*i] for i in 1:4]
ca = CellVectorFromDataAndStride(data,3)
test_index_cell_array( ca, a )

p = VectorValue(2.0,1.0)
data = [2,3,1,3,4,4,3,2,5,4,3,4]
ptrs = [1,4,4,7,13]
gid_to_val = [p, 2*p, 3*p, -p, p]
lid_to_gid = CellVectorFromDataAndPtrs(data,ptrs)
ca = CellVectorFromLocalToGlobal(lid_to_gid,gid_to_val)
a = [ gid_to_val[gids] for gids in lid_to_gid ]
test_index_cell_array( ca, a )

struct CartesianArray{T,N,A<:AbstractArray{T,N}} <:AbstractArray{T,N}
  data::A
end

Base.size(a::CartesianArray) = size(a.data)

function Base.getindex(
  a::CartesianArray{T,N},I::Vararg{Int,N}) where {T,N}
  a.data[I...]
end

b = fill([2,4],(3,4))
c = CartesianArray(b)
lid_to_gid = CellValueFromArray(c)
test_index_cell_array( lid_to_gid, c )
gid_to_val = [1, 2, 3, -1, 2]
ca = CellVectorFromLocalToGlobal(lid_to_gid,gid_to_val)
aa = fill([2,-1],(3,4))
a = CartesianArray(aa)
test_index_cell_array( ca, a )

data = [-2,3,-1,3,4,-4,-3,2,5,4,3,-4]
ptrs = [1,4,4,7,13]
p = 1
gid_to_val_pos = [p, 2*p, 3*p, -p, p]
gid_to_val_neg = [2*p, p, 5*p, 3*p, p]
lid_to_gid = CellVectorFromDataAndPtrs(data,ptrs)
ca = CellVectorFromLocalToGlobalPosAndNeg(
  lid_to_gid,gid_to_val_pos,gid_to_val_neg)
a = [ [1, 3, 2], Int64[], [3, -1, 3], [5, 2, 1, -1, 3, 3] ]
test_index_cell_array( ca, a )

b = fill([2,-4],(3,4))
c = CartesianArray(b)
lid_to_gid = CellValueFromArray(c)
test_index_cell_array( lid_to_gid, c )
gid_to_val_pos = [1, 2, 3, -1, 2]
gid_to_val_neg = [-1, -3, -3, 1, -4]
ca = CellVectorFromLocalToGlobalPosAndNeg(
  lid_to_gid,gid_to_val_pos,gid_to_val_neg)
aa = fill([2,1],(3,4))
a = CartesianArray(aa)
test_index_cell_array( ca, a )

cell_to_x_l = [2,3,1,3,4,4,3,2,5,4,3,4]
cell_to_x_p = [1,4,4,7,13]
cell_to_x = CellVectorFromDataAndPtrs(cell_to_x_l,cell_to_x_p)
x_to_vals_l = [5,4,1,2,3,6,7,8,9,10]
x_to_vals_p = [1,3,5,6,10,11]
x_to_vals = CellVectorFromDataAndPtrs(x_to_vals_l,x_to_vals_p)
cell_to_vals = CellVectorByComposition(cell_to_x, x_to_vals)
a = [
  [1, 2, 3, 5, 4],
  Int[],
  [3, 6, 7, 8, 9, 6, 7, 8, 9],
  [3, 1, 2, 10, 6, 7, 8, 9, 3, 6, 7, 8, 9]]
test_index_cell_array( cell_to_vals, a )

end # module