From 7091313fa8f4a6a0ebfb236033b2de11f3904d81 Mon Sep 17 00:00:00 2001
From: Damien Seguin <dmnsgn@gmail.com>
Date: Fri, 12 Apr 2024 23:24:41 +0100
Subject: [PATCH] feat: use correct component order for random.quat

- add link to Graphics Gems III source
---
 index.js      | 25 ++++++++++++++++---------
 test/index.js | 21 ++++++++++++++++-----
 2 files changed, 32 insertions(+), 14 deletions(-)

diff --git a/index.js b/index.js
index a4ef489..4a82662 100755
--- a/index.js
+++ b/index.js
@@ -144,21 +144,28 @@ class Random {
 
   /**
    * Get a random quaternion.
+   * @see [Graphics Gems III, Edited by David Kirk, III.6 UNIFORM RANDOM ROTATIONS]
    * @see [Steve LaValle]{@link https://web.archive.org/web/20211105205926/http://planning.cs.uiuc.edu/node198.html}
    * @returns {import("pex-math").quat}
    */
   quat() {
-    const u1 = this.rng();
-    const sqrt1MinU = Math.sqrt(1 - u1);
-    const sqrtU = Math.sqrt(u1);
-    const u2 = 2 * Math.PI * this.rng();
-    const u3 = 2 * Math.PI * this.rng();
+    // Let X0, X1, and X2 be three independent random variables that are uniformly distributed between 0 and 1.
+    const x0 = this.rng();
 
+    // Compute two uniformly distributed angles, θ1 = 2πX1 and θ2 = 2πX2, and their sines and cosines, s1, c1, s2, c2
+    const theta1 = 2 * Math.PI * this.rng();
+    const theta2 = 2 * Math.PI * this.rng();
+
+    // Also compute r1 = sqrt(1 – X0) and r2 = sqrt(X0)
+    const r1 = Math.sqrt(1 - x0);
+    const r2 = Math.sqrt(x0);
+
+    // Then return the unit quaternion with components [s1 r1, c1 r1, s2 r2, c2 r2]
     return [
-      sqrt1MinU * Math.cos(u2),
-      sqrtU * Math.sin(u3),
-      sqrtU * Math.cos(u3),
-      sqrt1MinU * Math.sin(u2),
+      Math.sin(theta1) * r1,
+      Math.cos(theta1) * r1,
+      Math.sin(theta2) * r2,
+      Math.cos(theta2) * r2,
     ];
   }
 
diff --git a/test/index.js b/test/index.js
index 26fad02..08c6d8f 100644
--- a/test/index.js
+++ b/test/index.js
@@ -2,7 +2,7 @@ import { describe, it } from "node:test";
 import assert from "node:assert";
 
 import { aabb, rect } from "pex-geom";
-import { vec2, vec3 } from "pex-math";
+import { vec2, vec3, quat } from "pex-math";
 
 import random from "../index.js";
 
@@ -110,7 +110,7 @@ describe("float(max)", () => {
     avg /= ITERATIONS;
     assert(
       avg > (max * 1) / 4 && avg < (max * 3) / 4,
-      `Avg:${avg} Min: 0 Max:${max}`
+      `Avg:${avg} Min: 0 Max:${max}`,
     );
   });
 });
@@ -163,7 +163,7 @@ describe("int(max)", () => {
     avg /= ITERATIONS;
     assert(
       avg > (max * 1) / 4 && avg < (max * 3) / 4,
-      `Avg:${avg} Min: 0 Max:${max}`
+      `Avg:${avg} Min: 0 Max:${max}`,
     );
   });
 });
@@ -250,12 +250,12 @@ describe("vec2InRect(rect)", () => {
 
     for (let i = 0; i < ITERATIONS; i++) {
       assert(
-        rect.containsPoint(POSITIVE_RECT, random.vec2InRect(POSITIVE_RECT))
+        rect.containsPoint(POSITIVE_RECT, random.vec2InRect(POSITIVE_RECT)),
       );
     }
     for (let i = 0; i < ITERATIONS; i++) {
       assert(
-        rect.containsPoint(NEGATIVE_RECT, random.vec2InRect(NEGATIVE_RECT))
+        rect.containsPoint(NEGATIVE_RECT, random.vec2InRect(NEGATIVE_RECT)),
       );
     }
     for (let i = 0; i < ITERATIONS; i++) {
@@ -306,6 +306,17 @@ describe("chance(probability)", () => {
     assert.equal(wins, ITERATIONS);
   });
 });
+describe("quat()", () => {
+  it("should randomize", () => {
+    const q = quat.create();
+
+    for (let i = 0; i < ITERATIONS; i++) {
+      quat.identity(q);
+      quat.set(q, random.quat());
+      assert(quat.length(q) > Number.EPSILON);
+    }
+  });
+});
 describe("element(list)", () => {
   it("should return element from the list", () => {
     const list = ["a", "b", "c", "d", "e", "f"];