Skip to content

Turning a vector, live

A vector turned by a body's orientation from the body's frame into the world's, R(q) v, q a unit quaternion (x, y, z, w): with u its vector part and t = 2 (u × v), R(q) v = v + w t + u × t (lusine.kinematics.rotate). The render turns a shape's layer gradients back to the world by it. Below, each step draws an orientation and a vector at random and the library's turn turns the one by the other on your GPU as you read; beside it the program turns the same vector by the rotation matrix the quaternion gives, entry by entry, a different sum of different products.

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows

  • The strip above the controls is green while every check holds, and turns red for good once one fails.
  • The turn: the chart draws, a point a step, how far the library's turn lands from the matrix's, over the vector's length, in binary32's unit roundoffs (2⁻²⁴), beside the bound it is held to, 32. Below it, the step, the turned vector's components beside the matrix's, and the worst miss so far.
  • The checks: the matrix's turn, the turn within 32 unit roundoffs of the matrix's; length kept, a turn keeping the vector's length within 32 unit roundoffs; q and -q alike, the orientation and its negative giving the same turn within a micrometre; a quarter turn, the turn by (0, 0, √½, √½) taking (x, y, z) to (−y, x, z) within 4 unit roundoffs of the length; the identity, (0, 0, 0, 1) leaving the vector as it is, bit for bit. Each is held at every step.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The seed the turns are drawn from is the program's setting, which the host may write; this page offers no control to change it yet.

The program

The step's turn, what the panels show and the checks, the panels the page draws, and the one procedure, main, run each step: an orientation and a vector drawn, turned by the law, by the matrix, by the negative, a quarter and the identity, and the checks.

// The turn's live page: a vector turned each step by an orientation drawn at random, by the library's turn, beside the
// same turn through the rotation matrix the orientation gives, and the turn held to keeping lengths, to a quaternion
// and its negative turning alike, to a quarter turn about z and to the identity, with the page's checks computed here
// on the GPU and shown by the panels the program declares.

namespace rotate;

lusine.version version = { { number = 1 } };

// The reader's settings: the seed the orientations and vectors are drawn from, and the quarter turn about +z, (0, 0,
// a, a) with a = sqrt(1/2), as its component's square.
host struct { integer seed = 2026; float half = 0.5; } settings;

// The step's orientation, a unit quaternion (x, y, z, w) drawn at random; the vector it turns; the library's turn of
// it, the matrix's and the turn by the orientation's negative; the quarter turn about +z of the vector, and the
// identity's.
struct {
  integer step; vector<float, 4> orientation; vector<float<m>, 3> value; vector<float<m>, 3> turned;
  vector<float<m>, 3> by_matrix; vector<float<m>, 3> negated; vector<float<m>, 3> quarter; vector<float<m>, 3> same;
} turn;

// What the panels show: the turned vector's components, the library's and the matrix's; the step's miss between the
// two over the vector's length, in binary32's unit roundoffs (2^-24), its worst so far, and the bound it is held to.
struct {
  float<m> x; float<m> y; float<m> z; float<m> matrix_x; float<m> matrix_y; float<m> matrix_z; float miss;
  float worst; float bound = 32.0;
} shown;

// The checks: the turn within 32 unit roundoffs of the matrix's, over the vector's length; the vector's length kept
// within 32 unit roundoffs; the orientation and its negative giving the same turn within 0.000001 m; a quarter turn
// about +z taking (x, y, z) to (-y, x, z) within 4 unit roundoffs of the length; and the identity leaving the vector
// as it is, bit for bit. Each at every step.
struct { bool matrix = true; bool length = true; bool negative = true; bool quarter = true; bool identity = true; }
  checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

float drawn(integer seed, integer step, integer index)
{
  return lusine.error.uniform(drawn = lusine.error.draw(seed = bits(seed), step = bits(step), key = fnv1a("turn"),
                                                        index = bits(index)));
}

// An orientation at random: four components each in [-1, 1), divided by their length.
vector<float, 4> orientation(integer seed, integer step)
{
  let raw = vector<float, 4>(2.0 * drawn(seed, step, index = 0) - 1.0, 2.0 * drawn(seed, step, index = 1) - 1.0,
                             2.0 * drawn(seed, step, index = 2) - 1.0, 2.0 * drawn(seed, step, index = 3) - 1.0);
  return raw / length(raw);
}

// A vector at random, each component in [-2, 2) m.
vector<float<m>, 3> value(integer seed, integer step)
{
  return vector<float<m>, 3>(4.0 m * (drawn(seed, step, index = 4) - 0.5), 4.0 m * (drawn(seed, step, index = 5) - 0.5),
                             4.0 m * (drawn(seed, step, index = 6) - 0.5));
}

