Skip to content

Turning a vector back, live

A vector turned back by a body's orientation, from the world's frame into the body's, R(q)ᵀ v: the turn by q's conjugate (lusine.kinematics.unrotate). The render evaluates a shape's layers at the hit point less its centre turned back so, and they turn with the ball. Below, each step draws an orientation and two vectors at random and the library's turn back runs on your GPU as you read; beside it the program turns the same vector back by the transpose of the rotation matrix the quaternion gives, entry by entry, and turns a vector into the world by the library's turn (lusine.kinematics.rotate) and back again.

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 back: the chart draws, a point a step, how far the turn back lands from the transpose's, and how far a vector turned there and back lands from itself, each over the vector's length in binary32's unit roundoffs (2⁻²⁴), beside the bounds they are held to, 32 and 64. Below it, the step, the turned-back vector's components beside the transpose's, and the two sides of the adjoint, R a · b and a · Rᵀ b.
  • The checks: the transpose's turn, the turn back within 32 unit roundoffs of the transpose's; undoes the turn, a vector turned into the world and back within 64 unit roundoffs of itself; the turn's adjoint, R a · b equal to a · Rᵀ b within 64 unit roundoffs of |a| |b|, which holds of the transpose and of nothing else; 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 turns, what the panels show and the checks, the panels the page draws, and the one procedure, main, run each step: an orientation and two vectors drawn, turned back by the law and by the transpose, there and back, and the checks.

// The turn back's live page: a vector turned back each step by an orientation drawn at random, by the library's turn
// back, beside the same vector through the transpose of the rotation matrix the orientation gives, and the turn back
// held to undoing the library's turn, to being the turn's adjoint and to the identity, with the page's checks computed
// here on the GPU and shown by the panels the program declares.

namespace unrotate;

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

// The reader's settings: the seed the orientations and vectors are drawn from.
host struct { integer seed = 2027; } settings;

// The step's orientation, a unit quaternion (x, y, z, w) drawn at random; two vectors, a and b; a turned back by the
// library, and by the matrix's transpose; a turned into the world by the library's turn and back again; b turned back
// and a turned into the world, for the adjoint; and a turned back by the identity.
struct {
  integer step; vector<float, 4> orientation; vector<float<m>, 3> a; vector<float<m>, 3> b;
  vector<float<m>, 3> back; vector<float<m>, 3> by_transpose; vector<float<m>, 3> there_and_back;
  vector<float<m>, 3> b_back; vector<float<m>, 3> a_there; vector<float<m>, 3> same;
} turn;

// What the panels show: the turned-back vector's components, the library's and the transpose's; the step's misses
// over the vector's length in binary32's unit roundoffs (2^-24), against the transpose and there and back; and the
// bounds they are held to.
struct {
  float<m> x; float<m> y; float<m> z; float<m> transpose_x; float<m> transpose_y; float<m> transpose_z;
  float transpose_miss; float round_trip_miss; float transpose_bound = 32.0; float round_trip_bound = 64.0;
  float<m2> adjoint_left; float<m2> adjoint_right;
} shown;

