diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 2830b81..a853593 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -6,12 +6,16 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, windows-latest]
+        nim-version: ['1.4.0', '1.4.x', 'stable']
 
     runs-on: ${{ matrix.os }}
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - uses: jiro4989/setup-nim-action@v1
+      with:
+        nim-version: ${{ matrix.nim-version }}
+        repo-token: ${{ secrets.GITHUB_TOKEN }}
     - run: nimble test -y
     - run: nimble test --gc:orc -y
     - run: nimble test -y -d:vmathObjBased
diff --git a/README.md b/README.md
index 2cf4f8c..57fd69b 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 <img src="docs/banner.png">
 
-# VMath - 2d and 3d vector math.
+# VMath - 2D and 3D vector math.
 
 `nimble install vmath`
 
@@ -10,6 +10,8 @@
 
 This library has no dependencies other than the Nim standard library.
 
+Supports c, cpp and js backend.
+
 ## About
 
 Your one stop shop for vector math routines for 2d and 3d graphics.
@@ -89,6 +91,44 @@ vmathObjArrayBased ................ 73.968 ms     74.292 ms    ±0.631   x100
 * [3d Ray Trace Benchmark](tests/bench_raytracer.nim)
 * [2d SVG Render Benchmark](https://github.com/treeform/pixie/blob/master/tests/bench_svg.nim)
 
+## Zmod - GLSL mod
+
+GLSL uses a different type of float point mod. Because mod is a Nim keyword please use `zmod` when you need GLSL `mod` behavior.
+
+## Coordinate System
+
+Right-hand z-forward coordinate system
+
+This is the same system used in the GLTF file format.
+
+> glTF uses a right-handed coordinate system.
+> glTF defines +Y as up, +Z as forward, and -X as right;
+> the front of a glTF asset faces +Z.
+
+[glTF Spec 2.0](https://registry.khronos.org/glTF/specs/2.0/glTF-2.0.html#coordinate-system-and-units)
+
+## OpenGL matrix column-major notation.
+
+> [9.005](https://www.opengl.org/archives/resources/faq/technical/transformations.htm) For programming purposes, OpenGL matrices are 16-value arrays with base vectors laid out contiguously in memory. The translation components occupy the 13th, 14th, and 15th elements of the 16-element matrix, where indices are numbered from 1 to 16 as described in section 2.11.2 of the [OpenGL 2.1 Specification](https://registry.khronos.org/OpenGL/specs/gl/glspec21.pdf).
+>
+> Sadly, the use of column-major format in the spec and blue book has resulted in endless confusion in the OpenGL programming community. Column-major notation suggests that matrices are not laid out in memory as a programmer would expect.
+
+OpenGL/GLSL/vmath vs Math/Specification notation:
+```
+        mat4([
+          a, b, c, 0,              | a d g x |
+          d, e, f, 0,              | b e h y |
+          g, h, i, 0,              | c f i z |
+          x, y, z, 1               | 0 0 0 1 |
+        ])
+```
+
+# 1.x.x to 2.0.0 vmath breaking changes:
+* New right-hand-z-forward cordate system and functions that care about
+coordinate system where moved there.
+* deprecated `lookAt()` please use `toAngles()`/`fromAngles()` instead.
+* deprecated `fractional()` use `frac()` instead.
+
 # 0.x.x to 1.0.0 vmath breaking changes:
 
 * `vec3(v)` no longer works please use `vec3(v.x, v.y, 0)` instead.
diff --git a/src/vmath.nim b/src/vmath.nim
index 26eaf0b..b1bb102 100644
--- a/src/vmath.nim
+++ b/src/vmath.nim
@@ -1,22 +1,20 @@
 ##[
 
-This library has no dependencies other than the Nim standard libarary.
-
 Your one stop shop for vector math routines for 2d and 3d graphics.
 
 * Pure Nim with no dependencies.
 * Very similar to GLSL Shader Language with extra stuff.
 * Extensively benchmarked.
 
-====== =========== ===================================================
+====== =========== =================================================
 Type   Constructor Description
-====== =========== ===================================================
-BVec#  bvec#       a vector of booleans
-IVec#  ivec#       a vector of signed integers
-UVec#  uvec#       a vector of unsigned integers
-Vec#   vec#        a vector of single-precision floating-point numbers
-DVec#  dvec#       a vector of double-precision floating-point numbers
-====== =========== ===================================================
+====== =========== =================================================
+BVec#  bvec#       vector of booleans
+IVec#  ivec#       vector of signed integers
+UVec#  uvec#       vector of unsigned integers
+Vec#   vec#        vector of single-precision floating-point numbers
+DVec#  dvec#       vector of double-precision floating-point numbers
+====== =========== =================================================
 
 You can use these constructors to make them:
 
@@ -406,13 +404,17 @@ proc quantize*[T: SomeFloat](v, n: T): T =
   ## Makes v be multiple of n. Rounding to integer quantize by 1.0.
   trunc(v / n) * n
 
-proc fractional*[T: SomeFloat](v: T): T =
+proc fract*[T: SomeFloat](v: T): T =
   ## Returns fractional part of a number.
   ## 3.14 -> 0.14
   ## -3.14 -> 0.14
   result = abs(v)
   result = result - trunc(result)
 
+proc fractional*[T: SomeFloat](v: T): T {.deprecated: "Use frac() insetad"} =
+  ## Returns fractional part of a number.
+  fract(v)
+
 proc inversesqrt*[T: float32|float64](v: T): T =
   ## Returns inverse square root.
   1/sqrt(v)
@@ -425,19 +427,21 @@ proc mix*[T: SomeFloat](a, b, v: T): T =
   v * (b - a) + a
 
 proc fixAngle*[T: SomeFloat](angle: T): T =
-  ## Make angle be from -PI to PI radians.
+  ## Normalize the angle to be from -PI to PI radians.
   result = angle
   while result > PI:
     result -= PI * 2
-  while result < -PI:
+  while result <= -PI:
     result += PI * 2
 
 proc angleBetween*[T: SomeFloat](a, b: T): T =
   ## Angle between angle a and angle b.
+  ## All angles assume radians.
   fixAngle(b - a)
 
 proc turnAngle*[T: SomeFloat](a, b, speed: T): T =
   ## Move from angle a to angle b with step of v.
+  ## All angles assume radians.
   var
     turn = fixAngle(b - a)
   if abs(turn) < speed:
@@ -725,6 +729,11 @@ genMathFn(sqrt)
 genMathFn(floor)
 genMathFn(ceil)
 genMathFn(abs)
+genMathFn(trunc)
+genMathFn(fract)
+genMathFn(quantize)
+genMathFn(toRadians)
+genMathFn(toDegrees)
 
 template genBoolFn(fn, op: untyped) =
   proc fn*[T](a, b: GVec2[T]): BVec2 =
@@ -970,10 +979,42 @@ proc `pos=`*[T](a: var GMat3[T], pos: GVec2[T]) =
   a[2, 0] = pos.x
   a[2, 1] = pos.y
 
+proc forward*[T](a: GMat4[T]): GVec3[T] {.inline.} =
+  ## Vector facing +Z.
+  result.x = a[2, 0]
+  result.y = a[2, 1]
+  result.z = a[2, 2]
+
+proc back*[T](a: GMat4[T]): GVec3[T] {.inline.} =
+  ## Vector facing -Z.
+  -a.forward()
+
+proc left*[T](a: GMat4[T]): GVec3[T] {.inline.} =
+  ## Vector facing +X.
+  result.x = a[0, 0]
+  result.y = a[0, 1]
+  result.z = a[0, 2]
+
+proc right*[T](a: GMat4[T]): GVec3[T] {.inline.} =
+  ## Vector facing -X.
+  -a.left()
+
+proc up*[T](a: GMat4[T]): GVec3[T] {.inline.} =
+  ## Vector facing +Y.
+  result.x = a[1, 0]
+  result.y = a[1, 1]
+  result.z = a[1, 2]
+
+proc down*[T](a: GMat4[T]): GVec3[T] {.inline.} =
+  ## Vector facing -X.
+  -a.up()
+
 proc pos*[T](a: GMat4[T]): GVec3[T] =
+  ## Position of the matrix.
   gvec3[T](a[3].x, a[3].y, a[3].z)
 
 proc `pos=`*[T](a: var GMat4[T], pos: GVec3[T]) =
+  ## See the position of the matrix.
   a[3, 0] = pos.x
   a[3, 1] = pos.y
   a[3, 2] = pos.z
@@ -1244,7 +1285,8 @@ proc translate*[T](v: GVec3[T]): GMat4[T] =
   )
 
 proc rotate*[T](angle: T): GMat3[T] =
-  ## Create a rotation matrix by an angle.
+  ## Create a 2D rotation matrix by an angle.
+  ## All angles assume radians.
   let
     sin = sin(angle)
     cos = cos(angle)
@@ -1254,22 +1296,126 @@ proc rotate*[T](angle: T): GMat3[T] =
     0, 0, 1
   )
 
-proc hrp*[T](m: GMat4[T]): GVec3[T] =
-  ## Return heading, rotation and pivot of a matrix.
-  var heading, pitch, roll: float32
-  if m[1] > 0.998: # singularity at north pole
-    heading = arctan2(m[2], m[10])
-    pitch = PI / 2
-    roll = 0
-  elif m[1] < -0.998: # singularity at south pole
-    heading = arctan2(m[2], m[10])
-    pitch = -PI / 2
-    roll = 0
+proc rotationOnly*[T](a: GMat4[T]): GMat4[T] {.inline.} =
+  ## Clears the positional component and returns rotation only.
+  ## Assumes matrix has not been scaled.
+  result = a
+  result.pos = gvec3(0, 0, 0)
+
+proc rotateX*[T](angle: T): GMat4[T] =
+  ## Return a rotation matrix around X with angle.
+  ## All angles assume radians.
+  result[0, 0] = 1
+  result[0, 1] = 0
+  result[0, 2] = 0
+  result[0, 3] = 0
+
+  result[1, 0] = 0
+  result[1, 1] = cos(angle)
+  result[1, 2] = -sin(angle)
+  result[1, 3] = 0
+
+  result[2, 0] = 0
+  result[2, 1] = sin(angle)
+  result[2, 2] = cos(angle)
+  result[2, 3] = 0
+
+  result[3, 0] = 0
+  result[3, 1] = 0
+  result[3, 2] = 0
+  result[3, 3] = 1
+
+proc rotateY*[T](angle: T): GMat4[T] =
+  ## Return a rotation matrix around Y with angle.
+  ## All angles assume radians.
+  result[0, 0] = cos(angle)
+  result[0, 1] = 0
+  result[0, 2] = sin(angle)
+  result[0, 3] = 0
+
+  result[1, 0] = 0
+  result[1, 1] = 1
+  result[1, 2] = 0
+  result[1, 3] = 0
+
+  result[2, 0] = -sin(angle)
+  result[2, 1] = 0
+  result[2, 2] = cos(angle)
+  result[2, 3] = 0
+
+  result[3, 0] = 0
+  result[3, 1] = 0
+  result[3, 2] = 0
+  result[3, 3] = 1
+
+proc rotateZ*[T](angle: T): GMat4[T] =
+  ## Return a rotation matrix around Z with angle.
+  ## All angles assume radians.
+  result[0, 0] = cos(angle)
+  result[0, 1] = -sin(angle)
+  result[0, 2] = 0
+  result[0, 3] = 0
+
+  result[1, 0] = sin(angle)
+  result[1, 1] = cos(angle)
+  result[1, 2] = 0
+  result[1, 3] = 0
+
+  result[2, 0] = 0
+  result[2, 1] = 0
+  result[2, 2] = 1
+  result[2, 3] = 0
+
+  result[3, 0] = 0
+  result[3, 1] = 0
+  result[3, 2] = 0
+  result[3, 3] = 1
+
+proc toAngles*[T](a: GVec3[T]): GVec3[T] =
+  ## Given a 3d vector, computes Euler angles: pitch and yaw
+  ##   pitch (x rotation)
+  ##   yaw (y rotation)
+  ##   roll (z rotation) - always 0 in vector case
+  ## All angles assume radians.
+  if a == gvec3[T](T(0), T(0), T(0)):
+    return
+  let
+    yaw = -arctan2(a.x, a.z)
+    pitch = -arctan2(sqrt(a.x*a.x + a.z*a.z), a.y) + T(PI/2)
+  result.x = pitch.fixAngle
+  result.y = yaw.fixAngle
+
+proc toAngles*[T](origin, target: GVec3[T]): GVec3[T] =
+  ## Gives Euler angles from origin to target
+  ##   pitch (x rotation)
+  ##   yaw (y rotation)
+  ##   roll (z rotation) - always 0 in vector case
+  ## All angles assume radians.
+  toAngles(target - origin)
+
+proc toAngles*[T](m: GMat4[T]): GVec3[T] =
+  ## Decomposes the matrix into Euler angles:
+  ##   pitch (x rotation)
+  ##   yaw (y rotation)
+  ##   roll (z rotation)
+  ## Assumes matrix has not been scaled.
+  ## All angles assume radians.
+  result.x = arcsin(m[2,1])
+  if result.x > PI/2:
+    # Degenerate case over north pole.
+    result.y = arctan2(m[0, 2], m[0, 0])
+  elif result.x < -PI/2:
+    # Degenerate case over south pole.
+    result.y = arctan2(m[0, 2], m[0, 0])
   else:
-    heading = arctan2(-m[8], m[0])
-    pitch = arctan2(-m[6], m[5])
-    roll = arcsin(m[4])
-  gvec3[T](heading, pitch, roll)
+    # Normal case.
+    result.y = -arctan2(m[2, 0], m[2, 2])
+    result.z = -arctan2(m[0, 1], m[1, 1])
+
+proc fromAngles*[T](a: GVec3[T]): GMat4[T] =
+  ## Takes a vector containing Euler angles and returns a matrix.
+  ## All angles assume radians.
+  rotateY(a.y) * rotateX(a.x) * rotateZ(a.z)
 
 proc frustum*[T](left, right, bottom, top, near, far: T): GMat4[T] =
   ## Create a frustum matrix.
@@ -1332,7 +1478,10 @@ proc ortho*[T](left, right, bottom, top, near, far: T): GMat4[T] =
   result[3, 2] = T(-(far + near) / fn)
   result[3, 3] = 1
 
-proc lookAt*[T](eye, center, up: GVec3[T]): GMat4[T] =
+proc lookAt*[T](eye, center, up: GVec3[T]): GMat4[T]
+  {.deprecated: "Wrong coordinate system. " &
+    "Use toAngles(eye, center).fromAngles() instead to get " &
+    "right-handed-z-forward coordinate system".} =
   ## Create a matrix that would convert eye pos to looking at center.
   let
     eyex = eye[0]
@@ -1413,7 +1562,10 @@ proc lookAt*[T](eye, center, up: GVec3[T]): GMat4[T] =
   result[3, 2] = -(z0 * eyex + z1 * eyey + z2 * eyez)
   result[3, 3] = 1
 
-proc lookAt*[T](eye, center: GVec3[T]): GMat4[T] =
+proc lookAt*[T](eye, center: GVec3[T]): GMat4[T]
+  {.deprecated: "Wrong coordinate system. " &
+    "Use toAngles(eye, center).fromAngles() instead to get " &
+    "right-handed-z-forward coordinate system".} =
   ## Look center from eye with default UP vector.
   lookAt(eye, center, gvec3(T(0), 0, 1))
 
@@ -1598,18 +1750,6 @@ proc rotate*[T](angle: T, axis: GVec3[T]): GMat4[T] =
   ## Return a rotation matrix with axis and angle.
   fromAxisAngle(axis, angle).mat4()
 
-proc rotateX*[T](angle: T): GMat4[T] =
-  ## Return a rotation matrix around X with angle.
-  fromAxisAngle(gvec3[T](1, 0, 0), angle).mat4()
-
-proc rotateY*[T](angle: T): GMat4[T] =
-  ## Return a rotation matrix around Y with angle.
-  fromAxisAngle(gvec3[T](0, 1, 0), angle).mat4()
-
-proc rotateZ*[T](angle: T): GMat4[T] =
-  ## Return a rotation matrix around Z with angle.
-  fromAxisAngle(gvec3[T](0, 0, 1), angle).mat4()
-
 when defined(release):
   {.pop.}
 {.pop.}
diff --git a/tests/test.nim b/tests/test.nim
index e7de428..ffba8db 100644
--- a/tests/test.nim
+++ b/tests/test.nim
@@ -53,18 +53,18 @@ block:
   doAssert quantize(1.23456789, 0.01) ~= 1.23
   doAssert quantize(-1.23456789, 0.01) ~= -1.23
 
-  doAssert fractional(0.0) ~= 0.0
-  doAssert fractional(3.14) ~= 0.14
-  doAssert fractional(-3.14) ~= 0.14
-  doAssert fractional(1.23456789) ~= 0.23456789
-  doAssert fractional(-1.23456789) ~= 0.23456789
+  doAssert fract(0.0) ~= 0.0
+  doAssert fract(3.14) ~= 0.14
+  doAssert fract(-3.14) ~= 0.14
+  doAssert fract(1.23456789) ~= 0.23456789
+  doAssert fract(-1.23456789) ~= 0.23456789
 
-  doAssert lerp(0.0, 1.0, 0.5) ~= 0.5
-  doAssert lerp(0.0, 10.0, 0.5) ~= 5.0
-  doAssert lerp(0.0, 100.0, 0.5) ~= 50.0
-  doAssert lerp(-1.0, 1.0, 0.25) ~= -0.5
-  doAssert lerp(-10.0, 10.0, 0.25) ~= -5.0
-  doAssert lerp(-100.0, 100.0, 0.25) ~= -50.0
+  doAssert mix(0.0, 1.0, 0.5) ~= 0.5
+  doAssert mix(0.0, 10.0, 0.5) ~= 5.0
+  doAssert mix(0.0, 100.0, 0.5) ~= 50.0
+  doAssert mix(-1.0, 1.0, 0.25) ~= -0.5
+  doAssert mix(-10.0, 10.0, 0.25) ~= -5.0
+  doAssert mix(-100.0, 100.0, 0.25) ~= -50.0
 
   doAssert mix(0.0, 1.0, 0.5) ~= 0.5
   doAssert mix(0.0, 10.0, 0.5) ~= 5.0
@@ -863,29 +863,7 @@ block:
     )
 
 block:
