Task: collision_detection
Category: Spatial--temporal algorithms
Difficulty: 3
Tags: broad_phase,narrow_phase,SAT,uniform_grid,AABB,spatial_hashing
Source: Box2D (separating-axis theorem reference): https://github.com/erincatto/box2d
Description: Batched 2D convex-polygon collision detection using bounding boxes and the separating-axis theorem.

Prompt template (YAML):
# 2D Convex Polygon Collision Detection -- Prompt Template

task_description: |
  Given N 2D convex polygons (with precomputed axis-aligned bounding boxes), determine
  all collision pairs and output the collision count per polygon.

  Each polygon is a convex hull with 3-8 vertices (counterclockwise order).
  Polygons are stored in a CSR-like format: `poly_offsets[i]` to `poly_offsets[i+1]`
  gives the vertex range in `vertices_x` and `vertices_y`.

  The classic two-phase approach:
  1. **Broad Phase**: Use a uniform spatial grid to quickly identify candidate pairs
     that might collide (based on AABB overlap within the same grid cell).
  2. **Narrow Phase**: For each candidate pair, use the Separating Axis Theorem (SAT)
     to determine exact collision.

  **Separating Axis Theorem (SAT)**:
  Two convex polygons overlap if and only if there is no separating axis.
  For each edge of both polygons, compute the edge normal, project all vertices of
  both polygons onto this axis, and check if the projections overlap.
  If any axis separates them, they don't collide.

  `counts[i]` = number of distinct polygons that polygon i collides with.

interface: |
  ```c
  extern "C" void solution_compute(
      int N,                      // number of polygons
      int total_verts,            // total vertices across all polygons
      int world_size_x100,        // world dimension * 100 (decode: / 100.0f)
      int cell_size_x100,         // grid cell size * 100 (decode: / 100.0f)
      const int* poly_offsets,    // host, [N+1], CSR offsets into vertex arrays
      const float* vertices_x,   // host, [total_verts], x-coordinates
      const float* vertices_y,   // host, [total_verts], y-coordinates
      const float* aabb,          // host, [N*4], per-polygon [min_x, min_y, max_x, max_y]
      int* counts                 // host, [N], output: collision count per polygon
  );

  ```

input_size_notes: |
  - **small**: 1,000 polygons, ~5,000 vertices
  - **medium**: 10,000 polygons, ~50,000 vertices
  - **large**: 50,000 polygons, ~250,000 vertices
  - Polygons are clustered (dense regions with many collisions, sparse regions with few)
  - Vertices are counterclockwise, 3-8 per polygon
  - `world_size` and `cell_size` are encoded as int * 100 (divide by 100.0f to decode)
  - AABB layout: `aabb[i*4+0]=min_x, aabb[i*4+1]=min_y, aabb[i*4+2]=max_x, aabb[i*4+3]=max_y`

output_contract: |
  - `counts[i]` = number of distinct polygons colliding with polygon i
  - Counts must match CPU reference exactly (integer values)
  - Each collision pair (i,j) increments both `counts[i]` and `counts[j]` by 1
  - Must not double-count: if polygons i and j collide, count it once per polygon
  - `solution_compute` must be **idempotent** and **synchronize** before returning

algorithm_background: |
  **Broad Phase — Uniform Grid**:
  1. Divide the world into cells of size `cell_size × cell_size`
  2. For each polygon, determine which grid cells its AABB overlaps
  3. Insert the polygon into each overlapping cell
  4. For each cell, generate candidate pairs from polygons sharing that cell
  5. Deduplicate candidate pairs (same pair may appear in multiple cells)

  **Narrow Phase — Separating Axis Theorem**:
  For two convex polygons A (na vertices) and B (nb vertices):
  1. For each edge of A and B (na + nb edges total):
     a. Compute edge normal: edge (dx, dy) → normal (-dy, dx)
     b. Project all vertices of A onto normal → interval [min_a, max_a]
     c. Project all vertices of B onto normal → interval [min_b, max_b]
     d. If max_a < min_b or max_b < min_a → separating axis found → NO collision
  2. If no separating axis found among all edges → collision

  Complexity: O(na × nb) per pair for SAT, O(C × V²) total for C candidate pairs.

hints_l2: |
  - Copy polygon data to GPU at the start of `solution_compute`
  - Broad phase: build grid on GPU, use spatial hashing for cell assignment
  - Narrow phase: one thread per candidate pair, each running SAT independently
  - Use `atomicAdd` to increment collision counts
  - Deduplicate candidate pairs to avoid double-counting

