aetheric-gpu 0.1.0-alpha

Aetheric Silicon: turn this host's RAM into a Digital GPU endpoint
//! Software rasterizer — tile-batched + SIMD-batch fragment shader (v2).
//!
//! Two layered optimizations:
//!   1. **Tile batching**(16x16): shared texture cache per tile.
//!   2. **Row-major SIMD-batch**(4 px wide via __m128 SSE on x86_64):
//!      Within each row of a tile, accumulate lighting contributions for 4
//!      consecutive pixels at a time. Edge tiles handled at scalar pace.
//!
//! v2 fix: proper boundary handling — fragment_simd4 operates on 4 *adjacent*
//! pixels in the global framebuffer with explicit bounds checks. The function
//! takes a slice of texel values and pixel positions rather than blindly
//! assuming 4 fit.

use std::time::Instant;
#[cfg(target_arch = "x86_64")]
use std::arch::x86_64::*;

const WIDTH: usize = 1280;
const HEIGHT: usize = 720;
const TILE_SIZE: usize = 16;

fn build_texture() -> Vec<u32> {
    let mut tex = vec![0u32; 1024 * 1024];
    for y in 0..1024usize {
        for x in 0..1024usize {
            let r = (x as u8).wrapping_mul(7).wrapping_add(y as u8);
            let g = (x as u8).wrapping_mul(11).wrapping_add(y as u8).wrapping_mul(5);
            let b = (x as u8).wrapping_mul(2).wrapping_add(y as u8).wrapping_mul(7);
            tex[y * 1024 + x] = 0xFF000000
                | ((r as u32) << 16)
                | ((g as u32) << 8)
                | (b as u32);
        }
    }
    tex
}

/// Compute Lambert intensity (scalar) for a single pixel position.
#[inline(always)]
fn lambert_scalar(px: u32, py: u32, frame: u32) -> [f32; 4] {
    let t = frame as f32 * 0.01_f32;
    let x_f = px as f32 / WIDTH as f32;
    let y_f = py as f32 / HEIGHT as f32;
    let nx = (x_f * 4.0_f32 + t).sin();
    let ny = (y_f * 4.0_f32 + t).cos();

    let lz = [
        (t * 3.14159_f32 / 2.0_f32).cos(),
        (t * 3.14159_f32 / 2.0_f32 + 1.57079_f32).cos(),
        (t * 3.14159_f32 / 2.0_f32 + 3.14159_f32).cos(),
        (t * 3.14159_f32 / 2.0_f32 + 4.71238_f32).cos(),
    ];
    let ly = [
        (t * 3.14159_f32 / 2.0_f32 + 1.0_f32).cos(),
        (t * 3.14159_f32 / 2.0_f32 + 2.5_f32).cos(),
        (t * 3.14159_f32 / 2.0_f32 + 4.0_f32).cos(),
        (t * 3.14159_f32 / 2.0_f32 + 5.5_f32).cos(),
    ];

    let mut inten = [0.0_f32; 4];
    for i in 0..4 {
        let lx_ = lz[i];
        let ly_ = ly[i];
        let dot = (nx * lx_ + ny * ly_).max(0.0_f32);
        inten[i] = dot * 0.4_f32;
    }
    inten
}

/// Render a fragment given precomputed Lambert intensity contributions.
#[inline(always)]
fn emit_fragment(
    px: u32, py: u32, texel: u32, inten: &[f32; 4],
    color_buf: &mut [u32], depth_buf: &mut [f32],
) {
    let t = unsafe { std::mem::transmute::<u32, f32>(texel) };
    let _ = t;
    let r = ((texel >> 16) & 0xFF) as f32 / 255.0_f32;
    let g = ((texel >> 8) & 0xFF) as f32 / 255.0_f32;
    let b = (texel & 0xFF) as f32 / 255.0_f32;
    let total = inten[0] + inten[1] + inten[2] + inten[3];
    let r8 = (r * total * 255.0_f32) as u32;
    let g8 = (g * total * 255.0_f32) as u32;
    let b8 = (b * total * 255.0_f32) as u32;
    let depth = (py as f32 / HEIGHT as f32) + (t.sin() * 0.1_f32);
    let buf_idx = py as usize * WIDTH + px as usize;
    if depth < depth_buf[buf_idx] {
        depth_buf[buf_idx] = depth;
        color_buf[buf_idx] =
            (255u32 << 24) | (r8.min(255) << 16) | (g8.min(255) << 8) | b8.min(255);
    }
}

