Simplify v*q, make q*q more accurate, add q*q and q/q unit tests

vector * quaternion: Simplified by precalculating the repeated products and removing math.fsum.

quaternion * quaternion: Add more F32's to get a result closer to SL's (still not there, but we're definitively closer now).

We left out quaternion/quaternion in the tests, and the q*q test was not general enough (had many zeros). Remedied that.
This commit is contained in:
Sei Lisa 2019-05-23 01:48:13 +02:00
parent ec509b3840
commit 00b103c8aa
3 changed files with 34 additions and 16 deletions

View file

@ -722,10 +722,11 @@ def mul(a, b, f32=True):
a = q2f(a)
b = q2f(b)
# quaternion product - product formula reversed
return Quaternion(F32((a[0] * b[3] + a[3] * b[0] + a[2] * b[1] - a[1] * b[2],
a[1] * b[3] - a[2] * b[0] + a[3] * b[1] + a[0] * b[2],
a[2] * b[3] + a[1] * b[0] - a[0] * b[1] + a[3] * b[2],
a[3] * b[3] - a[0] * b[0] - a[1] * b[1] - a[2] * b[2]), f32))
return Quaternion(F32((F32(a[0] * b[3]) + F32(a[3] * b[0]) + F32(a[2] * b[1]) - F32(a[1] * b[2]),
F32(a[1] * b[3]) - F32(a[2] * b[0]) + F32(a[3] * b[1]) + F32(a[0] * b[2]),
F32(a[2] * b[3]) + F32(a[1] * b[0]) - F32(a[0] * b[1]) + F32(a[3] * b[2]),
F32(a[3] * b[3]) - F32(a[0] * b[0]) - F32(a[1] * b[1]) - F32(a[2] * b[2])
), f32))
if ta != Vector:
raise ELSLInvalidType # Should never happen at this point
@ -745,21 +746,32 @@ def mul(a, b, f32=True):
raise ELSLInvalidType # Should never happen at this point
# vector * quaternion: perform conjugation
#v = mul(Quaternion((-b[0], -b[1], -b[2], b[3])), mul(Quaternion((a[0], a[1], a[2], 0.0)), b, f32=False))
#v = mul(Quaternion((-b[0], -b[1], -b[2], b[3])),
# mul(Quaternion((a[0], a[1], a[2], 0.0)), b, f32=False))
#return Vector((v[0], v[1], v[2]))
# this is more precise as it goes directly to the gist of it:
# this removes redundant calculations:
a = v2f(a)
b = q2f(b)
return Vector(F32((
math.fsum(( a[0]*(b[0]*b[0]-b[1]*b[1]-b[2]*b[2]+b[3]*b[3]),
a[1]*2*(b[0]*b[1]-b[2]*b[3]),
a[2]*2*(b[0]*b[2]+b[1]*b[3]))),
math.fsum(( a[0]*2*(b[0]*b[1]+b[2]*b[3]),
-a[1]*(b[0]*b[0]-b[1]*b[1]+b[2]*b[2]-b[3]*b[3]), # notice minus sign
a[2]*2*(b[1]*b[2]-b[0]*b[3]))),
math.fsum(( a[0]*2*(b[0]*b[2]-b[1]*b[3]),
a[1]*2*(b[1]*b[2]+b[0]*b[3]),
-a[2]*(b[0]*b[0]+b[1]*b[1]-b[2]*b[2]-b[3]*b[3]))) # notice minus sign
b0b0 = b[0] * b[0]
b0b1 = b[0] * b[1]
b0b2 = b[0] * b[2]
b0b3 = b[0] * b[3]
b1b1 = b[1] * b[1]
b1b2 = b[1] * b[2]
b1b3 = b[1] * b[3]
b2b2 = b[2] * b[2]
b2b3 = b[2] * b[3]
b3b3 = b[3] * b[3]
return Vector(F32(
( a[0] * (b0b0 - b1b1 - b2b2 + b3b3)
+ a[1] * (b0b1 - b2b3) * 2
+ a[2] * (b0b2 + b1b3) * 2
, a[0] * (b0b1 + b2b3) * 2
+ a[1] * (b1b1 - b0b0 - b2b2 + b3b3)
+ a[2] * (b1b2 - b0b3) * 2
, a[0] * (b0b2 - b1b3) * 2
+ a[1] * (b1b2 + b0b3) * 2
+ a[2] * (b2b2 - b0b0 - b1b1 + b3b3)
), f32))
def div(a, b, f32=True):