JuliaGeometry / Quaternions.jl

A Julia implementation of quaternions
https://juliageometry.github.io/Quaternions.jl
MIT License
116 stars 37 forks source link

Type instability in `axis` #109

Closed mikmoore closed 1 year ago

mikmoore commented 1 year ago

Return type is always Vector{Float64} when angle(q) is zero or NaN, leading to a type instability.

julia> code_warntype(axis,(QuaternionF32,))
MethodInstance for Quaternions.axis(::QuaternionF32)
  from axis(q::Quaternion) in Quaternions at /home/mike/.julia/packages/Quaternions/GVtxM/src/Quaternion.jl:205
Arguments
  #self#::Core.Const(Quaternions.axis)
  q@_2::QuaternionF32
Locals
  s::Float32
  q@_4::QuaternionF32
Body::Union{Vector{Float32}, Vector{Float64}}
1 ─       (q@_4 = q@_2)
│         (q@_4 = Quaternions.normalize(q@_4))
│   %3  = Quaternions.angle(q@_4)::Float32
│   %4  = (%3 / 2)::Float32
│         (s = Quaternions.sin(%4))
│   %6  = Quaternions.abs(s)::Float32
│   %7  = (%6 > 0)::Bool
└──       goto #3 if not %7
2 ─ %9  = Base.getproperty(q@_4, :v1)::Float32
│   %10 = Base.getproperty(q@_4, :v2)::Float32
│   %11 = Base.getproperty(q@_4, :v3)::Float32
│   %12 = Base.vect(%9, %10, %11)::Vector{Float32}
│   %13 = (%12 / s)::Vector{Float32}
└──       return %13
3 ─ %15 = Base.vect(1.0, 0.0, 0.0)::Vector{Float64}
└──       return %15

I'd have made a PR but am not at a machine where I can do so right now. I think something like

function axis(q::Quaternion)
    s = sin(angle(q) / 2) * abs(q)
    iszero(s) ?
        [one(q.v1/s), zero(q.v2/s), zero(q.v3/s)] : 
        [q.v1/s, q.v2/s, q.v3/s]
end

will fix the instability and also allow NaN propagation. This still doesn't fix the case where only one imaginary component isinf, where we could give a reasonable answer (assuming we care to have axis defined for infinite quaternions), but that was broken under the old version too.

EDIT: Alternatively,

function axis(q::Quaternion)
    s = abs_imag(q)
    iszero(s) ?
        [one(q.v1/s), zero(q.v2/s), zero(q.v3/s)] : 
        [q.v1/s, q.v2/s, q.v3/s]
end

would be more accurate, propagate NaN on only the imaginary part, and would return a well-defined vector for the zero quaternion. If we don't want zero (or any not-convincingly unit quaternion?) to be a valid input, then I'd suggest a DomainError.

hyrodium commented 1 year ago

I'm planning to remove the function in https://github.com/JuliaGeometry/Quaternions.jl/issues/86. Do you have any thoughts on this?

hyrodium commented 1 year ago

The function was removed in https://github.com/JuliaGeometry/Quaternions.jl/pull/111.