Intel Larrabee / KNF Vibe Coded 3D accelerated Driver?

[erek]

since the community has been unable to source the 3d drivers for the Larrabee

DirectXGfx "under FreeBSD it's just running an x86 program called DirectXGfx (248 threads of it). And it shares a file system with the host and you can telnet into it and give it other work to do and steal cores from your own graphics system - it was mind-bending!"
(reference: https://tomforsyth1000.github.io/blog.wiki.html#[[Why didn't Larrabee fail?]])

View: https://www.youtube.com/watch?v=0qw-z2SBpE4


so.. i think Agentic AI is powerful enough to just build us the software to run games on the a KNF / Xeon Phi co-processor, or even Real prototype Larrabee.

---

using chatgpt, glm-4.7 and Kimi K3:
Check out this conversation I had
https://chat.kimik3.pro/f/NzUzQXcJTICv?c=upfHyZ

Here is the complete, updated LRBni rasterizer driver with all corrections applied — correct 28.4 edge math, integer gather, safe framebuffer writes, and proper perspective interpolation.

```cpp
/*
* LRBni (Knights Ferry) Rasterizer — Final Corrected Version
*
* All vertex coordinates in 28.4 fixed-point (<< 4).
* Edge derivatives are per-pixel steps in 28.4 units.
* Pixel centers at (x*16 + 8, y*16 + 8).
* Edge values in 28.4^2 units (sign only matters for coverage).
*
* Intrinsic names are best guesses for ICC 12.x LRBni.
* Verify against actual zmmintrin.h before production use.
*/

#include <immintrin.h>
#include <stdint.h>
#include <math.h>

/* -----------------------------------------------------------------------
* Type definitions (unified register file, but compiler needs distinct types)
* ----------------------------------------------------------------------- */
typedef __m512 v16f; /* 16 floats */
typedef __m512 v16i; /* 16 integers in same register file */

/* -----------------------------------------------------------------------
* §1 Intrinsic wrappers — verify against zmmintrin.h
* ----------------------------------------------------------------------- */

/* Integer gather: load 32-bit values from base[index] with mask */
static inline v16i gather_epi32(const void* base, v16i indices, unsigned short mask)
{
/* LRBni: _mm512_mask_i32gather_epi32(src, base, index, mask, scale) */
v16i src = _mm512_setzero_epi32();
return _mm512_mask_i32gather_epi32(src, base, indices, mask, 4);
}

/* Extract 16-bit mask from vector of comparison results */
static inline unsigned short get_mask(v16i a)
{
/* If _mm512_movemask_epi32 exists, use it. Otherwise cast to float. */
return (unsigned short)_mm512_movemask_epi32(a);
}

/* Store 16 integers to memory */
static inline void store_epi32(void* addr, v16i a)
{
/* If _mm512_storeu_epi32 exists, use it. Otherwise cast to float store. */
_mm512_storeu_epi32(addr, a);
}

/* -----------------------------------------------------------------------
* §2 Triangle Setup
* ----------------------------------------------------------------------- */

typedef struct {
float x, y, z, w, u, v;
} Vertex;

typedef struct {
/* Edge data (28.4) */
int32_t e01_dx, e01_dy; /* per-pixel steps in 28.4 */
int32_t e12_dx, e12_dy;
int32_t e20_dx, e20_dy;
int32_t e01_row0, e12_row0, e20_row0; /* initial edge values (28.4^2) */

/* Bounding box (pixels) */
int xmin, ymin, xmax, ymax;

/* Perspective-correct UV (float) */
float diw_dx, diw_dy, duw_dx, duw_dy, dvw_dx, dvw_dy;
float iw_row0, uw_row0, vw_row0;

/* Texture */
int tex_w, tex_h, tex_w_mask, tex_h_mask, is_pow2;

/* Framebuffer */
int fb_stride;
uint32_t* fb_base;
} TriSetup;

static int triSetup(const Vertex* v0, const Vertex* v1, const Vertex* v2,
int fb_w, int fb_h, int tex_w, int tex_h,
uint32_t* fb, TriSetup* ts)
{
/* Convert vertices to 28.4 fixed-point */
int32_t v0x = (int32_t)(v0->x * 16.0f);
int32_t v0y = (int32_t)(v0->y * 16.0f);
int32_t v1x = (int32_t)(v1->x * 16.0f);
int32_t v1y = (int32_t)(v1->y * 16.0f);
int32_t v2x = (int32_t)(v2->x * 16.0f);
int32_t v2y = (int32_t)(v2->y * 16.0f);

/* Edge derivatives (per pixel step, in 28.4)
* e_dx = (y0 - y1) * 16 (since v0y = y0*16, v1y = y1*16)
* e_dy = (x1 - x0) * 16
*/
ts->e01_dx = (v0y - v1y) * 16;
ts->e01_dy = (v1x - v0x) * 16;
ts->e12_dx = (v1y - v2y) * 16;
ts->e12_dy = (v2x - v1x) * 16;
ts->e20_dx = (v2y - v0y) * 16;
ts->e20_dy = (v0x - v2x) * 16;

/* Area (64-bit to avoid overflow) — for backface culling */
int64_t area = (int64_t)(v0x - v2x) * (v1y - v2y)
- (int64_t)(v0y - v2y) * (v1x - v2x);
if (area <= 0) return 0; /* backface or degenerate */

/* Bounding box in pixel coordinates */
ts->xmin = (int)fminf(fminf(v0->x, v1->x), v2->x);
if (ts->xmin < 0) ts->xmin = 0;
ts->xmax = (int)fmaxf(fmaxf(v0->x, v1->x), v2->x);
if (ts->xmax >= fb_w) ts->xmax = fb_w - 1;
ts->ymin = (int)fminf(fminf(v0->y, v1->y), v2->y);
if (ts->ymin < 0) ts->ymin = 0;
ts->ymax = (int)fmaxf(fmaxf(v0->y, v1->y), v2->y);
if (ts->ymax >= fb_h) ts->ymax = fb_h - 1;
if (ts->xmin > ts->xmax || ts->ymin > ts->ymax) return 0;

/* Initial edge values at (xmin, ymin) pixel center
* Pixel center in 28.4: (x*16 + 8, y*16 + 8)
* Edge: e = (sy - v0y)*(v1x - v0x) - (sx - v0x)*(v1y - v0y)
*/
int32_t sx = ts->xmin * 16 + 8;
int32_t sy = ts->ymin * 16 + 8;

ts->e01_row0 = (sy - v0y) * (v1x - v0x) - (sx - v0x) * (v1y - v0y);
ts->e12_row0 = (sy - v1y) * (v2x - v1x) - (sx - v1x) * (v2y - v1y);
ts->e20_row0 = (sy - v2y) * (v0x - v2x) - (sx - v2x) * (v0y - v2y);

/* Perspective-correct UV interpolation setup (float) */
float inv_area = 1.0f / (float)(area >> 8); /* area / 256 */
float inv_w0 = 1.0f / v0->w;
float inv_w1 = 1.0f / v1->w;
float inv_w2 = 1.0f / v2->w;
float uw0 = v0->u * inv_w0, uw1 = v1->u * inv_w1, uw2 = v2->u * inv_w2;
float vw0 = v0->v * inv_w0, vw1 = v1->v * inv_w1, vw2 = v2->v * inv_w2;

/* Barycentric gradients in pixel units */
float w0_dx = (float)(v1y - v2y) * inv_area;
float w0_dy = (float)(v2x - v1x) * inv_area;
float w1_dx = (float)(v2y - v0y) * inv_area;
float w1_dy = (float)(v0x - v2x) * inv_area;
float w2_dx = (float)(v0y - v1y) * inv_area;
float w2_dy = (float)(v1x - v0x) * inv_area;

/* Barycentric values at (xmin, ymin) pixel center */
float w0_0 = (float)ts->e12_row0 * inv_area / 256.0f;
float w1_0 = (float)ts->e20_row0 * inv_area / 256.0f;
float w2_0 = (float)ts->e01_row0 * inv_area / 256.0f;

/* Perspective-correct attribute gradients */
ts->diw_dx = w0_dx * inv_w0 + w1_dx * inv_w1 + w2_dx * inv_w2;
ts->diw_dy = w0_dy * inv_w0 + w1_dy * inv_w1 + w2_dy * inv_w2;
ts->duw_dx = w0_dx * uw0 + w1_dx * uw1 + w2_dx * uw2;
ts->duw_dy = w0_dy * uw0 + w1_dy * uw1 + w2_dy * uw2;
ts->dvw_dx = w0_dx * vw0 + w1_dx * vw1 + w2_dx * vw2;
ts->dvw_dy = w0_dy * vw0 + w1_dy * vw1 + w2_dy * vw2;

ts->iw_row0 = w0_0 * inv_w0 + w1_0 * inv_w1 + w2_0 * inv_w2;
ts->uw_row0 = w0_0 * uw0 + w1_0 * uw1 + w2_0 * uw2;
ts->vw_row0 = w0_0 * vw0 + w1_0 * vw1 + w2_0 * vw2;

/* Texture info */
ts->tex_w = tex_w;
ts->tex_h = tex_h;
ts->is_pow2 = ((tex_w & (tex_w - 1)) == 0) && ((tex_h & (tex_h - 1)) == 0);
ts->tex_w_mask = tex_w - 1;
ts->tex_h_mask = tex_h - 1;
ts->fb_stride = fb_w;
ts->fb_base = fb;

return 1;
}

/* -----------------------------------------------------------------------
* §3 Pixel Pipeline (16-wide SIMD)
* ----------------------------------------------------------------------- */

/* Unpack RGBA from packed 32-bit integer */
static inline void unpack_rgba(v16i packed, v16i* r, v16i* g, v16i* b, v16i* a)
{
v16i mask_ff = _mm512_set1_epi32(0xFF);
*r = _mm512_and_epi32(packed, mask_ff);
*g = _mm512_and_epi32(_mm512_srli_epi32(packed, 8), mask_ff);
*b = _mm512_and_epi32(_mm512_srli_epi32(packed, 16), mask_ff);
*a = _mm512_srli_epi32(packed, 24);
}

/* Wrap or clamp texture coordinate */
static inline v16i wrap_coord(v16i coord, int size, int mask, int is_pow2)
{
if (is_pow2) {
return _mm512_and_epi32(coord, _mm512_set1_epi32(mask));
} else {
/* Clamp to [0, size-1] for non-power-of-two textures */
v16i zero = _mm512_setzero_epi32();
v16i max = _mm512_set1_epi32(size - 1);
return _mm512_min_epi32(_mm512_max_epi32(coord, zero), max);
}
}

/* Shade 16 pixels (8x2 tile) */
static void shadePixels16(
const TriSetup* ts, const uint32_t* texture,
v16i dx_pix, v16i dy_pix, v16i px_i, v16i py_i, unsigned short mask)
{
/* ---- Perspective interpolation (float) ---- */
v16f dx_f = _mm512_cvtepi32_ps(dx_pix);
v16f dy_f = _mm512_cvtepi32_ps(dy_pix);

v16f iw = _mm512_add_ps(
_mm512_add_ps(_mm512_set1_ps(ts->iw_row0),
_mm512_mul_ps(_mm512_set1_ps(ts->diw_dx), dx_f)),
_mm512_mul_ps(_mm512_set1_ps(ts->diw_dy), dy_f));

v16f uw = _mm512_add_ps(
_mm512_add_ps(_mm512_set1_ps(ts->uw_row0),
_mm512_mul_ps(_mm512_set1_ps(ts->duw_dx), dx_f)),
_mm512_mul_ps(_mm512_set1_ps(ts->duw_dy), dy_f));

v16f vw = _mm512_add_ps(
_mm512_add_ps(_mm512_set1_ps(ts->vw_row0),
_mm512_mul_ps(_mm512_set1_ps(ts->dvw_dx), dx_f)),
_mm512_mul_ps(_mm512_set1_ps(ts->dvw_dy), dy_f));

/* Perspective divide (Newton-Raphson refinement) */
v16f rcp_iw = _mm512_rcp_ps(iw);
rcp_iw = _mm512_mul_ps(rcp_iw,
_mm512_sub_ps(_mm512_set1_ps(2.0f),
_mm512_mul_ps(iw, rcp_iw)));

v16f u = _mm512_mul_ps(uw, rcp_iw);
v16f v = _mm512_mul_ps(vw, rcp_iw);

/* ---- Texel coordinates ---- */
v16f fu = _mm512_mul_ps(u, _mm512_set1_ps((float)ts->tex_w));
v16f fv = _mm512_mul_ps(v, _mm512_set1_ps((float)ts->tex_h));

v16i iu = _mm512_cvttps_epi32(fu);
v16i iv = _mm512_cvttps_epi32(fv);

v16f frac_u = _mm512_sub_ps(fu, _mm512_cvtepi32_ps(iu));
v16f frac_v = _mm512_sub_ps(fv, _mm512_cvtepi32_ps(iv));

/* ---- Bilinear tap coordinates ---- */
v16i iu0 = wrap_coord(iu, ts->tex_w, ts->tex_w_mask, ts->is_pow2);
v16i iv0 = wrap_coord(iv, ts->tex_h, ts->tex_h_mask, ts->is_pow2);
v16i iu1 = wrap_coord(_mm512_add_epi32(iu0, _mm512_set1_epi32(1)),
ts->tex_w, ts->tex_w_mask, ts->is_pow2);
v16i iv1 = wrap_coord(_mm512_add_epi32(iv0, _mm512_set1_epi32(1)),
ts->tex_h, ts->tex_h_mask, ts->is_pow2);

/* ---- Linear addresses: addr = y * width + x ---- */
v16i stride = _mm512_set1_epi32(ts->tex_w);
v16i addr00 = _mm512_add_epi32(_mm512_mullo_epi32(iv0, stride), iu0);
v16i addr10 = _mm512_add_epi32(_mm512_mullo_epi32(iv0, stride), iu1);
v16i addr01 = _mm512_add_epi32(_mm512_mullo_epi32(iv1, stride), iu0);
v16i addr11 = _mm512_add_epi32(_mm512_mullo_epi32(iv1, stride), iu1);

/* ---- Gather 4 texels (raw RGBA as 32-bit integers) ---- */
v16i c00 = gather_epi32(texture, addr00, mask);
v16i c10 = gather_epi32(texture, addr10, mask);
v16i c01 = gather_epi32(texture, addr01, mask);
v16i c11 = gather_epi32(texture, addr11, mask);

/* ---- Unpack RGBA channels ---- */
v16i r00, g00, b00, a00, r10, g10, b10, a10;
v16i r01, g01, b01, a01, r11, g11, b11, a11;
unpack_rgba(c00, &r00, &g00, &b00, &a00);
unpack_rgba(c10, &r10, &g10, &b10, &a10);
unpack_rgba(c01, &r01, &g01, &b01, &a01);
unpack_rgba(c11, &r11, &g11, &b11, &a11);

/* ---- Bilinear weights ---- */
v16f w0 = _mm512_mul_ps(
_mm512_sub_ps(_mm512_set1_ps(1.0f), frac_u),
_mm512_sub_ps(_mm512_set1_ps(1.0f), frac_v));
v16f w1 = _mm512_mul_ps(
frac_u,
_mm512_sub_ps(_mm512_set1_ps(1.0f), frac_v));
v16f w2 = _mm512_mul_ps(
_mm512_sub_ps(_mm512_set1_ps(1.0f), frac_u),
frac_v);
v16f w3 = _mm512_mul_ps(frac_u, frac_v);

/* ---- Blend channels ---- */
v16f r = _mm512_add_ps(
_mm512_add_ps(_mm512_mul_ps(w0, _mm512_cvtepi32_ps(r00)),
_mm512_mul_ps(w1, _mm512_cvtepi32_ps(r10))),
_mm512_add_ps(_mm512_mul_ps(w2, _mm512_cvtepi32_ps(r01)),
_mm512_mul_ps(w3, _mm512_cvtepi32_ps(r11))));
v16f g = _mm512_add_ps(
_mm512_add_ps(_mm512_mul_ps(w0, _mm512_cvtepi32_ps(g00)),
_mm512_mul_ps(w1, _mm512_cvtepi32_ps(g10))),
_mm512_add_ps(_mm512_mul_ps(w2, _mm512_cvtepi32_ps(g01)),
_mm512_mul_ps(w3, _mm512_cvtepi32_ps(g11))));
v16f b = _mm512_add_ps(
_mm512_add_ps(_mm512_mul_ps(w0, _mm512_cvtepi32_ps(b00)),
_mm512_mul_ps(w1, _mm512_cvtepi32_ps(b10))),
_mm512_add_ps(_mm512_mul_ps(w2, _mm512_cvtepi32_ps(b01)),
_mm512_mul_ps(w3, _mm512_cvtepi32_ps(b11))));
v16f a = _mm512_add_ps(
_mm512_add_ps(_mm512_mul_ps(w0, _mm512_cvtepi32_ps(a00)),
_mm512_mul_ps(w1, _mm512_cvtepi32_ps(a10))),
_mm512_add_ps(_mm512_mul_ps(w2, _mm512_cvtepi32_ps(a01)),
_mm512_mul_ps(w3, _mm512_cvtepi32_ps(a11))));

/* ---- Pack to uint32_t ---- */
v16i ir = _mm512_cvtps_epi32(r);
v16i ig = _mm512_slli_epi32(_mm512_cvtps_epi32(g), 8);
v16i ib = _mm512_slli_epi32(_mm512_cvtps_epi32(b), 16);
v16i ia = _mm512_slli_epi32(_mm512_cvtps_epi32(a), 24);
v16i packed = _mm512_or_epi32(
_mm512_or_epi32(ir, ig),
_mm512_or_epi32(ib, ia));

/* ---- Framebuffer write (scalar fallback for safety) ---- */
int32_t px[16], py[16], col[16];
store_epi32(px, px_i);
store_epi32(py, py_i);
store_epi32(col, packed);

unsigned short mask_row0 = mask & 0x00FF;
unsigned short mask_row1 = (mask >> 8) & 0x00FF;

/* Row 0 (lanes 0-7) */
if (mask_row0) {
uint32_t* fb0 = ts->fb_base + py[0] * ts->fb_stride + px[0];
for (int i = 0; i < 8; i++) {
if (mask_row0 & (1 << i)) {
fb0 = (uint32_t)col;
}
}
}

/* Row 1 (lanes 8-15) */
if (mask_row1) {
uint32_t* fb1 = ts->fb_base + py[8] * ts->fb_stride + px[8];
for (int i = 0; i < 8; i++) {
if (mask_row1 & (1 << i)) {
fb1 = (uint32_t)col[i + 8];
}
}
}
}

/* -----------------------------------------------------------------------
* §4 Tile Driver (8x2 pixel tiles)
* ----------------------------------------------------------------------- */

static void rasterizeTile16(const TriSetup* ts, const uint32_t* texture,
int tx, int ty)
{
int y0 = ty;
int y1 = ty + 1;

/* Early reject if tile is entirely outside triangle */
if (y0 > ts->ymax || y1 > ts->ymax || y0 < ts->ymin) return;

/* Build pixel coordinate vectors:
* px = [tx, tx+1, ..., tx+7, tx, tx+1, ..., tx+7]
* py = [y0, y0, ..., y0, y1, y1, ..., y1]
*/
int32_t px_arr[16], py_arr[16];
for (int i = 0; i < 8; i++) {
px_arr = tx + i;
px_arr[i + 8] = tx + i;
py_arr = y0;
py_arr[i + 8] = y1;
}

v16i px_i = _mm512_loadu_epi32(px_arr);
v16i py_i = _mm512_loadu_epi32(py_arr);

/* Offsets from bounding box origin (in pixels) */
v16i xmin_i = _mm512_set1_epi32(ts->xmin);
v16i ymin_i = _mm512_set1_epi32(ts->ymin);
v16i dx_pix = _mm512_sub_epi32(px_i, xmin_i);
v16i dy_pix = _mm512_sub_epi32(py_i, ymin_i);

/* ---- Edge evaluation (28.4 fixed-point) ----
* e = e_row0 + e_dx * dx_pix + e_dy * dy_pix
* All values in 28.4^2 units; sign determines coverage.
*/
v16i e01 = _mm512_add_epi32(
_mm512_set1_epi32(ts->e01_row0),
_mm512_add_epi32(
_mm512_mullo_epi32(_mm512_set1_epi32(ts->e01_dx), dx_pix),
_mm512_mullo_epi32(_mm512_set1_epi32(ts->e01_dy), dy_pix)));

v16i e12 = _mm512_add_epi32(
_mm512_set1_epi32(ts->e12_row0),
_mm512_add_epi32(
_mm512_mullo_epi32(_mm512_set1_epi32(ts->e12_dx), dx_pix),
_mm512_mullo_epi32(_mm512_set1_epi32(ts->e12_dy), dy_pix)));

v16i e20 = _mm512_add_epi32(
_mm512_set1_epi32(ts->e20_row0),
_mm512_add_epi32(
_mm512_mullo_epi32(_mm512_set1_epi32(ts->e20_dx), dx_pix),
_mm512_mullo_epi32(_mm512_set1_epi32(ts->e20_dy), dy_pix)));

/* ---- Coverage mask (all edges >= 0) ---- */
v16i zero = _mm512_setzero_epi32();
v16i inside = _mm512_and_epi32(
_mm512_and_epi32(
_mm512_cmpge_epi32(e01, zero),
_mm512_cmpge_epi32(e12, zero)),
_mm512_cmpge_epi32(e20, zero));

unsigned short mask = get_mask(inside);
if (mask == 0) return;

/* ---- Clamp pixel coordinates to screen bounds ---- */
v16i fb_w_max = _mm512_set1_epi32(ts->fb_stride - 1);
v16i ymax_i = _mm512_set1_epi32(ts->ymax);
px_i = _mm512_max_epi32(_mm512_min_epi32(px_i, fb_w_max), zero);
py_i = _mm512_max_epi32(_mm512_min_epi32(py_i, ymax_i), zero);

/* ---- Shade covered pixels ---- */
shadePixels16(ts, texture, dx_pix, dy_pix, px_i, py_i, mask);
}

/* -----------------------------------------------------------------------
* §5 Entry Point
* ----------------------------------------------------------------------- */

void rasterizeTriangle(const Vertex* v0, const Vertex* v1, const Vertex* v2,
uint32_t* fb, int fb_w, int fb_h,
const uint32_t* tex, int tex_w, int tex_h)
{
TriSetup ts;
if (!triSetup(v0, v1, v2, fb_w, fb_h, tex_w, tex_h, fb, &ts))
return;

/* Walk 8x2 tiles aligned to bounding box */
int tx0 = ts.xmin & ~7;
int ty0 = ts.ymin & ~1;
int tx1 = (ts.xmax + 7) & ~7;
int ty1 = (ts.ymax + 1) & ~1;

for (int ty = ty0; ty < ty1; ty += 2) {
for (int tx = tx0; tx < tx1; tx += 8) {
rasterizeTile16(&ts, tex, tx, ty);
}
}
}
```

## Summary of corrections applied

| Issue | Previous version | This version |
|-------|-----------------|--------------|
| Edge derivatives | `e_dx = (v0y - v1y)` (missing ×16) | `e_dx = (v0y - v1y) * 16` (correct per-pixel step) |
| Initial edge value | Mixed 28.4 and pixel units | Pure 28.4: `(sy - v0y)*(v1x - v0x) - (sx - v0x)*(v1y - v0y)` |
| Gather return type | Float with reinterpret cast | Integer gather `_mm512_mask_i32gather_epi32` |
| Framebuffer store | Two 128-bit stores (AVX-512F) | Scalar fallback (safe, verifiable) |
| Texture modulo | `_mm512_rem_epi32` (fantasy) | Power-of-two AND; non-pow2 clamp |
| Intrinsic naming | Mixed `_si512` and `_epi32` | Consistent `_epi32` suffix |
| Overflow handling | 32-bit area | 64-bit area in setup |

## Remaining work (requires actual ICC 12.x)

1. **Verify every intrinsic name** against `zmmintrin.h`
2. **Test edge rule** (top-left fill convention, shared edge handling)
3. **Benchmark scalar framebuffer store** vs masked store (if `_mm512_mask_storeu_epi32` exists)
4. **Test overflow behavior** for large triangles on high-res framebuffers
5. **Add sub-pixel precision** for vertex positions (currently truncates, should round)

The architecture is now correct. The implementation is as close to production-ready as possible without access to the actual compiler and hardware.

 

[erek]

/* ============================================================
 * lrb_knf_renderer_v3.c — Larrabee KNF Software Renderer (Rev 3)
 *
 * FUNDAMENTAL FIX: Separate lane assignments for coarse vs fine.
 *
 *   Coarse traversal: lanes = 16 TRIANGLES
 *     - Edge coefficients stored SoA across triangles
 *     - Evaluate all 16 tris' edges at a block corner in 3 FMA
 *     - 4 corners → fully_in / fully_out masks (1 bit per tri)
 *     - Amortizes hierarchical traversal across the packet
 *
 *   Fine rasterization: lanes = 16 PIXELS (4x4 block)
 *     - ONE triangle at a time, edge coeffs BROADCAST to all lanes
 *     - Pixel coords in lanes, evaluate e(lane) = A*px + B*py + C
 *     - Inside test → 16-bit coverage mask
 *     - Interpolate depth/color for 16 pixels via FMA
 *     - Depth test → predicated store to swizzled tile buffer
 *
 * Incremental stepping is PER-TRIANGLE SCALAR:
 *   e_base += A * 4  (scalar add per tri per step)
 *   Per-pixel: e(lane) = fma(bcast(A), lane_xoff, fma(bcast(B), lane_yoff, bcast(e_base)))
 *
 * This is the correct Larrabee programming model:
 *   - Packets amortize coarse traversal (lanes=tris)
 *   - Fine raster always uses full 16-wide pixel parallelism (lanes=pixels)
 *   - No scalar extraction + vector rebuild in the inner loop
 *
 * Build:
 *   gcc -O3 -march=knl -mavx512f -mavx512cd -mavx512pf -mavx512er \
 *       -mfma -lm -lpthread lrb_knf_renderer_v3.c -o lrb_knf_v3
 * ============================================================ */

#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdint.h>
#include <pthread.h>
#include <sched.h>
#include <immintrin.h>

/* ============================================================
 * 1. Architecture Constants
 ============================================================ */

#define LRB_CORES            32
#define LRB_THREADS_PER_CORE 4
#define LRB_TOTAL_THREADS    (LRB_CORES * LRB_THREADS_PER_CORE)
#define LRB_VPU_WIDTH       16
#define LRB_L2_BYTES         (256 * 1024)
#define LRB_TILE_SIZE        64
#define LRB_FINE_BLOCK       4       /* 4x4 px = 16 pixels = 1 vector */
#define LRB_COARSE_BLOCK     16      /* 16x16 px = 4x4 array of fine blocks */
#define LRB_ALIGN            64
#define LRB_PACKET_TRIS      16      /* triangles per coarse packet */

/* Tile working set: 64x64 = 4096 px
 *   Depth: 16 KB, Color: 64 KB, Total: 80 KB < 256 KB L2 */

/* ============================================================
 * 2. LRBni Intrinsic Layer
 * ============================================================ */

typedef __m512    lrb_v;
typedef __m512i   lrb_vi;
typedef __mmask16 lrb_k;

static inline lrb_v vadd(lrb_v a, lrb_v b) { return _mm512_add_ps(a, b); }
static inline lrb_v vsub(lrb_v a, lrb_v b) { return _mm512_sub_ps(a, b); }
static inline lrb_v vmul(lrb_v a, lrb_v b) { return _mm512_mul_ps(a, b); }
static inline lrb_v vfma(lrb_v a, lrb_v b, lrb_v c) {
    return _mm512_fmadd_ps(a, b, c);
}
static inline lrb_v vfnma(lrb_v a, lrb_v b, lrb_v c) {
    return _mm512_fnmadd_ps(a, b, c);
}
static inline lrb_v vbcast(float s) { return _mm512_set1_ps(s); }
static inline lrb_v vload(const float *p) { return _mm512_load_ps(p); }
static inline lrb_v vloadu(const float *p) { return _mm512_loadu_ps(p); }
static inline void vstore(float *p, lrb_v v) { _mm512_store_ps(p, v); }
static inline void vstore_k(float *p, lrb_v v, lrb_k k) {
    _mm512_mask_store_ps(p, k, v);
}
static inline lrb_k kge0(lrb_v a) {
    return _mm512_cmp_ps_mask(a, _mm512_setzero_ps(), _CMP_GE_OQ);
}
static inline lrb_k kle0(lrb_v a) {
    return _mm512_cmp_ps_mask(a, _mm512_setzero_ps(), _CMP_LE_OQ);
}
static inline lrb_k klt(lrb_v a, lrb_v b) {
    return _mm512_cmp_ps_mask(a, b, _CMP_LT_OQ);
}
static inline lrb_k kand(lrb_k a, lrb_k b) { return _kand_mask16(a, b); }
static inline lrb_k kor(lrb_k a, lrb_k b) { return _kor_mask16(a, b); }

static inline lrb_v vrcp(lrb_v x) {
    lrb_v r0 = _mm512_rcp14_ps(x);
    lrb_v nr = vfnma(r0, x, vbcast(2.0f));
    return vmul(r0, nr);
}

static inline void prefetch_l2(const void *p) {
    _mm_prefetch((const char *)p, _MM_HINT_T1);
}

/* Extract lane i from a vector (for transitioning coarse→fine) */
static inline float vlane(lrb_v v, int i) {
    return ((float *)&v)[i];
}

/* ============================================================
 * 3. Swizzled Tile-Local Framebuffer
 * ============================================================ */

#define TILE_BX     (LRB_TILE_SIZE / LRB_FINE_BLOCK)  /* 16 */
#define TILE_BY     (LRB_TILE_SIZE / LRB_FINE_BLOCK)  /* 16 */
#define TILE_BLOCKS (TILE_BX * TILE_BY)               /* 256 */
#define TILE_PIX    (LRB_TILE_SIZE * LRB_TILE_SIZE)   /* 4096 */

static inline int swizzle(int x, int y) {
    int bx = x >> 2, by = y >> 2;
    int lx = x & 3,  ly = y & 3;
    return (by * TILE_BX + bx) * 16 + ly * 4 + lx;
}

typedef struct {
    alignas(64) float depth[TILE_PIX];
    alignas(64) float color_r[TILE_PIX];
    alignas(64) float color_g[TILE_PIX];
    alignas(64) float color_b[TILE_PIX];
    alignas(64) float color_a[TILE_PIX];
} lrb_tile_buf;  /* 80 KB */

/* ============================================================
 * 4. Data Structures
 * ============================================================ */

typedef struct { float x, y, z, w; } vec4;
typedef struct { float u, v;         } vec2;

typedef struct {
    alignas(64) float x[16], y[16], z[16], w[16];
    alignas(64) float r[16], g[16], b[16], a[16];
    alignas(64) float u[16], v[16];
} lrb_vertex_soa;

typedef struct {
    vec4 pos;
    vec4 color;
    vec2 uv;
} lrb_vertex;

typedef struct { int v[3]; } lrb_tri;

/* ============================================================
 * 5. Triangle Packet — Coarse Traversal (lanes = triangles)
 * ============================================================
 *
 * Edge coefficients stored SoA across 16 triangles.
 * Used ONLY for coarse block testing.
 * Fine rasterization extracts per-triangle scalars and broadcasts
 * them to pixel lanes.
 * ============================================================ */

typedef struct {
    /* Edge coefficients: e(x,y) = A*x + B*y + C
     * Lane i = triangle i's coefficient */
    alignas(64) float A0[16], B0[16], C0[16];
    alignas(64) float A1[16], B1[16], C1[16];
    alignas(64) float A2[16], B2[16], C2[16];

    alignas(64) float inv_area[16];

    /* Per-triangle interpolation data (scalar, extracted for fine raster) */
    alignas(64) float z0[16], z1[16], z2[16];
    alignas(64) float r0[16], r1[16], r2[16];
    alignas(64) float g0[16], g1[16], g2[16];
    alignas(64) float b0[16], b1[16], b2[16];

    /* Bounding boxes (for tile binning) */
    alignas(64) float minx[16], maxx[16];
    alignas(64) float miny[16], maxy[16];

    /* Masks: 1 bit per triangle */
    lrb_k valid_mask;   /* non-degenerate */
    lrb_k ccw_mask;     /* CCW winding */

    int count;
} lrb_tri_packet;

/* ============================================================
 * 6. Per-Pixel Lane Offsets (for fine raster, lanes = pixels)
 * ============================================================
 * Lane = ly*4 + lx for a 4x4 block.
 * xoff[lane] = lx, yoff[lane] = ly.
 * Same for every block — computed once.
 * ============================================================ */

static lrb_v lane_xoff, lane_yoff;

static void init_lane_offsets(void) {
    alignas(64) float xo[16], yo[16];
    for (int ly = 0; ly < 4; ly++)
        for (int lx = 0; lx < 4; lx++) {
            xo[ly * 4 + lx] = (float)lx;
            yo[ly * 4 + lx] = (float)ly;
        }
    lane_xoff = vload(xo);
    lane_yoff = vload(yo);
}

/* ============================================================
 * 7. Matrix Math
 * ============================================================ */

typedef struct { float m[4][4]; } lrb_mat4;

static lrb_mat4 mat4_id(void) {
    lrb_mat4 r; memset(&r, 0, sizeof(r));
    r.m[0][0] = r.m[1][1] = r.m[2][2] = r.m[3][3] = 1.0f;
    return r;
}
static lrb_mat4 mat4_mul(lrb_mat4 a, lrb_mat4 b) {
    lrb_mat4 r;
    for (int i = 0; i < 4; i++)
        for (int j = 0; j < 4; j++) {
            r.m[i][j] = 0;
            for (int k = 0; k < 4; k++)
                r.m[i][j] += a.m[i][k] * b.m[k][j];
        }
    return r;
}
static lrb_mat4 mat4_persp(float fov, float asp, float zn, float zf) {
    lrb_mat4 r = mat4_id();
    float f = 1.0f / tanf(fov * 3.14159265f / 360.0f);
    r.m[0][0] = f / asp;  r.m[1][1] = f;
    r.m[2][2] = (zf + zn) / (zn - zf);
    r.m[2][3] = -1.0f;
    r.m[3][2] = (2.0f * zf * zn) / (zn - zf);
    r.m[3][3] = 0.0f;
    return r;
}
static lrb_mat4 mat4_trans(float x, float y, float z) {
    lrb_mat4 r = mat4_id();
    r.m[3][0] = x; r.m[3][1] = y; r.m[3][2] = z;
    return r;
}
static lrb_mat4 mat4_roty(float a) {
    lrb_mat4 r = mat4_id();
    float c = cosf(a), s = sinf(a);
    r.m[0][0] = c;  r.m[0][2] = s;
    r.m[2][0] = -s; r.m[2][2] = c;
    return r;
}

/* ============================================================
 * 8. AoS → SoA Transpose
 * ============================================================ */

static void aos_to_soa(const lrb_vertex *src, lrb_vertex_soa *dst, int n) {
    for (int i = 0; i < 16; i++) {
        if (i < n) {
            dst->x[i] = src[i].pos.x;  dst->y[i] = src[i].pos.y;
            dst->z[i] = src[i].pos.z;  dst->w[i] = src[i].pos.w;
            dst->r[i] = src[i].color.x; dst->g[i] = src[i].color.y;
            dst->b[i] = src[i].color.z; dst->a[i] = src[i].color.w;
            dst->u[i] = src[i].uv.u;    dst->v[i] = src[i].uv.v;
        } else {
            dst->x[i] = dst->y[i] = dst->z[i] = 0.0f;
            dst->w[i] = 1.0f;
            dst->r[i] = dst->g[i] = dst->b[i] = 0.0f;
            dst->a[i] = 1.0f;
            dst->u[i] = dst->v[i] = 0.0f;
        }
    }
}

/* ============================================================
 * 9. Vectorized Vertex Shader (lanes = vertices)
 * ============================================================ */

typedef struct {
    lrb_mat4 mvp;
    int fb_w, fb_h;
} lrb_shader_ctx;

static void vertex_shader(const lrb_vertex_soa *in, lrb_vertex_soa *out,
                          const lrb_shader_ctx *ctx)
{
    lrb_v vx = vload(in->x), vy = vload(in->y);
    lrb_v vz = vload(in->z), vw = vload(in->w);
    const lrb_mat4 *m = &ctx->mvp;

    lrb_v ox = vfma(vbcast(m->m[0][0]), vx, vbcast(0.0f));
    ox = vfma(vbcast(m->m[1][0]), vy, ox);
    ox = vfma(vbcast(m->m[2][0]), vz, ox);
    ox = vfma(vbcast(m->m[3][0]), vw, ox);

    lrb_v oy = vfma(vbcast(m->m[0][1]), vx, vbcast(0.0f));
    oy = vfma(vbcast(m->m[1][1]), vy, oy);
    oy = vfma(vbcast(m->m[2][1]), vz, oy);
    oy = vfma(vbcast(m->m[3][1]), vw, oy);

    lrb_v oz = vfma(vbcast(m->m[0][2]), vx, vbcast(0.0f));
    oz = vfma(vbcast(m->m[1][2]), vy, oz);
    oz = vfma(vbcast(m->m[2][2]), vz, oz);
    oz = vfma(vbcast(m->m[3][2]), vw, oz);

    lrb_v ow = vfma(vbcast(m->m[0][3]), vx, vbcast(0.0f));
    ow = vfma(vbcast(m->m[1][3]), vy, ow);
    ow = vfma(vbcast(m->m[2][3]), vz, ow);
    ow = vfma(vbcast(m->m[3][3]), vw, ow);

    lrb_v rcp_w = vrcp(ow);
    ox = vmul(ox, rcp_w);
    oy = vmul(oy, rcp_w);
    oz = vmul(oz, rcp_w);

    lrb_v half = vbcast(0.5f);
    ox = vmul(vfma(ox, half, half), vbcast((float)ctx->fb_w));
    oy = vmul(vfma(oy, half, half), vbcast((float)ctx->fb_h));

    vstore(out->x, ox);
    vstore(out->y, oy);
    vstore(out->z, oz);
    vstore(out->w, vbcast(1.0f));

    vstore(out->r, vload(in->r));
    vstore(out->g, vload(in->g));
    vstore(out->b, vload(in->b));
    vstore(out->a, vload(in->a));
    vstore(out->u, vload(in->u));
    vstore(out->v, vload(in->v));
}

/* ============================================================
 * 10. Triangle Setup → Packet
 * ============================================================ */

typedef struct {
    float x0, y0, z0, x1, y1, z1, x2, y2, z2;
    vec4 c0, c1, c2;
} lrb_setup_input;

static void build_packet(const lrb_setup_input *tris, int n,
                         lrb_tri_packet *pkt)
{
    memset(pkt, 0, sizeof(*pkt));
    pkt->count = (n < LRB_PACKET_TRIS) ? n : LRB_PACKET_TRIS;

    lrb_k valid = 0, ccw = 0;

    for (int i = 0; i < LRB_PACKET_TRIS; i++) {
        if (i >= n) {
            /* Pad with degenerate triangles */
            pkt->A0[i] = pkt->B0[i] = pkt->C0[i] = 0.0f;
            pkt->A1[i] = pkt->B1[i] = pkt->C1[i] = 0.0f;
            pkt->A2[i] = pkt->B2[i] = pkt->C2[i] = 0.0f;
            pkt->inv_area[i] = 0.0f;
            pkt->minx[i] = pkt->maxx[i] = 0.0f;
            pkt->miny[i] = pkt->maxy[i] = 0.0f;
            continue;
        }

        const lrb_setup_input *t = &tris[i];

        /* Edge 0: v0→v1 */
        pkt->A0[i] = t->y1 - t->y0;
        pkt->B0[i] = t->x0 - t->x1;
        pkt->C0[i] = t->x1 * t->y0 - t->x0 * t->y1;

        /* Edge 1: v1→v2 */
        pkt->A1[i] = t->y2 - t->y1;
        pkt->B1[i] = t->x1 - t->x2;
        pkt->C1[i] = t->x2 * t->y1 - t->x1 * t->y2;

        /* Edge 2: v2→v0 */
        pkt->A2[i] = t->y0 - t->y2;
        pkt->B2[i] = t->x2 - t->x0;
        pkt->C2[i] = t->x0 * t->y2 - t->x2 * t->y0;

        float area = pkt->A0[i] * t->x2 + pkt->B0[i] * t->y2 + pkt->C0[i];

        if (fabsf(area) > 1e-8f) {
            pkt->inv_area[i] = 1.0f / area;
            valid |= (1u << i);
            if (area > 0.0f) ccw |= (1u << i);
        } else {
            pkt->inv_area[i] = 0.0f;
        }

        pkt->z0[i] = t->z0; pkt->z1[i] = t->z1; pkt->z2[i] = t->z2;
        pkt->r0[i] = t->c0.x; pkt->r1[i] = t->c1.x; pkt->r2[i] = t->c2.x;
        pkt->g0[i] = t->c0.y; pkt->g1[i] = t->c1.y; pkt->g2[i] = t->c2.y;
        pkt->b0[i] = t->c0.z; pkt->b1[i] = t->c1.z; pkt->b2[i] = t->c2.z;

        pkt->minx[i] = fminf(fminf(t->x0, t->x1), t->x2);
        pkt->maxx[i] = fmaxf(fmaxf(t->x0, t->x1), t->x2);
        pkt->miny[i] = fminf(fminf(t->y0, t->y1), t->y2);
        pkt->maxy[i] = fmaxf(fmaxf(t->y0, t->y1), t->y2);
    }

    pkt->valid_mask = valid;
    pkt->ccw_mask = ccw;
}

/* ============================================================
 * 11. Coarse Traversal — lanes = TRIANGLES
 * ============================================================
 *
 * Evaluate all 16 triangles' edges at a scalar point (x,y).
 * Returns 3 vectors: ev0, ev1, ev2 (lane i = tri i's edge value).
 * 3 FMA per edge × 3 edges = 9 FMA for ALL 16 triangles.
 * ============================================================ */

static inline void eval_packet_at_point(
    lrb_v pA0, lrb_v pB0, lrb_v pC0,
    lrb_v pA1, lrb_v pB1, lrb_v pC1,
    lrb_v pA2, lrb_v pB2, lrb_v pC2,
    float x, float y,
    lrb_v *ev0, lrb_v *ev1, lrb_v *ev2)
{
    lrb_v vx = vbcast(x), vy = vbcast(y);
    *ev0 = vfma(pA0, vx, vfma(pB0, vy, pC0));
    *ev1 = vfma(pA1, vx, vfma(pB1, vy, pC1));
    *ev2 = vfma(pA2, vx, vfma(pB2, vy, pC2));
}

/* Coarse test: classify 16x16 block for all 16 tris.
 * Returns fully_in / fully_out masks (1 bit per tri). */
static void coarse_test_packet(
    lrb_v pA0, lrb_v pB0, lrb_v pC0,
    lrb_v pA1, lrb_v pB1, lrb_v pC1,
    lrb_v pA2, lrb_v pB2, lrb_v pC2,
    lrb_k ccw_mask, lrb_k valid_mask,
    int bx, int by,
    lrb_k *fully_in, lrb_k *fully_out)
{
    float x0 = (float)bx,         y0 = (float)by;
    float x1 = (float)(bx + 16),  y1 = (float)(by + 16);

    lrb_v e0_00, e1_00, e2_00, e0_10, e1_10, e2_10;
    lrb_v e0_01, e1_01, e2_01, e0_11, e1_11, e2_11;

    eval_packet_at_point(pA0,pB0,pC0, pA1,pB1,pC1, pA2,pB2,pC2,
                         x0, y0, &e0_00, &e1_00, &e2_00);
    eval_packet_at_point(pA0,pB0,pC0, pA1,pB1,pC1, pA2,pB2,pC2,
                         x1, y0, &e0_10, &e1_10, &e2_10);
    eval_packet_at_point(pA0,pB0,pC0, pA1,pB1,pC1, pA2,pB2,pC2,
                         x0, y1, &e0_01, &e1_01, &e2_01);
    eval_packet_at_point(pA0,pB0,pC0, pA1,pB1,pC1, pA2,pB2,pC2,
                         x1, y1, &e0_11, &e1_11, &e2_11);

    /* CCW: inside = all edges >= 0 at all 4 corners */
    lrb_k ccw = ccw_mask;
    lrb_k cw  = ~ccw & valid_mask;

    /* CCW fully-in: all 4 corners, all 3 edges >= 0 */
    lrb_k ccw_e0_all = kand(kand(kge0(e0_00), kge0(e0_10)),
                            kand(kge0(e0_01), kge0(e0_11)));
    lrb_k ccw_e1_all = kand(kand(kge0(e1_00), kge0(e1_10)),
                            kand(kge0(e1_01), kge0(e1_11)));
    lrb_k ccw_e2_all = kand(kand(kge0(e2_00), kge0(e2_10)),
                            kand(kge0(e2_01), kge0(e2_11)));
    lrb_k ccw_in = kand(kand(ccw_e0_all, ccw_e1_all), ccw_e2_all);
    ccw_in = kand(ccw_in, ccw);

    /* CW fully-in: all 4 corners, all 3 edges <= 0 */
    lrb_k cw_e0_all = kand(kand(kle0(e0_00), kle0(e0_10)),
                           kand(kle0(e0_01), kle0(e0_11)));
    lrb_k cw_e1_all = kand(kand(kle0(e1_00), kle0(e1_10)),
                           kand(kle0(e1_01), kle0(e1_11)));
    lrb_k cw_e2_all = kand(kand(kle0(e2_00), kle0(e2_10)),
                           kand(kle0(e2_01), kle0(e2_11)));
    lrb_k cw_in = kand(kand(cw_e0_all, cw_e1_all), cw_e2_all);
    cw_in = kand(cw_in, cw);

 

[erek]

    *fully_in = kor(ccw_in, cw_in);

    /* CCW fully-out: at least one edge < 0 at ALL 4 corners */
    lrb_k ccw_out_e0 = kand(kand(kle0(e0_00), kle0(e0_10)),
                            kand(kle0(e0_01), kle0(e0_11)));
    lrb_k ccw_out_e1 = kand(kand(kle0(e1_00), kle0(e1_10)),
                            kand(kle0(e1_01), kle0(e1_11)));
    lrb_k ccw_out_e2 = kand(kand(kle0(e2_00), kle0(e2_10)),
                            kand(kle0(e2_01), kle0(e2_11)));
    lrb_k ccw_out = kor(kor(ccw_out_e0, ccw_out_e1), ccw_out_e2);
    ccw_out = kand(ccw_out, ccw);

    /* CW fully-out: at least one edge > 0 at ALL 4 corners */
    lrb_k cw_out_e0 = kand(kand(kge0(e0_00), kge0(e0_10)),
                           kand(kge0(e0_01), kge0(e0_11)));
    lrb_k cw_out_e1 = kand(kand(kge0(e1_00), kge0(e1_10)),
                           kand(kge0(e1_01), kge0(e1_11)));
    lrb_k cw_out_e2 = kand(kand(kge0(e2_00), kge0(e2_10)),
                           kand(kge0(e2_01), kge0(e2_11)));
    lrb_k cw_out = kor(kor(cw_out_e0, cw_out_e1), cw_out_e2);
    cw_out = kand(cw_out, cw);

    *fully_out = kor(ccw_out, cw_out);
    *fully_out = kand(*fully_out, valid_mask);
}

/* ============================================================
 * 12. Fine Rasterization — lanes = PIXELS
 * ============================================================
 *
 * Process ONE triangle across a 4x4 block (16 pixels).
 * Edge coefficients are SCALARS extracted from the packet,
 * BROADCAST to all 16 pixel lanes.
 *
 * e(lane) = fma(bcast(A), pix_x, fma(bcast(B), pix_y, bcast(C)))
 *
 * Incremental: base edge value stepped as SCALAR per tri:
 *   e_base_x += A * 4   (scalar add)
 *   e_base_y += B * 4   (scalar add)
 * Then per-pixel: e(lane) = fma(bcast(A), lane_xoff,
 *                               fma(bcast(B), lane_yoff, bcast(e_base)))
 *
 * This keeps the VPU fully utilized on 16 pixels for each tri.
 * ============================================================ */

static inline void rasterize_tri_fine(
    lrb_tile_buf *tb,
    int swiz_off,
    /* This triangle's edge coefficients (scalars, broadcast to lanes) */
    float A0, float B0, float C0,
    float A1, float B1, float C1,
    float A2, float B2, float C2,
    float inv_area,
    float z0, float z1, float z2,
    float r0, float r1, float r2,
    float g0, float g1, float g2,
    float b0, float b1, float b2,
    /* Base edge values at block origin (scalars) */
    float e0_base, float e1_base, float e2_base,
    int is_ccw)
{
    /* Broadcast edge coefficients to all 16 pixel lanes */
    lrb_v vA0 = vbcast(A0), vB0 = vbcast(B0);
    lrb_v vA1 = vbcast(A1), vB1 = vbcast(B1);
    lrb_v vA2 = vbcast(A2), vB2 = vbcast(B2);

    /* Per-pixel edge values:
     * e(lane) = e_base + A*lane_xoff + B*lane_yoff
     *         = fma(A, lane_xoff, fma(B, lane_yoff, bcast(e_base))) */
    lrb_v ev0 = vfma(vA0, lane_xoff, vfma(vB0, lane_yoff, vbcast(e0_base)));
    lrb_v ev1 = vfma(vA1, lane_xoff, vfma(vB1, lane_yoff, vbcast(e1_base)));
    lrb_v ev2 = vfma(vA2, lane_xoff, vfma(vB2, lane_yoff, vbcast(e2_base)));

    /* Inside test → 16-bit coverage mask */
    lrb_k coverage;
    if (is_ccw)
        coverage = kand(kand(kge0(ev0), kge0(ev1)), kge0(ev2));
    else
        coverage = kand(kand(kle0(ev0), kle0(ev1)), kle0(ev2));

    if (coverage == 0) return;

    /* Barycentric: w0 = e1/area, w1 = e2/area, w2 = 1 - w0 - w1 */
    lrb_v v_inv = vbcast(inv_area);
    lrb_v bw0 = vmul(ev1, v_inv);
    lrb_v bw1 = vmul(ev2, v_inv);
    lrb_v bw2 = vsub(vsub(vbcast(1.0f), bw0), bw1);

    /* Interpolate depth & color via FMA (lanes = pixels) */
    lrb_v depth = vfma(bw0, vbcast(z0),
                  vfma(bw1, vbcast(z1), vmul(bw2, vbcast(z2))));
    lrb_v cr = vfma(bw0, vbcast(r0),
               vfma(bw1, vbcast(r1), vmul(bw2, vbcast(r2))));
    lrb_v cg = vfma(bw0, vbcast(g0),
               vfma(bw1, vbcast(g1), vmul(bw2, vbcast(g2))));
    lrb_v cb = vfma(bw0, vbcast(b0),
               vfma(bw1, vbcast(b1), vmul(bw2, vbcast(b2))));

    /* Depth test (lanes = pixels) */
    lrb_v db = vloadu(&tb->depth[swiz_off]);
    lrb_k write = kand(coverage, klt(depth, db));

    if (write == 0) return;

    /* Predicated stores to swizzled tile buffer */
    vstore_k(&tb->depth[swiz_off],   depth, write);
    vstore_k(&tb->color_r[swiz_off], cr,    write);
    vstore_k(&tb->color_g[swiz_off], cg,    write);
    vstore_k(&tb->color_b[swiz_off], cb,    write);
    vstore_k(&tb->color_a[swiz_off], vbcast(1.0f), write);
}

/* ============================================================
 * 13. Packet Rasterizer — Coarse (tris) → Fine (pixels)
 * ============================================================
 *
 * For each coarse block:
 *   1. Coarse test all 16 tris (lanes = tris) → fully_in/fully_out masks
 *   2. For each tri that is not fully_out:
 *      a. If fully_in: rasterize all 16 fine blocks at full coverage
 *      b. If mixed: descend, test each 4x4 block (lanes = pixels)
 *
 * Incremental stepping is PER-TRIANGLE SCALAR:
 *   e_base_x += A * 4, e_base_y += B * 4
 * ============================================================ */

static void rasterize_packet_tile(lrb_tile_buf *tb,
                                  const lrb_tri_packet *pkt,
                                  int tile_x, int tile_y,
                                  int fb_w, int fb_h)
{
    int tx0 = tile_x * LRB_TILE_SIZE;
    int ty0 = tile_y * LRB_TILE_SIZE;

    /* Packet bounding box clipped to tile */
    alignas(64) float minx_a[16], maxx_a[16], miny_a[16], maxy_a[16];
    vstore(minx_a, vload(pkt->minx));
    vstore(maxx_a, vload(pkt->maxx));
    vstore(miny_a, vload(pkt->miny));
    vstore(maxy_a, vload(pkt->maxy));

    float pb_minx = 1e30f, pb_maxx = -1e30f;
    float pb_miny = 1e30f, pb_maxy = -1e30f;
    for (int i = 0; i < pkt->count; i++) {
        if (pkt->valid_mask & (1u << i)) {
            pb_minx = fminf(pb_minx, minx_a[i]);
            pb_maxx = fmaxf(pb_maxx, maxx_a[i]);
            pb_miny = fminf(pb_miny, miny_a[i]);
            pb_maxy = fmaxf(pb_maxy, maxy_a[i]);
        }
    }
    if (pb_maxx < tx0 || pb_minx > tx0 + LRB_TILE_SIZE - 1 ||
        pb_maxy < ty0 || pb_miny > ty0 + LRB_TILE_SIZE - 1)
        return;

    int bx0 = (int)fmaxf(pb_minx, (float)tx0) & ~(LRB_COARSE_BLOCK - 1);
    int by0 = (int)fmaxf(pb_miny, (float)ty0) & ~(LRB_COARSE_BLOCK - 1);
    int bx1 = (int)fminf(pb_maxx, (float)(tx0 + LRB_TILE_SIZE - 1));
    int by1 = (int)fminf(pb_maxy, (float)(ty0 + LRB_TILE_SIZE - 1));

    /* Load packet edge coefficients (lanes = tris) — stay in registers */
    lrb_v pA0 = vload(pkt->A0), pB0 = vload(pkt->B0), pC0 = vload(pkt->C0);
    lrb_v pA1 = vload(pkt->A1), pB1 = vload(pkt->B1), pC1 = vload(pkt->C1);
    lrb_v pA2 = vload(pkt->A2), pB2 = vload(pkt->B2), pC2 = vload(pkt->C2);

    lrb_k valid = pkt->valid_mask;
    lrb_k ccw   = pkt->ccw_mask;

    /* --- Coarse traversal: 16x16 blocks --- */
    for (int cy = by0; cy <= by1; cy += LRB_COARSE_BLOCK) {
        for (int cx = bx0; cx <= bx1; cx += LRB_COARSE_BLOCK) {

            /* Coarse test: all 16 tris at once (lanes = tris) */
            lrb_k fully_in, fully_out;
            coarse_test_packet(pA0,pB0,pC0, pA1,pB1,pC1, pA2,pB2,pC2,
                               ccw, valid, cx, cy,
                               &fully_in, &fully_out);

            /* Skip block if ALL tris reject */
            if ((fully_out & valid) == valid) continue;

            /* Active tris = valid tris that are not fully out */
            lrb_k active = valid & ~fully_out;

            /* For each active triangle, do fine rasterization (lanes = pixels) */
            while (active != 0) {
                int ti = __builtin_ctz(active);
                active &= ~(1u << ti);

                int is_ccw = (ccw & (1u << ti)) != 0;

                /* Extract this triangle's scalar edge coefficients */
                float A0 = vlane(pA0, ti), B0 = vlane(pB0, ti), C0 = vlane(pC0, ti);
                float A1 = vlane(pA1, ti), B1 = vlane(pB1, ti), C1 = vlane(pC1, ti);
                float A2 = vlane(pA2, ti), B2 = vlane(pB2, ti), C2 = vlane(pC2, ti);
                float inv_a = pkt->inv_area[ti];

                float z0 = pkt->z0[ti], z1 = pkt->z1[ti], z2 = pkt->z2[ti];
                float r0 = pkt->r0[ti], r1 = pkt->r1[ti], r2 = pkt->r2[ti];
                float g0 = pkt->g0[ti], g1 = pkt->g1[ti], g2 = pkt->g2[ti];
                float b0 = pkt->b0[ti], b1 = pkt->b1[ti], b2 = pkt->b2[ti];

                /* Scalar incremental stepping deltas */
                float dA0x = A0 * 4.0f, dA1x = A1 * 4.0f, dA2x = A2 * 4.0f;
                float dB0y = B0 * 4.0f, dB1y = B1 * 4.0f, dB2y = B2 * 4.0f;

                /* Base edge values at coarse block origin (cx, cy) */
                float e0_row = A0 * ((float)cx + 0.5f) + B0 * ((float)cy + 0.5f) + C0;
                float e1_row = A1 * ((float)cx + 0.5f) + B1 * ((float)cy + 0.5f) + C1;
                float e2_row = A2 * ((float)cx + 0.5f) + B2 * ((float)cy + 0.5f) + C2;

                int tri_fully_in = (fully_in & (1u << ti)) != 0;

                /* --- Fine traversal: 4x4 blocks within 16x16 --- */
                for (int fy = 0; fy < LRB_COARSE_BLOCK; fy += LRB_FINE_BLOCK) {
                    float e0 = e0_row, e1 = e1_row, e2 = e2_row;

                    for (int fx = 0; fx < LRB_COARSE_BLOCK; fx += LRB_FINE_BLOCK) {
                        int px = cx + fx, py = cy + fy;
                        if (px >= fb_w || py >= fb_h) {
                            e0 += dA0x; e1 += dA1x; e2 += dA2x;
                            continue;
                        }

                        int lx = px - tx0, ly = py - ty0;
                        int swiz = swizzle(lx, ly);

                        if (tri_fully_in) {
                            /* Fully covered — skip edge test, use full mask.
                             * Build coverage = 0xFFFF directly, interpolate,
                             * depth test, and store. */
                            lrb_v vA0_ = vbcast(A0), vB0_ = vbcast(B0);
                            lrb_v vA1_ = vbcast(A1), vB1_ = vbcast(B1);
                            lrb_v vA2_ = vbcast(A2), vB2_ = vbcast(B2);

                            lrb_v ev0 = vfma(vA0_, lane_xoff,
                                       vfma(vB0_, lane_yoff, vbcast(e0)));
                            lrb_v ev1 = vfma(vA1_, lane_xoff,
                                       vfma(vB1_, lane_yoff, vbcast(e1)));
                            lrb_v ev2 = vfma(vA2_, lane_xoff,
                                       vfma(vB2_, lane_yoff, vbcast(e2)));

                            lrb_k coverage = 0xFFFF;  /* fully covered */

                            lrb_v v_inv = vbcast(inv_a);
                            lrb_v bw0 = vmul(ev1, v_inv);
                            lrb_v bw1 = vmul(ev2, v_inv);
                            lrb_v bw2 = vsub(vsub(vbcast(1.0f), bw0), bw1);

                            lrb_v depth = vfma(bw0, vbcast(z0),
                                          vfma(bw1, vbcast(z1),
                                               vmul(bw2, vbcast(z2))));
                            lrb_v cr = vfma(bw0, vbcast(r0),
                                       vfma(bw1, vbcast(r1),
                                            vmul(bw2, vbcast(r2))));
                            lrb_v cg = vfma(bw0, vbcast(g0),
                                       vfma(bw1, vbcast(g1),
                                            vmul(bw2, vbcast(g2))));
                            lrb_v cb = vfma(bw0, vbcast(b0),
                                       vfma(bw1, vbcast(b1),
                                            vmul(bw2, vbcast(b2))));

                            lrb_v db = vloadu(&tb->depth[swiz]);
                            lrb_k write = kand(coverage, klt(depth, db));

                            if (write != 0) {
                                vstore_k(&tb->depth[swiz],   depth, write);
                                vstore_k(&tb->color_r[swiz], cr,    write);
                                vstore_k(&tb->color_g[swiz], cg,    write);
                                vstore_k(&tb->color_b[swiz], cb,    write);
                                vstore_k(&tb->color_a[swiz], vbcast(1.0f), write);
                            }
                        } else {
                            /* Mixed — test edges per 4x4 block (lanes = pixels) */
                            rasterize_tri_fine(tb, swiz,
                                A0, B0, C0, A1, B1, C1, A2, B2, C2,
                                inv_a,
                                z0, z1, z2,
                                r0, r1, r2, g0, g1, g2, b0, b1, b2,
                                e0, e1, e2, is_ccw);
                        }

                        /* Prefetch next fine block row */
                        if (fy + LRB_FINE_BLOCK < LRB_COARSE_BLOCK) {
                            int next_swiz = swizzle(lx, ly + LRB_FINE_BLOCK);
                            prefetch_l2(&tb->depth[next_swiz]);
                        }

                        /* Incremental scalar step in x */
                        e0 += dA0x; e1 += dA1x; e2 += dA2x;
                    }
                    /* Incremental scalar step in y */
                    e0_row += dB0y; e1_row += dB1y; e2_row += dB2y;
                }
            }
        }
    }
}

/* ============================================================
 * 14. Core-Affine Thread Scheduling
 * ============================================================ */

typedef struct {
    int tile_x, tile_y;
    int *pkt_list;
    int num_pkts;
} lrb_tile_work;

typedef struct {
    lrb_tile_buf *tile_bufs;
    lrb_tri_packet *packets;
    int num_packets;
    lrb_tile_work *tiles;
    int num_tiles;
    atomic_int next_tile;
    int fb_w, fb_h;
    float *fb_color;
} lrb_render_ctx;

static void pin_to_core(int thread_id) {
    int core = thread_id / LRB_THREADS_PER_CORE;
    cpu_set_t cpuset;
    CPU_ZERO(&cpuset);
    CPU_SET(core, &cpuset);
    pthread_setaffinity_np(pthread_self(), sizeof(cpu_set_t), &cpuset);
}

static void resolve_tile(lrb_tile_buf *tb, int tile_x, int tile_y,
                         int fb_w, int fb_h, float *fb_color)
{
    int tx0 = tile_x * LRB_TILE_SIZE;
    int ty0 = tile_y * LRB_TILE_SIZE;

    for (int by = 0; by < TILE_BY && ty0 + by * 4 < fb_h; by++) {
        for (int bx = 0; bx < TILE_BX && tx0 + bx * 4 < fb_w; bx++) {
            int swiz_base = (by * TILE_BX + bx) * 16;

            lrb_v r = vload(&tb->color_r[swiz_base]);
            lrb_v g = vload(&tb->color_g[swiz_base]);
            lrb_v b = vload(&tb->color_b[swiz_base]);
            lrb_v a = vload(&tb->color_a[swiz_base]);

            alignas(64) float ra[16], ga[16], ba[16], aa[16];
            vstore(ra, r); vstore(ga, g);
            vstore(ba, b); vstore(aa, a);

            for (int ly = 0; ly < 4; ly++) {
                int sy = ty0 + by * 4 + ly;
                if (sy >= fb_h) break;
                for (int lx = 0; lx < 4; lx++) {
                    int sx = tx0 + bx * 4 + lx;
                    if (sx >= fb_w) break;
                    int lane = ly * 4 + lx;
                    int fi = (sy * fb_w + sx) * 4;
                    fb_color[fi + 0] = ra[lane];
                    fb_color[fi + 1] = ga[lane];
                    fb_color[fi + 2] = ba[lane];
                    fb_color[fi + 3] = aa[lane];
                }
            }
        }
    }
}

static atomic_int g_tid_counter = 0;

static void *render_worker(void *arg) {
    lrb_render_ctx *ctx = (lrb_render_ctx *)arg;

    int tid = atomic_fetch_add(&g_tid_counter, 1);
    pin_to_core(tid);

    lrb_tile_buf *tb = &ctx->tile_bufs[tid];

    for (;;) {
        int tile_idx = atomic_fetch_add(&ctx->next_tile, 1);
        if (tile_idx >= ctx->num_tiles) break;

        lrb_tile_work *tw = &ctx->tiles[tile_idx];

        /* Clear tile buffer */
        lrb_v cd = vbcast(1e30f), c0 = vbcast(0.0f), c1 = vbcast(1.0f);
        for (int i = 0; i < TILE_PIX; i += 16) {
            vstore(&tb->depth[i], cd);
            vstore(&tb->color_r[i], c0);
            vstore(&tb->color_g[i], c0);
            vstore(&tb->color_b[i], c0);
            vstore(&tb->color_a[i], c1);
        }

        /* Rasterize all packets binned to this tile */
        for (int p = 0; p < tw->num_pkts; p++) {
            lrb_tri_packet *pkt = &ctx->packets[tw->pkt_list[p]];
            rasterize_packet_tile(tb, pkt, tw->tile_x, tw->tile_y,
                                  ctx->fb_w, ctx->fb_h);
        }

        resolve_tile(tb, tw->tile_x, tw->tile_y,
                     ctx->fb_w, ctx->fb_h, ctx->fb_color);
    }
    return NULL;
}

/* ============================================================
 * 15. Full Pipeline
 * ============================================================ */

typedef struct {
    int fb_w, fb_h;
    float *fb_color;
    float *fb_depth;
    int tiles_x, tiles_y;
} lrb_device;

static void lrb_draw(lrb_device *dev, lrb_vertex *verts, int nv,
                     lrb_tri *tris, int nt, const lrb_shader_ctx *sctx)
{
    printf("[Draw] %d verts, %d tris\n", nv, nt);
    init_lane_offsets();

    /* Stage 1: Vertex shader */
    int nblocks = (nv + 15) / 16;
    lrb_vertex_soa *soa_in  = aligned_alloc(LRB_ALIGN, nblocks * sizeof(lrb_vertex_soa));
    lrb_vertex_soa *soa_out = aligned_alloc(LRB_ALIGN, nblocks * sizeof(lrb_vertex_soa));
    vec4 *xpos = calloc(nv, sizeof(vec4));
    vec4 *xcol = calloc(nv, sizeof(vec4));

    for (int b = 0; b < nblocks; b++) {
        int off = b * 16;
        int cnt = (off + 16 <= nv) ? 16 : (nv - off);
        aos_to_soa(&verts[off], &soa_in[b], cnt);
        vertex_shader(&soa_in[b], &soa_out[b], sctx);
        for (int i = 0; i < cnt; i++) {
            xpos[off+i].x = soa_out[b].x[i];  xpos[off+i].y = soa_out[b].y[i];
            xpos[off+i].z = soa_out[b].z[i];  xpos[off+i].w = soa_out[b].w[i];
            xcol[off+i].x = soa_out[b].r[i];  xcol[off+i].y = soa_out[b].g[i];
            xcol[off+i].z = soa_out[b].b[i];  xcol[off+i].w = soa_out[b].a[i];
        }
    }

    /* Stage 2: Build packets */
    int num_packets = (nt + LRB_PACKET_TRIS - 1) / LRB_PACKET_TRIS;
    lrb_tri_packet *packets = aligned_alloc(LRB_ALIGN,
                                            num_packets * sizeof(lrb_tri_packet));

    for (int p = 0; p < num_packets; p++) {
        lrb_setup_input inputs[LRB_PACKET_TRIS];
        int base = p * LRB_PACKET_TRIS;
        int cnt = (base + LRB_PACKET_TRIS <= nt) ? LRB_PACKET_TRIS : (nt - base);
        for (int i = 0; i < cnt; i++) {
            int t = base + i;
            int i0 = tris[t].v[0], i1 = tris[t].v[1], i2 = tris[t].v[2];
            inputs[i].x0 = xpos[i0].x; inputs[i].y0 = xpos[i0].y; inputs[i].z0 = xpos[i0].z;
            inputs[i].x1 = xpos[i1].x; inputs[i].y1 = xpos[i1].y; inputs[i].z1 = xpos[i1].z;
            inputs[i].x2 = xpos[i2].x; inputs[i].y2 = xpos[i2].y; inputs[i].z2 = xpos[i2].z;
            inputs[i].c0 = xcol[i0]; inputs[i].c1 = xcol[i1]; inputs[i].c2 = xcol[i2];
        }
        build_packet(inputs, cnt, &packets[p]);
    }

    /* Stage 3: Tile binning */
    int total_tiles = dev->tiles_x * dev->tiles_y;
    lrb_tile_work *tiles = calloc(total_tiles, sizeof(lrb_tile_work));
    for (int i = 0; i < total_tiles; i++) {
        tiles[i].tile_x = i % dev->tiles_x;
        tiles[i].tile_y = i / dev->tiles_x;
        tiles[i].pkt_list = calloc(num_packets, sizeof(int));
        tiles[i].num_pkts = 0;
    }

    for (int p = 0; p < num_packets; p++) {
        lrb_tri_packet *pkt = &packets[p];
        alignas(64) float minx_a[16], maxx_a[16], miny_a[16], maxy_a[16];
        vstore(minx_a, vload(pkt->minx)); vstore(maxx_a, vload(pkt->maxx));
        vstore(miny_a, vload(pkt->miny)); vstore(maxy_a, vload(pkt->maxy));

        float pb_minx = 1e30f, pb_maxx = -1e30f;
        float pb_miny = 1e30f, pb_maxy = -1e30f;
        for (int i = 0; i < pkt->count; i++) {
            if (pkt->valid_mask & (1u << i)) {
                pb_minx = fminf(pb_minx, minx_a[i]);
                pb_maxx = fmaxf(pb_maxx, maxx_a[i]);
                pb_miny = fminf(pb_miny, miny_a[i]);
                pb_maxy = fmaxf(pb_maxy, maxy_a[i]);
            }
        }

        int tx0 = (int)(pb_minx / LRB_TILE_SIZE);
        int tx1 = (int)(pb_maxx / LRB_TILE_SIZE);
        int ty0 = (int)(pb_miny / LRB_TILE_SIZE);
        int ty1 = (int)(pb_maxy / LRB_TILE_SIZE);
        if (tx0 < 0) tx0 = 0;
        if (ty0 < 0) ty0 = 0;
        if (tx1 >= dev->tiles_x) tx1 = dev->tiles_x - 1;
        if (ty1 >= dev->tiles_y) ty1 = dev->tiles_y - 1;

        for (int ty = ty0; ty <= ty1; ty++)
            for (int tx = tx0; tx <= tx1; tx++) {
                int ti = ty * dev->tiles_x + tx;
                tiles[ti].pkt_list[tiles[ti].num_pkts++] = p;
            }
    }

    /* Stage 4: Parallel rasterization */
    lrb_render_ctx ctx;
    ctx.tile_bufs = aligned_alloc(LRB_ALIGN,
                                  LRB_TOTAL_THREADS * sizeof(lrb_tile_buf));
    ctx.packets = packets;
    ctx.num_packets = num_packets;
    ctx.tiles = tiles;
    ctx.num_tiles = total_tiles;
    atomic_init(&ctx.next_tile, 0);
    atomic_init(&g_tid_counter, 0);
    ctx.fb_w = dev->fb_w;
    ctx.fb_h = dev->fb_h;
    ctx.fb_color = dev->fb_color;

    pthread_t threads[LRB_TOTAL_THREADS];
    printf("[Pipeline] %d threads → %d cores, %d tiles, %d packets\n",
           LRB_TOTAL_THREADS, LRB_CORES, total_tiles, num_packets);
    printf("[Pipeline] Coarse: lanes=tris (16 tri packet), "
           "Fine: lanes=pixels (4x4 block)\n");
    printf("[Pipeline] tile=%dx%d (80KB L2), incremental scalar edge stepping\n",
           LRB_TILE_SIZE, LRB_TILE_SIZE);

    for (int i = 0; i < LRB_TOTAL_THREADS; i++)
        pthread_create(&threads[i], NULL, render_worker, &ctx);
    for (int i = 0; i < LRB_TOTAL_THREADS; i++)
        pthread_join(threads[i], NULL);

    /* Cleanup */
    for (int i = 0; i < total_tiles; i++) free(tiles[i].pkt_list);
    free(tiles); free(packets); free(xpos); free(xcol);
    free(soa_in); free(soa_out); free(ctx.tile_bufs);
}

/* ============================================================
 * 16. Output & Main
 * ============================================================ */

static void present_ppm(lrb_device *dev, const char *fn) {
    FILE *f = fopen(fn, "wb");
    if (!f) return;
    fprintf(f, "P6\n%d %d\n255\n", dev->fb_w, dev->fb_h);
    for (int i = 0; i < dev->fb_w * dev->fb_h; i++) {
        unsigned char r = (unsigned char)(fminf(dev->fb_color[i*4],   1.0f) * 255);
        unsigned char g = (unsigned char)(fminf(dev->fb_color[i*4+1], 1.0f) * 255);
        unsigned char b = (unsigned char)(fminf(dev->fb_color[i*4+2], 1.0f) * 255);
        fputc(r, f); fputc(g, f); fputc(b, f);
    }
    fclose(f);
    printf("[Present] %s (%dx%d)\n", fn, dev->fb_w, dev->fb_h);
}

int main(void) {
    printf("=== Larrabee KNF Renderer v3 — Correct Lane Assignment ===\n\n");
    printf("Coarse traversal: lanes = 16 triangles (packet)\n");
    printf("Fine rasterization: lanes = 16 pixels (4x4 block)\n");
    printf("Incremental stepping: per-triangle scalar + vectorized per-pixel\n\n");

    lrb_device dev;
    dev.fb_w = 512; dev.fb_h = 512;
    dev.tiles_x = (512 + LRB_TILE_SIZE - 1) / LRB_TILE_SIZE;
    dev.tiles_y = (512 + LRB_TILE_SIZE - 1) / LRB_TILE_SIZE;
    posix_memalign((void **)&dev.fb_color, LRB_ALIGN, 512*512*4*sizeof(float));
    posix_memalign((void **)&dev.fb_depth, LRB_ALIGN, 512*512*sizeof(float));

    lrb_shader_ctx sctx;
    sctx.fb_w = 512; sctx.fb_h = 512;
    lrb_mat4 model = mat4_mul(mat4_trans(0, 0, -5), mat4_roty(0.7f));
    sctx.mvp = mat4_mul(mat4_persp(60, 1, 0.1f, 100), model);

    lrb_vertex verts[8];
    float p[8][3] = {{-1,-1,-1},{1,-1,-1},{1,1,-1},{-1,1,-1},
                     {-1,-1,1},{1,-1,1},{1,1,1},{-1,1,1}};
    vec4 col[8] = {{1,0,0,1},{0,1,0,1},{0,0,1,1},{1,1,0,1},
                   {1,0,1,1},{0,1,1,1},{1,1,1,1},{0.5f,0.5f,0.5f,1}};
    for (int i = 0; i < 8; i++) {
        verts[i].pos = (vec4){p[i][0], p[i][1], p[i][2], 1};
        verts[i].color = col[i];
        verts[i].uv = (vec2){0, 0};
    }
    int tris_arr[12][3] = {{0,1,2},{0,2,3},{5,4,7},{5,7,6},{4,0,3},{4,3,7},
                           {1,5,6},{1,6,2},{3,2,6},{3,6,7},{4,5,1},{4,1,0}};
    lrb_tri tris[12];
    for (int i = 0; i < 12; i++)
        tris[i] = (lrb_tri){{tris_arr[i][0], tris_arr[i][1], tris_arr[i][2]}};

    lrb_draw(&dev, verts, 8, tris, 12, &sctx);
    present_ppm(&dev, "lrb_knf_v3_output.ppm");

    free(dev.fb_color);
    free(dev.fb_depth);
    return 0;
}