Draw symbols on the GPU from the point table as it is

OpenGL ES 3.1: the file's points go up as one storage buffer, each symbol
gets 16 row bands listing the curves they reach (cut by integer division of
doubled coordinates, so CPU and shader agree), and every pixel of a symbol's
box sums the exact area of its bands' curves, as vega_coverage does.
vega-symbols gpu draws every symbol against the CPU reference and times
1080p pages.

On the Pi 5's V3D it matches the reference (mean 0.001-0.1 levels, worst
4) but takes 142-757 ms a page, 83% of it in the area function. A sorting
network in registers made it 2.3x faster yet computes wrong rows in the
fragment stage with this driver, while the same code is exact in a compute
shader; so drawing each glyph size once into an atlas by compute comes next.

3 files changed, 420 insertions, 2 deletions

diff --git a/src/symbols-gl.c b/src/symbols-gl.c
new file mode 100644
index 0000000..0e78efe
--- /dev/null
+++ b/src/symbols-gl.c
@@ -0,0 +1,260 @@
+#include "symbols-gl.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+enum { BANDS = 16 };
+
+/* Coordinates travel doubled, so implied points stay whole: 2x of a point
+ * is its raw i16 without the flag, and an implied point's is the sum of its
+ * two neighbours' coordinates. */
+#define SHARED \
+    "#version 310 es\n" \
+    "precision highp float;\n" \
+    "precision highp int;\n"
+
+static const char vertex_source[] = SHARED
+    "layout(location = 0) in vec4 place;   /* origin x, y in pixels; pixels per unit; unused */\n"
+    "layout(location = 1) in vec4 bounds;  /* left, bottom, right, top in units */\n"
+    "layout(location = 2) in uint bands;   /* its first band */\n"
+    "uniform vec2 target;\n"
+    "out vec2 units;\n"
+    "flat out float scale;\n"
+    "flat out ivec2 rows;                 /* bottom and top in units */\n"
+    "flat out uint base;\n"
+    "void main() {\n"
+    "    vec2 corner = vec2(float(gl_VertexID & 1), float(gl_VertexID >> 1));\n"
+    "    vec2 pad = vec2(1.0 / place.z);\n"
+    "    units = mix(bounds.xy - pad, bounds.zw + pad, corner);\n"
+    "    vec2 pixel = place.xy + vec2(units.x, -units.y) * place.z;\n"
+    "    gl_Position = vec4(pixel.x / target.x * 2.0 - 1.0, 1.0 - pixel.y / target.y * 2.0, 0.0, 1.0);\n"
+    "    scale = place.z;\n"
+    "    rows = ivec2(bounds.yw);\n"
+    "    base = bands;\n"
+    "}\n";
+
+static const char fragment_source[] = SHARED
+    "layout(std430, binding = 0) readonly buffer Points { uint points[]; };\n"
+    "layout(std430, binding = 1) readonly buffer Bands { uint bands[]; };\n"
+    "uniform uint listed; /* where the lists of curves start in bands */\n"
+    "in vec2 units;\n"
+    "flat in float scale;\n"
+    "flat in ivec2 rows;\n"
+    "flat in uint base;\n"
+    "out vec4 color;\n"
+    "\n"
+    "ivec3 point(uint i) { /* doubled x, doubled y, off */\n"
+    "    uint w = points[i];\n"
+    "    int x = int(w << 16u) >> 16, y = int(w) >> 16;\n"
+    "    return ivec3(x & ~1, y & ~1, x & 1);\n"
+    "}\n"
+    "int band(int doubled) { /* which of 16 bands a doubled y lies in */\n"
+    "    return clamp((doubled - 2 * rows.x) * 16 / (2 * (rows.y - rows.x) + 2), 0, 15);\n"
+    "}\n"
+    "int roots(float a, float b, float c, inout float ts[10], int n) {\n"
+    "    float r0 = -1.0, r1 = -1.0;\n"
+    "    if (a == 0.0) { if (b != 0.0) r0 = -c / b; }\n"
+    "    else {\n"
+    "        float d = b * b - 4.0 * a * c;\n"
+    "        if (d >= 0.0) {\n"
+    "            float q = -0.5 * (b + (b >= 0.0 ? sqrt(d) : -sqrt(d)));\n"
+    "            r0 = q / a;\n"
+    "            if (q != 0.0) r1 = c / q;\n"
+    "        }\n"
+    "    }\n"
+    "    if (r0 > 0.0 && r0 < 1.0) ts[n++] = r0;\n"
+    "    if (r1 > 0.0 && r1 < 1.0) ts[n++] = r1;\n"
+    "    return n;\n"
+    "}\n"
+    "/* As curve_area in symbols.c: the curve's share of the pixel [0, 1]^2. */\n"
+    "float area(vec2 p0, vec2 p1, vec2 p2) {\n"
+    "    if (max(max(p0.y, p1.y), p2.y) <= 0.0 || min(min(p0.y, p1.y), p2.y) >= 1.0 || min(min(p0.x, p1.x), p2.x) >= 1.0)\n"
+    "        return 0.0;\n"
+    "    if (max(max(p0.x, p1.x), p2.x) <= 0.0) return clamp(p2.y, 0.0, 1.0) - clamp(p0.y, 0.0, 1.0);\n"
+    "    vec2 a = p0 - 2.0 * p1 + p2, b = p0 - p1;\n"
+    "    float ts[10];\n"
+    "    ts[0] = 0.0;\n"
+    "    int n = 1;\n"
+    "    n = roots(a.y, -2.0 * b.y, p0.y, ts, n);\n"
+    "    n = roots(a.y, -2.0 * b.y, p0.y - 1.0, ts, n);\n"
+    "    n = roots(a.x, -2.0 * b.x, p0.x, ts, n);\n"
+    "    n = roots(a.x, -2.0 * b.x, p0.x - 1.0, ts, n);\n"
+    "    ts[n++] = 1.0;\n"
+    "    for (int i = 1; i < n; ++i)\n"
+    "        for (int j = i; j > 0 && ts[j] < ts[j - 1]; --j) { float t = ts[j]; ts[j] = ts[j - 1]; ts[j - 1] = t; }\n"
+    "    float u0 = 1.0 - p0.x, u1 = 2.0 * b.x, u2 = -a.x, v0 = -2.0 * b.y, v1 = 2.0 * a.y;\n"
+    "    float k1 = u0 * v0, k2 = (u0 * v1 + u1 * v0) / 2.0, k3 = (u1 * v1 + u2 * v0) / 3.0, k4 = u2 * v1 / 4.0, sum = 0.0;\n"
+    "    for (int i = 0; i + 1 < n; ++i) {\n"
+    "        float ta = ts[i], tb = ts[i + 1], t = (ta + tb) / 2.0;\n"
+    "        if (tb <= ta) continue;\n"
+    "        float y = (a.y * t - 2.0 * b.y) * t + p0.y, x = (a.x * t - 2.0 * b.x) * t + p0.x;\n"
+    "        if (y <= 0.0 || y >= 1.0 || x >= 1.0) continue;\n"
+    "        if (x <= 0.0) sum += (a.y * tb - 2.0 * b.y) * tb - (a.y * ta - 2.0 * b.y) * ta;\n"
+    "        else sum += (((k4 * tb + k3) * tb + k2) * tb + k1) * tb - (((k4 * ta + k3) * ta + k2) * ta + k1) * ta;\n"
+    "    }\n"
+    "    return sum;\n"
+    "}\n"
+    "void main() {\n"
+    "    vec2 corner = units - vec2(0.5 / scale); /* the pixel's lower left, in units */\n"
+    "    int b0 = band(2 * int(floor(corner.y))), b1 = band(2 * int(ceil(corner.y + 1.0 / scale)));\n"
+    "    float sum = 0.0;\n"
+    "    for (int b = b0; b <= b1; ++b) {\n"
+    "        uint at = bands[base + 2u * uint(b)], end = at + bands[base + 2u * uint(b) + 1u];\n"
+    "        for (; at < end; ++at) {\n"
+    "            uint k = bands[listed + at];\n"
+    "            ivec3 p = point(k), n = point(k + 1u);\n"
+    "            ivec2 s, c, e;\n"
+    "            if (p.z == 1) {\n"
+    "                ivec3 q = point(k - 1u);\n"
+    "                s = q.z == 1 ? (q.xy + p.xy) / 2 : q.xy;\n"
+    "                c = p.xy;\n"
+    "                e = n.z == 1 ? (p.xy + n.xy) / 2 : n.xy;\n"
+    "            } else {\n"
+    "                s = p.xy;\n"
+    "                c = (p.xy + n.xy) / 2;\n"
+    "                e = n.xy;\n"
+    "            }\n"
+    "            /* a curve over several of this pixel's bands counts in the first */\n"
+    "            if (max(band(min(min(s.y, c.y), e.y)), b0) != b) continue;\n"
+    "            sum += area((vec2(s) * 0.5 - corner) * scale, (vec2(c) * 0.5 - corner) * scale, (vec2(e) * 0.5 - corner) * scale);\n"
+    "        }\n"
+    "    }\n"
+    "    color = vec4(min(abs(sum), 1.0), 0.0, 0.0, 1.0);\n"
+    "}\n";
+
+static GLuint compile(GLenum type, const char *source, char *error, size_t error_size) {
+    GLuint shader = glCreateShader(type);
+    glShaderSource(shader, 1, &source, NULL);
+    glCompileShader(shader);
+    GLint ok;
+    glGetShaderiv(shader, GL_COMPILE_STATUS, &ok);
+    if (!ok) {
+        glGetShaderInfoLog(shader, (GLsizei)error_size, NULL, error);
+        glDeleteShader(shader);
+        return 0;
+    }
+    return shader;
+}
+
+/* Doubled coordinates make implied points whole; bands as the shader cuts them. */
+static int band(int doubled, int bottom, int top) {
+    int b = (doubled - 2 * bottom) * BANDS / (2 * (top - bottom) + 2);
+    return b < 0 ? 0 : b >= BANDS ? BANDS - 1 : b;
+}
+
+bool vega_gl_open(VegaGl *gl, const VegaSymbols *s, char *error, size_t error_size) {
+    *gl = (VegaGl){0};
+    GLuint vertex = compile(GL_VERTEX_SHADER, vertex_source, error, error_size);
+    GLuint fragment = vertex ? compile(GL_FRAGMENT_SHADER, fragment_source, error, error_size) : 0;
+    if (!fragment) return false;
+    gl->program = glCreateProgram();
+    glAttachShader(gl->program, vertex);
+    glAttachShader(gl->program, fragment);
+    glLinkProgram(gl->program);
+    glDeleteShader(vertex);
+    glDeleteShader(fragment);
+    GLint ok;
+    glGetProgramiv(gl->program, GL_LINK_STATUS, &ok);
+    if (!ok) {
+        glGetProgramInfoLog(gl->program, (GLsizei)error_size, NULL, error);
+        return false;
+    }
+    /* The point table, as the file has it. */
+    glGenBuffers(1, &gl->points);
+    glBindBuffer(GL_SHADER_STORAGE_BUFFER, gl->points);
+    glBufferData(GL_SHADER_STORAGE_BUFFER, (GLsizeiptr)(s->points ? s->points : 1) * VEGA_POINT_BYTES,
+                 s->bytes + VEGA_SYMBOLS_HEADER + (size_t)s->symbols * VEGA_SYMBOL_BYTES, GL_STATIC_DRAW);
+    /* Bands: per symbol, 16 (start, count) pairs into one list of curves,
+     * each curve (by the point that begins it) in every band it reaches. */
+    size_t table = (size_t)s->symbols * BANDS * 2, listed = 0;
+    uint32_t *counts = calloc(table + 1, sizeof(*counts));
+    for (int pass = 0; pass < 2; ++pass) {
+        uint32_t *bands = pass ? malloc((table + listed + 1) * sizeof(*bands)) : NULL;
+        size_t next = 0;
+        for (uint32_t i = 0; i < s->symbols; ++i) {
+            VegaSymbol symbol = vega_symbol_at(s, i);
+            uint32_t *mine = pass ? bands + (size_t)i * BANDS * 2 : counts + (size_t)i * BANDS * 2;
+            if (pass)
+                for (int b = 0; b < BANDS; ++b) {
+                    mine[2 * b] = (uint32_t)next;
+                    next += counts[(size_t)i * BANDS * 2 + 2 * b + 1];
+                    mine[2 * b + 1] = 0;
+                }
+            VegaCurve c;
+            for (uint32_t k = symbol.first; k < symbol.first + symbol.count; ++k) {
+                if (!vega_symbol_curve(s, &symbol, k, &c)) continue;
+                float low = c.y0 < c.y1 ? c.y0 : c.y1, high = c.y0 > c.y1 ? c.y0 : c.y1;
+                low = c.y2 < low ? c.y2 : low;
+                high = c.y2 > high ? c.y2 : high;
+                int from = band((int)(2 * low), symbol.bottom, symbol.top), to = band((int)(2 * high), symbol.bottom, symbol.top);
+                for (int b = from; b <= to; ++b) {
+                    if (pass) bands[table + mine[2 * b] + mine[2 * b + 1]] = k;
+                    ++mine[2 * b + 1];
+                    listed += !pass;
+                }
+            }
+        }
+        if (pass) {
+            glGenBuffers(1, &gl->bands);
+            glBindBuffer(GL_SHADER_STORAGE_BUFFER, gl->bands);
+            glBufferData(GL_SHADER_STORAGE_BUFFER, (GLsizeiptr)((table + listed + 1) * sizeof(*bands)), bands, GL_STATIC_DRAW);
+            free(bands);
+        }
+    }
+    free(counts);
+    glGenVertexArrays(1, &gl->array);
+    glGenBuffers(1, &gl->instances);
+    glBindVertexArray(gl->array);
+    glBindBuffer(GL_ARRAY_BUFFER, gl->instances);
+    glEnableVertexAttribArray(0);
+    glEnableVertexAttribArray(1);
+    glEnableVertexAttribArray(2);
+    glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 36, (void *)0);
+    glVertexAttribPointer(1, 4, GL_FLOAT, GL_FALSE, 36, (void *)16);
+    glVertexAttribIPointer(2, 1, GL_UNSIGNED_INT, 36, (void *)32);
+    glVertexAttribDivisor(0, 1);
+    glVertexAttribDivisor(1, 1);
+    glVertexAttribDivisor(2, 1);
+    glBindVertexArray(0);
+    glUseProgram(gl->program);
+    glUniform1ui(glGetUniformLocation(gl->program, "listed"), (GLuint)table);
+    gl->viewport = glGetUniformLocation(gl->program, "target");
+    if (glGetError() != GL_NO_ERROR) {
+        snprintf(error, error_size, "GL error while uploading");
+        return false;
+    }
+    return true;
+}
+
+void vega_gl_draw(VegaGl *gl, const VegaSymbols *s, const VegaGlyph *glyphs, uint32_t count, int width, int height) {
+    float *instances = malloc((size_t)(count ? count : 1) * 36);
+    for (uint32_t i = 0; i < count; ++i) {
+        VegaSymbol symbol = vega_symbol_at(s, glyphs[i].symbol);
+        float *f = instances + (size_t)i * 9;
+        f[0] = glyphs[i].x, f[1] = glyphs[i].y, f[2] = glyphs[i].scale, f[3] = 0;
+        f[4] = symbol.left, f[5] = symbol.bottom, f[6] = symbol.right, f[7] = symbol.top;
+        uint32_t base = glyphs[i].symbol * BANDS * 2;
+        memcpy(f + 8, &base, 4);
+    }
+    glUseProgram(gl->program);
+    glUniform2f(gl->viewport, (float)width, (float)height);
+    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, gl->points);
+    glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, gl->bands);
+    glBindVertexArray(gl->array);
+    glBindBuffer(GL_ARRAY_BUFFER, gl->instances);
+    glBufferData(GL_ARRAY_BUFFER, (GLsizeiptr)count * 36, instances, GL_STREAM_DRAW);
+    glEnable(GL_BLEND);
+    glBlendFunc(GL_ONE, GL_ONE);
+    glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, (GLsizei)count);
+    glBindVertexArray(0);
+    free(instances);
+}
+
+void vega_gl_close(VegaGl *gl) {
+    glDeleteBuffers(1, &gl->points);
+    glDeleteBuffers(1, &gl->bands);
+    glDeleteBuffers(1, &gl->instances);
+    glDeleteVertexArrays(1, &gl->array);
+    glDeleteProgram(gl->program);
+    *gl = (VegaGl){0};
+}
diff --git a/src/symbols-gl.h b/src/symbols-gl.h
new file mode 100644
index 0000000..4c6db67
--- /dev/null
+++ b/src/symbols-gl.h
@@ -0,0 +1,24 @@
+#ifndef VEGA_SYMBOLS_GL_H
+#define VEGA_SYMBOLS_GL_H
+#include "symbols.h"
+#include <GLES3/gl31.h>
+
+/* Symbols drawn by the GPU (OpenGL ES 3.1): the file's point table goes up
+ * as it is, one storage buffer, and each pixel of a symbol's box sums the
+ * exact area of the curves crossing its row, as vega_coverage does. Rows
+ * are cut into bands at upload, each listing its curves, so a pixel reads
+ * only those. Needs a current context. */
+typedef struct {
+    GLuint program, points, bands, instances, array;
+    GLint viewport;
+} VegaGl;
+
+/* A symbol placed: its origin in pixels from the target's top left, and
+ * pixels per unit. */
+typedef struct { float x, y, scale; uint32_t symbol; } VegaGlyph;
+
+bool vega_gl_open(VegaGl *, const VegaSymbols *, char *error, size_t error_size);
+/* Coverage added into the red channel of the bound framebuffer. */
+void vega_gl_draw(VegaGl *, const VegaSymbols *, const VegaGlyph *, uint32_t count, int width, int height);
+void vega_gl_close(VegaGl *);
+#endif
diff --git a/src/symbols-tool.c b/src/symbols-tool.c
index e957c9e..a94feb3 100644
--- a/src/symbols-tool.c
+++ b/src/symbols-tool.c
@@ -4,8 +4,12 @@
  *   vega-symbols build [--font FONT.ttf --codes 0x20-0x7e,...] [--svg CODE=ICON.svg ...] [--units N] -o OUT.vsym
  *   vega-symbols compare FONT.ttf OUT.vsym --sizes 12,16,32,64
  *   vega-symbols info OUT.vsym
