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

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>
#include <vector>
#include <algorithm>

struct Stats {
    float mean;
    float m2;
    int count;
};

__device__ __forceinline__ Stats merge_stats(Stats a, Stats b) {
    if (a.count == 0) return b;
    if (b.count == 0) return a;
    const int count = a.count + b.count;
    const float delta = b.mean - a.mean;
    const float ratio = float(b.count) / float(count);
    a.mean += delta * ratio;
    a.m2 += b.m2 + delta * delta * (float(a.count) * ratio);
    a.count = count;
    return a;
}

__device__ __forceinline__ void update_stats(Stats& s, float value) {
    ++s.count;
    const float delta = value - s.mean;
    s.mean += delta / float(s.count);
    s.m2 += delta * (value - s.mean);
}

__device__ __forceinline__ Stats warp_reduce(Stats s) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        Stats other;
        other.mean = __shfl_down_sync(0xffffffff, s.mean, offset);
        other.m2 = __shfl_down_sync(0xffffffff, s.m2, offset);
        other.count = __shfl_down_sync(0xffffffff, s.count, offset);
        if ((threadIdx.x & 31) + offset < 32)
            s = merge_stats(s, other);
    }
    return s;
}

__device__ __forceinline__ Stats block_reduce(Stats s, Stats* shared) {
    s = warp_reduce(s);
    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;
    if (lane == 0)
        shared[warp] = s;
    __syncthreads();

    Stats result = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8)
            result = shared[lane];
        result = warp_reduce(result);
    }
    return result;
}

__global__ void partial_stats_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial_mean,
    float* __restrict__ partial_m2,
    int n,
    int chunks
) {
    __shared__ Stats shared[8];
    const int row = blockIdx.x / chunks;
    const int chunk = blockIdx.x % chunks;
    const int first = chunk * 1024;
    const int vectors = n / 4;

    const float4* row_x =
        reinterpret_cast<const float4*>(x + int64_t(row) * n);

    Stats s = {0.0f, 0.0f, 0};
    #pragma unroll
    for (int step = 0; step < 4; ++step) {
        const int i = first + threadIdx.x + step * 256;
        if (i < vectors) {
            const float4 v = row_x[i];
            update_stats(s, v.x);
            update_stats(s, v.y);
            update_stats(s, v.z);
            update_stats(s, v.w);
        }
    }

    s = block_reduce(s, shared);
    if (threadIdx.x == 0) {
        partial_mean[blockIdx.x] = s.mean;
        partial_m2[blockIdx.x] = s.m2;
    }
}

__global__ void final_stats_kernel(
    const float* __restrict__ partial_mean,
    const float* __restrict__ partial_m2,
    float* __restrict__ mean,
    float* __restrict__ rstd,
    int n,
    int chunks,
    float eps
) {
    __shared__ Stats shared[8];
    const int row = blockIdx.x;
    Stats s = {0.0f, 0.0f, 0};

    for (int chunk = threadIdx.x; chunk < chunks; chunk += 256) {
        const int remaining = n - chunk * 4096;
        Stats other;
        other.mean = partial_mean[int64_t(row) * chunks + chunk];
        other.m2 = partial_m2[int64_t(row) * chunks + chunk];
        other.count = remaining < 4096 ? remaining : 4096;
        s = merge_stats(s, other);
    }

    s = block_reduce(s, shared);
    if (threadIdx.x == 0) {
        mean[row] = s.mean;
        rstd[row] = rsqrtf(fmaxf(s.m2 / float(n), 0.0f) + eps);
    }
}

__global__ void normalize_kernel(
    const float4* __restrict__ x,
    const float4* __restrict__ weight,
    const float4* __restrict__ bias,
    const float* __restrict__ mean,
    const float* __restrict__ rstd,
    float4* __restrict__ output,
    int vectors_per_row,
    int64_t total_vectors
) {
    for (int64_t i = int64_t(blockIdx.x) * blockDim.x + threadIdx.x;
         i < total_vectors;
         i += int64_t(blockDim.x) * gridDim.x) {
        const int row = i / vectors_per_row;
        const int column = i % vectors_per_row;

        const float m = mean[row];
        const float r = rstd[row];
        const float4 v = x[i];
        const float4 w = weight[column];
        const float4 b = bias[column];

        float4 y;
        y.x = ((v.x - m) * r) * w.x + b.x;
        y.y = ((v.y - m) * r) * w.y + b.y;
        y.z = ((v.z - m) * r) * w.z + b.z;
        y.w = ((v.w - m) * r) * w.w + b.w;
        output[i] = y;
    }
}

