On Sat, 19 Apr 2008, Chris Burke wrote:
See also:
open '~system/packages/math/quatern.ijs'
I just got around to looking at the quaternion routines
in j602/system/packages/math/quatern.ijs. By representing
the quaternions as 4x4 matrices, they are quite compact.
I wondered if I should use them. The verb for multiplying
quaternions given there is just
mp=. +/ . *
j=. ,: =i.4
j=. j, 0 1 0 0, _1 0 0 0, 0 0 0 _1,: 0 0 1 0
j=. j, 0 0 1 0, 0 0 0 1, _1 0 0 0,: 0 _1 0 0
j=. j, 0 0 0 1, 0 0 _1 0, 0 1 0 0 ,: _1 0 0 0
qnmat=. mp & j
qnmul=: (mp qnmat) f. "1
An array of quaternions is stored as a (n 4) dimension
element. I thought I should compare this with what I
had used for this mode of storage, which was
dot=: [: +/"1 * NB. dot products of vector arrays
pm1=: 1&|."1
pm2=: 2&|."1
pm12=: 4 : '(pm1 x)*pm2 y'
cross=: pm12 - pm12~ NB. cross products of vectors
NB. Multiply quaternions: q1 QM q2 --> q3
QM=: 4 : 0
s1=. 0{"1 x
v1=. }."1 x
s2=. 0{"1 y
v2=. }."1 y
s=. (s1*s2) - v1 dot v2
v=. (s1*v2)+(s2*v1)+(v1 cross v2)
s ,&.|: v NB. catenate under transpose
)
Since I'm interested in the speed of the operations for large
arrays, let's compare the verbs:
rand=: [: ? ] # 0"_
q1=. 200000 4 $ rand 800000
q2=. 200000 4 $ rand 800000
ts'q=. q1 qnmul q2'
0.313453 8.39142e6
ts'qq=. a QM b'
0.195947 6.50142e7
The output is the same:
12345{q
_0.949615 _0.673075 0.278588 0.426178
12345{qq
_0.949615 _0.673075 0.278588 0.426178
but qnmul is about 50% slower than QM, though it uses much
less space.
But what I normally do is store the quaternions as (4 n)
arrays, as I've found that faster. Thus I use
dot2=: [: +/ *
cross2=: ((1: |.[)*(_1: |. ]))-((_1: |.[)*(1:|.]))
QM2=: 4 : 0
s1=. 0{ x
v1=. }. x
s2=. 0{ y
v2=. }. y
s=. (s1*s2) - v1 dot2 v2
s,(s1*"1 v2)+(s2*"1 v1)+(v1 cross2 v2)
)
Comparing,
q1r=. |: q1
q2r=. |: q2
` ts'qq2=. q1r QM2 q2r'
0.139955 7.34037e7
12345{"1 qq2
_0.949615 _0.673075 0.278588 0.426178
Which is more than twice as fast as qnmul (at the expense
of 10 times the space). Could qnmul be tweaked to get
better speed on large arrays?
Patrick
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