hints_l1: |
  ### 1. Persistent GPU Data
  Copy `poly_offsets`, `vertices_x`, `vertices_y`, `aabb` to device at the start of `solution_compute`.
  Pre-allocate output `counts` buffer on device.

  ### 2. Grid Construction (GPU)
  - Kernel 1: For each polygon, compute cell range from AABB, emit (cell_hash, poly_id) pairs
  - Sort by cell_hash using `thrust::sort_by_key`
  - Kernel 2: Find cell boundaries (start/end of each cell in sorted array)

  ### 3. Candidate Pair Generation
  For each cell, generate all (i, j) pairs where i < j:
  - Option A: One thread per cell, iterate pairs (load-imbalanced for dense cells)
  - Option B: Flatten pairs with prefix sum, one thread per pair (better balance)

  ### 4. Deduplication
  Since the same pair can appear in multiple cells:
  - Sort candidate pairs by (i, j)
  - Remove consecutive duplicates (thrust::unique)

  ### 5. Parallel SAT
  One thread per unique candidate pair:
  - Load vertices of both polygons (can use shared memory for the pair)
  - Test all edge normals of both polygons
  - Early exit on first separating axis
  - If collision: `atomicAdd(&counts[i], 1); atomicAdd(&counts[j], 1);`

  ### 6. AABB Pre-filter
  Before running SAT, test AABB overlap (4 comparisons) as a cheap early exit.


CPU reference implementation (C):
// cpu_reference.c -- collision_detection CPU baseline
//
// Grid broad phase + SAT narrow phase for 2D convex polygon collision.
// NO file I/O. All I/O handled by task_io_cpu.c.

#include <stdlib.h>
#include <string.h>
#include <math.h>

// ===== Module-level state =====
static int g_N;
static int g_total_verts;
static float g_world_size;
static float g_cell_size;
static int g_grid_dim;

static const int*   g_poly_offsets;
static const float* g_vx;
static const float* g_vy;
static const float* g_aabb;  // N*4: min_x, min_y, max_x, max_y

// Grid: cell_heads[cell] = first polygon index, next[entry] = next polygon
// Using a flat hash map approach
#define MAX_GRID_CELLS 1000003  // prime for hashing
#define MAX_ENTRIES    8000000  // max polygon-to-cell entries

static int* g_cell_heads;    // MAX_GRID_CELLS, -1 = empty
static int* g_entry_poly;    // MAX_ENTRIES: polygon ID for each entry
static int* g_entry_next;    // MAX_ENTRIES: next entry in same cell
static int  g_num_entries;

// Candidate pairs (sorted unique)
#define MAX_PAIRS 20000000
static int* g_pair_i;
static int* g_pair_j;
static int  g_num_pairs;

// ===== Grid hash =====
static int grid_hash(int gx, int gy) {
    unsigned int h = (unsigned int)(gx * 73856093) ^ (unsigned int)(gy * 19349663);
    return (int)(h % MAX_GRID_CELLS);
}

// ===== AABB overlap =====
static int aabb_test(int i, int j) {
    return g_aabb[i*4+2] >= g_aabb[j*4+0] &&
           g_aabb[j*4+2] >= g_aabb[i*4+0] &&
           g_aabb[i*4+3] >= g_aabb[j*4+1] &&
           g_aabb[j*4+3] >= g_aabb[i*4+1];
}

// ===== SAT test =====
static int sat_test(int poly_a, int poly_b) {
    int sa = g_poly_offsets[poly_a], ea = g_poly_offsets[poly_a + 1];
    int sb = g_poly_offsets[poly_b], eb = g_poly_offsets[poly_b + 1];
    int na = ea - sa, nb = eb - sb;

    // Test edges of polygon A
    for (int i = 0; i < na; i++) {
        int j = (i + 1) % na;
        float nx = -(g_vy[sa + j] - g_vy[sa + i]);
        float ny =   g_vx[sa + j] - g_vx[sa + i];

        float min_a = 1e30f, max_a = -1e30f;
        for (int k = 0; k < na; k++) {
            float p = nx * g_vx[sa + k] + ny * g_vy[sa + k];
            if (p < min_a) min_a = p;
            if (p > max_a) max_a = p;
        }

        float min_b = 1e30f, max_b = -1e30f;
        for (int k = 0; k < nb; k++) {
            float p = nx * g_vx[sb + k] + ny * g_vy[sb + k];
            if (p < min_b) min_b = p;
            if (p > max_b) max_b = p;
        }

        if (max_a < min_b || max_b < min_a) return 0;
    }

    // Test edges of polygon B
    for (int i = 0; i < nb; i++) {
        int j = (i + 1) % nb;
        float nx = -(g_vy[sb + j] - g_vy[sb + i]);
        float ny =   g_vx[sb + j] - g_vx[sb + i];

        float min_a = 1e30f, max_a = -1e30f;
        for (int k = 0; k < na; k++) {
            float p = nx * g_vx[sa + k] + ny * g_vy[sa + k];
            if (p < min_a) min_a = p;
            if (p > max_a) max_a = p;
        }

        float min_b = 1e30f, max_b = -1e30f;
        for (int k = 0; k < nb; k++) {
            float p = nx * g_vx[sb + k] + ny * g_vy[sb + k];
            if (p < min_b) min_b = p;
            if (p > max_b) max_b = p;
        }

        if (max_a < min_b || max_b < min_a) return 0;
    }

    return 1;
}