// The rotation matrix of a unit quaternion (Shoemake, 1985) times the vector, row by row.
vector<float<m>, 3> by_matrix(vector<float, 4> q, vector<float<m>, 3> v)
{
  let x = q.x;
  let y = q.y;
  let z = q.z;
  let w = q.w;
  return vector<float<m>, 3>(
    (1.0 - 2.0 * (y * y + z * z)) * v.x + 2.0 * (x * y - z * w) * v.y + 2.0 * (x * z + y * w) * v.z,
    2.0 * (x * y + z * w) * v.x + (1.0 - 2.0 * (x * x + z * z)) * v.y + 2.0 * (y * z - x * w) * v.z,
    2.0 * (x * z - y * w) * v.x + 2.0 * (y * z + x * w) * v.y + (1.0 - 2.0 * (x * x + y * y)) * v.z);
}

vector<float, 4> quarter_z(float half) { return vector<float, 4>(0.0, 0.0, sqrt(half), sqrt(half)); }
vector<float, 4> identity() { return vector<float, 4>(0.0, 0.0, 0.0, 1.0); }

// The vector turned a quarter about +z by hand: (x, y, z) to (-y, x, z).
vector<float<m>, 3> quarter_by_hand(vector<float<m>, 3> v) { return vector<float<m>, 3>(-v.y, v.x, v.z); }

// A miss over a length, in unit roundoffs.
float roundoffs(vector<float<m>, 3> got, vector<float<m>, 3> want, vector<float<m>, 3> value)
{
  return length(got - want) / length(value) * 16777216.0;
}

bool kept(vector<float<m>, 3> turned, vector<float<m>, 3> value)
{
  return abs(length(turned) / length(value) - 1.0) <= 32.0 / 16777216.0;
}

bool same(vector<float<m>, 3> a, vector<float<m>, 3> b) { return a.x == b.x && a.y == b.y && a.z == b.z; }

output table<lusine.widget> panels = {
  lusine.series { { panel = "the turn", chart = "miss, unit roundoffs", quantity = shown.miss, steps = 200,
                    label = "turn against matrix" } },
  lusine.series { { panel = "the turn", chart = "miss, unit roundoffs", quantity = shown.bound, steps = 200,
                    label = "its bound" } },
  lusine.value { { panel = "the turn", quantity = turn.step, label = "step" } },
  lusine.value { { panel = "the turn", quantity = shown.x, label = "turned x" } },
  lusine.value { { panel = "the turn", quantity = shown.matrix_x, label = "matrix x" } },
  lusine.value { { panel = "the turn", quantity = shown.y, label = "turned y" } },
  lusine.value { { panel = "the turn", quantity = shown.matrix_y, label = "matrix y" } },
  lusine.value { { panel = "the turn", quantity = shown.z, label = "turned z" } },
  lusine.value { { panel = "the turn", quantity = shown.matrix_z, label = "matrix z" } },
  lusine.value { { panel = "the turn", quantity = shown.worst, label = "worst miss, roundoffs" } },
  lusine.value { { panel = "the checks", quantity = checks.matrix, label = "the matrix's turn" } },
  lusine.value { { panel = "the checks", quantity = checks.length, label = "length kept" } },
  lusine.value { { panel = "the checks", quantity = checks.negative, label = "q and -q alike" } },
  lusine.value { { panel = "the checks", quantity = checks.quarter, label = "a quarter turn" } },
  lusine.value { { panel = "the checks", quantity = checks.identity, label = "the identity" } },
};

void main()
{
  turn.step = turn.step + 1;
  turn.orientation = orientation(settings.seed, turn.step);
  turn.value = value(settings.seed, turn.step);
  turn.turned = lusine.kinematics.rotate(turn.orientation, turn.value);
  turn.by_matrix = by_matrix(q = turn.orientation, v = turn.value);
  turn.negated = lusine.kinematics.rotate(orientation = -turn.orientation, turn.value);
  turn.quarter = lusine.kinematics.rotate(orientation = quarter_z(settings.half), turn.value);
  turn.same = lusine.kinematics.rotate(orientation = identity(), turn.value);
  shown.x = turn.turned.x;
  shown.y = turn.turned.y;
  shown.z = turn.turned.z;
  shown.matrix_x = turn.by_matrix.x;
  shown.matrix_y = turn.by_matrix.y;
  shown.matrix_z = turn.by_matrix.z;
  shown.miss = roundoffs(got = turn.turned, want = turn.by_matrix, turn.value);
  shown.worst = shown.miss > shown.worst ? shown.miss : shown.worst;
  checks.matrix = checks.matrix && shown.miss <= shown.bound;
  checks.length = checks.length && kept(turn.turned, turn.value);
  checks.negative = checks.negative && length(turn.negated - turn.turned) <= 0.000001 m;
  checks.quarter = checks.quarter
                   && roundoffs(got = turn.quarter, want = quarter_by_hand(v = turn.value), turn.value) <= 4.0;
  checks.identity = checks.identity && same(a = turn.same, b = turn.value);
  lamps.colour = lamp_colour(holds = checks.matrix && checks.length && checks.negative && checks.quarter
                                     && checks.identity);
}

The turn's other proof (four hundred turns against a reference through the matrix, and the same invariants) is on its proof page.