Use dep to vendor things

vendor/github.com/go-gl/mathgl/AUTHORS

@@ -0,0 +1,14 @@
+# This is the official list of go-gl/mathgl authors for copyright purposes.
+# This file is distinct from the CONTRIBUTORS files.
+# See the latter for an explanation.
+
+# Names should be added to this file as
+#	Name or Organization <email address>
+# The email address is not required for organizations.
+
+# Please keep the list sorted.
+
+Alex Vasi <eee@someuser.com>
+Dmitri Shuralyov <shurcooL@gmail.com>
+GlenKelley <???>
+Zoe Juozapaitis <jragonmiris@gmail.com>

vendor/github.com/go-gl/mathgl/CONTRIBUTORS

@@ -0,0 +1,21 @@
+# This is the official list of people who can contribute
+# (and typically have contributed) code to the mathgl
+# repository.
+#
+# The AUTHORS file lists the copyright holders; this file
+# lists people.  For example, Google employees would be listed here
+# but not in AUTHORS, because Google would hold the copyright.
+#
+# When adding J Random Contributor's name to this file,
+# either J's name or J's organization's name should be
+# added to the AUTHORS file.
+#
+# Names should be added to this file like so:
+#     Name <email address>
+#
+# Please keep the list sorted.
+
+Alex Vasi <eee@someuser.com>
+Dmitri Shuralyov <shurcooL@gmail.com>
+GlenKelley <???>
+Zoe Juozapaitis <jragonmiris@gmail.com>

vendor/github.com/go-gl/mathgl/LICENSE

@@ -0,0 +1,23 @@
+Copyright ©2013 The go-gl Authors. All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    * Neither the name of the gonum project nor the names of its authors and
+      contributors may be used to endorse or promote products derived from this
+      software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

vendor/github.com/go-gl/mathgl/mgl32/codegen.go

@@ -0,0 +1,276 @@
+// Copyright 2014 The go-gl/mathgl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// +build ignore
+
+// codegen generates go code from templates. Intended to be
+// used with go generate; Also makes mgl64 from mgl32.
+// See the invocation in mgl32/util.go for details.
+// To use it, just run "go generate github.com/go-gl/mathgl/mgl32"
+// (or "go generate" in mgl32 directory).
+
+package main
+
+import (
+	"bytes"
+	"flag"
+	"fmt"
+	"io/ioutil"
+	"os"
+	"os/exec"
+	"path/filepath"
+	"strings"
+	"text/template"
+)
+
+type Context struct {
+	Comment      string
+	TemplateName string
+}
+
+type MatrixIter struct {
+	M     int // row
+	N     int // column
+	index int
+}
+
+var mgl64RewriteRules = []string{
+	"mgl32 -> mgl64",
+	"float32 -> float64",
+	"f32 -> f64",
+	"a.Float32 -> a.Float64",
+	"math.MaxFloat32 -> math.MaxFloat64",
+	"math.SmallestNonzeroFloat32 -> math.SmallestNonzeroFloat64",
+}
+
+func main() {
+	flag.Usage = func() {
+		fmt.Println("Usage: codegen -template file.tmpl -output file.go")
+		fmt.Println("Usage: codegen -mgl64 [-dir ../mgl64]")
+		flag.PrintDefaults()
+	}
+
+	tmplPath := flag.String("template", "file.tmpl", "template path")
+	oPath := flag.String("output", "file.go", "output path")
+	mgl64 := flag.Bool("mgl64", false, "make mgl64")
+	mgl64Path := flag.String("dir", "../mgl64", "path to mgl64 location")
+
+	flag.Parse()
+	if flag.NArg() > 0 || flag.NFlag() == 0 {
+		flag.Usage()
+		os.Exit(2)
+	}
+
+	if *mgl64 {
+		genMgl64(*mgl64Path)
+		return
+	}
+
+	tmpl := template.New("").Delims("<<", ">>").Funcs(template.FuncMap{
+		"typename":    typenameHelper,
+		"elementname": elementNameHelper,
+		"iter":        iterHelper,
+		"matiter":     matrixIterHelper,
+		"enum":        enumHelper,
+		"sep":         separatorHelper,
+		"repeat":      repeatHelper,
+		"add":         addHelper,
+		"mul":         mulHelper,
+	})
+	tmpl = template.Must(tmpl.ParseFiles(*tmplPath))
+	tmplName := filepath.Base(*tmplPath)
+
+	oFile, err := os.Create(*oPath)
+	if err != nil {
+		panic(err)
+	}
+
+	context := Context{
+		Comment:      "This file is generated by codegen.go; DO NOT EDIT",
+		TemplateName: tmplName,
+	}
+	if err = tmpl.ExecuteTemplate(oFile, tmplName, context); err != nil {
+		panic(err)
+	}
+	oFile.Close()
+
+	if err = rungofmt(*oPath, false, nil); err != nil {
+		panic(err)
+	}
+}
+
+func genMgl64(destPath string) {
+	HandleFile := func(source string, info os.FileInfo, err error) error {
+		if err != nil {
+			return err
+		}
+		dest := filepath.Join(destPath, source)
+
+		if info.IsDir() {
+			return os.MkdirAll(dest, info.Mode())
+		}
+		if !strings.HasSuffix(source, ".go") || info.Name() == "codegen.go" {
+			return nil
+		}
+		if !info.Mode().IsRegular() {
+			fmt.Println("Ignored, not a regular file:", source)
+			return nil
+		}
+
+		in, err := ioutil.ReadFile(source)
+		if err != nil {
+			return err
+		}
+
+		out, err := os.OpenFile(dest, os.O_RDWR|os.O_CREATE|os.O_TRUNC,
+			info.Mode())
+		if err != nil {
+			return err
+		}
+		defer out.Close()
+
+		comment := fmt.Sprintf(
+			"// This file is generated from mgl32/%s; DO NOT EDIT\n\n",
+			source)
+		if _, err = out.WriteString(comment); err != nil {
+			return err
+		}
+
+		r := strings.NewReplacer("//go:generate ", "//#go:generate ") // We don't want go generate directives in mgl64 package.
+
+		if _, err = r.WriteString(out, string(in)); err != nil {
+			return err
+		}
+
+		return rungofmt(dest, true, mgl64RewriteRules)
+	}
+
+	if err := filepath.Walk(".", HandleFile); err != nil {
+		panic(err)
+	}
+}
+
+func rungofmt(path string, fiximports bool, rewriteRules []string) error {
+	args := []string{"-w", path}
+	output, err := exec.Command("gofmt", args...).CombinedOutput()
+
+	for i := 0; err == nil && i < len(rewriteRules); i++ {
+		args = []string{"-w", "-r", rewriteRules[i], path}
+		output, err = exec.Command("gofmt", args...).CombinedOutput()
+	}
+
+	if fiximports && err == nil {
+		args = []string{"-w", path}
+		output, err = exec.Command("goimports", args...).CombinedOutput()
+	}
+
+	if err != nil {
+		fmt.Println("Error executing gofmt", strings.Join(args, " "))
+		os.Stdout.Write(output)
+	}
+
+	return err
+}
+
+func typenameHelper(m, n int) string {
+	if m == 1 {
+		return fmt.Sprintf("Vec%d", n)
+	}
+	if n == 1 {
+		return fmt.Sprintf("Vec%d", m)
+	}
+	if m == n {
+		return fmt.Sprintf("Mat%d", m)
+	}
+	return fmt.Sprintf("Mat%dx%d", m, n)
+}
+
+func elementNameHelper(m int) string {
+	switch m {
+	case 0:
+		return "X"
+	case 1:
+		return "Y"
+	case 2:
+		return "Z"
+	case 3:
+		return "W"
+	default:
+		panic("Can't generate element name")
+	}
+}
+
+func iterHelper(start, end int) []int {
+	iter := make([]int, end-start)
+	for i := start; i < end; i++ {
+		iter[i] = i
+	}
+	return iter
+}
+
+func matrixIterHelper(rows, cols int) []MatrixIter {
+	res := make([]MatrixIter, 0, rows*cols)
+
+	for n := 0; n < cols; n++ {
+		for m := 0; m < rows; m++ {
+			res = append(res, MatrixIter{
+				M:     m,
+				N:     n,
+				index: n*rows + m,
+			})
+		}
+	}
+
+	return res
+}
+
+// Template function that returns slice from its arguments. Indended to be used
+// in range loops.
+func enumHelper(args ...int) []int {
+	return args
+}
+
+// Template function to insert commas and '+' in range loops.
+func separatorHelper(sep string, iterCond int) string {
+	if iterCond > 0 {
+		return sep
+	}
+	return ""
+}
+
+// Template function to repeat string 'count' times. Inserting 'sep' between
+// repetitions. Also changes all occurrences of '%d' to repetition number.
+// For example, repeatHelper(3, "col%d", ",") will output "col0, col1, col2"
+func repeatHelper(count int, text string, sep string) string {
+	var res bytes.Buffer
+
+	for i := 0; i < count; i++ {
+		if i > 0 {
+			res.WriteString(sep)
+		}
+		res.WriteString(strings.Replace(text, "%d", fmt.Sprintf("%d", i), -1))
+	}
+
+	return res.String()
+}
+
+func addHelper(args ...int) int {
+	res := 0
+	for _, a := range args {
+		res += a
+	}
+	return res
+}
+
+func mulHelper(args ...int) int {
+	res := 1
+	for _, a := range args {
+		res *= a
+	}
+	return res
+}
+
+func (i MatrixIter) String() string {
+	return fmt.Sprintf("%d", i.index)
+}

vendor/github.com/go-gl/mathgl/mgl32/conv.go

@@ -0,0 +1,94 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import (
+	"math"
+)
+
+// Converts 3-dimensional cartesian coordinates (x,y,z) to spherical
+// coordinates with radius r, inclination theta, and azimuth phi.
+//
+// All angles are in radians.
+func CartesianToSpherical(coord Vec3) (r, theta, phi float32) {
+	r = coord.Len()
+	theta = float32(math.Acos(float64(coord[2] / r)))
+	phi = float32(math.Atan2(float64(coord[1]), float64(coord[0])))
+
+	return
+}
+
+// Converts 3-dimensional cartesian coordinates (x,y,z) to cylindrical
+// coordinates with radial distance r, azimuth phi, and height z.
+//
+// All angles are in radians.
+func CartesianToCylindical(coord Vec3) (rho, phi, z float32) {
+	rho = float32(math.Hypot(float64(coord[0]), float64(coord[1])))
+
+	phi = float32(math.Atan2(float64(coord[1]), float64(coord[0])))
+
+	z = coord[2]
+
+	return
+}
+
+// Converts spherical coordinates with radius r, inclination theta,
+// and azimuth phi to cartesian coordinates (x,y,z).
+//
+// Angles are in radians.
+func SphericalToCartesian(r, theta, phi float32) Vec3 {
+	st, ct := math.Sincos(float64(theta))
+	sp, cp := math.Sincos(float64(phi))
+
+	return Vec3{r * float32(st*cp), r * float32(st*sp), r * float32(ct)}
+}
+
+// Converts spherical coordinates with radius r, inclination theta,
+// and azimuth phi to cylindrical coordinates with radial distance r,
+// azimuth phi, and height z.
+//
+// Angles are in radians
+func SphericalToCylindrical(r, theta, phi float32) (rho, phi2, z float32) {
+	s, c := math.Sincos(float64(theta))
+
+	rho = r * float32(s)
+	z = r * float32(c)
+	phi2 = phi
+
+	return
+}
+
+// Converts cylindrical coordinates with radial distance r,
+// azimuth phi, and height z to spherical coordinates with radius r,
+// inclination theta, and azimuth phi.
+//
+// Angles are in radians
+func CylindircalToSpherical(rho, phi, z float32) (r, theta, phi2 float32) {
+	r = float32(math.Hypot(float64(rho), float64(z)))
+	phi2 = phi
+	theta = float32(math.Atan2(float64(rho), float64(z)))
+
+	return
+}
+
+// Converts cylindrical coordinates with radial distance r,
+// azimuth phi, and height z to cartesian coordinates (x,y,z)
+//
+// Angles are in radians.
+func CylindricalToCartesian(rho, phi, z float32) Vec3 {
+	s, c := math.Sincos(float64(phi))
+
+	return Vec3{rho * float32(c), rho * float32(s), z}
+}
+
+// Converts degrees to radians
+func DegToRad(angle float32) float32 {
+	return angle * float32(math.Pi) / 180
+}
+
+// Converts radians to degrees
+func RadToDeg(angle float32) float32 {
+	return angle * 180 / float32(math.Pi)
+}

vendor/github.com/go-gl/mathgl/mgl32/doc.go

@@ -0,0 +1,22 @@
+// Copyright 2014 The go-gl 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 mgl[32|64] (an abbreviation of mathgl since the packages were split between 32 and 64-bit versions)
+is a pure Go math package specialized for 3D math, with inspiration from GLM. It provides statically-sized vectors and matrices with
+compile-time generated calculations for most basic math operations. It also provides several basic graphics utilities such as bezier curves and surfaces,
+generation of basic primitives like circles, easy creation of common matrices such as perspective or rotation, and common operations like converting
+to/from screen/OpenGL coordinates or Projecting/Unprojecting from an MVP matrix. Quaternions are also supported.
+
+The basic vectors and matrices are written with code generation, so looking directly at the source will probably be a bit confusing. I recommend looking at the Godoc
+instead, as all basic functions are documented.
+
+This package is written in Column Major Order to make it easier with OpenGL. This means for uniform blocks you can use the default ordering, and when you call
+pass-in functions you can leave the "transpose" argument as false.
+
+The package now contains variable sized vectors and matrices. Using these is discouraged. They exist for corner cases where you need "small" matrices that are still
+bigger than 4x4. An example may be a Jacobean used for inverse kinematics. Things like computer vision or general linear algebra are best left to packages
+more directly suited for that task -- OpenCV, BLAS, LAPACK, numpy, gonum (if you want to stay in Go), and so on.
+*/
+package mgl32

vendor/github.com/go-gl/mathgl/mgl32/matmn.go