// Sort packed pairs
static int cmp_ll(const void* a, const void* b) {
    long long va = *(const long long*)a;
    long long vb = *(const long long*)b;
    return (va > vb) - (va < vb);
}

void solution_compute(int N, int total_verts,
                      int world_size_x100, int cell_size_x100,
                      const int* poly_offsets,
                      const float* vertices_x, const float* vertices_y,
                      const float* aabb,
                      int* counts) {
    g_N = N;
    g_total_verts = total_verts;
    g_world_size = (float)world_size_x100 / 100.0f;
    g_cell_size = (float)cell_size_x100 / 100.0f;
    g_grid_dim = (int)ceilf(g_world_size / g_cell_size) + 1;
    g_poly_offsets = poly_offsets;
    g_vx = vertices_x;
    g_vy = vertices_y;
    g_aabb = aabb;

    g_cell_heads = (int*)malloc(MAX_GRID_CELLS * sizeof(int));
    g_entry_poly = (int*)malloc(MAX_ENTRIES * sizeof(int));
    g_entry_next = (int*)malloc(MAX_ENTRIES * sizeof(int));
    g_pair_i = (int*)malloc(MAX_PAIRS * sizeof(int));
    g_pair_j = (int*)malloc(MAX_PAIRS * sizeof(int));

    memset(counts, 0, (size_t)N * sizeof(int));

    // === Phase 1: Build grid ===
    memset(g_cell_heads, -1, MAX_GRID_CELLS * sizeof(int));
    g_num_entries = 0;

    for (int i = 0; i < N; i++) {
        int gx0 = (int)(g_aabb[i*4+0] / g_cell_size);
        int gy0 = (int)(g_aabb[i*4+1] / g_cell_size);
        int gx1 = (int)(g_aabb[i*4+2] / g_cell_size);
        int gy1 = (int)(g_aabb[i*4+3] / g_cell_size);
        if (gx0 < 0) gx0 = 0;
        if (gy0 < 0) gy0 = 0;

        for (int gx = gx0; gx <= gx1; gx++) {
            for (int gy = gy0; gy <= gy1; gy++) {
                int h = grid_hash(gx, gy);
                if (g_num_entries < MAX_ENTRIES) {
                    int e = g_num_entries++;
                    g_entry_poly[e] = i;
                    g_entry_next[e] = g_cell_heads[h];
                    g_cell_heads[h] = e;
                }
            }
        }
    }

    // === Phase 2: Collect candidate pairs from grid cells ===
    g_num_pairs = 0;

    for (int c = 0; c < MAX_GRID_CELLS; c++) {
        if (g_cell_heads[c] < 0) continue;

        // Collect polygons in this cell
        int cell_polys[512];
        int cell_count = 0;
        int e = g_cell_heads[c];
        while (e >= 0 && cell_count < 512) {
            cell_polys[cell_count++] = g_entry_poly[e];
            e = g_entry_next[e];
        }

        // Generate pairs
        for (int a = 0; a < cell_count; a++) {
            for (int b = a + 1; b < cell_count; b++) {
                int pi = cell_polys[a], pj = cell_polys[b];
                if (pi > pj) { int tmp = pi; pi = pj; pj = tmp; }
                if (g_num_pairs < MAX_PAIRS) {
                    g_pair_i[g_num_pairs] = pi;
                    g_pair_j[g_num_pairs] = pj;
                    g_num_pairs++;
                }
            }
        }
    }

    // === Phase 3: Sort and deduplicate pairs ===
    long long* packed = (long long*)malloc((size_t)g_num_pairs * sizeof(long long));
    for (int k = 0; k < g_num_pairs; k++) {
        packed[k] = ((long long)g_pair_i[k] << 32) | (unsigned int)g_pair_j[k];
    }
    qsort(packed, (size_t)g_num_pairs, sizeof(long long), cmp_ll);

    // === Phase 4: Deduplicate + AABB test + SAT ===
    int prev = -1;
    for (int k = 0; k < g_num_pairs; k++) {
        // Skip duplicates
        if (k > 0 && packed[k] == packed[k - 1]) continue;

        int pi = (int)(packed[k] >> 32);
        int pj = (int)(packed[k] & 0xFFFFFFFF);

        // AABB test
        if (!aabb_test(pi, pj)) continue;

        // SAT test
        if (sat_test(pi, pj)) {
            counts[pi]++;
            counts[pj]++;
        }
    }

    free(packed);
    free(g_cell_heads);
    free(g_entry_poly);
    free(g_entry_next);
    free(g_pair_i);
    free(g_pair_j);
    g_cell_heads = NULL;
    g_entry_poly = NULL;
    g_entry_next = NULL;
    g_pair_i = NULL;
    g_pair_j = NULL;
}

void solution_free(void) {
    // All cleanup is done at the end of solution_compute
}