-  # test quat and matrix lookat
-  doAssert lookAt(vec3(1, 2, 3), vec3(0, 0, 0)).quat ~=
-    quat(
-      0.07232953608036041,
-      0.3063928484916687,
-      0.9237624406814575,
-      0.2180707305669785
-    )
-  doAssert lookAt(vec3(0, 0, 0), vec3(0, 0, 0)).quat ~= quat(0.0, 0.0, 0.0, 1.0)
-  doAssert lookAt(vec3(1, 0, 0), vec3(0, 0, 0)).quat ~= quat(0.5, 0.5, 0.5, 0.5)
-  doAssert lookAt(vec3(0, 1, 0), vec3(0, 0, 0)).quat ~=
-    quat(
-      0.0,
-      0.7071067690849304,
-      0.7071067690849304,
-      0.0
-    )
-  doAssert lookAt(vec3(0, 0, 1), vec3(0, 0, 0)).quat ~= quat(0.0, 0.0, 0.0, 1.0)
-
-  let
-    a = lookAt(vec3(1, 2, 3), vec3(0, 0, 0))
-    b = lookAt(dvec3(1, 2, 3), dvec3(0, 0, 0))
-
+  # test quat and matrix
   doAssert ortho[float32](-1, 1, 1, -1, -1000, 1000) ~= mat4(
     1.0, 0.0, 0.0, 0.0,
     0.0, -1.0, 0.0, 0.0,
@@ -1045,4 +1023,136 @@ block:
     b = a / 2
     when compiles(b = a div 2): doAssert false # type mismatch
 
+proc eq(a, b: Vec3): bool =
+  const epsilon = 0.001
+  return abs(angleBetween(a.x, b.x)) < epsilon and
+    abs(angleBetween(a.y, b.y)) < epsilon and
+    abs(angleBetween(a.z, b.z)) < epsilon
+
+const PI = PI.float32
+
+block:
+  # test Euler angles from a vector
+  doAssert vec3(0, 0, 0).toAngles.eq vec3(0f, 0f, 0f)
+  doAssert vec3(0, 0, 1).toAngles.eq vec3(0f, 0f, 0f) # forward
+  doAssert vec3(0, 0, -1).toAngles.eq vec3(0f, PI, 0f) # back
+  doAssert vec3(-1, 0, 0).toAngles.eq vec3(0f, PI/2, 0f) # right
+  doAssert vec3(1, 0, 0).toAngles.eq vec3(0f, -PI/2, 0f) # left
+  doAssert vec3(0, 1, 0).toAngles.eq vec3(PI/2, 0f, 0f) # up
+  doAssert vec3(0, -1, 0).toAngles.eq vec3(-PI/2, 0f, 0f) # down
+
+block:
+  # test Euler angles from a matrix
+  doAssert translate(vec3(0, 0, 0)).toAngles.eq vec3(0f, 0f, 0f)
+  doAssert rotateX(0f).toAngles.eq vec3(0f, 0f, 0f) # forward
+  doAssert rotateY(PI).toAngles.eq vec3(0f, -PI, 0f) # back
+  doAssert rotateY(PI/2).toAngles.eq vec3(0f, PI/2, 0f) # back
+  doAssert rotateY(-PI/2).toAngles.eq vec3(0f, -PI/2, 0f) # back
+  doAssert rotateX(PI/2).toAngles.eq vec3(PI/2, 0f, 0f) # up
+  doAssert rotateX(-PI/2).toAngles.eq vec3(-PI/2, 0f, 0f) # down
+  doAssert rotateZ(PI/2).toAngles.eq vec3(0f, 0f, PI/2) # tilt right
+  doAssert rotateZ(-PI/2).toAngles.eq vec3(0f, 0f, -PI/2) # tilt left
+
+  doAssert mat4().toAngles.eq vec3(0, 0, 0)
+
+  doAssert rotateX(10.toRadians()).toAngles.eq vec3(10.toRadians(), 0, 0)
+  doAssert rotateY(10.toRadians()).toAngles.eq vec3(0, 10.toRadians(), 0)
+  doAssert rotateZ(10.toRadians()).toAngles.eq vec3(0, 0, 10.toRadians())
+  doAssert rotateX(89.toRadians()).toAngles.eq vec3(89.toRadians(), 0, 0)
+  doAssert rotateY(89.toRadians()).toAngles.eq vec3(0, 89.toRadians(), 0)
+  doAssert rotateZ(89.toRadians()).toAngles.eq vec3(0, 0, 89.toRadians())
+  doAssert rotateX(90.toRadians()).toAngles.eq vec3(90.toRadians(), 0, 0)
+  doAssert rotateY(90.toRadians()).toAngles.eq vec3(0, 90.toRadians(), 0)
+  doAssert rotateZ(90.toRadians()).toAngles.eq vec3(0, 0, 90.toRadians())
+  doAssert rotateX(90.toRadians()).toAngles.eq vec3(90.toRadians(), 0, 0)
+  doAssert rotateY(90.toRadians()).toAngles.eq vec3(0, 90.toRadians(), 0)
+  doAssert rotateZ(-90.toRadians()).toAngles.eq vec3(0, 0, -90.toRadians())
+  doAssert rotateY(180.toRadians()).toAngles.eq vec3(0, -180.toRadians(), 0)
+  doAssert rotateZ(180.toRadians()).toAngles.eq vec3(0, 0, 180.toRadians())
+  doAssert rotateY(-180.toRadians()).toAngles.eq vec3(0, 180.toRadians(), 0)
+  doAssert rotateZ(-180.toRadians()).toAngles.eq vec3(0, 0, 180.toRadians())
+
+block:
+  # Euler angles fuzzing tests.
+
+  # Test fromAngles with and without roll have same forward
+  for i in 0 .. 1000:
+    let
+      xr = rand(-89.9f .. 89.9f).toRadians
+      yr = rand(-180 .. 180).toRadians
+      zr = rand(-180 .. 180).toRadians
+      a = vec3(xr, yr, zr)
+      b = vec3(xr, yr, 0f)
+      ma = fromAngles(a)
+      mb = fromAngles(b)
+
+    doAssert ma.forward() ~= mb.forward()
+
+  # Test forward/back, right/left, up/down combos
+  for i in 0 .. 1000:
+    let
+      xr = rand(-89.9f .. 89.9f).toRadians
+      yr = rand(-180 .. 180).toRadians
+      zr = rand(-180 .. 180).toRadians
+      b = vec3(xr, yr, zr)
+      m = fromAngles(b)
+
+    doAssert m.forward() ~= m * vec3(0, 0, 1)
+    doAssert m.back() ~= m * vec3(0, 0, -1)
+
+    doAssert m.right() ~= m * vec3(-1, 0, 0)
+    doAssert m.left() ~= m * vec3(1, 0, 0)
+
+    doAssert m.up() ~= m * vec3(0, 1, 0)
+    doAssert m.down() ~= m * vec3(0, -1, 0)
+
+  # Test non-polar and non-rotated cases
+  for i in 0 .. 1000:
+    let
+      xr = rand(-89.9f .. 89.9f).toRadians
+      yr = rand(-180 .. 180).toRadians
+      zr = 0f
+      b = vec3(xr, yr, zr)
+      m = fromAngles(b)
+      a = m.toAngles()
+    doAssert a.eq(b)
+
+  # Test non-polar cases
+  for i in 0 .. 1000:
+    let
+      xr = rand(-89.9f .. 89.9f).toRadians
+      yr = rand(-180 .. 180).toRadians
+      zr = rand(-180 .. 180).toRadians
+      b = vec3(xr, yr, zr)
+      m = fromAngles(b)
+      a = m.toAngles()
+    doAssert a.eq(b)
+
+  # Test polar and non-rotated cases
+  for i in 0 .. 1000:
+    let
+      xr = sample([-90, 90]).toRadians
+      yr = rand(-180 .. 180).toRadians
+      zr = 0f
+      b = vec3(xr, yr, zr)
+      m = fromAngles(b)
+      a = m.toAngles()
+    doAssert a.eq(b)
+
+  # Test polar and crazy rotated cases
+  for i in 0 .. 1000:
+    let
+      xr = sample([-90, 90]).toRadians
+      yr = rand(-180 .. 180).toRadians
+      zr = rand(-180 .. 180).toRadians
+      b = vec3(xr, yr, zr)
+      m = fromAngles(b)
+      a = m.toAngles()
+
+    doAssert abs(angleBetween(a.x, b.x)) < 0.001
+    if xr > 0:
+      doAssert abs(angleBetween(a.y, b.y + b.z)) < 0.001
+    else:
+      doAssert abs(angleBetween(a.y, b.y - b.z)) < 0.001
+
 echo "test finished successfully"
diff --git a/vmath.nimble b/vmath.nimble
index c414c6c..0198ad3 100644
--- a/vmath.nimble
+++ b/vmath.nimble
@@ -1,4 +1,4 @@
-version       = "1.2.0"
+version       = "2.0.0"
 author        = "Andre von Houck"
 description   = "Your single stop for vector math routines for 2d and 3d graphics."
 license       = "MIT"