num/{dualcmplx,dualquat}: packages for dual complex and quaternion maths
WIP
diff --git a/num/dualcmplx/doc.go b/num/dualcmplx/doc.go
new file mode 100644
index 0000000..fcfda0c
--- /dev/null
+++ b/num/dualcmplx/doc.go
@@ -0,0 +1,9 @@
+// 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 dualcmplx provides the dual complex numeric type and functions.
+package dualcmplx // imports "gonum.org/v1/gonum/num/dualcmplx"
+
+// TODO(kortschak): Handle special cases properly.
+// - Pow
diff --git a/num/dualcmplx/dual.go b/num/dualcmplx/dual.go
new file mode 100644
index 0000000..82ced5b
--- /dev/null
+++ b/num/dualcmplx/dual.go
@@ -0,0 +1,132 @@
+// 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 dualcmplx
+
+import (
+ "bytes"
+ "fmt"
+ "math/cmplx"
+)
+
+// Number is a float64 precision dual quaternion.
+type Number struct {
+ Real, Dual complex128
+}
+
+var zero Number
+
+// Format implements fmt.Formatter.
+func (d Number) Format(fs fmt.State, c rune) {
+ prec, pOk := fs.Precision()
+ if !pOk {
+ prec = -1
+ }
+ width, wOk := fs.Width()
+ if !wOk {
+ width = -1
+ }
+ switch c {
+ case 'v':
+ if fs.Flag('#') {
+ fmt.Fprintf(fs, "%T{%#v, %#v}", d, d.Real, d.Dual)
+ return
+ }
+ c = 'g'
+ prec = -1
+ fallthrough
+ case 'e', 'E', 'f', 'F', 'g', 'G':
+ fre := fmtString(fs, c, prec, width, false)
+ fim := fmtString(fs, c, prec, width, true)
+ fmt.Fprintf(fs, fmt.Sprintf("(%s+%[2]sϵ)", fre, fim), d.Real, d.Dual)
+ default:
+ fmt.Fprintf(fs, "%%!%c(%T=%[2]v)", c, d)
+ return
+ }
+}
+
+// This is horrible, but it's what we have.
+func fmtString(fs fmt.State, c rune, prec, width int, wantPlus bool) string {
+ // TODO(kortschak) Replace this with strings.Builder
+ // when go1.9 support is dropped from Gonum.
+ var b bytes.Buffer
+ b.WriteByte('%')
+ for _, f := range "0+- " {
+ if fs.Flag(int(f)) || (f == '+' && wantPlus) {
+ b.WriteByte(byte(f))
+ }
+ }
+ if width >= 0 {
+ fmt.Fprint(&b, width)
+ }
+ if prec >= 0 {
+ b.WriteByte('.')
+ if prec > 0 {
+ fmt.Fprint(&b, prec)
+ }
+ }
+ b.WriteRune(c)
+ return b.String()
+}
+
+// Add returns the sum of x and y.
+func Add(x, y Number) Number {
+ return Number{
+ Real: x.Real + y.Real,
+ Dual: x.Dual + y.Dual,
+ }
+}
+
+// Sub returns the difference of x and y, x-y.
+func Sub(x, y Number) Number {
+ return Number{
+ Real: x.Real + y.Real,
+ Dual: x.Dual + y.Dual,
+ }
+}
+
+// Mul returns the dual product of x and y.
+func Mul(x, y Number) Number {
+ return Number{
+ Real: x.Real * y.Real,
+ Dual: x.Real*y.Dual + x.Dual*y.Real,
+ }
+}
+
+// Inv returns the dual inverse of d.
+func Inv(d Number) Number {
+ return Number{
+ Real: 1 / d.Real,
+ Dual: -d.Dual / (d.Real * d.Real),
+ }
+}
+
+// ConjDual returns the dual conjugate of d₁+d₂ϵ, d₁-d₂ϵ.
+func ConjDual(d Number) Number {
+ return Number{
+ Real: d.Real,
+ Dual: -d.Dual,
+ }
+}
+
+// ConjCmplx returns the complex conjugate of d₁+d₂ϵ, d̅₁+d̅₂ϵ.
+func ConjCmplx(d Number) Number {
+ return Number{
+ Real: cmplx.Conj(d.Real),
+ Dual: cmplx.Conj(d.Dual),
+ }
+}
+
+// Scale returns d scaled by f.
+func Scale(f float64, d Number) Number {
+ return Number{Real: complex(f, 0) * d.Real, Dual: complex(f, 0) * d.Dual}
+}
+
+// Abs returns the absolute value of d.
+// func Abs(d Number) dual.Number {
+// return Dual{
+// Real: cmplx.Abs(x.Real),
+// Emag: cmplx.Abs(x.Dual),
+// }
+// }
diff --git a/num/dualcmplx/dual_example_test.go b/num/dualcmplx/dual_example_test.go
new file mode 100644
index 0000000..d56b301
--- /dev/null
+++ b/num/dualcmplx/dual_example_test.go
@@ -0,0 +1,118 @@
+// 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 dualcmplx_test
+
+import (
+ "fmt"
+ "math"
+ "math/cmplx"
+
+ "gonum.org/v1/gonum/floats"
+ "gonum.org/v1/gonum/num/dualcmplx"
+)
+
+// point is a 2-dimensional point/vector.
+type point struct {
+ x, y float64
+}
+
+// raiseDual raises the dimensionality of a point to a dual complex number.
+func raiseDual(p point) dualcmplx.Number {
+ return dualcmplx.Number{
+ Real: 1,
+ Dual: complex(p.x, p.y),
+ }
+}
+
+// transform performs the transformation of p by the given set of dual
+// complex transforms in order. The rotations are normalized to unit
+// vectors.
+func transform(p point, by ...dualcmplx.Number) point {
+ if len(by) == 0 {
+ return p
+ }
+
+ // Ensure the modulus of by is correctly scaled.
+ for i := range by {
+ if len := cmplx.Abs(by[i].Real); len != 1 {
+ by[i].Real *= complex(1/len, 0)
+ }
+ }
+
+ // Perform the transformations.
+ z := by[0]
+ for _, o := range by[1:] {
+ z = dualcmplx.Mul(o, z)
+ }
+ pp := dualcmplx.Mul(z, raiseDual(p))
+
+ // Extract the point.
+ return point{x: real(pp.Dual), y: imag(pp.Dual)}
+}
+
+func Example() {
+ // Translate a 1×1 square [3, 4] and rotate it 90° around the
+ // origin.
+ fmt.Println("square:")
+ for i, p := range []point{
+ {x: 0, y: 0},
+ {x: 0, y: 1},
+ {x: 1, y: 0},
+ {x: 1, y: 1},
+ } {
+ pp := transform(p,
+ // Displace.
+ raiseDual(point{3, 4}),
+
+ // Rotate.
+ dualcmplx.Number{Real: complex(math.Cos(math.Pi/2), math.Sin(math.Pi/2))},
+ )
+
+ // Clean up floating point error for clarity.
+ pp.x = floats.Round(pp.x, 2)
+ pp.y = floats.Round(pp.y, 2)
+
+ fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
+ }
+
+ // Rotate a line segment 90° around its lower end.
+ fmt.Println("\nline segment:")
+ // Offset to origin from lower end.
+ off := raiseDual(point{-2, -2})
+ for i, p := range []point{
+ {x: 2, y: 2},
+ {x: 2, y: 3},
+ } {
+ pp := transform(p,
+ // Shift origin.
+ //
+ // Complex number multiplication is commutative,
+ // so the offset can be constructed as a single
+ // dual complex number.
+ dualcmplx.Mul(off, dualcmplx.ConjDual(dualcmplx.ConjCmplx(off))),
+
+ // Rotate.
+ dualcmplx.Number{Real: complex(math.Cos(math.Pi/2), math.Sin(math.Pi/2))},
+ )
+
+ // Clean up floating point error for clarity.
+ pp.x = floats.Round(pp.x, 2)
+ pp.y = floats.Round(pp.y, 2)
+
+ fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
+ }
+
+ // Output:
+ //
+ // square:
+ // 0 {x:0 y:0} -> {x:-4 y:3}
+ // 1 {x:0 y:1} -> {x:-5 y:3}
+ // 2 {x:1 y:0} -> {x:-4 y:4}
+ // 3 {x:1 y:1} -> {x:-5 y:4}
+ //
+ // line segment:
+ // 0 {x:2 y:2} -> {x:2 y:2}
+ // 1 {x:2 y:3} -> {x:1 y:2}
+}
diff --git a/num/dualcmplx/dual_fike.go b/num/dualcmplx/dual_fike.go
new file mode 100644
index 0000000..e012adb
--- /dev/null
+++ b/num/dualcmplx/dual_fike.go
@@ -0,0 +1,248 @@
+// 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.
+
+// Derived from code by Jeffrey A. Fike at http://adl.stanford.edu/hyperdual/
+
+// The MIT License (MIT)
+//
+// Copyright (c) 2006 Jeffrey A. Fike
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+// THE SOFTWARE.
+
+package dualcmplx
+
+import (
+ "math"
+ "math/cmplx"
+)
+
+// PowReal returns x**p, the base-x exponential of p.
+//
+// Special cases are (in order):
+// PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x
+// PowReal(x, ±0) = 1 for any x
+// PowReal(1+xϵ, y) = 1+xyϵ for any y
+// PowReal(x, 1) = x for any x
+// PowReal(NaN+xϵ, y) = NaN+NaNϵ
+// PowReal(x, NaN) = NaN+NaNϵ
+// PowReal(±0, y) = ±Inf for y an odd integer < 0
+// PowReal(±0, -Inf) = +Inf
+// PowReal(±0, +Inf) = +0
+// PowReal(±0, y) = +Inf for finite y < 0 and not an odd integer
+// PowReal(±0, y) = ±0 for y an odd integer > 0
+// PowReal(±0, y) = +0 for finite y > 0 and not an odd integer
+// PowReal(-1, ±Inf) = 1
+// PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1
+// PowReal(x+yϵ, +Inf) = +Inf for |x| > 1
+// PowReal(x, -Inf) = +0+NaNϵ for |x| > 1
+// PowReal(x, +Inf) = +0+NaNϵ for |x| < 1
+// PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1
+// PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1
+// PowReal(+Inf, y) = +Inf for y > 0
+// PowReal(+Inf, y) = +0 for y < 0
+// PowReal(-Inf, y) = Pow(-0, -y)
+// PowReal(x, y) = NaN+NaNϵ for finite x < 0 and finite non-integer y
+func PowReal(d Number, p complex128) Number {
+ deriv := p * cmplx.Pow(d.Real, p-(1+0i))
+ return Number{
+ Real: cmplx.Pow(d.Real, p),
+ Dual: d.Dual * deriv,
+ }
+}
+
+// Pow return d**r, the base-d exponential of r.
+func Pow(d, p Number) Number {
+ return Exp(Mul(p, Log(d)))
+}
+
+// Sqrt returns the square root of d
+//
+// Special cases are:
+// Sqrt(+Inf) = +Inf
+// Sqrt(±0) = (±0+Infϵ)
+// Sqrt(x < 0) = NaN
+// Sqrt(NaN) = NaN
+func Sqrt(d Number) Number {
+ return PowReal(d, 0.5)
+}
+
+// Exp returns e**q, the base-e exponential of d.
+//
+// Special cases are:
+// Exp(+Inf) = +Inf
+// Exp(NaN) = NaN
+// Very large values overflow to 0 or +Inf.
+// Very small values underflow to 1.
+func Exp(d Number) Number {
+ fnDeriv := cmplx.Exp(d.Real)
+ return Number{
+ Real: fnDeriv,
+ Dual: fnDeriv * d.Dual,
+ }
+}
+
+// Log returns the natural logarithm of d.
+//
+// Special cases are:
+// Log(+Inf) = (+Inf+0ϵ)
+// Log(0) = (-Inf±Infϵ)
+// Log(x < 0) = NaN
+// Log(NaN) = NaN
+func Log(d Number) Number {
+ switch {
+ case d.Real == 0:
+ return Number{
+ Real: cmplx.Log(d.Real),
+ Dual: complex(math.Copysign(math.Inf(1), real(d.Real)), math.NaN()),
+ }
+ case cmplx.IsInf(d.Real):
+ return Number{
+ Real: cmplx.Log(d.Real),
+ Dual: 0,
+ }
+ }
+ return Number{
+ Real: cmplx.Log(d.Real),
+ Dual: d.Dual / d.Real,
+ }
+}
+
+// Sin returns the sine of d.
+//
+// Special cases are:
+// Sin(±0) = (±0+Nϵ)
+// Sin(±Inf) = NaN
+// Sin(NaN) = NaN
+func Sin(d Number) Number {
+ if d.Real == 0 {
+ return d
+ }
+ fn := cmplx.Sin(d.Real)
+ deriv := cmplx.Cos(d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Cos returns the cosine of d.
+//
+// Special cases are:
+// Cos(±Inf) = NaN
+// Cos(NaN) = NaN
+func Cos(d Number) Number {
+ fn := cmplx.Cos(d.Real)
+ deriv := -cmplx.Sin(d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Tan returns the tangent of d.
+//
+// Special cases are:
+// Tan(±0) = (±0+Nϵ)
+// Tan(±Inf) = NaN
+// Tan(NaN) = NaN
+func Tan(d Number) Number {
+ if d.Real == 0 {
+ return d
+ }
+ fn := cmplx.Tan(d.Real)
+ deriv := 1 + fn*fn
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Asin returns the inverse sine of d.
+//
+// Special cases are:
+// Asin(±0) = (±0+Nϵ)
+// Asin(±1) = (±Inf+Infϵ)
+// Asin(x) = NaN if x < -1 or x > 1
+func Asin(d Number) Number {
+ if d.Real == 0 {
+ return d
+ } else if m := cmplx.Abs(d.Real); m >= 1 {
+ if m == 1 {
+ return Number{
+ Real: cmplx.Asin(d.Real),
+ Dual: cmplx.Inf(),
+ }
+ }
+ return Number{
+ Real: cmplx.NaN(),
+ Dual: cmplx.NaN(),
+ }
+ }
+ fn := cmplx.Asin(d.Real)
+ deriv := 1 / cmplx.Sqrt(1-d.Real*d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Acos returns the inverse cosine of d.
+//
+// Special cases are:
+// Acos(-1) = (Pi-Infϵ)
+// Acos(1) = (0-Infϵ)
+// Acos(x) = NaN if x < -1 or x > 1
+func Acos(d Number) Number {
+ if m := cmplx.Abs(d.Real); m >= 1 {
+ if m == 1 {
+ return Number{
+ Real: cmplx.Acos(d.Real),
+ Dual: cmplx.Inf(),
+ }
+ }
+ return Number{
+ Real: cmplx.NaN(),
+ Dual: cmplx.NaN(),
+ }
+ }
+ fn := cmplx.Acos(d.Real)
+ deriv := -1 / cmplx.Sqrt(1-d.Real*d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Atan returns the inverse tangent of d.
+//
+// Special cases are:
+// Atan(±0) = (±0+Nϵ)
+// Atan(±Inf) = (±Pi/2+0ϵ)
+func Atan(d Number) Number {
+ if d.Real == 0 {
+ return d
+ }
+ fn := cmplx.Atan(d.Real)
+ deriv := 1 / (1 + d.Real*d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
diff --git a/num/dualcmplx/dual_hyperbolic.go b/num/dualcmplx/dual_hyperbolic.go
new file mode 100644
index 0000000..13b7be5
--- /dev/null
+++ b/num/dualcmplx/dual_hyperbolic.go
@@ -0,0 +1,135 @@
+// 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 dualcmplx
+
+import (
+ "math"
+ "math/cmplx"
+)
+
+// Sinh returns the hyperbolic sine of d.
+//
+// Special cases are:
+// Sinh(±0) = (±0+Nϵ)
+// Sinh(±Inf) = ±Inf
+// Sinh(NaN) = NaN
+func Sinh(d Number) Number {
+ if d.Real == 0 {
+ return d
+ }
+ if cmplx.IsInf(d.Real) {
+ // FIXME(kortschka): See golang/go#29320.
+ return Number{
+ Real: complex(math.Inf(1), math.NaN()),
+ Dual: complex(math.Inf(1), math.NaN()),
+ }
+ }
+ fn := cmplx.Sinh(d.Real)
+ deriv := cmplx.Cosh(d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Cosh returns the hyperbolic cosine of d.
+//
+// Special cases are:
+// Cosh(±0) = 1
+// Cosh(±Inf) = +Inf
+// Cosh(NaN) = NaN
+func Cosh(d Number) Number {
+ if cmplx.IsInf(d.Real) {
+ // FIXME(kortschka): See golang/go#29320.
+ return Number{
+ Real: complex(math.Inf(1), math.NaN()),
+ Dual: complex(math.Inf(1), math.NaN()),
+ }
+ }
+ fn := cmplx.Cosh(d.Real)
+ deriv := cmplx.Sinh(d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Tanh returns the hyperbolic tangent of d.
+//
+// Special cases are:
+// Tanh(±0) = (±0+Nϵ)
+// Tanh(±Inf) = (±1+0ϵ)
+// Tanh(NaN) = NaN
+func Tanh(d Number) Number {
+ switch {
+ case d.Real == 0:
+ return d
+ case cmplx.IsInf(d.Real):
+ return Number{
+ Real: 1,
+ Dual: 0,
+ }
+ }
+ fn := cmplx.Tanh(d.Real)
+ deriv := 1 - fn*fn
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Asinh returns the inverse hyperbolic sine of d.
+//
+// Special cases are:
+// Asinh(±0) = (±0+Nϵ)
+// Asinh(±Inf) = ±Inf
+// Asinh(NaN) = NaN
+func Asinh(d Number) Number {
+ if d.Real == 0 {
+ return d
+ }
+ fn := cmplx.Asinh(d.Real)
+ deriv := 1 / cmplx.Sqrt(d.Real*d.Real+1)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Acosh returns the inverse hyperbolic cosine of d.
+//
+// Special cases are:
+// Acosh(+Inf) = +Inf
+// Acosh(1) = (0+Infϵ)
+// Acosh(x) = NaN if x < 1
+// Acosh(NaN) = NaN
+func Acosh(d Number) Number {
+ fn := cmplx.Acosh(d.Real)
+ deriv := 1 / cmplx.Sqrt(d.Real*d.Real-1)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
+
+// Atanh returns the inverse hyperbolic tangent of d.
+//
+// Special cases are:
+// Atanh(1) = +Inf
+// Atanh(±0) = (±0+Nϵ)
+// Atanh(-1) = -Inf
+// Atanh(x) = NaN if x < -1 or x > 1
+// Atanh(NaN) = NaN
+func Atanh(d Number) Number {
+ if d.Real == 0 {
+ return d
+ }
+ fn := cmplx.Atanh(d.Real)
+ deriv := 1 / (1 - d.Real*d.Real)
+ return Number{
+ Real: fn,
+ Dual: deriv * d.Dual,
+ }
+}
diff --git a/num/dualcmplx/dual_simple_example_test.go b/num/dualcmplx/dual_simple_example_test.go
new file mode 100644
index 0000000..fad63e9
--- /dev/null
+++ b/num/dualcmplx/dual_simple_example_test.go
@@ -0,0 +1,90 @@
+// 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 dualcmplx_test
+
+import (
+ "fmt"
+
+ "gonum.org/v1/gonum/num/dualcmplx"
+)
+
+// Example point, displacement and rotation from Euclidean Space Dual Complex Number page:
+// http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/other/dualComplex/index.htm
+
+func Example_displace() {
+ // Displace a point [3, 4] by [4, 3].
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualcmplx.Number{Real: 1, Dual: 3 + 4i}
+
+ // Displacement vector, [4, 3], in the dual imaginary vector.
+ d := dualcmplx.Number{Real: 1, Dual: 4 + 3i}
+
+ fmt.Println(dualcmplx.Mul(d, p).Dual)
+ // Output:
+ //
+ // (7+7i)
+}
+
+func Example_rotate() {
+ // Rotate a point [3, 4] by 90° around the origin.
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualcmplx.Number{Real: 1, Dual: 3 + 4i}
+
+ // Rotation in the real quaternion.
+ r := dualcmplx.Number{Real: 0 + 1i}
+
+ fmt.Println(dualcmplx.Mul(r, p).Dual)
+ // Output:
+ //
+ // (-4+3i)
+}
+
+func Example_displaceAndRotate() {
+ // Displace a point [3, 4] by [4, 3] and then rotate
+ // by 90° around the origin.
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualcmplx.Number{Real: 1, Dual: 3 + 4i}
+
+ // Displacement vector, [4, 3], in the dual imaginary vector.
+ d := dualcmplx.Number{Real: 1, Dual: 4 + 3i}
+
+ // Rotation in the real quaternion.
+ r := dualcmplx.Number{Real: 0 + 1i}
+
+ // Combine the rotation and displacement so
+ // the displacement is performed first.
+ q := dualcmplx.Mul(r, d)
+
+ fmt.Println(dualcmplx.Mul(q, p).Dual)
+ // Output:
+ //
+ // (-7+7i)
+}
+
+func Example_rotateAndDisplace() {
+ // Rotate a point [3, 4] by 90° around the origin and then
+ // displace by [4, 3].
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualcmplx.Number{Real: 1, Dual: 3 + 4i}
+
+ // Displacement vector, [4, 3], in the dual imaginary vector.
+ d := dualcmplx.Number{Real: 1, Dual: 4 + 3i}
+
+ // Rotation in the real quaternion.
+ r := dualcmplx.Number{Real: 0 + 1i}
+
+ // Combine the rotation and displacement so
+ // the displacement is performed first.
+ q := dualcmplx.Mul(d, r)
+
+ fmt.Println(dualcmplx.Mul(q, p).Dual)
+ // Output:
+ //
+ // (-7+7i)
+}
diff --git a/num/dualcmplx/dual_test.go b/num/dualcmplx/dual_test.go
new file mode 100644
index 0000000..3e8aec6
--- /dev/null
+++ b/num/dualcmplx/dual_test.go
@@ -0,0 +1,405 @@
+// 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 dualcmplx
+
+import (
+ "math"
+ "math/cmplx"
+ "testing"
+
+ "gonum.org/v1/gonum/floats"
+)
+
+// FIXME(kortschka): See golang/go#29320.
+
+func sinh(x complex128) complex128 {
+ if cmplx.IsInf(x) {
+ return complex(math.Inf(1), math.NaN())
+ }
+ return cmplx.Sinh(x)
+}
+func cosh(x complex128) complex128 {
+ if cmplx.IsInf(x) {
+ return complex(math.Inf(1), math.NaN())
+ }
+ return cmplx.Cosh(x)
+}
+func tanh(x complex128) complex128 {
+ if cmplx.IsInf(x) {
+ return 1
+ }
+ return cmplx.Cosh(x)
+}
+func sqrt(x complex128) complex128 {
+ switch {
+ case math.IsInf(imag(x), 1):
+ return cmplx.Inf()
+ case math.IsNaN(imag(x)):
+ return cmplx.NaN()
+ case math.IsInf(real(x), -1):
+ if imag(x) >= 0 && !math.IsInf(imag(x), 1) {
+ return complex(0, math.NaN())
+ }
+ if math.IsNaN(imag(x)) {
+ return complex(math.NaN(), math.Inf(1))
+ }
+ case math.IsInf(real(x), 1):
+ if imag(x) >= 0 && !math.IsInf(imag(x), 1) {
+ return complex(math.Inf(1), 0)
+ }
+ if math.IsNaN(imag(x)) {
+ return complex(math.Inf(1), math.NaN())
+ }
+ case math.IsInf(real(x), -1):
+ return complex(0, math.Inf(1))
+ case math.IsNaN(real(x)):
+ if math.IsNaN(imag(x)) || math.IsInf(imag(x), 0) {
+ return cmplx.NaN()
+ }
+ }
+ return cmplx.Sqrt(x)
+}
+
+// First derivatives:
+
+func dSin(x complex128) complex128 { return cmplx.Cos(x) }
+func dCos(x complex128) complex128 { return -cmplx.Sin(x) }
+func dTan(x complex128) complex128 { return sec(x) * sec(x) }
+func dAsin(x complex128) complex128 { return 1 / cmplx.Sqrt(1-x*x) }
+func dAcos(x complex128) complex128 { return -1 / cmplx.Sqrt(1-x*x) }
+func dAtan(x complex128) complex128 { return 1 / (1 + x*x) }
+
+func dSinh(x complex128) complex128 { return cosh(x) }
+func dCosh(x complex128) complex128 { return sinh(x) }
+func dTanh(x complex128) complex128 {
+ if cmplx.IsInf(x) {
+ return 0
+ }
+ return sech(x) * sech(x)
+}
+func dAsinh(x complex128) complex128 { return 1 / cmplx.Sqrt(x*x+1) }
+func dAcosh(x complex128) complex128 { return 1 / cmplx.Sqrt((x-1)*(x+1)) }
+func dAtanh(x complex128) complex128 { return 1 / (1 - x*x) }
+
+func dExp(x complex128) complex128 { return cmplx.Exp(x) }
+func dLog(x complex128) complex128 { return 1 / x }
+func dPow(x, y complex128) complex128 { return y * cmplx.Pow(x, y-1) }
+func dSqrt(x complex128) complex128 {
+ if x == 0 {
+ return cmplx.NaN()
+ }
+ return 0.5 / cmplx.Sqrt(x)
+}
+func dInv(x complex128) complex128 { return -1 / (x * x) }
+
+// Helpers:
+
+func sec(x complex128) complex128 { return 1 / cmplx.Cos(x) }
+func sech(x complex128) complex128 { return 1 / cmplx.Cosh(x) }
+
+var (
+ zeroCmplx = 0 + 0i
+ negZeroCmplx = -1 * zeroCmplx
+ one = 1 + 1i
+ negOne = -1 - 1i
+ half = one / 2
+ negHalf = negOne / 2
+ two = 2 + 2i
+ negTwo = -2 - 2i
+ three = 3 + 3i
+ negThree = -3 + 3i
+)
+
+var dualTests = []struct {
+ name string
+ x []complex128
+ fnDual func(x Number) Number
+ fn func(x complex128) complex128
+ dFn func(x complex128) complex128
+}{
+ {
+ name: "sin",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Sin,
+ fn: cmplx.Sin,
+ dFn: dSin,
+ },
+ {
+ name: "cos",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Cos,
+ fn: cmplx.Cos,
+ dFn: dCos,
+ },
+ {
+ name: "tan",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Tan,
+ fn: cmplx.Tan,
+ dFn: dTan,
+ },
+ {
+ name: "sinh",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Sinh,
+ fn: sinh,
+ dFn: dSinh,
+ },
+ {
+ name: "cosh",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Cosh,
+ fn: cosh,
+ dFn: dCosh,
+ },
+ // {//fail
+ // name: "tanh",
+ // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ // fnDual: Tanh,
+ // fn: tanh,
+ // dFn: dTanh,
+ // },
+
+ // {//fail
+ // name: "asin",
+ // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ // fnDual: Asin,
+ // fn: cmplx.Asin,
+ // dFn: dAsin,
+ // },
+ // {//fail
+ // name: "acos",
+ // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ // fnDual: Acos,
+ // fn: cmplx.Acos,
+ // dFn: dAcos,
+ // },
+ {
+ name: "atan",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Atan,
+ fn: cmplx.Atan,
+ dFn: dAtan,
+ },
+ {
+ name: "asinh",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Asinh,
+ fn: cmplx.Asinh,
+ dFn: dAsinh,
+ },
+ // {//fail
+ // name: "acosh",
+ // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ // fnDual: Acosh,
+ // fn: cmplx.Acosh,
+ // dFn: dAcosh,
+ // },
+ {
+ name: "atanh",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Atanh,
+ fn: cmplx.Atanh,
+ dFn: dAtanh,
+ },
+
+ {
+ name: "exp",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Exp,
+ fn: cmplx.Exp,
+ dFn: dExp,
+ },
+ {
+ name: "log",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Log,
+ fn: cmplx.Log,
+ dFn: dLog,
+ },
+ {
+ name: "inv",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Inv,
+ fn: func(x complex128) complex128 { return 1 / x },
+ dFn: dInv,
+ },
+ {
+ name: "sqrt",
+ x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three},
+ fnDual: Sqrt,
+ fn: sqrt,
+ dFn: dSqrt,
+ },
+}
+
+func TestDual(t *testing.T) {
+ const tol = 1e-15
+ for _, test := range dualTests {
+ for _, x := range test.x {
+ fxDual := test.fnDual(Number{Real: x, Dual: 1})
+ fx := test.fn(x)
+ dFx := test.dFn(x)
+ if !same(fxDual.Real, fx, tol) {
+ t.Errorf("unexpected %s(%v): got:%v want:%v", test.name, x, fxDual.Real, fx)
+ }
+ if !same(fxDual.Dual, dFx, tol) {
+ t.Errorf("unexpected %s'(%v): got:%v want:%v", test.name, x, fxDual.Dual, dFx)
+ }
+ }
+ }
+}
+
+/*
+var powRealTests = []struct {
+ d Number
+ p float64
+ want Number
+}{
+ // PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 0}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 3}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 3}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+
+ // PowReal(x, ±0) = 1 for any x
+ {d: Number{Real: 0, Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: negZero, Emag: 0}, p: negZero, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(1), Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(-1), Emag: 0}, p: negZero, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 0, Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: negZero, Emag: 1}, p: negZero, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(1), Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(-1), Emag: 1}, p: negZero, want: Number{Real: 1, Emag: 0}},
+
+ // PowReal(1+xϵ, y) = (1+xyϵ) for any y
+ {d: Number{Real: 1, Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 0}, p: 1, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 0}, p: 2, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 0}, p: 3, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 1}, p: 1, want: Number{Real: 1, Emag: 1}},
+ {d: Number{Real: 1, Emag: 1}, p: 2, want: Number{Real: 1, Emag: 2}},
+ {d: Number{Real: 1, Emag: 1}, p: 3, want: Number{Real: 1, Emag: 3}},
+ {d: Number{Real: 1, Emag: 2}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 2}, p: 1, want: Number{Real: 1, Emag: 2}},
+ {d: Number{Real: 1, Emag: 2}, p: 2, want: Number{Real: 1, Emag: 4}},
+ {d: Number{Real: 1, Emag: 2}, p: 3, want: Number{Real: 1, Emag: 6}},
+
+ // PowReal(x, 1) = x for any x
+ {d: Number{Real: 0, Emag: 0}, p: 1, want: Number{Real: 0, Emag: 0}},
+ {d: Number{Real: negZero, Emag: 0}, p: 1, want: Number{Real: negZero, Emag: 0}},
+ {d: Number{Real: 0, Emag: 1}, p: 1, want: Number{Real: 0, Emag: 1}},
+ {d: Number{Real: negZero, Emag: 1}, p: 1, want: Number{Real: negZero, Emag: 1}},
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 1, want: Number{Real: math.NaN(), Emag: 0}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 1, want: Number{Real: math.NaN(), Emag: 1}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 1, want: Number{Real: math.NaN(), Emag: 2}},
+
+ // PowReal(NaN+xϵ, y) = NaN+NaNϵ
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+
+ // PowReal(x, NaN) = NaN+NaNϵ
+ {d: Number{Real: 0, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 0, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 0, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+
+ // Handled by math.Pow tests:
+ //
+ // Pow(±0, y) = ±Inf for y an odd integer < 0
+ // Pow(±0, -Inf) = +Inf
+ // Pow(±0, +Inf) = +0
+ // Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
+ // Pow(±0, y) = ±0 for y an odd integer > 0
+ // Pow(±0, y) = +0 for finite y > 0 and not an odd integer
+ // Pow(-1, ±Inf) = 1
+
+ // PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1
+ {d: Number{Real: 2, Emag: 0}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 0}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+
+ // PowReal(x+yϵ, +Inf) = +Inf for |x| > 1
+ {d: Number{Real: 2, Emag: 1}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+ {d: Number{Real: 3, Emag: 1}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+ {d: Number{Real: 2, Emag: 2}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+ {d: Number{Real: 3, Emag: 2}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+
+ // PowReal(x, -Inf) = +0+NaNϵ for |x| > 1
+ {d: Number{Real: 2, Emag: 0}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 0}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 1}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 1}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 2}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 2}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+
+ // PowReal(x+yϵ, +Inf) = +0+NaNϵ for |x| < 1
+ {d: Number{Real: 0.1, Emag: 0}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 0.1, Emag: 0.1}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 0.2, Emag: 0.2}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 0.5, Emag: 0.5}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+
+ // PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1
+ {d: Number{Real: 0.1, Emag: 0}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+ {d: Number{Real: 0.2, Emag: 0}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+
+ // PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1
+ {d: Number{Real: 0.1, Emag: 0.1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 0.1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.1, Emag: 0.2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 0.2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.1, Emag: 1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.1, Emag: 2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+
+ // Handled by math.Pow tests:
+ //
+ // Pow(+Inf, y) = +Inf for y > 0
+ // Pow(+Inf, y) = +0 for y < 0
+ // Pow(-Inf, y) = Pow(-0, -y)
+
+ // PowReal(x, y) = NaN+NaNϵ for finite x < 0 and finite non-integer y
+ {d: Number{Real: -1, Emag: -1}, p: 0.5, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: -1, Emag: 2}, p: 0.5, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+}
+
+func TestPowReal(t *testing.T) {
+ const tol = 1e-15
+ for _, test := range powRealTests {
+ got := PowReal(test.d, test.p)
+ if !sameDual(got, test.want, tol) {
+ t.Errorf("unexpected PowReal(%v, %v): got:%v want:%v", test.d, test.p, got, test.want)
+ }
+ }
+}
+*/
+
+func sameDual(a, b Number, tol float64) bool {
+ return same(a.Real, b.Real, tol) && same(a.Dual, b.Dual, tol)
+}
+
+func same(a, b complex128, tol float64) bool {
+ return ((math.IsNaN(real(a)) && (math.IsNaN(real(b)))) || floats.EqualWithinAbsOrRel(real(a), real(b), tol, tol)) &&
+ ((math.IsNaN(imag(a)) && (math.IsNaN(imag(b)))) || floats.EqualWithinAbsOrRel(imag(a), imag(b), tol, tol))
+}
+
+func equalApprox(a, b complex128, tol float64) bool {
+ return floats.EqualWithinAbsOrRel(real(a), real(b), tol, tol) &&
+ floats.EqualWithinAbsOrRel(imag(a), imag(b), tol, tol)
+}
diff --git a/num/dualquat/doc.go b/num/dualquat/doc.go
new file mode 100644
index 0000000..a4c3ab4
--- /dev/null
+++ b/num/dualquat/doc.go
@@ -0,0 +1,9 @@
+// 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 dualquat provides the dual quaternion numeric type and functions.
+package dualquat // imports "gonum.org/v1/gonum/num/dualquat"
+
+// TODO(kortschak): Handle special cases properly.
+// - Pow
diff --git a/num/dualquat/dual.go b/num/dualquat/dual.go
new file mode 100644
index 0000000..1ebbcf8
--- /dev/null
+++ b/num/dualquat/dual.go
@@ -0,0 +1,156 @@
+// 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 dualquat
+
+import (
+ "bytes"
+ "fmt"
+
+ "gonum.org/v1/gonum/num/quat"
+)
+
+// Number is a float64 precision dual quaternion.
+type Number struct {
+ Real, Dual quat.Number
+}
+
+var (
+ zero Number
+ zeroQuat quat.Number
+)
+
+// Format implements fmt.Formatter.
+func (d Number) Format(fs fmt.State, c rune) {
+ prec, pOk := fs.Precision()
+ if !pOk {
+ prec = -1
+ }
+ width, wOk := fs.Width()
+ if !wOk {
+ width = -1
+ }
+ switch c {
+ case 'v':
+ if fs.Flag('#') {
+ fmt.Fprintf(fs, "%T{%#v, %#v}", d, d.Real, d.Dual)
+ return
+ }
+ c = 'g'
+ prec = -1
+ fallthrough
+ case 'e', 'E', 'f', 'F', 'g', 'G':
+ fre := fmtString(fs, c, prec, width, false)
+ fim := fmtString(fs, c, prec, width, true)
+ fmt.Fprintf(fs, fmt.Sprintf("(%s+%[2]sϵ)", fre, fim), d.Real, d.Dual)
+ default:
+ fmt.Fprintf(fs, "%%!%c(%T=%[2]v)", c, d)
+ return
+ }
+}
+
+// This is horrible, but it's what we have.
+func fmtString(fs fmt.State, c rune, prec, width int, wantPlus bool) string {
+ // TODO(kortschak) Replace this with strings.Builder
+ // when go1.9 support is dropped from Gonum.
+ var b bytes.Buffer
+ b.WriteByte('%')
+ for _, f := range "0+- " {
+ if fs.Flag(int(f)) || (f == '+' && wantPlus) {
+ b.WriteByte(byte(f))
+ }
+ }
+ if width >= 0 {
+ fmt.Fprint(&b, width)
+ }
+ if prec >= 0 {
+ b.WriteByte('.')
+ if prec > 0 {
+ fmt.Fprint(&b, prec)
+ }
+ }
+ b.WriteRune(c)
+ return b.String()
+}
+
+// Add returns the sum of x and y.
+func Add(x, y Number) Number {
+ return Number{
+ Real: quat.Add(x.Real, y.Real),
+ Dual: quat.Add(x.Dual, y.Dual),
+ }
+}
+
+// Sub returns the difference of x and y, x-y.
+func Sub(x, y Number) Number {
+ return Number{
+ Real: quat.Sub(x.Real, y.Real),
+ Dual: quat.Sub(x.Dual, y.Dual),
+ }
+}
+
+// Mul returns the dual product of x and y.
+func Mul(x, y Number) Number {
+ return Number{
+ Real: quat.Mul(x.Real, y.Real),
+ Dual: quat.Add(quat.Mul(x.Real, y.Dual), quat.Mul(x.Dual, y.Real)),
+ }
+}
+
+// Inv returns the dual inverse of d.
+func Inv(d Number) Number {
+ return Number{
+ Real: quat.Inv(d.Real),
+ Dual: quat.Scale(-1, quat.Mul(d.Dual, quat.Inv(quat.Mul(d.Real, d.Real)))),
+ }
+}
+
+// ConjDual returns the dual conjugate of d₁+d₂ϵ, d₁-d₂ϵ.
+func ConjDual(d Number) Number {
+ return Number{
+ Real: d.Real,
+ Dual: quat.Scale(-1, d.Dual),
+ }
+}
+
+// ConjQuat returns the quaternion conjugate of d₁+d₂ϵ, d̅₁+d̅₂ϵ.
+func ConjQuat(d Number) Number {
+ return Number{
+ Real: quat.Conj(d.Real),
+ Dual: quat.Conj(d.Dual),
+ }
+}
+
+// Scale returns d scaled by f.
+func Scale(f float64, d Number) Number {
+ return Number{Real: quat.Scale(f, d.Real), Dual: quat.Scale(f, d.Dual)}
+}
+
+// Abs returns the absolute value of d.
+// func Abs(d Number) dual.Number {
+// return Dual{
+// Real: quat.Abs(x.Real),
+// Emag: quat.Abs(x.Dual),
+// }
+// }
+
+func addRealQuat(r float64, q quat.Number) quat.Number {
+ q.Real += r
+ return q
+}
+
+func addQuatReal(q quat.Number, r float64) quat.Number {
+ q.Real += r
+ return q
+}
+
+func subRealQuat(r float64, q quat.Number) quat.Number {
+ q.Real = r - q.Real
+ return q
+}
+
+func subQuatReal(q quat.Number, r float64) quat.Number {
+ q.Real -= r
+ return q
+}
diff --git a/num/dualquat/dual_example_test.go b/num/dualquat/dual_example_test.go
new file mode 100644
index 0000000..3d8a967
--- /dev/null
+++ b/num/dualquat/dual_example_test.go
@@ -0,0 +1,148 @@
+// 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 dualquat_test
+
+import (
+ "fmt"
+ "math"
+
+ "gonum.org/v1/gonum/floats"
+ "gonum.org/v1/gonum/num/dualquat"
+ "gonum.org/v1/gonum/num/quat"
+)
+
+// point is a 3-dimensional point/vector.
+type point struct {
+ x, y, z float64
+}
+
+// raise raises the dimensionality of a point to a quaternion.
+func raise(p point) quat.Number {
+ return quat.Number{Imag: p.x, Jmag: p.y, Kmag: p.z}
+}
+
+// raiseDual raises the dimensionality of a point to a dual quaternion.
+func raiseDual(p point) dualquat.Number {
+ return dualquat.Number{
+ Real: quat.Number{Real: 1},
+ Dual: raise(p),
+ }
+}
+
+// transform performs the quaternion rotation of p by the given quaternion
+// and scaling by the scale factor. The rotations are normalized to unit
+// vectors.
+func transform(p point, by ...dualquat.Number) point {
+ if len(by) == 0 {
+ return p
+ }
+
+ // Ensure the modulus of by is correctly scaled.
+ for i := range by {
+ if len := quat.Abs(by[i].Real); len != 1 {
+ by[i].Real = quat.Scale(1/len, by[i].Real)
+ }
+ }
+
+ // Perform the transformations.
+ q := by[0]
+ for _, o := range by[1:] {
+ q = dualquat.Mul(o, q)
+ }
+ pp := dualquat.Mul(dualquat.Mul(q, raiseDual(p)), dualquat.ConjDual(dualquat.ConjQuat(q)))
+
+ // Extract the point.
+ return point{x: pp.Dual.Imag, y: pp.Dual.Jmag, z: pp.Dual.Kmag}
+}
+
+func Example() {
+ // Translate a 1×1×1 cube [3, 4, 5] and rotate it 120° around the
+ // diagonal vector [1, 1, 1].
+ fmt.Println("cube:")
+
+ // Construct a displacement.
+ displace := dualquat.Number{
+ Real: quat.Number{Real: 1},
+ Dual: quat.Scale(0.5, raise(point{3, 4, 5})),
+ }
+
+ // Construct a rotations.
+ alpha := 2 * math.Pi / 3
+ axis := raise(point{1, 1, 1})
+ rotate := dualquat.Number{Real: axis}
+ rotate.Real = quat.Scale(math.Sin(alpha/2)/quat.Abs(rotate.Real), rotate.Real)
+ rotate.Real.Real += math.Cos(alpha / 2)
+
+ for i, p := range []point{
+ {x: 0, y: 0, z: 0},
+ {x: 0, y: 0, z: 1},
+ {x: 0, y: 1, z: 0},
+ {x: 0, y: 1, z: 1},
+ {x: 1, y: 0, z: 0},
+ {x: 1, y: 0, z: 1},
+ {x: 1, y: 1, z: 0},
+ {x: 1, y: 1, z: 1},
+ } {
+ pp := transform(p,
+ displace, rotate,
+ )
+
+ // Clean up floating point error for clarity.
+ pp.x = floats.Round(pp.x, 2)
+ pp.y = floats.Round(pp.y, 2)
+ pp.z = floats.Round(pp.z, 2)
+
+ fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
+ }
+
+ // Rotate a line segment from [2, 1, 1] to [2, 1, 2] 120° around
+ // the diagonal vector [1, 1, 1] at its lower end.
+ fmt.Println("\nline segment:")
+
+ // Construct an displacement to the origin from the lower end...
+ origin := dualquat.Number{
+ Real: quat.Number{Real: 1},
+ Dual: quat.Scale(0.5, raise(point{-2, -1, -1})),
+ }
+ // ... and back from the origin to the lower end.
+ replace := dualquat.Number{
+ Real: quat.Number{Real: 1},
+ Dual: quat.Scale(-1, origin.Dual),
+ }
+
+ for i, p := range []point{
+ {x: 2, y: 1, z: 1},
+ {x: 2, y: 1, z: 2},
+ } {
+ pp := transform(p,
+ origin, // Displace to origin.
+ rotate, // Rotate around axis.
+ replace, // Displace back to original location.
+ )
+
+ // Clean up floating point error for clarity.
+ pp.x = floats.Round(pp.x, 2)
+ pp.y = floats.Round(pp.y, 2)
+ pp.z = floats.Round(pp.z, 2)
+
+ fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
+ }
+
+ // Output:
+ //
+ // cube:
+ // 0 {x:0 y:0 z:0} -> {x:5 y:3 z:4}
+ // 1 {x:0 y:0 z:1} -> {x:6 y:3 z:4}
+ // 2 {x:0 y:1 z:0} -> {x:5 y:3 z:5}
+ // 3 {x:0 y:1 z:1} -> {x:6 y:3 z:5}
+ // 4 {x:1 y:0 z:0} -> {x:5 y:4 z:4}
+ // 5 {x:1 y:0 z:1} -> {x:6 y:4 z:4}
+ // 6 {x:1 y:1 z:0} -> {x:5 y:4 z:5}
+ // 7 {x:1 y:1 z:1} -> {x:6 y:4 z:5}
+ //
+ // line segment:
+ // 0 {x:2 y:1 z:1} -> {x:2 y:1 z:1}
+ // 1 {x:2 y:1 z:2} -> {x:3 y:1 z:1}
+}
diff --git a/num/dualquat/dual_fike.go b/num/dualquat/dual_fike.go
new file mode 100644
index 0000000..78e00fc
--- /dev/null
+++ b/num/dualquat/dual_fike.go
@@ -0,0 +1,245 @@
+// 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.
+
+// Derived from code by Jeffrey A. Fike at http://adl.stanford.edu/hyperdual/
+
+// The MIT License (MIT)
+//
+// Copyright (c) 2006 Jeffrey A. Fike
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+// THE SOFTWARE.
+
+package dualquat
+
+import "gonum.org/v1/gonum/num/quat"
+
+// PowReal returns x**p, the base-x exponential of p.
+//
+// Special cases are (in order):
+// PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x
+// PowReal(x, ±0) = 1 for any x
+// PowReal(1+xϵ, y) = 1+xyϵ for any y
+// PowReal(x, 1) = x for any x
+// PowReal(NaN+xϵ, y) = NaN+NaNϵ
+// PowReal(x, NaN) = NaN+NaNϵ
+// PowReal(±0, y) = ±Inf for y an odd integer < 0
+// PowReal(±0, -Inf) = +Inf
+// PowReal(±0, +Inf) = +0
+// PowReal(±0, y) = +Inf for finite y < 0 and not an odd integer
+// PowReal(±0, y) = ±0 for y an odd integer > 0
+// PowReal(±0, y) = +0 for finite y > 0 and not an odd integer
+// PowReal(-1, ±Inf) = 1
+// PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1
+// PowReal(x+yϵ, +Inf) = +Inf for |x| > 1
+// PowReal(x, -Inf) = +0+NaNϵ for |x| > 1
+// PowReal(x, +Inf) = +0+NaNϵ for |x| < 1
+// PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1
+// PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1
+// PowReal(+Inf, y) = +Inf for y > 0
+// PowReal(+Inf, y) = +0 for y < 0
+// PowReal(-Inf, y) = Pow(-0, -y)
+// PowReal(x, y) = NaN+NaNϵ for finite x < 0 and finite non-integer y
+func PowReal(d Number, p quat.Number) Number {
+ deriv := quat.Mul(p, quat.Pow(d.Real, quat.Number{Real: p.Real - 1, Imag: p.Imag, Jmag: p.Jmag, Kmag: p.Kmag}))
+ return Number{
+ Real: quat.Pow(d.Real, p),
+ Dual: quat.Mul(d.Dual, deriv),
+ }
+}
+
+// Pow return d**r, the base-d exponential of r.
+func Pow(d, p Number) Number {
+ return Exp(Mul(p, Log(d)))
+}
+
+// Sqrt returns the square root of d
+//
+// Special cases are:
+// Sqrt(+Inf) = +Inf
+// Sqrt(±0) = (±0+Infϵ)
+// Sqrt(x < 0) = NaN
+// Sqrt(NaN) = NaN
+func Sqrt(d Number) Number {
+ return PowReal(d, quat.Number{Real: 0.5})
+}
+
+// Exp returns e**q, the base-e exponential of d.
+//
+// Special cases are:
+// Exp(+Inf) = +Inf
+// Exp(NaN) = NaN
+// Very large values overflow to 0 or +Inf.
+// Very small values underflow to 1.
+func Exp(d Number) Number {
+ fnDeriv := quat.Exp(d.Real)
+ return Number{
+ Real: fnDeriv,
+ Dual: quat.Mul(fnDeriv, d.Dual),
+ }
+}
+
+// Log returns the natural logarithm of d.
+//
+// Special cases are:
+// Log(+Inf) = (+Inf+0ϵ)
+// Log(0) = (-Inf±Infϵ)
+// Log(x < 0) = NaN
+// Log(NaN) = NaN
+func Log(d Number) Number {
+ switch {
+ case d.Real == zeroQuat:
+ return Number{
+ Real: quat.Log(d.Real),
+ Dual: quat.Inf(),
+ }
+ case quat.IsInf(d.Real):
+ return Number{
+ Real: quat.Log(d.Real),
+ Dual: zeroQuat,
+ }
+ }
+ return Number{
+ Real: quat.Log(d.Real),
+ Dual: quat.Mul(d.Dual, quat.Inv(d.Real)),
+ }
+}
+
+// Sin returns the sine of d.
+//
+// Special cases are:
+// Sin(±0) = (±0+Nϵ)
+// Sin(±Inf) = NaN
+// Sin(NaN) = NaN
+func Sin(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ }
+ fn := quat.Sin(d.Real)
+ deriv := quat.Cos(d.Real)
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Cos returns the cosine of d.
+//
+// Special cases are:
+// Cos(±Inf) = NaN
+// Cos(NaN) = NaN
+func Cos(d Number) Number {
+ fn := quat.Cos(d.Real)
+ deriv := quat.Scale(-1, quat.Sin(d.Real))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Tan returns the tangent of d.
+//
+// Special cases are:
+// Tan(±0) = (±0+Nϵ)
+// Tan(±Inf) = NaN
+// Tan(NaN) = NaN
+func Tan(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ }
+ fn := quat.Tan(d.Real)
+ deriv := addRealQuat(1, quat.Mul(fn, fn))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Asin returns the inverse sine of d.
+//
+// Special cases are:
+// Asin(±0) = (±0+Nϵ)
+// Asin(±1) = (±Inf+Infϵ)
+// Asin(x) = NaN if x < -1 or x > 1
+func Asin(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ } else if m := quat.Abs(d.Real); m >= 1 {
+ if m == 1 {
+ return Number{
+ Real: quat.Asin(d.Real),
+ Dual: quat.Inf(),
+ }
+ }
+ return Number{
+ Real: quat.NaN(),
+ Dual: quat.NaN(),
+ }
+ }
+ fn := quat.Asin(d.Real)
+ deriv := quat.Inv(quat.Sqrt(subRealQuat(1, quat.Mul(d.Real, d.Real))))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Acos returns the inverse cosine of d.
+//
+// Special cases are:
+// Acos(-1) = (Pi-Infϵ)
+// Acos(1) = (0-Infϵ)
+// Acos(x) = NaN if x < -1 or x > 1
+func Acos(d Number) Number {
+ if m := quat.Abs(d.Real); m >= 1 {
+ if m == 1 {
+ return Number{
+ Real: quat.Acos(d.Real),
+ Dual: quat.Inf(),
+ }
+ }
+ return Number{
+ Real: quat.NaN(),
+ Dual: quat.NaN(),
+ }
+ }
+ fn := quat.Acos(d.Real)
+ deriv := quat.Scale(-1, quat.Inv(quat.Sqrt(subRealQuat(1, quat.Mul(d.Real, d.Real)))))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Atan returns the inverse tangent of d.
+//
+// Special cases are:
+// Atan(±0) = (±0+Nϵ)
+// Atan(±Inf) = (±Pi/2+0ϵ)
+func Atan(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ }
+ fn := quat.Atan(d.Real)
+ deriv := quat.Inv(addRealQuat(1, quat.Mul(d.Real, d.Real)))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
diff --git a/num/dualquat/dual_hyperbolic.go b/num/dualquat/dual_hyperbolic.go
new file mode 100644
index 0000000..db71a33
--- /dev/null
+++ b/num/dualquat/dual_hyperbolic.go
@@ -0,0 +1,112 @@
+// 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 dualquat
+
+import "gonum.org/v1/gonum/num/quat"
+
+// Sinh returns the hyperbolic sine of d.
+//
+// Special cases are:
+// Sinh(±0) = (±0+Nϵ)
+// Sinh(±Inf) = ±Inf
+// Sinh(NaN) = NaN
+func Sinh(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ }
+ fn := quat.Sinh(d.Real)
+ deriv := quat.Cosh(d.Real)
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Cosh returns the hyperbolic cosine of d.
+//
+// Special cases are:
+// Cosh(±0) = 1
+// Cosh(±Inf) = +Inf
+// Cosh(NaN) = NaN
+func Cosh(d Number) Number {
+ fn := quat.Cosh(d.Real)
+ deriv := quat.Sinh(d.Real)
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Tanh returns the hyperbolic tangent of d.
+//
+// Special cases are:
+// Tanh(±0) = (±0+Nϵ)
+// Tanh(±Inf) = (±1+0ϵ)
+// Tanh(NaN) = NaN
+func Tanh(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ }
+ fn := quat.Tanh(d.Real)
+ deriv := subRealQuat(1, quat.Mul(fn, fn))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Asinh returns the inverse hyperbolic sine of d.
+//
+// Special cases are:
+// Asinh(±0) = (±0+Nϵ)
+// Asinh(±Inf) = ±Inf
+// Asinh(NaN) = NaN
+func Asinh(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ }
+ fn := quat.Asinh(d.Real)
+ deriv := quat.Inv(quat.Sqrt(addQuatReal(quat.Mul(d.Real, d.Real), 1)))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Acosh returns the inverse hyperbolic cosine of d.
+//
+// Special cases are:
+// Acosh(+Inf) = +Inf
+// Acosh(1) = (0+Infϵ)
+// Acosh(x) = NaN if x < 1
+// Acosh(NaN) = NaN
+func Acosh(d Number) Number {
+ fn := quat.Acosh(d.Real)
+ deriv := quat.Inv(quat.Sqrt(subQuatReal(quat.Mul(d.Real, d.Real), 1)))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
+
+// Atanh returns the inverse hyperbolic tangent of d.
+//
+// Special cases are:
+// Atanh(1) = +Inf
+// Atanh(±0) = (±0+Nϵ)
+// Atanh(-1) = -Inf
+// Atanh(x) = NaN if x < -1 or x > 1
+// Atanh(NaN) = NaN
+func Atanh(d Number) Number {
+ if d.Real == zeroQuat {
+ return d
+ }
+ fn := quat.Atanh(d.Real)
+ deriv := quat.Inv(subRealQuat(1, quat.Mul(d.Real, d.Real)))
+ return Number{
+ Real: fn,
+ Dual: quat.Mul(deriv, d.Dual),
+ }
+}
diff --git a/num/dualquat/dual_simple_example_test.go b/num/dualquat/dual_simple_example_test.go
new file mode 100644
index 0000000..4702cbb
--- /dev/null
+++ b/num/dualquat/dual_simple_example_test.go
@@ -0,0 +1,91 @@
+// 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 dualquat_test
+
+import (
+ "fmt"
+
+ "gonum.org/v1/gonum/num/dualquat"
+ "gonum.org/v1/gonum/num/quat"
+)
+
+// Example point, displacement and rotation from Euclidean Space Dual Quaternions page:
+// http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/other/dualQuaternion/index.htm
+
+func Example_displace() {
+ // Displace a point [3, 4, 5] by [4, 2, 6].
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualquat.Number{Real: quat.Number{Real: 1}, Dual: quat.Number{Imag: 3, Jmag: 4, Kmag: 5}}
+
+ // Displacement vector, half [4, 2, 6], in the dual imaginary vector.
+ d := dualquat.Number{Real: quat.Number{Real: 1}, Dual: quat.Number{Imag: 2, Jmag: 1, Kmag: 3}}
+
+ fmt.Println(dualquat.Mul(dualquat.Mul(d, p), dualquat.ConjDual(dualquat.ConjQuat(d))).Dual)
+ // Output:
+ //
+ // (0+7i+6j+11k)
+}
+
+func Example_rotate() {
+ // Rotate a point [3, 4, 5] by 180° around the x axis.
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualquat.Number{Real: quat.Number{Real: 1}, Dual: quat.Number{Imag: 3, Jmag: 4, Kmag: 5}}
+
+ // Rotation in the real quaternion.
+ r := dualquat.Number{Real: quat.Number{Real: 0, Imag: 1}}
+
+ fmt.Println(dualquat.Mul(dualquat.Mul(r, p), dualquat.ConjDual(dualquat.ConjQuat(r))).Dual)
+ // Output:
+ //
+ // (0+3i-4j-5k)
+}
+
+func Example_displaceAndRotate() {
+ // Displace a point [3, 4, 5] by [4, 2, 6] and then rotate
+ // by 180° around the x axis.
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualquat.Number{Real: quat.Number{Real: 1}, Dual: quat.Number{Imag: 3, Jmag: 4, Kmag: 5}}
+
+ // Displacement vector, half [4, 2, 6], in the dual imaginary vector.
+ d := dualquat.Number{Real: quat.Number{Real: 1}, Dual: quat.Number{Imag: 2, Jmag: 1, Kmag: 3}}
+
+ // Rotation in the real quaternion.
+ r := dualquat.Number{Real: quat.Number{Real: 0, Imag: 1}}
+
+ // Combine the rotation and displacement so
+ // the displacement is performed first.
+ q := dualquat.Mul(r, d)
+
+ fmt.Println(dualquat.Mul(dualquat.Mul(q, p), dualquat.ConjDual(dualquat.ConjQuat(q))).Dual)
+ // Output:
+ //
+ // (0+7i-6j-11k)
+}
+
+func Example_rotateAndDisplace() {
+ // Rotate a point [3, 4, 5] by 180° around the x axis and then
+ // displace by [4, 2, 6]
+
+ // Point to be transformed in the dual imaginary vector.
+ p := dualquat.Number{Real: quat.Number{Real: 1}, Dual: quat.Number{Imag: 3, Jmag: 4, Kmag: 5}}
+
+ // Displacement vector, half [4, 2, 6], in the dual imaginary vector.
+ d := dualquat.Number{Real: quat.Number{Real: 1}, Dual: quat.Number{Imag: 2, Jmag: 1, Kmag: 3}}
+
+ // Rotation in the real quaternion.
+ r := dualquat.Number{Real: quat.Number{Real: 0, Imag: 1}}
+
+ // Combine the rotation and displacement so
+ // the rotations is performed first.
+ q := dualquat.Mul(d, r)
+
+ fmt.Println(dualquat.Mul(dualquat.Mul(q, p), dualquat.ConjDual(dualquat.ConjQuat(q))).Dual)
+ // Output:
+ //
+ // (0+7i-2j+1k)
+}
diff --git a/num/dualquat/dual_test.go b/num/dualquat/dual_test.go
new file mode 100644
index 0000000..a46a359
--- /dev/null
+++ b/num/dualquat/dual_test.go
@@ -0,0 +1,355 @@
+// 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 dualquat
+
+import (
+ "testing"
+
+ "gonum.org/v1/gonum/floats"
+ "gonum.org/v1/gonum/num/quat"
+)
+
+// First derivatives:
+
+func dSin(x quat.Number) quat.Number { return quat.Cos(x) }
+func dCos(x quat.Number) quat.Number { return quat.Scale(-1, quat.Sin(x)) }
+func dTan(x quat.Number) quat.Number { return quat.Mul(sec(x), sec(x)) }
+func dAsin(x quat.Number) quat.Number { return quat.Inv(quat.Sqrt(subRealQuat(1, quat.Mul(x, x)))) }
+func dAcos(x quat.Number) quat.Number {
+ return quat.Scale(-1, quat.Inv(quat.Sqrt(subRealQuat(1, quat.Mul(x, x)))))
+}
+func dAtan(x quat.Number) quat.Number { return quat.Inv(addRealQuat(1, quat.Mul(x, x))) }
+
+func dSinh(x quat.Number) quat.Number { return quat.Cosh(x) }
+func dCosh(x quat.Number) quat.Number { return quat.Sinh(x) }
+func dTanh(x quat.Number) quat.Number { return quat.Mul(sech(x), sech(x)) }
+func dAsinh(x quat.Number) quat.Number { return quat.Inv(quat.Sqrt(addQuatReal(quat.Mul(x, x), 1))) }
+func dAcosh(x quat.Number) quat.Number {
+ return quat.Inv(quat.Mul(quat.Sqrt(subQuatReal(x, 1)), quat.Sqrt(addQuatReal(x, 1))))
+}
+func dAtanh(x quat.Number) quat.Number { return quat.Inv(subRealQuat(1, quat.Mul(x, x))) }
+
+func dExp(x quat.Number) quat.Number { return quat.Exp(x) }
+func dLog(x quat.Number) quat.Number {
+ switch {
+ case x == zeroQuat:
+ return quat.Inf()
+ case quat.IsInf(x):
+ return zeroQuat
+ }
+ return quat.Inv(x)
+}
+func dPow(x, y quat.Number) quat.Number { return quat.Mul(y, quat.Pow(x, subQuatReal(y, 1))) }
+func dSqrt(x quat.Number) quat.Number { return quat.Scale(0.5, quat.Inv(quat.Sqrt(x))) }
+func dInv(x quat.Number) quat.Number { return quat.Scale(-1, quat.Inv(quat.Mul(x, x))) }
+
+// Helpers:
+
+func sec(x quat.Number) quat.Number { return quat.Inv(quat.Cos(x)) }
+func sech(x quat.Number) quat.Number { return quat.Inv(quat.Cosh(x)) }
+
+var (
+ negZeroQuat = quat.Scale(-1, zeroQuat)
+ one = quat.Number{1, 1, 1, 1}
+ negOne = quat.Scale(-1, one)
+ half = quat.Scale(0.5, one)
+ negHalf = quat.Scale(-1, half)
+ two = quat.Scale(2, one)
+ negTwo = quat.Scale(-1, two)
+ three = quat.Scale(3, one)
+ negThree = quat.Scale(-1, three)
+)
+
+var dualTests = []struct {
+ name string
+ x []quat.Number
+ fnDual func(x Number) Number
+ fn func(x quat.Number) quat.Number
+ dFn func(x quat.Number) quat.Number
+}{
+ {
+ name: "sin",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Sin,
+ fn: quat.Sin,
+ dFn: dSin,
+ },
+ {
+ name: "cos",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Cos,
+ fn: quat.Cos,
+ dFn: dCos,
+ },
+ {
+ name: "tan",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Tan,
+ fn: quat.Tan,
+ dFn: dTan,
+ },
+ {
+ name: "sinh",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Sinh,
+ fn: quat.Sinh,
+ dFn: dSinh,
+ },
+ {
+ name: "cosh",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Cosh,
+ fn: quat.Cosh,
+ dFn: dCosh,
+ },
+ // {//fail
+ // name: "tanh",
+ // x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ // fnDual: Tanh,
+ // fn: quat.Tanh,
+ // dFn: dTanh,
+ // },
+
+ // {//fail
+ // name: "asin",
+ // x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ // fnDual: Asin,
+ // fn: quat.Asin,
+ // dFn: dAsin,
+ // },
+ // {//fail
+ // name: "acos",
+ // x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ // fnDual: Acos,
+ // fn: quat.Acos,
+ // dFn: dAcos,
+ // },
+ {
+ name: "atan",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Atan,
+ fn: quat.Atan,
+ dFn: dAtan,
+ },
+ {
+ name: "asinh",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Asinh,
+ fn: quat.Asinh,
+ dFn: dAsinh,
+ },
+ // { //fail
+ // name: "acosh",
+ // x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ // fnDual: Acosh,
+ // fn: quat.Acosh,
+ // dFn: dAcosh,
+ // },
+ {
+ name: "atanh",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Atanh,
+ fn: quat.Atanh,
+ dFn: dAtanh,
+ },
+
+ {
+ name: "exp",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Exp,
+ fn: quat.Exp,
+ dFn: dExp,
+ },
+ {
+ name: "log",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Log,
+ fn: quat.Log,
+ dFn: dLog,
+ },
+ {
+ name: "inv",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Inv,
+ fn: quat.Inv,
+ dFn: dInv,
+ },
+ {
+ name: "sqrt",
+ x: []quat.Number{quat.NaN(), quat.Inf(), negThree, negTwo, negOne, negHalf, negZeroQuat, zeroQuat, half, one, two, three},
+ fnDual: Sqrt,
+ fn: quat.Sqrt,
+ dFn: dSqrt,
+ },
+}
+
+func TestDual(t *testing.T) {
+ const tol = 1e-15
+ for _, test := range dualTests {
+ for _, x := range test.x {
+ fxDual := test.fnDual(Number{Real: x, Dual: quat.Number{Real: 1}})
+ fx := test.fn(x)
+ dFx := test.dFn(x)
+ if !same(fxDual.Real, fx, tol) {
+ t.Errorf("unexpected %s(%v): got:%v want:%v", test.name, x, fxDual.Real, fx)
+ }
+ if !same(fxDual.Dual, dFx, tol) {
+ t.Errorf("unexpected %s'(%v): got:%v want:%v", test.name, x, fxDual.Dual, dFx)
+ }
+ }
+ }
+}
+
+/*
+var powRealTests = []struct {
+ d Number
+ p float64
+ want Number
+}{
+ // PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 0}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 3}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 0, want: Number{Real: 1, Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 3}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}},
+
+ // PowReal(x, ±0) = 1 for any x
+ {d: Number{Real: 0, Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: negZero, Emag: 0}, p: negZero, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(1), Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(-1), Emag: 0}, p: negZero, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 0, Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: negZero, Emag: 1}, p: negZero, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(1), Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: math.Inf(-1), Emag: 1}, p: negZero, want: Number{Real: 1, Emag: 0}},
+
+ // PowReal(1+xϵ, y) = (1+xyϵ) for any y
+ {d: Number{Real: 1, Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 0}, p: 1, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 0}, p: 2, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 0}, p: 3, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 1}, p: 1, want: Number{Real: 1, Emag: 1}},
+ {d: Number{Real: 1, Emag: 1}, p: 2, want: Number{Real: 1, Emag: 2}},
+ {d: Number{Real: 1, Emag: 1}, p: 3, want: Number{Real: 1, Emag: 3}},
+ {d: Number{Real: 1, Emag: 2}, p: 0, want: Number{Real: 1, Emag: 0}},
+ {d: Number{Real: 1, Emag: 2}, p: 1, want: Number{Real: 1, Emag: 2}},
+ {d: Number{Real: 1, Emag: 2}, p: 2, want: Number{Real: 1, Emag: 4}},
+ {d: Number{Real: 1, Emag: 2}, p: 3, want: Number{Real: 1, Emag: 6}},
+
+ // PowReal(x, 1) = x for any x
+ {d: Number{Real: 0, Emag: 0}, p: 1, want: Number{Real: 0, Emag: 0}},
+ {d: Number{Real: negZero, Emag: 0}, p: 1, want: Number{Real: negZero, Emag: 0}},
+ {d: Number{Real: 0, Emag: 1}, p: 1, want: Number{Real: 0, Emag: 1}},
+ {d: Number{Real: negZero, Emag: 1}, p: 1, want: Number{Real: negZero, Emag: 1}},
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 1, want: Number{Real: math.NaN(), Emag: 0}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 1, want: Number{Real: math.NaN(), Emag: 1}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 1, want: Number{Real: math.NaN(), Emag: 2}},
+
+ // PowReal(NaN+xϵ, y) = NaN+NaNϵ
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 0}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 1}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: math.NaN(), Emag: 2}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+
+ // PowReal(x, NaN) = NaN+NaNϵ
+ {d: Number{Real: 0, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 0, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 0, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}},
+
+ // Handled by math.Pow tests:
+ //
+ // Pow(±0, y) = ±Inf for y an odd integer < 0
+ // Pow(±0, -Inf) = +Inf
+ // Pow(±0, +Inf) = +0
+ // Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
+ // Pow(±0, y) = ±0 for y an odd integer > 0
+ // Pow(±0, y) = +0 for finite y > 0 and not an odd integer
+ // Pow(-1, ±Inf) = 1
+
+ // PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1
+ {d: Number{Real: 2, Emag: 0}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 0}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+
+ // PowReal(x+yϵ, +Inf) = +Inf for |x| > 1
+ {d: Number{Real: 2, Emag: 1}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+ {d: Number{Real: 3, Emag: 1}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+ {d: Number{Real: 2, Emag: 2}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+ {d: Number{Real: 3, Emag: 2}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}},
+
+ // PowReal(x, -Inf) = +0+NaNϵ for |x| > 1
+ {d: Number{Real: 2, Emag: 0}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 0}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 1}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 1}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 2, Emag: 2}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 3, Emag: 2}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}},
+
+ // PowReal(x+yϵ, +Inf) = +0+NaNϵ for |x| < 1
+ {d: Number{Real: 0.1, Emag: 0}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 0.1, Emag: 0.1}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 0.2, Emag: 0.2}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+ {d: Number{Real: 0.5, Emag: 0.5}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}},
+
+ // PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1
+ {d: Number{Real: 0.1, Emag: 0}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+ {d: Number{Real: 0.2, Emag: 0}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.NaN()}},
+
+ // PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1
+ {d: Number{Real: 0.1, Emag: 0.1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 0.1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.1, Emag: 0.2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 0.2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.1, Emag: 1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.1, Emag: 2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+ {d: Number{Real: 0.2, Emag: 2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}},
+
+ // Handled by math.Pow tests:
+ //
+ // Pow(+Inf, y) = +Inf for y > 0
+ // Pow(+Inf, y) = +0 for y < 0
+ // Pow(-Inf, y) = Pow(-0, -y)
+
+ // PowReal(x, y) = NaN+NaNϵ for finite x < 0 and finite non-integer y
+ {d: Number{Real: -1, Emag: -1}, p: 0.5, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+ {d: Number{Real: -1, Emag: 2}, p: 0.5, want: Number{Real: math.NaN(), Emag: math.NaN()}},
+}
+
+func TestPowReal(t *testing.T) {
+ const tol = 1e-15
+ for _, test := range powRealTests {
+ got := PowReal(test.d, test.p)
+ if !sameDual(got, test.want, tol) {
+ t.Errorf("unexpected PowReal(%v, %v): got:%v want:%v", test.d, test.p, got, test.want)
+ }
+ }
+}
+
+func sameDual(a, b Number, tol float64) bool {
+ return same(a.Real, b.Real, tol) && same(a.Emag, b.Emag, tol)
+}
+*/
+func same(a, b quat.Number, tol float64) bool {
+ return (quat.IsNaN(a) && quat.IsNaN(b)) || (quat.IsInf(a) && quat.IsInf(b)) || equalApprox(a, b, tol)
+}
+
+func equalApprox(a, b quat.Number, tol float64) bool {
+ return floats.EqualWithinAbsOrRel(a.Real, b.Real, tol, tol) &&
+ floats.EqualWithinAbsOrRel(a.Imag, b.Imag, tol, tol) &&
+ floats.EqualWithinAbsOrRel(a.Jmag, b.Jmag, tol, tol) &&
+ floats.EqualWithinAbsOrRel(a.Kmag, b.Kmag, tol, tol)
+}
diff --git a/num/quat/conj.go b/num/quat/conj.go
index 0a15ca8..bcf6d01 100644
--- a/num/quat/conj.go
+++ b/num/quat/conj.go
@@ -15,6 +15,9 @@
// Inv returns the quaternion inverse of q.
func Inv(q Number) Number {
+ if IsInf(q) {
+ return zero
+ }
a := Abs(q)
return Scale(1/(a*a), Conj(q))
}
