Inconsistent behavior of `jacobian(::Type{RotMatrix}, ::QuatRotation)`?

The method jacobian(::Type{RotMatrix}, ::QuatRotation) is supposed to implement the jacobian of vec(RotMatrix(q)) with respect to the parameters of q, where q isa QuatRotation.
The docstring says:

jacobian(::Type{output_param}, R::input_param)
Returns the jacobian for transforming from the input rotation parameterization
to the output parameterization, centered at the value of R.

Since QuatRotation can be constructed with or without normalization of the input quaternion, it is not entirely clear to me what this means.
I think there are essentially two possibilities:

• Either it is the 9x4 matrix of derivatives of vec(RotMatrix(q)) with respect to the parameters of a unit quaternion, i.e. the jacobian of vec(RotMatrix(QuatRotation(p, false))), where p is a SVector{4} with norm(p) == 1.
• Or it is the 9x4 matrix of derivatives of vec(RotMatrix(q)) with respect to the parameters of a general quaternion, i.e. the jacobian of vec(RotMatrix(QuatRotation(p, true))), where p is a SVector{4} with arbitrary norm. In this case the jacobian should be the product of the jacobian of the above case, times the jacobian of the normalization operation.

Now what jacobian(::Type{output_param}, R::input_param) tries to implement is apparently the second case, i.e. the jacobian including the normalization operation.
However even though it includes the normalization operation, the implementation still assumes that the input quaternion is normalized. I think this is inconsistent (but probably even a bug?)

There is a test for this method but it only tests for the case where the QuatRotation is constructed with an already normalized unit quaternion:

@testset "Jacobian (QuatRotation -> RotMatrix)" begin
    for i in 1:10    # do some repeats
        q = rand(QuatRotation)  # a random quaternion

        # test jacobian to a Rotation matrix
        R_jac = Rotations.jacobian(RotMatrix, q)
        FD_jac = ForwardDiff.jacobian(x -> SVector{9}(QuatRotation(x[1], x[2], x[3], x[4])),
                                      Rotations.params(q))

        # compare
        @test FD_jac ≈ R_jac
    end
end

If one would change line 3 st q = rand(QuatRotation) --> q = QuatRotation(rand(QuaternionF64), false), the test would actually fail.

So in summary I think this method should either assume a unit quaternion and thus exclude the jacobian of the normalization operation, or assume a general (not necessarily unit-) quaternion and include the jacobian of the normalization operation.

One suggestion would be to add a jacobian(::Type{RotMatrix}, ::Quaternion) (note the Quaternion instead of QuatRotation), which assumes a general quaternion and returns the jacobian including the normalization term.
And change the behavior of the existing jacobian(::Type{RotMatrix}, ::QuatRotation) to return the jacobian without the normalization term.
I could provide a PR for that. I've proposed a PR: #291 which relaxes above test in said way and provides the suggested change (as a base for discussion).

Start with test/derivative_tests.jl and the jacobian(::Type{RotMatrix}, ::QuatRotation) entry point described in the issue. Review the proposed semantics and PR #291, then verify the behavior for both normalized and non-normalized quaternion inputs. Done means the API behavior, implementation, and derivative tests consistently reflect the chosen normalization rules.