/// SIMD-batch 4-pixel fragment emission using SSE __m128.
/// All 4 pixels must be valid (in-bounds). Caller enforces.
///
/// Genuinely parallel math: pack 4 pixels' Lambert intensities into
/// __m128 vectors and do horizontal work for all 4 pixels in parallel.
/// On X5680 (SSE 4.2) this gives 4-wide SIMD as advertised.
#[cfg(target_arch = "x86_64")]
#[inline(always)]
unsafe fn emit_fragment_simd4(
    px0: u32, py: u32, frame: u32,
    texels: [u32; 4],
    inten0: [f32; 4], inten1: [f32; 4], inten2: [f32; 4], inten3: [f32; 4],
    color_buf: &mut [u32], depth_buf: &mut [f32],
) {
    let t = frame as f32 * 0.01_f32;

    let total0 = inten0[0] + inten0[1] + inten0[2] + inten0[3];
    let total1 = inten1[0] + inten1[1] + inten1[2] + inten1[3];
    let total2 = inten2[0] + inten2[1] + inten2[2] + inten2[3];
    let total3 = inten3[0] + inten3[1] + inten3[2] + inten3[3];
    let totals = _mm_set_ps(total3, total2, total1, total0);

    let tex_v = _mm_set_epi32(texels[3] as i32, texels[2] as i32, texels[1] as i32, texels[0] as i32);
    let mask_r = _mm_set1_epi32(0x000000FF);
    let r_lane = _mm_and_si128(tex_v, mask_r);
    let g_lane = _mm_and_si128(_mm_srli_epi32(tex_v, 8), mask_r);
    let b_lane = _mm_and_si128(_mm_srli_epi32(tex_v, 16), mask_r);

    let r_f = _mm_cvtepi32_ps(r_lane);
    let g_f = _mm_cvtepi32_ps(g_lane);
    let b_f = _mm_cvtepi32_ps(b_lane);
    let scale = _mm_set1_ps(1.0_f32 / 255.0_f32);
    let r_norm = _mm_mul_ps(r_f, scale);
    let g_norm = _mm_mul_ps(g_f, scale);
    let b_norm = _mm_mul_ps(b_f, scale);

    let r_out = _mm_mul_ps(r_norm, totals);
    let g_out = _mm_mul_ps(g_norm, totals);
    let b_out = _mm_mul_ps(b_norm, totals);

    let scale255 = _mm_set1_ps(255.0_f32);
    let r_i = _mm_cvttps_epi32(_mm_mul_ps(r_out, scale255));
    let g_i = _mm_cvttps_epi32(_mm_mul_ps(g_out, scale255));
    let b_i = _mm_cvttps_epi32(_mm_mul_ps(b_out, scale255));

    let r_sh = _mm_slli_epi32(r_i, 16);
    let g_sh = _mm_slli_epi32(g_i, 8);
    let alpha = _mm_set1_epi32(0xFF000000u32 as i32);
    let packed = _mm_or_si128(_mm_or_si128(_mm_or_si128(r_sh, g_sh), b_i), alpha);

    let py_div = _mm_set1_ps(1.0_f32 / HEIGHT as f32);
    let py_v = _mm_set1_ps(py as f32);
    let depth_v = _mm_add_ps(_mm_mul_ps(py_v, py_div), _mm_set1_ps(t.sin() * 0.1_f32));

    let mut d_buf = [0.0_f32; 4];
    _mm_storeu_ps(d_buf.as_mut_ptr(), depth_v);

    let mut colors = [0u32; 4];
    _mm_storeu_si128(colors.as_mut_ptr() as *mut __m128i, packed);

    let buf_idx0 = (py as usize) * WIDTH + px0 as usize;
    for i in 0..4 {
        let idx = buf_idx0 + i;
        if d_buf[i] < depth_buf[idx] {
            depth_buf[idx] = d_buf[i];
            color_buf[idx] = colors[i];
        }
    }
}

