quat: new package for quaternions

num/quat/abs.go

@@ -0,0 +1,52 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Copyright 2017 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package quat
+
+import "math"
+
+// Abs returns the absolute value (also called the modulus) of q.
+func Abs(q Quat) float64 {
+	// Special cases.
+	switch {
+	case IsInf(q):
+		return math.Inf(1)
+	case IsNaN(q):
+		return math.NaN()
+	}
+
+	r, i, j, k := q.Real, q.Imag, q.Jmag, q.Kmag
+	if r < 0 {
+		r = -r
+	}
+	if i < 0 {
+		i = -i
+	}
+	if j < 0 {
+		j = -j
+	}
+	if k < 0 {
+		k = -k
+	}
+	if r < i {
+		r, i = i, r
+	}
+	if r < j {
+		r, j = j, r
+	}
+	if r < k {
+		r, k = k, r
+	}
+	if r == 0 {
+		return 0
+	}
+	i /= r
+	j /= r
+	k /= r
+	return r * math.Sqrt(1+i*i+j*j+k*k)
+}

num/quat/abs_test.go

@@ -0,0 +1,38 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package quat
+
+import (
+	"math"
+	"testing"
+)
+
+var absTests = []struct {
+	q    Quat
+	want float64
+}{
+	{q: Quat{}, want: 0},
+	{q: NaN(), want: nan},
+	{q: Inf(), want: inf},
+	{q: Quat{Real: 1, Imag: 1, Jmag: 1, Kmag: 1}, want: 2},
+	{q: Quat{Real: -1, Imag: 1, Jmag: -1, Kmag: 1}, want: 2},
+	{q: Quat{Real: 1, Imag: 2, Jmag: 3, Kmag: 4}, want: math.Sqrt(1 + 4 + 9 + 16)},
+	{q: Quat{Real: -1, Imag: -2, Jmag: -3, Kmag: -4}, want: math.Sqrt(1 + 4 + 9 + 16)},
+}
+
+func TestAbs(t *testing.T) {
+	for _, test := range absTests {
+		got := Abs(test.q)
+		if math.IsNaN(test.want) {
+			if !math.IsNaN(got) {
+				t.Errorf("unexpected result for Abs(%v): got:%v want:%v", test.q, got, test.want)
+			}
+			continue
+		}
+		if got != test.want {
+			t.Errorf("unexpected result for Abs(%v): got:%v want:%v", test.q, got, test.want)
+		}
+	}
+}

num/quat/conj.go

@@ -0,0 +1,20 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Copyright 2017 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package quat
+
+// Conj returns the quaternion conjugate of q.
+func Conj(q Quat) Quat {
+	return Quat{Real: q.Real, Imag: -q.Imag, Jmag: -q.Jmag, Kmag: -q.Kmag}
+}
+
+// Inv returns the quaternion inverse of q.
+func Inv(q Quat) Quat {
+	a := Abs(q)
+	return Scale(1/(a*a), Conj(q))
+}

num/quat/conj_test.go

@@ -0,0 +1,41 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package quat
+
+import (
+	"testing"
+
+	"gonum.org/v1/gonum/floats"
+)
+
+var invTests = []struct {
+	q       Quat
+	wantNaN bool
+}{
+	{q: Quat{Real: 1, Imag: 1, Jmag: 1, Kmag: 1}},
+	{q: Quat{Real: 3, Imag: -1, Jmag: 5, Kmag: -40}},
+	{q: Quat{Real: 1e6, Imag: -1e5, Jmag: 4, Kmag: -10}},
+	{q: Quat{Real: 0, Imag: 1, Jmag: 1, Kmag: 1}},
+	{q: Quat{Real: 1, Imag: 0, Jmag: 1, Kmag: 1}},
+	{q: Quat{Real: 1, Imag: 1, Jmag: 0, Kmag: 1}},
+	{q: Quat{Real: 1, Imag: 1, Jmag: 1, Kmag: 0}},
+	{q: Quat{}, wantNaN: true},
+}
+
+func TestInv(t *testing.T) {
+	const tol = 1e-14
+	for _, test := range invTests {
+		got := Mul(test.q, Inv(test.q))
+		if test.wantNaN {
+			if !IsNaN(got) {
+				t.Errorf("unexpected result for Mul(%v, Inv(%[1]v)): got:%v want:%v", test.q, got, NaN())
+			}
+			continue
+		}
+		if !(floats.EqualWithinAbsOrRel(got.Real, 1, tol, tol) && floats.EqualWithinAbsOrRel(Abs(got), 1, tol, tol)) {
+			t.Errorf("unexpected result for Mul(%v, Inv(%[1]v)): got:%v want:%v", test.q, got, Quat{Real: 1})
+		}
+	}
+}

num/quat/doc.go

@@ -0,0 +1,10 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package quat provides the quaternion numeric type and functions.
+//
+// For a good treatment of uses and behaviors of quaternions, see
+// the interactive videos by Ben Eater and Grant Sanderson here
+// https://eater.net/quaternions.
+package quat // imports "gonum.org/v1/gonum/num/quat"

num/quat/exp.go

@@ -0,0 +1,64 @@
+// Copyright ©2018 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Copyright 2017 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package quat
+
+import "math"
+
+// Exp returns e**q, the base-e exponential of q.
+func Exp(q Quat) Quat {
+	w, uv := split(q)
+	if uv == zero {
+		return lift(math.Exp(w))
+	}
+	v := Abs(uv)
+	e := math.Exp(w)
+	s, c := math.Sincos(v)
+	return join(e*c, Scale(e*s/v, uv))
+}
+
+// Log returns the natural logarithm of q.
+func Log(q Quat) Quat {
+	w, uv := split(q)
+	if uv == zero {
+		return lift(math.Log(w))
+	}
+	v := Abs(uv)
+	return join(math.Log(Abs(q)), Scale(math.Atan2(v, w)/v, uv))
+}
+
+// Pow return q**r, the base-q exponential of r.
+// For generalized compatibility with math.Pow:
+//      Pow(0, ±0) returns 1+0i+0j+0k
+//      Pow(0, c) for real(c)<0 returns Inf+0i+0j+0k if imag(c), jmag(c), kmag(c) are zero,
+//          otherwise Inf+Inf i+Inf j+Inf k.
+func Pow(q, r Quat) Quat {
+	if q == zero {
+		w, uv := split(r)
+		switch {
+		case w == 0:
+			return Quat{Real: 1}
+		case w < 0:
+			if uv == zero {
+				return Quat{Real: math.Inf(1)}
+			}
+			return Inf()
+		case w > 0:
+			return zero
+		}
+	}
+	return Exp(Mul(Log(q), r))
+}
+
+// Sqrt returns the square root of q.
+func Sqrt(q Quat) Quat {
+	if q == zero {
+		return zero
+	}
+	return Pow(q, Quat{Real: 0.5})
+}