- *   vega-symbols render OUT.vsym CODE PIXELS OUT.pgm */
+ *   vega-symbols render OUT.vsym CODE PIXELS OUT.pgm
+ *   vega-symbols gpu OUT.vsym --sizes 12,16,32,64 [--frames N] */
+#include "symbols-gl.h"
 #include "symbols.h"
+#include <EGL/egl.h>
+#include <EGL/eglext.h>
 #include <ft2build.h>
 #include FT_FREETYPE_H
 #include FT_OUTLINE_H
@@ -737,11 +741,141 @@ static int render(int argc, char **argv) {
     return 0;
 }
 
+/* ---- The GPU: a context of its own, no window (EGL, Mesa's surfaceless platform). */
+
+static void gl_context(void) {
+    PFNEGLGETPLATFORMDISPLAYEXTPROC platform = (PFNEGLGETPLATFORMDISPLAYEXTPROC)eglGetProcAddress("eglGetPlatformDisplayEXT");
+    EGLDisplay display = platform ? platform(EGL_PLATFORM_SURFACELESS_MESA, EGL_DEFAULT_DISPLAY, NULL) : EGL_NO_DISPLAY;
+    EGLint major, minor, found;
+    EGLConfig config;
+    EGLint wanted[] = {EGL_RENDERABLE_TYPE, EGL_OPENGL_ES3_BIT, EGL_NONE};
+    EGLint version[] = {EGL_CONTEXT_MAJOR_VERSION, 3, EGL_CONTEXT_MINOR_VERSION, 1, EGL_NONE};
+    if (display == EGL_NO_DISPLAY || !eglInitialize(display, &major, &minor) || !eglBindAPI(EGL_OPENGL_ES_API) ||
+        !eglChooseConfig(display, wanted, &config, 1, &found))
+        fail("no EGL display for OpenGL ES 3.1");
+    /* With no surface, a context needs no config (EGL_KHR_no_config_context). */
+    EGLContext context = eglCreateContext(display, found ? config : EGL_NO_CONFIG_KHR, EGL_NO_CONTEXT, version);
+    if (context == EGL_NO_CONTEXT || !eglMakeCurrent(display, EGL_NO_SURFACE, EGL_NO_SURFACE, context))
+        fail("no OpenGL ES 3.1 context");
+}
+
+/* A width x height target of 8-bit red, bound. */
+static GLuint gl_target(int width, int height) {
+    GLuint texture, framebuffer;
+    glGenTextures(1, &texture);
+    glBindTexture(GL_TEXTURE_2D, texture);
+    glTexStorage2D(GL_TEXTURE_2D, 1, GL_R8, width, height);
+    glGenFramebuffers(1, &framebuffer);
+    glBindFramebuffer(GL_FRAMEBUFFER, framebuffer);
+    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, texture, 0);
+    if (glCheckFramebufferStatus(GL_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE) fail("no framebuffer");
+    glViewport(0, 0, width, height);
+    glClearColor(0, 0, 0, 0);
+    glClear(GL_COLOR_BUFFER_BIT);
+    return texture;
+}
+
+/* Every symbol in a grid of cells, drawn by the GPU and by vega_coverage:
+ * how far apart; then pages of text (glyphs cycling through the file)
+ * drawn again and again: how long a page takes. */
+static int gpu(int argc, char **argv) {
+    if (argc < 5 || strcmp(argv[3], "--sizes")) fail("usage: vega-symbols gpu SYMBOLS.vsym --sizes 12,16,32,64 [--frames N]");
+    int frames = argc > 6 && !strcmp(argv[5], "--frames") ? atoi(argv[6]) : 20;
+    if (frames < 1) fail("--frames");
+    char *bytes, error[4096];
+    VegaSymbols s = open_symbols(argv[2], &bytes);
+    gl_context();
+    printf("{\"renderer\":\"%s\"}\n", (const char *)glGetString(GL_RENDERER));
+    VegaGl gl;
+    double started = seconds();
+    if (!vega_gl_open(&gl, &s, error, sizeof(error))) fail(error);
+    glFinish();
+    double open_ms = (seconds() - started) * 1e3;
+    for (const char *p = argv[4]; *p;) {
+        char *end;
+        long size = strtol(p, &end, 10);
+        p = *end == ',' ? end + 1 : end;
+        if (size <= 0 || size > 512) fail("--sizes");
+        double scale = (double)size / s.units;
+        /* Cells that hold any symbol, so none reaches into another's. */
+        int left = 0, bottom = 0, right = 0, top = 0;
+        for (uint32_t i = 0; i < s.symbols; ++i) {
+            VegaSymbol symbol = vega_symbol_at(&s, i);
+            left = symbol.left < left ? symbol.left : left;
+            bottom = symbol.bottom < bottom ? symbol.bottom : bottom;
+            right = symbol.right > right ? symbol.right : right;
+            top = symbol.top > top ? symbol.top : top;
+        }
+        int cell = (int)ceil((right - left > top - bottom ? right - left : top - bottom) * scale) + 4;
+        int across = 4096 / cell, rows = (int)((s.symbols + (uint32_t)across - 1) / (uint32_t)across);
+        if (rows * cell > 4096) rows = 4096 / cell;
+        int width = across * cell, height = rows * cell;
+        uint32_t shown = s.symbols < (uint32_t)(across * rows) ? s.symbols : (uint32_t)(across * rows);
+        VegaGlyph *glyphs = malloc(shown * sizeof(*glyphs));
+        for (uint32_t i = 0; i < shown; ++i)
+            glyphs[i] = (VegaGlyph){(float)(i % (uint32_t)across * (uint32_t)cell) + 2 - (float)(left * scale),
+                                    (float)(i / (uint32_t)across * (uint32_t)cell) + 2 + (float)(top * scale), (float)scale, i};
+        GLuint texture = gl_target(width, height);
+        vega_gl_draw(&gl, &s, glyphs, shown, width, height);
+        unsigned char *drawn = malloc((size_t)width * (size_t)height), *reference = malloc((size_t)cell * (size_t)cell);
+        glPixelStorei(GL_PACK_ALIGNMENT, 1);
+        glReadPixels(0, 0, width, height, GL_RED, GL_UNSIGNED_BYTE, drawn);
+        double total = 0;
+        uint64_t pixels = 0;
+        int worst = 0;
+        for (uint32_t i = 0; i < shown; ++i) {
+            VegaSymbol symbol = vega_symbol_at(&s, i);
+            int x0 = (int)(i % (uint32_t)across) * cell, y0 = (int)(i / (uint32_t)across) * cell;
+            vega_symbol_raster(&s, &symbol, scale, (x0 - glyphs[i].x) / scale, (glyphs[i].y - y0) / scale, cell, cell, reference);
+            for (int j = 0; j < cell; ++j)
+                for (int k = 0; k < cell; ++k) {
+                    /* GL's rows run bottom up */
+                    int gpu_value = drawn[(size_t)(height - 1 - (y0 + j)) * (size_t)width + (size_t)(x0 + k)];
+                    int difference = abs(gpu_value - reference[j * cell + k]);
+                    total += difference;
+                    worst = difference > worst ? difference : worst;
+                    ++pixels;
+                }
+        }
+        /* Pages: glyphs at their advances in lines, a full target's worth. */
+        int page_width = 1920, page_height = 1080;
+        uint32_t page = 0, capacity = 1 << 16;
+        VegaGlyph *text = malloc(capacity * sizeof(*text));
+        float x = 0, line = (float)size * 1.2f;
+        for (float y = line; y + (float)size * 0.3f < (float)page_height && page < capacity; ) {
+            VegaSymbol symbol = vega_symbol_at(&s, page % s.symbols);
+            float advance = (float)(symbol.advance * scale);
+            if (x + advance > (float)page_width) { x = 0; y += line; continue; }
+            text[page++] = (VegaGlyph){x, y, (float)scale, page % s.symbols};
+            x += advance;
+        }
+        glDeleteTextures(1, &texture);
+        texture = gl_target(page_width, page_height);
+        vega_gl_draw(&gl, &s, text, page, page_width, page_height);
+        glFinish();
+        started = seconds();
+        for (int f = 0; f < frames; ++f) vega_gl_draw(&gl, &s, text, page, page_width, page_height);
+        glFinish();
+        double page_ms = (seconds() - started) / frames * 1e3;
+        printf("{\"pixels_per_em\":%ld,\"symbols\":%u,\"mean_difference\":%.3f,\"worst\":%d,\"page_glyphs\":%u,"
+               "\"page_ms\":%.2f,\"upload_ms\":%.1f}\n", size, shown, pixels ? total / (double)pixels : 0, worst, page, page_ms, open_ms);
+        glDeleteTextures(1, &texture);
+        free(text);
+        free(glyphs);
+        free(drawn);
+        free(reference);
+    }
+    vega_gl_close(&gl);
+    free(bytes);
+    return 0;
+}
+
 int main(int argc, char **argv) {
     if (argc > 1 && !strcmp(argv[1], "build")) return build(argc, argv);
     if (argc > 1 && !strcmp(argv[1], "compare")) return compare(argc, argv);
     if (argc > 1 && !strcmp(argv[1], "info")) return info(argc, argv);
     if (argc > 1 && !strcmp(argv[1], "render")) return render(argc, argv);
-    fail("vega-symbols build|compare|info|render ...");
+    if (argc > 1 && !strcmp(argv[1], "gpu")) return gpu(argc, argv);
+    fail("vega-symbols build|compare|info|render|gpu ...");
     return 2;
 }