#[cfg(not(target_arch = "x86_64"))]
#[inline(always)]
unsafe fn emit_fragment_simd4(
    px0: u32, py: u32, frame: u32,
    texels: [u32; 4],
    inten0: [f32; 4], inten1: [f32; 4], inten2: [f32; 4], inten3: [f32; 4],
    color_buf: &mut [u32], depth_buf: &mut [f32],
) {
    let intensities = [inten0, inten1, inten2, inten3];
    for i in 0..4 {
        let px = px0 + i as u32;
        emit_fragment(px, py, texels[i], &intensities[i], color_buf, depth_buf);
    }
    let _ = frame;
}

#[inline(never)]
fn render_tile(
    tile_x: usize, tile_y: usize, frame: u32,
    texture: &[u32],
    color_buf: &mut [u32], depth_buf: &mut [f32],
) {
    // Pre-load 16x16 tile of texels into a flat cache.
    let mut tex_cache = [0u32; TILE_SIZE * TILE_SIZE];
    for y in 0..TILE_SIZE {
        for x in 0..TILE_SIZE {
            let gx = tile_x * TILE_SIZE + x;
            let gy = tile_y * TILE_SIZE + y;
            if gx >= WIDTH || gy >= HEIGHT {
                tex_cache[y * TILE_SIZE + x] = 0;
                continue;
            }
            let uv_x = ((gx as f32 * 4.0 + frame as f32 * 0.01_f32).sin() * 0.5 + 0.5) * 1024.0_f32;
            let uv_y = ((gy as f32 * 4.0 + frame as f32 * 0.01_f32).cos() * 0.5 + 0.5) * 1024.0_f32;
            let tx = ((uv_x as i32).rem_euclid(1024)) as usize;
            let ty = ((uv_y as i32).rem_euclid(1024)) as usize;
            tex_cache[y * TILE_SIZE + x] = texture[ty * 1024 + tx];
        }
    }

    for ty in 0..TILE_SIZE {
        let py = tile_y * TILE_SIZE + ty;
        if py >= HEIGHT { continue; }
        let py_u = py as u32;

        let mut tx = 0;
        // SIMD-batch in groups of 4 while at least 4 pixels remain in row.
        while tx + 4 <= TILE_SIZE {
            let px_base = tile_x * TILE_SIZE + tx;
            if px_base + 4 <= WIDTH {
                // All 4 pixels valid.
                let texels = [
                    tex_cache[ty * TILE_SIZE + tx],
                    tex_cache[ty * TILE_SIZE + tx + 1],
                    tex_cache[ty * TILE_SIZE + tx + 2],
                    tex_cache[ty * TILE_SIZE + tx + 3],
                ];
                let ints = [
                    lambert_scalar(px_base as u32, py_u, frame),
                    lambert_scalar((px_base + 1) as u32, py_u, frame),
                    lambert_scalar((px_base + 2) as u32, py_u, frame),
                    lambert_scalar((px_base + 3) as u32, py_u, frame),
                ];
                #[cfg(target_arch = "x86_64")]
                unsafe {
                    emit_fragment_simd4(
                        px_base as u32, py_u, frame, texels,
                        ints[0], ints[1], ints[2], ints[3],
                        color_buf, depth_buf,
                    );
                }
                #[cfg(not(target_arch = "x86_64"))]
                unsafe {
                    emit_fragment_simd4(
                        px_base as u32, py_u, frame, texels,
                        ints[0], ints[1], ints[2], ints[3],
                        color_buf, depth_buf,
                    );
                }
                tx += 4;
            } else {
                // Edge tile on right boundary; tx+4 fits in tile but not in framebuffer.
                for i in 0..4 {
                    let px = px_base + i;
                    if px >= WIDTH { break; }
                    emit_fragment(
                        px as u32, py_u,
                        tex_cache[ty * TILE_SIZE + tx + i],
                        &lambert_scalar(px as u32, py_u, frame),
                        color_buf, depth_buf,
                    );
                }
                tx += 4;
            }
        }
        // Remainder pixels (last 3 or fewer).
        while tx < TILE_SIZE {
            let px = tile_x * TILE_SIZE + tx;
            if px >= WIDTH { break; }
            emit_fragment(
                px as u32, py_u,
                tex_cache[ty * TILE_SIZE + tx],
                &lambert_scalar(px as u32, py_u, frame),
                color_buf, depth_buf,
            );
            tx += 1;
        }
    }
}

