Add more `geomean` methods for efficient calculation

[hyrodium]

See Note section in #177.
The mean method returns SMatrix (<: AbstractMatrix), but it's better to have a method for mean of rotation which returns Rotation.

The draft of this method is in src/mean.jl. (See also #90)

[hyrodium]

src/mean.jl will be removed in #255. The original file can be viewed at https://github.com/JuliaGeometry/Rotations.jl/blob/v1.4.0/src/mean.jl.

[hyrodium]

We already have StatsBase.geomean that calculates the mean of rotations.

julia> import StatsBase

julia> r1 = QuatRotation(5,1,0,1)
3×3 QuatRotation{Float64} with indices SOneTo(3)×SOneTo(3)(QuaternionF64(0.96225, 0.19245, 0.0, 0.19245)):
 0.925926   -0.37037    0.0740741
 0.37037     0.851852  -0.37037
 0.0740741   0.37037    0.925926

julia> r2 = QuatRotation(6,1,0.2,0)
3×3 QuatRotation{Float64} with indices SOneTo(3)×SOneTo(3)(QuaternionF64(0.985861, 0.16431, 0.032862, 0.0)):
  0.99784    0.0107991   0.0647948
  0.0107991  0.946004   -0.323974
 -0.0647948  0.323974    0.943844

julia> r3 = QuatRotation(8,0,0.1,2)
3×3 QuatRotation{Float64} with indices SOneTo(3)×SOneTo(3)(QuaternionF64(0.970071, 0.0, 0.0121259, 0.242518)):
  0.882076  -0.470519    0.023526
  0.470519   0.88237     0.00588149
 -0.023526   0.00588149  0.999706

julia> StatsBase.geomean([r1,r2,r3])
3×3 RotationVec{Float64} with indices SOneTo(3)×SOneTo(3)(0.239999, 0.0301774, 0.29325):
 0.957068    -0.282661   0.064224
 0.289816     0.929063  -0.229888
 0.00531226   0.238632   0.971096

We can implement more efficient methods into this function as a package extension.