@@ -0,0 +1,494 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import (
+	"math"
+)
+
+// An arbitrary mxn matrix backed by a slice of floats.
+//
+// This is emphatically not recommended for hardcore n-dimensional
+// linear algebra. For that purpose I recommend github.com/gonum/matrix or
+// well-tested C libraries such as BLAS or LAPACK.
+//
+// This is meant to complement future algorithms that may require matrices larger than
+// 4x4, but still relatively small (e.g. Jacobeans for inverse kinematics).
+//
+// It makes use of the same memory sync.Pool set that VecN does, with the same sizing rules.
+//
+// MatMN will always check if the receiver is nil on any method. Meaning MathMN(nil).Add(dst,m2)
+// should always work. Except for the Reshape function, the semantics of this is to "propogate" nils
+// forward, so if an invalid operation occurs in a long chain of matrix operations, the overall result will be nil.
+type MatMxN struct {
+	m, n int
+	dat  []float32
+}
+
+// Creates a matrix backed by a new slice of size m*n
+func NewMatrix(m, n int) (mat *MatMxN) {
+	if shouldPool {
+		return &MatMxN{m: m, n: n, dat: grabFromPool(m * n)}
+	} else {
+		return &MatMxN{m: m, n: n, dat: make([]float32, m*n)}
+	}
+}
+
+// Returns a matrix with data specified by the data in src
+//
+// For instance, to create a 3x3 MatMN from a Mat3
+//
+//    m1 := mgl32.Rotate3DX(3.14159)
+//    mat := mgl32.NewBackedMatrix(m1[:],3,3)
+//
+// will create an MN matrix matching the data in the original
+// rotation matrix. This matrix is NOT backed by the initial slice;
+// it's a copy of the data
+//
+// If m*n > cap(src), this function will panic.
+func NewMatrixFromData(src []float32, m, n int) *MatMxN {
+	var internal []float32
+	if shouldPool {
+		internal = grabFromPool(m * n)
+	} else {
+		internal = make([]float32, m*n)
+	}
+	copy(internal, src[:m*n])
+
+	return &MatMxN{m: m, n: n, dat: internal}
+}
+
+// Copies src into dst. This Reshapes dst
+// to the same size as src.
+//
+// If dst or src is nil, this is a no-op
+func CopyMatMN(dst, src *MatMxN) {
+	if dst == nil || src == nil {
+		return
+	}
+	dst.Reshape(src.m, src.n)
+	copy(dst.dat, src.dat)
+}
+
+// Stores the NxN identity matrix in dst, reallocating as necessary.
+func IdentN(dst *MatMxN, n int) *MatMxN {
+	dst = dst.Reshape(n, n)
+
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			if i == j {
+				dst.Set(i, j, 1)
+			} else {
+				dst.Set(i, j, 0)
+			}
+		}
+	}
+
+	return dst
+}
+
+// Creates an NxN diagonal matrix seeded by the diagonal vector
+// diag. Meaning: for all entries, where i==j, dst.At(i,j) = diag[i]. Otherwise
+// dst.At(i,j) = 0
+//
+// This reshapes dst to the correct size, returning/grabbing from the memory pool as necessary.
+func DiagN(dst *MatMxN, diag *VecN) *MatMxN {
+	dst = dst.Reshape(len(diag.vec), len(diag.vec))
+	n := len(diag.vec)
+
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			if i == j {
+				dst.Set(i, j, diag.vec[i])
+			} else {
+				dst.Set(i, j, 0)
+			}
+		}
+	}
+
+	return dst
+}
+
+// Reshapes the matrix to m by n and zeroes out all
+// elements.
+func (mat *MatMxN) Zero(m, n int) {
+	if mat == nil {
+		return
+	}
+
+	mat.Reshape(m, n)
+	for i := range mat.dat {
+		mat.dat[i] = 0
+	}
+}
+
+// Returns the underlying matrix slice to the memory pool
+func (mat *MatMxN) destroy() {
+	if mat == nil {
+		return
+	}
+
+	if shouldPool && mat.dat != nil {
+		returnToPool(mat.dat)
+	}
+	mat.m, mat.n = 0, 0
+	mat.dat = nil
+}
+
+// Reshapes the matrix to the desired dimensions.
+// If the overall size of the new matrix (m*n) is bigger
+// than the current size, the underlying slice will
+// be grown, sending the current slice to the memory pool
+// and grabbing a bigger one if necessary
+//
+// If the caller is a nil pointer, the return value will be a new
+// matrix, as if NewMatrix(m,n) had been called. Otherwise it's
+// simply the caller.
+func (mat *MatMxN) Reshape(m, n int) *MatMxN {
+	if mat == nil {
+		return NewMatrix(m, n)
+	}
+
+	if m*n <= cap(mat.dat) {
+		if mat.dat != nil {
+			mat.dat = mat.dat[:m*n]
+		} else {
+			mat.dat = []float32{}
+		}
+		mat.m, mat.n = m, n
+		return mat
+	}
+
+	if shouldPool && mat.dat != nil {
+		returnToPool(mat.dat)
+	}
+	(*mat) = (*NewMatrix(m, n))
+
+	return mat
+}
+
+// Infers an MxN matrix from a constant matrix from this package. For instance,
+// a Mat2x3 inferred with this function will work just like NewMatrixFromData(m[:],2,3)
+// where m is the Mat2x3. This uses a type switch.
+//
+// I personally recommend using NewMatrixFromData, because it avoids a potentially costly type switch.
+// However, this is also more robust and less error prone if you change the size of your matrix somewhere.
+//
+// If the value passed in is not recognized, it returns an InferMatrixError.
+func (mat *MatMxN) InferMatrix(m interface{}) (*MatMxN, error) {
+	switch raw := m.(type) {
+	case Mat2:
+		return NewMatrixFromData(raw[:], 2, 2), nil
+	case Mat2x3:
+		return NewMatrixFromData(raw[:], 2, 3), nil
+	case Mat2x4:
+		return NewMatrixFromData(raw[:], 2, 4), nil
+	case Mat3:
+		return NewMatrixFromData(raw[:], 3, 3), nil
+	case Mat3x2:
+		return NewMatrixFromData(raw[:], 3, 2), nil
+	case Mat3x4:
+		return NewMatrixFromData(raw[:], 3, 4), nil
+	case Mat4:
+		return NewMatrixFromData(raw[:], 4, 4), nil
+	case Mat4x2:
+		return NewMatrixFromData(raw[:], 4, 2), nil
+	case Mat4x3:
+		return NewMatrixFromData(raw[:], 4, 3), nil
+	default:
+		return nil, InferMatrixError{}
+	}
+}
+
+// Returns the trace of a square matrix (sum of all diagonal elements). If the matrix
+// is nil, or not square, the result will be NaN.
+func (mat *MatMxN) Trace() float32 {
+	if mat == nil || mat.m != mat.n {
+		return float32(math.NaN())
+	}
+
+	var out float32
+	for i := 0; i < mat.m; i++ {
+		out += mat.At(i, i)
+	}
+
+	return out
+}
+
+// Takes the transpose of mat and puts it in dst.
+//
+// If dst is not of the correct dimensions, it will be Reshaped,
+// if dst and mat are the same, a temporary matrix of the correct size will
+// be allocated; these resources will be released via the memory pool.
+//
+// This should be improved in the future.
+func (mat *MatMxN) Transpose(dst *MatMxN) (t *MatMxN) {
+	if mat == nil {
+		return nil
+	}
+
+	if dst == mat {
+		dst = NewMatrix(mat.n, mat.m)
+
+		// Copy data to correct matrix,
+		// delete temporary buffer,
+		// and set the return value to the
+		// correct one
+		defer func() {
+			copy(mat.dat, dst.dat)
+
+			mat.m, mat.n = mat.n, mat.m
+
+			dst.destroy()
+			t = mat
+		}()
+
+		return mat
+	} else {
+		dst = dst.Reshape(mat.n, mat.m)
+	}
+
+	for r := 0; r < mat.m; r++ {
+		for c := 0; c < mat.n; c++ {
+			dst.dat[r*dst.m+c] = mat.dat[c*mat.m+r]
+		}
+	}
+
+	return dst
+}
+
+// Returns the raw slice backing this matrix
+func (mat *MatMxN) Raw() []float32 {
+	if mat == nil {
+		return nil
+	}
+
+	return mat.dat
+}
+
+// Returns the number of rows in this matrix
+func (mat *MatMxN) NumRows() int {
+	return mat.m
+}
+
+// Returns the number of columns in this matrix
+func (mat *MatMxN) NumCols() int {
+	return mat.n
+}
+
+// Returns the number of rows and columns in this matrix
+// as a single operation
+func (mat *MatMxN) NumRowCols() (rows, cols int) {
+	return mat.m, mat.n
+}
+
+// Returns the element at the given row and column.
+// This is garbage in/garbage out and does no bounds
+// checking. If the computation happens to lead to an invalid
+// element, it will be returned; or it may panic.
+func (mat *MatMxN) At(row, col int) float32 {
+	return mat.dat[col*mat.m+row]
+}
+
+// Sets the element at the given row and column.
+// This is garbage in/garbage out and does no bounds
+// checking. If the computation happens to lead to an invalid
+// element, it will be set; or it may panic.
+func (mat *MatMxN) Set(row, col int, val float32) {
+	mat.dat[col*mat.m+row] = val
+}
+
+func (mat *MatMxN) Add(dst *MatMxN, addend *MatMxN) *MatMxN {
+	if mat == nil || addend == nil || mat.m != addend.m || mat.n != addend.n {
+		return nil
+	}
+
+	dst = dst.Reshape(mat.m, mat.n)
+
+	// No need to care about rows and columns
+	// since it's element-wise anyway
+	for i, el := range mat.dat {
+		dst.dat[i] = el + addend.dat[i]
+	}
+
+	return dst
+}
+
+func (mat *MatMxN) Sub(dst *MatMxN, subtrahend *MatMxN) *MatMxN {
+	if mat == nil || subtrahend == nil || mat.m != subtrahend.m || mat.n != subtrahend.n {
+		return nil
+	}
+
+	dst = dst.Reshape(mat.m, mat.n)
+
+	// No need to care about rows and columns
+	// since it's element-wise anyway
+	for i, el := range mat.dat {
+		dst.dat[i] = el - subtrahend.dat[i]
+	}
+
+	return dst
+}
+
+// Performs matrix multiplication on MxN matrix mat and NxO matrix mul, storing the result in dst.
+// This returns dst, or nil if the operation is not able to be performed.
+//
+// If mat == dst, or mul == dst a temporary matrix will be used.
+//
+// This uses the naive algorithm (though on smaller matrices,
+// this can actually be faster; about len(mat)+len(mul) < ~100)
+func (mat *MatMxN) MulMxN(dst *MatMxN, mul *MatMxN) *MatMxN {
+	if mat == nil || mul == nil || mat.n != mul.m {
+		return nil
+	}
+
+	if dst == mul {
+		mul = NewMatrix(mul.m, mul.n)
+		copy(mul.dat, dst.dat)
+
+		// If mat==dst==mul, we need to change
+		// mat too or we have a bug
+		if mat == dst {
+			mat = mul
+		}
+
+		defer mul.destroy()
+	} else if dst == mat {
+		mat = NewMatrix(mat.m, mat.n)
+		copy(mat.dat, dst.dat)
+
+		defer mat.destroy()
+	}
+
+	dst = dst.Reshape(mat.m, mul.n)
+	for r1 := 0; r1 < mat.m; r1++ {
+		for c2 := 0; c2 < mul.n; c2++ {
+
+			dst.dat[c2*mat.m+r1] = 0
+			for i := 0; i < mat.n; i++ {
+				dst.dat[c2*mat.m+r1] += mat.dat[i*mat.m+r1] * mul.dat[c2*mul.m+i]
+			}
+
+		}
+	}
+
+	return dst
+}
+
+// Performs a scalar multiplication between mat and some constant c,
+// storing the result in dst. Mat and dst can be equal. If dst is not the
+// correct size, a Reshape will occur.
+func (mat *MatMxN) Mul(dst *MatMxN, c float32) *MatMxN {
+	if mat == nil {
+		return nil
+	}
+
+	dst = dst.Reshape(mat.m, mat.n)
+
+	for i, el := range mat.dat {
+		dst.dat[i] = el * c
+	}
+
+	return dst
+}
+
+// Multiplies the matrix by a vector of size n. If mat or v is
+// nil, this returns nil. If the number of columns in mat does not match
+// the Size of v, this also returns nil.
+//
+// Dst will be resized if it's not big enough. If dst == v; a temporary
+// vector will be allocated and returned via the realloc callback when complete.
+func (mat *MatMxN) MulNx1(dst, v *VecN) *VecN {
+	if mat == nil || v == nil || mat.n != len(v.vec) {
+		return nil
+	}
+	if dst == v {
+		v = NewVecN(len(v.vec))
+		copy(v.vec, dst.vec)
+
+		defer v.destroy()
+	}
+
+	dst = dst.Resize(mat.m)
+
+	for r := 0; r < mat.m; r++ {
+		dst.vec[r] = 0
+
+		for c := 0; c < mat.n; c++ {
+			dst.vec[r] += mat.At(r, c) * v.vec[c]
+		}
+	}
+
+	return dst
+}
+
+func (mat *MatMxN) ApproxEqual(m2 *MatMxN) bool {
+	if mat == m2 {
+		return true
+	}
+	if mat.m != m2.m || mat.n != m2.n {
+		return false
+	}
+
+	for i, el := range mat.dat {
+		if !FloatEqual(el, m2.dat[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (mat *MatMxN) ApproxEqualThreshold(m2 *MatMxN, epsilon float32) bool {
+	if mat == m2 {
+		return true
+	}
+	if mat.m != m2.m || mat.n != m2.n {
+		return false
+	}
+
+	for i, el := range mat.dat {
+		if !FloatEqualThreshold(el, m2.dat[i], epsilon) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (mat *MatMxN) ApproxEqualFunc(m2 *MatMxN, comp func(float32, float32) bool) bool {
+	if mat == m2 {
+		return true
+	}
+	if mat.m != m2.m || mat.n != m2.n {
+		return false
+	}
+
+	for i, el := range mat.dat {
+		if !comp(el, m2.dat[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+type InferMatrixError struct{}
+
+func (me InferMatrixError) Error() string {
+	return "could not infer matrix. Make sure you're using a constant matrix such as Mat3 from within the same package (meaning: mgl32.MatMxN can't handle a mgl64.Mat2x3)."
+}
+
+type RectangularMatrixError struct{}
+
+func (mse RectangularMatrixError) Error() string {
+	return "the matrix was the wrong shape, needed a square matrix."
+}
+
+type NilMatrixError struct{}
+
+func (me NilMatrixError) Error() string {
+	return "the matrix is nil"
+}

vendor/github.com/go-gl/mathgl/mgl32/matrix.tmpl

@@ -0,0 +1,370 @@
+// Copyright 2014 The go-gl/mathgl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// <<.Comment>>
+// Edit <<.TemplateName>> and run "go generate" to make changes.
+
+package mgl32
+
+import (
+	"bytes"
+	"fmt"
+	"golang.org/x/image/math/f32"
+	"text/tabwriter"
+)
+
+type Mat2 [4]float32
+type Mat2x3 [6]float32
+type Mat2x4 [8]float32
+type Mat3x2 [6]float32
+type Mat3 f32.Mat3
+type Mat3x4 [12]float32
+type Mat4x2 [8]float32
+type Mat4x3 [12]float32
+type Mat4 f32.Mat4
+
+func (m Mat2) Mat3() Mat3 {
+	col0, col1 := m.Cols()
+	return Mat3FromCols(
+		col0.Vec3(0),
+		col1.Vec3(0),
+		Vec3{0, 0, 1},
+	)
+}
+
+func (m Mat2) Mat4() Mat4 {
+	col0, col1 := m.Cols()
+	return Mat4FromCols(
+		col0.Vec4(0, 0),
+		col1.Vec4(0, 0),
+		Vec4{0, 0, 1, 0},
+		Vec4{0, 0, 0, 1},
+	)
+}
+
+func (m Mat3) Mat2() Mat2 {
+	col0, col1, _ := m.Cols()
+	return Mat2FromCols(
+		col0.Vec2(),
+		col1.Vec2(),
+	)
+}
+
+func (m Mat3) Mat4() Mat4 {
+	col0, col1, col2 := m.Cols()
+	return Mat4FromCols(
+		col0.Vec4(0),
+		col1.Vec4(0),
+		col2.Vec4(0),
+		Vec4{0, 0, 0, 1},
+	)
+}
+
+func (m Mat4) Mat2() Mat2 {
+	col0, col1, _, _ := m.Cols()
+	return Mat2FromCols(
+		col0.Vec2(),
+		col1.Vec2(),
+	)
+}
+
+func (m Mat4) Mat3() Mat3 {
+	col0, col1, col2, _ := m.Cols()
+	return Mat3FromCols(
+		col0.Vec3(),
+		col1.Vec3(),
+		col2.Vec3(),
+	)
+}
+
+
+<</* Common functions for all matrices */>>
+<<range $m := enum 2 3 4>><<range $n := enum 2 3 4>>
+<<$type := typename $m $n>>
+
+// Sets a Column within the Matrix, so it mutates the calling matrix.
+func (m *<<$type>>) SetCol(col int, v <<typename $m 1>>) {
+	<<range $i := iter 0 $m>><<sep "," $i>>m[col*<<$m>>+<<$i>>]<<end>> = <<repeat $m "v[%d]" ",">>
+}
+
+// Sets a Row within the Matrix, so it mutates the calling matrix.
+func (m *<<$type>>) SetRow(row int, v <<typename $n 1>>) {
+	<<range $i := iter 0 $n>><<sep "," $i>>m[row+<<mul $m $i>>]<<end>> = <<repeat $n "v[%d]" ",">>
+}
+
+<<if eq $m $n>>
+// Diag is a basic operation on a square matrix that simply
+// returns main diagonal (meaning all elements such that row==col).
+func (m <<$type>>) Diag() <<typename $m 1>> {
+	return <<typename $m 1>>{<<range $i := iter 0 $m>>m[<<mul $i $m | add $i>>],<<end>>}
+}
+<<end>>
+
+<<if eq $m $n>>
+// Ident<<$m>> returns the <<$m>>x<<$n>> identity matrix.
+// The identity matrix is a square matrix with the value 1 on its
+// diagonals. The characteristic property of the identity matrix is that
+// any matrix multiplied by it is itself. (MI = M; IN = N)
+func Ident<<$m>>() <<$type>> {
+	return <<$type>>{<<range $i := matiter $m $n>><<if eq $i.M $i.N>>1<<else>>0<<end>>,<<end>>}
+}
+<<end>>
+
+<<if eq $m $n>>
+// Diag<<$m>> creates a diagonal matrix from the entries of the input vector.
+// That is, for each pointer for row==col, vector[row] is the entry. Otherwise it's 0.
+//
+// Another way to think about it is that the identity is this function where the every vector element is 1.
+func Diag<<$m>>(v <<typename $m 1>>) <<$type>> {
+	var m <<$type>>
+	<<range $i := iter 0 $m>><<sep "," $i>>m[<<mul $i $m | add $i>>]<<end>> = <<repeat $m "v[%d]" ",">>
+	return m
+}
+<<end>>
+
+// <<$type>>FromRows builds a new matrix from row vectors.
+// The resulting matrix will still be in column major order, but this can be
+// good for hand-building matrices.
+
+func <<$type>>FromRows(<<range $i := iter 0 $m>><<sep "," $i>>row<<$i>><<end>> <<typename $n 1>>) <<$type>> {
+	return <<$type>>{<<range $i := matiter $m $n>>row<<$i.M>>[<<$i.N>>],<<end>>}
+}
+
+// <<$type>>FromCols builds a new matrix from column vectors.
+func <<$type>>FromCols(<<repeat $n "col%d" ",">> <<typename $m 1>>) <<$type>> {
+	return <<$type>>{<<range $i := matiter $m $n>>col<<$i.N>>[<<$i.M>>], <<end>>}
+}
+
+// Add performs an element-wise addition of two matrices, this is
+// equivalent to iterating over every element of m1 and adding the corresponding value of m2.
+func (m1 <<$type>>) Add(m2 <<$type>>) <<$type>> {
+	return <<$type>>{<< range $i := matiter $m $n>>m1[<<$i>>] + m2[<<$i>>], <<end>>}
+}
+
+// Sub performs an element-wise subtraction of two matrices, this is
+// equivalent to iterating over every element of m1 and subtracting the corresponding value of m2.
+func (m1 <<$type>>) Sub(m2 <<$type>>) <<$type>> {
+	return <<$type>>{<< range $i := matiter $m $n>>m1[<<$i>>] - m2[<<$i>>], <<end>>}
+}
+
+// Mul performs a scalar multiplcation of the matrix. This is equivalent to iterating
+// over every element of the matrix and multiply it by c.
+func (m1 <<$type>>) Mul(c float32) <<$type>> {
+	return <<$type>>{<< range $i := matiter $m $n>>m1[<<$i>>] * c, <<end>>}
+}
+
+<<range $o := enum 1 2 3 4>>
+// Mul<<$n>><<if ne $n $o>>x<<$o>><<end>> performs a "matrix product" between this matrix
+// and another of the given dimension. For any two matrices of dimensionality
+// MxN and NxO, the result will be MxO. For instance, Mat4 multiplied using
+// Mul4x2 will result in a Mat4x2.
+func (m1 <<$type>>) Mul<<$n>><<if ne $n $o>>x<<$o>><<end>>(m2 <<typename $n $o>>) <<typename $m $o>> {
+	return <<typename $m $o>>{<<range $i := matiter $m $o>>
+		<<range $k := iter 0 $n>><<sep "+" $k>>m1[<<mul $k $m | add $i.M>>]*m2[<<mul $i.N $n| add $k>>]<<end>>,<<end>>
+	}
+}
+<<end>>
+
+// Transpose produces the transpose of this matrix. For any MxN matrix
+// the transpose is an NxM matrix with the rows swapped with the columns. For instance
+// the transpose of the Mat3x2 is a Mat2x3 like so:
+//
+//    [[a b]]    [[a c e]]
+//    [[c d]] =  [[b d f]]
+//    [[e f]]
+func (m1 <<$type>>) Transpose() <<typename $n $m>> {
+	return <<typename $n $m>>{<<range $i := matiter $n $m>>m1[<<mul $m $i.M | add $i.N>>], <<end>>}
+}
+
+<<if eq $m $n>>
+// The determinant of a matrix is a measure of a square matrix's
+// singularity and invertability, among other things. In this library, the
+// determinant is hard coded based on pre-computed cofactor expansion, and uses
+// no loops. Of course, the addition and multiplication must still be done.
+func (m <<$type>>) Det() float32 {
+	<<if eq $m 2 ->>
+	return m[0]*m[3] - m[1]*m[2]
+	<<else if eq $m 3 ->>
+	return m[0]*m[4]*m[8] + m[3]*m[7]*m[2] + m[6]*m[1]*m[5] - m[6]*m[4]*m[2] - m[3]*m[1]*m[8] - m[0]*m[7]*m[5]
+	<<else if eq $m 4 ->>
+	return m[0]*m[5]*m[10]*m[15] - m[0]*m[5]*m[11]*m[14] - m[0]*m[6]*m[9]*m[15] + m[0]*m[6]*m[11]*m[13] + m[0]*m[7]*m[9]*m[14] - m[0]*m[7]*m[10]*m[13] - m[1]*m[4]*m[10]*m[15] + m[1]*m[4]*m[11]*m[14] + m[1]*m[6]*m[8]*m[15] - m[1]*m[6]*m[11]*m[12] - m[1]*m[7]*m[8]*m[14] + m[1]*m[7]*m[10]*m[12] + m[2]*m[4]*m[9]*m[15] - m[2]*m[4]*m[11]*m[13] - m[2]*m[5]*m[8]*m[15] + m[2]*m[5]*m[11]*m[12] + m[2]*m[7]*m[8]*m[13] - m[2]*m[7]*m[9]*m[12] - m[3]*m[4]*m[9]*m[14] + m[3]*m[4]*m[10]*m[13] + m[3]*m[5]*m[8]*m[14] - m[3]*m[5]*m[10]*m[12] - m[3]*m[6]*m[8]*m[13] + m[3]*m[6]*m[9]*m[12]
+	<<end ->>
+}
+<<end>>
+
+<<if eq $m $n>>
+// Inv computes the inverse of a square matrix. An inverse is a square matrix such that when multiplied by the
+// original, yields the identity.
+//
+// M_inv * M = M * M_inv = I
+//
+// In this library, the math is precomputed, and uses no loops, though the multiplications, additions, determinant calculation, and scaling
+// are still done. This can still be (relatively) expensive for a 4x4.
+//
+// This function checks the determinant to see if the matrix is invertible.
+// If the determinant is 0.0, this function returns the zero matrix. However, due to floating point errors, it is
+// entirely plausible to get a false positive or negative.
+// In the future, an alternate function may be written which takes in a pre-computed determinant.
+func (m <<$type>>) Inv() <<$type>> {
+	det := m.Det()
+	if FloatEqual(det, float32(0.0)) {
+		return <<$type>>{}
+	}
+	<<if eq $m 2>>
+	retMat := Mat2{m[3], -m[1], -m[2], m[0]}
+	<<else if eq $m 3>>
+	retMat := Mat3{
+		m[4]*m[8] - m[5]*m[7],
+		m[2]*m[7] - m[1]*m[8],
+		m[1]*m[5] - m[2]*m[4],
+		m[5]*m[6] - m[3]*m[8],
+		m[0]*m[8] - m[2]*m[6],
+		m[2]*m[3] - m[0]*m[5],
+		m[3]*m[7] - m[4]*m[6],
+		m[1]*m[6] - m[0]*m[7],
+		m[0]*m[4] - m[1]*m[3],
+	}
+	<<else if eq $m 4>>
+	retMat := Mat4{
+		-m[7]*m[10]*m[13] + m[6]*m[11]*m[13] + m[7]*m[9]*m[14] - m[5]*m[11]*m[14] - m[6]*m[9]*m[15] + m[5]*m[10]*m[15],
+		 m[3]*m[10]*m[13] - m[2]*m[11]*m[13] - m[3]*m[9]*m[14] + m[1]*m[11]*m[14] + m[2]*m[9]*m[15] - m[1]*m[10]*m[15],
+		 -m[3]*m[6]*m[13] + m[2]*m[7]*m[13] + m[3]*m[5]*m[14] - m[1]*m[7]*m[14] - m[2]*m[5]*m[15] + m[1]*m[6]*m[15],
+		 m[3]*m[6]*m[9] - m[2]*m[7]*m[9] - m[3]*m[5]*m[10] + m[1]*m[7]*m[10] + m[2]*m[5]*m[11] - m[1]*m[6]*m[11],
+		 m[7]*m[10]*m[12] - m[6]*m[11]*m[12] - m[7]*m[8]*m[14] + m[4]*m[11]*m[14] + m[6]*m[8]*m[15] - m[4]*m[10]*m[15],
+		 -m[3]*m[10]*m[12] + m[2]*m[11]*m[12] + m[3]*m[8]*m[14] - m[0]*m[11]*m[14] - m[2]*m[8]*m[15] + m[0]*m[10]*m[15],
+		 m[3]*m[6]*m[12] - m[2]*m[7]*m[12] - m[3]*m[4]*m[14] + m[0]*m[7]*m[14] + m[2]*m[4]*m[15] - m[0]*m[6]*m[15],
+		 -m[3]*m[6]*m[8] + m[2]*m[7]*m[8] + m[3]*m[4]*m[10] - m[0]*m[7]*m[10] - m[2]*m[4]*m[11] + m[0]*m[6]*m[11],
+		 -m[7]*m[9]*m[12] + m[5]*m[11]*m[12] + m[7]*m[8]*m[13] - m[4]*m[11]*m[13] - m[5]*m[8]*m[15] + m[4]*m[9]*m[15],
+		 m[3]*m[9]*m[12] - m[1]*m[11]*m[12] - m[3]*m[8]*m[13] + m[0]*m[11]*m[13] + m[1]*m[8]*m[15] - m[0]*m[9]*m[15],
+		 -m[3]*m[5]*m[12] + m[1]*m[7]*m[12] + m[3]*m[4]*m[13] - m[0]*m[7]*m[13] - m[1]*m[4]*m[15] + m[0]*m[5]*m[15],
+		 m[3]*m[5]*m[8] - m[1]*m[7]*m[8] - m[3]*m[4]*m[9] + m[0]*m[7]*m[9] + m[1]*m[4]*m[11] - m[0]*m[5]*m[11],
+		 m[6]*m[9]*m[12] - m[5]*m[10]*m[12] - m[6]*m[8]*m[13] + m[4]*m[10]*m[13] + m[5]*m[8]*m[14] - m[4]*m[9]*m[14],
+		 -m[2]*m[9]*m[12] + m[1]*m[10]*m[12] + m[2]*m[8]*m[13] - m[0]*m[10]*m[13] - m[1]*m[8]*m[14] + m[0]*m[9]*m[14],
+		 m[2]*m[5]*m[12] - m[1]*m[6]*m[12] - m[2]*m[4]*m[13] + m[0]*m[6]*m[13] + m[1]*m[4]*m[14] - m[0]*m[5]*m[14],
+		-m[2]*m[5]*m[8] + m[1]*m[6]*m[8] + m[2]*m[4]*m[9] - m[0]*m[6]*m[9] - m[1]*m[4]*m[10] + m[0]*m[5]*m[10],
+	}
+	<<end>>
+	return retMat.Mul(1 / det)
+}
+<<end>>
+
+// ApproxEqual performs an element-wise approximate equality test between two matrices,
+// as if FloatEqual had been used.
+func (m1 <<$type>>) ApproxEqual(m2 <<$type>>) bool {
+	for i := range m1 {
+		if !FloatEqual(m1[i], m2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxEqualThreshold performs an element-wise approximate equality test between two matrices
+// with a given epsilon threshold, as if FloatEqualThreshold had been used.
+func (m1 <<$type>>) ApproxEqualThreshold(m2 <<$type>>, threshold float32) bool {
+	for i := range m1 {
+		if !FloatEqualThreshold(m1[i], m2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxEqualFunc performs an element-wise approximate equality test between two matrices
+// with a given equality functions, intended to be used with FloatEqualFunc; although and comparison
+// function may be used in practice.
+func (m1 <<$type>>) ApproxFuncEqual(m2 <<$type>>, eq func(float32, float32) bool) bool {
+	for i := range m1 {
+		if !eq(m1[i], m2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// At returns the matrix element at the given row and column.
+// This is equivalent to mat[col * numRow + row] where numRow is constant
+// (E.G. for a Mat3x2 it's equal to 3)
+//
+// This method is garbage-in garbage-out. For instance, on a Mat4 asking for
+// At(5,0) will work just like At(1,1). Or it may panic if it's out of bounds.
+func (m <<$type>>) At(row, col int) float32 {
+	return m[col*<<$m>>+row]
+}
+
+// Set sets the corresponding matrix element at the given row and column.
+// This has a pointer receiver because it mutates the matrix.
+//
+// This method is garbage-in garbage-out. For instance, on a Mat4 asking for
+// Set(5,0,val) will work just like Set(1,1,val). Or it may panic if it's out of bounds.
+func (m *<<$type>>) Set(row, col int, value float32) {
+	m[col*<<$m>>+row] = value
+}
+
+// Index returns the index of the given row and column, to be used with direct
+// access. E.G. Index(0,0) = 0.
+//
+// This is a garbage-in garbage-out method. For instance, on a Mat4 asking for the index of
+// (5,0) will work the same as asking for (1,1). Or it may give you a value that will cause
+// a panic if you try to access the array with it if it's truly out of bounds.
+func (m <<$type>>) Index(row, col int) int {
+	return col*<<$m>> + row
+}
+
+// Row returns a vector representing the corresponding row (starting at row 0).
+// This package makes no distinction between row and column vectors, so it
+// will be a normal VecM for a MxN matrix.
+func (m <<$type>>) Row(row int) <<typename $n 1>> {
+	return <<typename $n 1>>{<<range $i := iter 0 $n>>m[row+<<mul $m $i>>],<<end>>}
+}
+
+// Rows decomposes a matrix into its corresponding row vectors.
+// This is equivalent to calling mat.Row for each row.
+func (m <<$type>>) Rows() (<<repeat $m "row%d" ",">> <<typename $n 1>>) {
+	return <<repeat $m "m.Row(%d)" ",">>
+}
+
+// Col returns a vector representing the corresponding column (starting at col 0).
+// This package makes no distinction between row and column vectors, so it
+// will be a normal VecN for a MxN matrix.
+func (m <<$type>>) Col(col int) <<typename $m 1>> {
+	return <<typename $m 1>>{<<range $i := iter 0 $m>>m[col*<<$m>>+<<$i>>],<<end>>}
+}
+
+// Cols decomposes a matrix into its corresponding column vectors.
+// This is equivalent to calling mat.Col for each column.
+func (m <<$type>>) Cols() (<<repeat $n "col%d" ",">> <<typename $m 1>>) {
+	return <<repeat $n "m.Col(%d)" ",">>
+}
+
+<<if eq $m $n>>
+// Trace is a basic operation on a square matrix that simply
+// sums up all elements on the main diagonal (meaning all elements such that row==col).
+func (m <<$type>>) Trace() float32 {
+	return <<range $i := iter 0 $m>><<sep "+" $i>> m[<<mul $i $m | add $i>>]<<end>>
+}
+<<end>>
+
+// Abs returns the element-wise absolute value of this matrix
+func (m <<$type>>) Abs() <<$type>> {
+	return <<$type>>{<<repeat (mul $m $n) "Abs(m[%d])" ",">>}
+}
+
+// Pretty prints the matrix
+func (m <<$type>>) String() string {
+	buf := new(bytes.Buffer)
+	w := tabwriter.NewWriter(buf, 4, 4, 1, ' ', tabwriter.AlignRight)
+	for i := 0; i < <<$m>>; i++ {
+		for _, col := range m.Row(i) {
+			fmt.Fprintf(w, "%f\t", col)
+		}
+
+		fmt.Fprintln(w, "")
+	}
+	w.Flush()
+
+	return buf.String()
+}
+
+<<end>><<end>> <</* range $m range $n */>>

vendor/github.com/go-gl/mathgl/mgl32/mempool.go

@@ -0,0 +1,121 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import (
+	"sync"
+)
+
+var (
+	slicePools []*sync.Pool
+	listLock   sync.RWMutex
+)
+
+var shouldPool = true
+
+func DisableMemoryPooling() {
+	shouldPool = false
+}
+
+// Returns the given memory pool. If the pool doesn't exist, it will
+// create all pools up to element i. The number "i" corresponds to "p"
+// in most other comments. That is, it's Ceil(log_2(whatever)). So i=0
+// means you'll get the pool for slices of size 1, i=1 for size 2, i=2 for size 4,
+// and so on.
+//
+// This is concurrency safe and uses an RWMutex to protect the list expansion.
+func getPool(i int) *sync.Pool {
+	listLock.RLock()
+	if i >= len(slicePools) {
+
+		// Promote to a write lock because we now
+		// need to mutate the pool
+		listLock.RUnlock()
+		listLock.Lock()
+		defer listLock.Unlock()
+
+		for n := i - len(slicePools); n >= 0; n-- {
+			newFunc := genPoolNew(1 << uint(len(slicePools)))
+			slicePools = append(slicePools, &sync.Pool{New: newFunc})
+		}
+	} else {
+		defer listLock.RUnlock()
+	}
+
+	return slicePools[i]
+}
+
+func genPoolNew(i int) func() interface{} {
+	return func() interface{} {
+		return make([]float32, 0, i)
+	}
+}
+
+// Grabs a slice from the memory pool, such that its cap
+// is 2^p where p is Ceil(log_2(size)). It will be downsliced
+// such that the len is size.
+func grabFromPool(size int) []float32 {
+	pool, exact := binLog(size)
+
+	// Tried to grab something of size
+	// zero or less
+	if pool == -1 {
+		return nil
+	}
+
+	// If the log is not exact, we
+	// need to "overallocate" so we have
+	// log+1
+	if !exact {
+		pool++
+	}
+
+	slice := getPool(pool).Get().([]float32)
+	slice = slice[:size]
+	return slice
+}
+
+// Returns a slice to the appropriate pool. If the slice does not have a cap that's precisely
+// a power of 2, this will panic.
+func returnToPool(slice []float32) {
+	if cap(slice) == 0 {
+		return
+	}
+
+	pool, exact := binLog(cap(slice))
+
+	if !exact {
+		panic("attempt to pool slice with non-exact cap. If you're a user, please file an issue with github.com/go-gl/mathgl about this bug. This should never happen.")
+	}
+
+	getPool(pool).Put(slice)
+}
+
+// This returns the integer base 2 log of the value
+// and whether the log is exact or rounded down.
+//
+// This is only for positive integers.
+//
+// There are faster ways to do this, I'm open to suggestions. Most rely on knowing system endianness
+// which Go makes hard to do. I'm hesistant to use float conversions and the math package because of off-by-one errors.
+func binLog(val int) (int, bool) {
+	if val <= 0 {
+		return -1, false
+	}
+
+	exact := true
+	l := 0
+	for ; val > 1; val = val >> 1 {
+		// If the current lsb is 1 and the number
+		// is not equal to 1, this is not an exact
+		// log, but rather a rounding of it
+		if val&1 != 0 {
+			exact = false
+		}
+		l++
+	}
+
+	return l, exact
+}

vendor/github.com/go-gl/mathgl/mgl32/project.go

@@ -0,0 +1,98 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import (
+	"errors"
+	"math"
+)
+
+func Ortho(left, right, bottom, top, near, far float32) Mat4 {
+	rml, tmb, fmn := (right - left), (top - bottom), (far - near)
+
+	return Mat4{float32(2. / rml), 0, 0, 0, 0, float32(2. / tmb), 0, 0, 0, 0, float32(-2. / fmn), 0, float32(-(right + left) / rml), float32(-(top + bottom) / tmb), float32(-(far + near) / fmn), 1}
+}
+
+// Equivalent to Ortho with the near and far planes being -1 and 1, respectively
+func Ortho2D(left, right, bottom, top float32) Mat4 {
+	return Ortho(left, right, bottom, top, -1, 1)
+}
+
+func Perspective(fovy, aspect, near, far float32) Mat4 {
+	// fovy = (fovy * math.Pi) / 180.0 // convert from degrees to radians
+	nmf, f := near-far, float32(1./math.Tan(float64(fovy)/2.0))
+
+	return Mat4{float32(f / aspect), 0, 0, 0, 0, float32(f), 0, 0, 0, 0, float32((near + far) / nmf), -1, 0, 0, float32((2. * far * near) / nmf), 0}
+}
+
+func Frustum(left, right, bottom, top, near, far float32) Mat4 {
+	rml, tmb, fmn := (right - left), (top - bottom), (far - near)
+	A, B, C, D := (right+left)/rml, (top+bottom)/tmb, -(far+near)/fmn, -(2*far*near)/fmn
+
+	return Mat4{float32((2. * near) / rml), 0, 0, 0, 0, float32((2. * near) / tmb), 0, 0, float32(A), float32(B), float32(C), -1, 0, 0, float32(D), 0}
+}
+
+func LookAt(eyeX, eyeY, eyeZ, centerX, centerY, centerZ, upX, upY, upZ float32) Mat4 {
+	return LookAtV(Vec3{eyeX, eyeY, eyeZ}, Vec3{centerX, centerY, centerZ}, Vec3{upX, upY, upZ})
+}
+
+//  LookAtV generates a transform matrix from world space into the specific eye space
+func LookAtV(eye, center, up Vec3) Mat4 {
+	f := center.Sub(eye).Normalize()
+	s := f.Cross(up.Normalize()).Normalize()
+	u := s.Cross(f)
+
+	M := Mat4{
+		s[0], u[0], -f[0], 0,
+		s[1], u[1], -f[1], 0,
+		s[2], u[2], -f[2], 0,
+		0, 0, 0, 1,
+	}
+
+	return M.Mul4(Translate3D(float32(-eye[0]), float32(-eye[1]), float32(-eye[2])))
+}
+
+// Transform a set of coordinates from object space (in obj) to window coordinates (with depth)
+//
+// Window coordinates are continuous, not discrete (well, as continuous as an IEEE Floating Point can be), so you won't get exact pixel locations
+// without rounding or similar
+func Project(obj Vec3, modelview, projection Mat4, initialX, initialY, width, height int) (win Vec3) {
+	obj4 := obj.Vec4(1)
+
+	vpp := projection.Mul4(modelview).Mul4x1(obj4)
+	vpp = vpp.Mul(1 / vpp.W())
+	win[0] = float32(initialX) + (float32(width)*(vpp[0]+1))/2
+	win[1] = float32(initialY) + (float32(height)*(vpp[1]+1))/2
+	win[2] = (vpp[2] + 1) / 2
+
+	return win
+}
+
+// Transform a set of window coordinates to object space. If your MVP (projection.Mul(modelview) matrix is not invertible, this will return an error
+//
+// Note that the projection may not be perfect if you use strict pixel locations rather than the exact values given by Projectf.
+// (It's still unlikely to be perfect due to precision errors, but it will be closer)
+func UnProject(win Vec3, modelview, projection Mat4, initialX, initialY, width, height int) (obj Vec3, err error) {
+	inv := projection.Mul4(modelview).Inv()
+	var blank Mat4
+	if inv == blank {
+		return Vec3{}, errors.New("Could not find matrix inverse (projection times modelview is probably non-singular)")
+	}
+
+	obj4 := inv.Mul4x1(Vec4{
+		(2 * (win[0] - float32(initialX)) / float32(width)) - 1,
+		(2 * (win[1] - float32(initialY)) / float32(height)) - 1,
+		2*win[2] - 1,
+		1.0,
+	})
+	obj = obj4.Vec3()
+
+	//if obj4[3] > MinValue {}
+	obj[0] /= obj4[3]
+	obj[1] /= obj4[3]
+	obj[2] /= obj4[3]
+
+	return obj, nil
+}

vendor/github.com/go-gl/mathgl/mgl32/quat.go

@@ -0,0 +1,458 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import (
+	"math"
+)
+
+// A rotation order is the order in which
+// rotations will be transformed for the purposes of AnglesToQuat
+type RotationOrder int
+
+const (
+	XYX RotationOrder = iota
+	XYZ
+	XZX
+	XZY
+	YXY
+	YXZ
+	YZY
+	YZX
+	ZYZ
+	ZYX
+	ZXZ
+	ZXY
+)
+
+// A Quaternion is an extension of the imaginary numbers; there's all sorts of
+// interesting theory behind it. In 3D graphics we mostly use it as a cheap way of
+// representing rotation since quaternions are cheaper to multiply by, and easier to
+// interpolate than matrices.
+//
+// A Quaternion has two parts: W, the so-called scalar component,
+// and "V", the vector component. The vector component is considered to
+// be the part in 3D space, while W (loosely interpreted) is its 4D coordinate.
+type Quat struct {
+	W float32
+	V Vec3
+}
+
+// The quaternion identity: W=1; V=(0,0,0).
+//
+// As with all identities, multiplying any quaternion by this will yield the same
+// quaternion you started with.
+func QuatIdent() Quat {
+	return Quat{1., Vec3{0, 0, 0}}
+}
+
+// Creates an angle from an axis and an angle relative to that axis.
+//
+// This is cheaper than HomogRotate3D.
+func QuatRotate(angle float32, axis Vec3) Quat {
+	// angle = (float32(math.Pi) * angle) / 180.0
+
+	c, s := float32(math.Cos(float64(angle/2))), float32(math.Sin(float64(angle/2)))
+
+	return Quat{c, axis.Mul(s)}
+}
+
+// A convenient alias for q.V[0]
+func (q Quat) X() float32 {
+	return q.V[0]
+}
+
+// A convenient alias for q.V[1]
+func (q Quat) Y() float32 {
+	return q.V[1]
+}
+
+// A convenient alias for q.V[2]
+func (q Quat) Z() float32 {
+	return q.V[2]
+}
+
+// Adds two quaternions. It's no more complicated than
+// adding their W and V components.
+func (q1 Quat) Add(q2 Quat) Quat {
+	return Quat{q1.W + q2.W, q1.V.Add(q2.V)}
+}
+
+// Subtracts two quaternions. It's no more complicated than
+// subtracting their W and V components.
+func (q1 Quat) Sub(q2 Quat) Quat {
+	return Quat{q1.W - q2.W, q1.V.Sub(q2.V)}
+}
+
+// Multiplies two quaternions. This can be seen as a rotation. Note that
+// Multiplication is NOT commutative, meaning q1.Mul(q2) does not necessarily
+// equal q2.Mul(q1).
+func (q1 Quat) Mul(q2 Quat) Quat {
+	return Quat{q1.W*q2.W - q1.V.Dot(q2.V), q1.V.Cross(q2.V).Add(q2.V.Mul(q1.W)).Add(q1.V.Mul(q2.W))}
+}
+
+// Scales every element of the quaternion by some constant factor.
+func (q1 Quat) Scale(c float32) Quat {
+	return Quat{q1.W * c, Vec3{q1.V[0] * c, q1.V[1] * c, q1.V[2] * c}}
+}
+
+// Returns the conjugate of a quaternion. Equivalent to
+// Quat{q1.W, q1.V.Mul(-1)}
+func (q1 Quat) Conjugate() Quat {
+	return Quat{q1.W, q1.V.Mul(-1)}
+}
+
+// Returns the Length of the quaternion, also known as its Norm. This is the same thing as
+// the Len of a Vec4
+func (q1 Quat) Len() float32 {
+	return float32(math.Sqrt(float64(q1.W*q1.W + q1.V[0]*q1.V[0] + q1.V[1]*q1.V[1] + q1.V[2]*q1.V[2])))
+}
+
+// Norm() is an alias for Len() since both are very common terms.
+func (q1 Quat) Norm() float32 {
+	return q1.Len()
+}
+
+// Normalizes the quaternion, returning its versor (unit quaternion).
+//
+// This is the same as normalizing it as a Vec4.
+func (q1 Quat) Normalize() Quat {
+	length := q1.Len()
+
+	if FloatEqual(1, length) {
+		return q1
+	}
+	if length == 0 {
+		return QuatIdent()
+	}
+	if length == InfPos {
+		length = MaxValue
+	}
+
+	return Quat{q1.W * 1 / length, q1.V.Mul(1 / length)}
+}
+
+// The inverse of a quaternion. The inverse is equivalent
+// to the conjugate divided by the square of the length.
+//
+// This method computes the square norm by directly adding the sum
+// of the squares of all terms instead of actually squaring q1.Len(),
+// both for performance and precision.
+func (q1 Quat) Inverse() Quat {
+	return q1.Conjugate().Scale(1 / q1.Dot(q1))
+}
+
+// Rotates a vector by the rotation this quaternion represents.
+// This will result in a 3D vector. Strictly speaking, this is
+// equivalent to q1.v.q* where the "."" is quaternion multiplication and v is interpreted
+// as a quaternion with W 0 and V v. In code:
+// q1.Mul(Quat{0,v}).Mul(q1.Conjugate()), and
+// then retrieving the imaginary (vector) part.
+//
+// In practice, we hand-compute this in the general case and simplify
+// to save a few operations.
+func (q1 Quat) Rotate(v Vec3) Vec3 {
+	cross := q1.V.Cross(v)
+	// v + 2q_w * (q_v x v) + 2q_v x (q_v x v)
+	return v.Add(cross.Mul(2 * q1.W)).Add(q1.V.Mul(2).Cross(cross))
+}
+
+// Returns the homogeneous 3D rotation matrix corresponding to the quaternion.
+func (q1 Quat) Mat4() Mat4 {
+	w, x, y, z := q1.W, q1.V[0], q1.V[1], q1.V[2]
+	return Mat4{
+		1 - 2*y*y - 2*z*z, 2*x*y + 2*w*z, 2*x*z - 2*w*y, 0,
+		2*x*y - 2*w*z, 1 - 2*x*x - 2*z*z, 2*y*z + 2*w*x, 0,
+		2*x*z + 2*w*y, 2*y*z - 2*w*x, 1 - 2*x*x - 2*y*y, 0,
+		0, 0, 0, 1,
+	}
+}
+
+// The dot product between two quaternions, equivalent to if this was a Vec4
+func (q1 Quat) Dot(q2 Quat) float32 {
+	return q1.W*q2.W + q1.V[0]*q2.V[0] + q1.V[1]*q2.V[1] + q1.V[2]*q2.V[2]
+}
+
+// Returns whether the quaternions are approximately equal, as if
+// FloatEqual was called on each matching element
+func (q1 Quat) ApproxEqual(q2 Quat) bool {
+	return FloatEqual(q1.W, q2.W) && q1.V.ApproxEqual(q2.V)
+}
+
+// Returns whether the quaternions are approximately equal with a given tolerence, as if
+// FloatEqualThreshold was called on each matching element with the given epsilon
+func (q1 Quat) ApproxEqualThreshold(q2 Quat, epsilon float32) bool {
+	return FloatEqualThreshold(q1.W, q2.W, epsilon) && q1.V.ApproxEqualThreshold(q2.V, epsilon)
+}
+
+// Returns whether the quaternions are approximately equal using the given comparison function, as if
+// the function had been called on each individual element
+func (q1 Quat) ApproxEqualFunc(q2 Quat, f func(float32, float32) bool) bool {
+	return f(q1.W, q2.W) && q1.V.ApproxFuncEqual(q2.V, f)
+}
+
+// Returns whether the quaternions represents the same orientation
+//
+// Different values can represent the same orientation (q == -q) because quaternions avoid singularities
+// and discontinuities involved with rotation in 3 dimensions by adding extra dimensions
+func (q1 Quat) OrientationEqual(q2 Quat) bool {
+	return q1.OrientationEqualThreshold(q2, Epsilon)
+}
+
+// Returns whether the quaternions represents the same orientation with a given tolerence
+func (q1 Quat) OrientationEqualThreshold(q2 Quat, epsilon float32) bool {
+	return Abs(q1.Normalize().Dot(q2.Normalize())) > 1-epsilon
+}
+
+// Slerp is *S*pherical *L*inear Int*erp*olation, a method of interpolating
+// between two quaternions. This always takes the straightest path on the sphere between
+// the two quaternions, and maintains constant velocity.
+//
+// However, it's expensive and QuatSlerp(q1,q2) is not the same as QuatSlerp(q2,q1)
+func QuatSlerp(q1, q2 Quat, amount float32) Quat {
+	q1, q2 = q1.Normalize(), q2.Normalize()
+	dot := q1.Dot(q2)
+
+	// If the inputs are too close for comfort, linearly interpolate and normalize the result.
+	if dot > 0.9995 {
+		return QuatNlerp(q1, q2, amount)
+	}
+
+	// This is here for precision errors, I'm perfectly aware that *technically* the dot is bound [-1,1], but since Acos will freak out if it's not (even if it's just a liiiiitle bit over due to normal error) we need to clamp it
+	dot = Clamp(dot, -1, 1)
+
+	theta := float32(math.Acos(float64(dot))) * amount
+	c, s := float32(math.Cos(float64(theta))), float32(math.Sin(float64(theta)))
+	rel := q2.Sub(q1.Scale(dot)).Normalize()
+
+	return q1.Scale(c).Add(rel.Scale(s))
+}
+
+// *L*inear Int*erp*olation between two Quaternions, cheap and simple.
+//
+// Not excessively useful, but uses can be found.
+func QuatLerp(q1, q2 Quat, amount float32) Quat {
+	return q1.Add(q2.Sub(q1).Scale(amount))
+}
+
+// *Normalized* *L*inear Int*erp*olation between two Quaternions. Cheaper than Slerp
+// and usually just as good. This is literally Lerp with Normalize() called on it.
+//
+// Unlike Slerp, constant velocity isn't maintained, but it's much faster and
+// Nlerp(q1,q2) and Nlerp(q2,q1) return the same path. You should probably
+// use this more often unless you're suffering from choppiness due to the
+// non-constant velocity problem.
+func QuatNlerp(q1, q2 Quat, amount float32) Quat {
+	return QuatLerp(q1, q2, amount).Normalize()
+}
+
+// Performs a rotation in the specified order. If the order is not
+// a valid RotationOrder, this function will panic
+//
+// The rotation "order" is more of an axis descriptor. For instance XZX would
+// tell the function to interpret angle1 as a rotation about the X axis, angle2 about
+// the Z axis, and angle3 about the X axis again.
+//
+// Based off the code for the Matlab function "angle2quat", though this implementation
+// only supports 3 single angles as opposed to multiple angles.
+func AnglesToQuat(angle1, angle2, angle3 float32, order RotationOrder) Quat {
+	var s [3]float64
+	var c [3]float64
+
+	s[0], c[0] = math.Sincos(float64(angle1 / 2))
+	s[1], c[1] = math.Sincos(float64(angle2 / 2))
+	s[2], c[2] = math.Sincos(float64(angle3 / 2))
+
+	ret := Quat{}
+	switch order {
+	case ZYX:
+		ret.W = float32(c[0]*c[1]*c[2] + s[0]*s[1]*s[2])
+		ret.V = Vec3{float32(c[0]*c[1]*s[2] - s[0]*s[1]*c[2]),
+			float32(c[0]*s[1]*c[2] + s[0]*c[1]*s[2]),
+			float32(s[0]*c[1]*c[2] - c[0]*s[1]*s[2]),
+		}
+	case ZYZ:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float32(c[0]*s[1]*s[2] - s[0]*s[1]*c[2]),
+			float32(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float32(s[0]*c[1]*c[2] + c[0]*c[1]*s[2]),
+		}
+	case ZXY:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*s[1]*s[2])
+		ret.V = Vec3{float32(c[0]*s[1]*c[2] - s[0]*c[1]*s[2]),
+			float32(c[0]*c[1]*s[2] + s[0]*s[1]*c[2]),
+			float32(c[0]*s[1]*s[2] + s[0]*c[1]*c[2]),
+		}
+	case ZXZ:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float32(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float32(s[0]*s[1]*c[2] - c[0]*s[1]*s[2]),
+			float32(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+		}
+	case YXZ:
+		ret.W = float32(c[0]*c[1]*c[2] + s[0]*s[1]*s[2])
+		ret.V = Vec3{float32(c[0]*s[1]*c[2] + s[0]*c[1]*s[2]),
+			float32(s[0]*c[1]*c[2] - c[0]*s[1]*s[2]),
+			float32(c[0]*c[1]*s[2] - s[0]*s[1]*c[2]),
+		}
+	case YXY:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float32(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float32(s[0]*c[1]*c[2] + c[0]*c[1]*s[2]),
+			float32(c[0]*s[1]*s[2] - s[0]*s[1]*c[2]),
+		}
+	case YZX:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*s[1]*s[2])
+		ret.V = Vec3{float32(c[0]*c[1]*s[2] + s[0]*s[1]*c[2]),
+			float32(c[0]*s[1]*s[2] + s[0]*c[1]*c[2]),
+			float32(c[0]*s[1]*c[2] - s[0]*c[1]*s[2]),
+		}
+	case YZY:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float32(s[0]*s[1]*c[2] - c[0]*s[1]*s[2]),
+			float32(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+			float32(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+		}
+	case XYZ:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*s[1]*s[2])
+		ret.V = Vec3{float32(c[0]*s[1]*s[2] + s[0]*c[1]*c[2]),
+			float32(c[0]*s[1]*c[2] - s[0]*c[1]*s[2]),
+			float32(c[0]*c[1]*s[2] + s[0]*s[1]*c[2]),
+		}
+	case XYX:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float32(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+			float32(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+			float32(s[0]*s[1]*c[2] - c[0]*s[1]*s[2]),
+		}
+	case XZY:
+		ret.W = float32(c[0]*c[1]*c[2] + s[0]*s[1]*s[2])
+		ret.V = Vec3{float32(s[0]*c[1]*c[2] - c[0]*s[1]*s[2]),
+			float32(c[0]*c[1]*s[2] - s[0]*s[1]*c[2]),
+			float32(c[0]*s[1]*c[2] + s[0]*c[1]*s[2]),
+		}
+	case XZX:
+		ret.W = float32(c[0]*c[1]*c[2] - s[0]*c[1]*s[2])
+		ret.V = Vec3{float32(c[0]*c[1]*s[2] + s[0]*c[1]*c[2]),
+			float32(c[0]*s[1]*s[2] - s[0]*s[1]*c[2]),
+			float32(c[0]*s[1]*c[2] + s[0]*s[1]*s[2]),
+		}
+	default:
+		panic("Unsupported rotation order")
+	}
+	return ret
+}
+
+// Mat4ToQuat converts a pure rotation matrix into a quaternion
+func Mat4ToQuat(m Mat4) Quat {
+	// http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
+
+	if tr := m[0] + m[5] + m[10]; tr > 0 {
+		s := float32(0.5 / math.Sqrt(float64(tr+1.0)))
+		return Quat{
+			0.25 / s,
+			Vec3{
+				(m[6] - m[9]) * s,
+				(m[8] - m[2]) * s,
+				(m[1] - m[4]) * s,
+			},
+		}
+	}
+
+	if (m[0] > m[5]) && (m[0] > m[10]) {
+		s := float32(2.0 * math.Sqrt(float64(1.0+m[0]-m[5]-m[10])))
+		return Quat{
+			(m[6] - m[9]) / s,
+			Vec3{
+				0.25 * s,
+				(m[4] + m[1]) / s,
+				(m[8] + m[2]) / s,
+			},
+		}
+	}
+
+	if m[5] > m[10] {
+		s := float32(2.0 * math.Sqrt(float64(1.0+m[5]-m[0]-m[10])))
+		return Quat{
+			(m[8] - m[2]) / s,
+			Vec3{
+				(m[4] + m[1]) / s,
+				0.25 * s,
+				(m[9] + m[6]) / s,
+			},
+		}
+
+	}
+
+	s := float32(2.0 * math.Sqrt(float64(1.0+m[10]-m[0]-m[5])))
+	return Quat{
+		(m[1] - m[4]) / s,
+		Vec3{
+			(m[8] + m[2]) / s,
+			(m[9] + m[6]) / s,
+			0.25 * s,
+		},
+	}
+}
+
+// QuatLookAtV creates a rotation from an eye vector to a center vector
+//
+// It assumes the front of the rotated object at Z- and up at Y+
+func QuatLookAtV(eye, center, up Vec3) Quat {
+	// http://www.opengl-tutorial.org/intermediate-tutorials/tutorial-17-quaternions/#I_need_an_equivalent_of_gluLookAt__How_do_I_orient_an_object_towards_a_point__
+	// https://bitbucket.org/sinbad/ogre/src/d2ef494c4a2f5d6e2f0f17d3bfb9fd936d5423bb/OgreMain/src/OgreCamera.cpp?at=default#cl-161
+
+	direction := center.Sub(eye).Normalize()
+
+	// Find the rotation between the front of the object (that we assume towards Z-,
+	// but this depends on your model) and the desired direction
+	rotDir := QuatBetweenVectors(Vec3{0, 0, -1}, direction)
+
+	// Recompute up so that it's perpendicular to the direction
+	// You can skip that part if you really want to force up
+	//right := direction.Cross(up)
+	//up = right.Cross(direction)
+
+	// Because of the 1rst rotation, the up is probably completely screwed up.
+	// Find the rotation between the "up" of the rotated object, and the desired up
+	upCur := rotDir.Rotate(Vec3{0, 1, 0})
+	rotUp := QuatBetweenVectors(upCur, up)
+
+	rotTarget := rotUp.Mul(rotDir) // remember, in reverse order.
+	return rotTarget.Inverse()     // camera rotation should be inversed!
+}
+
+// QuatBetweenVectors calculates the rotation between two vectors
+func QuatBetweenVectors(start, dest Vec3) Quat {
+	// http://www.opengl-tutorial.org/intermediate-tutorials/tutorial-17-quaternions/#I_need_an_equivalent_of_gluLookAt__How_do_I_orient_an_object_towards_a_point__
+	// https://github.com/g-truc/glm/blob/0.9.5/glm/gtx/quaternion.inl#L225
+	// https://bitbucket.org/sinbad/ogre/src/d2ef494c4a2f5d6e2f0f17d3bfb9fd936d5423bb/OgreMain/include/OgreVector3.h?at=default#cl-654
+
+	start = start.Normalize()
+	dest = dest.Normalize()
+	epsilon := float32(0.001)
+
+	cosTheta := start.Dot(dest)
+	if cosTheta < -1.0+epsilon {
+		// special case when vectors in opposite directions:
+		// there is no "ideal" rotation axis
+		// So guess one; any will do as long as it's perpendicular to start
+		axis := Vec3{1, 0, 0}.Cross(start)
+		if axis.Dot(axis) < epsilon {
+			// bad luck, they were parallel, try again!
+			axis = Vec3{0, 1, 0}.Cross(start)
+		}
+
+		return QuatRotate(math.Pi, axis.Normalize())
+	}
+
+	axis := start.Cross(dest)
+	s := float32(math.Sqrt(float64(1.0+cosTheta) * 2.0))
+
+	return Quat{
+		s * 0.5,
+		axis.Mul(1.0 / s),
+	}
+}

vendor/github.com/go-gl/mathgl/mgl32/shapes.go

@@ -0,0 +1,306 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import (
+	"fmt"
+	"math"
+)
+
+// Generates a circle centered at (0,0) with a given radius.
+// The radii are assumed to be in GL's coordinate sizing.
+//
+// Technically this draws an ellipse with two axes that match with the X and Y axes, the reason it has a radiusX and radiusY is because GL's coordinate system
+// is proportional to screen width and screen height. So if you have a non-square viewport, a single radius will appear
+// to "squash" the circle in one direction (usually the Y direction), so the X and Y radius allow for a circle to be made.
+// A good way to get the correct radii is with mathgl.ScreenToGLCoords(radius, radius, screenWidth, screenHeight) which will get you the correct
+// proportional GL coords.
+//
+// The numSlices argument specifies how many triangles you want your circle divided into, setting this
+// number to too low a value may cause problem (and too high will cause it to take a lot of memory and time to compute
+// without much gain in resolution).
+//
+// This uses discrete triangles, not a triangle fan
+func Circle(radiusX, radiusY float32, numSlices int) []Vec2 {
+	twoPi := float32(2.0 * math.Pi)
+
+	circlePoints := make([]Vec2, 0, numSlices*3)
+	center := Vec2{0.0, 0.0}
+	previous := Vec2{radiusX, 0.0}
+
+	for theta := twoPi / float32(numSlices); !FloatEqual(theta, twoPi); theta = Clamp(theta+twoPi/float32(numSlices), 0.0, twoPi) {
+		sin, cos := math.Sincos(float64(theta))
+		curr := Vec2{float32(cos) * radiusX, float32(sin) * radiusY}
+
+		circlePoints = append(circlePoints, center, previous, curr)
+		previous = curr
+	}
+
+	// Now add the final point at theta=2pi
+	circlePoints = append(circlePoints, center, previous, Vec2{radiusX, 0.0})
+	return circlePoints
+}
+
+// Generates a 2-triangle rectangle for use with GL_TRIANGLES. The width and height should use GL's proportions (that is, where a width of 1.0
+// is equivalent to half of the width of the render target); however, the y-coordinates grow downwards, not upwards. That is, it
+// assumes you want the origin of the rectangle with the top-left corner at (0.0,0.0).
+//
+// Keep in mind that GL's coordinate system is proportional, so width=height will not result in a square unless your viewport is square.
+// If you want to maintain proportionality regardless of screen size, use the results of w,h := ScreenToGLCoordsf(absoluteWidth, absoluteHeight, screenWidth, screenHeight);
+// w,h=w+1,h-1 in the call to this function. (The w+1,h-1 step maps the coordinates to start at 0.0 rather than -1.0)
+func Rect(width, height float32) []Vec2 {
+	return []Vec2{
+		{0.0, 0.0},
+		{0.0, -height},
+		{width, -height},
+
+		{0.0, 0.0},
+		{width, -height},
+		{width, 0.0},
+	}
+}
+
+func QuadraticBezierCurve2D(t float32, cPoint1, cPoint2, cPoint3 Vec2) Vec2 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1.0 - t) * (1.0 - t)).Add(cPoint2.Mul(2 * (1 - t) * t)).Add(cPoint3.Mul(t * t))
+}
+
+func QuadraticBezierCurve3D(t float32, cPoint1, cPoint2, cPoint3 Vec3) Vec3 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1.0 - t) * (1.0 - t)).Add(cPoint2.Mul(2 * (1 - t) * t)).Add(cPoint3.Mul(t * t))
+}
+
+func CubicBezierCurve2D(t float32, cPoint1, cPoint2, cPoint3, cPoint4 Vec2) Vec2 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1 - t) * (1 - t) * (1 - t)).Add(cPoint2.Mul(3 * (1 - t) * (1 - t) * t)).Add(cPoint3.Mul(3 * (1 - t) * t * t)).Add(cPoint4.Mul(t * t * t))
+}
+
+func CubicBezierCurve3D(t float32, cPoint1, cPoint2, cPoint3, cPoint4 Vec3) Vec3 {
+	if t < 0.0 || t > 1.0 {
+		panic("Can't interpolate on bezier curve with t out of range [0.0,1.0]")
+	}
+
+	return cPoint1.Mul((1 - t) * (1 - t) * (1 - t)).Add(cPoint2.Mul(3 * (1 - t) * (1 - t) * t)).Add(cPoint3.Mul(3 * (1 - t) * t * t)).Add(cPoint4.Mul(t * t * t))
+}
+
+// Returns the point at point t along an n-control point Bezier curve
+//
+// t must be in the range 0.0 and 1.0 or this function will panic. Consider [0.0,1.0] to be similar to a percentage,
+// 0.0 is first control point, and the point at 1.0 is the last control point. Any point in between is how far along the path you are between 0 and 1.
+//
+// This function is not sensitive to the coordinate system of the control points. It will correctly interpolate regardless of whether they're in screen coords,
+// gl coords, or something else entirely
+func BezierCurve2D(t float32, cPoints []Vec2) Vec2 {
+	if t < 0.0 || t > 1.0 {
+		panic("Input to bezier has t not in range [0,1]. If you think this is a precision error, use mathgl.Clamp[f|d] before calling this function")
+	}
+
+	n := len(cPoints) - 1
+	point := cPoints[0].Mul(float32(math.Pow(float64(1.0-t), float64(n))))
+
+	for i := 1; i <= n; i++ {
+		point = point.Add(cPoints[i].Mul(float32(float64(choose(n, i)) * math.Pow(float64(1-t), float64(n-i)) * math.Pow(float64(t), float64(i))))) // P += P_i * nCi * (1-t)^(n-i) * t^i
+	}
+
+	return point
+}
+
+// Same as the 2D version, except the line is in 3D space
+func BezierCurve3D(t float32, cPoints []Vec3) Vec3 {
+	if t < 0.0 || t > 1.0 {
+		panic("Input to bezier has t not in range [0,1]. If you think this is a precision error, use mathgl.Clamp[f|d] before calling this function")
+	}
+
+	n := len(cPoints) - 1
+	point := cPoints[0].Mul(float32(math.Pow(float64(1.0-t), float64(n))))
+
+	for i := 1; i <= n; i++ {
+		point = point.Add(cPoints[i].Mul(float32(float64(choose(n, i)) * math.Pow(float64(1-t), float64(n-i)) * math.Pow(float64(t), float64(i))))) // P += P_i * nCi * (1-t)^(n-i) * t^i
+	}
+
+	return point
+}
+
+// Generates a bezier curve with controlPoints cPoints. The numPoints argument
+// determines how many "samples" it makes along the line. For instance, a
+// call to this with numPoints 2 will have exactly two points: the start and end points
+// For any points above 2 it will divide it into numPoints-1 chunks (which means it will generate numPoints-2 vertices other than the beginning and end).
+// So for 3 points it will divide it in half, 4 points into thirds, and so on.
+//
+// This is likely to get rather expensive for anything over perhaps a cubic curve.
+func MakeBezierCurve2D(numPoints int, cPoints []Vec2) (line []Vec2) {
+	line = make([]Vec2, numPoints)
+	if numPoints == 0 {
+		return
+	} else if numPoints == 1 {
+		line[0] = cPoints[0]
+		return
+	} else if numPoints == 2 {
+		line[0] = cPoints[0]
+		line[1] = cPoints[len(cPoints)-1]
+		return
+	}
+
+	line[0] = cPoints[0]
+	for i := 1; i < numPoints-1; i++ {
+		line[i] = BezierCurve2D(Clamp(float32(i)/float32(numPoints-1), 0.0, 1.0), cPoints)
+	}
+	line[numPoints-1] = cPoints[len(cPoints)-1]
+
+	return
+}
+
+// Same as the 2D version, except with the line in 3D space
+func MakeBezierCurve3D(numPoints int, cPoints []Vec3) (line []Vec3) {
+	line = make([]Vec3, numPoints)
+	if numPoints == 0 {
+		return
+	} else if numPoints == 1 {
+		line[0] = cPoints[0]
+		return
+	} else if numPoints == 2 {
+		line[0] = cPoints[0]
+		line[1] = cPoints[len(cPoints)-1]
+		return
+	}
+
+	line[0] = cPoints[0]
+	for i := 1; i < numPoints-1; i++ {
+		line[i] = BezierCurve3D(Clamp(float32(i)/float32(numPoints-1), 0.0, 1.0), cPoints)
+	}
+	line[numPoints-1] = cPoints[len(cPoints)-1]
+
+	return
+}
+
+// Creates a 2-dimensional Bezier surface of arbitrary degree in 3D Space
+// Like the curve functions, if u or v are not in the range [0.0,1.0] the function will panic, use Clamp[f|d]
+// to ensure it is correct.
+//
+// The control point matrix must not be jagged, or this will end up panicking from an index out of bounds exception
+func BezierSurface(u, v float32, cPoints [][]Vec3) Vec3 {
+	if u < 0.0 || u > 1.0 || v < 1.0 || v > 1.0 {
+		panic("u or v not in range [0.0,1.0] in BezierSurface")
+	}
+
+	n := len(cPoints) - 1
+	m := len(cPoints[0]) - 1
+
+	point := cPoints[0][0].Mul(float32(math.Pow(float64(1.0-u), float64(n)) * math.Pow(float64(1.0-v), float64(m))))
+
+	for i := 0; i <= n; i++ {
+		for j := 0; j <= m; j++ {
+			if i == 0 && j == 0 {
+				continue
+			}
+
+			point = point.Add(cPoints[i][j].Mul(float32(float64(choose(n, i)) * math.Pow(float64(u), float64(i)) * math.Pow(float64(1.0-u), float64(n-i)) * float64(choose(m, j)) * math.Pow(float64(v), float64(j)) * math.Pow(float64(1.0-v), float64(m-j)))))
+		}
+	}
+
+	return point
+}
+
+// Does interpolation over a spline of several bezier curves. Each bezier curve must have a finite range,
+// though the spline may be disjoint. The bezier curves are not required to be in any particular order.
+//
+// If t is out of the range of all given curves, this function will panic
+func BezierSplineInterpolate2D(t float32, ranges [][2]float32, cPoints [][]Vec2) Vec2 {
+	if len(ranges) != len(cPoints) {
+		panic("Each bezier curve needs a range")
+	}
+
+	for i, curveRange := range ranges {
+		if t >= curveRange[0] && t <= curveRange[1] {
+			return BezierCurve2D((t-curveRange[0])/(curveRange[1]-curveRange[0]), cPoints[i])
+		}
+	}
+
+	panic("t is out of the range of all bezier curves in this spline")
+}
+
+// Does interpolation over a spline of several bezier curves. Each bezier curve must have a finite range,
+// though the spline may be disjoint. The bezier curves are not required to be in any particular order.
+//
+// If t is out of the range of all given curves, this function will panic
+func BezierSplineInterpolate3D(t float32, ranges [][2]float32, cPoints [][]Vec3) Vec3 {
+	if len(ranges) != len(cPoints) {
+		panic("Each bezier curve needs a range")
+	}
+
+	for i, curveRange := range ranges {
+		if t >= curveRange[0] && t <= curveRange[1] {
+			return BezierCurve3D((t-curveRange[0])/(curveRange[1]-curveRange[0]), cPoints[i])
+		}
+	}
+
+	panic("t is out of the range of all bezier curves in this spline")
+}
+
+// Reticulates ALL the Splines
+//
+// For the overly serious: the function is just for fun. It does nothing except prints a Maxis reference. Technically you could "reticulate splines"
+// by joining a bunch of splines together, but that ruins the joke.
+func ReticulateSplines(ranges [][][2]float32, cPoints [][][]Vec2, withLlamas bool) {
+	if !withLlamas {
+		fmt.Println("You can't reticulate splines without llamas, silly.")
+	} else {
+		fmt.Println("Actually, you can't even reticulate splines WITH llamas")
+	}
+}
+
+// Transform from pixel coordinates to GL coordinates.
+//
+// This assumes that your pixel coordinate system considers its origin to be in the top left corner (GL's is in the bottom left).
+// The coordinates x and y may be out of the range [0,screenWidth-1] and [0,screeneHeight-1].
+//
+// GL's coordinate system maps [screenWidth-1,0] to [1.0,1.0] and [0,screenHeight-1] to [-1.0,-1.0]. If x and y are out of the range, they'll still
+// be mapped correctly, just off the screen. (e.g. if y = 2*(screenHeight-1) you'll get -3.0 for yOut)
+//
+// This is similar to Unproject, except for 2D cases and much simpler (especially since an inverse may always be found)
+func ScreenToGLCoords(x, y int, screenWidth, screenHeight int) (xOut, yOut float32) {
+	xOut = 2.0*float32(x)/float32(screenWidth-1) - 1.0
+	yOut = -2.0*float32(y)/float32(screenHeight-1) + 1.0
+
+	return
+}
+
+// Transform from GL's proportional system to pixel coordinates.
+//
+// Assumes the pixel coordinate system has its origin in the top left corner. (GL's is in the bottom left)
+//
+// GL's coordinate system maps [screenWidth-1,0] to [1.0,1.0] and [0,screenHeight-1] to [-1.0,-1.0]. If x and y are out of the range, they'll still
+// be mapped correctly, just off the screen. (e.g. if y=-3.0, you'll get 2*(screenHeight-1) for yOut)
+//
+// This is similar to Project, except for 2D cases and much simpler
+func GLToScreenCoords(x, y float32, screenWidth, screenHeight int) (xOut, yOut int) {
+	xOut = int((x + 1.0) * float32(screenWidth-1) / 2.0)
+	yOut = int((1.0 - y) * float32(screenHeight-1) / 2.0)
+
+	return
+}
+
+func choose(n, k int) (result int) {
+	if k == 0 {
+		return 1
+	} else if n == 0 {
+		return 0
+	}
+	result = (n - (k - 1))
+	for i := 2; i <= k; i++ {
+		result *= (n - (k - i)) / i
+	}
+
+	return result
+}

vendor/github.com/go-gl/mathgl/mgl32/transform.go

@@ -0,0 +1,223 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import "math"
+
+// Rotate2D returns a rotation Matrix about a angle in 2-D space. Specifically about the origin.
+// It is a 2x2 matrix, if you need a 3x3 for Homogeneous math (e.g. composition with a Translation matrix)
+// see HomogRotate2D
+func Rotate2D(angle float32) Mat2 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+	return Mat2{cos, sin, -sin, cos}
+}
+
+// Rotate3DX returns a 3x3 (non-homogeneous) Matrix that rotates by angle about the X-axis
+//
+// Where c is cos(angle) and s is sin(angle)
+//    [1 0 0]
+//    [0 c -s]
+//    [0 s c ]
+func Rotate3DX(angle float32) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+
+	return Mat3{1, 0, 0, 0, cos, sin, 0, -sin, cos}
+}
+
+// Rotate3DY returns a 3x3 (non-homogeneous) Matrix that rotates by angle about the Y-axis
+//
+// Where c is cos(angle) and s is sin(angle)
+//    [c 0 s]
+//    [0 1 0]
+//    [s 0 c ]
+func Rotate3DY(angle float32) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+
+	return Mat3{cos, 0, -sin, 0, 1, 0, sin, 0, cos}
+}
+
+// Rotate3DZ returns a 3x3 (non-homogeneous) Matrix that rotates by angle about the Z-axis
+//
+// Where c is cos(angle) and s is sin(angle)
+//    [c -s 0]
+//    [s c 0]
+//    [0 0 1 ]
+func Rotate3DZ(angle float32) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+
+	return Mat3{cos, sin, 0, -sin, cos, 0, 0, 0, 1}
+}
+
+// Translate2D returns a homogeneous (3x3 for 2D-space) Translation matrix that moves a point by Tx units in the x-direction and Ty units in the y-direction
+//
+//    [[1, 0, Tx]]
+//    [[0, 1, Ty]]
+//    [[0, 0, 1 ]]
+func Translate2D(Tx, Ty float32) Mat3 {
+	return Mat3{1, 0, 0, 0, 1, 0, float32(Tx), float32(Ty), 1}
+}
+
+// Translate3D returns a homogeneous (4x4 for 3D-space) Translation matrix that moves a point by Tx units in the x-direction, Ty units in the y-direction,
+// and Tz units in the z-direction
+//
+//    [[1, 0, 0, Tx]]
+//    [[0, 1, 0, Ty]]
+//    [[0, 0, 1, Tz]]
+//    [[0, 0, 0, 1 ]]
+func Translate3D(Tx, Ty, Tz float32) Mat4 {
+	return Mat4{1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, float32(Tx), float32(Ty), float32(Tz), 1}
+}
+
+// Same as Rotate2D, except homogeneous (3x3 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate2D(angle float32) Mat3 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+	return Mat3{cos, sin, 0, -sin, cos, 0, 0, 0, 1}
+}
+
+// Same as Rotate3DX, except homogeneous (4x4 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate3DX(angle float32) Mat4 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+
+	return Mat4{1, 0, 0, 0, 0, cos, sin, 0, 0, -sin, cos, 0, 0, 0, 0, 1}
+}
+
+// Same as Rotate3DY, except homogeneous (4x4 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate3DY(angle float32) Mat4 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+	return Mat4{cos, 0, -sin, 0, 0, 1, 0, 0, sin, 0, cos, 0, 0, 0, 0, 1}
+}
+
+// Same as Rotate3DZ, except homogeneous (4x4 with the extra row/col being all zeroes with a one in the bottom right)
+func HomogRotate3DZ(angle float32) Mat4 {
+	//angle = (angle * math.Pi) / 180.0
+	sin, cos := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+	return Mat4{cos, sin, 0, 0, -sin, cos, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}
+}
+
+// Scale3D creates a homogeneous 3D scaling matrix
+// [[ scaleX, 0     , 0     , 0 ]]
+// [[ 0     , scaleY, 0     , 0 ]]
+// [[ 0     , 0     , scaleZ, 0 ]]
+// [[ 0     , 0     , 0     , 1 ]]
+func Scale3D(scaleX, scaleY, scaleZ float32) Mat4 {
+
+	return Mat4{float32(scaleX), 0, 0, 0, 0, float32(scaleY), 0, 0, 0, 0, float32(scaleZ), 0, 0, 0, 0, 1}
+}
+
+// Scale2D creates a homogeneous 2D scaling matrix
+// [[ scaleX, 0     , 0 ]]
+// [[ 0     , scaleY, 0 ]]
+// [[ 0     , 0     , 1 ]]
+func Scale2D(scaleX, scaleY float32) Mat3 {
+	return Mat3{float32(scaleX), 0, 0, 0, float32(scaleY), 0, 0, 0, 1}
+}
+
+// ShearX2D creates a homogeneous 2D shear matrix along the X-axis
+func ShearX2D(shear float32) Mat3 {
+	return Mat3{1, 0, 0, float32(shear), 1, 0, 0, 0, 1}
+}
+
+// ShearY2D creates a homogeneous 2D shear matrix along the Y-axis
+func ShearY2D(shear float32) Mat3 {
+	return Mat3{1, float32(shear), 0, 0, 1, 0, 0, 0, 1}
+}
+
+// ShearX3D creates a homogeneous 3D shear matrix along the X-axis
+func ShearX3D(shearY, shearZ float32) Mat4 {
+
+	return Mat4{1, float32(shearY), float32(shearZ), 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1}
+}
+
+// ShearY3D creates a homogeneous 3D shear matrix along the Y-axis
+func ShearY3D(shearX, shearZ float32) Mat4 {
+	return Mat4{1, 0, 0, 0, float32(shearX), 1, float32(shearZ), 0, 0, 0, 1, 0, 0, 0, 0, 1}
+}
+
+// ShearZ3D creates a homogeneous 3D shear matrix along the Z-axis
+func ShearZ3D(shearX, shearY float32) Mat4 {
+	return Mat4{1, 0, 0, 0, 0, 1, 0, 0, float32(shearX), float32(shearY), 1, 0, 0, 0, 0, 1}
+}
+
+// HomogRotate3D creates a 3D rotation Matrix that rotates by (radian) angle about some arbitrary axis given by a Vector.
+// It produces a homogeneous matrix (4x4)
+//
+// Where c is cos(angle) and s is sin(angle), and x, y, and z are the first, second, and third elements of the axis vector (respectively):
+//
+//    [[ x^2(1-c)+c, xy(1-c)-zs, xz(1-c)+ys, 0 ]]
+//    [[ xy(1-c)+zs, y^2(1-c)+c, yz(1-c)-xs, 0 ]]
+//    [[ xz(1-c)-ys, yz(1-c)+xs, z^2(1-c)+c, 0 ]]
+//    [[ 0         , 0         , 0         , 1 ]]
+func HomogRotate3D(angle float32, axis Vec3) Mat4 {
+	x, y, z := axis[0], axis[1], axis[2]
+	s, c := float32(math.Sin(float64(angle))), float32(math.Cos(float64(angle)))
+	k := 1 - c
+
+	return Mat4{x*x*k + c, x*y*k + z*s, x*z*k - y*s, 0, x*y*k - z*s, y*y*k + c, y*z*k + x*s, 0, x*z*k + y*s, y*z*k - x*s, z*z*k + c, 0, 0, 0, 0, 1}
+}
+
+// Extracts the 3d scaling from a homogeneous matrix
+func Extract3DScale(m Mat4) (x, y, z float32) {
+	return float32(math.Sqrt(float64(m[0]*m[0] + m[1]*m[1] + m[2]*m[2]))),
+		float32(math.Sqrt(float64(m[4]*m[4] + m[5]*m[5] + m[6]*m[6]))),
+		float32(math.Sqrt(float64(m[8]*m[8] + m[9]*m[9] + m[10]*m[10])))
+}
+
+// Extracts the maximum scaling from a homogeneous matrix
+func ExtractMaxScale(m Mat4) float32 {
+	scaleX := float64(m[0]*m[0] + m[1]*m[1] + m[2]*m[2])
+	scaleY := float64(m[4]*m[4] + m[5]*m[5] + m[6]*m[6])
+	scaleZ := float64(m[8]*m[8] + m[9]*m[9] + m[10]*m[10])
+
+	return float32(math.Sqrt(math.Max(scaleX, math.Max(scaleY, scaleZ))))
+}
+
+// Calculates the Normal of the Matrix (aka the inverse transpose)
+func Mat4Normal(m Mat4) Mat3 {
+	n := m.Inv().Transpose()
+	return Mat3{n[0], n[1], n[2], n[4], n[5], n[6], n[8], n[9], n[10]}
+}
+
+// Multiplies a 3D vector by a transformation given by
+// the homogeneous 4D matrix m, applying any translation.
+// If this transformation is non-affine, it will project this
+// vector onto the plane w=1 before returning the result.
+//
+// This is similar to saying you're transforming and projecting a point.
+//
+// This is effectively equivalent to the GLSL code
+//     vec4 r = (m * vec4(v,1.));
+//     r = r/r.w;
+//     vec3 newV = r.xyz;
+func TransformCoordinate(v Vec3, m Mat4) Vec3 {
+	t := v.Vec4(1)
+	t = m.Mul4x1(t)
+	t = t.Mul(1 / t[3])
+
+	return t.Vec3()
+}
+
+// Multiplies a 3D vector by a transformation given by
+// the homogeneous 4D matrix m, NOT applying any translations.
+//
+// This is similar to saying you're applying a transformation
+// to a direction or normal. Rotation still applies (as does scaling),
+// but translating a direction or normal is meaningless.
+//
+// This is effectively equivalent to the GLSL code
+//    vec4 r = (m * vec4(v,0.));
+//    vec3 newV = r.xyz
+func TransformNormal(v Vec3, m Mat4) Vec3 {
+	t := v.Vec4(0)
+	t = m.Mul4x1(t)
+
+	return t.Vec3()
+}

vendor/github.com/go-gl/mathgl/mgl32/util.go

@@ -0,0 +1,142 @@
+// Copyright 2014 The go-gl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+//go:generate go run codegen.go -template vector.tmpl -output vector.go
+//go:generate go run codegen.go -template matrix.tmpl -output matrix.go
+//go:generate go run codegen.go -mgl64
+
+package mgl32
+
+import (
+	"math"
+)
+
+// Epsilon is some tiny value that determines how precisely equal we want our floats to be
+// This is exported and left as a variable in case you want to change the default threshold for the
+// purposes of certain methods (e.g. Unproject uses the default epsilon when determining
+// if the determinant is "close enough" to zero to mean there's no inverse).
+//
+// This is, obviously, not mutex protected so be **absolutely sure** that no functions using Epsilon
+// are being executed when you change this.
+var Epsilon float32 = 1e-10
+
+// A direct copy of the math package's Abs. This is here for the mgl32
+// package, to prevent rampant type conversions during equality tests.
+func Abs(a float32) float32 {
+	if a < 0 {
+		return -a
+	} else if a == 0 {
+		return 0
+	}
+
+	return a
+}
+
+// FloatEqual is a safe utility function to compare floats.
+// It's Taken from http://floating-point-gui.de/errors/comparison/
+//
+// It is slightly altered to not call Abs when not needed.
+func FloatEqual(a, b float32) bool {
+	return FloatEqualThreshold(a, b, Epsilon)
+}
+
+// FloatEqualFunc is a utility closure that will generate a function that
+// always approximately compares floats like FloatEqualThreshold with a different
+// threshold.
+func FloatEqualFunc(epsilon float32) func(float32, float32) bool {
+	return func(a, b float32) bool {
+		return FloatEqualThreshold(a, b, epsilon)
+	}
+}
+
+var (
+	MinNormal = float32(1.1754943508222875e-38) // 1 / 2**(127 - 1)
+	MinValue  = float32(math.SmallestNonzeroFloat32)
+	MaxValue  = float32(math.MaxFloat32)
+
+	InfPos = float32(math.Inf(1))
+	InfNeg = float32(math.Inf(-1))
+	NaN    = float32(math.NaN())
+)
+
+// FloatEqualThreshold is a utility function to compare floats.
+// It's Taken from http://floating-point-gui.de/errors/comparison/
+//
+// It is slightly altered to not call Abs when not needed.
+//
+// This differs from FloatEqual in that it lets you pass in your comparison threshold, so that you can adjust the comparison value to your specific needs
+func FloatEqualThreshold(a, b, epsilon float32) bool {
+	if a == b { // Handles the case of inf or shortcuts the loop when no significant error has accumulated
+		return true
+	}
+
+	diff := Abs(a - b)
+	if a*b == 0 || diff < MinNormal { // If a or b are 0 or both are extremely close to it
+		return diff < epsilon*epsilon
+	}
+
+	// Else compare difference
+	return diff/(Abs(a)+Abs(b)) < epsilon
+}
+
+// Clamp takes in a value and two thresholds. If the value is smaller than the low
+// threshold, it returns the low threshold. If it's bigger than the high threshold
+// it returns the high threshold. Otherwise it returns the value.
+//
+// Useful to prevent some functions from freaking out because a value was
+// teeeeechnically out of range.
+func Clamp(a, low, high float32) float32 {
+	if a < low {
+		return low
+	} else if a > high {
+		return high
+	}
+
+	return a
+}
+
+// ClampFunc generates a closure that returns its parameter
+// clamped to the range [low,high].
+func ClampFunc(low, high float32) func(float32) float32 {
+	return func(a float32) float32 {
+		return Clamp(a, low, high)
+	}
+}
+
+/* The IsClamped functions use strict equality (meaning: not the FloatEqual function)
+there shouldn't be any major issues with this since clamp is often used to fix minor errors*/
+
+// Checks if a is clamped between low and high as if
+// Clamp(a, low, high) had been called.
+//
+// In most cases it's probably better to just call Clamp
+// without checking this since it's relatively cheap.
+func IsClamped(a, low, high float32) bool {
+	return a >= low && a <= high
+}
+
+// If a > b, then a will be set to the value of b.
+func SetMin(a, b *float32) {
+	if *b < *a {
+		*a = *b
+	}
+}
+
+// If a < b, then a will be set to the value of b.
+func SetMax(a, b *float32) {
+	if *a < *b {
+		*a = *b
+	}
+}
+
+// Round shortens a float32 value to a specified precision (number of digits after the decimal point)
+// with "round half up" tie-braking rule. Half-way values (23.5) are always rounded up (24).
+func Round(v float32, precision int) float32 {
+	p := float64(precision)
+	t := float64(v) * math.Pow(10, p)
+	if t > 0 {
+		return float32(math.Floor(t+0.5) / math.Pow(10, p))
+	}
+	return float32(math.Ceil(t-0.5) / math.Pow(10, p))
+}

vendor/github.com/go-gl/mathgl/mgl32/vecn.go

@@ -0,0 +1,350 @@
+// Copyright 2014 The go-gl 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 mgl32
+
+import (
+	"math"
+)
+
+// A vector of N elements backed by a slice
+//
+// As with MatMxN, this is not for hardcore linear algebra with large dimensions. Use github.com/gonum/matrix
+// or something like BLAS/LAPACK for that. This is for corner cases in 3D math where you require
+// something a little bigger that 4D, but still relatively small.
+//
+// This VecN uses several sync.Pool objects as a memory pool. The rule is that for any sized vector, the backing slice
+// has CAPACITY (not length) of 2^p where p is Ceil(log_2(N)) -- or in other words, rounding up the base-2
+// log of the size of the vector. E.G. a VecN of size 17 will have a backing slice of Cap 32.
+type VecN struct {
+	vec []float32
+}
+
+// Creates a new vector with a backing slice filled with the contents
+// of initial. It is NOT backed by initial, but rather a slice with cap
+// 2^p where p is Ceil(log_2(len(initial))), with the data from initial copied into
+// it.
+func NewVecNFromData(initial []float32) *VecN {
+	if initial == nil {
+		return &VecN{}
+	}
+	var internal []float32
+	if shouldPool {
+		internal = grabFromPool(len(initial))
+	} else {
+		internal = make([]float32, len(initial))
+	}
+	copy(internal, initial)
+	return &VecN{vec: internal}
+}
+
+// Creates a new vector with a backing slice of
+// 2^p where p = Ceil(log_2(n))
+func NewVecN(n int) *VecN {
+	if shouldPool {
+		return &VecN{vec: grabFromPool(n)}
+	} else {
+		return &VecN{vec: make([]float32, n)}
+	}
+}
+
+// Returns the raw slice backing the VecN
+//
+// This may be sent back to the memory pool at any time
+// and you aren't advised to rely on this value
+func (vn VecN) Raw() []float32 {
+	return vn.vec
+}
+
+// Gets the element at index i from the vector.
+// This does not bounds check, and will panic if i is
+// out of range.
+func (vn VecN) Get(i int) float32 {
+	return vn.vec[i]
+}
+
+func (vn *VecN) Set(i int, val float32) {
+	vn.vec[i] = val
+}
+
+// Sends the allocated memory through the callback if it exists
+func (vn *VecN) destroy() {
+	if vn == nil || vn.vec == nil {
+		return
+	}
+
+	if shouldPool {
+		returnToPool(vn.vec)
+	}
+	vn.vec = nil
+}
+
+// Resizes the underlying slice to the desired amount, reallocating or retrieving from the pool
+// if necessary. The values after a Resize cannot be expected to be related to the values before a Resize.
+//
+// If the caller is a nil pointer, this returns a value as if NewVecN(n) had been called,
+// otherwise it simply returns the caller.
+func (vn *VecN) Resize(n int) *VecN {
+	if vn == nil {
+		return NewVecN(n)
+	}
+
+	if n <= cap(vn.vec) {
+		if vn.vec != nil {
+			vn.vec = vn.vec[:n]
+		} else {
+			vn.vec = []float32{}
+		}
+		return vn
+	}
+
+	if shouldPool && vn.vec != nil {
+		returnToPool(vn.vec)
+	}
+	*vn = (*NewVecN(n))
+
+	return vn
+}
+
+// Sets the vector's backing slice to the given
+// new one.
+func (vn *VecN) SetBackingSlice(newSlice []float32) {
+	vn.vec = newSlice
+}
+
+// Return the len of the vector's underlying slice.
+// This is not titled Len because it conflicts the package's
+// convention of calling the Norm the Len.
+func (vn *VecN) Size() int {
+	return len(vn.vec)
+}
+
+// Returns the cap of the vector's underlying slice.
+func (vn *VecN) Cap() int {
+	return cap(vn.vec)
+}
+
+// Sets the vector's size to n and zeroes out the vector.
+// If n is bigger than the vector's size, it will realloc.
+func (vn *VecN) Zero(n int) {
+	vn.Resize(n)
+	for i := range vn.vec {
+		vn.vec[i] = 0
+	}
+}
+
+// Adds vn and addend, storing the result in dst.
+// If dst does not have sufficient size it will be resized
+// Dst may be one of the other arguments. If dst is nil, it will be allocated.
+// The value returned is dst, for easier method chaining
+//
+// If vn and addend are not the same size, this function will add min(vn.Size(), addend.Size())
+// elements.
+func (vn *VecN) Add(dst *VecN, subtrahend *VecN) *VecN {
+	if vn == nil || subtrahend == nil {
+		return nil
+	}
+	size := intMin(len(vn.vec), len(subtrahend.vec))
+	dst = dst.Resize(size)
+
+	for i := 0; i < size; i++ {
+		dst.vec[i] = vn.vec[i] + subtrahend.vec[i]
+	}
+
+	return dst
+}
+
+// Subtracts addend from vn, storing the result in dst.
+// If dst does not have sufficient size it will be resized
+// Dst may be one of the other arguments. If dst is nil, it will be allocated.
+// The value returned is dst, for easier method chaining
+//
+// If vn and addend are not the same size, this function will add min(vn.Size(), addend.Size())
+// elements.
+func (vn *VecN) Sub(dst *VecN, addend *VecN) *VecN {
+	if vn == nil || addend == nil {
+		return nil
+	}
+	size := intMin(len(vn.vec), len(addend.vec))
+	dst = dst.Resize(size)
+
+	for i := 0; i < size; i++ {
+		dst.vec[i] = vn.vec[i] - addend.vec[i]
+	}
+
+	return dst
+}
+
+// Takes the binary cross product of vn and other, and stores it in dst.
+// If either vn or other are not of size 3 this function will panic
+//
+// If dst is not of sufficient size, or is nil, a new slice is allocated.
+// Dst is permitted to be one of the other arguments
+func (vn *VecN) Cross(dst *VecN, other *VecN) *VecN {
+	if vn == nil || other == nil {
+		return nil
+	}
+	if len(vn.vec) != 3 || len(other.vec) != 3 {
+		panic("Cannot take binary cross product of non-3D elements (7D cross product not implemented)")
+	}
+
+	dst = dst.Resize(3)
+	dst.vec[0], dst.vec[1], dst.vec[2] = vn.vec[1]*other.vec[2]-vn.vec[2]*other.vec[1], vn.vec[2]*other.vec[0]-vn.vec[0]*other.vec[2], vn.vec[0]*other.vec[1]-vn.vec[1]*other.vec[0]
+
+	return dst
+}
+
+func intMin(a, b int) int {
+	if a < b {
+		return a
+	}
+
+	return b
+}
+
+// Computes the dot product of two VecNs, if
+// the two vectors are not of the same length -- this
+// will return NaN.
+func (vn *VecN) Dot(other *VecN) float32 {
+	if vn == nil || other == nil || len(vn.vec) != len(other.vec) {
+		return float32(math.NaN())
+	}
+
+	var result float32 = 0.0
+	for i, el := range vn.vec {
+		result += el * other.vec[i]
+	}
+
+	return result
+}
+
+// Computes the vector length (also called the Norm) of the
+// vector. Equivalent to math.Sqrt(vn.Dot(vn)) with the appropriate
+// type conversions.
+//
+// If vn is nil, this returns NaN
+func (vn *VecN) Len() float32 {
+	if vn == nil {
+		return float32(math.NaN())
+	}
+	if len(vn.vec) == 0 {
+		return 0
+	}
+
+	return float32(math.Sqrt(float64(vn.Dot(vn))))
+}
+
+// Normalizes the vector and stores the result in dst, which
+// will be returned. Dst will be appropraitely resized to the
+// size of vn.
+//
+// The destination can be vn itself and nothing will go wrong.
+//
+// This is equivalent to vn.Mul(dst, 1/vn.Len())
+func (vn *VecN) Normalize(dst *VecN) *VecN {
+	if vn == nil {
+		return nil
+	}
+
+	return vn.Mul(dst, 1/vn.Len())
+}
+
+// Multiplied the vector by some scalar value and stores the result in dst, which
+// will be returned. Dst will be appropraitely resized to the
+// size of vn.
+//
+// The destination can be vn itself and nothing will go wrong.
+func (vn *VecN) Mul(dst *VecN, c float32) *VecN {
+	if vn == nil {
+		return nil
+	}
+	dst = dst.Resize(len(vn.vec))
+
+	for i, el := range vn.vec {
+		dst.vec[i] = el * c
+	}
+
+	return dst
+}
+
+// Performs the vector outer product between vn and v2.
+// The outer product is like a "reverse" dot product. Where the dot product
+// aligns both vectors with the "sized" part facing "inward" (Vec3*Vec3=Mat1x3*Mat3x1=Mat1x1=Scalar).
+// The outer product multiplied them with it facing "outward"
+// (Vec3*Vec3=Mat3x1*Mat1x3=Mat3x3).
+//
+// The matrix dst will be Reshaped to the correct size, if vn or v2 are nil,
+// this returns nil.
+func (vn *VecN) OuterProd(dst *MatMxN, v2 *VecN) *MatMxN {
+	if vn == nil || v2 == nil {
+		return nil
+	}
+
+	dst = dst.Reshape(len(vn.vec), len(v2.vec))
+
+	for c, el1 := range v2.vec {
+		for r, el2 := range vn.vec {
+			dst.Set(r, c, el1*el2)
+		}
+	}
+
+	return dst
+}
+
+func (vn *VecN) ApproxEqual(vn2 *VecN) bool {
+	if vn == nil || vn2 == nil || len(vn.vec) != len(vn2.vec) {
+		return false
+	}
+
+	for i, el := range vn.vec {
+		if !FloatEqual(el, vn2.vec[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (vn *VecN) ApproxEqualThreshold(vn2 *VecN, epsilon float32) bool {
+	if vn == nil || vn2 == nil || len(vn.vec) != len(vn2.vec) {
+		return false
+	}
+
+	for i, el := range vn.vec {
+		if !FloatEqualThreshold(el, vn2.vec[i], epsilon) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (vn *VecN) ApproxEqualFunc(vn2 *VecN, comp func(float32, float32) bool) bool {
+	if vn == nil || vn2 == nil || len(vn.vec) != len(vn2.vec) {
+		return false
+	}
+
+	for i, el := range vn.vec {
+		if !comp(el, vn2.vec[i]) {
+			return false
+		}
+	}
+
+	return true
+}
+
+func (vn *VecN) Vec2() Vec2 {
+	raw := vn.Raw()
+	return Vec2{raw[0], raw[1]}
+}
+
+func (vn *VecN) Vec3() Vec3 {
+	raw := vn.Raw()
+	return Vec3{raw[0], raw[1], raw[2]}
+}
+
+func (vn *VecN) Vec4() Vec4 {
+	raw := vn.Raw()
+	return Vec4{raw[0], raw[1], raw[2], raw[3]}
+}

vendor/github.com/go-gl/mathgl/mgl32/vector.go

@@ -0,0 +1,562 @@
+// Copyright 2014 The go-gl/mathgl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// This file is generated by codegen.go; DO NOT EDIT
+// Edit vector.tmpl and run "go generate" to make changes.
+
+package mgl32
+
+import (
+	"golang.org/x/image/math/f32"
+	"math"
+)
+
+type Vec2 f32.Vec2
+type Vec3 f32.Vec3
+type Vec4 f32.Vec4
+
+func (v Vec2) Vec3(z float32) Vec3 {
+	return Vec3{v[0], v[1], z}
+}
+
+func (v Vec2) Vec4(z, w float32) Vec4 {
+	return Vec4{v[0], v[1], z, w}
+}
+
+func (v Vec3) Vec2() Vec2 {
+	return Vec2{v[0], v[1]}
+}
+
+func (v Vec3) Vec4(w float32) Vec4 {
+	return Vec4{v[0], v[1], v[2], w}
+}
+
+func (v Vec4) Vec2() Vec2 {
+	return Vec2{v[0], v[1]}
+}
+
+func (v Vec4) Vec3() Vec3 {
+	return Vec3{v[0], v[1], v[2]}
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec2) Elem() (x, y float32) {
+	return v[0], v[1]
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec3) Elem() (x, y, z float32) {
+	return v[0], v[1], v[2]
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec4) Elem() (x, y, z, w float32) {
+	return v[0], v[1], v[2], v[3]
+}
+
+// The vector cross product is an operation only defined on 3D vectors. It is equivalent to
+// Vec3{v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0]}.
+// Another interpretation is that it's the vector whose magnitude is |v1||v2|sin(theta)
+// where theta is the angle between v1 and v2.
+//
+// The cross product is most often used for finding surface normals. The cross product of vectors
+// will generate a vector that is perpendicular to the plane they form.
+//
+// Technically, a generalized cross product exists as an "(N-1)ary" operation
+// (that is, the 4D cross product requires 3 4D vectors). But the binary
+// 3D (and 7D) cross product is the most important. It can be considered
+// the area of a parallelogram with sides v1 and v2.
+//
+// Like the dot product, the cross product is roughly a measure of directionality.
+// Two normalized perpendicular vectors will return a vector with a magnitude of
+// 1.0 or -1.0 and two parallel vectors will return a vector with magnitude 0.0.
+// The cross product is "anticommutative" meaning v1.Cross(v2) = -v2.Cross(v1),
+// this property can be useful to know when finding normals,
+// as taking the wrong cross product can lead to the opposite normal of the one you want.
+func (v1 Vec3) Cross(v2 Vec3) Vec3 {
+	return Vec3{v1[1]*v2[2] - v1[2]*v2[1], v1[2]*v2[0] - v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0]}
+}
+
+// Add performs element-wise addition between two vectors. It is equivalent to iterating
+// over every element of v1 and adding the corresponding element of v2 to it.
+func (v1 Vec2) Add(v2 Vec2) Vec2 {
+	return Vec2{v1[0] + v2[0], v1[1] + v2[1]}
+}
+
+// Sub performs element-wise subtraction between two vectors. It is equivalent to iterating
+// over every element of v1 and subtracting the corresponding element of v2 from it.
+func (v1 Vec2) Sub(v2 Vec2) Vec2 {
+	return Vec2{v1[0] - v2[0], v1[1] - v2[1]}
+}
+
+// Mul performs a scalar multiplication between the vector and some constant value
+// c. This is equivalent to iterating over every vector element and multiplying by c.
+func (v1 Vec2) Mul(c float32) Vec2 {
+	return Vec2{v1[0] * c, v1[1] * c}
+}
+
+// Dot returns the dot product of this vector with another. There are multiple ways
+// to describe this value. One is the multiplication of their lengths and cos(theta) where
+// theta is the angle between the vectors: v1.v2 = |v1||v2|cos(theta).
+//
+// The other (and what is actually done) is the sum of the element-wise multiplication of all
+// elements. So for instance, two Vec3s would yield v1.x * v2.x + v1.y * v2.y + v1.z * v2.z.
+//
+// This means that the dot product of a vector and itself is the square of its Len (within
+// the bounds of floating points error).
+//
+// The dot product is roughly a measure of how closely two vectors are to pointing in the same
+// direction. If both vectors are normalized, the value will be -1 for opposite pointing,
+// one for same pointing, and 0 for perpendicular vectors.
+func (v1 Vec2) Dot(v2 Vec2) float32 {
+	return v1[0]*v2[0] + v1[1]*v2[1]
+}
+
+// Len returns the vector's length. Note that this is NOT the dimension of
+// the vector (len(v)), but the mathematical length. This is equivalent to the square
+// root of the sum of the squares of all elements. E.G. for a Vec2 it's
+// math.Hypot(v[0], v[1]).
+func (v1 Vec2) Len() float32 {
+
+	return float32(math.Hypot(float64(v1[0]), float64(v1[1])))
+
+}
+
+// Normalize normalizes the vector. Normalization is (1/|v|)*v,
+// making this equivalent to v.Scale(1/v.Len()). If the len is 0.0,
+// this function will return an infinite value for all elements due
+// to how floating point division works in Go (n/0.0 = math.Inf(Sign(n))).
+//
+// Normalization makes a vector's Len become 1.0 (within the margin of floating point error),
+// while maintaining its directionality.
+//
+// (Can be seen here: http://play.golang.org/p/Aaj7SnbqIp )
+func (v1 Vec2) Normalize() Vec2 {
+	l := 1.0 / v1.Len()
+	return Vec2{v1[0] * l, v1[1] * l}
+}
+
+// ApproxEqual takes in a vector and does an element-wise
+// approximate float comparison as if FloatEqual had been used
+func (v1 Vec2) ApproxEqual(v2 Vec2) bool {
+	for i := range v1 {
+		if !FloatEqual(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxThresholdEq takes in a threshold for comparing two floats, and uses it to do an
+// element-wise comparison of the vector to another.
+func (v1 Vec2) ApproxEqualThreshold(v2 Vec2, threshold float32) bool {
+	for i := range v1 {
+		if !FloatEqualThreshold(v1[i], v2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxFuncEq takes in a func that compares two floats, and uses it to do an element-wise
+// comparison of the vector to another. This is intended to be used with FloatEqualFunc
+func (v1 Vec2) ApproxFuncEqual(v2 Vec2, eq func(float32, float32) bool) bool {
+	for i := range v1 {
+		if !eq(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec2) X() float32 {
+	return v[0]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec2) Y() float32 {
+	return v[1]
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 2x2 matrix. E.G. a Vec2 * Vec2 = Mat2.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec2) OuterProd2(v2 Vec2) Mat2 {
+	return Mat2{v1[0] * v2[0], v1[1] * v2[0], v1[0] * v2[1], v1[1] * v2[1]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 2x3 matrix. E.G. a Vec2 * Vec3 = Mat2x3.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec2) OuterProd3(v2 Vec3) Mat2x3 {
+	return Mat2x3{v1[0] * v2[0], v1[1] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[0] * v2[2], v1[1] * v2[2]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 2x4 matrix. E.G. a Vec2 * Vec4 = Mat2x4.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec2) OuterProd4(v2 Vec4) Mat2x4 {
+	return Mat2x4{v1[0] * v2[0], v1[1] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[0] * v2[3], v1[1] * v2[3]}
+}
+
+// Add performs element-wise addition between two vectors. It is equivalent to iterating
+// over every element of v1 and adding the corresponding element of v2 to it.
+func (v1 Vec3) Add(v2 Vec3) Vec3 {
+	return Vec3{v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2]}
+}
+
+// Sub performs element-wise subtraction between two vectors. It is equivalent to iterating
+// over every element of v1 and subtracting the corresponding element of v2 from it.
+func (v1 Vec3) Sub(v2 Vec3) Vec3 {
+	return Vec3{v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2]}
+}
+
+// Mul performs a scalar multiplication between the vector and some constant value
+// c. This is equivalent to iterating over every vector element and multiplying by c.
+func (v1 Vec3) Mul(c float32) Vec3 {
+	return Vec3{v1[0] * c, v1[1] * c, v1[2] * c}
+}
+
+// Dot returns the dot product of this vector with another. There are multiple ways
+// to describe this value. One is the multiplication of their lengths and cos(theta) where
+// theta is the angle between the vectors: v1.v2 = |v1||v2|cos(theta).
+//
+// The other (and what is actually done) is the sum of the element-wise multiplication of all
+// elements. So for instance, two Vec3s would yield v1.x * v2.x + v1.y * v2.y + v1.z * v2.z.
+//
+// This means that the dot product of a vector and itself is the square of its Len (within
+// the bounds of floating points error).
+//
+// The dot product is roughly a measure of how closely two vectors are to pointing in the same
+// direction. If both vectors are normalized, the value will be -1 for opposite pointing,
+// one for same pointing, and 0 for perpendicular vectors.
+func (v1 Vec3) Dot(v2 Vec3) float32 {
+	return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
+}
+
+// Len returns the vector's length. Note that this is NOT the dimension of
+// the vector (len(v)), but the mathematical length. This is equivalent to the square
+// root of the sum of the squares of all elements. E.G. for a Vec2 it's
+// math.Hypot(v[0], v[1]).
+func (v1 Vec3) Len() float32 {
+
+	return float32(math.Sqrt(float64(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2])))
+
+}
+
+// Normalize normalizes the vector. Normalization is (1/|v|)*v,
+// making this equivalent to v.Scale(1/v.Len()). If the len is 0.0,
+// this function will return an infinite value for all elements due
+// to how floating point division works in Go (n/0.0 = math.Inf(Sign(n))).
+//
+// Normalization makes a vector's Len become 1.0 (within the margin of floating point error),
+// while maintaining its directionality.
+//
+// (Can be seen here: http://play.golang.org/p/Aaj7SnbqIp )
+func (v1 Vec3) Normalize() Vec3 {
+	l := 1.0 / v1.Len()
+	return Vec3{v1[0] * l, v1[1] * l, v1[2] * l}
+}
+
+// ApproxEqual takes in a vector and does an element-wise
+// approximate float comparison as if FloatEqual had been used
+func (v1 Vec3) ApproxEqual(v2 Vec3) bool {
+	for i := range v1 {
+		if !FloatEqual(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxThresholdEq takes in a threshold for comparing two floats, and uses it to do an
+// element-wise comparison of the vector to another.
+func (v1 Vec3) ApproxEqualThreshold(v2 Vec3, threshold float32) bool {
+	for i := range v1 {
+		if !FloatEqualThreshold(v1[i], v2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxFuncEq takes in a func that compares two floats, and uses it to do an element-wise
+// comparison of the vector to another. This is intended to be used with FloatEqualFunc
+func (v1 Vec3) ApproxFuncEqual(v2 Vec3, eq func(float32, float32) bool) bool {
+	for i := range v1 {
+		if !eq(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec3) X() float32 {
+	return v[0]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec3) Y() float32 {
+	return v[1]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec3) Z() float32 {
+	return v[2]
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 3x2 matrix. E.G. a Vec3 * Vec2 = Mat3x2.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec3) OuterProd2(v2 Vec2) Mat3x2 {
+	return Mat3x2{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 3x3 matrix. E.G. a Vec3 * Vec3 = Mat3.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec3) OuterProd3(v2 Vec3) Mat3 {
+	return Mat3{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 3x4 matrix. E.G. a Vec3 * Vec4 = Mat3x4.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec3) OuterProd4(v2 Vec4) Mat3x4 {
+	return Mat3x4{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2], v1[0] * v2[3], v1[1] * v2[3], v1[2] * v2[3]}
+}
+
+// Add performs element-wise addition between two vectors. It is equivalent to iterating
+// over every element of v1 and adding the corresponding element of v2 to it.
+func (v1 Vec4) Add(v2 Vec4) Vec4 {
+	return Vec4{v1[0] + v2[0], v1[1] + v2[1], v1[2] + v2[2], v1[3] + v2[3]}
+}
+
+// Sub performs element-wise subtraction between two vectors. It is equivalent to iterating
+// over every element of v1 and subtracting the corresponding element of v2 from it.
+func (v1 Vec4) Sub(v2 Vec4) Vec4 {
+	return Vec4{v1[0] - v2[0], v1[1] - v2[1], v1[2] - v2[2], v1[3] - v2[3]}
+}
+
+// Mul performs a scalar multiplication between the vector and some constant value
+// c. This is equivalent to iterating over every vector element and multiplying by c.
+func (v1 Vec4) Mul(c float32) Vec4 {
+	return Vec4{v1[0] * c, v1[1] * c, v1[2] * c, v1[3] * c}
+}
+
+// Dot returns the dot product of this vector with another. There are multiple ways
+// to describe this value. One is the multiplication of their lengths and cos(theta) where
+// theta is the angle between the vectors: v1.v2 = |v1||v2|cos(theta).
+//
+// The other (and what is actually done) is the sum of the element-wise multiplication of all
+// elements. So for instance, two Vec3s would yield v1.x * v2.x + v1.y * v2.y + v1.z * v2.z.
+//
+// This means that the dot product of a vector and itself is the square of its Len (within
+// the bounds of floating points error).
+//
+// The dot product is roughly a measure of how closely two vectors are to pointing in the same
+// direction. If both vectors are normalized, the value will be -1 for opposite pointing,
+// one for same pointing, and 0 for perpendicular vectors.
+func (v1 Vec4) Dot(v2 Vec4) float32 {
+	return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] + v1[3]*v2[3]
+}
+
+// Len returns the vector's length. Note that this is NOT the dimension of
+// the vector (len(v)), but the mathematical length. This is equivalent to the square
+// root of the sum of the squares of all elements. E.G. for a Vec2 it's
+// math.Hypot(v[0], v[1]).
+func (v1 Vec4) Len() float32 {
+
+	return float32(math.Sqrt(float64(v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2] + v1[3]*v1[3])))
+
+}
+
+// Normalize normalizes the vector. Normalization is (1/|v|)*v,
+// making this equivalent to v.Scale(1/v.Len()). If the len is 0.0,
+// this function will return an infinite value for all elements due
+// to how floating point division works in Go (n/0.0 = math.Inf(Sign(n))).
+//
+// Normalization makes a vector's Len become 1.0 (within the margin of floating point error),
+// while maintaining its directionality.
+//
+// (Can be seen here: http://play.golang.org/p/Aaj7SnbqIp )
+func (v1 Vec4) Normalize() Vec4 {
+	l := 1.0 / v1.Len()
+	return Vec4{v1[0] * l, v1[1] * l, v1[2] * l, v1[3] * l}
+}
+
+// ApproxEqual takes in a vector and does an element-wise
+// approximate float comparison as if FloatEqual had been used
+func (v1 Vec4) ApproxEqual(v2 Vec4) bool {
+	for i := range v1 {
+		if !FloatEqual(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxThresholdEq takes in a threshold for comparing two floats, and uses it to do an
+// element-wise comparison of the vector to another.
+func (v1 Vec4) ApproxEqualThreshold(v2 Vec4, threshold float32) bool {
+	for i := range v1 {
+		if !FloatEqualThreshold(v1[i], v2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxFuncEq takes in a func that compares two floats, and uses it to do an element-wise
+// comparison of the vector to another. This is intended to be used with FloatEqualFunc
+func (v1 Vec4) ApproxFuncEqual(v2 Vec4, eq func(float32, float32) bool) bool {
+	for i := range v1 {
+		if !eq(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) X() float32 {
+	return v[0]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) Y() float32 {
+	return v[1]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) Z() float32 {
+	return v[2]
+}
+
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v Vec4) W() float32 {
+	return v[3]
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 4x2 matrix. E.G. a Vec4 * Vec2 = Mat4x2.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec4) OuterProd2(v2 Vec2) Mat4x2 {
+	return Mat4x2{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[3] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[3] * v2[1]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 4x3 matrix. E.G. a Vec4 * Vec3 = Mat4x3.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec4) OuterProd3(v2 Vec3) Mat4x3 {
+	return Mat4x3{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[3] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[3] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2], v1[3] * v2[2]}
+}
+
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// 4x4 matrix. E.G. a Vec4 * Vec4 = Mat4.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 Vec4) OuterProd4(v2 Vec4) Mat4 {
+	return Mat4{v1[0] * v2[0], v1[1] * v2[0], v1[2] * v2[0], v1[3] * v2[0], v1[0] * v2[1], v1[1] * v2[1], v1[2] * v2[1], v1[3] * v2[1], v1[0] * v2[2], v1[1] * v2[2], v1[2] * v2[2], v1[3] * v2[2], v1[0] * v2[3], v1[1] * v2[3], v1[2] * v2[3], v1[3] * v2[3]}
+}

vendor/github.com/go-gl/mathgl/mgl32/vector.tmpl

@@ -0,0 +1,207 @@
+// Copyright 2014 The go-gl/mathgl Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// <<.Comment>>
+// Edit <<.TemplateName>> and run "go generate" to make changes.
+
+package mgl32
+
+import (
+	"golang.org/x/image/math/f32"
+	"math"
+)
+
+type Vec2 f32.Vec2
+type Vec3 f32.Vec3
+type Vec4 f32.Vec4
+
+func (v Vec2) Vec3(z float32) Vec3 {
+	return Vec3{v[0], v[1], z}
+}
+
+func (v Vec2) Vec4(z, w float32) Vec4 {
+	return Vec4{v[0], v[1], z, w}
+}
+
+func (v Vec3) Vec2() Vec2 {
+	return Vec2{v[0], v[1]}
+}
+
+func (v Vec3) Vec4(w float32) Vec4 {
+	return Vec4{v[0], v[1], v[2], w}
+}
+
+func (v Vec4) Vec2() Vec2 {
+	return Vec2{v[0], v[1]}
+}
+
+func (v Vec4) Vec3() Vec3 {
+	return Vec3{v[0], v[1], v[2]}
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec2) Elem() (x, y float32) {
+	return v[0], v[1]
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec3) Elem() (x, y, z float32) {
+	return v[0], v[1], v[2]
+}
+
+// Elem extracts the elements of the vector for direct value assignment.
+func (v Vec4) Elem() (x, y, z, w float32) {
+	return v[0], v[1], v[2], v[3]
+}
+
+// The vector cross product is an operation only defined on 3D vectors. It is equivalent to
+// Vec3{v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0]}.
+// Another interpretation is that it's the vector whose magnitude is |v1||v2|sin(theta)
+// where theta is the angle between v1 and v2.
+//
+// The cross product is most often used for finding surface normals. The cross product of vectors
+// will generate a vector that is perpendicular to the plane they form.
+//
+// Technically, a generalized cross product exists as an "(N-1)ary" operation
+// (that is, the 4D cross product requires 3 4D vectors). But the binary
+// 3D (and 7D) cross product is the most important. It can be considered
+// the area of a parallelogram with sides v1 and v2.
+//
+// Like the dot product, the cross product is roughly a measure of directionality.
+// Two normalized perpendicular vectors will return a vector with a magnitude of
+// 1.0 or -1.0 and two parallel vectors will return a vector with magnitude 0.0.
+// The cross product is "anticommutative" meaning v1.Cross(v2) = -v2.Cross(v1),
+// this property can be useful to know when finding normals,
+// as taking the wrong cross product can lead to the opposite normal of the one you want.
+func (v1 Vec3) Cross(v2 Vec3) Vec3 {
+	return Vec3{v1[1]*v2[2] - v1[2]*v2[1], v1[2]*v2[0] - v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0]}
+}
+
+
+<</* Common functions for all vectors */>>
+<<range $m := enum 2 3 4>>
+<<$type := typename $m 1>>
+
+// Add performs element-wise addition between two vectors. It is equivalent to iterating
+// over every element of v1 and adding the corresponding element of v2 to it.
+func (v1 <<$type>>) Add(v2 <<$type>>) <<$type>> {
+	return <<$type>>{<<range $i := iter 0 $m>> v1[<<$i>>] + v2[<<$i>>], <<end>>}
+}
+
+// Sub performs element-wise subtraction between two vectors. It is equivalent to iterating
+// over every element of v1 and subtracting the corresponding element of v2 from it.
+func (v1 <<$type>>) Sub(v2 <<$type>>) <<$type>> {
+	return <<$type>>{<<range $i := iter 0 $m>> v1[<<$i>>] - v2[<<$i>>], <<end>>}
+}
+
+// Mul performs a scalar multiplication between the vector and some constant value
+// c. This is equivalent to iterating over every vector element and multiplying by c.
+func (v1 <<$type>>) Mul(c float32) <<$type>> {
+	return <<$type>>{<<range $i := iter 0 $m>> v1[<<$i>>] * c, <<end>>}
+}
+
+// Dot returns the dot product of this vector with another. There are multiple ways
+// to describe this value. One is the multiplication of their lengths and cos(theta) where
+// theta is the angle between the vectors: v1.v2 = |v1||v2|cos(theta).
+//
+// The other (and what is actually done) is the sum of the element-wise multiplication of all
+// elements. So for instance, two Vec3s would yield v1.x * v2.x + v1.y * v2.y + v1.z * v2.z.
+//
+// This means that the dot product of a vector and itself is the square of its Len (within
+// the bounds of floating points error).
+//
+// The dot product is roughly a measure of how closely two vectors are to pointing in the same
+// direction. If both vectors are normalized, the value will be -1 for opposite pointing,
+// one for same pointing, and 0 for perpendicular vectors.
+func (v1 <<$type>>) Dot(v2 <<$type>>) float32 {
+	return <<range $i := iter 0 $m>><<sep "+" $i>> v1[<<$i>>]*v2[<<$i>>] <<end>>
+}
+
+// Len returns the vector's length. Note that this is NOT the dimension of
+// the vector (len(v)), but the mathematical length. This is equivalent to the square
+// root of the sum of the squares of all elements. E.G. for a Vec2 it's
+// math.Hypot(v[0], v[1]).
+func (v1 <<$type>>) Len() float32 {
+	<<if eq $m 2 >>
+	return float32(math.Hypot(float64(v1[0]), float64(v1[1])))
+	<<else>>
+	return float32(math.Sqrt(float64(<<repeat $m "v1[%d]*v1[%d]" "+">>)))
+	<<end>>
+}
+
+// Normalize normalizes the vector. Normalization is (1/|v|)*v,
+// making this equivalent to v.Scale(1/v.Len()). If the len is 0.0,
+// this function will return an infinite value for all elements due
+// to how floating point division works in Go (n/0.0 = math.Inf(Sign(n))).
+//
+// Normalization makes a vector's Len become 1.0 (within the margin of floating point error),
+// while maintaining its directionality.
+//
+// (Can be seen here: http://play.golang.org/p/Aaj7SnbqIp )
+func (v1 <<$type>>) Normalize() <<$type>> {
+	l := 1.0 / v1.Len()
+	return <<$type>>{<<range $i := iter 0 $m>>v1[<<$i>>] * l,<<end>>}
+}
+
+// ApproxEqual takes in a vector and does an element-wise
+// approximate float comparison as if FloatEqual had been used
+func (v1 <<$type>>) ApproxEqual(v2 <<$type>>) bool {
+	for i := range v1 {
+		if !FloatEqual(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxThresholdEq takes in a threshold for comparing two floats, and uses it to do an
+// element-wise comparison of the vector to another.
+func (v1 <<$type>>) ApproxEqualThreshold(v2 <<$type>>, threshold float32) bool {
+	for i := range v1 {
+		if !FloatEqualThreshold(v1[i], v2[i], threshold) {
+			return false
+		}
+	}
+	return true
+}
+
+// ApproxFuncEq takes in a func that compares two floats, and uses it to do an element-wise
+// comparison of the vector to another. This is intended to be used with FloatEqualFunc
+func (v1 <<$type>>) ApproxFuncEqual(v2 <<$type>>, eq func(float32, float32) bool) bool {
+	for i := range v1 {
+		if !eq(v1[i], v2[i]) {
+			return false
+		}
+	}
+	return true
+}
+
+<<range $i := iter 0 $m>>
+// This is an element access func, it is equivalent to v[n] where
+// n is some valid index. The mappings are XYZW (X=0, Y=1 etc). Benchmarks
+// show that this is more or less as fast as direct acces, probably due to
+// inlining, so use v[0] or v.X() depending on personal preference.
+func (v <<$type>>) <<elementname $i>>() float32 {
+	return v[<<$i>>]
+}
+<<end>>
+
+<<range $n := enum 2 3 4>>
+// Does the vector outer product
+// of two vectors. The outer product produces an
+// <<$m>>x<<$n>> matrix. E.G. a Vec<<$m>> * Vec<<$n>> = <<typename $m $n>>.
+//
+// The outer product can be thought of as the "opposite"
+// of the Dot product. The Dot product treats both vectors like matrices
+// oriented such that the left one has N columns and the right has N rows.
+// So Vec3.Vec3 = Mat1x3*Mat3x1 = Mat1 = Scalar.
+//
+// The outer product orients it so they're facing "outward": Vec2*Vec3
+// = Mat2x1*Mat1x3 = Mat2x3.
+func (v1 <<$type>>) OuterProd<<$n>>(v2 <<typename $n 1>>) <<typename $m $n>> {
+	return <<typename $m $n>>{<<range $i := matiter $m $n>>v1[<<$i.M>>] * v2[<<$i.N>>], <<end>>}
+}
+<<end>>
+
+<<end>> <</* range $m */>>