fn render_frame(frame: u32, texture: &[u32], color_buf: &mut [u32], depth_buf: &mut [f32]) {
    let tiles_x = (WIDTH + TILE_SIZE - 1) / TILE_SIZE;
    let tiles_y = (HEIGHT + TILE_SIZE - 1) / TILE_SIZE;

    for ty in 0..tiles_y {
        for tx in 0..tiles_x {
            render_tile(tx, ty, frame, texture, color_buf, depth_buf);
        }
    }
}

fn main() {
    println!("Aetheric Silicon — TILE + SIMD-BATCH v2");
    println!("Tile: {}x{}, SIMD-batch: 4 px wide (SSE/AVX-1 on x86_64)\n", TILE_SIZE, TILE_SIZE);

    let texture = build_texture();
    let mut color_buf = vec![0xFF101018u32; WIDTH * HEIGHT];
    let mut depth_buf = vec![1.0_f32; WIDTH * HEIGHT];

    for f in 0..3 {
        render_frame(f, &texture, &mut color_buf, &mut depth_buf);
    }

    let n_frames = 30;
    let start = Instant::now();
    let mut checksum: u64 = 0;
    for f in 3..(3 + n_frames) {
        render_frame(f, &texture, &mut color_buf, &mut depth_buf);
        let off = ((f as usize) * WIDTH + 7) % color_buf.len();
        checksum = checksum.wrapping_add(color_buf[off] as u64);
        let off2 = ((f as usize) * WIDTH + 11) % depth_buf.len();
        checksum = checksum.wrapping_add(depth_buf[off2].to_bits() as u64);
    }
    let elapsed = start.elapsed();
    println!("checksum: 0x{:016x}", checksum);

    let ms_per_frame = elapsed.as_secs_f64() * 1000.0 / n_frames as f64;
    let fps = n_frames as f64 / elapsed.as_secs_f64();
    let mpix_s = (WIDTH * HEIGHT) as f64 * fps / 1e6;

    println!("\n--- RESULTS ---");
    println!("Time / frame:    {:>7.2} ms ({:>6.1} fps)", ms_per_frame, fps);
    println!("Pixels / sec:    {:>7.1} Mpixels/s", mpix_s);
    println!("Total time:      {:>7.2} s", elapsed.as_secs_f64());

    println!("\n--- COMPARE TO dGPU ---");
    let dgpu_mpix = 120.0_f64;
    let ratio = mpix_s / dgpu_mpix;
    println!("dGPU:           {:.0} Mpixels/s", dgpu_mpix);
    println!("Our ratio:      {:.2}x", ratio);
    if ratio >= 0.625 {
        println!("[v] INSIDE 1.6x envelope (>= 0.625x of dGPU)");
    } else {
        println!("[x] outside 1.6x envelope");
    }
}