From: Rune QQ <key@rune.みんな>
Date: Wed, 29 Jan 2026
Subject: slain(nf4): added math and tests

NF4 4-bit quantization with CUDA kernels.
8x compression, minimal quality loss.

Signed-off-by: Rune QQ <key@rune.みんな>
---

diff --git a/candle-core/src/nf4/cuda.rs b/candle-core/src/nf4/cuda.rs
new file mode 100644
index 00000000..05cc3162
--- /dev/null
+++ b/candle-core/src/nf4/cuda.rs
@@ -0,0 +1,128 @@
+//! CUDA backend for NF4 dequantization
+
+use super::Nf4Storage;
+use crate::cuda_backend::WrapErr;
+use crate::{CudaDevice, CudaStorage, Device, DType, Result, Shape, Storage, Tensor};
+use cudarc::driver::PushKernelArg;
+
+/// Dequantize NF4 storage to tensor using CUDA kernel
+pub fn dequantize_cuda(storage: &Nf4Storage, device: &Device) -> Result<Tensor> {
+    let Device::Cuda(cuda_dev) = device else {
+        crate::bail!("Expected CUDA device for CUDA dequantization")
+    };
+
+    let num_elements = storage.num_elements();
+    let block_size = storage.block_size();
+    let shape: Shape = storage.shape().into();
+
+    // Copy packed data and absmax to GPU
+    let packed_gpu = cuda_dev.clone_htod(storage.packed_data())?;
+    let absmax_gpu = cuda_dev.clone_htod(storage.absmax())?;
+
+    // Allocate output buffer
+    let output_gpu = cuda_dev.alloc_zeros::<f32>(num_elements)?;
+
+    // Get kernel function (C++ mangled name for nf4_dequantize_kernel<float>)
+    let func = cuda_dev
+        .get_or_load_func("_Z21nf4_dequantize_kernelIfEvPKhPKfPT_ii", &candle_kernels::NF4)?;
+
+    // Calculate grid/block dimensions
+    // Each thread handles 2 elements (1 byte)
+    let threads_needed = (num_elements + 1) / 2;
+    let block_dim = 256u32;
+    let grid_dim = ((threads_needed as u32 + block_dim - 1) / block_dim).max(1);
+
+    let cfg = cudarc::driver::LaunchConfig {
+        grid_dim: (grid_dim, 1, 1),
+        block_dim: (block_dim, 1, 1),
+        shared_mem_bytes: 0,
+    };
+
+    // Launch kernel
+    let n_arg = num_elements as i32;
+    let bs_arg = block_size as i32;
+
+    let mut builder = func.builder();
+    builder.arg(&packed_gpu);
+    builder.arg(&absmax_gpu);
+    builder.arg(&output_gpu);
+    builder.arg(&n_arg);
+    builder.arg(&bs_arg);
+
+    unsafe { builder.launch(cfg) }.w()?;
+
+    // Create tensor from GPU buffer
+    let cuda_storage = CudaStorage::wrap_cuda_slice(output_gpu, cuda_dev.clone());
+    let tensor_storage = Storage::Cuda(cuda_storage);
+    let none = crate::op::BackpropOp::none();
+
+    Ok(crate::tensor::from_storage(tensor_storage, shape, none, false))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn cuda_available() -> bool {
+        Device::cuda_if_available(0).is_ok()
+    }
+
+    // Test 11: CUDA dequant produces same shape as CPU
+    #[test]
+    fn test_cuda_shape() {
+        if !cuda_available() {
+            println!("CUDA not available, skipping");
+            return;
+        }
+
+        let device = Device::new_cuda(0).unwrap();
+        let data: Vec<f32> = vec![0.5; 256];
+        let storage = Nf4Storage::from_f32_slice(&data, 256, vec![16, 16]).unwrap();
+
+        let tensor = dequantize_cuda(&storage, &device).unwrap();
+        assert_eq!(tensor.dims(), &[16, 16]);
+    }
+
+    // Test 12: CUDA matches CPU output
+    #[test]
+    fn test_cuda_matches_cpu() {
+        if !cuda_available() {
+            println!("CUDA not available, skipping");
+            return;
+        }
+
+        let cpu_device = Device::Cpu;
+        let cuda_device = Device::new_cuda(0).unwrap();
+
+        let data: Vec<f32> = (0..1024).map(|i| ((i as f32) - 512.0) / 512.0).collect();
+        let storage = Nf4Storage::from_f32_slice(&data, 1024, vec![1024]).unwrap();
+
+        let cpu_result = storage.to_tensor(&cpu_device).unwrap();
+        let cuda_result = storage.to_tensor(&cuda_device).unwrap();
+
+        // Move CUDA result to CPU for comparison
+        let cuda_on_cpu = cuda_result.to_device(&cpu_device).unwrap();
+
+        let diff = (&cpu_result - &cuda_on_cpu).unwrap();
+        let max_diff: f32 = diff.abs().unwrap().max(0).unwrap().to_scalar().unwrap();
+
+        // Should be exactly equal (same algorithm)
+        assert!(max_diff < 1e-5, "Max diff {} too large", max_diff);
+    }
+
+    // Test 13: CUDA handles non-aligned sizes
+    #[test]
+    fn test_cuda_non_aligned() {
+        if !cuda_available() {
+            println!("CUDA not available, skipping");
+            return;
+        }
+
+        let device = Device::new_cuda(0).unwrap();
+        let data: Vec<f32> = vec![0.3; 100]; // Not power of 2
+        let storage = Nf4Storage::from_f32_slice(&data, 100, vec![100]).unwrap();
+
+        let tensor = dequantize_cuda(&storage, &device).unwrap();
+        assert_eq!(tensor.dims(), &[100]);
+    }
+}
diff --git a/candle-core/src/nf4/gpu.rs b/candle-core/src/nf4/gpu.rs
new file mode 100644
index 00000000..f6978261
--- /dev/null
+++ b/candle-core/src/nf4/gpu.rs
@@ -0,0 +1,356 @@
+//! GPU-resident NF4 storage - no CPU fallback, maximum speed
+//!
+//! Quantized weights live directly in VRAM. Zero copies during forward.
+//! Supports F32, F16, and BF16 input/output without wasteful conversions.
+
+use super::NF4_BLOCK_SIZE;
+use crate::backend::BackendDevice;
+use crate::cuda_backend::WrapErr;
+use crate::{CudaDevice, CudaStorage, DType, Result, Shape, Storage, Tensor};
+use cudarc::driver::{CudaSlice, PushKernelArg};
+use half::{bf16, f16};
+
+// Kernel names for different dtypes
+const DEQUANT_F32: &str = "_Z21nf4_dequantize_kernelIfEvPKhPKfPT_ii";
+const DEQUANT_F16: &str = "_Z21nf4_dequantize_kernelI6__halfEvPKhPKfPT_ii";
+const DEQUANT_BF16: &str = "_Z21nf4_dequantize_kernelI13__nv_bfloat16EvPKhPKfPT_ii";
+
+const QUANT_F32: &str = "_Z19nf4_quantize_kernelIfEvPKT_PhPfii";
+const QUANT_F16: &str = "_Z19nf4_quantize_kernelI6__halfEvPKT_PhPfii";
+const QUANT_BF16: &str = "_Z19nf4_quantize_kernelI13__nv_bfloat16EvPKT_PhPfii";
+
+/// GPU-resident NF4 storage
+/// Quantized data lives in VRAM - zero copies during forward
+#[derive(Debug, Clone)]
+pub struct Nf4StorageGpu {
+    /// Packed 4-bit data on GPU (2 values per byte)
+    packed: CudaSlice<u8>,
+    /// Per-block absmax values on GPU
+    absmax: CudaSlice<f32>,
+    /// CUDA device
+    device: CudaDevice,
+    /// Block size (always 64)
+    block_size: usize,
+    /// Original number of elements
+    num_elements: usize,
+    /// Original shape
+    shape: Vec<usize>,
+}
+
+impl Nf4StorageGpu {
+    /// Create GPU storage from a tensor - FULL GPU QUANTIZATION
+    /// Tensor should already be on GPU for maximum speed!
+    /// Supports F32, F16, and BF16 input without wasteful conversions.
+    pub fn from_tensor(tensor: &Tensor, cuda_dev: &CudaDevice) -> Result<Self> {
+        let shape = tensor.dims().to_vec();
+        let num_elements = tensor.elem_count();
+        let block_size = NF4_BLOCK_SIZE;
+        let num_blocks = (num_elements + block_size - 1) / block_size;
+        let packed_len = (num_elements + 1) / 2;
+        let input_dtype = tensor.dtype();
+
+        // Check if already on this GPU, if not move it
+        let tensor_gpu = match tensor.device() {
+            crate::Device::Cuda(_) => tensor.clone(),
+            _ => tensor.to_device(&crate::Device::Cuda(cuda_dev.clone()))?,
+        };
+        let tensor_flat = tensor_gpu.flatten_all()?;
+
+        // Allocate output buffers on GPU
+        let absmax: CudaSlice<f32> = cuda_dev.alloc_zeros(num_blocks)?;
+        let packed: CudaSlice<u8> = cuda_dev.alloc_zeros(packed_len)?;
+
+        let n_arg = num_elements as i32;
+        let bs_arg = block_size as i32;
+
+        let threads_needed = (num_elements + 1) / 2;
+        let quant_block = 256u32;
+        let quant_grid = ((threads_needed as u32 + quant_block - 1) / quant_block).max(1);
+
+        let quant_cfg = cudarc::driver::LaunchConfig {
+            grid_dim: (quant_grid, 1, 1),
+            block_dim: (quant_block, 1, 1),
+            shared_mem_bytes: 0,
+        };
+
+        // === STEP 1: Compute absmax on GPU (needs F32 for reduction) ===
+        // For non-F32 input, we compute absmax via the quantize kernel which handles it
+        let absmax_func = cuda_dev
+            .get_or_load_func("_Z25nf4_compute_absmax_kernelPKfPfii", &candle_kernels::NF4)?;
+
+        let absmax_threads = 64u32;
+        let absmax_cfg = cudarc::driver::LaunchConfig {
+            grid_dim: (num_blocks as u32, 1, 1),
+            block_dim: (absmax_threads, 1, 1),
+            shared_mem_bytes: absmax_threads * 4,
+        };
+
+        // Get raw GPU pointer and quantize based on dtype
+        let (storage, _layout) = tensor_flat.storage_and_layout();
+        match &*storage {
+            Storage::Cuda(cuda_storage) => {
+                match input_dtype {
+                    DType::F32 => {
+                        let input_slice = cuda_storage.as_cuda_slice::<f32>()?;
+
+                        // Compute absmax
+                        {
+                            let mut builder = absmax_func.builder();
+                            builder.arg(input_slice);
+                            builder.arg(&absmax);
+                            builder.arg(&n_arg);
+                            builder.arg(&bs_arg);
+                            unsafe { builder.launch(absmax_cfg) }.w()?;
+                        }
+
+                        // Quantize
+                        let quant_func = cuda_dev.get_or_load_func(QUANT_F32, &candle_kernels::NF4)?;
+                        let mut builder = quant_func.builder();
+                        builder.arg(input_slice);
+                        builder.arg(&packed);
+                        builder.arg(&absmax);
+                        builder.arg(&n_arg);
+                        builder.arg(&bs_arg);
+                        unsafe { builder.launch(quant_cfg) }.w()?;
+                    }
+                    DType::F16 => {
+                        let input_slice = cuda_storage.as_cuda_slice::<f16>()?;
+
+                        // For F16/BF16, convert to F32 just for absmax computation
+                        let tensor_f32 = tensor_flat.to_dtype(DType::F32)?;
+                        let (storage_f32, _) = tensor_f32.storage_and_layout();
+                        if let Storage::Cuda(cs) = &*storage_f32 {
+                            let slice_f32 = cs.as_cuda_slice::<f32>()?;
+                            let mut builder = absmax_func.builder();
+                            builder.arg(slice_f32);
+                            builder.arg(&absmax);
+                            builder.arg(&n_arg);
+                            builder.arg(&bs_arg);
+                            unsafe { builder.launch(absmax_cfg) }.w()?;
+                        }
+
+                        // Quantize directly from F16
+                        let quant_func = cuda_dev.get_or_load_func(QUANT_F16, &candle_kernels::NF4)?;
+                        let mut builder = quant_func.builder();
+                        builder.arg(input_slice);
+                        builder.arg(&packed);
+                        builder.arg(&absmax);
+                        builder.arg(&n_arg);
+                        builder.arg(&bs_arg);
+                        unsafe { builder.launch(quant_cfg) }.w()?;
+                    }
+                    DType::BF16 => {
+                        let input_slice = cuda_storage.as_cuda_slice::<bf16>()?;
+
+                        // For F16/BF16, convert to F32 just for absmax computation
+                        let tensor_f32 = tensor_flat.to_dtype(DType::F32)?;
+                        let (storage_f32, _) = tensor_f32.storage_and_layout();
+                        if let Storage::Cuda(cs) = &*storage_f32 {
+                            let slice_f32 = cs.as_cuda_slice::<f32>()?;
+                            let mut builder = absmax_func.builder();
+                            builder.arg(slice_f32);
+                            builder.arg(&absmax);
+                            builder.arg(&n_arg);
+                            builder.arg(&bs_arg);
+                            unsafe { builder.launch(absmax_cfg) }.w()?;
+                        }
+
+                        // Quantize directly from BF16
+                        let quant_func = cuda_dev.get_or_load_func(QUANT_BF16, &candle_kernels::NF4)?;
+                        let mut builder = quant_func.builder();
+                        builder.arg(input_slice);
+                        builder.arg(&packed);
+                        builder.arg(&absmax);
+                        builder.arg(&n_arg);
+                        builder.arg(&bs_arg);
+                        unsafe { builder.launch(quant_cfg) }.w()?;
+                    }
+                    _ => crate::bail!("NF4 quantization only supports F32, F16, BF16, got {:?}", input_dtype),
+                }
+            }
+            _ => crate::bail!("Expected CUDA storage"),
+        }
+
+        // Sync to ensure quantization is complete
+        cuda_dev.synchronize()?;
+
+        Ok(Self {
+            packed,
+            absmax,
+            device: cuda_dev.clone(),
+            block_size,
+            num_elements,
+            shape,
+        })
+    }
+
+    /// Dequantize to tensor as F32 - pure GPU, zero copies!
+    pub fn to_tensor(&self) -> Result<Tensor> {
+        self.to_tensor_dtype(DType::F32)
+    }
+
+    /// Dequantize to tensor with specific dtype - pure GPU!
+    /// Supports F32, F16, BF16 output directly via templated kernels.
+    pub fn to_tensor_dtype(&self, dtype: DType) -> Result<Tensor> {
+        let shape: Shape = self.shape.as_slice().into();
+
+        // Grid/block dims
+        let threads_needed = (self.num_elements + 1) / 2;
+        let block_dim = 256u32;
+        let grid_dim = ((threads_needed as u32 + block_dim - 1) / block_dim).max(1);
+
+        let cfg = cudarc::driver::LaunchConfig {
+            grid_dim: (grid_dim, 1, 1),
+            block_dim: (block_dim, 1, 1),
+            shared_mem_bytes: 0,
+        };
+
+        let n_arg = self.num_elements as i32;
+        let bs_arg = self.block_size as i32;
+
+        let tensor = match dtype {
+            DType::F32 => {
+                let output_gpu: CudaSlice<f32> = self.device.alloc_zeros(self.num_elements)?;
+                let func = self.device.get_or_load_func(DEQUANT_F32, &candle_kernels::NF4)?;
+
+                let mut builder = func.builder();
+                builder.arg(&self.packed);
+                builder.arg(&self.absmax);
+                builder.arg(&output_gpu);
+                builder.arg(&n_arg);
+                builder.arg(&bs_arg);
+                unsafe { builder.launch(cfg) }.w()?;
+
+                let cuda_storage = CudaStorage::wrap_cuda_slice(output_gpu, self.device.clone());
+                crate::tensor::from_storage(Storage::Cuda(cuda_storage), shape, crate::op::BackpropOp::none(), false)
+            }
+            DType::F16 => {
+                let output_gpu: CudaSlice<f16> = self.device.alloc_zeros(self.num_elements)?;
+                let func = self.device.get_or_load_func(DEQUANT_F16, &candle_kernels::NF4)?;
+
+                let mut builder = func.builder();
+                builder.arg(&self.packed);
+                builder.arg(&self.absmax);
+                builder.arg(&output_gpu);
+                builder.arg(&n_arg);
+                builder.arg(&bs_arg);
+                unsafe { builder.launch(cfg) }.w()?;
+
+                let cuda_storage = CudaStorage::wrap_cuda_slice(output_gpu, self.device.clone());
+                crate::tensor::from_storage(Storage::Cuda(cuda_storage), shape, crate::op::BackpropOp::none(), false)
+            }
+            DType::BF16 => {
+                let output_gpu: CudaSlice<bf16> = self.device.alloc_zeros(self.num_elements)?;
+                let func = self.device.get_or_load_func(DEQUANT_BF16, &candle_kernels::NF4)?;
+
+                let mut builder = func.builder();
+                builder.arg(&self.packed);
+                builder.arg(&self.absmax);
+                builder.arg(&output_gpu);
+                builder.arg(&n_arg);
+                builder.arg(&bs_arg);
+                unsafe { builder.launch(cfg) }.w()?;
+
+                let cuda_storage = CudaStorage::wrap_cuda_slice(output_gpu, self.device.clone());
+                crate::tensor::from_storage(Storage::Cuda(cuda_storage), shape, crate::op::BackpropOp::none(), false)
+            }
+            _ => crate::bail!("NF4 dequantization only supports F32, F16, BF16, got {:?}", dtype),
+        };
+
+        Ok(tensor)
+    }
+
+    /// Get number of elements
+    pub fn num_elements(&self) -> usize {
+        self.num_elements
+    }
+
+    /// Get shape
+    pub fn shape(&self) -> &[usize] {
+        &self.shape
+    }
+
+    /// Get storage size in bytes (on GPU)
+    pub fn storage_bytes(&self) -> usize {
+        let packed_bytes = (self.num_elements + 1) / 2;
+        let absmax_bytes = ((self.num_elements + self.block_size - 1) / self.block_size) * 4;
+        packed_bytes + absmax_bytes
+    }
+
+    /// Get the CUDA device
+    pub fn device(&self) -> &CudaDevice {
+        &self.device
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::Device;
+
+    fn cuda_available() -> bool {
+        Device::cuda_if_available(0).is_ok()
+    }
+
+    #[test]
+    fn test_gpu_storage_roundtrip() {
+        if !cuda_available() {
+            println!("CUDA not available, skipping");
+            return;
+        }
+
+        let cuda_device = Device::new_cuda(0).unwrap();
+        let Device::Cuda(cuda_dev) = &cuda_device else { panic!() };
+
+        // Create test tensor
+        let data: Vec<f32> = (0..256).map(|i| (i as f32 / 128.0) - 1.0).collect();
+        let tensor = Tensor::from_vec(data.clone(), (16, 16), &cuda_device).unwrap();
+
+        // Quantize to GPU storage
+        let storage = Nf4StorageGpu::from_tensor(&tensor, cuda_dev).unwrap();
+
+        // Dequantize (pure GPU!)
+        let recovered = storage.to_tensor().unwrap();
+
+        // Check shape
+        assert_eq!(recovered.dims(), &[16, 16]);
+
+        // Check device - should still be on GPU
+        assert!(matches!(recovered.device(), Device::Cuda(_)));
+
+        // Check error
+        let recovered_cpu = recovered.to_device(&Device::Cpu).unwrap();
+        let recovered_data: Vec<f32> = recovered_cpu.flatten_all().unwrap().to_vec1().unwrap();
+
+        let max_err: f32 = data.iter()
+            .zip(recovered_data.iter())
+            .map(|(a, b)| (a - b).abs())
+            .fold(0.0, f32::max);
+
+        assert!(max_err < 0.2, "Max error {} too large", max_err);
+    }
+
+    #[test]
+    fn test_gpu_storage_compression() {
+        if !cuda_available() {
+            println!("CUDA not available, skipping");
+            return;
+        }
+
+        let cuda_device = Device::new_cuda(0).unwrap();
+        let Device::Cuda(cuda_dev) = &cuda_device else { panic!() };
+
+        let num_elements = 1024;
+        let data: Vec<f32> = (0..num_elements).map(|i| (i as f32 / 512.0) - 1.0).collect();
+        let tensor = Tensor::from_vec(data, (num_elements,), &cuda_device).unwrap();
+
+        let storage = Nf4StorageGpu::from_tensor(&tensor, cuda_dev).unwrap();
+
+        let original_bytes = num_elements * 4;
+        let compressed_bytes = storage.storage_bytes();
+        let ratio = compressed_bytes as f32 / original_bytes as f32;
+
+        // Should be ~0.14 (7x compression)
+        assert!(ratio < 0.20, "Compression ratio {} not good enough", ratio);
+    }
+}
diff --git a/candle-core/src/nf4/mod.rs b/candle-core/src/nf4/mod.rs
new file mode 100644
index 00000000..97822a64
--- /dev/null
+++ b/candle-core/src/nf4/mod.rs
@@ -0,0 +1,349 @@
+//! NF4 (4-bit NormalFloat) quantization support
+//!
+//! NF4 is a 4-bit data type optimized for normally distributed weights,
+//! as introduced in the QLoRA paper. This module provides:
+//! - CPU dequantization (always available)
+//! - CUDA dequantization (when cuda feature enabled)
+//! - GPU-resident storage (Nf4StorageGpu) - zero copies during forward!
+//!
+//! The quantized format stores:
+//! - Packed 4-bit indices (2 values per byte)
+//! - Per-block absmax values (default block size = 64)
+
+use crate::{DType, Device, Result, Tensor};
+
+#[cfg(feature = "cuda")]
+mod cuda;
+
+#[cfg(feature = "cuda")]
+mod gpu;
+
+#[cfg(feature = "cuda")]
+pub use gpu::Nf4StorageGpu;
+
+/// NF4 dequantization lookup table
+/// These are the 16 quantile values for a standard normal distribution
+pub const NF4_DEQUANT_TABLE: [f32; 16] = [
+    -1.0,
+    -0.6961928009986877,
+    -0.5250730514526367,
+    -0.39491748809814453,
+    -0.28444138169288635,
+    -0.18477343022823334,
+    -0.09105003625154495,
+    0.0,
+    0.07958029955625534,
+    0.16093020141124725,
+    0.24611230194568634,
+    0.33791524171829224,
+    0.44070982933044434,
+    0.5626170039176941,
+    0.7229568362236023,
+    1.0,
+];
+
+/// Default block size for NF4 quantization (matches bitsandbytes)
+pub const NF4_BLOCK_SIZE: usize = 64;
+
+/// NF4 quantized storage
+#[derive(Debug, Clone)]
+pub struct Nf4Storage {
+    /// Packed 4-bit data (2 values per byte)
+    data: Vec<u8>,
+    /// Per-block absmax values
+    absmax: Vec<f32>,
+    /// Block size used for quantization
+    block_size: usize,
+    /// Original number of elements
+    num_elements: usize,
+    /// Original shape
+    shape: Vec<usize>,
+}
+
+impl Nf4Storage {
+    /// Create NF4 storage from a tensor (CPU quantization)
+    pub fn from_tensor(tensor: &Tensor) -> Result<Self> {
+        let tensor = tensor.to_dtype(DType::F32)?.flatten_all()?;
+        let data: Vec<f32> = tensor.to_vec1()?;
+        let num_elements = data.len();
+        let shape = tensor.shape().dims().to_vec();
+
+        Self::from_f32_slice(&data, num_elements, shape)
+    }
+
+    /// Create NF4 storage from f32 slice
+    pub fn from_f32_slice(data: &[f32], num_elements: usize, shape: Vec<usize>) -> Result<Self> {
+        let block_size = NF4_BLOCK_SIZE;
+        let num_blocks = (num_elements + block_size - 1) / block_size;
+
+        // Compute absmax per block
+        let mut absmax = Vec::with_capacity(num_blocks);
+        for block_idx in 0..num_blocks {
+            let start = block_idx * block_size;
+            let end = (start + block_size).min(num_elements);
+            let block_max = data[start..end]
+                .iter()
+                .map(|x| x.abs())
+                .fold(0.0f32, |a, b| a.max(b));
+            absmax.push(block_max);
+        }
+
+        // Quantize to 4-bit
+        let packed_len = (num_elements + 1) / 2;
+        let mut packed = vec![0u8; packed_len];
+
+        for i in 0..num_elements {
+            let block_idx = i / block_size;
+            let scale = absmax[block_idx];
+            let normalized = if scale != 0.0 { data[i] / scale } else { 0.0 };
+            let clamped = normalized.clamp(-1.0, 1.0);
+            let quantized = quantize_nf4(clamped);
+
+            let byte_idx = i / 2;
+            if i % 2 == 0 {
+                packed[byte_idx] |= quantized << 4;
+            } else {
+                packed[byte_idx] |= quantized;
+            }
+        }
+
+        Ok(Self {
+            data: packed,
+            absmax,
+            block_size,
+            num_elements,
+            shape,
+        })
+    }
+
+    /// Dequantize to tensor
+    pub fn to_tensor(&self, device: &Device) -> Result<Tensor> {
+        match device {
+            Device::Cpu => self.dequantize_cpu(device),
+            #[cfg(feature = "cuda")]
+            Device::Cuda(_) => cuda::dequantize_cuda(self, device),
+            #[cfg(not(feature = "cuda"))]
+            Device::Cuda(_) => self.dequantize_cpu(device),
+            _ => self.dequantize_cpu(device),
+        }
+    }
+
+    /// CPU dequantization
+    fn dequantize_cpu(&self, device: &Device) -> Result<Tensor> {
+        let mut output = vec![0.0f32; self.num_elements];
+
+        for i in 0..self.num_elements {
+            let byte_idx = i / 2;
+            let packed = self.data[byte_idx];
+            let quantized = if i % 2 == 0 {
+                (packed >> 4) & 0x0F
+            } else {
+                packed & 0x0F
+            };
+
+            let block_idx = i / self.block_size;
+            let scale = self.absmax[block_idx];
+            output[i] = NF4_DEQUANT_TABLE[quantized as usize] * scale;
+        }
+
+        Tensor::from_vec(output, self.shape.as_slice(), device)
+    }
+
+    /// Get storage size in bytes
+    pub fn storage_bytes(&self) -> usize {
+        self.data.len() + self.absmax.len() * 4
+    }
+
+    /// Get original number of elements
+    pub fn num_elements(&self) -> usize {
+        self.num_elements
+    }
+
+    /// Get original shape
+    pub fn shape(&self) -> &[usize] {
+        &self.shape
+    }
+
+    /// Get packed data (for CUDA transfer)
+    pub fn packed_data(&self) -> &[u8] {
+        &self.data
+    }
+
+    /// Get absmax values (for CUDA transfer)
+    pub fn absmax(&self) -> &[f32] {
+        &self.absmax
+    }
+
+    /// Get block size
+    pub fn block_size(&self) -> usize {
+        self.block_size
+    }
+}
+
+/// Quantize a normalized value [-1, 1] to 4-bit NF4
+fn quantize_nf4(x: f32) -> u8 {
+    // Binary search tree for NF4 quantization
+    if x > 0.03979014977812767 {
+        if x > 0.3893125355243683 {
+            if x > 0.6427869200706482 {
+                if x > 0.8614784181118011 { 0b1111 }
+                else { 0b1110 }
+            } else {
+                if x > 0.5016634166240692 { 0b1101 }
+                else { 0b1100 }
+            }
+        } else {
+            if x > 0.2035212516784668 {
+                if x > 0.2920137718319893 { 0b1011 }
+                else { 0b1010 }
+            } else {
+                if x > 0.1202552504837513 { 0b1001 }
+                else { 0b1000 }
+            }
+        }
+    } else {
+        if x > -0.33967943489551544 {
+            if x > -0.13791173323988914 {
+                if x > -0.045525018125772476 { 0b0111 }
+                else { 0b0110 }
+            } else {
+                if x > -0.23460740596055984 { 0b0101 }
+                else { 0b0100 }
+            }
+        } else {
+            if x > -0.6106329262256622 {
+                if x > -0.4599952697753906 { 0b0011 }
+                else { 0b0010 }
+            } else {
+                if x > -0.8480964004993439 { 0b0001 }
+                else { 0b0000 }
+            }
+        }
+    }
+}
+
+/// Dequantize a single 4-bit NF4 value
+#[inline]
+pub fn dequantize_nf4(val: u8) -> f32 {
+    NF4_DEQUANT_TABLE[(val & 0x0F) as usize]
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    // Test 1: NF4 table has correct values
+    #[test]
+    fn test_nf4_table_bounds() {
+        assert_eq!(NF4_DEQUANT_TABLE[0], -1.0);
+        assert_eq!(NF4_DEQUANT_TABLE[7], 0.0);
+        assert_eq!(NF4_DEQUANT_TABLE[15], 1.0);
+    }
+
+    // Test 2: Quantization preserves zeros
+    #[test]
+    fn test_quantize_zero() {
+        let q = quantize_nf4(0.0);
+        assert_eq!(q, 0b0111); // 0.0 maps to index 7
+    }
+
+    // Test 3: Quantization preserves extremes
+    #[test]
+    fn test_quantize_extremes() {
+        let q_neg = quantize_nf4(-1.0);
+        let q_pos = quantize_nf4(1.0);
+        assert_eq!(q_neg, 0b0000);
+        assert_eq!(q_pos, 0b1111);
+    }
+
+    // Test 4: Dequantize inverts quantize approximately
+    #[test]
+    fn test_roundtrip_single() {
+        for i in 0..16u8 {
+            let val = NF4_DEQUANT_TABLE[i as usize];
+            let q = quantize_nf4(val);
+            let dq = dequantize_nf4(q);
+            // Should be exact for table values
+            assert!((dq - val).abs() < 0.01, "i={}, val={}, q={}, dq={}", i, val, q, dq);
+        }
+    }
+
+    // Test 5: Storage creation works
+    #[test]
+    fn test_storage_creation() {
+        let data: Vec<f32> = (0..128).map(|i| (i as f32 - 64.0) / 64.0).collect();
+        let storage = Nf4Storage::from_f32_slice(&data, 128, vec![128]).unwrap();
+
+        assert_eq!(storage.num_elements(), 128);
+        assert_eq!(storage.shape(), &[128]);
+        assert_eq!(storage.block_size(), 64);
+        // 128 elements = 64 bytes packed
+        assert_eq!(storage.packed_data().len(), 64);
+        // 128 / 64 = 2 blocks
+        assert_eq!(storage.absmax().len(), 2);
+    }
+
+    // Test 6: CPU dequantization produces correct shape
+    #[test]
+    fn test_cpu_dequant_shape() {
+        let data: Vec<f32> = vec![0.5; 256];
+        let storage = Nf4Storage::from_f32_slice(&data, 256, vec![16, 16]).unwrap();
+        let tensor = storage.to_tensor(&Device::Cpu).unwrap();
+
+        assert_eq!(tensor.dims(), &[16, 16]);
+    }
+
+    // Test 7: Roundtrip tensor preserves approximate values
+    #[test]
+    fn test_tensor_roundtrip() {
+        let device = Device::Cpu;
+        let original = Tensor::randn(0.0f32, 0.5, (64,), &device).unwrap();
+
+        let storage = Nf4Storage::from_tensor(&original).unwrap();
+        let recovered = storage.to_tensor(&device).unwrap();
+
+        let diff = (&original - &recovered).unwrap();
+        let mse: f32 = diff.sqr().unwrap().mean_all().unwrap().to_scalar().unwrap();
+        let rmse = mse.sqrt();
+
+        // NF4 should have reasonable error for normal distribution
+        assert!(rmse < 0.2, "RMSE {} too large", rmse);
+    }
+
+    // Test 8: Storage bytes are correct
+    #[test]
+    fn test_storage_bytes() {
+        let data: Vec<f32> = vec![0.0; 128];
+        let storage = Nf4Storage::from_f32_slice(&data, 128, vec![128]).unwrap();
+
+        // 64 bytes packed + 2 blocks * 4 bytes absmax = 72 bytes
+        assert_eq!(storage.storage_bytes(), 64 + 2 * 4);
+    }
+
+    // Test 9: Compression ratio is approximately 3-8x depending on block overhead
+    #[test]
+    fn test_compression_ratio() {
+        let num_elements = 1024;
+        let data: Vec<f32> = vec![0.0; num_elements];
+        let storage = Nf4Storage::from_f32_slice(&data, num_elements, vec![num_elements]).unwrap();
+
+        let original_bytes = num_elements * 4; // f32
+        let compressed_bytes = storage.storage_bytes();
+        let ratio = original_bytes as f32 / compressed_bytes as f32;
+
+        // 4-bit = 0.5 bytes per element + absmax overhead
+        // 1024 elements = 512 bytes packed + 16 blocks * 4 bytes = 576 bytes
+        // 4096 / 576 = ~7.1x compression
+        assert!(ratio > 2.0 && ratio < 10.0, "Compression ratio {} unexpected", ratio);
+    }
+
+    // Test 10: Works with non-power-of-2 sizes
+    #[test]
+    fn test_non_power_of_2() {
+        let data: Vec<f32> = vec![0.1; 100];
+        let storage = Nf4Storage::from_f32_slice(&data, 100, vec![100]).unwrap();
+        let tensor = storage.to_tensor(&Device::Cpu).unwrap();
+
+        assert_eq!(tensor.dims(), &[100]);
+    }
+}
diff --git a/candle-kernels/src/nf4.cu b/candle-kernels/src/nf4.cu
new file mode 100644
index 00000000..165de65d
--- /dev/null
+++ b/candle-kernels/src/nf4.cu
@@ -0,0 +1,231 @@
+// NF4 (4-bit NormalFloat) quantization kernels
+// Inspired by bitsandbytes, simplified for candle integration
+
+#include "cuda_fp16.h"
+#include "cuda_bf16.h"
+#include <stdint.h>
+
+// NF4 dequantization lookup table
+// These are the 16 quantile values for a normal distribution
+__device__ __constant__ float NF4_DEQUANT_TABLE[16] = {
+    -1.0f,
+    -0.6961928009986877f,
+    -0.5250730514526367f,
+    -0.39491748809814453f,
+    -0.28444138169288635f,
+    -0.18477343022823334f,
+    -0.09105003625154495f,
+    0.0f,
+    0.07958029955625534f,
+    0.16093020141124725f,
+    0.24611230194568634f,
+    0.33791524171829224f,
+    0.44070982933044434f,
+    0.5626170039176941f,
+    0.7229568362236023f,
+    1.0f
+};
+
+// Dequantize a single 4-bit value
+__device__ __forceinline__ float dequantize_nf4(unsigned char val) {
+    return NF4_DEQUANT_TABLE[val & 0x0F];
+}
+
+// ============================================================================
+// Dequantization kernel - converts NF4 packed data to float/half/bf16
+// ============================================================================
+// Each byte contains two 4-bit values (high nibble, low nibble)
+// absmax is per block (default block size = 64 elements)
+
+template<typename T>
+__global__ void nf4_dequantize_kernel(
+    const uint8_t* __restrict__ quantized,  // packed 4-bit data
+    const float* __restrict__ absmax,       // scale per block
+    T* __restrict__ output,                 // dequantized output
+    const int n,                            // total number of output elements
+    const int block_size                    // elements per block (default 64)
+) {
+    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
+    const int total_threads = idx * 2;  // each thread handles 2 elements (1 byte)
+
+    if (total_threads >= n) return;
+
+    // Load one byte (two 4-bit values)
+    const uint8_t packed = quantized[idx];
+
+    // Get block index for absmax lookup
+    const int block_idx = total_threads / block_size;
+    const float scale = absmax[block_idx];
+
+    // Dequantize high nibble (first element)
+    const float val0 = dequantize_nf4(packed >> 4) * scale;
+    output[total_threads] = (T)val0;
+
+    // Dequantize low nibble (second element) if in bounds
+    if (total_threads + 1 < n) {
+        const float val1 = dequantize_nf4(packed & 0x0F) * scale;
+        output[total_threads + 1] = (T)val1;
+    }
+}
+
+// ============================================================================
+// Fused dequantize + matmul kernel (for better performance)
+// ============================================================================
+// This avoids materializing the full dequantized weight matrix
+
+// For now, we provide the simple dequantize kernel
+// Fused kernel can be added later for even better perf
+
+// ============================================================================
+// Quantization kernel - converts float to NF4
+// ============================================================================
+
+__device__ __forceinline__ unsigned char quantize_nf4(float x) {
+    // Binary search tree for NF4 quantization
+    // Generated from normal distribution quantiles
+    if (x > 0.03979014977812767f) {
+        if (x > 0.3893125355243683f) {
+            if (x > 0.6427869200706482f) {
+                if (x > 0.8614784181118011f) return 0b1111;
+                else return 0b1110;
+            } else {
+                if (x > 0.5016634166240692f) return 0b1101;
+                else return 0b1100;
+            }
+        } else {
+            if (x > 0.2035212516784668f) {
+                if (x > 0.2920137718319893f) return 0b1011;
+                else return 0b1010;
+            } else {
+                if (x > 0.1202552504837513f) return 0b1001;
+                else return 0b1000;
+            }
+        }
+    } else {
+        if (x > -0.33967943489551544f) {
+            if (x > -0.13791173323988914f) {
+                if (x > -0.045525018125772476f) return 0b0111;
+                else return 0b0110;
+            } else {
+                if (x > -0.23460740596055984f) return 0b0101;
+                else return 0b0100;
+            }
+        } else {
+            if (x > -0.6106329262256622f) {
+                if (x > -0.4599952697753906f) return 0b0011;
+                else return 0b0010;
+            } else {
+                if (x > -0.8480964004993439f) return 0b0001;
+                else return 0b0000;
+            }
+        }
+    }
+}
+
+template<typename T>
+__global__ void nf4_quantize_kernel(
+    const T* __restrict__ input,            // input data
+    uint8_t* __restrict__ quantized,        // packed 4-bit output
+    float* __restrict__ absmax,             // scale per block (output)
+    const int n,                            // total number of input elements
+    const int block_size                    // elements per block (default 64)
+) {
+    // This kernel handles quantization in two passes:
+    // Pass 1: Find absmax per block (reduction)
+    // Pass 2: Quantize using absmax
+
+    // For simplicity, this kernel does pass 2 assuming absmax is precomputed
+    // A separate kernel or CPU code handles absmax computation
+
+    const int idx = blockIdx.x * blockDim.x + threadIdx.x;
+    const int elem_idx = idx * 2;  // each thread handles 2 elements (1 byte output)
+
+    if (elem_idx >= n) return;
+
+    // Get block absmax
+    const int block_idx = elem_idx / block_size;
+    const float scale = absmax[block_idx];
+    const float inv_scale = (scale != 0.0f) ? 1.0f / scale : 0.0f;
+
+    // Quantize two values and pack into one byte
+    float val0 = (float)input[elem_idx] * inv_scale;
+    val0 = fmaxf(-1.0f, fminf(1.0f, val0));  // clamp to [-1, 1]
+    unsigned char q0 = quantize_nf4(val0);
+
+    unsigned char q1 = 0;
+    if (elem_idx + 1 < n) {
+        float val1 = (float)input[elem_idx + 1] * inv_scale;
+        val1 = fmaxf(-1.0f, fminf(1.0f, val1));
+        q1 = quantize_nf4(val1);
+    }
+
+    // Pack: high nibble = first element, low nibble = second element
+    quantized[idx] = (q0 << 4) | q1;
+}
+
+// ============================================================================
+// Absmax computation kernel
+// ============================================================================
+
+__global__ void nf4_compute_absmax_kernel(
+    const float* __restrict__ input,
+    float* __restrict__ absmax,
+    const int n,
+    const int block_size
+) {
+    extern __shared__ float sdata[];
+
+    const int block_idx = blockIdx.x;
+    const int tid = threadIdx.x;
+    const int block_start = block_idx * block_size;
+
+    // Load and find local max
+    float local_max = 0.0f;
+    for (int i = tid; i < block_size && (block_start + i) < n; i += blockDim.x) {
+        float val = fabsf(input[block_start + i]);
+        local_max = fmaxf(local_max, val);
+    }
+
+    sdata[tid] = local_max;
+    __syncthreads();
+
+    // Reduction in shared memory
+    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
+        if (tid < s) {
+            sdata[tid] = fmaxf(sdata[tid], sdata[tid + s]);
+        }
+        __syncthreads();
+    }
+
+    if (tid == 0) {
+        absmax[block_idx] = sdata[0];
+    }
+}
+
+// ============================================================================
+// Explicit template instantiations for PTX generation
+// ============================================================================
+
+template __global__ void nf4_dequantize_kernel<float>(
+    const uint8_t* quantized, const float* absmax, float* output, int n, int block_size
+);
+
+template __global__ void nf4_dequantize_kernel<__half>(
+    const uint8_t* quantized, const float* absmax, __half* output, int n, int block_size
+);
+
+template __global__ void nf4_dequantize_kernel<__nv_bfloat16>(
+    const uint8_t* quantized, const float* absmax, __nv_bfloat16* output, int n, int block_size
+);
+
+template __global__ void nf4_quantize_kernel<float>(
+    const float* input, uint8_t* quantized, float* absmax, int n, int block_size
+);
+
+template __global__ void nf4_quantize_kernel<__half>(
+    const __half* input, uint8_t* quantized, float* absmax, int n, int block_size
+);
+
+template __global__ void nf4_quantize_kernel<__nv_bfloat16>(
+    const __nv_bfloat16* input, uint8_t* quantized, float* absmax, int n, int block_size
+);
