import torch
import torch.nn as nn

class Model(nn.Module):
    """
    Simple model that performs RMS Normalization.
    """
    def __init__(self, num_features: int, eps: float = 1e-5):
        """
        Initializes the RMSNorm layer.

        Args:
            num_features (int): Number of features in the input tensor.
            eps (float, optional): A small value added to the denominator to avoid division by zero. Defaults to 1e-5.
        """
        super(Model, self).__init__()
        self.num_features = num_features
        self.eps = eps

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        """
        Applies RMS Normalization to the input tensor.

        Args:
            x (torch.Tensor): Input tensor of shape (batch_size, num_features, *).

        Returns:
            torch.Tensor: Output tensor with RMS Normalization applied, same shape as input.
        """
        # Calculate the RMS along the feature dimension
        rms = torch.sqrt(torch.mean(x ** 2, dim=1, keepdim=True) + self.eps)

        # Normalize the input by dividing by the RMS
        return x / rms

batch_size = 112
features = 64
dim1 = 512
dim2 = 512

def get_inputs():
    x = torch.rand(batch_size, features, dim1, dim2)
    return [x]

def get_init_inputs():
    return [features]

import torch
import torch.nn as nn
from torch.utils.cpp_extension import load_inline

_cpp_src = r"""
#include <torch/extension.h>

void rmsnorm_forward_cuda(torch::Tensor x, torch::Tensor y, double eps);

torch::Tensor rmsnorm_forward(torch::Tensor x, double eps) {
    auto y = torch::empty_like(x);
    rmsnorm_forward_cuda(x, y, eps);
    return y;
}
"""

_cuda_src = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda.h>
#include <cuda_runtime.h>

__global__ void rmsnorm_c64_kernel(const float* __restrict__ x, float* __restrict__ y,
                                   int N, int S, float eps) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int total = N * S;
    if (col >= total) return;
    int n = col / S;
    int s = col - n * S;
    long long base = ((long long)n) * 64LL * (long long)S + s;

    float s0 = 0.0f, s1 = 0.0f, s2 = 0.0f, s3 = 0.0f;
    #pragma unroll
    for (int c = 0; c < 64; c += 4) {
        float v0 = x[base + (long long)(c + 0) * S];
        float v1 = x[base + (long long)(c + 1) * S];
        float v2 = x[base + (long long)(c + 2) * S];
        float v3 = x[base + (long long)(c + 3) * S];
        s0 = fmaf(v0, v0, s0);
        s1 = fmaf(v1, v1, s1);
        s2 = fmaf(v2, v2, s2);
        s3 = fmaf(v3, v3, s3);
    }
    float scale = rsqrtf(((s0 + s1) + (s2 + s3)) * (1.0f / 64.0f) + eps);
    #pragma unroll
    for (int c = 0; c < 64; c += 4) {
        long long i0 = base + (long long)(c + 0) * S;
        long long i1 = base + (long long)(c + 1) * S;
        long long i2 = base + (long long)(c + 2) * S;
        long long i3 = base + (long long)(c + 3) * S;
        y[i0] = x[i0] * scale;
        y[i1] = x[i1] * scale;
        y[i2] = x[i2] * scale;
        y[i3] = x[i3] * scale;
    }
}

__global__ void rmsnorm_generic_kernel(const float* __restrict__ x, float* __restrict__ y,
                                       int N, int C, int S, float eps) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int total = N * S;
    if (col >= total) return;
    int n = col / S;
    int s = col - n * S;
    long long base = ((long long)n) * C * (long long)S + s;
    float sum = 0.0f;
    for (int c = 0; c < C; ++c) {
        float v = x[base + (long long)c * S];
        sum = fmaf(v, v, sum);
    }
    float scale = rsqrtf(sum / (float)C + eps);
    for (int c = 0; c < C; ++c) {
        long long idx = base + (long long)c * S;
        y[idx] = x[idx] * scale;
    }
}

void rmsnorm_forward_cuda(torch::Tensor x, torch::Tensor y, double eps) {
    int N = (int)x.size(0);
    int C = (int)x.size(1);
    long long Sll = 1;
    for (int i = 2; i < x.dim(); ++i) Sll *= x.size(i);
    int S = (int)Sll;
    int total = N * S;
    const int threads = 256;
    int blocks = (total + threads - 1) / threads;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();
    if (C == 64) {
        rmsnorm_c64_kernel<<<blocks, threads, 0, stream>>>(x.data_ptr<float>(), y.data_ptr<float>(), N, S, (float)eps);
    } else {
        rmsnorm_generic_kernel<<<blocks, threads, 0, stream>>>(x.data_ptr<float>(), y.data_ptr<float>(), N, C, S, (float)eps);
    }
}
"""

_ext = load_inline(
    name="rmsnorm_ext_final",
    cpp_sources=[_cpp_src],
    cuda_sources=[_cuda_src],
    functions=["rmsnorm_forward"],
    extra_cuda_cflags=["-O3", "--use_fast_math"],
    extra_cflags=["-O3"],
    verbose=False,
)

class ModelNew(nn.Module):
    def __init__(self, num_features: int, eps: float = 1e-5):
        super(ModelNew, self).__init__()
        self.num_features = num_features
        self.eps = eps

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        return _ext.rmsnorm_forward(x, float(self.eps))