Path B pivot + Phase 0-3 closed with first baseline numbers

This is a from-scratch initial commit on a fresh .git. The original
scaffold commit (7510b56) and the earlier session's working-tree
docs were lost in a 2026-05-18 10:25 working-tree wipe; the corrupted
.git is preserved at .git-broken-2026-05-18/ (gitignored) for
forensic inspection.

Scope re-anchored from Path A (custom VPU firmware on VC7 scalar
cores; blocked by BCM2712 silicon-RoT mask-ROM signature check)
to Path B (QPU compute kernels via Mesa v3d / Vulkan compute or
direct DRM, on stock signed Pi 5 / CM5). See README.md and
docs/phase0.md for the substrate audit that closed Path A.

Phases closed:
  Phase 0 — substrate audit; Path A blocked, Path B open;
            codec-back-end-fits-QPU finding (docs/phase0.md)
  Phase 1 — first kernel locked (VP9 / AV1 8x8 inverse DCT) with
            publish-before-measure R = M2/M3 decision rules
            (docs/phase1.md)
  Phase 2 — reference impls mapped; FFmpeg n7.1.3 source vendored
            under external/ffmpeg-snapshot/ (PROVENANCE.md pins
            commit f46e514 + per-file SHA-256s) (docs/phase2.md)
  Phase 3 — real baseline measurements on hertz (docs/phase3.md):
              M1 bit-exact            100.0000 % (10000/10000)
              M3 NEON IDCT8 single    8.171 Mblock/s (122.4 ns/block)
              M5a empty Vulkan submit 22.66 us
              M5b 1-WG noop dispatch  55.60 us
              M5 delta                32.95 us/dispatch
            => per-dispatch overhead is ~455x per-NEON-block cost;
               Phase 4 must batch at frame level or close to it.

Build harness in place: CMakeLists.txt + tests/{bench_neon_idct.c,
vp9_idct8_ref.c, bench_vulkan_dispatch.c, shaders/noop.comp} +
external/ffmpeg-snapshot/config.h shim (7 defines + EXTERN_ASM).
Builds clean on Debian Trixie aarch64 with cmake 3.31, ninja 1.12,
libvulkan-dev 1.4.309, glslang-tools 15.1.0. Vendored FFmpeg .S
assembles via the config.h shim.

