diff --git a/.changeset/tidy-otters-listen.md b/.changeset/tidy-otters-listen.md
new file mode 100644
index 0000000..4600df4
--- /dev/null
+++ b/.changeset/tidy-otters-listen.md
@@ -0,0 +1,5 @@
+---
+"math": minor
+---
+
+accepts readonly inputs across core, shapes, color, geometry, time, ik and random operations
diff --git a/src/color/color.ts b/src/color/color.ts
index 732e982..589b0b6 100644
--- a/src/color/color.ts
+++ b/src/color/color.ts
@@ -5,6 +5,9 @@ export * from './parse';
 /** A linear-sRGB color: [r, g, b] floats in [0, 1]. */
 export type Color = [r: number, g: number, b: number];
 
+/** A read-only linear sRGB color */
+export type RColor = Readonly<Color>;
+
 /** Accepted input types for creating or parsing a Color. */
 export type ColorInput =
     | string // '#f00', '#ff0000', 'red', 'rgb(255,0,0)', 'hsl(0,100%,50%)'
@@ -22,12 +25,12 @@ export function fromValues(r: number, g: number, b: number): Color {
 }
 
 /** Create a new Color that is a copy of `c`. */
-export function clone(c: Color): Color {
+export function clone(c: RColor): Color {
     return [c[0], c[1], c[2]];
 }
 
 /** Copy the values from `src` into `out`. Returns `out`. */
-export function copy(out: Color, src: Color): Color {
+export function copy(out: Color, src: RColor): Color {
     out[0] = src[0];
     out[1] = src[1];
     out[2] = src[2];
@@ -54,7 +57,7 @@ export function setScalar(out: Color, s: number): Color {
  * Set `out` from an sRGB gamma-encoded [r, g, b] array with values in [0, 1].
  * Converts from sRGB gamma space to linear. Returns `out`.
  */
-export function setFromSRGB(out: Color, srgb: [number, number, number]): Color {
+export function setFromSRGB(out: Color, srgb: readonly [number, number, number]): Color {
     out[0] = srgbToLinear(srgb[0]);
     out[1] = srgbToLinear(srgb[1]);
     out[2] = srgbToLinear(srgb[2]);
@@ -62,12 +65,12 @@ export function setFromSRGB(out: Color, srgb: [number, number, number]): Color {
 }
 
 /** Create a new Color from an sRGB gamma-encoded [r, g, b] array with values in [0, 1]. */
-export function fromSRGB(srgb: [number, number, number]): Color {
+export function fromSRGB(srgb: readonly [number, number, number]): Color {
     return setFromSRGB(create(), srgb);
 }
 
 /** Write the sRGB gamma-encoded [r, g, b] of a linear Color into `out` (values [0, 1]). */
-export function toSRGB(out: [number, number, number], c: Color): [number, number, number] {
+export function toSRGB(out: [number, number, number], c: RColor): [number, number, number] {
     out[0] = linearToSrgb(c[0]);
     out[1] = linearToSrgb(c[1]);
     out[2] = linearToSrgb(c[2]);
@@ -75,22 +78,22 @@ export function toSRGB(out: [number, number, number], c: Color): [number, number
 }
 
 /** Create a CSS `rgb(...)` string in sRGB gamma space (for HTML/canvas use). */
-export function toCSS(c: Color): string {
+export function toCSS(c: RColor): string {
     return `rgb(${to255(c[0])}, ${to255(c[1])}, ${to255(c[2])})`;
 }
 
 /** Convert to a 0xRRGGBB integer in sRGB gamma space. */
-export function toHex(c: Color): number {
+export function toHex(c: RColor): number {
     return (to255(c[0]) << 16) | (to255(c[1]) << 8) | to255(c[2]);
 }
 
 /** Convert to a 6-digit sRGB hex string without a leading '#', e.g. 'ff8800'. */
-export function toHexString(c: Color): string {
+export function toHexString(c: RColor): string {
     return toHex(c).toString(16).padStart(6, '0');
 }
 
 /** Add `a + b` component-wise into `out`. Returns `out`. */
-export function add(out: Color, a: Color, b: Color): Color {
+export function add(out: Color, a: RColor, b: RColor): Color {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -98,7 +101,7 @@ export function add(out: Color, a: Color, b: Color): Color {
 }
 
 /** Add scalar `s` to each channel of `a` into `out`. Returns `out`. */
-export function addScalar(out: Color, a: Color, s: number): Color {
+export function addScalar(out: Color, a: RColor, s: number): Color {
     out[0] = a[0] + s;
     out[1] = a[1] + s;
     out[2] = a[2] + s;
@@ -106,7 +109,7 @@ export function addScalar(out: Color, a: Color, s: number): Color {
 }
 
 /** Subtract `a - b` component-wise into `out`. Returns `out`. */
-export function sub(out: Color, a: Color, b: Color): Color {
+export function sub(out: Color, a: RColor, b: RColor): Color {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     out[2] = a[2] - b[2];
@@ -114,7 +117,7 @@ export function sub(out: Color, a: Color, b: Color): Color {
 }
 
 /** Multiply `a * b` component-wise into `out` (tinting). Returns `out`. */
-export function multiply(out: Color, a: Color, b: Color): Color {
+export function multiply(out: Color, a: RColor, b: RColor): Color {
     out[0] = a[0] * b[0];
     out[1] = a[1] * b[1];
     out[2] = a[2] * b[2];
@@ -122,7 +125,7 @@ export function multiply(out: Color, a: Color, b: Color): Color {
 }
 
 /** Scale each channel of `a` by `s` into `out` (brightness). Returns `out`. */
-export function multiplyScalar(out: Color, a: Color, s: number): Color {
+export function multiplyScalar(out: Color, a: RColor, s: number): Color {
     out[0] = a[0] * s;
     out[1] = a[1] * s;
     out[2] = a[2] * s;
@@ -130,7 +133,7 @@ export function multiplyScalar(out: Color, a: Color, s: number): Color {
 }
 
 /** Linearly interpolate from `a` to `b` by `t` into `out` (physically-correct blend). Returns `out`. */
-export function lerp(out: Color, a: Color, b: Color, t: number): Color {
+export function lerp(out: Color, a: RColor, b: RColor, t: number): Color {
     out[0] = a[0] + (b[0] - a[0]) * t;
     out[1] = a[1] + (b[1] - a[1]) * t;
     out[2] = a[2] + (b[2] - a[2]) * t;
@@ -138,7 +141,7 @@ export function lerp(out: Color, a: Color, b: Color, t: number): Color {
 }
 
 /** Clamp each channel of `c` to [0, 1] into `out`. Returns `out`. */
-export function clamp(out: Color, c: Color): Color {
+export function clamp(out: Color, c: RColor): Color {
     out[0] = clamp01(c[0]);
     out[1] = clamp01(c[1]);
     out[2] = clamp01(c[2]);
@@ -146,12 +149,12 @@ export function clamp(out: Color, c: Color): Color {
 }
 
 /** Whether `a` and `b` are equal, within an optional per-channel `epsilon` (default exact). */
-export function equals(a: Color, b: Color, epsilon = 0): boolean {
+export function equals(a: RColor, b: RColor, epsilon = 0): boolean {
     return Math.abs(a[0] - b[0]) <= epsilon && Math.abs(a[1] - b[1]) <= epsilon && Math.abs(a[2] - b[2]) <= epsilon;
 }
 
 /** Relative luminance in [0, 1] (Rec. 709 weights, on linear light). */
-export function luminance(c: Color): number {
+export function luminance(c: RColor): number {
     return 0.2126 * c[0] + 0.7152 * c[1] + 0.0722 * c[2];
 }
 
diff --git a/src/color/colorspace.ts b/src/color/colorspace.ts
index 617c7e4..8350400 100644
--- a/src/color/colorspace.ts
+++ b/src/color/colorspace.ts
@@ -1,4 +1,4 @@
-import type { Color } from './color';
+import type { Color, RColor } from './color';
 
 // Color-space conversions (pure functions — no global working-space state).
 //
@@ -21,7 +21,7 @@ export function linearToSrgb(c: number): number {
  * Convert a linear-sRGB Color to linear Display-P3 primaries, into `out`. Returns `out`.
  * (Both spaces share the sRGB transfer curve; this changes only the primaries.)
  */
-export function linearSrgbToLinearDisplayP3(out: Color, c: Color): Color {
+export function linearSrgbToLinearDisplayP3(out: Color, c: RColor): Color {
     const r = c[0];
     const g = c[1];
     const b = c[2];
@@ -35,7 +35,7 @@ export function linearSrgbToLinearDisplayP3(out: Color, c: Color): Color {
  * Convert a linear Display-P3 Color to linear-sRGB primaries, into `out`. Returns `out`.
  * Colors outside the sRGB gamut yield channels outside [0, 1] — clamp if needed.
  */
-export function linearDisplayP3ToLinearSrgb(out: Color, c: Color): Color {
+export function linearDisplayP3ToLinearSrgb(out: Color, c: RColor): Color {
     const r = c[0];
     const g = c[1];
     const b = c[2];
diff --git a/src/color/hsl.ts b/src/color/hsl.ts
index 505cbd7..ea74de4 100644
--- a/src/color/hsl.ts
+++ b/src/color/hsl.ts
@@ -1,9 +1,12 @@
-import type { Color } from './color';
+import type { Color, RColor } from './color';
 import { linearToSrgb, srgbToLinear } from './colorspace';
 
 /** A hue-saturation-lightness color: [h, s, l], all in [0, 1] (hue wraps). */
 export type HSL = [hue: number, saturation: number, lightness: number];
 
+/** A read-only HSL color */
+export type RHSL = Readonly<HSL>;
+
 /** Create a new HSL initialized to [0, 0, 0] (black). */
 export function create(): HSL {
     return [0, 0, 0];
@@ -15,12 +18,12 @@ export function fromValues(h: number, s: number, l: number): HSL {
 }
 
 /** Create a new HSL that is a copy of `a`. */
-export function clone(a: HSL): HSL {
+export function clone(a: RHSL): HSL {
     return [a[0], a[1], a[2]];
 }
 
 /** Copy the values from `src` into `out`. Returns `out`. */
-export function copy(out: HSL, src: HSL): HSL {
+export function copy(out: HSL, src: RHSL): HSL {
     out[0] = src[0];
     out[1] = src[1];
     out[2] = src[2];
@@ -36,7 +39,7 @@ export function set(out: HSL, h: number, s: number, l: number): HSL {
 }
 
 /** Write the HSL of a linear Color into `out`. Returns `out`. */
-export function fromColor(out: HSL, c: Color): HSL {
+export function fromColor(out: HSL, c: RColor): HSL {
     // linear -> sRGB gamma; HSL is defined on gamma-encoded sRGB
     const r = linearToSrgb(c[0]);
     const g = linearToSrgb(c[1]);
@@ -64,7 +67,7 @@ export function fromColor(out: HSL, c: Color): HSL {
 }
 
 /** Write the linear Color of an HSL into `out`. Returns `out`. */
-export function toColor(out: Color, a: HSL): Color {
+export function toColor(out: Color, a: RHSL): Color {
     const h = a[0];
     const s = a[1];
     const l = a[2];
@@ -90,7 +93,7 @@ export function toColor(out: Color, a: HSL): Color {
  * the hue wheel (so e.g. 350°→10° passes through 0°, not all the way back).
  * Returns `out`.
  */
-export function lerp(out: HSL, a: HSL, b: HSL, t: number): HSL {
+export function lerp(out: HSL, a: RHSL, b: RHSL, t: number): HSL {
     let dh = b[0] - a[0];
     if (dh > 0.5) dh -= 1;
     else if (dh < -0.5) dh += 1;
@@ -108,7 +111,7 @@ export function lerp(out: HSL, a: HSL, b: HSL, t: number): HSL {
  * Offset `a` by (dh, ds, dl) into `out`: hue wraps into [0, 1), saturation and
  * lightness are clamped to [0, 1]. Returns `out`.
  */
-export function offset(out: HSL, a: HSL, dh: number, ds: number, dl: number): HSL {
+export function offset(out: HSL, a: RHSL, dh: number, ds: number, dl: number): HSL {
     let h = a[0] + dh;
     h -= Math.floor(h);
     out[0] = h;
diff --git a/src/color/index.ts b/src/color/index.ts
index 1ee674b..904d0ff 100644
--- a/src/color/index.ts
+++ b/src/color/index.ts
@@ -2,5 +2,5 @@ export * as color from './color';
 export * as colorspace from './colorspace';
 export * as hsl from './hsl';
 
-export type { Color, ColorInput } from './color';
-export type { HSL } from './hsl';
+export type { Color, ColorInput, RColor } from './color';
+export type { HSL, RHSL } from './hsl';
diff --git a/src/core/euler.ts b/src/core/euler.ts
index 7e8594f..27886dd 100644
--- a/src/core/euler.ts
+++ b/src/core/euler.ts
@@ -1,5 +1,5 @@
-import type { Mat4 } from './mat4';
-import type { Quat } from './quat';
+import type { RMat4 } from './mat4';
+import type { RQuat } from './quat';
 import * as quat from './quat';
 import { clamp, EPSILON } from './scalar';
 
@@ -11,6 +11,9 @@ export type EulerOrder = 'xyz' | 'xzy' | 'yxz' | 'yzx' | 'zxy' | 'zyx';
 /** A Euler in 3D space, with an optional order (default is 'xyz') */
 export type Euler = [x: number, y: number, z: number, order?: EulerOrder];
 
+/** A read-only set of Euler angles */
+export type REuler = Readonly<Euler>;
+
 /**
  * Creates a new Euler with default values (0, 0, 0, 'xyz').
  */
@@ -71,7 +74,7 @@ export function fromDegrees(out: Euler, x: number, y: number, z: number, order:
  * @param order The order of the Euler angles.
  * @returns The output Euler.
  */
-export function fromRotationMat4(out: Euler, rotationMatrix: Mat4, order: EulerOrder = out[3] || 'xyz'): Euler {
+export function fromRotationMat4(out: Euler, rotationMatrix: RMat4, order: EulerOrder = out[3] || 'xyz'): Euler {
     return fromRotationMatrixValues(
         out,
         rotationMatrix[0],
@@ -199,7 +202,7 @@ function fromRotationMatrixValues(
  * @param b The second euler.
  * @returns True if the euler angles are equal, false otherwise.
  */
-export function exactEquals(a: Euler, b: Euler): boolean {
+export function exactEquals(a: REuler, b: REuler): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2] && a[3] === b[3];
 }
 
@@ -210,7 +213,7 @@ export function exactEquals(a: Euler, b: Euler): boolean {
  * @param b The second euler.
  * @returns True if the euler angles are equal, false otherwise.
  */
-export function equals(a: Euler, b: Euler): boolean {
+export function equals(a: REuler, b: REuler): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -232,7 +235,7 @@ export function equals(a: Euler, b: Euler): boolean {
  * @param order The order of the Euler.
  * @returns The output Euler
  */
-export function fromQuat(out: Euler, q: Quat, order: EulerOrder): Euler {
+export function fromQuat(out: Euler, q: RQuat, order: EulerOrder): Euler {
     // compute the rotation matrix elements directly from the quaternion
     const x = q[0];
     const y = q[1];
@@ -276,7 +279,7 @@ const _reorderQuaternion = /*@__PURE__*/ quat.create();
  * @param order The order of the Euler.
  * @returns The output Euler.
  */
-export function reorder(out: Euler, a: Euler, order: EulerOrder): Euler {
+export function reorder(out: Euler, a: REuler, order: EulerOrder): Euler {
     quat.fromEuler(_reorderQuaternion, a);
     fromQuat(out, _reorderQuaternion, order);
     return out;
diff --git a/src/core/index.ts b/src/core/index.ts
index b56a212..34173a2 100644
--- a/src/core/index.ts
+++ b/src/core/index.ts
@@ -1,39 +1,40 @@
+export type { DeepReadonly } from './readonly';
 export * from './scalar';
 export * from './angle';
 export type { MutableArrayLike } from './arrays';
 
 export * as vec2 from './vec2';
-export type { Vec2 } from './vec2';
+export type { RVec2, Vec2 } from './vec2';
 
 export * as vec3 from './vec3';
-export type { Vec3 } from './vec3';
+export type { RVec3, Vec3 } from './vec3';
 
 export * as vec4 from './vec4';
-export type { Vec4 } from './vec4';
+export type { RVec4, Vec4 } from './vec4';
 
 export * as euler from './euler';
-export type { Euler, EulerOrder } from './euler';
+export type { Euler, EulerOrder, REuler } from './euler';
 
 export * as quat from './quat';
-export type { Quat } from './quat';
+export type { Quat, RQuat } from './quat';
 
 export * as quat2 from './quat2';
-export type { Quat2 } from './quat2';
+export type { Quat2, RQuat2 } from './quat2';
 
 export * as mat2 from './mat2';
-export type { Mat2 } from './mat2';
+export type { Mat2, RMat2 } from './mat2';
 
 export * as mat2d from './mat2d';
-export type { Mat2d } from './mat2d';
+export type { Mat2d, RMat2d } from './mat2d';
 
 export * as mat3 from './mat3';
-export type { Mat3 } from './mat3';
+export type { Mat3, RMat3 } from './mat3';
 
 export * as mat4 from './mat4';
-export type { Mat4 } from './mat4';
+export type { Mat4, RMat4 } from './mat4';
 
 export * as spherical from './spherical';
-export type { Spherical } from './spherical';
+export type { RSpherical, Spherical } from './spherical';
 
 export * as polar from './polar';
-export type { Polar } from './polar';
+export type { Polar, RPolar } from './polar';
diff --git a/src/core/mat2.ts b/src/core/mat2.ts
index d4e26a8..34e2119 100644
--- a/src/core/mat2.ts
+++ b/src/core/mat2.ts
@@ -1,9 +1,12 @@
 import { EPSILON } from './scalar';
-import type { Vec2 } from './vec2';
+import type { RVec2 } from './vec2';
 
 /** A 2x2 matrix */
 export type Mat2 = [e1: number, e2: number, e3: number, e4: number];
 
+/** A read-only 2x2 matrix */
+export type RMat2 = Readonly<Mat2>;
+
 /**
  * Creates a new identity mat2
  *
@@ -19,7 +22,7 @@ export function create(): Mat2 {
  * @param a matrix to clone
  * @returns a new 2x2 matrix
  */
-export function clone(a: Mat2): Mat2 {
+export function clone(a: RMat2): Mat2 {
     return [a[0], a[1], a[2], a[3]];
 }
 
@@ -30,7 +33,7 @@ export function clone(a: Mat2): Mat2 {
  * @param a the source matrix
  * @returns out
  */
-export function copy(out: Mat2, a: Mat2): Mat2 {
+export function copy(out: Mat2, a: RMat2): Mat2 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -90,7 +93,7 @@ export function set(out: Mat2, m00: number, m01: number, m10: number, m11: numbe
  * @param a the source matrix
  * @returns out
  */
-export function transpose(out: Mat2, a: Mat2): Mat2 {
+export function transpose(out: Mat2, a: RMat2): Mat2 {
     // If we are transposing ourselves we can skip a few steps but have to cache
     // some values
     if (out === a) {
@@ -114,7 +117,7 @@ export function transpose(out: Mat2, a: Mat2): Mat2 {
  * @param a the source matrix
  * @returns out, or null if source matrix is not invertible
  */
-export function invert(out: Mat2, a: Mat2): Mat2 | null {
+export function invert(out: Mat2, a: RMat2): Mat2 | null {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -143,7 +146,7 @@ export function invert(out: Mat2, a: Mat2): Mat2 | null {
  * @param a the source matrix
  * @returns out
  */
-export function adjoint(out: Mat2, a: Mat2): Mat2 {
+export function adjoint(out: Mat2, a: RMat2): Mat2 {
     // Caching this value is necessary if out == a
     const a0 = a[0];
     out[0] = a[3];
@@ -160,7 +163,7 @@ export function adjoint(out: Mat2, a: Mat2): Mat2 {
  * @param a the source matrix
  * @returns determinant of a
  */
-export function determinant(a: Mat2): number {
+export function determinant(a: RMat2): number {
     return a[0] * a[3] - a[2] * a[1];
 }
 
@@ -172,7 +175,7 @@ export function determinant(a: Mat2): number {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Mat2, a: Mat2, b: Mat2): Mat2 {
+export function multiply(out: Mat2, a: RMat2, b: RMat2): Mat2 {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -196,7 +199,7 @@ export function multiply(out: Mat2, a: Mat2, b: Mat2): Mat2 {
  * @param rad the angle to rotate the matrix by
  * @returns out
  */
-export function rotate(out: Mat2, a: Mat2, rad: number): Mat2 {
+export function rotate(out: Mat2, a: RMat2, rad: number): Mat2 {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -218,7 +221,7 @@ export function rotate(out: Mat2, a: Mat2, rad: number): Mat2 {
  * @param v the vec2 to scale the matrix by
  * @returns out
  **/
-export function scale(out: Mat2, a: Mat2, v: Vec2): Mat2 {
+export function scale(out: Mat2, a: RMat2, v: RVec2): Mat2 {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -264,7 +267,7 @@ export function fromRotation(out: Mat2, rad: number): Mat2 {
  * @param v Scaling vector
  * @returns out
  */
-export function fromScaling(out: Mat2, v: Vec2): Mat2 {
+export function fromScaling(out: Mat2, v: RVec2): Mat2 {
     out[0] = v[0];
     out[1] = 0;
     out[2] = 0;
@@ -278,7 +281,7 @@ export function fromScaling(out: Mat2, v: Vec2): Mat2 {
  * @param a matrix to represent as a string
  * @returns string representation of the matrix
  */
-export function str(a: Mat2): string {
+export function str(a: RMat2): string {
     return `mat2(${a[0]}, ${a[1]}, ${a[2]}, ${a[3]})`;
 }
 
@@ -288,7 +291,7 @@ export function str(a: Mat2): string {
  * @param a the matrix to calculate Frobenius norm of
  * @returns Frobenius norm
  */
-export function frob(a: Mat2): number {
+export function frob(a: RMat2): number {
     return Math.sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3]);
 }
 
@@ -300,7 +303,7 @@ export function frob(a: Mat2): number {
  * @param a the input matrix to factorize
  */
 
-export function LDU(L: Mat2, D: Mat2, U: Mat2, a: Mat2): [Mat2, Mat2, Mat2] {
+export function LDU(L: Mat2, D: Mat2, U: Mat2, a: RMat2): [Mat2, Mat2, Mat2] {
     L[2] = a[2] / a[0];
     U[0] = a[0];
     U[1] = a[1];
@@ -316,7 +319,7 @@ export function LDU(L: Mat2, D: Mat2, U: Mat2, a: Mat2): [Mat2, Mat2, Mat2] {
  * @param b the second operand
  * @returns out
  */
-export function add(out: Mat2, a: Mat2, b: Mat2): Mat2 {
+export function add(out: Mat2, a: RMat2, b: RMat2): Mat2 {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -332,7 +335,7 @@ export function add(out: Mat2, a: Mat2, b: Mat2): Mat2 {
  * @param b the second operand
  * @returns out
  */
-export function subtract(out: Mat2, a: Mat2, b: Mat2): Mat2 {
+export function subtract(out: Mat2, a: RMat2, b: RMat2): Mat2 {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     out[2] = a[2] - b[2];
@@ -347,7 +350,7 @@ export function subtract(out: Mat2, a: Mat2, b: Mat2): Mat2 {
  * @param b The second matrix.
  * @returns True if the matrices are equal, false otherwise.
  */
-export function exactEquals(a: Mat2, b: Mat2): boolean {
+export function exactEquals(a: RMat2, b: RMat2): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2] && a[3] === b[3];
 }
 
@@ -358,7 +361,7 @@ export function exactEquals(a: Mat2, b: Mat2): boolean {
  * @param b The second matrix.
  * @returns True if the matrices are equal, false otherwise.
  */
-export function equals(a: Mat2, b: Mat2): boolean {
+export function equals(a: RMat2, b: RMat2): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -383,7 +386,7 @@ export function equals(a: Mat2, b: Mat2): boolean {
  * @param b amount to scale the matrix's elements by
  * @returns out
  */
-export function multiplyScalar(out: Mat2, a: Mat2, b: number): Mat2 {
+export function multiplyScalar(out: Mat2, a: RMat2, b: number): Mat2 {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     out[2] = a[2] * b;
@@ -400,7 +403,7 @@ export function multiplyScalar(out: Mat2, a: Mat2, b: number): Mat2 {
  * @param scale the amount to scale b's elements by before adding
  * @returns out
  */
-export function multiplyScalarAndAdd(out: Mat2, a: Mat2, b: Mat2, scale: number): Mat2 {
+export function multiplyScalarAndAdd(out: Mat2, a: RMat2, b: RMat2, scale: number): Mat2 {
     out[0] = a[0] + b[0] * scale;
     out[1] = a[1] + b[1] * scale;
     out[2] = a[2] + b[2] * scale;
diff --git a/src/core/mat2d.ts b/src/core/mat2d.ts
index bcc1a67..2d1cacf 100644
--- a/src/core/mat2d.ts
+++ b/src/core/mat2d.ts
@@ -1,9 +1,12 @@
 import { EPSILON } from './scalar';
-import type { Vec2 } from './vec2';
+import type { RVec2 } from './vec2';
 
 /** A 2D affine transform matrix */
 export type Mat2d = [e1: number, e2: number, e3: number, e4: number, e5: number, e6: number];
 
+/** A read-only 2x3 matrix */
+export type RMat2d = Readonly<Mat2d>;
+
 /**
  * Creates a new identity mat2d
  *
@@ -19,7 +22,7 @@ export function create(): Mat2d {
  * @param a matrix to clone
  * @returns a new 2x3 matrix
  */
-export function clone(a: Mat2d): Mat2d {
+export function clone(a: RMat2d): Mat2d {
     return [a[0], a[1], a[2], a[3], a[4], a[5]];
 }
 
@@ -30,7 +33,7 @@ export function clone(a: Mat2d): Mat2d {
  * @param a the source matrix
  * @returns out
  */
-export function copy(out: Mat2d, a: Mat2d): Mat2d {
+export function copy(out: Mat2d, a: RMat2d): Mat2d {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -100,7 +103,7 @@ export function set(out: Mat2d, a: number, b: number, c: number, d: number, tx:
  * @param a the source matrix
  * @returns out, or null if source matrix is not invertible
  */
-export function invert(out: Mat2d, a: Mat2d): Mat2d | null {
+export function invert(out: Mat2d, a: RMat2d): Mat2d | null {
     const aa = a[0];
     const ab = a[1];
     const ac = a[2];
@@ -129,7 +132,7 @@ export function invert(out: Mat2d, a: Mat2d): Mat2d | null {
  * @param a the source matrix
  * @returns determinant of a
  */
-export function determinant(a: Mat2d): number {
+export function determinant(a: RMat2d): number {
     return a[0] * a[3] - a[1] * a[2];
 }
 
@@ -141,7 +144,7 @@ export function determinant(a: Mat2d): number {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Mat2d, a: Mat2d, b: Mat2d): Mat2d {
+export function multiply(out: Mat2d, a: RMat2d, b: RMat2d): Mat2d {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -171,7 +174,7 @@ export function multiply(out: Mat2d, a: Mat2d, b: Mat2d): Mat2d {
  * @param rad the angle to rotate the matrix by
  * @returns out
  */
-export function rotate(out: Mat2d, a: Mat2d, rad: number): Mat2d {
+export function rotate(out: Mat2d, a: RMat2d, rad: number): Mat2d {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -197,7 +200,7 @@ export function rotate(out: Mat2d, a: Mat2d, rad: number): Mat2d {
  * @param v the vec2 to scale the matrix by
  * @returns out
  **/
-export function scale(out: Mat2d, a: Mat2d, v: Vec2): Mat2d {
+export function scale(out: Mat2d, a: RMat2d, v: RVec2): Mat2d {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -223,7 +226,7 @@ export function scale(out: Mat2d, a: Mat2d, v: Vec2): Mat2d {
  * @param v the vec2 to translate the matrix by
  * @returns out
  **/
-export function translate(out: Mat2d, a: Mat2d, v: Vec2): Mat2d {
+export function translate(out: Mat2d, a: RMat2d, v: RVec2): Mat2d {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -275,7 +278,7 @@ export function fromRotation(out: Mat2d, rad: number): Mat2d {
  * @param v Scaling vector
  * @returns out
  */
-export function fromScaling(out: Mat2d, v: Vec2): Mat2d {
+export function fromScaling(out: Mat2d, v: RVec2): Mat2d {
     out[0] = v[0];
     out[1] = 0;
     out[2] = 0;
@@ -296,7 +299,7 @@ export function fromScaling(out: Mat2d, v: Vec2): Mat2d {
  * @param v Translation vector
  * @returns out
  */
-export function fromTranslation(out: Mat2d, v: Vec2): Mat2d {
+export function fromTranslation(out: Mat2d, v: RVec2): Mat2d {
     out[0] = 1;
     out[1] = 0;
     out[2] = 0;
@@ -312,7 +315,7 @@ export function fromTranslation(out: Mat2d, v: Vec2): Mat2d {
  * @param a matrix to represent as a string
  * @returns string representation of the matrix
  */
-export function str(a: Mat2d): string {
+export function str(a: RMat2d): string {
     return `mat2d(${a[0]}, ${a[1]}, ${a[2]}, ${a[3]}, ${a[4]}, ${a[5]})`;
 }
 
@@ -322,7 +325,7 @@ export function str(a: Mat2d): string {
  * @param a the matrix to calculate Frobenius norm of
  * @returns Frobenius norm
  */
-export function frob(a: Mat2d): number {
+export function frob(a: RMat2d): number {
     return Math.sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3] + a[4] * a[4] + a[5] * a[5] + 1);
 }
 
@@ -334,7 +337,7 @@ export function frob(a: Mat2d): number {
  * @param b the second operand
  * @returns out
  */
-export function add(out: Mat2d, a: Mat2d, b: Mat2d): Mat2d {
+export function add(out: Mat2d, a: RMat2d, b: RMat2d): Mat2d {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -352,7 +355,7 @@ export function add(out: Mat2d, a: Mat2d, b: Mat2d): Mat2d {
  * @param b the second operand
  * @returns out
  */
-export function subtract(out: Mat2d, a: Mat2d, b: Mat2d): Mat2d {
+export function subtract(out: Mat2d, a: RMat2d, b: RMat2d): Mat2d {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     out[2] = a[2] - b[2];
@@ -370,7 +373,7 @@ export function subtract(out: Mat2d, a: Mat2d, b: Mat2d): Mat2d {
  * @param b amount to scale the matrix's elements by
  * @returns out
  */
-export function multiplyScalar(out: Mat2d, a: Mat2d, b: number): Mat2d {
+export function multiplyScalar(out: Mat2d, a: RMat2d, b: number): Mat2d {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     out[2] = a[2] * b;
@@ -389,7 +392,7 @@ export function multiplyScalar(out: Mat2d, a: Mat2d, b: number): Mat2d {
  * @param scale the amount to scale b's elements by before adding
  * @returns out
  */
-export function multiplyScalarAndAdd(out: Mat2d, a: Mat2d, b: Mat2d, scale: number): Mat2d {
+export function multiplyScalarAndAdd(out: Mat2d, a: RMat2d, b: RMat2d, scale: number): Mat2d {
     out[0] = a[0] + b[0] * scale;
     out[1] = a[1] + b[1] * scale;
     out[2] = a[2] + b[2] * scale;
@@ -406,7 +409,7 @@ export function multiplyScalarAndAdd(out: Mat2d, a: Mat2d, b: Mat2d, scale: numb
  * @param b The second matrix.
  * @returns True if the matrices are equal, false otherwise.
  */
-export function exactEquals(a: Mat2d, b: Mat2d): boolean {
+export function exactEquals(a: RMat2d, b: RMat2d): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2] && a[3] === b[3] && a[4] === b[4] && a[5] === b[5];
 }
 
@@ -417,7 +420,7 @@ export function exactEquals(a: Mat2d, b: Mat2d): boolean {
  * @param b The second matrix.
  * @returns True if the matrices are equal, false otherwise.
  */
-export function equals(a: Mat2d, b: Mat2d): boolean {
+export function equals(a: RMat2d, b: RMat2d): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
diff --git a/src/core/mat3.ts b/src/core/mat3.ts
index 1719310..a8cfef9 100644
--- a/src/core/mat3.ts
+++ b/src/core/mat3.ts
@@ -1,12 +1,15 @@
 import { EPSILON } from './scalar';
-import type { Mat2d } from './mat2d';
-import type { Mat4 } from './mat4';
-import type { Quat } from './quat';
-import type { Vec2 } from './vec2';
+import type { RMat2d } from './mat2d';
+import type { RMat4 } from './mat4';
+import type { RQuat } from './quat';
+import type { RVec2 } from './vec2';
 
 /** A 3x3 matrix */
 export type Mat3 = [e1: number, e2: number, e3: number, e4: number, e5: number, e6: number, e7: number, e8: number, e9: number];
 
+/** A read-only 3x3 matrix */
+export type RMat3 = Readonly<Mat3>;
+
 /**
  * Creates a new identity mat3
  *
@@ -23,7 +26,7 @@ export function create(): Mat3 {
  * @param a   the source 4x4 matrix
  * @returns out
  */
-export function fromMat4(out: Mat3, a: Mat4): Mat3 {
+export function fromMat4(out: Mat3, a: RMat4): Mat3 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -42,7 +45,7 @@ export function fromMat4(out: Mat3, a: Mat4): Mat3 {
  * @param a matrix to clone
  * @returns a new 3x3 matrix
  */
-export function clone(a: Mat3): Mat3 {
+export function clone(a: RMat3): Mat3 {
     return [a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8]];
 }
 
@@ -53,7 +56,7 @@ export function clone(a: Mat3): Mat3 {
  * @param a the source matrix
  * @returns out
  */
-export function copy(out: Mat3, a: Mat3): Mat3 {
+export function copy(out: Mat3, a: RMat3): Mat3 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -178,7 +181,7 @@ export function zero(out: Mat3): Mat3 {
  * @param a the source matrix
  * @returns out
  */
-export function transpose(out: Mat3, a: Mat3): Mat3 {
+export function transpose(out: Mat3, a: RMat3): Mat3 {
     // If we are transposing ourselves we can skip a few steps but have to cache some values
     if (out === a) {
         const a01 = a[1];
@@ -212,7 +215,7 @@ export function transpose(out: Mat3, a: Mat3): Mat3 {
  * @param a the source matrix
  * @returns out
  */
-export function invert(out: Mat3, a: Mat3): Mat3 | null {
+export function invert(out: Mat3, a: RMat3): Mat3 | null {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -254,7 +257,7 @@ export function invert(out: Mat3, a: Mat3): Mat3 | null {
  * @param a the source matrix
  * @returns out
  */
-export function adjoint(out: Mat3, a: Mat3): Mat3 {
+export function adjoint(out: Mat3, a: RMat3): Mat3 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -283,7 +286,7 @@ export function adjoint(out: Mat3, a: Mat3): Mat3 {
  * @param a the source matrix
  * @returns determinant of a
  */
-export function determinant(a: Mat3): number {
+export function determinant(a: RMat3): number {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -305,7 +308,7 @@ export function determinant(a: Mat3): number {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Mat3, a: Mat3, b: Mat3): Mat3 {
+export function multiply(out: Mat3, a: RMat3, b: RMat3): Mat3 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -348,7 +351,7 @@ export function multiply(out: Mat3, a: Mat3, b: Mat3): Mat3 {
  * @param v vector to translate by
  * @returns out
  */
-export function translate(out: Mat3, a: Mat3, v: Vec2): Mat3 {
+export function translate(out: Mat3, a: RMat3, v: RVec2): Mat3 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -383,7 +386,7 @@ export function translate(out: Mat3, a: Mat3, v: Vec2): Mat3 {
  * @param rad the angle to rotate the matrix by
  * @returns out
  */
-export function rotate(out: Mat3, a: Mat3, rad: number): Mat3 {
+export function rotate(out: Mat3, a: RMat3, rad: number): Mat3 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -418,7 +421,7 @@ export function rotate(out: Mat3, a: Mat3, rad: number): Mat3 {
  * @param v the vec2 to scale the matrix by
  * @returns out
  **/
-export function scale(out: Mat3, a: Mat3, v: Vec2): Mat3 {
+export function scale(out: Mat3, a: RMat3, v: RVec2): Mat3 {
     const x = v[0];
     const y = v[1];
 
@@ -447,7 +450,7 @@ export function scale(out: Mat3, a: Mat3, v: Vec2): Mat3 {
  * @param v Translation vector
  * @returns out
  */
-export function fromTranslation(out: Mat3, v: Vec2): Mat3 {
+export function fromTranslation(out: Mat3, v: RVec2): Mat3 {
     out[0] = 1;
     out[1] = 0;
     out[2] = 0;
@@ -500,7 +503,7 @@ export function fromRotation(out: Mat3, rad: number): Mat3 {
  * @param v Scaling vector
  * @returns out
  */
-export function fromScaling(out: Mat3, v: Vec2): Mat3 {
+export function fromScaling(out: Mat3, v: RVec2): Mat3 {
     out[0] = v[0];
     out[1] = 0;
     out[2] = 0;
@@ -522,7 +525,7 @@ export function fromScaling(out: Mat3, v: Vec2): Mat3 {
  * @param a the matrix to copy
  * @returns out
  **/
-export function fromMat2d(out: Mat3, a: Mat2d): Mat3 {
+export function fromMat2d(out: Mat3, a: RMat2d): Mat3 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = 0;
@@ -545,7 +548,7 @@ export function fromMat2d(out: Mat3, a: Mat2d): Mat3 {
  *
  * @returns out
  */
-export function fromQuat(out: Mat3, q: Quat): Mat3 {
+export function fromQuat(out: Mat3, q: RQuat): Mat3 {
     const x = q[0];
     const y = q[1];
     const z = q[2];
@@ -587,7 +590,7 @@ export function fromQuat(out: Mat3, q: Quat): Mat3 {
  *
  * @returns out
  */
-export function normalFromMat4(out: Mat3, a: Mat4): Mat3 | null {
+export function normalFromMat4(out: Mat3, a: RMat4): Mat3 | null {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -668,7 +671,7 @@ export function projection(out: Mat3, width: number, height: number): Mat3 {
  * @param a matrix to represent as a string
  * @returns string representation of the matrix
  */
-export function str(a: Mat3): string {
+export function str(a: RMat3): string {
     return `mat3(${a[0]}, ${a[1]}, ${a[2]}, ${a[3]}, ${a[4]}, ${a[5]}, ${a[6]}, ${a[7]}, ${a[8]})`;
 }
 
@@ -678,7 +681,7 @@ export function str(a: Mat3): string {
  * @param a the matrix to calculate Frobenius norm of
  * @returns Frobenius norm
  */
-export function frob(a: Mat3): number {
+export function frob(a: RMat3): number {
     return Math.sqrt(
         a[0] * a[0] +
             a[1] * a[1] +
@@ -700,7 +703,7 @@ export function frob(a: Mat3): number {
  * @param b the second operand
  * @returns out
  */
-export function add(out: Mat3, a: Mat3, b: Mat3): Mat3 {
+export function add(out: Mat3, a: RMat3, b: RMat3): Mat3 {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -721,7 +724,7 @@ export function add(out: Mat3, a: Mat3, b: Mat3): Mat3 {
  * @param b the second operand
  * @returns out
  */
-export function subtract(out: Mat3, a: Mat3, b: Mat3): Mat3 {
+export function subtract(out: Mat3, a: RMat3, b: RMat3): Mat3 {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     out[2] = a[2] - b[2];
@@ -742,7 +745,7 @@ export function subtract(out: Mat3, a: Mat3, b: Mat3): Mat3 {
  * @param b amount to scale the matrix's elements by
  * @returns out
  */
-export function multiplyScalar(out: Mat3, a: Mat3, b: number): Mat3 {
+export function multiplyScalar(out: Mat3, a: RMat3, b: number): Mat3 {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     out[2] = a[2] * b;
@@ -764,7 +767,7 @@ export function multiplyScalar(out: Mat3, a: Mat3, b: number): Mat3 {
  * @param scale the amount to scale b's elements by before adding
  * @returns out
  */
-export function multiplyScalarAndAdd(out: Mat3, a: Mat3, b: Mat3, scale: number): Mat3 {
+export function multiplyScalarAndAdd(out: Mat3, a: RMat3, b: RMat3, scale: number): Mat3 {
     out[0] = a[0] + b[0] * scale;
     out[1] = a[1] + b[1] * scale;
     out[2] = a[2] + b[2] * scale;
@@ -784,7 +787,7 @@ export function multiplyScalarAndAdd(out: Mat3, a: Mat3, b: Mat3, scale: number)
  * @param b The second matrix.
  * @returns True if the matrices are equal, false otherwise.
  */
-export function exactEquals(a: Mat3, b: Mat3): boolean {
+export function exactEquals(a: RMat3, b: RMat3): boolean {
     return (
         a[0] === b[0] &&
         a[1] === b[1] &&
@@ -805,7 +808,7 @@ export function exactEquals(a: Mat3, b: Mat3): boolean {
  * @param b The second matrix.
  * @returns True if the matrices are equal, false otherwise.
  */
-export function equals(a: Mat3, b: Mat3): boolean {
+export function equals(a: RMat3, b: RMat3): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
diff --git a/src/core/mat4.ts b/src/core/mat4.ts
index 21b2a9a..3d1541b 100644
--- a/src/core/mat4.ts
+++ b/src/core/mat4.ts
@@ -1,7 +1,7 @@
 import { EPSILON } from './scalar';
-import type { Quat } from './quat';
-import type { Quat2 } from './quat2';
-import type { Vec3 } from './vec3';
+import type { Quat, RQuat } from './quat';
+import type { RQuat2 } from './quat2';
+import type { RVec3, Vec3 } from './vec3';
 
 /** A 4x4 matrix */
 export type Mat4 = [
@@ -23,6 +23,9 @@ export type Mat4 = [
     e16: number,
 ];
 
+/** A read-only 4x4 matrix */
+export type RMat4 = Readonly<Mat4>;
+
 /**
  * Creates a new identity mat4
  *
@@ -38,7 +41,7 @@ export function create(): Mat4 {
  * @param a matrix to clone
  * @returns a new 4x4 matrix
  */
-export function clone(a: Mat4): Mat4 {
+export function clone(a: RMat4): Mat4 {
     return [a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10], a[11], a[12], a[13], a[14], a[15]];
 }
 
@@ -49,7 +52,7 @@ export function clone(a: Mat4): Mat4 {
  * @param a the source matrix
  * @returns out
  */
-export function copy(out: Mat4, a: Mat4): Mat4 {
+export function copy(out: Mat4, a: RMat4): Mat4 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -230,7 +233,7 @@ export function zero(out: Mat4): Mat4 {
  * @param a the source matrix
  * @returns out
  */
-export function transpose(out: Mat4, a: Mat4): Mat4 {
+export function transpose(out: Mat4, a: RMat4): Mat4 {
     // If we are transposing ourselves we can skip a few steps but have to cache some values
     if (out === a) {
         const a01 = a[1];
@@ -281,7 +284,7 @@ export function transpose(out: Mat4, a: Mat4): Mat4 {
  * @param a the source matrix
  * @returns out, or null if source matrix is not invertible
  */
-export function invert(out: Mat4, a: Mat4): Mat4 | null {
+export function invert(out: Mat4, a: RMat4): Mat4 | null {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -349,7 +352,7 @@ export function invert(out: Mat4, a: Mat4): Mat4 | null {
  * @param a the source matrix
  * @returns out, or null if the 3x3 part is not invertible
  */
-export function invert3x3(out: Mat4, a: Mat4): Mat4 | null {
+export function invert3x3(out: Mat4, a: RMat4): Mat4 | null {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -398,7 +401,7 @@ export function invert3x3(out: Mat4, a: Mat4): Mat4 | null {
  * @param a the source matrix
  * @returns out
  */
-export function adjoint(out: Mat4, a: Mat4): Mat4 {
+export function adjoint(out: Mat4, a: RMat4): Mat4 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -454,7 +457,7 @@ export function adjoint(out: Mat4, a: Mat4): Mat4 {
  * @param a the source matrix
  * @returns determinant of a
  */
-export function determinant(a: Mat4): number {
+export function determinant(a: RMat4): number {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -495,7 +498,7 @@ export function determinant(a: Mat4): number {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Mat4, a: Mat4, b: Mat4): Mat4 {
+export function multiply(out: Mat4, a: RMat4, b: RMat4): Mat4 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -562,7 +565,7 @@ export function multiply(out: Mat4, a: Mat4, b: Mat4): Mat4 {
  * @param b the second operand
  * @returns out
  */
-export function multiply3x3(out: Mat4, a: Mat4, b: Mat4): Mat4 {
+export function multiply3x3(out: Mat4, a: RMat4, b: RMat4): Mat4 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -616,7 +619,7 @@ export function multiply3x3(out: Mat4, a: Mat4, b: Mat4): Mat4 {
  * @param b the second operand (will be transposed)
  * @returns out
  */
-export function multiply3x3RightTransposed(out: Mat4, a: Mat4, b: Mat4): Mat4 {
+export function multiply3x3RightTransposed(out: Mat4, a: RMat4, b: RMat4): Mat4 {
     const a00 = a[0];
     const a01 = a[1];
     const a02 = a[2];
@@ -667,7 +670,7 @@ export function multiply3x3RightTransposed(out: Mat4, a: Mat4, b: Mat4): Mat4 {
  * @param vec the vector to transform
  * @returns out
  */
-export function multiply3x3TransposedVec(out: Vec3, mat: Mat4, vec: Vec3): Vec3 {
+export function multiply3x3TransposedVec(out: Vec3, mat: RMat4, vec: RVec3): Vec3 {
     const x = vec[0];
     const y = vec[1];
     const z = vec[2];
@@ -687,7 +690,7 @@ export function multiply3x3TransposedVec(out: Vec3, mat: Mat4, vec: Vec3): Vec3
  * @param vec the vector to transform
  * @returns out
  */
-export function multiply3x3Vec(out: Vec3, mat: Mat4, vec: Vec3): Vec3 {
+export function multiply3x3Vec(out: Vec3, mat: RMat4, vec: RVec3): Vec3 {
     const x = vec[0];
     const y = vec[1];
     const z = vec[2];
@@ -707,7 +710,7 @@ export function multiply3x3Vec(out: Vec3, mat: Mat4, vec: Vec3): Vec3 {
  * @param v the vector to create the cross product matrix from
  * @returns out
  */
-export function crossProductMatrix(out: Mat4, v: Vec3): Mat4 {
+export function crossProductMatrix(out: Mat4, v: RVec3): Mat4 {
     const x = v[0];
     const y = v[1];
     const z = v[2];
@@ -741,7 +744,7 @@ export function crossProductMatrix(out: Mat4, v: Vec3): Mat4 {
  * @param v vector to translate by
  * @returns out
  */
-export function translate(out: Mat4, a: Mat4, v: Vec3): Mat4 {
+export function translate(out: Mat4, a: RMat4, v: RVec3): Mat4 {
     const x = v[0];
     const y = v[1];
     const z = v[2];
@@ -807,7 +810,7 @@ export function translate(out: Mat4, a: Mat4, v: Vec3): Mat4 {
  * @param v the vec3 to scale the matrix by
  * @returns out
  **/
-export function scale(out: Mat4, a: Mat4, v: Vec3): Mat4 {
+export function scale(out: Mat4, a: RMat4, v: RVec3): Mat4 {
     const x = v[0];
     const y = v[1];
     const z = v[2];
@@ -840,7 +843,7 @@ export function scale(out: Mat4, a: Mat4, v: Vec3): Mat4 {
  * @param axis the axis to rotate around
  * @returns out
  */
-export function rotate(out: Mat4, a: Mat4, rad: number, axis: Vec3): Mat4 | null {
+export function rotate(out: Mat4, a: RMat4, rad: number, axis: RVec3): Mat4 | null {
     let x = axis[0];
     let y = axis[1];
     let z = axis[2];
@@ -915,7 +918,7 @@ export function rotate(out: Mat4, a: Mat4, rad: number, axis: Vec3): Mat4 | null
  * @param rad the angle to rotate the matrix by
  * @returns out
  */
-export function rotateX(out: Mat4, a: Mat4, rad: number): Mat4 {
+export function rotateX(out: Mat4, a: RMat4, rad: number): Mat4 {
     const s = Math.sin(rad);
     const c = Math.cos(rad);
     const a10 = a[4];
@@ -959,7 +962,7 @@ export function rotateX(out: Mat4, a: Mat4, rad: number): Mat4 {
  * @param rad the angle to rotate the matrix by
  * @returns out
  */
-export function rotateY(out: Mat4, a: Mat4, rad: number): Mat4 {
+export function rotateY(out: Mat4, a: RMat4, rad: number): Mat4 {
     const s = Math.sin(rad);
     const c = Math.cos(rad);
     const a00 = a[0];
@@ -1003,7 +1006,7 @@ export function rotateY(out: Mat4, a: Mat4, rad: number): Mat4 {
  * @param rad the angle to rotate the matrix by
  * @returns out
  */
-export function rotateZ(out: Mat4, a: Mat4, rad: number): Mat4 {
+export function rotateZ(out: Mat4, a: RMat4, rad: number): Mat4 {
     const s = Math.sin(rad);
     const c = Math.cos(rad);
     const a00 = a[0];
@@ -1050,7 +1053,7 @@ export function rotateZ(out: Mat4, a: Mat4, rad: number): Mat4 {
  * @param v Translation vector
  * @returns out
  */
-export function fromTranslation(out: Mat4, v: Vec3): Mat4 {
+export function fromTranslation(out: Mat4, v: RVec3): Mat4 {
     out[0] = 1;
     out[1] = 0;
     out[2] = 0;
@@ -1081,7 +1084,7 @@ export function fromTranslation(out: Mat4, v: Vec3): Mat4 {
  * @param v Scaling vector
  * @returns out
  */
-export function fromScaling(out: Mat4, v: Vec3): Mat4 {
+export function fromScaling(out: Mat4, v: RVec3): Mat4 {
     out[0] = v[0];
     out[1] = 0;
     out[2] = 0;
@@ -1113,7 +1116,7 @@ export function fromScaling(out: Mat4, v: Vec3): Mat4 {
  * @param axis the axis to rotate around
  * @returns out
  */
-export function fromRotation(out: Mat4, rad: number, axis: Vec3): Mat4 | null {
+export function fromRotation(out: Mat4, rad: number, axis: RVec3): Mat4 | null {
     let x = axis[0];
     let y = axis[1];
     let z = axis[2];
@@ -1272,7 +1275,7 @@ export function fromZRotation(out: Mat4, rad: number): Mat4 {
  * @param v Translation vector
  * @returns out
  */
-export function fromRotationTranslation(out: Mat4, q: Quat | Quat2, v: Vec3): Mat4 {
+export function fromRotationTranslation(out: Mat4, q: RQuat | RQuat2, v: RVec3): Mat4 {
     // Quaternion math
     const x = q[0];
     const y = q[1];
@@ -1319,7 +1322,7 @@ export function fromRotationTranslation(out: Mat4, q: Quat | Quat2, v: Vec3): Ma
  * @param a Dual Quaternion
  * @returns mat4 receiving operation result
  */
-export function fromQuat2(out: Mat4, a: Quat2): Mat4 {
+export function fromQuat2(out: Mat4, a: RQuat2): Mat4 {
     const translation = [0, 0, 0] as Vec3;
     const bx = -a[0];
     const by = -a[1];
@@ -1354,7 +1357,7 @@ export function fromQuat2(out: Mat4, a: Quat2): Mat4 {
  * @param mat Matrix to be decomposed (input)
  * @return out
  */
-export function getTranslation(out: Vec3, mat: Mat4): Vec3 {
+export function getTranslation(out: Vec3, mat: RMat4): Vec3 {
     out[0] = mat[12];
     out[1] = mat[13];
     out[2] = mat[14];
@@ -1372,7 +1375,7 @@ export function getTranslation(out: Vec3, mat: Mat4): Vec3 {
  * @param mat Matrix to be decomposed (input)
  * @return out
  */
-export function getScaling(out: Vec3, mat: Mat4): Vec3 {
+export function getScaling(out: Vec3, mat: RMat4): Vec3 {
     const m11 = mat[0];
     const m12 = mat[1];
     const m13 = mat[2];
@@ -1399,7 +1402,7 @@ export function getScaling(out: Vec3, mat: Mat4): Vec3 {
  * @param mat Matrix to be decomposed (input)
  * @return out
  */
-export function getRotation(out: Quat, mat: Mat4): Quat {
+export function getRotation(out: Quat, mat: RMat4): Quat {
     const scaling = [0, 0, 0] as Vec3;
     getScaling(scaling, mat);
 
@@ -1458,7 +1461,7 @@ export function getRotation(out: Quat, mat: Mat4): Quat {
  * @param mat Matrix to be decomposed (input)
  * @returns out_r
  */
-export function decompose(out_r: Quat, out_t: Vec3, out_s: Vec3, mat: Mat4): Quat {
+export function decompose(out_r: Quat, out_t: Vec3, out_s: Vec3, mat: RMat4): Quat {
     out_t[0] = mat[12];
     out_t[1] = mat[13];
     out_t[2] = mat[14];
@@ -1540,7 +1543,7 @@ export function decompose(out_r: Quat, out_t: Vec3, out_s: Vec3, mat: Mat4): Qua
  * @param s Scaling vector
  * @returns out
  */
-export function fromRotationTranslationScale(out: Mat4, q: Quat, v: Vec3, s: Vec3): Mat4 {
+export function fromRotationTranslationScale(out: Mat4, q: RQuat, v: RVec3, s: RVec3): Mat4 {
     // Quaternion math
     const x = q[0];
     const y = q[1];
@@ -1603,7 +1606,7 @@ export function fromRotationTranslationScale(out: Mat4, q: Quat, v: Vec3, s: Vec
  * @param o The origin vector around which to scale and rotate
  * @returns out
  */
-export function fromRotationTranslationScaleOrigin(out: Mat4, q: Quat, v: Vec3, s: Vec3, o: Vec3): Mat4 {
+export function fromRotationTranslationScaleOrigin(out: Mat4, q: RQuat, v: RVec3, s: RVec3, o: RVec3): Mat4 {
     // Quaternion math
     const x = q[0];
     const y = q[1];
@@ -1669,7 +1672,7 @@ export function fromRotationTranslationScaleOrigin(out: Mat4, q: Quat, v: Vec3,
  *
  * @returns out
  */
-export function fromQuat(out: Mat4, q: Quat): Mat4 {
+export function fromQuat(out: Mat4, q: RQuat): Mat4 {
     const x = q[0];
     const y = q[1];
     const z = q[2];
@@ -2040,7 +2043,7 @@ export function orthoZO(out: Mat4, left: number, right: number, bottom: number,
  * @param up vec3 pointing up
  * @returns out
  */
-export function lookAt(out: Mat4, eye: Vec3, center: Vec3, up: Vec3): Mat4 {
+export function lookAt(out: Mat4, eye: RVec3, center: RVec3, up: RVec3): Mat4 {
     let x0: number;
     let x1: number;
     let x2: number;
@@ -2134,7 +2137,7 @@ export function lookAt(out: Mat4, eye: Vec3, center: Vec3, up: Vec3): Mat4 {
  * @param up vec3 pointing up
  * @returns out
  */
-export function targetTo(out: Mat4, eye: Vec3, target: Vec3, up: Vec3): Mat4 {
+export function targetTo(out: Mat4, eye: RVec3, target: RVec3, up: RVec3): Mat4 {
     const eyex = eye[0];
     const eyey = eye[1];
     const eyez = eye[2];
@@ -2191,7 +2194,7 @@ export function targetTo(out: Mat4, eye: Vec3, target: Vec3, up: Vec3): Mat4 {
  * @param a matrix to represent as a string
  * @returns {String} string representation of the matrix
  */
-export function str(a: Mat4): string {
+export function str(a: RMat4): string {
     return `mat4(${a[0]}, ${a[1]}, ${a[2]}, ${a[3]}, ${a[4]}, ${a[5]}, ${a[6]}, ${a[7]}, ${a[8]}, ${a[9]}, ${a[10]}, ${a[11]}, ${a[12]}, ${a[13]}, ${a[14]}, ${a[15]})`;
 }
 
@@ -2201,7 +2204,7 @@ export function str(a: Mat4): string {
  * @param a the matrix to calculate Frobenius norm of
  * @returns Frobenius norm
  */
-export function frob(a: Mat4): number {
+export function frob(a: RMat4): number {
     return Math.sqrt(
         a[0] * a[0] +
             a[1] * a[1] +
@@ -2230,7 +2233,7 @@ export function frob(a: Mat4): number {
  * @param b the second operand
  * @returns out
  */
-export function add(out: Mat4, a: Mat4, b: Mat4): Mat4 {
+export function add(out: Mat4, a: RMat4, b: RMat4): Mat4 {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -2258,7 +2261,7 @@ export function add(out: Mat4, a: Mat4, b: Mat4): Mat4 {
  * @param b the second operand
  * @returns out
  */
-export function subtract(out: Mat4, a: Mat4, b: Mat4): Mat4 {
+export function subtract(out: Mat4, a: RMat4, b: RMat4): Mat4 {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     out[2] = a[2] - b[2];
@@ -2286,7 +2289,7 @@ export function subtract(out: Mat4, a: Mat4, b: Mat4): Mat4 {
  * @param b amount to scale the matrix's elements by
  * @returns out
  */
-export function multiplyScalar(out: Mat4, a: Mat4, b: number): Mat4 {
+export function multiplyScalar(out: Mat4, a: RMat4, b: number): Mat4 {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     out[2] = a[2] * b;
@@ -2315,7 +2318,7 @@ export function multiplyScalar(out: Mat4, a: Mat4, b: number): Mat4 {
  * @param scale the amount to scale b's elements by before adding
  * @returns out
  */
-export function multiplyScalarAndAdd(out: Mat4, a: Mat4, b: Mat4, scale: number): Mat4 {
+export function multiplyScalarAndAdd(out: Mat4, a: RMat4, b: RMat4, scale: number): Mat4 {
     out[0] = a[0] + b[0] * scale;
     out[1] = a[1] + b[1] * scale;
     out[2] = a[2] + b[2] * scale;
@@ -2342,7 +2345,7 @@ export function multiplyScalarAndAdd(out: Mat4, a: Mat4, b: Mat4, scale: number)
  * @param b The second matrix.
  * @returns {Boolean} True if the matrices are equal, false otherwise.
  */
-export function exactEquals(a: Mat4, b: Mat4): boolean {
+export function exactEquals(a: RMat4, b: RMat4): boolean {
     return (
         a[0] === b[0] &&
         a[1] === b[1] &&
@@ -2370,7 +2373,7 @@ export function exactEquals(a: Mat4, b: Mat4): boolean {
  * @param b The second matrix.
  * @returns {Boolean} True if the matrices are equal, false otherwise.
  */
-export function equals(a: Mat4, b: Mat4): boolean {
+export function equals(a: RMat4, b: RMat4): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
diff --git a/src/core/polar.ts b/src/core/polar.ts
index aad176b..e42082e 100644
--- a/src/core/polar.ts
+++ b/src/core/polar.ts
@@ -1,6 +1,6 @@
 import * as scalar from './scalar';
 import { wrapAngle } from './angle';
-import type { Vec2 } from './vec2';
+import type { RVec2, Vec2 } from './vec2';
 
 /**
  * A point in polar coordinates [r, theta]
@@ -14,6 +14,9 @@ import type { Vec2 } from './vec2';
  */
 export type Polar = [r: number, theta: number];
 
+/** A read-only polar coordinate */
+export type RPolar = Readonly<Polar>;
+
 /**
  * Creates a new polar coordinate at r=1, theta=0
  *
@@ -40,7 +43,7 @@ export function fromValues(r: number, theta: number): Polar {
  * @param a the source Polar
  * @returns a new Polar
  */
-export function clone(a: Polar): Polar {
+export function clone(a: RPolar): Polar {
     return [a[0], a[1]];
 }
 
@@ -51,7 +54,7 @@ export function clone(a: Polar): Polar {
  * @param a the source Polar
  * @returns out
  */
-export function copy(out: Polar, a: Polar): Polar {
+export function copy(out: Polar, a: RPolar): Polar {
     out[0] = a[0];
     out[1] = a[1];
     return out;
@@ -78,7 +81,7 @@ export function set(out: Polar, r: number, theta: number): Polar {
  * @param a the source Polar
  * @returns out
  */
-export function normalize(out: Polar, a: Polar): Polar {
+export function normalize(out: Polar, a: RPolar): Polar {
     out[0] = 1;
     out[1] = a[1];
     return out;
@@ -92,7 +95,7 @@ export function normalize(out: Polar, a: Polar): Polar {
  * @param s scalar to multiply r by
  * @returns out
  */
-export function scale(out: Polar, a: Polar, s: number): Polar {
+export function scale(out: Polar, a: RPolar, s: number): Polar {
     out[0] = a[0] * s;
     out[1] = a[1];
     return out;
@@ -106,7 +109,7 @@ export function scale(out: Polar, a: Polar, s: number): Polar {
  * @param rad the angle to add to theta
  * @returns out
  */
-export function rotate(out: Polar, a: Polar, rad: number): Polar {
+export function rotate(out: Polar, a: RPolar, rad: number): Polar {
     out[0] = a[0];
     out[1] = wrapAngle(a[1] + rad);
     return out;
@@ -122,7 +125,7 @@ export function rotate(out: Polar, a: Polar, rad: number): Polar {
  * @param t interpolation factor in [0, 1]
  * @returns out
  */
-export function lerp(out: Polar, a: Polar, b: Polar, t: number): Polar {
+export function lerp(out: Polar, a: RPolar, b: RPolar, t: number): Polar {
     out[0] = scalar.lerp(a[0], b[0], t);
     out[1] = a[1] + wrapAngle(b[1] - a[1]) * t;
     return out;
@@ -137,7 +140,7 @@ export function lerp(out: Polar, a: Polar, b: Polar, t: number): Polar {
  * @param v the source Vec2
  * @returns out
  */
-export function setFromVec2(out: Polar, v: Vec2): Polar {
+export function setFromVec2(out: Polar, v: RVec2): Polar {
     const x = v[0];
     const y = v[1];
     out[0] = Math.sqrt(x * x + y * y);
@@ -157,7 +160,7 @@ export const fromVec2 = setFromVec2;
  * @param a the source Polar
  * @returns out
  */
-export function toVec2(out: Vec2, a: Polar): Vec2 {
+export function toVec2(out: Vec2, a: RPolar): Vec2 {
     const r = a[0];
     const theta = a[1];
     out[0] = r * Math.cos(theta);
@@ -173,7 +176,7 @@ export function toVec2(out: Vec2, a: Polar): Vec2 {
  * @param b the second Polar
  * @returns angle in radians in [0, pi]
  */
-export function angleTo(a: Polar, b: Polar): number {
+export function angleTo(a: RPolar, b: RPolar): number {
     return Math.abs(wrapAngle(b[1] - a[1]));
 }
 
@@ -185,7 +188,7 @@ export function angleTo(a: Polar, b: Polar): number {
  * @param b the second Polar
  * @returns the Euclidean distance between the two points
  */
-export function distance(a: Polar, b: Polar): number {
+export function distance(a: RPolar, b: RPolar): number {
     const ra = a[0];
     const rb = b[0];
     const d = ra * ra + rb * rb - 2 * ra * rb * Math.cos(b[1] - a[1]);
@@ -200,7 +203,7 @@ export function distance(a: Polar, b: Polar): number {
  * @param b the second Polar
  * @returns true if approximately equal
  */
-export function equals(a: Polar, b: Polar): boolean {
+export function equals(a: RPolar, b: RPolar): boolean {
     return scalar.equals(a[0], b[0]) && scalar.equals(a[1], b[1]);
 }
 
@@ -211,7 +214,7 @@ export function equals(a: Polar, b: Polar): boolean {
  * @param b the second Polar
  * @returns true if exactly equal
  */
-export function exactEquals(a: Polar, b: Polar): boolean {
+export function exactEquals(a: RPolar, b: RPolar): boolean {
     return a[0] === b[0] && a[1] === b[1];
 }
 
@@ -221,6 +224,6 @@ export function exactEquals(a: Polar, b: Polar): boolean {
  * @param a the source Polar
  * @returns string representation
  */
-export function str(a: Polar): string {
+export function str(a: RPolar): string {
     return `Polar(${a[0]}, ${a[1]})`;
 }
diff --git a/src/core/quat.ts b/src/core/quat.ts
index 3b0d86e..83f5bd3 100644
--- a/src/core/quat.ts
+++ b/src/core/quat.ts
@@ -1,16 +1,19 @@
 import type { MutableArrayLike } from './arrays';
-import type { Euler, EulerOrder } from './euler';
-import type { Mat3 } from './mat3';
+import type { Euler, EulerOrder, REuler } from './euler';
+import type { RMat3 } from './mat3';
 import * as mat3 from './mat3';
-import type { Mat4 } from './mat4';
+import type { RMat4 } from './mat4';
 import { EPSILON } from './scalar';
-import type { Vec3 } from './vec3';
+import type { RVec3, Vec3 } from './vec3';
 import * as vec3 from './vec3';
 import * as vec4 from './vec4';
 
 /** A quaternion that represents rotation */
 export type Quat = [x: number, y: number, z: number, w: number];
 
+/** A read-only quaternion */
+export type RQuat = Readonly<Quat>;
+
 /**
  * Creates a new identity quat
  *
@@ -42,7 +45,7 @@ export function fromBuffer(out: Quat, buffer: ArrayLike<number>, startIndex: num
  * @param startIndex The starting index in the buffer
  * @returns The output buffer
  */
-export function toBuffer(outBuffer: MutableArrayLike<number>, q: Quat, startIndex: number): MutableArrayLike<number> {
+export function toBuffer(outBuffer: MutableArrayLike<number>, q: RQuat, startIndex: number): MutableArrayLike<number> {
     outBuffer[startIndex] = q[0];
     outBuffer[startIndex + 1] = q[1];
     outBuffer[startIndex + 2] = q[2];
@@ -73,7 +76,7 @@ export function identity(out: Quat): Quat {
  * @param rad the angle in radians
  * @returns out
  **/
-export function setAxisAngle(out: Quat, axis: Vec3, rad: number): Quat {
+export function setAxisAngle(out: Quat, axis: RVec3, rad: number): Quat {
     rad *= 0.5;
     const s = Math.sin(rad);
     out[0] = s * axis[0];
@@ -96,7 +99,7 @@ export function setAxisAngle(out: Quat, axis: Vec3, rad: number): Quat {
  * @param  q     Quaternion to be decomposed
  * @return     Angle, in radians, of the rotation
  */
-export function getAxisAngle(out_axis: Vec3, q: Quat): number {
+export function getAxisAngle(out_axis: Vec3, q: RQuat): number {
     const rad = Math.acos(q[3]) * 2.0;
     const s = Math.sin(rad / 2.0);
     if (s > EPSILON) {
@@ -119,7 +122,7 @@ export function getAxisAngle(out_axis: Vec3, q: Quat): number {
  * @param  b     Destination unit quaternion
  * @return     Angle, in radians, between the two quaternions
  */
-export function getAngle(a: Quat, b: Quat): number {
+export function getAngle(a: RQuat, b: RQuat): number {
     const dotproduct = dot(a, b);
 
     return Math.acos(2 * dotproduct * dotproduct - 1);
@@ -133,7 +136,7 @@ export function getAngle(a: Quat, b: Quat): number {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Quat, a: Quat, b: Quat): Quat {
+export function multiply(out: Quat, a: RQuat, b: RQuat): Quat {
     const ax = a[0];
     const ay = a[1];
     const az = a[2];
@@ -158,7 +161,7 @@ export function multiply(out: Quat, a: Quat, b: Quat): Quat {
  * @param rad angle (in radians) to rotate
  * @returns out
  */
-export function rotateX(out: Quat, a: Quat, rad: number): Quat {
+export function rotateX(out: Quat, a: RQuat, rad: number): Quat {
     rad *= 0.5;
 
     const ax = a[0];
@@ -183,7 +186,7 @@ export function rotateX(out: Quat, a: Quat, rad: number): Quat {
  * @param rad angle (in radians) to rotate
  * @returns out
  */
-export function rotateY(out: Quat, a: Quat, rad: number): Quat {
+export function rotateY(out: Quat, a: RQuat, rad: number): Quat {
     rad *= 0.5;
 
     const ax = a[0];
@@ -208,7 +211,7 @@ export function rotateY(out: Quat, a: Quat, rad: number): Quat {
  * @param rad angle (in radians) to rotate
  * @returns out
  */
-export function rotateZ(out: Quat, a: Quat, rad: number): Quat {
+export function rotateZ(out: Quat, a: RQuat, rad: number): Quat {
     rad *= 0.5;
 
     const ax = a[0];
@@ -234,7 +237,7 @@ export function rotateZ(out: Quat, a: Quat, rad: number): Quat {
  * @param a quat to calculate W component of
  * @returns out
  */
-export function calculateW(out: Quat, a: Quat): Quat {
+export function calculateW(out: Quat, a: RQuat): Quat {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -253,7 +256,7 @@ export function calculateW(out: Quat, a: Quat): Quat {
  * @param a quat to calculate the exponential of
  * @returns out
  */
-export function exp(out: Quat, a: Quat): Quat {
+export function exp(out: Quat, a: RQuat): Quat {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -278,7 +281,7 @@ export function exp(out: Quat, a: Quat): Quat {
  * @param a quat to calculate the exponential of
  * @returns out
  */
-export function ln(out: Quat, a: Quat): Quat {
+export function ln(out: Quat, a: RQuat): Quat {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -303,7 +306,7 @@ export function ln(out: Quat, a: Quat): Quat {
  * @param b amount to scale the quaternion by
  * @returns out
  */
-export function pow(out: Quat, a: Quat, b: number): Quat {
+export function pow(out: Quat, a: RQuat, b: number): Quat {
     ln(out, a);
     scale(out, out, b);
     exp(out, out);
@@ -319,7 +322,7 @@ export function pow(out: Quat, a: Quat, b: number): Quat {
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function slerp(out: Quat, a: Quat, b: Quat, t: number): Quat {
+export function slerp(out: Quat, a: RQuat, b: RQuat, t: number): Quat {
     // benchmarks:
     //    http://jsperf.com/quaternion-slerp-implementations
     const ax = a[0];
@@ -376,7 +379,7 @@ export function slerp(out: Quat, a: Quat, b: Quat, t: number): Quat {
  * @param a quat to calculate inverse of
  * @returns out
  */
-export function invert(out: Quat, a: Quat): Quat {
+export function invert(out: Quat, a: RQuat): Quat {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -407,7 +410,7 @@ export function invert(out: Quat, a: Quat): Quat {
  * @param a quat to calculate conjugate of
  * @returns out
  */
-export function conjugate(out: Quat, a: Quat): Quat {
+export function conjugate(out: Quat, a: RQuat): Quat {
     out[0] = -a[0];
     out[1] = -a[1];
     out[2] = -a[2];
@@ -425,7 +428,7 @@ export function conjugate(out: Quat, a: Quat): Quat {
  * @param m rotation matrix
  * @returns out
  */
-export function fromMat3(out: Quat, m: Mat3): Quat {
+export function fromMat3(out: Quat, m: RMat3): Quat {
     // Algorithm in Ken Shoemake's article in 1987 SIGGRAPH course notes
     // article "Quaternion Calculus and Fast Animation".
     const fTrace = m[0] + m[4] + m[8];
@@ -466,7 +469,7 @@ export function fromMat3(out: Quat, m: Mat3): Quat {
  * @param m rotation matrix
  * @returns out
  */
-export function fromMat4(out: Quat, m: Mat4): Quat {
+export function fromMat4(out: Quat, m: RMat4): Quat {
     const m3 = mat3.create();
     mat3.fromMat4(m3, m);
     return fromMat3(out, m3);
@@ -478,7 +481,7 @@ export function fromMat4(out: Quat, m: Mat4): Quat {
  * @param euler the euler to create the quaternion from
  * @returns out
  */
-export function fromEuler(out: Quat, euler: Euler): Quat {
+export function fromEuler(out: Quat, euler: REuler): Quat {
     const x = euler[0];
     const y = euler[1];
     const z = euler[2];
@@ -573,7 +576,7 @@ export function fromDegrees(out: Quat, x: number, y: number, z: number, order: E
  * @param a vector to represent as a string
  * @returns string representation of the vector
  */
-export function str(a: Quat): string {
+export function str(a: RQuat): string {
     return `quat(${a[0]}, ${a[1]}, ${a[2]}, ${a[3]})`;
 }
 
@@ -713,7 +716,7 @@ export const exactEquals = vec4.exactEquals;
  * @param b The second quaternion.
  * @returns True if the quaternions are equal, false otherwise.
  */
-export function equals(a: Quat, b: Quat): boolean {
+export function equals(a: RQuat, b: RQuat): boolean {
     return Math.abs(vec4.dot(a, b)) >= 1 - EPSILON;
 }
 
@@ -733,7 +736,7 @@ export const rotationTo = /* @__PURE__ */ (() => {
     const xUnitVec3 = vec3.fromValues(1, 0, 0);
     const yUnitVec3 = vec3.fromValues(0, 1, 0);
 
-    return (out: Quat, a: Vec3, b: Vec3): Quat => {
+    return (out: Quat, a: RVec3, b: RVec3): Quat => {
         const dot = vec3.dot(a, b);
 
         if (dot < -0.999999) {
@@ -776,7 +779,7 @@ export const sqlerp = /* @__PURE__ */ (() => {
     const temp1 = create();
     const temp2 = create();
 
-    return (out: Quat, a: Quat, b: Quat, c: Quat, d: Quat, t: number): Quat => {
+    return (out: Quat, a: RQuat, b: RQuat, c: RQuat, d: RQuat, t: number): Quat => {
         slerp(temp1, a, d, t);
         slerp(temp2, b, c, t);
         slerp(out, temp1, temp2, 2 * t * (1 - t));
@@ -798,7 +801,7 @@ export const sqlerp = /* @__PURE__ */ (() => {
 export const setAxes = /* @__PURE__ */ (() => {
     const matr = mat3.create();
 
-    return (out: Quat, view: Vec3, right: Vec3, up: Vec3): Quat => {
+    return (out: Quat, view: RVec3, right: RVec3, up: RVec3): Quat => {
         matr[0] = right[0];
         matr[3] = right[1];
         matr[6] = right[2];
diff --git a/src/core/quat2.ts b/src/core/quat2.ts
index ec4283e..76c035d 100644
--- a/src/core/quat2.ts
+++ b/src/core/quat2.ts
@@ -1,11 +1,14 @@
-import type { Mat4 } from './mat4';
-import type { Quat } from './quat';
+import type { RMat4 } from './mat4';
+import type { Quat, RQuat } from './quat';
 import { EPSILON } from './scalar';
-import type { Vec3 } from './vec3';
+import type { RVec3, Vec3 } from './vec3';
 
 /** A dual quaternion that represents both rotation and translation */
 export type Quat2 = [x: number, y: number, z: number, w: number, x2: number, y2: number, z2: number, w2: number];
 
+/** A read-only dual quaternion */
+export type RQuat2 = Readonly<Quat2>;
+
 /**
  * Creates a new identity dual quat
  *
@@ -26,7 +29,7 @@ export function create(): Quat2 {
  * @returns new dual quaternion
  * @function
  */
-export function clone(a: Quat2): Quat2 {
+export function clone(a: RQuat2): Quat2 {
     return [a[0], a[1], a[2], a[3], a[4], a[5], a[6], a[7]];
 }
 
@@ -103,7 +106,7 @@ export function fromRotationTranslationValues(
  * @returns dual quaternion receiving operation result
  * @function
  */
-export function fromRotationTranslation(out: Quat2, q: Quat, t: Vec3): Quat2 {
+export function fromRotationTranslation(out: Quat2, q: RQuat, t: RVec3): Quat2 {
     const ax = t[0] * 0.5;
     const ay = t[1] * 0.5;
     const az = t[2] * 0.5;
@@ -130,7 +133,7 @@ export function fromRotationTranslation(out: Quat2, q: Quat, t: Vec3): Quat2 {
  * @returns dual quaternion receiving operation result
  * @function
  */
-export function fromTranslation(out: Quat2, t: Vec3): Quat2 {
+export function fromTranslation(out: Quat2, t: RVec3): Quat2 {
     out[0] = 0;
     out[1] = 0;
     out[2] = 0;
@@ -150,7 +153,7 @@ export function fromTranslation(out: Quat2, t: Vec3): Quat2 {
  * @returns dual quaternion receiving operation result
  * @function
  */
-export function fromRotation(out: Quat2, q: Quat): Quat2 {
+export function fromRotation(out: Quat2, q: RQuat): Quat2 {
     out[0] = q[0];
     out[1] = q[1];
     out[2] = q[2];
@@ -169,7 +172,7 @@ export function fromRotation(out: Quat2, q: Quat): Quat2 {
  * @param a the matrix
  * @returns dual quat receiving operation result
  */
-export function fromMat4(out: Quat2, a: Mat4): Quat2 {
+export function fromMat4(out: Quat2, a: RMat4): Quat2 {
     // Rotation: extract the quaternion from the (possibly scaled) upper-3x3.
     // Inlined from mat4.getRotation/getScaling so no scratch quat or Vec3 is allocated.
     const is1 = 1 / Math.sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
@@ -242,7 +245,7 @@ export function fromMat4(out: Quat2, a: Mat4): Quat2 {
  * @returns out
  * @function
  */
-export function copy(out: Quat2, a: Quat2): Quat2 {
+export function copy(out: Quat2, a: RQuat2): Quat2 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -316,7 +319,7 @@ export function set(
  * @param  a Dual Quaternion
  * @return real part
  */
-export function getReal(out: Quat, a: Quat2): Quat {
+export function getReal(out: Quat, a: RQuat2): Quat {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -330,7 +333,7 @@ export function getReal(out: Quat, a: Quat2): Quat {
  * @param  a Dual Quaternion
  * @return dual part
  */
-export function getDual(out: Quat, a: Quat2): Quat {
+export function getDual(out: Quat, a: RQuat2): Quat {
     out[0] = a[4];
     out[1] = a[5];
     out[2] = a[6];
@@ -345,7 +348,7 @@ export function getDual(out: Quat, a: Quat2): Quat {
  * @param q a quaternion representing the real part
  * @returns out
  */
-export function setReal(out: Quat2, q: Quat): Quat2 {
+export function setReal(out: Quat2, q: RQuat): Quat2 {
     out[0] = q[0];
     out[1] = q[1];
     out[2] = q[2];
@@ -361,7 +364,7 @@ export function setReal(out: Quat2, q: Quat): Quat2 {
  * @returns out
  * @function
  */
-export function setDual(out: Quat2, q: Quat): Quat2 {
+export function setDual(out: Quat2, q: RQuat): Quat2 {
     out[4] = q[0];
     out[5] = q[1];
     out[6] = q[2];
@@ -375,7 +378,7 @@ export function setDual(out: Quat2, q: Quat): Quat2 {
  * @param  a Dual Quaternion to be decomposed
  * @return translation
  */
-export function getTranslation(out: Vec3, a: Quat2): Vec3 {
+export function getTranslation(out: Vec3, a: RQuat2): Vec3 {
     const ax = a[4];
     const ay = a[5];
     const az = a[6];
@@ -398,7 +401,7 @@ export function getTranslation(out: Vec3, a: Quat2): Vec3 {
  * @param v vector to translate by
  * @returns out
  */
-export function translate(out: Quat2, a: Quat2, v: Vec3): Quat2 {
+export function translate(out: Quat2, a: RQuat2, v: RVec3): Quat2 {
     const ax1 = a[0];
     const ay1 = a[1];
     const az1 = a[2];
@@ -429,7 +432,7 @@ export function translate(out: Quat2, a: Quat2, v: Vec3): Quat2 {
  * @param rad how far should the rotation be
  * @returns out
  */
-export function rotateX(out: Quat2, a: Quat2, rad: number): Quat2 {
+export function rotateX(out: Quat2, a: RQuat2, rad: number): Quat2 {
     let bx = -a[0];
     let by = -a[1];
     let bz = -a[2];
@@ -473,7 +476,7 @@ export function rotateX(out: Quat2, a: Quat2, rad: number): Quat2 {
  * @param rad how far should the rotation be
  * @returns out
  */
-export function rotateY(out: Quat2, a: Quat2, rad: number): Quat2 {
+export function rotateY(out: Quat2, a: RQuat2, rad: number): Quat2 {
     let bx = -a[0];
     let by = -a[1];
     let bz = -a[2];
@@ -517,7 +520,7 @@ export function rotateY(out: Quat2, a: Quat2, rad: number): Quat2 {
  * @param rad how far should the rotation be
  * @returns out
  */
-export function rotateZ(out: Quat2, a: Quat2, rad: number): Quat2 {
+export function rotateZ(out: Quat2, a: RQuat2, rad: number): Quat2 {
     let bx = -a[0];
     let by = -a[1];
     let bz = -a[2];
@@ -561,7 +564,7 @@ export function rotateZ(out: Quat2, a: Quat2, rad: number): Quat2 {
  * @param q quaternion to rotate by
  * @returns out
  */
-export function rotateByQuatAppend(out: Quat2, a: Quat2, q: Quat): Quat2 {
+export function rotateByQuatAppend(out: Quat2, a: RQuat2, q: RQuat): Quat2 {
     const qx = q[0];
     const qy = q[1];
     const qz = q[2];
@@ -594,7 +597,7 @@ export function rotateByQuatAppend(out: Quat2, a: Quat2, q: Quat): Quat2 {
  * @param a the dual quaternion to rotate
  * @returns out
  */
-export function rotateByQuatPrepend(out: Quat2, q: Quat, a: Quat2): Quat2 {
+export function rotateByQuatPrepend(out: Quat2, q: RQuat, a: RQuat2): Quat2 {
     const qx = q[0];
     const qy = q[1];
     const qz = q[2];
@@ -628,7 +631,7 @@ export function rotateByQuatPrepend(out: Quat2, q: Quat, a: Quat2): Quat2 {
  * @param rad how far the rotation should be
  * @returns out
  */
-export function rotateAroundAxis(out: Quat2, a: Quat2, axis: Vec3, rad: number): Quat2 {
+export function rotateAroundAxis(out: Quat2, a: RQuat2, axis: RVec3, rad: number): Quat2 {
     //Special case for rad = 0
     if (Math.abs(rad) < EPSILON) {
         return copy(out, a);
@@ -672,7 +675,7 @@ export function rotateAroundAxis(out: Quat2, a: Quat2, axis: Vec3, rad: number):
  * @returns out
  * @function
  */
-export function add(out: Quat2, a: Quat2, b: Quat2): Quat2 {
+export function add(out: Quat2, a: RQuat2, b: RQuat2): Quat2 {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -692,7 +695,7 @@ export function add(out: Quat2, a: Quat2, b: Quat2): Quat2 {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Quat2, a: Quat2, b: Quat2): Quat2 {
+export function multiply(out: Quat2, a: RQuat2, b: RQuat2): Quat2 {
     const ax0 = a[0];
     const ay0 = a[1];
     const az0 = a[2];
@@ -735,7 +738,7 @@ export const mul = multiply;
  * @returns out
  * @function
  */
-export function scale(out: Quat2, a: Quat2, b: number): Quat2 {
+export function scale(out: Quat2, a: RQuat2, b: number): Quat2 {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     out[2] = a[2] * b;
@@ -754,7 +757,7 @@ export function scale(out: Quat2, a: Quat2, b: number): Quat2 {
  * @param b the second operand
  * @returns dot product of a and b
  */
-export function dot(a: Quat2, b: Quat2): number {
+export function dot(a: RQuat2, b: RQuat2): number {
     return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
 }
 
@@ -768,7 +771,7 @@ export function dot(a: Quat2, b: Quat2): number {
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function lerp(out: Quat2, a: Quat2, b: Quat2, t: number): Quat2 {
+export function lerp(out: Quat2, a: RQuat2, b: RQuat2, t: number): Quat2 {
     const mt = 1 - t;
     // dot of the real (rotation) parts, matching quat2.dot
     if (a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3] < 0) t = -t;
@@ -792,7 +795,7 @@ export function lerp(out: Quat2, a: Quat2, b: Quat2, t: number): Quat2 {
  * @param a dual quat to calculate inverse of
  * @returns out
  */
-export function invert(out: Quat2, a: Quat2): Quat2 {
+export function invert(out: Quat2, a: RQuat2): Quat2 {
     const sqlen = a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3];
     out[0] = -a[0] / sqlen;
     out[1] = -a[1] / sqlen;
@@ -813,7 +816,7 @@ export function invert(out: Quat2, a: Quat2): Quat2 {
  * @param a quat to calculate conjugate of
  * @returns out
  */
-export function conjugate(out: Quat2, a: Quat2): Quat2 {
+export function conjugate(out: Quat2, a: RQuat2): Quat2 {
     out[0] = -a[0];
     out[1] = -a[1];
     out[2] = -a[2];
@@ -831,7 +834,7 @@ export function conjugate(out: Quat2, a: Quat2): Quat2 {
  * @param a dual quat to calculate length of
  * @returns length of a
  */
-export function length(a: Quat2): number {
+export function length(a: RQuat2): number {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -851,7 +854,7 @@ export const len = length;
  * @param a dual quat to calculate squared length of
  * @returns squared length of a
  */
-export function squaredLength(a: Quat2): number {
+export function squaredLength(a: RQuat2): number {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -873,7 +876,7 @@ export const sqrLen = squaredLength;
  * @returns out
  * @function
  */
-export function normalize(out: Quat2, a: Quat2): Quat2 {
+export function normalize(out: Quat2, a: RQuat2): Quat2 {
     let magnitude = a[0] * a[0] + a[1] * a[1] + a[2] * a[2] + a[3] * a[3];
     if (magnitude > 0) {
         magnitude = Math.sqrt(magnitude);
@@ -909,7 +912,7 @@ export function normalize(out: Quat2, a: Quat2): Quat2 {
  * @param a dual quaternion to represent as a string
  * @returns string representation of the dual quat
  */
-export function str(a: Quat2): string {
+export function str(a: RQuat2): string {
     return `quat2(${a[0]}, ${a[1]}, ${a[2]}, ${a[3]}, ${a[4]}, ${a[5]}, ${a[6]}, ${a[7]})`;
 }
 
@@ -920,7 +923,7 @@ export function str(a: Quat2): string {
  * @param b the second dual quaternion.
  * @returns true if the dual quaternions are equal, false otherwise.
  */
-export function exactEquals(a: Quat2, b: Quat2): boolean {
+export function exactEquals(a: RQuat2, b: RQuat2): boolean {
     return (
         a[0] === b[0] &&
         a[1] === b[1] &&
@@ -940,7 +943,7 @@ export function exactEquals(a: Quat2, b: Quat2): boolean {
  * @param b the second dual quat.
  * @returns true if the dual quats are equal, false otherwise.
  */
-export function equals(a: Quat2, b: Quat2): boolean {
+export function equals(a: RQuat2, b: RQuat2): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
diff --git a/src/core/readonly.ts b/src/core/readonly.ts
new file mode 100644
index 0000000..0a93503
--- /dev/null
+++ b/src/core/readonly.ts
@@ -0,0 +1,2 @@
+/** Recursively makes properties and elements of plain objects, arrays, and tuples read-only. */
+export type DeepReadonly<T> = { readonly [K in keyof T]: DeepReadonly<T[K]> };
diff --git a/src/core/spherical.ts b/src/core/spherical.ts
index 3baed09..fd2312f 100644
--- a/src/core/spherical.ts
+++ b/src/core/spherical.ts
@@ -1,7 +1,7 @@
 import { wrapAngle } from './angle';
 import * as scalar from './scalar';
-import type { Vec2 } from './vec2';
-import type { Vec3 } from './vec3';
+import type { RVec2, Vec2 } from './vec2';
+import type { RVec3, Vec3 } from './vec3';
 
 /**
  * A point in spherical coordinates [r, theta, phi] (Three.js / OpenGL convention)
@@ -11,6 +11,9 @@ import type { Vec3 } from './vec3';
  */
 export type Spherical = [r: number, theta: number, phi: number];
 
+/** A read-only spherical coordinate */
+export type RSpherical = Readonly<Spherical>;
+
 /**
  * Creates a new spherical coordinate at r=1, theta=0, phi=0
  *
@@ -38,7 +41,7 @@ export function fromValues(r: number, theta: number, phi: number): Spherical {
  * @param a the source Spherical
  * @returns a new Spherical
  */
-export function clone(a: Spherical): Spherical {
+export function clone(a: RSpherical): Spherical {
     return [a[0], a[1], a[2]];
 }
 
@@ -49,7 +52,7 @@ export function clone(a: Spherical): Spherical {
  * @param a the source Spherical
  * @returns out
  */
-export function copy(out: Spherical, a: Spherical): Spherical {
+export function copy(out: Spherical, a: RSpherical): Spherical {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -79,7 +82,7 @@ export function set(out: Spherical, r: number, theta: number, phi: number): Sphe
  * @param a the source Spherical
  * @returns out
  */
-export function normalize(out: Spherical, a: Spherical): Spherical {
+export function normalize(out: Spherical, a: RSpherical): Spherical {
     out[0] = 1;
     out[1] = a[1];
     out[2] = a[2];
@@ -94,7 +97,7 @@ export function normalize(out: Spherical, a: Spherical): Spherical {
  * @param s scalar to multiply r by
  * @returns out
  */
-export function scale(out: Spherical, a: Spherical, s: number): Spherical {
+export function scale(out: Spherical, a: RSpherical, s: number): Spherical {
     out[0] = a[0] * s;
     out[1] = a[1];
     out[2] = a[2];
@@ -111,7 +114,7 @@ export function scale(out: Spherical, a: Spherical, s: number): Spherical {
  * @param t interpolation factor in [0, 1]
  * @returns out
  */
-export function lerp(out: Spherical, a: Spherical, b: Spherical, t: number): Spherical {
+export function lerp(out: Spherical, a: RSpherical, b: RSpherical, t: number): Spherical {
     out[0] = scalar.lerp(a[0], b[0], t);
     out[1] = a[1] + wrapAngle(b[1] - a[1]) * t;
     out[2] = a[2] + wrapAngle(b[2] - a[2]) * t;
@@ -128,7 +131,7 @@ export function lerp(out: Spherical, a: Spherical, b: Spherical, t: number): Sph
  * @param v the source Vec3
  * @returns out
  */
-export function setFromVec3(out: Spherical, v: Vec3): Spherical {
+export function setFromVec3(out: Spherical, v: RVec3): Spherical {
     const x = v[0];
     const y = v[1];
     const z = v[2];
@@ -151,7 +154,7 @@ export const fromVec3 = setFromVec3;
  * @param a the source Spherical
  * @returns out
  */
-export function makeSafe(out: Spherical, a: Spherical): Spherical {
+export function makeSafe(out: Spherical, a: RSpherical): Spherical {
     const EPS = scalar.EPSILON;
     out[0] = a[0];
     out[1] = a[1];
@@ -169,7 +172,7 @@ export function makeSafe(out: Spherical, a: Spherical): Spherical {
  * @param a the source Spherical
  * @returns out
  */
-export function toVec3(out: Vec3, a: Spherical): Vec3 {
+export function toVec3(out: Vec3, a: RSpherical): Vec3 {
     const r = a[0];
     const theta = a[1];
     const phi = a[2];
@@ -188,7 +191,7 @@ export function toVec3(out: Vec3, a: Spherical): Vec3 {
  * @param v the source Vec2 interpreted as (x, z)
  * @returns out
  */
-export function fromVec2(out: Spherical, v: Vec2): Spherical {
+export function fromVec2(out: Spherical, v: RVec2): Spherical {
     const x = v[0];
     const z = v[1];
     const r = Math.sqrt(x * x + z * z);
@@ -206,7 +209,7 @@ export function fromVec2(out: Spherical, v: Vec2): Spherical {
  * @param a the source Spherical
  * @returns out
  */
-export function toVec2(out: Vec2, a: Spherical): Vec2 {
+export function toVec2(out: Vec2, a: RSpherical): Vec2 {
     const r = a[0];
     const theta = a[1];
     const phi = a[2];
@@ -224,7 +227,7 @@ export function toVec2(out: Vec2, a: Spherical): Vec2 {
  * @param b the second Spherical
  * @returns true if approximately equal
  */
-export function equals(a: Spherical, b: Spherical): boolean {
+export function equals(a: RSpherical, b: RSpherical): boolean {
     return scalar.equals(a[0], b[0]) && scalar.equals(a[1], b[1]) && scalar.equals(a[2], b[2]);
 }
 
@@ -235,7 +238,7 @@ export function equals(a: Spherical, b: Spherical): boolean {
  * @param b the second Spherical
  * @returns true if exactly equal
  */
-export function exactEquals(a: Spherical, b: Spherical): boolean {
+export function exactEquals(a: RSpherical, b: RSpherical): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2];
 }
 
@@ -245,7 +248,7 @@ export function exactEquals(a: Spherical, b: Spherical): boolean {
  * @param a the source Spherical
  * @returns string representation
  */
-export function str(a: Spherical): string {
+export function str(a: RSpherical): string {
     return `Spherical(${a[0]}, ${a[1]}, ${a[2]})`;
 }
 
@@ -260,7 +263,7 @@ export function str(a: Spherical): string {
  * @param b the second Spherical
  * @returns angle in radians in [0, π]
  */
-export function angleTo(a: Spherical, b: Spherical): number {
+export function angleTo(a: RSpherical, b: RSpherical): number {
     const phiA = a[2];
     const phiB = b[2];
     const dTheta = b[1] - a[1];
diff --git a/src/core/vec2.ts b/src/core/vec2.ts
index f22b9ad..2cb2031 100644
--- a/src/core/vec2.ts
+++ b/src/core/vec2.ts
@@ -1,14 +1,17 @@
 import type { MutableArrayLike } from './arrays';
-import type { Mat2 } from './mat2';
-import type { Mat2d } from './mat2d';
-import type { Mat3 } from './mat3';
-import type { Mat4 } from './mat4';
+import type { RMat2 } from './mat2';
+import type { RMat2d } from './mat2d';
+import type { RMat3 } from './mat3';
+import type { RMat4 } from './mat4';
 import * as scalar from './scalar';
 import type { Vec3 } from './vec3';
 
 /** A 2D vector */
 export type Vec2 = [x: number, y: number];
 
+/** A read-only 2D vector */
+export type RVec2 = Readonly<Vec2>;
+
 /**
  * Creates a new, empty vec2
  *
@@ -24,7 +27,7 @@ export function create(): Vec2 {
  * @param a vector to clone
  * @returns a new 2D vector
  */
-export function clone(a: Vec2): Vec2 {
+export function clone(a: RVec2): Vec2 {
     return [a[0], a[1]];
 }
 
@@ -46,7 +49,7 @@ export function fromValues(x: number, y: number): Vec2 {
  * @param a the source vector
  * @returns out
  */
-export function copy(out: Vec2, a: Vec2): Vec2 {
+export function copy(out: Vec2, a: RVec2): Vec2 {
     out[0] = a[0];
     out[1] = a[1];
     return out;
@@ -86,7 +89,7 @@ export function fromBuffer(out: Vec2, buffer: ArrayLike<number>, startIndex: num
  * @param startIndex The starting index in the buffer
  * @returns The output buffer
  */
-export function toBuffer(outBuffer: MutableArrayLike<number>, vec: Vec2, startIndex: number): MutableArrayLike<number> {
+export function toBuffer(outBuffer: MutableArrayLike<number>, vec: RVec2, startIndex: number): MutableArrayLike<number> {
     outBuffer[startIndex] = vec[0];
     outBuffer[startIndex + 1] = vec[1];
     return outBuffer;
@@ -100,7 +103,7 @@ export function toBuffer(outBuffer: MutableArrayLike<number>, vec: Vec2, startIn
  * @param b the second operand
  * @returns out
  */
-export function add(out: Vec2, a: Vec2, b: Vec2): Vec2 {
+export function add(out: Vec2, a: RVec2, b: RVec2): Vec2 {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     return out;
@@ -114,7 +117,7 @@ export function add(out: Vec2, a: Vec2, b: Vec2): Vec2 {
  * @param b the scalar value to add
  * @returns out
  */
-export function addScalar(out: Vec2, a: Vec2, b: number): Vec2 {
+export function addScalar(out: Vec2, a: RVec2, b: number): Vec2 {
     out[0] = a[0] + b;
     out[1] = a[1] + b;
     return out;
@@ -128,7 +131,7 @@ export function addScalar(out: Vec2, a: Vec2, b: number): Vec2 {
  * @param b the second operand
  * @returns out
  */
-export function subtract(out: Vec2, a: Vec2, b: Vec2): Vec2 {
+export function subtract(out: Vec2, a: RVec2, b: RVec2): Vec2 {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     return out;
@@ -142,7 +145,7 @@ export function subtract(out: Vec2, a: Vec2, b: Vec2): Vec2 {
  * @param b the scalar value to subtract
  * @returns out
  */
-export function subtractScalar(out: Vec2, a: Vec2, b: number): Vec2 {
+export function subtractScalar(out: Vec2, a: RVec2, b: number): Vec2 {
     out[0] = a[0] - b;
     out[1] = a[1] - b;
     return out;
@@ -156,7 +159,7 @@ export function subtractScalar(out: Vec2, a: Vec2, b: number): Vec2 {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Vec2, a: Vec2, b: Vec2): Vec2 {
+export function multiply(out: Vec2, a: RVec2, b: RVec2): Vec2 {
     out[0] = a[0] * b[0];
     out[1] = a[1] * b[1];
     return out;
@@ -170,7 +173,7 @@ export function multiply(out: Vec2, a: Vec2, b: Vec2): Vec2 {
  * @param b the second operand
  * @returns out
  */
-export function divide(out: Vec2, a: Vec2, b: Vec2): Vec2 {
+export function divide(out: Vec2, a: RVec2, b: RVec2): Vec2 {
     out[0] = a[0] / b[0];
     out[1] = a[1] / b[1];
     return out;
@@ -183,7 +186,7 @@ export function divide(out: Vec2, a: Vec2, b: Vec2): Vec2 {
  * @param a vector to ceil
  * @returns out
  */
-export function ceil(out: Vec2, a: Vec2): Vec2 {
+export function ceil(out: Vec2, a: RVec2): Vec2 {
     out[0] = Math.ceil(a[0]);
     out[1] = Math.ceil(a[1]);
     return out;
@@ -196,7 +199,7 @@ export function ceil(out: Vec2, a: Vec2): Vec2 {
  * @param a vector to floor
  * @returns out
  */
-export function floor(out: Vec2, a: Vec2): Vec2 {
+export function floor(out: Vec2, a: RVec2): Vec2 {
     out[0] = Math.floor(a[0]);
     out[1] = Math.floor(a[1]);
     return out;
@@ -210,7 +213,7 @@ export function floor(out: Vec2, a: Vec2): Vec2 {
  * @param b the second operand
  * @returns out
  */
-export function min(out: Vec2, a: Vec2, b: Vec2): Vec2 {
+export function min(out: Vec2, a: RVec2, b: RVec2): Vec2 {
     out[0] = Math.min(a[0], b[0]);
     out[1] = Math.min(a[1], b[1]);
     return out;
@@ -224,7 +227,7 @@ export function min(out: Vec2, a: Vec2, b: Vec2): Vec2 {
  * @param b the second operand
  * @returns out
  */
-export function max(out: Vec2, a: Vec2, b: Vec2): Vec2 {
+export function max(out: Vec2, a: RVec2, b: RVec2): Vec2 {
     out[0] = Math.max(a[0], b[0]);
     out[1] = Math.max(a[1], b[1]);
     return out;
@@ -237,7 +240,7 @@ export function max(out: Vec2, a: Vec2, b: Vec2): Vec2 {
  * @param a vector to round
  * @returns out
  */
-export function round(out: Vec2, a: Vec2): Vec2 {
+export function round(out: Vec2, a: RVec2): Vec2 {
     out[0] = scalar.round(a[0]);
     out[1] = scalar.round(a[1]);
     return out;
@@ -251,7 +254,7 @@ export function round(out: Vec2, a: Vec2): Vec2 {
  * @param b amount to scale the vector by
  * @returns out
  */
-export function scale(out: Vec2, a: Vec2, b: number): Vec2 {
+export function scale(out: Vec2, a: RVec2, b: number): Vec2 {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     return out;
@@ -266,7 +269,7 @@ export function scale(out: Vec2, a: Vec2, b: number): Vec2 {
  * @param scale the amount to scale b by before adding
  * @returns out
  */
-export function scaleAndAdd(out: Vec2, a: Vec2, b: Vec2, scale: number): Vec2 {
+export function scaleAndAdd(out: Vec2, a: RVec2, b: RVec2, scale: number): Vec2 {
     out[0] = a[0] + b[0] * scale;
     out[1] = a[1] + b[1] * scale;
     return out;
@@ -279,7 +282,7 @@ export function scaleAndAdd(out: Vec2, a: Vec2, b: Vec2, scale: number): Vec2 {
  * @param b the second operand
  * @returns distance between a and b
  */
-export function distance(a: Vec2, b: Vec2): number {
+export function distance(a: RVec2, b: RVec2): number {
     const x = b[0] - a[0];
     const y = b[1] - a[1];
     return Math.sqrt(x * x + y * y);
@@ -292,7 +295,7 @@ export function distance(a: Vec2, b: Vec2): number {
  * @param b the second operand
  * @returns squared distance between a and b
  */
-export function squaredDistance(a: Vec2, b: Vec2): number {
+export function squaredDistance(a: RVec2, b: RVec2): number {
     const x = b[0] - a[0];
     const y = b[1] - a[1];
     return x * x + y * y;
@@ -304,7 +307,7 @@ export function squaredDistance(a: Vec2, b: Vec2): number {
  * @param a vector to calculate length of
  * @returns length of a
  */
-export function length(a: Vec2): number {
+export function length(a: RVec2): number {
     const x = a[0];
     const y = a[1];
     return Math.sqrt(x * x + y * y);
@@ -316,7 +319,7 @@ export function length(a: Vec2): number {
  * @param a vector to calculate squared length of
  * @returns squared length of a
  */
-export function squaredLength(a: Vec2): number {
+export function squaredLength(a: RVec2): number {
     const x = a[0];
     const y = a[1];
     return x * x + y * y;
@@ -329,7 +332,7 @@ export function squaredLength(a: Vec2): number {
  * @param a vector to negate
  * @returns out
  */
-export function negate(out: Vec2, a: Vec2): Vec2 {
+export function negate(out: Vec2, a: RVec2): Vec2 {
     out[0] = -a[0];
     out[1] = -a[1];
     return out;
@@ -342,7 +345,7 @@ export function negate(out: Vec2, a: Vec2): Vec2 {
  * @param a vector to invert
  * @returns out
  */
-export function inverse(out: Vec2, a: Vec2): Vec2 {
+export function inverse(out: Vec2, a: RVec2): Vec2 {
     out[0] = 1.0 / a[0];
     out[1] = 1.0 / a[1];
     return out;
@@ -355,7 +358,7 @@ export function inverse(out: Vec2, a: Vec2): Vec2 {
  * @param a vector to normalize
  * @returns out
  */
-export function normalize(out: Vec2, a: Vec2): Vec2 {
+export function normalize(out: Vec2, a: RVec2): Vec2 {
     const x = a[0];
     const y = a[1];
     let len = x * x + y * y;
@@ -374,7 +377,7 @@ export function normalize(out: Vec2, a: Vec2): Vec2 {
  * @param b the second operand
  * @returns dot product of a and b
  */
-export function dot(a: Vec2, b: Vec2): number {
+export function dot(a: RVec2, b: RVec2): number {
     return a[0] * b[0] + a[1] * b[1];
 }
 
@@ -387,7 +390,7 @@ export function dot(a: Vec2, b: Vec2): number {
  * @param b the second operand
  * @returns out
  */
-export function cross(out: Vec3, a: Vec2, b: Vec2): Vec3 {
+export function cross(out: Vec3, a: RVec2, b: RVec2): Vec3 {
     const z = a[0] * b[1] - a[1] * b[0];
     out[0] = out[1] = 0;
     out[2] = z;
@@ -403,7 +406,7 @@ export function cross(out: Vec3, a: Vec2, b: Vec2): Vec3 {
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function lerp(out: Vec2, a: Vec2, b: Vec2, t: number): Vec2 {
+export function lerp(out: Vec2, a: RVec2, b: RVec2, t: number): Vec2 {
     const ax = a[0];
     const ay = a[1];
     out[0] = ax + t * (b[0] - ax);
@@ -422,7 +425,7 @@ export function lerp(out: Vec2, a: Vec2, b: Vec2, t: number): Vec2 {
  * @param t interpolation amount
  * @returns out
  */
-export function lagrange(out: Vec2, a: Vec2, b: Vec2, c: Vec2, t: number): Vec2 {
+export function lagrange(out: Vec2, a: RVec2, b: RVec2, c: RVec2, t: number): Vec2 {
     const c0 = 2 * (t - 1) * (t - 0.5);
     const c1 = -4 * (t - 1) * t;
     const c2 = 2 * (t - 0.5) * t;
@@ -439,7 +442,7 @@ export function lagrange(out: Vec2, a: Vec2, b: Vec2, c: Vec2, t: number): Vec2
  * @param m matrix to transform with
  * @returns out
  */
-export function transformMat2(out: Vec2, a: Vec2, m: Mat2): Vec2 {
+export function transformMat2(out: Vec2, a: RVec2, m: RMat2): Vec2 {
     const x = a[0];
     const y = a[1];
     out[0] = m[0] * x + m[2] * y;
@@ -455,7 +458,7 @@ export function transformMat2(out: Vec2, a: Vec2, m: Mat2): Vec2 {
  * @param m matrix to transform with
  * @returns out
  */
-export function transformMat2d(out: Vec2, a: Vec2, m: Mat2d): Vec2 {
+export function transformMat2d(out: Vec2, a: RVec2, m: RMat2d): Vec2 {
     const x = a[0];
     const y = a[1];
     out[0] = m[0] * x + m[2] * y + m[4];
@@ -472,7 +475,7 @@ export function transformMat2d(out: Vec2, a: Vec2, m: Mat2d): Vec2 {
  * @param m matrix to transform with
  * @returns out
  */
-export function transformMat3(out: Vec2, a: Vec2, m: Mat3): Vec2 {
+export function transformMat3(out: Vec2, a: RVec2, m: RMat3): Vec2 {
     const x = a[0];
     const y = a[1];
     out[0] = m[0] * x + m[3] * y + m[6];
@@ -490,7 +493,7 @@ export function transformMat3(out: Vec2, a: Vec2, m: Mat3): Vec2 {
  * @param m matrix to transform with
  * @returns out
  */
-export function transformMat4(out: Vec2, a: Vec2, m: Mat4): Vec2 {
+export function transformMat4(out: Vec2, a: RVec2, m: RMat4): Vec2 {
     const x = a[0];
     const y = a[1];
     out[0] = m[0] * x + m[4] * y + m[12];
@@ -506,7 +509,7 @@ export function transformMat4(out: Vec2, a: Vec2, m: Mat4): Vec2 {
  * @param rad The angle of rotation in radians
  * @returns out
  */
-export function rotate(out: Vec2, a: Vec2, b: Vec2, rad: number): Vec2 {
+export function rotate(out: Vec2, a: RVec2, b: RVec2, rad: number): Vec2 {
     //Translate point to the origin
     const p0 = a[0] - b[0];
     const p1 = a[1] - b[1];
@@ -526,7 +529,7 @@ export function rotate(out: Vec2, a: Vec2, b: Vec2, rad: number): Vec2 {
  * @param b The second operand
  * @returns The angle in radians
  */
-export function angle(a: Vec2, b: Vec2): number {
+export function angle(a: RVec2, b: RVec2): number {
     const x1 = a[0];
     const y1 = a[1];
     const x2 = b[0];
@@ -551,7 +554,7 @@ export function angle(a: Vec2, b: Vec2): number {
  * @param b the second operand
  * @returns the signed angle in radians
  */
-export function signedAngle(a: Vec2, b: Vec2): number {
+export function signedAngle(a: RVec2, b: RVec2): number {
     const ax = a[0];
     const ay = a[1];
     const bx = b[0];
@@ -579,7 +582,7 @@ export function zero(out: Vec2): Vec2 {
  * @param a vector to represent as a string
  * @returns string representation of the vector
  */
-export function str(a: Vec2): string {
+export function str(a: RVec2): string {
     return `vec2(${a[0]}, ${a[1]})`;
 }
 
@@ -590,7 +593,7 @@ export function str(a: Vec2): string {
  * @param b The second vector.
  * @returns True if the vectors are equal, false otherwise.
  */
-export function exactEquals(a: Vec2, b: Vec2): boolean {
+export function exactEquals(a: RVec2, b: RVec2): boolean {
     return a[0] === b[0] && a[1] === b[1];
 }
 
@@ -601,7 +604,7 @@ export function exactEquals(a: Vec2, b: Vec2): boolean {
  * @param b The second vector.
  * @returns True if the vectors are equal, false otherwise.
  */
-export function equals(a: Vec2, b: Vec2): boolean {
+export function equals(a: RVec2, b: RVec2): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const b0 = b[0];
@@ -617,7 +620,7 @@ export function equals(a: Vec2, b: Vec2): boolean {
  * @param a vector to test
  * @returns whether or not the vector is finite
  */
-export function finite(a: Vec2): boolean {
+export function finite(a: RVec2): boolean {
     return Number.isFinite(a[0]) && Number.isFinite(a[1]);
 }
 
diff --git a/src/core/vec3.ts b/src/core/vec3.ts
index efbf24e..1f20479 100644
--- a/src/core/vec3.ts
+++ b/src/core/vec3.ts
@@ -1,12 +1,15 @@
 import type { MutableArrayLike } from './arrays';
-import type { Mat3 } from './mat3';
-import type { Mat4 } from './mat4';
-import type { Quat } from './quat';
+import type { RMat3 } from './mat3';
+import type { RMat4 } from './mat4';
+import type { RQuat } from './quat';
 import * as scalar from './scalar';
 
 /** A 3D vector */
 export type Vec3 = [x: number, y: number, z: number];
 
+/** A read-only 3D vector */
+export type RVec3 = Readonly<Vec3>;
+
 /**
  * Creates a new, empty vec3
  *
@@ -22,7 +25,7 @@ export function create(): Vec3 {
  * @param a vector to clone
  * @returns a new 3D vector
  */
-export function clone(a: Vec3): Vec3 {
+export function clone(a: RVec3): Vec3 {
     return [a[0], a[1], a[2]];
 }
 
@@ -44,7 +47,7 @@ export function fromValues(x: number, y: number, z: number): Vec3 {
  * @param a vector to calculate length of
  * @returns length of a
  */
-export function length(a: Vec3): number {
+export function length(a: RVec3): number {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -58,7 +61,7 @@ export function length(a: Vec3): number {
  * @param a the source vector
  * @returns out
  */
-export function copy(out: Vec3, a: Vec3): Vec3 {
+export function copy(out: Vec3, a: RVec3): Vec3 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -116,7 +119,7 @@ export function fromBuffer(out: Vec3, buffer: ArrayLike<number>, startIndex: num
  * @param startIndex The starting index in the buffer
  * @returns The output buffer
  */
-export function toBuffer(outBuffer: MutableArrayLike<number>, vec: Vec3, startIndex: number): MutableArrayLike<number> {
+export function toBuffer(outBuffer: MutableArrayLike<number>, vec: RVec3, startIndex: number): MutableArrayLike<number> {
     outBuffer[startIndex] = vec[0];
     outBuffer[startIndex + 1] = vec[1];
     outBuffer[startIndex + 2] = vec[2];
@@ -131,7 +134,7 @@ export function toBuffer(outBuffer: MutableArrayLike<number>, vec: Vec3, startIn
  * @param b the second operand
  * @returns out
  */
-export function add(out: Vec3, a: Vec3, b: Vec3): Vec3 {
+export function add(out: Vec3, a: RVec3, b: RVec3): Vec3 {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -146,7 +149,7 @@ export function add(out: Vec3, a: Vec3, b: Vec3): Vec3 {
  * @param b the scalar value to add
  * @returns out
  */
-export function addScalar(out: Vec3, a: Vec3, b: number): Vec3 {
+export function addScalar(out: Vec3, a: RVec3, b: number): Vec3 {
     out[0] = a[0] + b;
     out[1] = a[1] + b;
     out[2] = a[2] + b;
@@ -161,7 +164,7 @@ export function addScalar(out: Vec3, a: Vec3, b: number): Vec3 {
  * @param b the second operand
  * @returns out
  */
-export function subtract(out: Vec3, a: Vec3, b: Vec3): Vec3 {
+export function subtract(out: Vec3, a: RVec3, b: RVec3): Vec3 {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     out[2] = a[2] - b[2];
@@ -176,7 +179,7 @@ export function subtract(out: Vec3, a: Vec3, b: Vec3): Vec3 {
  * @param b the scalar value to subtract
  * @returns out
  */
-export function subtractScalar(out: Vec3, a: Vec3, b: number): Vec3 {
+export function subtractScalar(out: Vec3, a: RVec3, b: number): Vec3 {
     out[0] = a[0] - b;
     out[1] = a[1] - b;
     out[2] = a[2] - b;
@@ -190,7 +193,7 @@ export function subtractScalar(out: Vec3, a: Vec3, b: number): Vec3 {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Vec3, a: Vec3, b: Vec3): Vec3 {
+export function multiply(out: Vec3, a: RVec3, b: RVec3): Vec3 {
     out[0] = a[0] * b[0];
     out[1] = a[1] * b[1];
     out[2] = a[2] * b[2];
@@ -205,7 +208,7 @@ export function multiply(out: Vec3, a: Vec3, b: Vec3): Vec3 {
  * @param b the second operand
  * @returns out
  */
-export function divide(out: Vec3, a: Vec3, b: Vec3): Vec3 {
+export function divide(out: Vec3, a: RVec3, b: RVec3): Vec3 {
     out[0] = a[0] / b[0];
     out[1] = a[1] / b[1];
     out[2] = a[2] / b[2];
@@ -219,7 +222,7 @@ export function divide(out: Vec3, a: Vec3, b: Vec3): Vec3 {
  * @param a vector to ceil
  * @returns out
  */
-export function ceil(out: Vec3, a: Vec3): Vec3 {
+export function ceil(out: Vec3, a: RVec3): Vec3 {
     out[0] = Math.ceil(a[0]);
     out[1] = Math.ceil(a[1]);
     out[2] = Math.ceil(a[2]);
@@ -233,7 +236,7 @@ export function ceil(out: Vec3, a: Vec3): Vec3 {
  * @param a vector to floor
  * @returns out
  */
-export function floor(out: Vec3, a: Vec3): Vec3 {
+export function floor(out: Vec3, a: RVec3): Vec3 {
     out[0] = Math.floor(a[0]);
     out[1] = Math.floor(a[1]);
     out[2] = Math.floor(a[2]);
@@ -248,7 +251,7 @@ export function floor(out: Vec3, a: Vec3): Vec3 {
  * @param b the second operand
  * @returns out
  */
-export function min(out: Vec3, a: Vec3, b: Vec3): Vec3 {
+export function min(out: Vec3, a: RVec3, b: RVec3): Vec3 {
     out[0] = Math.min(a[0], b[0]);
     out[1] = Math.min(a[1], b[1]);
     out[2] = Math.min(a[2], b[2]);
@@ -263,7 +266,7 @@ export function min(out: Vec3, a: Vec3, b: Vec3): Vec3 {
  * @param b the second operand
  * @returns out
  */
-export function max(out: Vec3, a: Vec3, b: Vec3): Vec3 {
+export function max(out: Vec3, a: RVec3, b: RVec3): Vec3 {
     out[0] = Math.max(a[0], b[0]);
     out[1] = Math.max(a[1], b[1]);
     out[2] = Math.max(a[2], b[2]);
@@ -277,7 +280,7 @@ export function max(out: Vec3, a: Vec3, b: Vec3): Vec3 {
  * @param a vector to round
  * @returns out
  */
-export function round(out: Vec3, a: Vec3): Vec3 {
+export function round(out: Vec3, a: RVec3): Vec3 {
     out[0] = scalar.round(a[0]);
     out[1] = scalar.round(a[1]);
     out[2] = scalar.round(a[2]);
@@ -292,7 +295,7 @@ export function round(out: Vec3, a: Vec3): Vec3 {
  * @param b amount to scale the vector by
  * @returns out
  */
-export function scale(out: Vec3, a: Vec3, b: number): Vec3 {
+export function scale(out: Vec3, a: RVec3, b: number): Vec3 {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     out[2] = a[2] * b;
@@ -308,7 +311,7 @@ export function scale(out: Vec3, a: Vec3, b: number): Vec3 {
  * @param scale the amount to scale b by before adding
  * @returns out
  */
-export function scaleAndAdd(out: Vec3, a: Vec3, b: Vec3, scale: number): Vec3 {
+export function scaleAndAdd(out: Vec3, a: RVec3, b: RVec3, scale: number): Vec3 {
     out[0] = a[0] + b[0] * scale;
     out[1] = a[1] + b[1] * scale;
     out[2] = a[2] + b[2] * scale;
@@ -322,7 +325,7 @@ export function scaleAndAdd(out: Vec3, a: Vec3, b: Vec3, scale: number): Vec3 {
  * @param b the second operand
  * @returns distance between a and b
  */
-export function distance(a: Vec3, b: Vec3): number {
+export function distance(a: RVec3, b: RVec3): number {
     const x = b[0] - a[0];
     const y = b[1] - a[1];
     const z = b[2] - a[2];
@@ -336,7 +339,7 @@ export function distance(a: Vec3, b: Vec3): number {
  * @param b the second operand
  * @returns squared distance between a and b
  */
-export function squaredDistance(a: Vec3, b: Vec3): number {
+export function squaredDistance(a: RVec3, b: RVec3): number {
     const x = b[0] - a[0];
     const y = b[1] - a[1];
     const z = b[2] - a[2];
@@ -349,7 +352,7 @@ export function squaredDistance(a: Vec3, b: Vec3): number {
  * @param a vector to calculate squared length of
  * @returns squared length of a
  */
-export function squaredLength(a: Vec3): number {
+export function squaredLength(a: RVec3): number {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -363,7 +366,7 @@ export function squaredLength(a: Vec3): number {
  * @param a vector to negate
  * @returns out
  */
-export function negate(out: Vec3, a: Vec3): Vec3 {
+export function negate(out: Vec3, a: RVec3): Vec3 {
     out[0] = -a[0];
     out[1] = -a[1];
     out[2] = -a[2];
@@ -377,7 +380,7 @@ export function negate(out: Vec3, a: Vec3): Vec3 {
  * @param a vector to invert
  * @returns out
  */
-export function inverse(out: Vec3, a: Vec3): Vec3 {
+export function inverse(out: Vec3, a: RVec3): Vec3 {
     out[0] = 1.0 / a[0];
     out[1] = 1.0 / a[1];
     out[2] = 1.0 / a[2];
@@ -391,7 +394,7 @@ export function inverse(out: Vec3, a: Vec3): Vec3 {
  * @param a vector to normalize
  * @returns out
  */
-export function normalize(out: Vec3, a: Vec3): Vec3 {
+export function normalize(out: Vec3, a: RVec3): Vec3 {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -412,7 +415,7 @@ export function normalize(out: Vec3, a: Vec3): Vec3 {
  * @param b the second operand
  * @returns dot product of a and b
  */
-export function dot(a: Vec3, b: Vec3): number {
+export function dot(a: RVec3, b: RVec3): number {
     return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
 }
 
@@ -424,7 +427,7 @@ export function dot(a: Vec3, b: Vec3): number {
  * @param b the second operand
  * @returns out
  */
-export function cross(out: Vec3, a: Vec3, b: Vec3): Vec3 {
+export function cross(out: Vec3, a: RVec3, b: RVec3): Vec3 {
     const ax = a[0];
     const ay = a[1];
     const az = a[2];
@@ -449,7 +452,7 @@ export function cross(out: Vec3, a: Vec3, b: Vec3): Vec3 {
  * @param a the source vector
  * @returns the out vector
  */
-export function perpendicular(out: Vec3, a: Vec3): Vec3 {
+export function perpendicular(out: Vec3, a: RVec3): Vec3 {
     if (Math.abs(a[0]) > Math.abs(a[1])) {
         const len = Math.sqrt(a[0] * a[0] + a[2] * a[2]);
         const invLen = 1.0 / len;
@@ -485,7 +488,7 @@ export function perpendicular(out: Vec3, a: Vec3): Vec3 {
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function lerp(out: Vec3, a: Vec3, b: Vec3, t: number): Vec3 {
+export function lerp(out: Vec3, a: RVec3, b: RVec3, t: number): Vec3 {
     const ax = a[0];
     const ay = a[1];
     const az = a[2];
@@ -506,7 +509,7 @@ export function lerp(out: Vec3, a: Vec3, b: Vec3, t: number): Vec3 {
  * @param t interpolation amount
  * @returns out
  */
-export function lagrange(out: Vec3, a: Vec3, b: Vec3, c: Vec3, t: number): Vec3 {
+export function lagrange(out: Vec3, a: RVec3, b: RVec3, c: RVec3, t: number): Vec3 {
     const c0 = 2 * (t - 1) * (t - 0.5);
     const c1 = -4 * (t - 1) * t;
     const c2 = 2 * (t - 0.5) * t;
@@ -525,7 +528,7 @@ export function lagrange(out: Vec3, a: Vec3, b: Vec3, c: Vec3, t: number): Vec3
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function slerp(out: Vec3, a: Vec3, b: Vec3, t: number): Vec3 {
+export function slerp(out: Vec3, a: RVec3, b: RVec3, t: number): Vec3 {
     const angle = Math.acos(Math.min(Math.max(dot(a, b), -1), 1));
     const sinTotal = Math.sin(angle);
 
@@ -549,7 +552,7 @@ export function slerp(out: Vec3, a: Vec3, b: Vec3, t: number): Vec3 {
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function hermite(out: Vec3, a: Vec3, b: Vec3, c: Vec3, d: Vec3, t: number): Vec3 {
+export function hermite(out: Vec3, a: RVec3, b: RVec3, c: RVec3, d: RVec3, t: number): Vec3 {
     const factorTimes2 = t * t;
     const factor1 = factorTimes2 * (2 * t - 3) + 1;
     const factor2 = factorTimes2 * (t - 2) + t;
@@ -574,7 +577,7 @@ export function hermite(out: Vec3, a: Vec3, b: Vec3, c: Vec3, d: Vec3, t: number
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function bezier(out: Vec3, a: Vec3, b: Vec3, c: Vec3, d: Vec3, t: number): Vec3 {
+export function bezier(out: Vec3, a: RVec3, b: RVec3, c: RVec3, d: RVec3, t: number): Vec3 {
     const inverseFactor = 1 - t;
     const inverseFactorTimesTwo = inverseFactor * inverseFactor;
     const factorTimes2 = t * t;
@@ -599,7 +602,7 @@ export function bezier(out: Vec3, a: Vec3, b: Vec3, c: Vec3, d: Vec3, t: number)
  * @param m matrix to transform with
  * @returns out
  */
-export function transformMat4(out: Vec3, a: Vec3, m: Mat4): Vec3 {
+export function transformMat4(out: Vec3, a: RVec3, m: RMat4): Vec3 {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -619,7 +622,7 @@ export function transformMat4(out: Vec3, a: Vec3, m: Mat4): Vec3 {
  * @param m the 3x3 matrix to transform with
  * @returns out
  */
-export function transformMat3(out: Vec3, a: Vec3, m: Mat3): Vec3 {
+export function transformMat3(out: Vec3, a: RVec3, m: RMat3): Vec3 {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -638,7 +641,7 @@ export function transformMat3(out: Vec3, a: Vec3, m: Mat3): Vec3 {
  * @param q quaternion to transform with
  * @returns out
  */
-export function transformQuat(out: Vec3, a: Vec3, q: Quat): Vec3 {
+export function transformQuat(out: Vec3, a: RVec3, q: RQuat): Vec3 {
     // benchmarks: https://jsperf.com/quaternion-transform-vec3-implementations-fixed
     const qx = q[0];
     const qy = q[1];
@@ -680,7 +683,7 @@ export function transformQuat(out: Vec3, a: Vec3, q: Quat): Vec3 {
  * @param rad The angle of rotation in radians
  * @returns out
  */
-export function rotateX(out: Vec3, a: Vec3, b: Vec3, rad: number): Vec3 {
+export function rotateX(out: Vec3, a: RVec3, b: RVec3, rad: number): Vec3 {
     const p: number[] = [];
     const r: number[] = [];
     //Translate point to the origin
@@ -709,7 +712,7 @@ export function rotateX(out: Vec3, a: Vec3, b: Vec3, rad: number): Vec3 {
  * @param rad The angle of rotation in radians
  * @returns out
  */
-export function rotateY(out: Vec3, a: Vec3, b: Vec3, rad: number): Vec3 {
+export function rotateY(out: Vec3, a: RVec3, b: RVec3, rad: number): Vec3 {
     const p: number[] = [];
     const r: number[] = [];
 
@@ -739,7 +742,7 @@ export function rotateY(out: Vec3, a: Vec3, b: Vec3, rad: number): Vec3 {
  * @param rad The angle of rotation in radians
  * @returns out
  */
-export function rotateZ(out: Vec3, a: Vec3, b: Vec3, rad: number): Vec3 {
+export function rotateZ(out: Vec3, a: RVec3, b: RVec3, rad: number): Vec3 {
     const p: number[] = [];
     const r: number[] = [];
     // translate point to the origin
@@ -766,7 +769,7 @@ export function rotateZ(out: Vec3, a: Vec3, b: Vec3, rad: number): Vec3 {
  * @param b The second operand
  * @returns The angle in radians
  */
-export function angle(a: Vec3, b: Vec3): number {
+export function angle(a: RVec3, b: RVec3): number {
     const ax = a[0];
     const ay = a[1];
     const az = a[2];
@@ -793,7 +796,7 @@ export function angle(a: Vec3, b: Vec3): number {
  * @param axis the axis to measure the rotation about, assumed to be unit length
  * @returns the signed angle in radians
  */
-export function signedAngle(a: Vec3, b: Vec3, axis: Vec3): number {
+export function signedAngle(a: RVec3, b: RVec3, axis: RVec3): number {
     const nx = axis[0];
     const ny = axis[1];
     const nz = axis[2];
@@ -836,7 +839,7 @@ const _rotateTowards_axis: Vec3 = [0, 0, 0];
  * @param maxAngle the maximum rotation, in radians
  * @returns out
  */
-export function rotateTowards(out: Vec3, from: Vec3, to: Vec3, maxAngle: number): Vec3 {
+export function rotateTowards(out: Vec3, from: RVec3, to: RVec3, maxAngle: number): Vec3 {
     const fx = from[0];
     const fy = from[1];
     const fz = from[2];
@@ -905,7 +908,7 @@ export function zero(out: Vec3): Vec3 {
  * @param a vector to represent as a string
  * @returns string representation of the vector
  */
-export function str(a: Vec3): string {
+export function str(a: RVec3): string {
     return `vec3(${a[0]}, ${a[1]}, ${a[2]})`;
 }
 
@@ -916,7 +919,7 @@ export function str(a: Vec3): string {
  * @param b The second vector.
  * @returns True if the vectors are equal, false otherwise.
  */
-export function exactEquals(a: Vec3, b: Vec3): boolean {
+export function exactEquals(a: RVec3, b: RVec3): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2];
 }
 
@@ -927,7 +930,7 @@ export function exactEquals(a: Vec3, b: Vec3): boolean {
  * @param b The second vector.
  * @returns True if the vectors are equal, false otherwise.
  */
-export function equals(a: Vec3, b: Vec3): boolean {
+export function equals(a: RVec3, b: RVec3): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -946,7 +949,7 @@ export function equals(a: Vec3, b: Vec3): boolean {
  * @param a vector to test
  * @returns whether or not the vector is finite
  */
-export function finite(a: Vec3): boolean {
+export function finite(a: RVec3): boolean {
     return Number.isFinite(a[0]) && Number.isFinite(a[1]) && Number.isFinite(a[2]);
 }
 
@@ -957,7 +960,7 @@ export function finite(a: Vec3): boolean {
  * @param scale The scale vector to test
  * @returns true if the scale represents a reflection (odd number of negative components)
  */
-export function isScaleInsideOut(scale: Vec3): boolean {
+export function isScaleInsideOut(scale: RVec3): boolean {
     // create a bitmask of which components are negative
     // each component that is < 0 contributes a bit (1, 2, or 4)
     const mask = (scale[0] < 0 ? 1 : 0) | (scale[1] < 0 ? 2 : 0) | (scale[2] < 0 ? 4 : 0);
diff --git a/src/core/vec4.ts b/src/core/vec4.ts
index 9ccc4ff..7de79b3 100644
--- a/src/core/vec4.ts
+++ b/src/core/vec4.ts
@@ -1,11 +1,14 @@
 import type { MutableArrayLike } from './arrays';
-import type { Mat4 } from './mat4';
-import type { Quat } from './quat';
+import type { RMat4 } from './mat4';
+import type { RQuat } from './quat';
 import * as scalar from './scalar';
 
 /** A 4D vector */
 export type Vec4 = [x: number, y: number, z: number, w: number];
 
+/** A read-only 4D vector */
+export type RVec4 = Readonly<Vec4>;
+
 /**
  * Creates a new, empty vec4
  *
@@ -21,7 +24,7 @@ export function create(): Vec4 {
  * @param a vector to clone
  * @returns a new 4D vector
  */
-export function clone(a: Vec4): Vec4 {
+export function clone(a: RVec4): Vec4 {
     return [a[0], a[1], a[2], a[3]];
 }
 
@@ -45,7 +48,7 @@ export function fromValues(x: number, y: number, z: number, w: number): Vec4 {
  * @param a the source vector
  * @returns out
  */
-export function copy(out: Vec4, a: Vec4): Vec4 {
+export function copy(out: Vec4, a: RVec4): Vec4 {
     out[0] = a[0];
     out[1] = a[1];
     out[2] = a[2];
@@ -93,7 +96,7 @@ export function fromBuffer(out: Vec4, buffer: ArrayLike<number>, startIndex: num
  * @param startIndex The starting index in the buffer
  * @returns The output buffer
  */
-export function toBuffer(outBuffer: MutableArrayLike<number>, vec: Vec4, startIndex: number): MutableArrayLike<number> {
+export function toBuffer(outBuffer: MutableArrayLike<number>, vec: RVec4, startIndex: number): MutableArrayLike<number> {
     outBuffer[startIndex] = vec[0];
     outBuffer[startIndex + 1] = vec[1];
     outBuffer[startIndex + 2] = vec[2];
@@ -109,7 +112,7 @@ export function toBuffer(outBuffer: MutableArrayLike<number>, vec: Vec4, startIn
  * @param b the second operand
  * @returns out
  */
-export function add(out: Vec4, a: Vec4, b: Vec4): Vec4 {
+export function add(out: Vec4, a: RVec4, b: RVec4): Vec4 {
     out[0] = a[0] + b[0];
     out[1] = a[1] + b[1];
     out[2] = a[2] + b[2];
@@ -125,7 +128,7 @@ export function add(out: Vec4, a: Vec4, b: Vec4): Vec4 {
  * @param b the second operand
  * @returns out
  */
-export function subtract(out: Vec4, a: Vec4, b: Vec4): Vec4 {
+export function subtract(out: Vec4, a: RVec4, b: RVec4): Vec4 {
     out[0] = a[0] - b[0];
     out[1] = a[1] - b[1];
     out[2] = a[2] - b[2];
@@ -141,7 +144,7 @@ export function subtract(out: Vec4, a: Vec4, b: Vec4): Vec4 {
  * @param b the second operand
  * @returns out
  */
-export function multiply(out: Vec4, a: Vec4, b: Vec4): Vec4 {
+export function multiply(out: Vec4, a: RVec4, b: RVec4): Vec4 {
     out[0] = a[0] * b[0];
     out[1] = a[1] * b[1];
     out[2] = a[2] * b[2];
@@ -157,7 +160,7 @@ export function multiply(out: Vec4, a: Vec4, b: Vec4): Vec4 {
  * @param b the second operand
  * @returns out
  */
-export function divide(out: Vec4, a: Vec4, b: Vec4): Vec4 {
+export function divide(out: Vec4, a: RVec4, b: RVec4): Vec4 {
     out[0] = a[0] / b[0];
     out[1] = a[1] / b[1];
     out[2] = a[2] / b[2];
@@ -172,7 +175,7 @@ export function divide(out: Vec4, a: Vec4, b: Vec4): Vec4 {
  * @param a vector to ceil
  * @returns out
  */
-export function ceil(out: Vec4, a: Vec4): Vec4 {
+export function ceil(out: Vec4, a: RVec4): Vec4 {
     out[0] = Math.ceil(a[0]);
     out[1] = Math.ceil(a[1]);
     out[2] = Math.ceil(a[2]);
@@ -187,7 +190,7 @@ export function ceil(out: Vec4, a: Vec4): Vec4 {
  * @param a vector to floor
  * @returns out
  */
-export function floor(out: Vec4, a: Vec4): Vec4 {
+export function floor(out: Vec4, a: RVec4): Vec4 {
     out[0] = Math.floor(a[0]);
     out[1] = Math.floor(a[1]);
     out[2] = Math.floor(a[2]);
@@ -203,7 +206,7 @@ export function floor(out: Vec4, a: Vec4): Vec4 {
  * @param b the second operand
  * @returns out
  */
-export function min(out: Vec4, a: Vec4, b: Vec4): Vec4 {
+export function min(out: Vec4, a: RVec4, b: RVec4): Vec4 {
     out[0] = Math.min(a[0], b[0]);
     out[1] = Math.min(a[1], b[1]);
     out[2] = Math.min(a[2], b[2]);
@@ -219,7 +222,7 @@ export function min(out: Vec4, a: Vec4, b: Vec4): Vec4 {
  * @param b the second operand
  * @returns out
  */
-export function max(out: Vec4, a: Vec4, b: Vec4): Vec4 {
+export function max(out: Vec4, a: RVec4, b: RVec4): Vec4 {
     out[0] = Math.max(a[0], b[0]);
     out[1] = Math.max(a[1], b[1]);
     out[2] = Math.max(a[2], b[2]);
@@ -234,7 +237,7 @@ export function max(out: Vec4, a: Vec4, b: Vec4): Vec4 {
  * @param a vector to round
  * @returns out
  */
-export function round(out: Vec4, a: Vec4): Vec4 {
+export function round(out: Vec4, a: RVec4): Vec4 {
     out[0] = scalar.round(a[0]);
     out[1] = scalar.round(a[1]);
     out[2] = scalar.round(a[2]);
@@ -250,7 +253,7 @@ export function round(out: Vec4, a: Vec4): Vec4 {
  * @param b amount to scale the vector by
  * @returns out
  */
-export function scale(out: Vec4, a: Vec4, b: number): Vec4 {
+export function scale(out: Vec4, a: RVec4, b: number): Vec4 {
     out[0] = a[0] * b;
     out[1] = a[1] * b;
     out[2] = a[2] * b;
@@ -267,7 +270,7 @@ export function scale(out: Vec4, a: Vec4, b: number): Vec4 {
  * @param scale the amount to scale b by before adding
  * @returns out
  */
-export function scaleAndAdd(out: Vec4, a: Vec4, b: Vec4, scale: number): Vec4 {
+export function scaleAndAdd(out: Vec4, a: RVec4, b: RVec4, scale: number): Vec4 {
     out[0] = a[0] + b[0] * scale;
     out[1] = a[1] + b[1] * scale;
     out[2] = a[2] + b[2] * scale;
@@ -282,7 +285,7 @@ export function scaleAndAdd(out: Vec4, a: Vec4, b: Vec4, scale: number): Vec4 {
  * @param b the second operand
  * @returns distance between a and b
  */
-export function distance(a: Vec4, b: Vec4): number {
+export function distance(a: RVec4, b: RVec4): number {
     const x = b[0] - a[0];
     const y = b[1] - a[1];
     const z = b[2] - a[2];
@@ -297,7 +300,7 @@ export function distance(a: Vec4, b: Vec4): number {
  * @param b the second operand
  * @returns squared distance between a and b
  */
-export function squaredDistance(a: Vec4, b: Vec4): number {
+export function squaredDistance(a: RVec4, b: RVec4): number {
     const x = b[0] - a[0];
     const y = b[1] - a[1];
     const z = b[2] - a[2];
@@ -311,7 +314,7 @@ export function squaredDistance(a: Vec4, b: Vec4): number {
  * @param a vector to calculate length of
  * @returns length of a
  */
-export function length(a: Vec4): number {
+export function length(a: RVec4): number {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -325,7 +328,7 @@ export function length(a: Vec4): number {
  * @param a vector to calculate squared length of
  * @returns squared length of a
  */
-export function squaredLength(a: Vec4): number {
+export function squaredLength(a: RVec4): number {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -340,7 +343,7 @@ export function squaredLength(a: Vec4): number {
  * @param a vector to negate
  * @returns out
  */
-export function negate(out: Vec4, a: Vec4): Vec4 {
+export function negate(out: Vec4, a: RVec4): Vec4 {
     out[0] = -a[0];
     out[1] = -a[1];
     out[2] = -a[2];
@@ -355,7 +358,7 @@ export function negate(out: Vec4, a: Vec4): Vec4 {
  * @param a vector to invert
  * @returns out
  */
-export function inverse(out: Vec4, a: Vec4): Vec4 {
+export function inverse(out: Vec4, a: RVec4): Vec4 {
     out[0] = 1.0 / a[0];
     out[1] = 1.0 / a[1];
     out[2] = 1.0 / a[2];
@@ -370,7 +373,7 @@ export function inverse(out: Vec4, a: Vec4): Vec4 {
  * @param a vector to normalize
  * @returns out
  */
-export function normalize(out: Vec4, a: Vec4): Vec4 {
+export function normalize(out: Vec4, a: RVec4): Vec4 {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -393,7 +396,7 @@ export function normalize(out: Vec4, a: Vec4): Vec4 {
  * @param b the second operand
  * @returns dot product of a and b
  */
-export function dot(a: Vec4, b: Vec4): number {
+export function dot(a: RVec4, b: RVec4): number {
     return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
 }
 
@@ -406,7 +409,7 @@ export function dot(a: Vec4, b: Vec4): number {
  * @param w the third vector
  * @returns result
  */
-export function cross(out: Vec4, u: Vec4, v: Vec4, w: Vec4): Vec4 {
+export function cross(out: Vec4, u: RVec4, v: RVec4, w: RVec4): Vec4 {
     const A = v[0] * w[1] - v[1] * w[0];
     const B = v[0] * w[2] - v[2] * w[0];
     const C = v[0] * w[3] - v[3] * w[0];
@@ -435,7 +438,7 @@ export function cross(out: Vec4, u: Vec4, v: Vec4, w: Vec4): Vec4 {
  * @param t interpolation amount, in the range [0-1], between the two inputs
  * @returns out
  */
-export function lerp(out: Vec4, a: Vec4, b: Vec4, t: number): Vec4 {
+export function lerp(out: Vec4, a: RVec4, b: RVec4, t: number): Vec4 {
     const ax = a[0];
     const ay = a[1];
     const az = a[2];
@@ -458,7 +461,7 @@ export function lerp(out: Vec4, a: Vec4, b: Vec4, t: number): Vec4 {
  * @param t interpolation amount
  * @returns out
  */
-export function lagrange(out: Vec4, a: Vec4, b: Vec4, c: Vec4, t: number): Vec4 {
+export function lagrange(out: Vec4, a: RVec4, b: RVec4, c: RVec4, t: number): Vec4 {
     const c0 = 2 * (t - 1) * (t - 0.5);
     const c1 = -4 * (t - 1) * t;
     const c2 = 2 * (t - 0.5) * t;
@@ -477,7 +480,7 @@ export function lagrange(out: Vec4, a: Vec4, b: Vec4, c: Vec4, t: number): Vec4
  * @param m matrix to transform with
  * @returns out
  */
-export function transformMat4(out: Vec4, a: Vec4, m: Mat4): Vec4 {
+export function transformMat4(out: Vec4, a: RVec4, m: RMat4): Vec4 {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -497,7 +500,7 @@ export function transformMat4(out: Vec4, a: Vec4, m: Mat4): Vec4 {
  * @param q quaternion to transform with
  * @returns out
  */
-export function transformQuat(out: Vec4, a: Vec4, q: Quat): Vec4 {
+export function transformQuat(out: Vec4, a: RVec4, q: RQuat): Vec4 {
     const x = a[0];
     const y = a[1];
     const z = a[2];
@@ -540,7 +543,7 @@ export function zero(out: Vec4): Vec4 {
  * @param a vector to represent as a string
  * @returns string representation of the vector
  */
-export function str(a: Vec4): string {
+export function str(a: RVec4): string {
     return `vec4(${a[0]}, ${a[1]}, ${a[2]}, ${a[3]})`;
 }
 
@@ -551,7 +554,7 @@ export function str(a: Vec4): string {
  * @param b The second vector.
  * @returns True if the vectors are equal, false otherwise.
  */
-export function exactEquals(a: Vec4, b: Vec4): boolean {
+export function exactEquals(a: RVec4, b: RVec4): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2] && a[3] === b[3];
 }
 
@@ -562,7 +565,7 @@ export function exactEquals(a: Vec4, b: Vec4): boolean {
  * @param b The second vector.
  * @returns True if the vectors are equal, false otherwise.
  */
-export function equals(a: Vec4, b: Vec4): boolean {
+export function equals(a: RVec4, b: RVec4): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -584,7 +587,7 @@ export function equals(a: Vec4, b: Vec4): boolean {
  * @param a vector to test
  * @returns whether or not the vector is finite
  */
-export function finite(a: Vec4): boolean {
+export function finite(a: RVec4): boolean {
     return Number.isFinite(a[0]) && Number.isFinite(a[1]) && Number.isFinite(a[2]) && Number.isFinite(a[3]);
 }
 
diff --git a/src/geometry/circumcircle.ts b/src/geometry/circumcircle.ts
index d3358eb..4b1ffda 100644
--- a/src/geometry/circumcircle.ts
+++ b/src/geometry/circumcircle.ts
@@ -1,5 +1,5 @@
 import { EPSILON } from '../core/scalar';
-import type { Vec2 } from '../core/vec2';
+import type { RVec2 } from '../core/vec2';
 import type { Circle } from '../shapes/circle';
 
 /**
@@ -8,7 +8,7 @@ import type { Circle } from '../shapes/circle';
  * @param triangle The triangle defined by three points
  * @returns
  */
-export function circumcircle(out: Circle, a: Vec2, b: Vec2, c: Vec2): Circle {
+export function circumcircle(out: Circle, a: RVec2, b: RVec2, c: RVec2): Circle {
     // work relative to `a` at the origin, which collapses the circumcenter formula
     const ax = a[0];
     const ay = a[1];
diff --git a/src/geometry/polygon2-decompose.ts b/src/geometry/polygon2-decompose.ts
index 81b3af4..48d1ff8 100644
--- a/src/geometry/polygon2-decompose.ts
+++ b/src/geometry/polygon2-decompose.ts
@@ -27,7 +27,7 @@ function triArea(ax: number, ay: number, bx: number, by: number, cx: number, cy:
 }
 
 /** True if vertex `i` is reflex; delegates to the shared polygon2 primitive. */
-function isReflex(poly: number[], i: number): boolean {
+function isReflex(poly: readonly number[], i: number): boolean {
     return isReflexVertex(poly, poly.length >> 1, i);
 }
 
@@ -108,14 +108,14 @@ function segmentsIntersect(
 }
 
 /** Appends vertices `from..to-1` of `src` onto `dst` (both flat arrays). */
-function appendRange(dst: number[], src: number[], from: number, to: number): void {
+function appendRange(dst: number[], src: readonly number[], from: number, to: number): void {
     for (let k = from; k < to; k++) {
         dst.push(src[k * 2], src[k * 2 + 1]);
     }
 }
 
 /** Copies vertices `i..j` (cyclic) of `poly` into a new flat polygon. */
-function polygonCopy(poly: number[], i: number, j: number): number[] {
+function polygonCopy(poly: readonly number[], i: number, j: number): number[] {
     const s = poly.length >> 1;
     const out: number[] = [];
     if (i < j) {
@@ -128,7 +128,7 @@ function polygonCopy(poly: number[], i: number, j: number): number[] {
 }
 
 /** Writes a CCW copy of the first `n` vertices of `vertices` into `out`, reversing if needed. */
-function toCCW(out: number[], vertices: number[], n: number): number[] {
+function toCCW(out: number[], vertices: readonly number[], n: number): number[] {
     for (let k = 0; k < n * 2; k++) out[k] = vertices[k];
     if (signedArea(out, n) < 0) reverse(out, out, n);
     return out;
@@ -137,7 +137,7 @@ function toCCW(out: number[], vertices: number[], n: number): number[] {
 const QUICK_DECOMP_MAX_LEVEL = 100;
 
 /** True if vertices `a` and `b` can see each other without any edge blocking (segment test). */
-function canSeeSegment(poly: number[], a: number, b: number): boolean {
+function canSeeSegment(poly: readonly number[], a: number, b: number): boolean {
     const s = poly.length >> 1;
     const ax = poly[a * 2];
     const ay = poly[a * 2 + 1];
@@ -167,7 +167,7 @@ function canSeeSegment(poly: number[], a: number, b: number): boolean {
  * @param n number of vertices to read from `vertices`
  * @returns an array of convex sub-polygons, each a flat `[x0, y0, ...]` array (CCW)
  */
-export function decomposePolygon2Quick(vertices: number[], n: number): number[][] {
+export function decomposePolygon2Quick(vertices: readonly number[], n: number): number[][] {
     if (n < 3) return [];
 
     const out: number[][] = [];
@@ -338,7 +338,7 @@ export function decomposePolygon2Quick(vertices: number[], n: number): number[][
 }
 
 /** True if vertices `a` and `b` can see each other (visibility test, used by the quality decomposition). */
-function canSeeVisibility(poly: number[], a: number, b: number): boolean {
+function canSeeVisibility(poly: readonly number[], a: number, b: number): boolean {
     const s = poly.length >> 1;
     const ax = poly[a * 2];
     const ay = poly[a * 2 + 1];
@@ -377,7 +377,7 @@ function canSeeVisibility(poly: number[], a: number, b: number): boolean {
 }
 
 /** Finds the minimal set of cut edges (as [ax, ay, bx, by]) that convex-partition the polygon. */
-function getCutEdges(poly: number[]): number[][] {
+function getCutEdges(poly: readonly number[]): number[][] {
     const s = poly.length >> 1;
     let min: number[][] = [];
     let nDiags = Number.MAX_VALUE;
@@ -403,7 +403,7 @@ function getCutEdges(poly: number[]): number[][] {
 }
 
 /** Index of the vertex in `poly` with the exact coordinates (x, y), or -1. */
-function indexOfVertex(poly: number[], x: number, y: number): number {
+function indexOfVertex(poly: readonly number[], x: number, y: number): number {
     const s = poly.length >> 1;
     for (let k = 0; k < s; k++) {
         if (poly[k * 2] === x && poly[k * 2 + 1] === y) return k;
@@ -443,7 +443,7 @@ function sliceByEdges(poly: number[], cutEdges: number[][]): number[][] {
  * @param n number of vertices to read from `vertices`
  * @returns an array of convex sub-polygons, each a flat `[x0, y0, ...]` array (CCW)
  */
-export function decomposePolygon2Quality(vertices: number[], n: number): number[][] {
+export function decomposePolygon2Quality(vertices: readonly number[], n: number): number[][] {
     if (n < 3) return [];
     const poly = toCCW([], vertices, n);
     const edges = getCutEdges(poly);
diff --git a/src/geometry/polygon2-triangulate.ts b/src/geometry/polygon2-triangulate.ts
index c2ff230..4cc75d6 100644
--- a/src/geometry/polygon2-triangulate.ts
+++ b/src/geometry/polygon2-triangulate.ts
@@ -64,16 +64,24 @@ function intersectSeg(ax: number, ay: number, bx: number, by: number, cx: number
 }
 
 /** X of the original vertex referenced by working slot `s`. */
-const vX = (vertices: number[], indices: number[], s: number): number => vertices[(indices[s] & IDX_MASK) * 2];
+const vX = (vertices: readonly number[], indices: readonly number[], s: number): number => vertices[(indices[s] & IDX_MASK) * 2];
 /** Y of the original vertex referenced by working slot `s`. */
-const vY = (vertices: number[], indices: number[], s: number): number => vertices[(indices[s] & IDX_MASK) * 2 + 1];
+const vY = (vertices: readonly number[], indices: readonly number[], s: number): number =>
+    vertices[(indices[s] & IDX_MASK) * 2 + 1];
 
 /**
  * True if the segment between working slots `i` and `j` does not intersect any
  * polygon edge. `loose` uses proper-intersection only (tolerating collinear
  * contact), used as a fallback when no strict ear can be found.
  */
-function diagonalie(i: number, j: number, n: number, vertices: number[], indices: number[], loose: boolean): boolean {
+function diagonalie(
+    i: number,
+    j: number,
+    n: number,
+    vertices: readonly number[],
+    indices: readonly number[],
+    loose: boolean,
+): boolean {
     const d0x = vX(vertices, indices, i);
     const d0y = vY(vertices, indices, i);
     const d1x = vX(vertices, indices, j);
@@ -107,7 +115,14 @@ function diagonalie(i: number, j: number, n: number, vertices: number[], indices
 }
 
 /** True if the diagonal from slot `i` to slot `j` stays inside the cone at vertex `i`. */
-function inCone(i: number, j: number, n: number, vertices: number[], indices: number[], loose: boolean): boolean {
+function inCone(
+    i: number,
+    j: number,
+    n: number,
+    vertices: readonly number[],
+    indices: readonly number[],
+    loose: boolean,
+): boolean {
     const ax = vX(vertices, indices, i);
     const ay = vY(vertices, indices, i);
     const bx = vX(vertices, indices, j);
@@ -130,7 +145,14 @@ function inCone(i: number, j: number, n: number, vertices: number[], indices: nu
     return !(area2(ax, ay, bx, by, nx, ny) <= 0 && area2(bx, by, ax, ay, px, py) <= 0);
 }
 
-function diagonal(i: number, j: number, n: number, vertices: number[], indices: number[], loose: boolean): boolean {
+function diagonal(
+    i: number,
+    j: number,
+    n: number,
+    vertices: readonly number[],
+    indices: readonly number[],
+    loose: boolean,
+): boolean {
     return inCone(i, j, n, vertices, indices, loose) && diagonalie(i, j, n, vertices, indices, loose);
 }
 
@@ -148,7 +170,7 @@ function diagonal(i: number, j: number, n: number, vertices: number[], indices:
  * @param n number of vertices to read from `vertices`
  * @returns the number of triangles written
  */
-export function triangulatePolygon2(out: number[], vertices: number[], n: number): number {
+export function triangulatePolygon2(out: number[], vertices: readonly number[], n: number): number {
     if (n < 3) return 0;
 
     // Order the working list so the resolved polygon winds clockwise, which is
diff --git a/src/geometry/quickhull2.ts b/src/geometry/quickhull2.ts
index 42dd699..7ba4983 100644
--- a/src/geometry/quickhull2.ts
+++ b/src/geometry/quickhull2.ts
@@ -9,7 +9,7 @@ const EPSILON = 1e-10;
  * @param points flat array of 2D points: [x0, y0, x1, y1, ...]
  * @returns indices of hull vertices in ccw order
  */
-export function quickhull2(points: number[]): number[] {
+export function quickhull2(points: readonly number[]): number[] {
     const n = Math.floor(points.length / 2);
     if (n < 3) return Array.from({ length: n }, (_, i) => i);
 
@@ -69,7 +69,7 @@ export function quickhull2(points: number[]): number[] {
  * Finds points on convex hull from set Sk that are on the right side of oriented line from P to Q.
  * Points are inserted into hull array at the end (before the final endpoint).
  */
-function findHull(points: number[], sk: number[], p: number, q: number, hull: number[]): void {
+function findHull(points: readonly number[], sk: readonly number[], p: number, q: number, hull: number[]): void {
     if (sk.length === 0) return;
 
     // find farthest point C from segment PQ
@@ -123,7 +123,7 @@ function findHull(points: number[], sk: number[], p: number, q: number, hull: nu
  * < 0: p3 is on the right of line p1→p2 (clockwise)
  * = 0: collinear
  */
-function crossProduct(points: number[], p1: number, p2: number, p3: number): number {
+function crossProduct(points: readonly number[], p1: number, p2: number, p3: number): number {
     const x1 = points[p1 * 2];
     const y1 = points[p1 * 2 + 1];
     const x2 = points[p2 * 2];
diff --git a/src/geometry/quickhull3.ts b/src/geometry/quickhull3.ts
index d70fd9f..d4440b5 100644
--- a/src/geometry/quickhull3.ts
+++ b/src/geometry/quickhull3.ts
@@ -54,7 +54,7 @@ type VertexList = {
 };
 
 type HullState = {
-    points: number[];
+    points: readonly number[];
     tolerance: number;
     faces: Face[];
     newFaces: Face[];
@@ -69,7 +69,7 @@ type HullState = {
  * @param points An array of numbers representing the 3D points (x1, y1, z1, x2, y2, z2, ...)
  * @returns An array of indices representing the triangles of the convex hull (i1, j1, k1, i2, j2, k2, ...).
  */
-export function quickhull3(points: number[]): number[] {
+export function quickhull3(points: readonly number[]): number[] {
     const n = points.length / 3;
     if (n < 4) return [];
 
@@ -91,7 +91,7 @@ export function quickhull3(points: number[]): number[] {
 
 // Hull state management
 
-function createHullState(points: number[], n: number): HullState {
+function createHullState(points: readonly number[], n: number): HullState {
     const vertices: VertexNode[] = [];
     for (let i = 0; i < n; i++) {
         vertices.push(createVertexNode(i));
@@ -294,7 +294,7 @@ function faceGetEdge(face: Face, i: number): HalfEdge | null {
     return edge;
 }
 
-function faceCompute(face: Face, points: number[]): void {
+function faceCompute(face: Face, points: readonly number[]): void {
     const a = halfEdgeTail(face.edge!)!;
     const b = halfEdgeHead(face.edge!);
     const c = halfEdgeHead(face.edge!.next!);
@@ -346,7 +346,7 @@ function faceCompute(face: Face, points: number[]): void {
     face.constant = face.normal[0] * face.midpoint[0] + face.normal[1] * face.midpoint[1] + face.normal[2] * face.midpoint[2];
 }
 
-function faceDistanceToPoint(face: Face, points: number[], vertexIndex: number): number {
+function faceDistanceToPoint(face: Face, points: readonly number[], vertexIndex: number): number {
     const idx = vertexIndex * 3;
     return face.normal[0] * points[idx] + face.normal[1] * points[idx + 1] + face.normal[2] * points[idx + 2] - face.constant;
 }
@@ -720,7 +720,13 @@ function reindexFaces(state: HullState): void {
 
 // Helper functions
 
-function computePlane(points: number[], v0: number, v1: number, v2: number, outNormal: [number, number, number]): number {
+function computePlane(
+    points: readonly number[],
+    v0: number,
+    v1: number,
+    v2: number,
+    outNormal: [number, number, number],
+): number {
     const p0x = points[v0 * 3];
     const p0y = points[v0 * 3 + 1];
     const p0z = points[v0 * 3 + 2];
@@ -758,14 +764,19 @@ function computePlane(points: number[], v0: number, v1: number, v2: number, outN
     return -(outNormal[0] * p0x + outNormal[1] * p0y + outNormal[2] * p0z);
 }
 
-function distanceToPlane(points: number[], idx: number, normal: [number, number, number], offset: number): number {
+function distanceToPlane(
+    points: readonly number[],
+    idx: number,
+    normal: readonly [number, number, number],
+    offset: number,
+): number {
     const x = points[idx * 3];
     const y = points[idx * 3 + 1];
     const z = points[idx * 3 + 2];
     return normal[0] * x + normal[1] * y + normal[2] * z + offset;
 }
 
-function distanceToLineSquared(points: number[], idx: number, v0: number, v1: number): number {
+function distanceToLineSquared(points: readonly number[], idx: number, v0: number, v1: number): number {
     const px = points[idx * 3];
     const py = points[idx * 3 + 1];
     const pz = points[idx * 3 + 2];
diff --git a/src/ik/fabrik2.ts b/src/ik/fabrik2.ts
index 14f7f15..8e73b6d 100644
--- a/src/ik/fabrik2.ts
+++ b/src/ik/fabrik2.ts
@@ -1,4 +1,4 @@
-import { type Vec2, vec2 } from '../core';
+import { type RVec2, type Vec2, vec2 } from '../core';
 
 // FABRIK (Forward And Backward Reaching Inverse Kinematics) for 2D chains.
 //
@@ -179,7 +179,7 @@ const DEGENERATE_SQUARED_LENGTH = 1e-24;
  * A zero axis cannot be normalized, and storing one turns every constraint that reads it into NaN.
  * `out` always starts as a valid unit vector, so keeping it is the safe fallback.
  */
-function setUnitAxis(out: Vec2, axis: Vec2): Vec2 {
+function setUnitAxis(out: Vec2, axis: RVec2): Vec2 {
     if (!hasDirection(vec2.squaredLength(axis))) return out;
     return vec2.normalize(out, axis);
 }
@@ -239,7 +239,7 @@ export function createChain2(): Chain2 {
  * @param joint the bone's joint, or a fresh unconstrained one if omitted
  * @returns the appended bone
  */
-export function addBone(chain: Chain2, start: Vec2, end: Vec2, joint: Joint2 = createJoint2()): Bone2 {
+export function addBone(chain: Chain2, start: RVec2, end: RVec2, joint: Joint2 = createJoint2()): Bone2 {
     const bone: Bone2 = {
         start: [start[0], start[1]],
         end: [end[0], end[1]],
@@ -272,7 +272,7 @@ const _addConsecutive_end: Vec2 = [0, 0];
  * @param joint the bone's joint, or a fresh unconstrained one if omitted
  * @returns the appended bone
  */
-export function addConsecutiveBone(chain: Chain2, direction: Vec2, length: number, joint: Joint2 = createJoint2()): Bone2 {
+export function addConsecutiveBone(chain: Chain2, direction: RVec2, length: number, joint: Joint2 = createJoint2()): Bone2 {
     const previous = chain.bones[chain.bones.length - 1];
 
     _addConsecutive_end[0] = previous.end[0] + direction[0] * length;
@@ -295,7 +295,7 @@ export function addConsecutiveBone(chain: Chain2, direction: Vec2, length: numbe
  * @param joint the joint for the junction this creates, or a fresh unconstrained one if omitted
  * @returns the prepended bone
  */
-export function addBoneAtBase(chain: Chain2, direction: Vec2, length: number, joint: Joint2 = createJoint2()): Bone2 {
+export function addBoneAtBase(chain: Chain2, direction: RVec2, length: number, joint: Joint2 = createJoint2()): Bone2 {
     const first = chain.bones[0];
 
     // the bone that was first now sits at index 1, so it is the one whose joint governs the new
@@ -345,7 +345,7 @@ export function setLocalJoint(joint: Joint2, clockwise: number, anticlockwise: n
  * @param anticlockwise how far it may swing anticlockwise, in radians, clamped to [0, PI]
  * @returns the joint
  */
-export function setGlobalJoint(joint: Joint2, axis: Vec2, clockwise: number, anticlockwise: number): Joint2 {
+export function setGlobalJoint(joint: Joint2, axis: RVec2, clockwise: number, anticlockwise: number): Joint2 {
     joint.coordinateSystem = ConstraintCoordinateSystem.GLOBAL;
     joint.clockwise = clampAngle(clockwise);
     joint.anticlockwise = clampAngle(anticlockwise);
@@ -366,7 +366,7 @@ export function setGlobalJoint(joint: Joint2, axis: Vec2, clockwise: number, ant
 export function setBaseboneConstraint(
     chain: Chain2,
     type: BaseboneConstraintType,
-    axis: Vec2,
+    axis: RVec2,
     clockwise: number,
     anticlockwise: number,
 ): Chain2 {
@@ -385,7 +385,7 @@ export function setBaseboneConstraint(
  *
  * The next {@link backward} or {@link solve} pulls the chain to it.
  */
-export function setBaseLocation(chain: Chain2, base: Vec2): Chain2 {
+export function setBaseLocation(chain: Chain2, base: RVec2): Chain2 {
     chain.base[0] = base[0];
     chain.base[1] = base[1];
     return chain;
@@ -397,7 +397,7 @@ export function setBaseLocation(chain: Chain2, base: Vec2): Chain2 {
  * A dead-straight chain is the worst starting pose for {@link solve} - see the note there. Bend
  * `direction` slightly between bones instead if the chain will be solved cold.
  */
-export function straighten(chain: Chain2, direction: Vec2): Chain2 {
+export function straighten(chain: Chain2, direction: RVec2): Chain2 {
     const bones = chain.bones;
 
     let x = chain.base[0];
@@ -450,7 +450,7 @@ export function getBoneAngle(chain: Chain2, index: number): number {
 }
 
 /** Whether `target` is within reach of the chain's base, so a solve can place the effector exactly on it. */
-export function isReachable(chain: Chain2, target: Vec2): boolean {
+export function isReachable(chain: Chain2, target: RVec2): boolean {
     return vec2.squaredDistance(chain.base, target) <= chain.length * chain.length;
 }
 
@@ -467,7 +467,7 @@ export function isReachable(chain: Chain2, target: Vec2): boolean {
  * @param target where the end effector should go
  * @returns the chain
  */
-export function forward(chain: Chain2, target: Vec2): Chain2 {
+export function forward(chain: Chain2, target: RVec2): Chain2 {
     const bones = chain.bones;
     const count = bones.length;
 
@@ -540,7 +540,7 @@ export function forward(chain: Chain2, target: Vec2): Chain2 {
  * @param base where the base should go, used only when `chain.fixedBase` is set
  * @returns the chain
  */
-export function backward(chain: Chain2, base: Vec2): Chain2 {
+export function backward(chain: Chain2, base: RVec2): Chain2 {
     const bones = chain.bones;
     const count = bones.length;
 
@@ -612,7 +612,7 @@ export function backward(chain: Chain2, base: Vec2): Chain2 {
  * @param target where the end effector should go
  * @returns the distance from the effector to `target` afterwards
  */
-export function iterate(chain: Chain2, target: Vec2): number {
+export function iterate(chain: Chain2, target: RVec2): number {
     if (chain.bones.length === 0) return Number.POSITIVE_INFINITY;
 
     forward(chain, target);
@@ -634,7 +634,7 @@ export function iterate(chain: Chain2, target: Vec2): number {
  * @param target where the end effector should go
  * @returns the distance from the effector to `target`, also stored as `chain.solveDistance`
  */
-export function solve(chain: Chain2, target: Vec2): number {
+export function solve(chain: Chain2, target: RVec2): number {
     const count = chain.bones.length;
 
     if (count === 0) {
@@ -727,7 +727,7 @@ export function connectChain(
  * @param structure the structure to solve, mutated in place
  * @param target the target for every chain that does not use an embedded target
  */
-export function solveStructure(structure: Structure2, target: Vec2): void {
+export function solveStructure(structure: Structure2, target: RVec2): void {
     const chains = structure.chains;
 
     for (let i = 0; i < chains.length; i++) {
@@ -783,7 +783,7 @@ function clampAngle(radians: number): number {
  * Clamps the direction `(x, y)` into the wedge reaching `clockwise` one way and `anticlockwise` the
  * other from `baseline`. Writes `_pass_direction`.
  */
-function constrainToWedge(x: number, y: number, baseline: Vec2, clockwise: number, anticlockwise: number): void {
+function constrainToWedge(x: number, y: number, baseline: RVec2, clockwise: number, anticlockwise: number): void {
     _pass_direction[0] = x;
     _pass_direction[1] = y;
 
diff --git a/src/ik/fabrik3.ts b/src/ik/fabrik3.ts
index 2702d66..708b2ab 100644
--- a/src/ik/fabrik3.ts
+++ b/src/ik/fabrik3.ts
@@ -1,4 +1,4 @@
-import { type Mat3, mat3, type Quat, quat, type Vec3, vec3 } from '../core';
+import { type Mat3, mat3, type Quat, quat, type RVec3, type Vec3, vec3 } from '../core';
 
 // FABRIK (Forward And Backward Reaching Inverse Kinematics) for 3D chains.
 //
@@ -210,7 +210,7 @@ const DEGENERATE_SQUARED_LENGTH = 1e-24;
  * A zero axis cannot be normalized, and storing one turns every constraint that reads it into NaN.
  * `out` always starts as a valid unit vector, so keeping it is the safe fallback.
  */
-function setUnitAxis(out: Vec3, axis: Vec3): Vec3 {
+function setUnitAxis(out: Vec3, axis: RVec3): Vec3 {
     if (!hasDirection(vec3.squaredLength(axis))) return out;
     return vec3.normalize(out, axis);
 }
@@ -275,7 +275,7 @@ export function createChain3(): Chain3 {
  * @param joint the bone's joint, or a fresh unconstrained one if omitted
  * @returns the appended bone
  */
-export function addBone(chain: Chain3, start: Vec3, end: Vec3, joint: Joint3 = createJoint3()): Bone3 {
+export function addBone(chain: Chain3, start: RVec3, end: RVec3, joint: Joint3 = createJoint3()): Bone3 {
     const bone: Bone3 = {
         start: [start[0], start[1], start[2]],
         end: [end[0], end[1], end[2]],
@@ -307,7 +307,7 @@ export function addBone(chain: Chain3, start: Vec3, end: Vec3, joint: Joint3 = c
  * @param joint the bone's joint, or a fresh unconstrained one if omitted
  * @returns the appended bone
  */
-export function addConsecutiveBone(chain: Chain3, direction: Vec3, length: number, joint: Joint3 = createJoint3()): Bone3 {
+export function addConsecutiveBone(chain: Chain3, direction: RVec3, length: number, joint: Joint3 = createJoint3()): Bone3 {
     const previous = chain.bones[chain.bones.length - 1];
 
     _addConsecutive_end[0] = previous.end[0] + direction[0] * length;
@@ -333,7 +333,7 @@ const _addConsecutive_end: Vec3 = [0, 0, 0];
  * @param joint the joint for the junction this creates, or a fresh unconstrained one if omitted
  * @returns the prepended bone
  */
-export function addBoneAtBase(chain: Chain3, direction: Vec3, length: number, joint: Joint3 = createJoint3()): Bone3 {
+export function addBoneAtBase(chain: Chain3, direction: RVec3, length: number, joint: Joint3 = createJoint3()): Bone3 {
     const first = chain.bones[0];
 
     // the bone that was first now sits at index 1, so it is the one whose joint governs the new
@@ -395,10 +395,10 @@ export function setBallJoint(joint: Joint3, rotor: number): Joint3 {
 export function setHingeJoint(
     joint: Joint3,
     type: JointType.GLOBAL_HINGE | JointType.LOCAL_HINGE,
-    rotationAxis: Vec3,
+    rotationAxis: RVec3,
     clockwise: number,
     anticlockwise: number,
-    referenceAxis: Vec3,
+    referenceAxis: RVec3,
 ): Joint3 {
     joint.type = type;
     joint.clockwise = clampAngle(clockwise);
@@ -422,7 +422,7 @@ export function setHingeJoint(
 export function setBaseboneRotorConstraint(
     chain: Chain3,
     type: BaseboneConstraintType.GLOBAL_ROTOR | BaseboneConstraintType.LOCAL_ROTOR,
-    axis: Vec3,
+    axis: RVec3,
     rotor: number,
 ): Chain3 {
     chain.baseboneConstraintType = type;
@@ -449,10 +449,10 @@ export function setBaseboneRotorConstraint(
 export function setBaseboneHingeConstraint(
     chain: Chain3,
     type: BaseboneConstraintType.GLOBAL_HINGE | BaseboneConstraintType.LOCAL_HINGE,
-    rotationAxis: Vec3,
+    rotationAxis: RVec3,
     clockwise: number,
     anticlockwise: number,
-    referenceAxis: Vec3,
+    referenceAxis: RVec3,
 ): Chain3 {
     chain.baseboneConstraintType = type;
     chain.baseboneClockwise = clampAngle(clockwise);
@@ -472,7 +472,7 @@ export function setBaseboneHingeConstraint(
  *
  * The next {@link backward} or {@link solve} pulls the chain to it.
  */
-export function setBaseLocation(chain: Chain3, base: Vec3): Chain3 {
+export function setBaseLocation(chain: Chain3, base: RVec3): Chain3 {
     chain.base[0] = base[0];
     chain.base[1] = base[1];
     chain.base[2] = base[2];
@@ -485,7 +485,7 @@ export function setBaseLocation(chain: Chain3, base: Vec3): Chain3 {
  * A dead-straight chain is the worst starting pose for {@link solve} - see the note there. Bend
  * `direction` slightly between bones instead if the chain will be solved cold.
  */
-export function straighten(chain: Chain3, direction: Vec3): Chain3 {
+export function straighten(chain: Chain3, direction: RVec3): Chain3 {
     const bones = chain.bones;
 
     let x = chain.base[0];
@@ -535,7 +535,7 @@ export function getBoneDirection(out: Vec3, chain: Chain3, index: number): Vec3
  * Use it to orient a mesh along a bone, passing whichever axis the mesh is modelled along - `up`
  * is `[0, 1, 0]` for a cylinder or capsule built along Y. The roll about the bone is arbitrary.
  */
-export function getBoneRotation(out: Quat, chain: Chain3, index: number, up: Vec3): Quat {
+export function getBoneRotation(out: Quat, chain: Chain3, index: number, up: RVec3): Quat {
     getBoneDirection(_boneRotation_direction, chain, index);
     return quat.rotationTo(out, up, _boneRotation_direction);
 }
@@ -543,7 +543,7 @@ export function getBoneRotation(out: Quat, chain: Chain3, index: number, up: Vec
 const _boneRotation_direction: Vec3 = [0, 0, 0];
 
 /** Whether `target` is within reach of the chain's base, so a solve can place the effector exactly on it. */
-export function isReachable(chain: Chain3, target: Vec3): boolean {
+export function isReachable(chain: Chain3, target: RVec3): boolean {
     return vec3.squaredDistance(chain.base, target) <= chain.length * chain.length;
 }
 
@@ -560,7 +560,7 @@ export function isReachable(chain: Chain3, target: Vec3): boolean {
  * @param target where the end effector should go
  * @returns the chain
  */
-export function forward(chain: Chain3, target: Vec3): Chain3 {
+export function forward(chain: Chain3, target: RVec3): Chain3 {
     const bones = chain.bones;
     const count = bones.length;
 
@@ -644,7 +644,7 @@ export function forward(chain: Chain3, target: Vec3): Chain3 {
  * @param base where the base should go, used only when `chain.fixedBase` is set
  * @returns the chain
  */
-export function backward(chain: Chain3, base: Vec3): Chain3 {
+export function backward(chain: Chain3, base: RVec3): Chain3 {
     const bones = chain.bones;
     const count = bones.length;
 
@@ -727,7 +727,7 @@ export function backward(chain: Chain3, base: Vec3): Chain3 {
  * @param target where the end effector should go
  * @returns the distance from the effector to `target` afterwards
  */
-export function iterate(chain: Chain3, target: Vec3): number {
+export function iterate(chain: Chain3, target: RVec3): number {
     if (chain.bones.length === 0) return Number.POSITIVE_INFINITY;
 
     forward(chain, target);
@@ -749,7 +749,7 @@ export function iterate(chain: Chain3, target: Vec3): number {
  * @param target where the end effector should go
  * @returns the distance from the effector to `target`, also stored as `chain.solveDistance`
  */
-export function solve(chain: Chain3, target: Vec3): number {
+export function solve(chain: Chain3, target: RVec3): number {
     const count = chain.bones.length;
 
     if (count === 0) {
@@ -843,7 +843,7 @@ export function connectChain(
  * @param structure the structure to solve, mutated in place
  * @param target the target for every chain that does not use an embedded target
  */
-export function solveStructure(structure: Structure3, target: Vec3): void {
+export function solveStructure(structure: Structure3, target: RVec3): void {
     const chains = structure.chains;
 
     for (let i = 0; i < chains.length; i++) {
@@ -902,7 +902,7 @@ function clampAngle(radians: number): number {
 }
 
 /** Writes the component of `a` perpendicular to the unit vector `axis`, normalized. */
-function orthonormalize(out: Vec3, a: Vec3, axis: Vec3): Vec3 {
+function orthonormalize(out: Vec3, a: RVec3, axis: RVec3): Vec3 {
     // the component of `a` perpendicular to `axis`
     vec3.scaleAndAdd(out, a, axis, -vec3.dot(a, axis));
 
@@ -915,7 +915,7 @@ function orthonormalize(out: Vec3, a: Vec3, axis: Vec3): Vec3 {
 }
 
 /** Normalizes `out` in place, falling back to `fallback` when it has no length. */
-function normalizeOr(out: Vec3, fallback: Vec3): Vec3 {
+function normalizeOr(out: Vec3, fallback: RVec3): Vec3 {
     if (!hasDirection(vec3.squaredLength(out))) {
         return vec3.copy(out, fallback);
     }
@@ -928,7 +928,7 @@ function normalizeOr(out: Vec3, fallback: Vec3): Vec3 {
  * A direction parallel to the hinge axis projects to nothing, leaving no in-plane direction to
  * normalize. Fall back to the hinge's reference axis, which lies in the plane by construction.
  */
-function projectOntoHinge(out: Vec3, x: number, y: number, z: number, axis: Vec3, referenceAxis: Vec3): Vec3 {
+function projectOntoHinge(out: Vec3, x: number, y: number, z: number, axis: RVec3, referenceAxis: RVec3): Vec3 {
     const d = x * axis[0] + y * axis[1] + z * axis[2];
 
     out[0] = x - axis[0] * d;
@@ -943,7 +943,7 @@ function projectOntoHinge(out: Vec3, x: number, y: number, z: number, axis: Vec3
 }
 
 /** Rotates the unit vector `a` about the unit vector `axis` by `radians` (Rodrigues). */
-function rotateAboutAxis(out: Vec3, a: Vec3, axis: Vec3, radians: number): Vec3 {
+function rotateAboutAxis(out: Vec3, a: RVec3, axis: RVec3, radians: number): Vec3 {
     const ax = a[0];
     const ay = a[1];
     const az = a[2];
@@ -967,8 +967,8 @@ function constrainHinge(
     x: number,
     y: number,
     z: number,
-    axis: Vec3,
-    referenceAxis: Vec3,
+    axis: RVec3,
+    referenceAxis: RVec3,
     clockwise: number,
     anticlockwise: number,
 ): void {
@@ -998,7 +998,7 @@ function constrainHinge(
  * except at `direction` = (0, 0, -1), where the basis flips, so a local hinge whose parent swings
  * through there will pop.
  */
-function basisFromDirection(out: Mat3, direction: Vec3): Mat3 {
+function basisFromDirection(out: Mat3, direction: RVec3): Mat3 {
     const x = direction[0];
     const y = direction[1];
     const z = direction[2];
diff --git a/src/random/random.ts b/src/random/random.ts
index a816423..be19e73 100644
--- a/src/random/random.ts
+++ b/src/random/random.ts
@@ -50,7 +50,7 @@ export function sign(random: RandomGenerator, plusChance = 0.5): number {
  * @param items the array to choose from
  * @throws if the array is empty
  */
-export function choice<T>(random: RandomGenerator, items: T[]): T {
+export function choice<T>(random: RandomGenerator, items: readonly T[]): T {
     if (items.length === 0) {
         throw new Error('cannot choose from an empty array');
     }
diff --git a/src/shapes/box2.ts b/src/shapes/box2.ts
index eee017b..66f9780 100644
--- a/src/shapes/box2.ts
+++ b/src/shapes/box2.ts
@@ -1,10 +1,13 @@
 import { EPSILON } from '../core/scalar';
-import type { Vec2 } from '../core/vec2';
-import type { Circle } from './circle';
+import type { RVec2, Vec2 } from '../core/vec2';
+import type { RCircle } from './circle';
 
 /** An axis-aligned box in 2D space, as [minX, minY, maxX, maxY] */
 export type Box2 = [minX: number, minY: number, maxX: number, maxY: number];
 
+/** A read-only 2D axis-aligned bounding box */
+export type RBox2 = Readonly<Box2>;
+
 /**
  * Create a new empty Box2 with "min" set to positive infinity and "max" set to negative infinity
  * @returns A new Box2
@@ -18,7 +21,7 @@ export function create(): Box2 {
  * @param box - A Box2 to clone
  * @returns a clone of box
  */
-export function clone(box: Box2): Box2 {
+export function clone(box: RBox2): Box2 {
     return [box[0], box[1], box[2], box[3]];
 }
 
@@ -28,7 +31,7 @@ export function clone(box: Box2): Box2 {
  * @param box the input Box2
  * @returns the output Box2
  */
-export function copy(out: Box2, box: Box2): Box2 {
+export function copy(out: Box2, box: RBox2): Box2 {
     out[0] = box[0];
     out[1] = box[1];
     out[2] = box[2];
@@ -60,7 +63,7 @@ export function set(out: Box2, minX: number, minY: number, maxX: number, maxY: n
  * @param max - The maximum corner
  * @returns The updated Box2
  */
-export function setFromVectors(out: Box2, min: Vec2, max: Vec2): Box2 {
+export function setFromVectors(out: Box2, min: RVec2, max: RVec2): Box2 {
     out[0] = min[0];
     out[1] = min[1];
     out[2] = max[0];
@@ -74,7 +77,7 @@ export function setFromVectors(out: Box2, min: Vec2, max: Vec2): Box2 {
  * @param box - The input Box2
  * @returns The minimum corner
  */
-export function min(out: Vec2, box: Box2): Vec2 {
+export function min(out: Vec2, box: RBox2): Vec2 {
     out[0] = box[0];
     out[1] = box[1];
     return out;
@@ -86,7 +89,7 @@ export function min(out: Vec2, box: Box2): Vec2 {
  * @param box - The input Box2
  * @returns The maximum corner
  */
-export function max(out: Vec2, box: Box2): Vec2 {
+export function max(out: Vec2, box: RBox2): Vec2 {
     out[0] = box[2];
     out[1] = box[3];
     return out;
@@ -111,7 +114,7 @@ export function empty(out: Box2): Box2 {
  * @param b - The second box
  * @returns True if the boxes are equal, false otherwise
  */
-export function exactEquals(a: Box2, b: Box2): boolean {
+export function exactEquals(a: RBox2, b: RBox2): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2] && a[3] === b[3];
 }
 
@@ -121,7 +124,7 @@ export function exactEquals(a: Box2, b: Box2): boolean {
  * @param b - The second box
  * @returns True if the boxes are equal, false otherwise
  */
-export function equals(a: Box2, b: Box2): boolean {
+export function equals(a: RBox2, b: RBox2): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -145,7 +148,7 @@ export function equals(a: Box2, b: Box2): boolean {
  * @param size - The size of the box
  * @returns The updated Box2
  */
-export function setFromCenterAndSize(out: Box2, center: Vec2, size: Vec2): Box2 {
+export function setFromCenterAndSize(out: Box2, center: RVec2, size: RVec2): Box2 {
     const hx = size[0] * 0.5;
     const hy = size[1] * 0.5;
     out[0] = center[0] - hx;
@@ -162,7 +165,7 @@ export function setFromCenterAndSize(out: Box2, center: Vec2, size: Vec2): Box2
  * @param point - The point to include
  * @returns The expanded Box2
  */
-export function expandByPoint(out: Box2, box: Box2, point: Vec2): Box2 {
+export function expandByPoint(out: Box2, box: RBox2, point: RVec2): Box2 {
     out[0] = Math.min(box[0], point[0]);
     out[1] = Math.min(box[1], point[1]);
     out[2] = Math.max(box[2], point[0]);
@@ -178,7 +181,7 @@ export function expandByPoint(out: Box2, box: Box2, point: Vec2): Box2 {
  * @param vector - The vector to expand by
  * @returns The expanded Box2
  */
-export function expandByExtents(out: Box2, box: Box2, vector: Vec2): Box2 {
+export function expandByExtents(out: Box2, box: RBox2, vector: RVec2): Box2 {
     out[0] = box[0] - vector[0];
     out[1] = box[1] - vector[1];
     out[2] = box[2] + vector[0];
@@ -194,7 +197,7 @@ export function expandByExtents(out: Box2, box: Box2, vector: Vec2): Box2 {
  * @param margin - The uniform margin to expand by
  * @returns The expanded Box2
  */
-export function expandByMargin(out: Box2, box: Box2, margin: number): Box2 {
+export function expandByMargin(out: Box2, box: RBox2, margin: number): Box2 {
     out[0] = box[0] - margin;
     out[1] = box[1] - margin;
     out[2] = box[2] + margin;
@@ -210,7 +213,7 @@ export function expandByMargin(out: Box2, box: Box2, margin: number): Box2 {
  * @param boxB - The second Box2
  * @returns The union Box2
  */
-export function union(out: Box2, boxA: Box2, boxB: Box2): Box2 {
+export function union(out: Box2, boxA: RBox2, boxB: RBox2): Box2 {
     out[0] = Math.min(boxA[0], boxB[0]);
     out[1] = Math.min(boxA[1], boxB[1]);
     out[2] = Math.max(boxA[2], boxB[2]);
@@ -224,7 +227,7 @@ export function union(out: Box2, boxA: Box2, boxB: Box2): Box2 {
  * @param box - The input Box2
  * @returns The center point
  */
-export function center(out: Vec2, box: Box2): Vec2 {
+export function center(out: Vec2, box: RBox2): Vec2 {
     out[0] = (box[0] + box[2]) * 0.5;
     out[1] = (box[1] + box[3]) * 0.5;
     return out;
@@ -236,7 +239,7 @@ export function center(out: Vec2, box: Box2): Vec2 {
  * @param box - The input Box2
  * @returns The extents (distance from center to each edge)
  */
-export function extents(out: Vec2, box: Box2): Vec2 {
+export function extents(out: Vec2, box: RBox2): Vec2 {
     out[0] = (box[2] - box[0]) * 0.5;
     out[1] = (box[3] - box[1]) * 0.5;
     return out;
@@ -248,7 +251,7 @@ export function extents(out: Vec2, box: Box2): Vec2 {
  * @param box - The input Box2
  * @returns The size (width, height)
  */
-export function size(out: Vec2, box: Box2): Vec2 {
+export function size(out: Vec2, box: RBox2): Vec2 {
     out[0] = box[2] - box[0];
     out[1] = box[3] - box[1];
     return out;
@@ -259,7 +262,7 @@ export function size(out: Vec2, box: Box2): Vec2 {
  * @param box - The input Box2
  * @returns The area (width * height)
  */
-export function area(box: Box2): number {
+export function area(box: RBox2): number {
     return (box[2] - box[0]) * (box[3] - box[1]);
 }
 
@@ -270,7 +273,7 @@ export function area(box: Box2): number {
  * @param scale - The scale to apply (as a Vec2)
  * @returns The scaled Box2
  */
-export function scale(out: Box2, box: Box2, scale: Vec2): Box2 {
+export function scale(out: Box2, box: RBox2, scale: RVec2): Box2 {
     const minX = box[0] * scale[0];
     const maxX = box[2] * scale[0];
     const minY = box[1] * scale[1];
@@ -291,7 +294,7 @@ export function scale(out: Box2, box: Box2, scale: Vec2): Box2 {
  * @param point - The point to test
  * @returns true if the point is inside or on the boundary of the box
  */
-export function containsPoint(box: Box2, point: Vec2): boolean {
+export function containsPoint(box: RBox2, point: RVec2): boolean {
     return point[0] >= box[0] && point[0] <= box[2] && point[1] >= box[1] && point[1] <= box[3];
 }
 
@@ -301,7 +304,7 @@ export function containsPoint(box: Box2, point: Vec2): boolean {
  * @param contained - The Box2 that might be contained
  * @returns true if the container Box2 completely contains the contained Box2
  */
-export function containsBox2(container: Box2, contained: Box2): boolean {
+export function containsBox2(container: RBox2, contained: RBox2): boolean {
     return (
         contained[0] >= container[0] &&
         contained[2] <= container[2] &&
@@ -313,14 +316,14 @@ export function containsBox2(container: Box2, contained: Box2): boolean {
 /**
  * Check whether two bounding boxes intersect
  */
-export function intersectsBox2(boxA: Box2, boxB: Box2): boolean {
+export function intersectsBox2(boxA: RBox2, boxB: RBox2): boolean {
     return boxA[0] <= boxB[2] && boxA[2] >= boxB[0] && boxA[1] <= boxB[3] && boxA[3] >= boxB[1];
 }
 
 /**
  * Test intersection between an axis-aligned bounding box and a circle.
  */
-export function intersectsCircle(box: Box2, circle: Circle): boolean {
+export function intersectsCircle(box: RBox2, circle: RCircle): boolean {
     const { center, radius } = circle;
     const cx = center[0];
     const cy = center[1];
diff --git a/src/shapes/box3.ts b/src/shapes/box3.ts
index ff3d349..be7f283 100644
--- a/src/shapes/box3.ts
+++ b/src/shapes/box3.ts
@@ -1,12 +1,15 @@
 import { EPSILON } from '../core/scalar';
-import type { Mat4 } from '../core/mat4';
-import type { Vec3 } from '../core/vec3';
-import type { Plane3 } from './plane3';
-import type { Sphere } from './sphere';
+import type { RMat4 } from '../core/mat4';
+import type { RVec3, Vec3 } from '../core/vec3';
+import type { RPlane3 } from './plane3';
+import type { RSphere } from './sphere';
 
 /** A box in 3D space */
 export type Box3 = [minX: number, minY: number, minZ: number, maxX: number, maxY: number, maxZ: number];
 
+/** A read-only 3D axis-aligned bounding box */
+export type RBox3 = Readonly<Box3>;
+
 /**
  * Create a new empty Box3 with "min" set to positive infinity and "max" set to negative infinity
  * @returns A new Box3
@@ -27,7 +30,7 @@ export function create(): Box3 {
  * @param box - A Box3 to clone
  * @returns a clone of box
  */
-export function clone(box: Box3): Box3 {
+export function clone(box: RBox3): Box3 {
     return [box[0], box[1], box[2], box[3], box[4], box[5]];
 }
 
@@ -37,7 +40,7 @@ export function clone(box: Box3): Box3 {
  * @param box the input Box3
  * @returns the output Box3
  */
-export function copy(out: Box3, box: Box3): Box3 {
+export function copy(out: Box3, box: RBox3): Box3 {
     out[0] = box[0];
     out[1] = box[1];
     out[2] = box[2];
@@ -75,7 +78,7 @@ export function set(out: Box3, minX: number, minY: number, minZ: number, maxX: n
  * @param max - The maximum corner
  * @returns The updated Box3
  */
-export function setFromVectors(out: Box3, min: Vec3, max: Vec3): Box3 {
+export function setFromVectors(out: Box3, min: RVec3, max: RVec3): Box3 {
     out[0] = min[0];
     out[1] = min[1];
     out[2] = min[2];
@@ -91,7 +94,7 @@ export function setFromVectors(out: Box3, min: Vec3, max: Vec3): Box3 {
  * @param box - The input Box3
  * @returns The minimum corner
  */
-export function min(out: Vec3, box: Box3): Vec3 {
+export function min(out: Vec3, box: RBox3): Vec3 {
     out[0] = box[0];
     out[1] = box[1];
     out[2] = box[2];
@@ -104,7 +107,7 @@ export function min(out: Vec3, box: Box3): Vec3 {
  * @param box - The input Box3
  * @returns The maximum corner
  */
-export function max(out: Vec3, box: Box3): Vec3 {
+export function max(out: Vec3, box: RBox3): Vec3 {
     out[0] = box[3];
     out[1] = box[4];
     out[2] = box[5];
@@ -132,7 +135,7 @@ export function empty(out: Box3): Box3 {
  * @param b - The second box
  * @returns True if the boxes are equal, false otherwise
  */
-export function exactEquals(a: Box3, b: Box3): boolean {
+export function exactEquals(a: RBox3, b: RBox3): boolean {
     return a[0] === b[0] && a[1] === b[1] && a[2] === b[2] && a[3] === b[3] && a[4] === b[4] && a[5] === b[5];
 }
 
@@ -142,7 +145,7 @@ export function exactEquals(a: Box3, b: Box3): boolean {
  * @param b - The second box
  * @returns True if the boxes are equal, false otherwise
  */
-export function equals(a: Box3, b: Box3): boolean {
+export function equals(a: RBox3, b: RBox3): boolean {
     const a0 = a[0];
     const a1 = a[1];
     const a2 = a[2];
@@ -172,7 +175,7 @@ export function equals(a: Box3, b: Box3): boolean {
  * @param size - The size of the box
  * @returns The updated Box3
  */
-export function setFromCenterAndSize(out: Box3, center: Vec3, size: Vec3): Box3 {
+export function setFromCenterAndSize(out: Box3, center: RVec3, size: RVec3): Box3 {
     const hx = size[0] * 0.5;
     const hy = size[1] * 0.5;
     const hz = size[2] * 0.5;
@@ -192,7 +195,7 @@ export function setFromCenterAndSize(out: Box3, center: Vec3, size: Vec3): Box3
  * @param point - The point to include
  * @returns The expanded Box3
  */
-export function expandByPoint(out: Box3, box: Box3, point: Vec3): Box3 {
+export function expandByPoint(out: Box3, box: RBox3, point: RVec3): Box3 {
     out[0] = Math.min(box[0], point[0]);
     out[1] = Math.min(box[1], point[1]);
     out[2] = Math.min(box[2], point[2]);
@@ -210,7 +213,7 @@ export function expandByPoint(out: Box3, box: Box3, point: Vec3): Box3 {
  * @param vector - The vector to expand by
  * @returns The expanded Box3
  */
-export function expandByExtents(out: Box3, box: Box3, vector: Vec3): Box3 {
+export function expandByExtents(out: Box3, box: RBox3, vector: RVec3): Box3 {
     out[0] = box[0] - vector[0];
     out[1] = box[1] - vector[1];
     out[2] = box[2] - vector[2];
@@ -228,7 +231,7 @@ export function expandByExtents(out: Box3, box: Box3, vector: Vec3): Box3 {
  * @param margin - The uniform margin to expand by
  * @returns The expanded Box3
  */
-export function expandByMargin(out: Box3, box: Box3, margin: number): Box3 {
+export function expandByMargin(out: Box3, box: RBox3, margin: number): Box3 {
     out[0] = box[0] - margin;
     out[1] = box[1] - margin;
     out[2] = box[2] - margin;
@@ -246,7 +249,7 @@ export function expandByMargin(out: Box3, box: Box3, margin: number): Box3 {
  * @param boxB - The second Box3
  * @returns The union Box3
  */
-export function union(out: Box3, boxA: Box3, boxB: Box3): Box3 {
+export function union(out: Box3, boxA: RBox3, boxB: RBox3): Box3 {
     out[0] = Math.min(boxA[0], boxB[0]);
     out[1] = Math.min(boxA[1], boxB[1]);
     out[2] = Math.min(boxA[2], boxB[2]);
@@ -262,7 +265,7 @@ export function union(out: Box3, boxA: Box3, boxB: Box3): Box3 {
  * @param box - The input Box3
  * @returns The center point
  */
-export function center(out: Vec3, box: Box3): Vec3 {
+export function center(out: Vec3, box: RBox3): Vec3 {
     out[0] = (box[0] + box[3]) * 0.5;
     out[1] = (box[1] + box[4]) * 0.5;
     out[2] = (box[2] + box[5]) * 0.5;
@@ -275,7 +278,7 @@ export function center(out: Vec3, box: Box3): Vec3 {
  * @param box - The input Box3
  * @returns The extents (distance from center to each face)
  */
-export function extents(out: Vec3, box: Box3): Vec3 {
+export function extents(out: Vec3, box: RBox3): Vec3 {
     out[0] = (box[3] - box[0]) * 0.5;
     out[1] = (box[4] - box[1]) * 0.5;
     out[2] = (box[5] - box[2]) * 0.5;
@@ -288,7 +291,7 @@ export function extents(out: Vec3, box: Box3): Vec3 {
  * @param box - The input Box3
  * @returns The size (width, height, depth)
  */
-export function size(out: Vec3, box: Box3): Vec3 {
+export function size(out: Vec3, box: RBox3): Vec3 {
     out[0] = box[3] - box[0];
     out[1] = box[4] - box[1];
     out[2] = box[5] - box[2];
@@ -300,7 +303,7 @@ export function size(out: Vec3, box: Box3): Vec3 {
  * @param box - The input Box3
  * @returns The surface area
  */
-export function surfaceArea(box: Box3): number {
+export function surfaceArea(box: RBox3): number {
     const width = box[3] - box[0];
     const height = box[4] - box[1];
     const depth = box[5] - box[2];
@@ -314,7 +317,7 @@ export function surfaceArea(box: Box3): number {
  * @param scale - The scale to apply (as a Vec3)
  * @returns The scaled Box3
  */
-export function scale(out: Box3, box: Box3, scale: Vec3): Box3 {
+export function scale(out: Box3, box: RBox3, scale: RVec3): Box3 {
     const minX = box[0] * scale[0];
     const maxX = box[3] * scale[0];
     const minY = box[1] * scale[1];
@@ -349,7 +352,7 @@ export function scale(out: Box3, box: Box3, scale: Vec3): Box3 {
  * @param mat - The 4x4 transformation matrix
  * @returns The transformed Box3
  */
-export function transformMat4(out: Box3, box: Box3, mat: Mat4): Box3 {
+export function transformMat4(out: Box3, box: RBox3, mat: RMat4): Box3 {
     const bMinX = box[0];
     const bMinY = box[1];
     const bMinZ = box[2];
@@ -410,7 +413,7 @@ export function transformMat4(out: Box3, box: Box3, mat: Mat4): Box3 {
  * @param point - The point to test
  * @returns true if the point is inside or on the boundary of the box
  */
-export function containsPoint(box: Box3, point: Vec3): boolean {
+export function containsPoint(box: RBox3, point: RVec3): boolean {
     return (
         point[0] >= box[0] &&
         point[0] <= box[3] &&
@@ -427,7 +430,7 @@ export function containsPoint(box: Box3, point: Vec3): boolean {
  * @param contained - The Box3 that might be contained
  * @returns true if the container Box3 completely contains the contained Box3
  */
-export function containsBox3(container: Box3, contained: Box3): boolean {
+export function containsBox3(container: RBox3, contained: RBox3): boolean {
     return (
         contained[0] >= container[0] &&
         contained[3] <= container[3] &&
@@ -441,7 +444,7 @@ export function containsBox3(container: Box3, contained: Box3): boolean {
 /**
  * Check whether two bounding boxes intersect
  */
-export function intersectsBox3(boxA: Box3, boxB: Box3): boolean {
+export function intersectsBox3(boxA: RBox3, boxB: RBox3): boolean {
     return (
         boxA[0] <= boxB[3] &&
         boxA[3] >= boxB[0] &&
@@ -465,7 +468,7 @@ export function intersectsBox3(boxA: Box3, boxB: Box3): boolean {
  * projections are needed. An all-zero cross axis (edge parallel to a box axis)
  * collapses every projection and the radius to 0, passing automatically.
  */
-export function intersectsTriangle3(box: Box3, a: Vec3, b: Vec3, c: Vec3): boolean {
+export function intersectsTriangle3(box: RBox3, a: RVec3, b: RVec3, c: RVec3): boolean {
     // Empty box quick reject
     if (box[0] > box[3] || box[1] > box[4] || box[2] > box[5]) return false;
 
@@ -574,7 +577,7 @@ export function intersectsTriangle3(box: Box3, a: Vec3, b: Vec3, c: Vec3): boole
 /**
  * Test intersection between axis-aligned bounding box and a sphere.
  */
-export function intersectsSphere(box: Box3, sphere: Sphere): boolean {
+export function intersectsSphere(box: RBox3, sphere: RSphere): boolean {
     const { center, radius } = sphere;
     const cx = center[0];
     const cy = center[1];
@@ -589,7 +592,7 @@ export function intersectsSphere(box: Box3, sphere: Sphere): boolean {
 /**
  * Test intersection between axis-aligned bounding box and plane.
  */
-export function intersectsPlane3(box: Box3, plane: Plane3): boolean {
+export function intersectsPlane3(box: RBox3, plane: RPlane3): boolean {
     const { normal, constant } = plane;
     const nx = normal[0];
     const ny = normal[1];
diff --git a/src/shapes/circle.ts b/src/shapes/circle.ts
index 5ad3459..d801e5a 100644
--- a/src/shapes/circle.ts
+++ b/src/shapes/circle.ts
@@ -1,8 +1,12 @@
+import type { DeepReadonly } from '../core/readonly';
 import type { Vec2 } from '../core/vec2';
 
 /** A circle in 2D space */
 export type Circle = { center: Vec2; radius: number };
 
+/** A read-only circle in 2D space */
+export type RCircle = DeepReadonly<Circle>;
+
 export function create(): Circle {
     return { center: [0, 0], radius: 0 };
 }
diff --git a/src/shapes/frustum.ts b/src/shapes/frustum.ts
index 737ece1..ffcf5cb 100644
--- a/src/shapes/frustum.ts
+++ b/src/shapes/frustum.ts
@@ -1,9 +1,10 @@
-import type { Mat4 } from '../core/mat4';
-import type { Vec3 } from '../core/vec3';
-import type { Box3 } from './box3';
+import type { RMat4 } from '../core/mat4';
+import type { DeepReadonly } from '../core/readonly';
+import type { RVec3, Vec3 } from '../core/vec3';
+import type { RBox3 } from './box3';
 import type { Plane3 } from './plane3';
 import * as plane3 from './plane3';
-import type { Sphere } from './sphere';
+import type { RSphere } from './sphere';
 
 /**
  * A view frustum, represented as the six bounding planes of a camera's view volume.
@@ -12,6 +13,9 @@ import type { Sphere } from './sphere';
  */
 export type Frustum = [Plane3, Plane3, Plane3, Plane3, Plane3, Plane3];
 
+/** A read-only view frustum */
+export type RFrustum = DeepReadonly<Frustum>;
+
 /**
  * The eight corners of a frustum, as returned by [[corners]].
  * Ordered near bottom-left, near top-left, near top-right, near bottom-right,
@@ -19,6 +23,9 @@ export type Frustum = [Plane3, Plane3, Plane3, Plane3, Plane3, Plane3];
  */
 export type FrustumCorners = [Vec3, Vec3, Vec3, Vec3, Vec3, Vec3, Vec3, Vec3];
 
+/** The eight read-only corners of a frustum */
+export type RFrustumCorners = DeepReadonly<FrustumCorners>;
+
 /**
  * Creates a new frustum of zeroed planes.
  * @returns A new frustum
@@ -39,7 +46,7 @@ export function create(): Frustum {
  * @param f - The frustum to clone
  * @returns A new frustum
  */
-export function clone(f: Frustum): Frustum {
+export function clone(f: RFrustum): Frustum {
     const p0 = f[0];
     const p1 = f[1];
     const p2 = f[2];
@@ -62,7 +69,7 @@ export function clone(f: Frustum): Frustum {
  * @param f - The source frustum
  * @returns The output frustum
  */
-export function copy(out: Frustum, f: Frustum): Frustum {
+export function copy(out: Frustum, f: RFrustum): Frustum {
     plane3.copy(out[0], f[0]);
     plane3.copy(out[1], f[1]);
     plane3.copy(out[2], f[2]);
@@ -83,7 +90,7 @@ export function copy(out: Frustum, f: Frustum): Frustum {
  * @param view - The view matrix
  * @returns The output frustum
  */
-export function setFromViewProjectionMatrixNO(out: Frustum, proj: Mat4, view: Mat4): Frustum {
+export function setFromViewProjectionMatrixNO(out: Frustum, proj: RMat4, view: RMat4): Frustum {
     const p0 = proj[0];
     const p1 = proj[1];
     const p2 = proj[2];
@@ -225,7 +232,7 @@ export function setFromViewProjectionMatrixNO(out: Frustum, proj: Mat4, view: Ma
  * @param view - The view matrix
  * @returns The output frustum
  */
-export function setFromViewProjectionMatrixZO(out: Frustum, proj: Mat4, view: Mat4): Frustum {
+export function setFromViewProjectionMatrixZO(out: Frustum, proj: RMat4, view: RMat4): Frustum {
     const p0 = proj[0];
     const p1 = proj[1];
     const p2 = proj[2];
@@ -370,7 +377,7 @@ export function setFromViewProjectionMatrixZO(out: Frustum, proj: Mat4, view: Ma
  * @param view - The view matrix
  * @returns The output frustum
  */
-export function setFromViewProjectionMatrixSides(out: Frustum, proj: Mat4, view: Mat4): Frustum {
+export function setFromViewProjectionMatrixSides(out: Frustum, proj: RMat4, view: RMat4): Frustum {
     // row2 (near/far) coefficients are not needed, so skip p2, p6, p10, p14
     const p0 = proj[0];
     const p1 = proj[1];
@@ -475,7 +482,7 @@ export function setFromViewProjectionMatrixSides(out: Frustum, proj: Mat4, view:
  * @param s - The sphere
  * @returns True if the sphere intersects or is inside the frustum
  */
-export function intersectsSphere(f: Frustum, s: Sphere): boolean {
+export function intersectsSphere(f: RFrustum, s: RSphere): boolean {
     const cx = s.center[0];
     const cy = s.center[1];
     const cz = s.center[2];
@@ -493,7 +500,7 @@ export function intersectsSphere(f: Frustum, s: Sphere): boolean {
  * @param s - The sphere
  * @returns True if the sphere intersects or is inside the frustum's sides
  */
-export function sidesIntersectsSphere(f: Frustum, s: Sphere): boolean {
+export function sidesIntersectsSphere(f: RFrustum, s: RSphere): boolean {
     const cx = s.center[0];
     const cy = s.center[1];
     const cz = s.center[2];
@@ -511,7 +518,7 @@ export function sidesIntersectsSphere(f: Frustum, s: Sphere): boolean {
  * @param box - The box
  * @returns True if the box intersects or is inside the frustum
  */
-export function intersectsBox3(f: Frustum, box: Box3): boolean {
+export function intersectsBox3(f: RFrustum, box: RBox3): boolean {
     const minX = box[0];
     const minY = box[1];
     const minZ = box[2];
@@ -537,7 +544,7 @@ export function intersectsBox3(f: Frustum, box: Box3): boolean {
  * @param box - The box
  * @returns True if the box intersects or is inside the frustum's sides
  */
-export function sidesIntersectsBox3(f: Frustum, box: Box3): boolean {
+export function sidesIntersectsBox3(f: RFrustum, box: RBox3): boolean {
     const minX = box[0];
     const minY = box[1];
     const minZ = box[2];
@@ -562,7 +569,7 @@ export function sidesIntersectsBox3(f: Frustum, box: Box3): boolean {
  * @param p - The point
  * @returns True if the point is inside or on the boundary of the frustum
  */
-export function containsPoint(f: Frustum, p: Vec3): boolean {
+export function containsPoint(f: RFrustum, p: RVec3): boolean {
     const x = p[0];
     const y = p[1];
     const z = p[2];
@@ -579,7 +586,7 @@ export function containsPoint(f: Frustum, p: Vec3): boolean {
  * @param p - The point
  * @returns True if the point is inside or on the boundary of the frustum's sides
  */
-export function sidesContainsPoint(f: Frustum, p: Vec3): boolean {
+export function sidesContainsPoint(f: RFrustum, p: RVec3): boolean {
     const x = p[0];
     const y = p[1];
     const z = p[2];
@@ -598,7 +605,7 @@ export function sidesContainsPoint(f: Frustum, p: Vec3): boolean {
  * @param direction - Ray direction (need not be normalized)
  * @returns True if the ray intersects the frustum
  */
-export function intersectsRay(f: Frustum, origin: Vec3, direction: Vec3): boolean {
+export function intersectsRay(f: RFrustum, origin: RVec3, direction: RVec3): boolean {
     const ox = origin[0];
     const oy = origin[1];
     const oz = origin[2];
@@ -636,7 +643,7 @@ export function intersectsRay(f: Frustum, origin: Vec3, direction: Vec3): boolea
  * @param direction - Ray direction (need not be normalized)
  * @returns True if the ray intersects the frustum's sides
  */
-export function sidesIntersectsRay(f: Frustum, origin: Vec3, direction: Vec3): boolean {
+export function sidesIntersectsRay(f: RFrustum, origin: RVec3, direction: RVec3): boolean {
     const ox = origin[0];
     const oy = origin[1];
     const oz = origin[2];
@@ -677,7 +684,7 @@ export function sidesIntersectsRay(f: Frustum, origin: Vec3, direction: Vec3): b
  * @param f - The frustum
  * @returns The output corners
  */
-export function corners(out: FrustumCorners, f: Frustum): FrustumCorners {
+export function corners(out: FrustumCorners, f: RFrustum): FrustumCorners {
     // near = f[4], far = f[5], left = f[0], right = f[1], bottom = f[2], top = f[3]
     plane3.intersect(out[0], f[4], f[0], f[2]);
     plane3.intersect(out[1], f[4], f[0], f[3]);
diff --git a/src/shapes/index.ts b/src/shapes/index.ts
index a81fe65..4b1f252 100644
--- a/src/shapes/index.ts
+++ b/src/shapes/index.ts
@@ -3,22 +3,22 @@
 export type * from '../core';
 
 export * as box2 from './box2';
-export type { Box2 } from './box2';
+export type { Box2, RBox2 } from './box2';
 
 export * as box3 from './box3';
-export type { Box3 } from './box3';
+export type { Box3, RBox3 } from './box3';
 
 export * as obb3 from './obb3';
-export type { OBB3 } from './obb3';
+export type { OBB3, ROBB3 } from './obb3';
 
 export * as plane3 from './plane3';
-export type { Plane3 } from './plane3';
+export type { Plane3, RPlane3 } from './plane3';
 
 export * as sphere from './sphere';
-export type { Sphere } from './sphere';
+export type { RSphere, Sphere } from './sphere';
 
 export * as circle from './circle';
-export type { Circle } from './circle';
+export type { Circle, RCircle } from './circle';
 
 export * as segment2 from './segment2';
 
@@ -31,4 +31,4 @@ export * as triangle3 from './triangle3';
 export * as raycast3 from './raycast3';
 
 export * as frustum from './frustum';
-export type { Frustum, FrustumCorners } from './frustum';
+export type { Frustum, FrustumCorners, RFrustum, RFrustumCorners } from './frustum';
diff --git a/src/shapes/obb3.ts b/src/shapes/obb3.ts
index 4184c3e..e054720 100644
--- a/src/shapes/obb3.ts
+++ b/src/shapes/obb3.ts
@@ -1,19 +1,23 @@
-import type { Mat3 } from '../core/mat3';
+import type { Mat3, RMat3 } from '../core/mat3';
 import * as mat3 from '../core/mat3';
-import type { Mat4 } from '../core/mat4';
-import type { Quat } from '../core/quat';
+import type { RMat4 } from '../core/mat4';
+import type { RQuat } from '../core/quat';
+import type { DeepReadonly } from '../core/readonly';
 import { EPSILON } from '../core/scalar';
-import type { Vec3 } from '../core/vec3';
-import type { Box3 } from './box3';
+import type { RVec3, Vec3 } from '../core/vec3';
+import type { RBox3 } from './box3';
 
 /** An oriented bounding box in 3D space */
 export type OBB3 = { center: Vec3; halfExtents: Vec3; rotation: Mat3 };
 
+/** A read-only oriented bounding box in 3D space */
+export type ROBB3 = DeepReadonly<OBB3>;
+
 export function create(): OBB3 {
     return { center: [0, 0, 0], halfExtents: [1, 1, 1], rotation: mat3.create() };
 }
 
-export function clone(a: OBB3): OBB3 {
+export function clone(a: ROBB3): OBB3 {
     return {
         center: [a.center[0], a.center[1], a.center[2]],
         halfExtents: [a.halfExtents[0], a.halfExtents[1], a.halfExtents[2]],
@@ -21,7 +25,7 @@ export function clone(a: OBB3): OBB3 {
     };
 }
 
-export function copy(out: OBB3, a: OBB3): OBB3 {
+export function copy(out: OBB3, a: ROBB3): OBB3 {
     out.center[0] = a.center[0];
     out.center[1] = a.center[1];
     out.center[2] = a.center[2];
@@ -48,7 +52,7 @@ export function copy(out: OBB3, a: OBB3): OBB3 {
  * @param rotation the Mat3 rotation matrix
  * @returns the OBB with the given center, half extents, and rotation
  */
-export function set(out: OBB3, center: Vec3, halfExtents: Vec3, rotation: Mat3): OBB3 {
+export function set(out: OBB3, center: RVec3, halfExtents: RVec3, rotation: RMat3): OBB3 {
     out.center[0] = center[0];
     out.center[1] = center[1];
     out.center[2] = center[2];
@@ -77,7 +81,7 @@ export function set(out: OBB3, center: Vec3, halfExtents: Vec3, rotation: Mat3):
  * @param q - The quaternion representing the OBB's orientation
  * @returns out
  */
-export function setFromCenterHalfExtentsQuaternion(out: OBB3, center: Vec3, halfExtents: Vec3, q: Quat): OBB3 {
+export function setFromCenterHalfExtentsQuaternion(out: OBB3, center: RVec3, halfExtents: RVec3, q: RQuat): OBB3 {
     out.center[0] = center[0];
     out.center[1] = center[1];
     out.center[2] = center[2];
@@ -97,7 +101,7 @@ export function setFromCenterHalfExtentsQuaternion(out: OBB3, center: Vec3, half
  * @param aabb - The AABB (min and max corners)
  * @returns out
  */
-export function setFromBox3(out: OBB3, aabb: Box3): OBB3 {
+export function setFromBox3(out: OBB3, aabb: RBox3): OBB3 {
     // Center = (min + max) / 2
     out.center[0] = (aabb[0] + aabb[3]) * 0.5;
     out.center[1] = (aabb[1] + aabb[4]) * 0.5;
@@ -121,7 +125,7 @@ export function setFromBox3(out: OBB3, aabb: Box3): OBB3 {
  * @param point - The point to test
  * @returns true if the point is inside the OBB
  */
-export function containsPoint(obb: OBB3, point: Vec3): boolean {
+export function containsPoint(obb: ROBB3, point: RVec3): boolean {
     // Vector from center to point
     const dx = point[0] - obb.center[0];
     const dy = point[1] - obb.center[1];
@@ -147,7 +151,7 @@ export function containsPoint(obb: OBB3, point: Vec3): boolean {
  * @param point - The point to clamp
  * @returns out
  */
-export function clampPoint(out: Vec3, obb: OBB3, point: Vec3): Vec3 {
+export function clampPoint(out: Vec3, obb: ROBB3, point: RVec3): Vec3 {
     // OBB axes are the columns of the rotation matrix, read directly from r[].
     const r = obb.rotation;
 
@@ -202,7 +206,7 @@ export function clampPoint(out: Vec3, obb: OBB3, point: Vec3): Vec3 {
  * @param epsilon - Squared-sine threshold below which near-parallel edge axes are skipped
  * @returns true if the OBBs intersect
  */
-export function intersectsOBB3(a: OBB3, b: OBB3, epsilon = EPSILON): boolean {
+export function intersectsOBB3(a: ROBB3, b: ROBB3, epsilon = EPSILON): boolean {
     const rotA = a.rotation;
     const rotB = b.rotation;
 
@@ -368,7 +372,7 @@ export function intersectsOBB3(a: OBB3, b: OBB3, epsilon = EPSILON): boolean {
  * @param aabb - The AABB (axis-aligned bounding box)
  * @returns true if they intersect
  */
-export function intersectsBox3(obb: OBB3, aabb: Box3): boolean {
+export function intersectsBox3(obb: ROBB3, aabb: RBox3): boolean {
     const rotA = obb.rotation;
     const epsilon = EPSILON;
 
@@ -514,7 +518,7 @@ export function intersectsBox3(obb: OBB3, aabb: Box3): boolean {
  * @param matrix - The 4x4 transformation matrix
  * @returns out
  */
-export function applyMatrix4(out: OBB3, obb: OBB3, matrix: Mat4): OBB3 {
+export function applyMatrix4(out: OBB3, obb: ROBB3, matrix: RMat4): OBB3 {
     // read the upper-left 3x3 (the affine linear part) into locals once. Columns
     // m0*, m1*, m2* correspond to mat4 columns 0, 1, 2.
     const m00 = matrix[0];
diff --git a/src/shapes/plane3.ts b/src/shapes/plane3.ts
index 009a28b..f538c9b 100644
--- a/src/shapes/plane3.ts
+++ b/src/shapes/plane3.ts
@@ -1,7 +1,8 @@
-import type { Mat4 } from '../core/mat4';
-import type { Vec3 } from '../core/vec3';
+import type { RMat4 } from '../core/mat4';
+import type { DeepReadonly } from '../core/readonly';
+import type { RVec3, Vec3 } from '../core/vec3';
 import * as vec3 from '../core/vec3';
-import type { Sphere } from './sphere';
+import type { RSphere } from './sphere';
 
 /**
  * A plane in 3D space
@@ -10,6 +11,9 @@ import type { Sphere } from './sphere';
  */
 export type Plane3 = { normal: Vec3; constant: number };
 
+/** A read-only plane in 3D space */
+export type RPlane3 = DeepReadonly<Plane3>;
+
 /**
  * Creates a new plane with normal (0, 1, 0) and constant 0
  * @returns A new plane
@@ -25,7 +29,7 @@ export function create(): Plane3 {
  * @param constant - The signed distance from origin
  * @returns The output plane
  */
-export function fromNormalAndConstant(out: Plane3, normal: Vec3, constant: number): Plane3 {
+export function fromNormalAndConstant(out: Plane3, normal: RVec3, constant: number): Plane3 {
     vec3.copy(out.normal, normal);
     out.constant = constant;
     return out;
@@ -38,7 +42,7 @@ export function fromNormalAndConstant(out: Plane3, normal: Vec3, constant: numbe
  * @param point - A point on the plane
  * @returns The output plane
  */
-export function fromNormalAndPoint(out: Plane3, normal: Vec3, point: Vec3): Plane3 {
+export function fromNormalAndPoint(out: Plane3, normal: RVec3, point: RVec3): Plane3 {
     vec3.copy(out.normal, normal);
     out.constant = -vec3.dot(normal, point);
     return out;
@@ -52,7 +56,7 @@ export function fromNormalAndPoint(out: Plane3, normal: Vec3, point: Vec3): Plan
  * @param c - Third point
  * @returns The output plane
  */
-export function fromCoplanarPoints(out: Plane3, a: Vec3, b: Vec3, c: Vec3): Plane3 {
+export function fromCoplanarPoints(out: Plane3, a: RVec3, b: RVec3, c: RVec3): Plane3 {
     const ax = a[0];
     const ay = a[1];
     const az = a[2];
@@ -90,7 +94,7 @@ export function fromCoplanarPoints(out: Plane3, a: Vec3, b: Vec3, c: Vec3): Plan
  * @param plane - The plane to clone
  * @returns A new plane
  */
-export function clone(plane: Plane3): Plane3 {
+export function clone(plane: RPlane3): Plane3 {
     return {
         normal: vec3.clone(plane.normal),
         constant: plane.constant,
@@ -103,7 +107,7 @@ export function clone(plane: Plane3): Plane3 {
  * @param plane - The source plane
  * @returns The output plane
  */
-export function copy(out: Plane3, plane: Plane3): Plane3 {
+export function copy(out: Plane3, plane: RPlane3): Plane3 {
     vec3.copy(out.normal, plane.normal);
     out.constant = plane.constant;
     return out;
@@ -115,7 +119,7 @@ export function copy(out: Plane3, plane: Plane3): Plane3 {
  * @param plane - The input plane
  * @returns The normalized plane
  */
-export function normalize(out: Plane3, plane: Plane3): Plane3 {
+export function normalize(out: Plane3, plane: RPlane3): Plane3 {
     const invMagnitude = 1.0 / vec3.length(plane.normal);
     vec3.scale(out.normal, plane.normal, invMagnitude);
     out.constant = plane.constant * invMagnitude;
@@ -128,7 +132,7 @@ export function normalize(out: Plane3, plane: Plane3): Plane3 {
  * @param plane - The input plane
  * @returns The negated plane
  */
-export function negate(out: Plane3, plane: Plane3): Plane3 {
+export function negate(out: Plane3, plane: RPlane3): Plane3 {
     vec3.negate(out.normal, plane.normal);
     out.constant = -plane.constant;
     return out;
@@ -141,7 +145,7 @@ export function negate(out: Plane3, plane: Plane3): Plane3 {
  * @param distance - The distance to offset (positive = in direction of normal)
  * @returns The offset plane
  */
-export function offset(out: Plane3, plane: Plane3, distance: number): Plane3 {
+export function offset(out: Plane3, plane: RPlane3, distance: number): Plane3 {
     vec3.copy(out.normal, plane.normal);
     out.constant = plane.constant - distance;
     return out;
@@ -153,7 +157,7 @@ export function offset(out: Plane3, plane: Plane3, distance: number): Plane3 {
  * @param point - The point
  * @returns The signed distance (positive = in direction of normal)
  */
-export function distanceToPoint(plane: Plane3, point: Vec3): number {
+export function distanceToPoint(plane: RPlane3, point: RVec3): number {
     return vec3.dot(plane.normal, point) + plane.constant;
 }
 
@@ -164,7 +168,7 @@ export function distanceToPoint(plane: Plane3, point: Vec3): number {
  * @param point - The point to project
  * @returns The projected point
  */
-export function projectPoint(out: Vec3, plane: Plane3, point: Vec3): Vec3 {
+export function projectPoint(out: Vec3, plane: RPlane3, point: RVec3): Vec3 {
     const distance = distanceToPoint(plane, point);
     return vec3.scaleAndAdd(out, point, plane.normal, -distance);
 }
@@ -176,7 +180,7 @@ export function projectPoint(out: Vec3, plane: Plane3, point: Vec3): Vec3 {
  * @param matrix - The transformation matrix
  * @returns The transformed plane
  */
-export function transform(out: Plane3, plane: Plane3, matrix: Mat4): Plane3 {
+export function transform(out: Plane3, plane: RPlane3, matrix: RMat4): Plane3 {
     // Transform the normal by the inverse transpose of the matrix
     // For a proper implementation, you'd need mat4.invert and proper normal transformation
     // This is a simplified version (rotation-only normal transform). fully scalar, no allocations.
@@ -224,7 +228,7 @@ export function transform(out: Plane3, plane: Plane3, matrix: Mat4): Plane3 {
  * @param sphere - The sphere
  * @returns True if they intersect
  */
-export function intersectsSphere(plane: Plane3, sphere: Sphere): boolean {
+export function intersectsSphere(plane: RPlane3, sphere: RSphere): boolean {
     const distance = Math.abs(distanceToPoint(plane, sphere.center));
     return distance <= sphere.radius;
 }
@@ -235,7 +239,7 @@ export function intersectsSphere(plane: Plane3, sphere: Sphere): boolean {
  * @param b - Second plane
  * @returns True if planes are exactly equal
  */
-export function exactEquals(a: Plane3, b: Plane3): boolean {
+export function exactEquals(a: RPlane3, b: RPlane3): boolean {
     return vec3.exactEquals(a.normal, b.normal) && a.constant === b.constant;
 }
 
@@ -247,7 +251,7 @@ export function exactEquals(a: Plane3, b: Plane3): boolean {
  * @param p3 - Third plane
  * @returns True if intersection exists, false if planes are degenerate or parallel
  */
-export function intersect(out: Vec3, p1: Plane3, p2: Plane3, p3: Plane3): boolean {
+export function intersect(out: Vec3, p1: RPlane3, p2: RPlane3, p3: RPlane3): boolean {
     // point = -(d1*(N2×N3) + d2*(N3×N1) + d3*(N1×N2)) / (N1·(N2×N3))
     // Cramer's rule: the three cross products are the columns of adj(M) and are reused
     // between the determinant and the numerator. fully scalar, zero allocations.
@@ -302,6 +306,6 @@ export function intersect(out: Vec3, p1: Plane3, p2: Plane3, p3: Plane3): boolea
  * @param b - Second plane
  * @returns True if planes are equal
  */
-export function equals(a: Plane3, b: Plane3): boolean {
+export function equals(a: RPlane3, b: RPlane3): boolean {
     return vec3.equals(a.normal, b.normal) && Math.abs(a.constant - b.constant) < 0.000001;
 }
diff --git a/src/shapes/polygon2.ts b/src/shapes/polygon2.ts
index 0369f7e..2aab075 100644
--- a/src/shapes/polygon2.ts
+++ b/src/shapes/polygon2.ts
@@ -1,4 +1,4 @@
-import type { Vec2 } from '../core/vec2';
+import type { RVec2, Vec2 } from '../core/vec2';
 import type { Box2 } from './box2';
 
 /**
@@ -20,7 +20,7 @@ import type { Box2 } from './box2';
  * @param n number of vertices to read from `vertices`
  * @returns the signed area
  */
-export function signedArea(vertices: number[], n: number): number {
+export function signedArea(vertices: readonly number[], n: number): number {
     let area = 0;
     for (let i = 0, j = n - 1; i < n; j = i++) {
         const xi = vertices[i * 2];
@@ -39,7 +39,7 @@ export function signedArea(vertices: number[], n: number): number {
  * @param n number of vertices to read from `vertices`
  * @returns the absolute area
  */
-export function area(vertices: number[], n: number): number {
+export function area(vertices: readonly number[], n: number): number {
     return Math.abs(signedArea(vertices, n));
 }
 
@@ -52,7 +52,7 @@ export function area(vertices: number[], n: number): number {
  * @param point the point to test
  * @returns true if the point is inside (or on the boundary of) the polygon
  */
-export function containsPoint(vertices: number[], n: number, point: Vec2): boolean {
+export function containsPoint(vertices: readonly number[], n: number, point: RVec2): boolean {
     let inside = false;
     const x = point[0];
     const y = point[1];
@@ -112,7 +112,7 @@ export function containsPoint(vertices: number[], n: number, point: Vec2): boole
  * @param n number of vertices to read from `vertices`
  * @returns out
  */
-export function centroid(out: Vec2, vertices: number[], n: number): Vec2 {
+export function centroid(out: Vec2, vertices: readonly number[], n: number): Vec2 {
     let cx = 0;
     let cy = 0;
     let a2 = 0; // twice the signed area
@@ -154,7 +154,7 @@ export function centroid(out: Vec2, vertices: number[], n: number): Vec2 {
  * @param n number of vertices to read from `vertices`
  * @returns the perimeter
  */
-export function perimeter(vertices: number[], n: number): number {
+export function perimeter(vertices: readonly number[], n: number): number {
     let total = 0;
     for (let i = 0, j = n - 1; i < n; j = i++) {
         const dx = vertices[i * 2] - vertices[j * 2];
@@ -173,7 +173,7 @@ export function perimeter(vertices: number[], n: number): number {
  * @param n number of vertices to read from `vertices`
  * @returns 1 (CCW), -1 (CW), or 0 (degenerate)
  */
-export function winding(vertices: number[], n: number): number {
+export function winding(vertices: readonly number[], n: number): number {
     const a = signedArea(vertices, n);
     if (a > 0) return 1;
     if (a < 0) return -1;
@@ -188,7 +188,7 @@ export function winding(vertices: number[], n: number): number {
  * @param n number of vertices to read from `vertices`
  * @returns true if the polygon is convex
  */
-export function isConvex(vertices: number[], n: number): boolean {
+export function isConvex(vertices: readonly number[], n: number): boolean {
     if (n < 3) return false;
 
     let sign = 0;
@@ -224,7 +224,7 @@ export function isConvex(vertices: number[], n: number): boolean {
  * @param i index of the vertex to test
  * @returns true if vertex `i` is reflex
  */
-export function isReflexVertex(vertices: number[], n: number, i: number): boolean {
+export function isReflexVertex(vertices: readonly number[], n: number, i: number): boolean {
     const p = ((i - 1 + n) % n) * 2;
     const c = i * 2;
     const q = ((i + 1) % n) * 2;
@@ -246,7 +246,7 @@ export function isReflexVertex(vertices: number[], n: number, i: number): boolea
  * @param n number of vertices to read from `vertices`
  * @returns out
  */
-export function reverse(out: number[], vertices: number[], n: number): number[] {
+export function reverse(out: number[], vertices: readonly number[], n: number): number[] {
     for (let lo = 0, hi = n - 1; lo < hi; lo++, hi--) {
         const x = vertices[lo * 2];
         const y = vertices[lo * 2 + 1];
@@ -273,7 +273,7 @@ export function reverse(out: number[], vertices: number[], n: number): number[]
  * @param n number of vertices to read from `vertices`
  * @returns out
  */
-export function bounds(out: Box2, vertices: number[], n: number): Box2 {
+export function bounds(out: Box2, vertices: readonly number[], n: number): Box2 {
     let minX = Number.POSITIVE_INFINITY;
     let minY = Number.POSITIVE_INFINITY;
     let maxX = Number.NEGATIVE_INFINITY;
@@ -306,7 +306,7 @@ export function bounds(out: Box2, vertices: number[], n: number): Box2 {
  * @param point the query point
  * @returns out
  */
-export function closestPoint(out: Vec2, vertices: number[], n: number, point: Vec2): Vec2 {
+export function closestPoint(out: Vec2, vertices: readonly number[], n: number, point: RVec2): Vec2 {
     const px = point[0];
     const py = point[1];
     let bestDistSq = Number.POSITIVE_INFINITY;
@@ -352,7 +352,7 @@ export function closestPoint(out: Vec2, vertices: number[], n: number, point: Ve
  * @param point the query point
  * @returns the signed distance (negative inside, positive outside)
  */
-export function signedDistance(vertices: number[], n: number, point: Vec2): number {
+export function signedDistance(vertices: readonly number[], n: number, point: RVec2): number {
     const px = point[0];
     const py = point[1];
     let bestDistSq = Number.POSITIVE_INFINITY;
@@ -382,7 +382,7 @@ const _projA: [number, number] = [0, 0];
 const _projB: [number, number] = [0, 0];
 
 /** Projects a polygon onto the axis (nx, ny), storing [min, max] in `out`. */
-function projectOntoAxis(out: [number, number], nx: number, ny: number, vertices: number[], n: number): void {
+function projectOntoAxis(out: [number, number], nx: number, ny: number, vertices: readonly number[], n: number): void {
     let min = nx * vertices[0] + ny * vertices[1];
     let max = min;
     for (let i = 1; i < n; i++) {
@@ -395,7 +395,7 @@ function projectOntoAxis(out: [number, number], nx: number, ny: number, vertices
 }
 
 /** Returns true if the two polygons are separated along any edge normal of A. */
-function separatedByEdgesOf(verticesA: number[], numA: number, verticesB: number[], numB: number): boolean {
+function separatedByEdgesOf(verticesA: readonly number[], numA: number, verticesB: readonly number[], numB: number): boolean {
     for (let i = 0, j = numA - 1; i < numA; j = i++) {
         // Normal of edge (a -> b); orientation does not matter for separation.
         const nx = verticesA[i * 2 + 1] - verticesA[j * 2 + 1];
@@ -420,7 +420,7 @@ function separatedByEdgesOf(verticesA: number[], numA: number, verticesB: number
  * @param numB number of vertices in the second polygon
  * @returns true if the polygons overlap
  */
-export function overlapConvex(verticesA: number[], numA: number, verticesB: number[], numB: number): boolean {
+export function overlapConvex(verticesA: readonly number[], numA: number, verticesB: readonly number[], numB: number): boolean {
     if (separatedByEdgesOf(verticesA, numA, verticesB, numB)) return false;
     if (separatedByEdgesOf(verticesB, numB, verticesA, numA)) return false;
     return true;
@@ -437,7 +437,7 @@ export function overlapConvex(verticesA: number[], numA: number, verticesB: numb
  * @param b end of the segment
  * @returns true if the segment intersects the polygon
  */
-export function intersectsSegment(vertices: number[], n: number, a: Vec2, b: Vec2): boolean {
+export function intersectsSegment(vertices: readonly number[], n: number, a: RVec2, b: RVec2): boolean {
     if (containsPoint(vertices, n, a) || containsPoint(vertices, n, b)) return true;
 
     const ax = a[0];
diff --git a/src/shapes/raycast3.ts b/src/shapes/raycast3.ts
index 10560b8..91a1498 100644
--- a/src/shapes/raycast3.ts
+++ b/src/shapes/raycast3.ts
@@ -1,5 +1,5 @@
-import type { Vec3 } from '../core/vec3';
-import type { Box3 } from './box3';
+import type { RVec3 } from '../core/vec3';
+import type { RBox3 } from './box3';
 
 /**
  * Result of a ray-triangle intersection test
@@ -39,12 +39,12 @@ export function createIntersectsTriangleResult(): IntersectsTriangleResult {
  */
 export function intersectsTriangle(
     out: IntersectsTriangleResult,
-    origin: Vec3,
-    direction: Vec3,
+    origin: RVec3,
+    direction: RVec3,
     length: number,
-    a: Vec3,
-    b: Vec3,
-    c: Vec3,
+    a: RVec3,
+    b: RVec3,
+    c: RVec3,
     backfaceCulling: boolean,
 ): void {
     // compute edge1 = b - a
@@ -165,7 +165,7 @@ export function intersectsTriangle(
  * @param aabb AABB to test against
  * @returns true if ray intersects the AABB, false otherwise
  */
-export function intersectsBox3(origin: Vec3, direction: Vec3, length: number, aabb: Box3): boolean {
+export function intersectsBox3(origin: RVec3, direction: RVec3, length: number, aabb: RBox3): boolean {
     let tmin = 0;
     let tmax = length;
 
diff --git a/src/shapes/segment2.ts b/src/shapes/segment2.ts
index 21970eb..088c4d6 100644
--- a/src/shapes/segment2.ts
+++ b/src/shapes/segment2.ts
@@ -1,4 +1,4 @@
-import type { Vec2 } from '../core/vec2';
+import type { RVec2, Vec2 } from '../core/vec2';
 
 /**
  * Calculates the closest point on a line segment to a given point
@@ -7,7 +7,7 @@ import type { Vec2 } from '../core/vec2';
  * @param a First endpoint of the segment
  * @param b Second endpoint of the segment
  */
-export function closestPoint(out: Vec2, point: Vec2, a: Vec2, b: Vec2): Vec2 {
+export function closestPoint(out: Vec2, point: RVec2, a: RVec2, b: RVec2): Vec2 {
     const pqx = b[0] - a[0];
     const pqz = b[1] - a[1];
     const dx = point[0] - a[0];
@@ -36,7 +36,7 @@ export function closestPoint(out: Vec2, point: Vec2, a: Vec2, b: Vec2): Vec2 {
  * @param d second endpoint of the second segment
  * @returns true if the segments intersect
  */
-export function intersects(a: Vec2, b: Vec2, c: Vec2, d: Vec2): boolean {
+export function intersects(a: RVec2, b: RVec2, c: RVec2, d: RVec2): boolean {
     const rx = b[0] - a[0];
     const ry = b[1] - a[1];
     const ex = d[0] - c[0];
@@ -65,7 +65,7 @@ export function intersects(a: Vec2, b: Vec2, c: Vec2, d: Vec2): boolean {
  * @param d second endpoint of the second segment
  * @returns out if the segments intersect, otherwise null
  */
-export function intersection(out: Vec2, a: Vec2, b: Vec2, c: Vec2, d: Vec2): Vec2 | null {
+export function intersection(out: Vec2, a: RVec2, b: RVec2, c: RVec2, d: RVec2): Vec2 | null {
     const rx = b[0] - a[0];
     const ry = b[1] - a[1];
     const ex = d[0] - c[0];
diff --git a/src/shapes/sphere.ts b/src/shapes/sphere.ts
index b426ac9..9f15b3d 100644
--- a/src/shapes/sphere.ts
+++ b/src/shapes/sphere.ts
@@ -1,8 +1,12 @@
-import type { Vec3 } from '../core/vec3';
+import type { DeepReadonly } from '../core/readonly';
+import type { RVec3, Vec3 } from '../core/vec3';
 
 /** A sphere in 3D space */
 export type Sphere = { center: Vec3; radius: number };
 
+/** A read-only sphere in 3D space */
+export type RSphere = DeepReadonly<Sphere>;
+
 /**
  * Creates a new sphere with a default center 0,0,0 and radius 1
  * @returns A new sphere.
@@ -18,7 +22,7 @@ export function create(): Sphere {
  * @param point the point to test
  * @returns true if the point is within the sphere's radius
  */
-export function containsPoint(sphere: Sphere, point: Vec3): boolean {
+export function containsPoint(sphere: RSphere, point: RVec3): boolean {
     const dx = point[0] - sphere.center[0];
     const dy = point[1] - sphere.center[1];
     const dz = point[2] - sphere.center[2];
diff --git a/src/shapes/triangle2.ts b/src/shapes/triangle2.ts
index 35ee55c..adbc5c7 100644
--- a/src/shapes/triangle2.ts
+++ b/src/shapes/triangle2.ts
@@ -1,4 +1,4 @@
-import type { Vec2 } from '../core/vec2';
+import type { RVec2, Vec2 } from '../core/vec2';
 import type { Box2 } from './box2';
 
 /**
@@ -11,7 +11,7 @@ import type { Box2 } from './box2';
  * @param c the third vertex of the triangle.
  * @returns the signed area.
  */
-export function signedArea(a: Vec2, b: Vec2, c: Vec2): number {
+export function signedArea(a: RVec2, b: RVec2, c: RVec2): number {
     return ((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) / 2;
 }
 
@@ -23,7 +23,7 @@ export function signedArea(a: Vec2, b: Vec2, c: Vec2): number {
  * @param c the third vertex of the triangle.
  * @returns the absolute area.
  */
-export function area(a: Vec2, b: Vec2, c: Vec2): number {
+export function area(a: RVec2, b: RVec2, c: RVec2): number {
     return Math.abs(signedArea(a, b, c));
 }
 
@@ -36,7 +36,7 @@ export function area(a: Vec2, b: Vec2, c: Vec2): number {
  * @param c the third vertex of the triangle.
  * @returns out.
  */
-export function centroid(out: Vec2, a: Vec2, b: Vec2, c: Vec2): Vec2 {
+export function centroid(out: Vec2, a: RVec2, b: RVec2, c: RVec2): Vec2 {
     out[0] = (a[0] + b[0] + c[0]) / 3;
     out[1] = (a[1] + b[1] + c[1]) / 3;
     return out;
@@ -51,7 +51,7 @@ export function centroid(out: Vec2, a: Vec2, b: Vec2, c: Vec2): Vec2 {
  * @param c the third vertex of the triangle.
  * @returns out.
  */
-export function bounds(out: Box2, a: Vec2, b: Vec2, c: Vec2): Box2 {
+export function bounds(out: Box2, a: RVec2, b: RVec2, c: RVec2): Box2 {
     out[0] = Math.min(a[0], b[0], c[0]);
     out[1] = Math.min(a[1], b[1], c[1]);
     out[2] = Math.max(a[0], b[0], c[0]);
@@ -69,7 +69,7 @@ export function bounds(out: Box2, a: Vec2, b: Vec2, c: Vec2): Box2 {
  * @param point the point to test.
  * @returns true if the point is inside (or on the boundary of) the triangle.
  */
-export function containsPoint(a: Vec2, b: Vec2, c: Vec2, point: Vec2): boolean {
+export function containsPoint(a: RVec2, b: RVec2, c: RVec2, point: RVec2): boolean {
     const px = point[0];
     const py = point[1];
 
diff --git a/src/shapes/triangle3.ts b/src/shapes/triangle3.ts
index e3ad375..42343cf 100644
--- a/src/shapes/triangle3.ts
+++ b/src/shapes/triangle3.ts
@@ -1,4 +1,4 @@
-import type { Vec3 } from '../core/vec3';
+import type { RVec3, Vec3 } from '../core/vec3';
 import type { Box3 } from './box3';
 
 /**
@@ -9,7 +9,7 @@ import type { Box3 } from './box3';
  * @param c the third vertex of the triangle.
  * @returns the output box containing the axis-aligned bounding box of the triangle.
  */
-export function bounds(out: Box3, a: Vec3, b: Vec3, c: Vec3): Box3 {
+export function bounds(out: Box3, a: RVec3, b: RVec3, c: RVec3): Box3 {
     out[0] = Math.min(a[0], b[0], c[0]);
     out[1] = Math.min(a[1], b[1], c[1]);
     out[2] = Math.min(a[2], b[2], c[2]);
@@ -29,7 +29,7 @@ export function bounds(out: Box3, a: Vec3, b: Vec3, c: Vec3): Box3 {
  * @param c the third vertex of the triangle.
  * @returns the output vector containing the normal of the triangle.
  */
-export function normal(out: Vec3, a: Vec3, b: Vec3, c: Vec3): Vec3 {
+export function normal(out: Vec3, a: RVec3, b: RVec3, c: RVec3): Vec3 {
     const abx = b[0] - a[0];
     const aby = b[1] - a[1];
     const abz = b[2] - a[2];
@@ -61,7 +61,7 @@ export function normal(out: Vec3, a: Vec3, b: Vec3, c: Vec3): Vec3 {
  * @param c the third vertex of the triangle.
  * @returns the output vector containing the centroid of the triangle.
  */
-export function centroid(out: Vec3, a: Vec3, b: Vec3, c: Vec3): Vec3 {
+export function centroid(out: Vec3, a: RVec3, b: RVec3, c: RVec3): Vec3 {
     out[0] = (a[0] + b[0] + c[0]) / 3;
     out[1] = (a[1] + b[1] + c[1]) / 3;
     out[2] = (a[2] + b[2] + c[2]) / 3;
diff --git a/src/time/spring2.ts b/src/time/spring2.ts
index 2a3224d..e3b3bfd 100644
--- a/src/time/spring2.ts
+++ b/src/time/spring2.ts
@@ -1,17 +1,23 @@
-import type { Vec2 } from '../core';
+import type { RVec2, Vec2 } from '../core';
 import { coef, coefficients, type Spring } from './spring-core';
 
 // A Vec2 spring. See `spring` (scalar) for the model; the same four coefficients
 // are applied to each component. `damp` is `update` pinned to dampingRatio = 1.
 
 /** Creates a Vec2 spring at `value` (copied), at rest. */
-export const create = (value: Vec2 = [0, 0]): Spring<Vec2> => ({ value: [value[0], value[1]], velocity: [0, 0] });
+export const create = (value: RVec2 = [0, 0]): Spring<Vec2> => ({ value: [value[0], value[1]], velocity: [0, 0] });
 
 /**
  * Springs `state.value` toward `target`, mutating `state` in place. Returns it.
  * @param dampingRatio 1 = critically damped (no overshoot), <1 bouncy, >1 sluggish
  */
-export function update(state: Spring<Vec2>, target: Vec2, smoothTime: number, dampingRatio: number, delta: number): Spring<Vec2> {
+export function update(
+    state: Spring<Vec2>,
+    target: RVec2,
+    smoothTime: number,
+    dampingRatio: number,
+    delta: number,
+): Spring<Vec2> {
     coefficients(smoothTime, dampingRatio, delta);
     const val = state.value;
     const vel = state.velocity;
@@ -27,6 +33,6 @@ export function update(state: Spring<Vec2>, target: Vec2, smoothTime: number, da
 }
 
 /** Critically-damped Vec2 spring (dampingRatio = 1). See {@link update}. */
-export function damp(state: Spring<Vec2>, target: Vec2, smoothTime: number, delta: number): Spring<Vec2> {
+export function damp(state: Spring<Vec2>, target: RVec2, smoothTime: number, delta: number): Spring<Vec2> {
     return update(state, target, smoothTime, 1, delta);
 }
diff --git a/src/time/spring3.ts b/src/time/spring3.ts
index 8b7211a..906c7a5 100644
--- a/src/time/spring3.ts
+++ b/src/time/spring3.ts
@@ -1,11 +1,11 @@
-import type { Vec3 } from '../core';
+import type { RVec3, Vec3 } from '../core';
 import { coef, coefficients, type Spring } from './spring-core';
 
 // A Vec3 spring. See `spring` (scalar) for the model. `damp` is `update` pinned
 // to dampingRatio = 1.
 
 /** Creates a Vec3 spring at `value` (copied), at rest. */
-export const create = (value: Vec3 = [0, 0, 0]): Spring<Vec3> => ({
+export const create = (value: RVec3 = [0, 0, 0]): Spring<Vec3> => ({
     value: [value[0], value[1], value[2]],
     velocity: [0, 0, 0],
 });
@@ -14,7 +14,13 @@ export const create = (value: Vec3 = [0, 0, 0]): Spring<Vec3> => ({
  * Springs `state.value` toward `target`, mutating `state` in place. Returns it.
  * @param dampingRatio 1 = critically damped (no overshoot), <1 bouncy, >1 sluggish
  */
-export function update(state: Spring<Vec3>, target: Vec3, smoothTime: number, dampingRatio: number, delta: number): Spring<Vec3> {
+export function update(
+    state: Spring<Vec3>,
+    target: RVec3,
+    smoothTime: number,
+    dampingRatio: number,
+    delta: number,
+): Spring<Vec3> {
     coefficients(smoothTime, dampingRatio, delta);
     const val = state.value;
     const vel = state.velocity;
@@ -28,6 +34,6 @@ export function update(state: Spring<Vec3>, target: Vec3, smoothTime: number, da
 }
 
 /** Critically-damped Vec3 spring (dampingRatio = 1). See {@link update}. */
-export function damp(state: Spring<Vec3>, target: Vec3, smoothTime: number, delta: number): Spring<Vec3> {
+export function damp(state: Spring<Vec3>, target: RVec3, smoothTime: number, delta: number): Spring<Vec3> {
     return update(state, target, smoothTime, 1, delta);
 }
diff --git a/src/time/spring4.ts b/src/time/spring4.ts
index c7a6401..264c781 100644
--- a/src/time/spring4.ts
+++ b/src/time/spring4.ts
@@ -1,11 +1,11 @@
-import type { Vec4 } from '../core';
+import type { RVec4, Vec4 } from '../core';
 import { coef, coefficients, type Spring } from './spring-core';
 
 // A Vec4 spring. See `spring` (scalar) for the model. `damp` is `update` pinned
 // to dampingRatio = 1.
 
 /** Creates a Vec4 spring at `value` (copied), at rest. */
-export const create = (value: Vec4 = [0, 0, 0, 0]): Spring<Vec4> => ({
+export const create = (value: RVec4 = [0, 0, 0, 0]): Spring<Vec4> => ({
     value: [value[0], value[1], value[2], value[3]],
     velocity: [0, 0, 0, 0],
 });
@@ -14,7 +14,13 @@ export const create = (value: Vec4 = [0, 0, 0, 0]): Spring<Vec4> => ({
  * Springs `state.value` toward `target`, mutating `state` in place. Returns it.
  * @param dampingRatio 1 = critically damped (no overshoot), <1 bouncy, >1 sluggish
  */
-export function update(state: Spring<Vec4>, target: Vec4, smoothTime: number, dampingRatio: number, delta: number): Spring<Vec4> {
+export function update(
+    state: Spring<Vec4>,
+    target: RVec4,
+    smoothTime: number,
+    dampingRatio: number,
+    delta: number,
+): Spring<Vec4> {
     coefficients(smoothTime, dampingRatio, delta);
     const val = state.value;
     const vel = state.velocity;
@@ -28,6 +34,6 @@ export function update(state: Spring<Vec4>, target: Vec4, smoothTime: number, da
 }
 
 /** Critically-damped Vec4 spring (dampingRatio = 1). See {@link update}. */
-export function damp(state: Spring<Vec4>, target: Vec4, smoothTime: number, delta: number): Spring<Vec4> {
+export function damp(state: Spring<Vec4>, target: RVec4, smoothTime: number, delta: number): Spring<Vec4> {
     return update(state, target, smoothTime, 1, delta);
 }