// The checks: the turn back within 32 unit roundoffs of the transpose's, over the vector's length; a vector turned into
// the world by the library's turn and back again within 64 unit roundoffs of itself; the turn back the turn's
// adjoint, R a . b equal to a . R^T b within 64 unit roundoffs of |a| |b|; and the identity leaving the vector as it
// is, bit for bit. Each at every step.
struct { bool transpose = true; bool round_trip = true; bool adjoint = 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, from three draws from the first index on.
vector<float<m>, 3> value(integer seed, integer step, integer first)
{
  return vector<float<m>, 3>(4.0 m * (drawn(seed, step, index = first) - 0.5),
                             4.0 m * (drawn(seed, step, index = first + 1) - 0.5),
                             4.0 m * (drawn(seed, step, index = first + 2) - 0.5));
}

// The transpose of a unit quaternion's rotation matrix (Shoemake, 1985) times the vector, column by column.
vector<float<m>, 3> by_transpose(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> identity() { return vector<float, 4>(0.0, 0.0, 0.0, 1.0); }

// 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 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 back", chart = "miss, unit roundoffs", quantity = shown.transpose_miss,
                    steps = 200, label = "against the transpose" } },
  lusine.series { { panel = "the turn back", chart = "miss, unit roundoffs", quantity = shown.round_trip_miss,
                    steps = 200, label = "there and back" } },
  lusine.series { { panel = "the turn back", chart = "miss, unit roundoffs", quantity = shown.transpose_bound,
                    steps = 200, label = "bound, transpose" } },
  lusine.series { { panel = "the turn back", chart = "miss, unit roundoffs", quantity = shown.round_trip_bound,
                    steps = 200, label = "bound, there and back" } },
  lusine.value { { panel = "the turn back", quantity = turn.step, label = "step" } },
  lusine.value { { panel = "the turn back", quantity = shown.x, label = "turned back x" } },
  lusine.value { { panel = "the turn back", quantity = shown.transpose_x, label = "transpose x" } },
  lusine.value { { panel = "the turn back", quantity = shown.y, label = "turned back y" } },
  lusine.value { { panel = "the turn back", quantity = shown.transpose_y, label = "transpose y" } },
  lusine.value { { panel = "the turn back", quantity = shown.z, label = "turned back z" } },
  lusine.value { { panel = "the turn back", quantity = shown.transpose_z, label = "transpose z" } },
  lusine.value { { panel = "the turn back", quantity = shown.adjoint_left, label = "R a . b" } },
  lusine.value { { panel = "the turn back", quantity = shown.adjoint_right, label = "a . R^T b" } },
  lusine.value { { panel = "the checks", quantity = checks.transpose, label = "the transpose's turn" } },
  lusine.value { { panel = "the checks", quantity = checks.round_trip, label = "undoes the turn" } },
  lusine.value { { panel = "the checks", quantity = checks.adjoint, label = "the turn's adjoint" } },
  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.a = value(settings.seed, turn.step, first = 4);
  turn.b = value(settings.seed, turn.step, first = 7);
  turn.back = lusine.kinematics.unrotate(turn.orientation, value = turn.a);
  turn.by_transpose = by_transpose(q = turn.orientation, v = turn.a);
  turn.a_there = lusine.kinematics.rotate(turn.orientation, value = turn.a);
  turn.there_and_back = lusine.kinematics.unrotate(turn.orientation, value = turn.a_there);
  turn.b_back = lusine.kinematics.unrotate(turn.orientation, value = turn.b);
  turn.same = lusine.kinematics.unrotate(orientation = identity(), value = turn.a);
  shown.x = turn.back.x;
  shown.y = turn.back.y;
  shown.z = turn.back.z;
  shown.transpose_x = turn.by_transpose.x;
  shown.transpose_y = turn.by_transpose.y;
  shown.transpose_z = turn.by_transpose.z;
  shown.transpose_miss = roundoffs(got = turn.back, want = turn.by_transpose, value = turn.a);
  shown.round_trip_miss = roundoffs(got = turn.there_and_back, want = turn.a, value = turn.a);
  shown.adjoint_left = dot(turn.a_there, turn.b);
  shown.adjoint_right = dot(turn.a, turn.b_back);
  checks.transpose = checks.transpose && shown.transpose_miss <= shown.transpose_bound;
  checks.round_trip = checks.round_trip && shown.round_trip_miss <= shown.round_trip_bound;
  checks.adjoint = checks.adjoint
                   && abs(shown.adjoint_left - shown.adjoint_right)
                      <= 64.0 / 16777216.0 * length(turn.a) * length(turn.b);
  checks.identity = checks.identity && same(a = turn.same, b = turn.a);
  lamps.colour = lamp_colour(holds = checks.transpose && checks.round_trip && checks.adjoint && checks.identity);
}

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