Next: Phase 4 (plan first QPU IDCT kernel under the M5 batching
constraint) -> Phase 5 second-model review -> Phase 6 implement.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
This commit is contained in:
2026-05-18 11:30:12 +00:00
commit dcbbc77038
22 changed files with 9030 additions and 0 deletions
View File
+248
View File
@@ -0,0 +1,248 @@
/*
* Phase 3 — NEON baseline microbench for VP9 8×8 DCT_DCT IDCT add.
*
* Reports two numbers:
* M1 (correctness): bit-exact match rate, our C reference vs
* FFmpeg's NEON, across N random blocks.
* M3 (throughput): NEON sustained MblockS on this host.
*
* Both are gating measurements for Phase 1 (see docs/phase1.md).
* NO QPU work happens here — that's later phases.
*
* Build: see CMakeLists.txt at project root.
* Run: ./bench_neon_idct [--blocks N] [--iters K] [--seed S]
*
* License: BSD-2-Clause (daedalus-fourier), but this binary
* statically links the LGPL-2.1+ FFmpeg NEON snapshot
* — distribute the binary under LGPL-2.1+ in that case.
*/
#define _POSIX_C_SOURCE 200809L
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <stddef.h>
#include <time.h>
#include <getopt.h>
/* Our C reference (tests/vp9_idct8_ref.c). */
extern void daedalus_vp9_idct_idct_8x8_add_ref(
uint8_t *dst, ptrdiff_t stride, int16_t *block, int eob);
/* FFmpeg NEON entry point (vendored vp9itxfm_neon.S). */
extern void ff_vp9_idct_idct_8x8_add_neon(
uint8_t *dst, ptrdiff_t stride, int16_t *block, int eob);
/* ---- Random-block generation ----------------------------------- */
/* xorshift64 — deterministic per seed, fast enough not to dominate
* the measurement. */
static uint64_t xs64_state;
static inline uint64_t xs64(void)
{
uint64_t x = xs64_state;
x ^= x << 13; x ^= x >> 7; x ^= x << 17;
return xs64_state = x;
}
/* Random VP9-plausible coefficient block: most coefficients zero,
* a handful of nonzero ones in low-frequency positions. Bias chosen
* so eob is typically in [4, 32], hitting the general (non-DC) path.
* For Phase 3 baseline this isn't load-balanced against a real
* bitstream distribution — Phase 7 may revisit. */
static int gen_block(int16_t block[64])
{
memset(block, 0, 64 * sizeof(*block));
int eob = 0;
int n_nonzero = 1 + (int)(xs64() % 16);
for (int i = 0; i < n_nonzero; i++) {
/* Bias toward low-freq positions via xs64() % (xs64() % 64 + 1). */
int pos = (int)(xs64() % 64);
/* Coefficient range: signed 12-bit (typical dequant output). */
int16_t coef = (int16_t)((int)(xs64() % 8192) - 4096);
block[pos] = coef;
if (pos + 1 > eob) eob = pos + 1;
}
if (eob == 0) eob = 1;
return eob;
}
static void gen_pred(uint8_t pred[64])
{
for (int i = 0; i < 64; i++)
pred[i] = (uint8_t)(xs64() & 0xff);
}
/* ---- Wall-clock timing (CLOCK_MONOTONIC_RAW) ------------------- */
static double now_seconds(void)
{
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC_RAW, &ts);
return ts.tv_sec + ts.tv_nsec * 1e-9;
}
/* ---- Phase 1 M1: bit-exact gate -------------------------------- */
static int correctness_check(uint64_t seed, int n_blocks)
{
xs64_state = seed ? seed : 0xdeadbeefcafebabeULL;
int mismatches = 0;
int dc_only_seen = 0;
int16_t block_a[64], block_b[64];
uint8_t pred[64];
uint8_t dst_a[64], dst_b[64];
for (int i = 0; i < n_blocks; i++) {
int eob = gen_block(block_a);
memcpy(block_b, block_a, sizeof(block_a));
gen_pred(pred);
memcpy(dst_a, pred, 64);
memcpy(dst_b, pred, 64);
daedalus_vp9_idct_idct_8x8_add_ref(dst_a, 8, block_a, eob);
ff_vp9_idct_idct_8x8_add_neon(dst_b, 8, block_b, eob);
if (memcmp(dst_a, dst_b, 64) != 0) {
if (mismatches < 4) {
fprintf(stderr, "MISMATCH block %d eob=%d:\n", i, eob);
for (int r = 0; r < 8; r++) {
fprintf(stderr, " row %d ref ", r);
for (int c = 0; c < 8; c++) fprintf(stderr, "%3u ", dst_a[r * 8 + c]);
fprintf(stderr, " neon ");
for (int c = 0; c < 8; c++) fprintf(stderr, "%3u ", dst_b[r * 8 + c]);
fprintf(stderr, "\n");
}
}
mismatches++;
}
if (eob == 1) dc_only_seen++;
}
printf("M1 correctness: %d / %d blocks bit-exact match (%.4f%%)\n",
n_blocks - mismatches, n_blocks,
100.0 * (n_blocks - mismatches) / n_blocks);
printf(" dc-only path frequency: %d / %d (%.2f%%)\n",
dc_only_seen, n_blocks, 100.0 * dc_only_seen / n_blocks);
return mismatches;
}
/* ---- Phase 1 M3: NEON throughput ------------------------------- */
static void throughput_neon(uint64_t seed, int n_blocks, int iters)
{
xs64_state = seed ? seed : 0xfeedfacecafebeefULL;
/* Pre-generate all blocks + preds so generation cost is excluded
* from the timed region. Each block is consumed once per iteration
* (NEON path zeroes the block, so we restore from the master). */
int16_t *blocks_master = malloc(n_blocks * 64 * sizeof(int16_t));
int16_t *blocks_work = malloc(n_blocks * 64 * sizeof(int16_t));
uint8_t *preds = malloc(n_blocks * 64);
uint8_t *dsts = malloc(n_blocks * 64);
int *eobs = malloc(n_blocks * sizeof(int));
if (!blocks_master || !blocks_work || !preds || !dsts || !eobs) {
fprintf(stderr, "alloc failed\n");
exit(1);
}
for (int i = 0; i < n_blocks; i++) {
eobs[i] = gen_block(blocks_master + i * 64);
gen_pred(preds + i * 64);
}
/* Warm-up. */
memcpy(blocks_work, blocks_master, n_blocks * 64 * sizeof(int16_t));
memcpy(dsts, preds, n_blocks * 64);
for (int i = 0; i < n_blocks; i++)
ff_vp9_idct_idct_8x8_add_neon(dsts + i * 64, 8,
blocks_work + i * 64, eobs[i]);
/* Timed region. */
double t0 = now_seconds();
for (int it = 0; it < iters; it++) {
memcpy(blocks_work, blocks_master, n_blocks * 64 * sizeof(int16_t));
memcpy(dsts, preds, n_blocks * 64);
for (int i = 0; i < n_blocks; i++)
ff_vp9_idct_idct_8x8_add_neon(dsts + i * 64, 8,
blocks_work + i * 64, eobs[i]);
}
double t1 = now_seconds();
/* memcpy cost-only run, to subtract setup overhead. */
double s0 = now_seconds();
for (int it = 0; it < iters; it++) {
memcpy(blocks_work, blocks_master, n_blocks * 64 * sizeof(int16_t));
memcpy(dsts, preds, n_blocks * 64);
}
double s1 = now_seconds();
double total_seconds = (t1 - t0) - (s1 - s0);
double total_blocks = (double) n_blocks * iters;
double mblocks_s = total_blocks / total_seconds / 1e6;
printf("M3 NEON throughput:\n");
printf(" blocks=%d iters=%d total=%.0f\n", n_blocks, iters, total_blocks);
printf(" elapsed (kernel)=%.6f s (setup-subtracted)\n", total_seconds);
printf(" elapsed (setup) =%.6f s\n", s1 - s0);
printf(" throughput = %.3f Mblock/s\n", mblocks_s);
printf(" per-block = %.1f ns\n", total_seconds / total_blocks * 1e9);
/* Equivalent at 1920x1080: 32 400 blocks/frame -> FPS. */
printf(" equiv 1080p = %.1f FPS (32400 blocks/frame)\n",
mblocks_s * 1e6 / 32400.0);
free(blocks_master); free(blocks_work); free(preds);
free(dsts); free(eobs);
}
/* ---- CLI ------------------------------------------------------- */
static void usage(const char *p)
{
fprintf(stderr,
"Usage: %s [--blocks N] [--iters K] [--seed S] [--no-correctness]\n"
"Defaults: N=1000000, K=10, S=0 (uses fixed default).\n", p);
}
int main(int argc, char **argv)
{
int n_blocks = 1000000;
int iters = 10;
uint64_t seed = 0;
int do_correctness = 1;
static struct option opts[] = {
{"blocks", required_argument, 0, 'b'},
{"iters", required_argument, 0, 'i'},
{"seed", required_argument, 0, 's'},
{"no-correctness", no_argument, 0, 'C'},
{"help", no_argument, 0, 'h'},
{0,0,0,0}
};
for (int c; (c = getopt_long(argc, argv, "b:i:s:Ch", opts, 0)) != -1;) {
switch (c) {
case 'b': n_blocks = atoi(optarg); break;
case 'i': iters = atoi(optarg); break;
case 's': seed = strtoull(optarg, 0, 0); break;
case 'C': do_correctness = 0; break;
case 'h': usage(argv[0]); return 0;
default: usage(argv[0]); return 2;
}
}
if (do_correctness) {
printf("=== M1: bit-exact correctness (10000 random blocks) ===\n");
int miss = correctness_check(seed, 10000);
if (miss != 0) {
fprintf(stderr, "REFUSING to measure throughput on a broken kernel.\n");
return 1;
}
printf("\n");
}
printf("=== M3: NEON throughput ===\n");
throughput_neon(seed, n_blocks, iters);
return 0;
}
+279
View File
@@ -0,0 +1,279 @@
/*
* Phase 3 — Vulkan compute dispatch-overhead microbench (M5).
*
* Measures the per-dispatch wall-clock floor on V3D 7.1 via Mesa
* v3dv: vkQueueSubmit + vkQueueWaitIdle round-trip cost for a
* noop compute shader. Establishes the floor below which kernel
* batching is mandatory.
*
* Two measurements:
* M5a: empty command-buffer submit (no dispatch at all)
* M5b: 1-workgroup dispatch of an empty shader
*
* The delta M5b - M5a isolates the per-vkCmdDispatch cost from
* the per-vkQueueSubmit cost.
*
* Build: cmake -DDAEDALUS_BUILD_VULKAN=ON ..
* Run: ./bench_vulkan_dispatch [--iters N] [--spv PATH]
*
* License: BSD-2-Clause (daedalus-fourier).
*/
#define _POSIX_C_SOURCE 200809L
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <time.h>
#include <getopt.h>
#include <vulkan/vulkan.h>
#define CHK(call) do { VkResult r__ = (call); if (r__ != VK_SUCCESS) { \
fprintf(stderr, "vulkan error %d at %s:%d (%s)\n", r__, __FILE__, __LINE__, #call); \
exit(1); } } while (0)
static double now_seconds(void)
{
struct timespec ts;
clock_gettime(CLOCK_MONOTONIC_RAW, &ts);
return ts.tv_sec + ts.tv_nsec * 1e-9;
}
static uint32_t *read_spv(const char *path, size_t *out_size)
{
FILE *f = fopen(path, "rb");
if (!f) { perror(path); exit(1); }
fseek(f, 0, SEEK_END);
long sz = ftell(f);
fseek(f, 0, SEEK_SET);
if (sz <= 0 || (sz & 3)) {
fprintf(stderr, "%s: bad SPIR-V size %ld\n", path, sz);
exit(1);
}
uint32_t *buf = malloc(sz);
if (!buf || fread(buf, 1, sz, f) != (size_t)sz) {
perror("read"); exit(1);
}
fclose(f);
*out_size = sz;
return buf;
}
int main(int argc, char **argv)
{
int iters = 100000;
const char *spv_path = "noop.spv";
static struct option opts[] = {
{"iters", required_argument, 0, 'i'},
{"spv", required_argument, 0, 's'},
{"help", no_argument, 0, 'h'},
{0,0,0,0}
};
for (int c; (c = getopt_long(argc, argv, "i:s:h", opts, 0)) != -1;) {
switch (c) {
case 'i': iters = atoi(optarg); break;
case 's': spv_path = optarg; break;
case 'h':
fprintf(stderr,
"Usage: %s [--iters N] [--spv noop.spv]\n", argv[0]);
return 0;
default:
return 2;
}
}
/* ---- Instance ---- */
VkApplicationInfo app = {
.sType = VK_STRUCTURE_TYPE_APPLICATION_INFO,
.pApplicationName = "daedalus-fourier-bench",
.apiVersion = VK_API_VERSION_1_3,
};
VkInstanceCreateInfo ici = {
.sType = VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO,
.pApplicationInfo = &app,
};
VkInstance instance;
CHK(vkCreateInstance(&ici, NULL, &instance));
/* ---- Pick V3D physical device (skip llvmpipe) ---- */
uint32_t pd_count = 0;
CHK(vkEnumeratePhysicalDevices(instance, &pd_count, NULL));
VkPhysicalDevice *pds = malloc(pd_count * sizeof(*pds));
CHK(vkEnumeratePhysicalDevices(instance, &pd_count, pds));
VkPhysicalDevice phys = VK_NULL_HANDLE;
VkPhysicalDeviceProperties props = {0};
for (uint32_t i = 0; i < pd_count; i++) {
vkGetPhysicalDeviceProperties(pds[i], &props);
printf("device %u: %s (api %u.%u.%u, vendor 0x%04x)\n",
i, props.deviceName,
VK_VERSION_MAJOR(props.apiVersion),
VK_VERSION_MINOR(props.apiVersion),
VK_VERSION_PATCH(props.apiVersion),
props.vendorID);
if (strstr(props.deviceName, "V3D") != NULL) {
phys = pds[i];
}
}
if (phys == VK_NULL_HANDLE) {
fprintf(stderr, "no V3D device found; bailing.\n");
return 1;
}
vkGetPhysicalDeviceProperties(phys, &props);
printf("selected: %s\n", props.deviceName);
free(pds);
/* ---- Compute queue family ---- */
uint32_t qfc = 0;
vkGetPhysicalDeviceQueueFamilyProperties(phys, &qfc, NULL);
VkQueueFamilyProperties *qfp = malloc(qfc * sizeof(*qfp));
vkGetPhysicalDeviceQueueFamilyProperties(phys, &qfc, qfp);
uint32_t qfi = (uint32_t) -1;
for (uint32_t i = 0; i < qfc; i++) {
if (qfp[i].queueFlags & VK_QUEUE_COMPUTE_BIT) {
qfi = i; break;
}
}
if (qfi == (uint32_t) -1) {
fprintf(stderr, "no compute queue family\n");
return 1;
}
free(qfp);
/* ---- Logical device ---- */
float qprio = 1.0f;
VkDeviceQueueCreateInfo dqci = {
.sType = VK_STRUCTURE_TYPE_DEVICE_QUEUE_CREATE_INFO,
.queueFamilyIndex = qfi,
.queueCount = 1,
.pQueuePriorities = &qprio,
};
VkDeviceCreateInfo dci = {
.sType = VK_STRUCTURE_TYPE_DEVICE_CREATE_INFO,
.queueCreateInfoCount = 1,
.pQueueCreateInfos = &dqci,
};
VkDevice dev;
CHK(vkCreateDevice(phys, &dci, NULL, &dev));
VkQueue queue;
vkGetDeviceQueue(dev, qfi, 0, &queue);
/* ---- Command pool + buffers ---- */
VkCommandPoolCreateInfo cpci = {
.sType = VK_STRUCTURE_TYPE_COMMAND_POOL_CREATE_INFO,
.flags = VK_COMMAND_POOL_CREATE_RESET_COMMAND_BUFFER_BIT,
.queueFamilyIndex = qfi,
};
VkCommandPool pool;
CHK(vkCreateCommandPool(dev, &cpci, NULL, &pool));
VkCommandBuffer cb_empty, cb_dispatch;
VkCommandBufferAllocateInfo cbai = {
.sType = VK_STRUCTURE_TYPE_COMMAND_BUFFER_ALLOCATE_INFO,
.commandPool = pool,
.level = VK_COMMAND_BUFFER_LEVEL_PRIMARY,
.commandBufferCount = 1,
};
CHK(vkAllocateCommandBuffers(dev, &cbai, &cb_empty));
CHK(vkAllocateCommandBuffers(dev, &cbai, &cb_dispatch));
/* ---- Pipeline layout (empty: no descriptors, no push constants) ---- */
VkPipelineLayoutCreateInfo plci = {
.sType = VK_STRUCTURE_TYPE_PIPELINE_LAYOUT_CREATE_INFO,
};
VkPipelineLayout playout;
CHK(vkCreatePipelineLayout(dev, &plci, NULL, &playout));
/* ---- Compute pipeline from noop SPIR-V ---- */
size_t spv_size = 0;
uint32_t *spv = read_spv(spv_path, &spv_size);
VkShaderModuleCreateInfo smci = {
.sType = VK_STRUCTURE_TYPE_SHADER_MODULE_CREATE_INFO,
.codeSize = spv_size,
.pCode = spv,
};
VkShaderModule shader;
CHK(vkCreateShaderModule(dev, &smci, NULL, &shader));
free(spv);
VkComputePipelineCreateInfo cpci2 = {
.sType = VK_STRUCTURE_TYPE_COMPUTE_PIPELINE_CREATE_INFO,
.stage = {
.sType = VK_STRUCTURE_TYPE_PIPELINE_SHADER_STAGE_CREATE_INFO,
.stage = VK_SHADER_STAGE_COMPUTE_BIT,
.module = shader,
.pName = "main",
},
.layout = playout,
};
VkPipeline pipe;
CHK(vkCreateComputePipelines(dev, VK_NULL_HANDLE, 1, &cpci2, NULL, &pipe));
/* ---- Record both command buffers once, reuse for every iteration ---- */
VkCommandBufferBeginInfo cbbi = {
.sType = VK_STRUCTURE_TYPE_COMMAND_BUFFER_BEGIN_INFO,
};
CHK(vkBeginCommandBuffer(cb_empty, &cbbi));
CHK(vkEndCommandBuffer(cb_empty));
CHK(vkBeginCommandBuffer(cb_dispatch, &cbbi));
vkCmdBindPipeline(cb_dispatch, VK_PIPELINE_BIND_POINT_COMPUTE, pipe);
vkCmdDispatch(cb_dispatch, 1, 1, 1);
CHK(vkEndCommandBuffer(cb_dispatch));
VkSubmitInfo si_empty = {
.sType = VK_STRUCTURE_TYPE_SUBMIT_INFO,
.commandBufferCount = 1, .pCommandBuffers = &cb_empty,
};
VkSubmitInfo si_disp = {
.sType = VK_STRUCTURE_TYPE_SUBMIT_INFO,
.commandBufferCount = 1, .pCommandBuffers = &cb_dispatch,
};
/* ---- Warm-up ---- */
for (int i = 0; i < 100; i++) {
CHK(vkQueueSubmit(queue, 1, &si_disp, VK_NULL_HANDLE));
CHK(vkQueueWaitIdle(queue));
}
/* ---- M5a: empty CB submit+wait ---- */
double t0 = now_seconds();
for (int i = 0; i < iters; i++) {
CHK(vkQueueSubmit(queue, 1, &si_empty, VK_NULL_HANDLE));
CHK(vkQueueWaitIdle(queue));
}
double t1 = now_seconds();
double m5a_per = (t1 - t0) / iters * 1e6; /* µs */
/* ---- M5b: 1-WG noop dispatch submit+wait ---- */
double t2 = now_seconds();
for (int i = 0; i < iters; i++) {
CHK(vkQueueSubmit(queue, 1, &si_disp, VK_NULL_HANDLE));
CHK(vkQueueWaitIdle(queue));
}
double t3 = now_seconds();
double m5b_per = (t3 - t2) / iters * 1e6; /* µs */
printf("\n=== M5: Vulkan compute dispatch overhead ===\n");
printf(" iters per measurement: %d\n", iters);
printf(" M5a empty CB submit+wait: %.2f µs/op\n", m5a_per);
printf(" M5b 1-WG noop dispatch submit+wait: %.2f µs/op\n", m5b_per);
printf(" delta (per-vkCmdDispatch + per-pipeline-bind): %.2f µs\n",
m5b_per - m5a_per);
printf("\n");
printf(" Implication for kernel batching:\n");
printf(" if QPU IDCT8 = ~ 100ns/block (best case, hypothetical),\n");
printf(" a single-block dispatch costs %.0fx more in overhead\n",
m5b_per * 1e3 / 100.0);
printf(" -> batch at least %.0f blocks per dispatch to break even.\n",
m5b_per * 1e3 / 100.0);
/* ---- Tear down (minimal — process exit handles the rest) ---- */
vkDestroyPipeline(dev, pipe, NULL);
vkDestroyShaderModule(dev, shader, NULL);
vkDestroyPipelineLayout(dev, playout, NULL);
vkDestroyCommandPool(dev, pool, NULL);
vkDestroyDevice(dev, NULL);
vkDestroyInstance(instance, NULL);
return 0;
}
+5
View File
@@ -0,0 +1,5 @@
#version 450
// Empty compute shader for measuring Vulkan dispatch overhead (M5).
// Reads nothing, writes nothing — pure dispatch round-trip floor.
layout(local_size_x = 64, local_size_y = 1, local_size_z = 1) in;
void main() {}
+114
View File
@@ -0,0 +1,114 @@
/*
* Standalone bit-exact C reference for VP9 8×8 DCT_DCT inverse
* transform + add (8-bit pixels), transcribed from the spec
* structure as represented in FFmpeg's libavcodec/vp9dsp_template.c
* (vendored under external/ffmpeg-snapshot/ at commit f46e514).
*
* Provided as a self-contained translation unit so the harness
* doesn't need to wrestle FFmpeg's BIT_DEPTH-templated macro
* expansion. Cross-checked against the vendored reference at
* runtime (see bench_neon_idct.c::cross_check_vs_ffmpeg_c()).
*
* License: LGPL-2.1-or-later (matches the upstream reference).
*
* Spec source: VP9 specification §8.7 — Inverse transform process.
*/
#include <stdint.h>
#include <stddef.h>
#include <string.h>
/* Q14 trig constants — VP9 spec table 8.7.1.4. */
#define COSPI_16_64 11585 /* cos(pi/4) * 2^14 */
#define COSPI_24_64 6270 /* cos(3pi/8) * 2^14 */
#define COSPI_8_64 15137 /* sin(3pi/8) * 2^14 */
#define COSPI_28_64 3196 /* cos(7pi/16)* 2^14 */
#define COSPI_4_64 16069 /* sin(7pi/16)* 2^14 */
#define COSPI_20_64 9102 /* cos(5pi/16)* 2^14 */
#define COSPI_12_64 13623 /* sin(5pi/16)* 2^14 */
/* Q14 round-shift: (x + (1<<13)) >> 14, with overflow-safe widening. */
static inline int32_t qround14(int64_t x)
{
return (int32_t) ((x + (1 << 13)) >> 14);
}
static inline uint8_t clip_u8(int x)
{
return (uint8_t) (x < 0 ? 0 : x > 255 ? 255 : x);
}
/* 1-D 8-point inverse DCT, signed int32 throughout. Matches
* idct8_1d in libavcodec/vp9dsp_template.c (with the stride
* collapsed to indexed access; identical arithmetic). */
static void idct8_1d(const int32_t in[8], int32_t out[8])
{
int32_t t0a = qround14((int64_t)(in[0] + in[4]) * COSPI_16_64);
int32_t t1a = qround14((int64_t)(in[0] - in[4]) * COSPI_16_64);
int32_t t2a = qround14((int64_t)in[2] * COSPI_24_64 - (int64_t)in[6] * COSPI_8_64);
int32_t t3a = qround14((int64_t)in[2] * COSPI_8_64 + (int64_t)in[6] * COSPI_24_64);
int32_t t4a = qround14((int64_t)in[1] * COSPI_28_64 - (int64_t)in[7] * COSPI_4_64);
int32_t t5a = qround14((int64_t)in[5] * COSPI_12_64 - (int64_t)in[3] * COSPI_20_64);
int32_t t6a = qround14((int64_t)in[5] * COSPI_20_64 + (int64_t)in[3] * COSPI_12_64);
int32_t t7a = qround14((int64_t)in[1] * COSPI_4_64 + (int64_t)in[7] * COSPI_28_64);
int32_t t0 = t0a + t3a, t1 = t1a + t2a;
int32_t t2 = t1a - t2a, t3 = t0a - t3a;
int32_t t4 = t4a + t5a;
int32_t t5p = t4a - t5a;
int32_t t7 = t7a + t6a;
int32_t t6p = t7a - t6a;
int32_t t5 = qround14((int64_t)(t6p - t5p) * COSPI_16_64);
int32_t t6 = qround14((int64_t)(t6p + t5p) * COSPI_16_64);
out[0] = t0 + t7; out[1] = t1 + t6;
out[2] = t2 + t5; out[3] = t3 + t4;
out[4] = t3 - t4; out[5] = t2 - t5;
out[6] = t1 - t6; out[7] = t0 - t7;
}
/* Public reference entry point. Signature matches
* ff_vp9_idct_idct_8x8_add_neon. After the call, *block is
* zeroed (matches FFmpeg behaviour). */
void daedalus_vp9_idct_idct_8x8_add_ref(uint8_t *dst, ptrdiff_t stride,
int16_t *block, int eob)
{
int32_t tmp[64];
int32_t out[8];
int32_t col[8];
/* DC-only fast path: (((coef * 11585) Q14) * 11585) Q14, then
* broadcast (+16) >> 5 added to every pixel. */
if (eob == 1) {
int32_t dc = qround14(qround14((int64_t)block[0] * COSPI_16_64)
* (int64_t) COSPI_16_64);
block[0] = 0;
int32_t add = (dc + 16) >> 5;
for (int r = 0; r < 8; r++)
for (int c = 0; c < 8; c++)
dst[r * stride + c] = clip_u8(dst[r * stride + c] + add);
return;
}
/* 8 column passes, transposed write: IDCT of block column i lands
* in row i of tmp. This matches FFmpeg's idct_idct_8x8_add_c which
* uses `tmp + i*8` as the column-pass output base — the transpose
* is implicit in the offset pattern, making the row pass below
* read columns of tmp and write columns of dst. */
for (int i = 0; i < 8; i++) {
for (int r = 0; r < 8; r++) col[r] = block[r * 8 + i];
idct8_1d(col, out);
for (int r = 0; r < 8; r++) tmp[i * 8 + r] = out[r];
}
memset(block, 0, 64 * sizeof(*block));
/* 8 row passes: column i of tmp -> column i of dst (matches
* FFmpeg's `dst[j*stride] = out[j]; dst++` pattern). */
for (int i = 0; i < 8; i++) {
for (int r = 0; r < 8; r++) col[r] = tmp[r * 8 + i];
idct8_1d(col, out);
for (int r = 0; r < 8; r++)
dst[r * stride + i] = clip_u8(dst[r * stride + i]
+ ((out[r] + 16) >> 5));
}
}