std::vector<torch::Tensor> large_layer_norm_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    double eps
) {
    c10::cuda::CUDAGuard guard(x.device());

    const int n = static_cast<int>(weight.numel());
    const int rows = static_cast<int>(x.numel() / n);
    const int chunks = (n - 1) / 4096 + 1;

    auto output = torch::empty_like(x);
    auto partial = torch::empty({2, rows, chunks}, x.options());
    auto mean = torch::empty({rows}, x.options());
    auto rstd = torch::empty({rows}, x.options());

    float* partial_mean = partial.data_ptr<float>();
    float* partial_m2 = partial_mean + int64_t(rows) * chunks;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    partial_stats_kernel<<<rows * chunks, 256, 0, stream>>>(
        x.data_ptr<float>(), partial_mean, partial_m2, n, chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    final_stats_kernel<<<rows, 256, 0, stream>>>(
        partial_mean, partial_m2,
        mean.data_ptr<float>(), rstd.data_ptr<float>(),
        n, chunks, static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    const int64_t total_vectors = x.numel() / 4;
    const int blocks = static_cast<int>(
        std::min<int64_t>((total_vectors + 255) / 256, 65535));

    normalize_kernel<<<blocks, 256, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        reinterpret_cast<const float4*>(weight.data_ptr<float>()),
        reinterpret_cast<const float4*>(bias.data_ptr<float>()),
        mean.data_ptr<float>(), rstd.data_ptr<float>(),
        reinterpret_cast<float4*>(output.data_ptr<float>()),
        n / 4, total_vectors);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return {output, mean, rstd};
}
"""

_extension = load_inline(
    name="large_layer_norm_welford_cuda",
    cpp_sources=r"""
#include <torch/extension.h>
#include <vector>
std::vector<torch::Tensor> large_layer_norm_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    double eps);
""",
    cuda_sources=_cuda_source,
    functions=["large_layer_norm_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class _LargeLayerNorm(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, eps):
        output, mean, rstd = _extension.large_layer_norm_cuda(
            x, weight, bias, eps
        )
        ctx.save_for_backward(x, weight, mean, rstd)
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, mean, rstd = ctx.saved_tensors
        n = weight.numel()
        grad = grad_output.reshape(-1, n)
        inv_std = rstd[:, None]

        dx = dw = db = None
        normalized = None
        if ctx.needs_input_grad[0] or ctx.needs_input_grad[1]:
            normalized = (x.reshape(-1, n) - mean[:, None]) * inv_std

        if ctx.needs_input_grad[0]:
            scaled_grad = grad * weight.reshape(1, n)
            dx = (
                scaled_grad
                - scaled_grad.mean(dim=1, keepdim=True)
                - normalized
                * (scaled_grad * normalized).mean(dim=1, keepdim=True)
            ) * inv_std
            dx = dx.reshape_as(x)

        if ctx.needs_input_grad[1]:
            dw = (grad * normalized).sum(dim=0).reshape_as(weight)

        if ctx.needs_input_grad[2]:
            db = grad.sum(dim=0).reshape_as(weight)

        return dx, dw, db, None


class ModelNew(nn.Module):
    def __init__(self, normalized_shape: tuple):
        super().__init__()
        self.ln = nn.LayerNorm(normalized_shape=normalized_shape)

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        weight = self.ln.weight
        bias = self.ln.bias
        shape = self.ln.normalized_shape
        n = weight.numel()

        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or weight.dtype != torch.float32
            or bias.dtype != torch.float32
            or weight.device != x.device
            or bias.device != x.device
            or not x.is_contiguous()
            or not weight.is_contiguous()
            or not bias.is_contiguous()
            or x.ndim < len(shape)
            or tuple(x.shape[-len(shape):]) != tuple(shape)
            or x.numel() == 0
            or n < 4096
            or n % 4 != 0
            or n > 2147483647
            or x.numel() // n * ((n - 1) // 4096 + 1) > 2147483647
            or x.data_ptr() % 16 != 0
            or weight.data_ptr() % 16 != 0
            or bias.data_ptr() % 16 != 0
        ):
            return self.ln(x)

        return _LargeLayerNorm.apply(x, weight, bias, self.ln.eps)
```