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


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

torch::Tensor padded_channels_last(torch::Tensor x, int64_t pd,
                                   int64_t ph, int64_t pw);
torch::Tensor reverse_weight(torch::Tensor weight);
torch::Tensor relu_group_norm(torch::Tensor x, torch::Tensor weight,
                              torch::Tensor bias, double eps);
"""

_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>

namespace {

constexpr int THREADS = 256;
constexpr int CHANNELS = 128;
constexpr int GROUPS = 8;
constexpr int VECTORS = 32;
constexpr int SPATIAL_TILE = 128;
constexpr int OUTPUT_TILE = 2048;

__device__ __forceinline__ float warp_sum(float x) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        x += __shfl_down_sync(0xffffffff, x, offset);
    return x;
}

__device__ __forceinline__ float4 relu4(float4 x) {
    x.x = x.x < 0.0f ? 0.0f : x.x;
    x.y = x.y < 0.0f ? 0.0f : x.y;
    x.z = x.z < 0.0f ? 0.0f : x.z;
    x.w = x.w < 0.0f ? 0.0f : x.w;
    return x;
}

// Fuse NCDHW -> NDHWC conversion with the halo required by the
// equivalent forward convolution.
__global__ __launch_bounds__(THREADS)
void padded_channels_last_kernel(
    const float* __restrict__ input,
    float* __restrict__ output,
    int depth, int height, int width,
    int padded_height, int padded_width,
    int pd, int ph, int pw,
    int input_spatial, int output_spatial)
{
    __shared__ float tile[64][33];

    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;
    const int start = blockIdx.x * 32;
    const int s = start + lane;

    const int ox = s % padded_width;
    const int q = s / padded_width;
    const int oy = q % padded_height;
    const int oz = q / padded_height;

    const int ix = ox - pw;
    const int iy = oy - ph;
    const int iz = oz - pd;

    const bool valid =
        s < output_spatial &&
        static_cast<unsigned int>(ix) < static_cast<unsigned int>(width) &&
        static_cast<unsigned int>(iy) < static_cast<unsigned int>(height) &&
        static_cast<unsigned int>(iz) < static_cast<unsigned int>(depth);

    const int source = (iz * height + iy) * width + ix;
    const int64_t input_base =
        static_cast<int64_t>(blockIdx.y) * 64 * input_spatial;
    const int64_t output_base =
        static_cast<int64_t>(blockIdx.y) * 64 * output_spatial;

    #pragma unroll
    for (int j = 0; j < 8; ++j) {
        const int c = warp + j * 8;
        tile[c][lane] = valid
            ? input[input_base + static_cast<int64_t>(c) * input_spatial
                    + source]
            : 0.0f;
    }
    __syncthreads();

    #pragma unroll
    for (int j = 0; j < 4; ++j) {
        const int row = warp + j * 8;
        if (start + row < output_spatial) {
            const int64_t index =
                output_base + static_cast<int64_t>(start + row) * 64 + lane;
            output[index] = tile[lane][row];
            output[index + 32] = tile[lane + 32][row];
        }
    }
}

__global__ __launch_bounds__(THREADS)
void reverse_weight_kernel(
    const float* __restrict__ input,
    float* __restrict__ output,
    int kernel_volume)
{
    __shared__ float tile[32][33];

    const int lane = threadIdx.x & 31;
    const int warp = threadIdx.x >> 5;
    const int k = blockIdx.x;
    const int input_start = blockIdx.y * 32;
    const int output_start = blockIdx.z * 32;

    #pragma unroll
    for (int j = 0; j < 4; ++j) {
        const int ic = input_start + warp + j * 8;
        const int oc = output_start + lane;
        tile[warp + j * 8][lane] =
            input[(static_cast<int64_t>(ic) * kernel_volume + k) * 128 + oc];
    }
    __syncthreads();

    const int reversed_k = kernel_volume - 1 - k;

    #pragma unroll
    for (int j = 0; j < 4; ++j) {
        const int oc = output_start + warp + j * 8;
        const int ic = input_start + lane;
        output[(static_cast<int64_t>(oc) * kernel_volume + reversed_k) * 64
               + ic] = tile[lane][warp + j * 8];
    }
}

__global__ __launch_bounds__(THREADS)
void partial_stats_kernel(
    const float4* __restrict__ x,
    float2* __restrict__ partial,
    int spatial,
    int chunks)
{
    __shared__ float shared_sum[64];
    __shared__ float shared_square[64];

    const int tid = threadIdx.x;
    const int channel_vector = tid & 31;
    const int spatial_lane = tid >> 5;
    const int n = blockIdx.y;
    const int chunk = blockIdx.x;
    const int64_t base = static_cast<int64_t>(n) * spatial * VECTORS;

    float sum = 0.0f;
    float square = 0.0f;

    #pragma unroll
    for (int j = 0; j < SPATIAL_TILE / 8; ++j) {
        const int s = chunk * SPATIAL_TILE + spatial_lane + j * 8;
        if (s < spatial) {
            const float4 v = relu4(
                x[base + static_cast<int64_t>(s) * VECTORS + channel_vector]);
            sum += (v.x + v.y) + (v.z + v.w);
            square += (v.x * v.x + v.y * v.y)
                    + (v.z * v.z + v.w * v.w);
        }
    }

    #pragma unroll
    for (int offset = 2; offset > 0; offset >>= 1) {
        sum += __shfl_down_sync(0xffffffff, sum, offset, 4);
        square += __shfl_down_sync(0xffffffff, square, offset, 4);
    }

    if ((tid & 3) == 0) {
        shared_sum[tid >> 2] = sum;
        shared_square[tid >> 2] = square;
    }
    __syncthreads();

    if (tid < 64) {
        const int group = tid >> 3;
        const int lane = tid & 7;

        sum = shared_sum[lane * GROUPS + group];
        square = shared_square[lane * GROUPS + group];

        #pragma unroll
        for (int offset = 4; offset > 0; offset >>= 1) {
            sum += __shfl_down_sync(0xffffffff, sum, offset, 8);
            square += __shfl_down_sync(0xffffffff, square, offset, 8);
        }

        if (lane == 0)
            partial[(static_cast<int64_t>(n) * GROUPS + group) * chunks
                    + chunk] = make_float2(sum, square);
    }
}

__global__ __launch_bounds__(THREADS)
void finish_stats_kernel(
    const float2* __restrict__ partial,
    const float* __restrict__ weight,
    const float* __restrict__ bias,
    float* __restrict__ scale,
    float* __restrict__ shift,
    int chunks,
    int spatial,
    float eps)
{
    __shared__ float shared_sum[8];
    __shared__ float shared_square[8];

    const int ng = blockIdx.x;
    const int tid = threadIdx.x;
    const int lane = tid & 31;
    const int warp = tid >> 5;

    float sum = 0.0f;
    float square = 0.0f;

    for (int i = tid; i < chunks; i += THREADS) {
        const float2 p = partial[static_cast<int64_t>(ng) * chunks + i];
        sum += p.x;
        square += p.y;
    }

    sum = warp_sum(sum);
    square = warp_sum(square);

    if (lane == 0) {
        shared_sum[warp] = sum;
        shared_square[warp] = square;
    }
    __syncthreads();

    if (warp == 0) {
        sum = warp_sum(lane < 8 ? shared_sum[lane] : 0.0f);
        square = warp_sum(lane < 8 ? shared_square[lane] : 0.0f);
        sum = __shfl_sync(0xffffffff, sum, 0);
        square = __shfl_sync(0xffffffff, square, 0);

        if (lane < 16) {
            const float reciprocal =
                1.0f / (static_cast<float>(spatial) * 16.0f);
            const float mean = sum * reciprocal;
            const float variance =
                fmaxf(square * reciprocal - mean * mean, 0.0f);
            const float rstd = rsqrtf(variance + eps);
            const int channel = (ng % GROUPS) * 16 + lane;
            const float a = rstd * weight[channel];

            scale[ng * 16 + lane] = a;
            shift[ng * 16 + lane] = fmaf(-mean, a, bias[channel]);
        }
    }
}

__global__ __launch_bounds__(THREADS)
void normalize_kernel(
    float4* __restrict__ x,
    const float4* __restrict__ scale,
    const float4* __restrict__ shift,
    int vectors_per_batch)
{
    const int tid = threadIdx.x;
    const int n = blockIdx.y;
    const int channel_vector = tid & 31;

    const float4 a = scale[n * VECTORS + channel_vector];
    const float4 b = shift[n * VECTORS + channel_vector];
    const int64_t base = static_cast<int64_t>(n) * vectors_per_batch;
    const int start = blockIdx.x * OUTPUT_TILE + tid;

    #pragma unroll
    for (int j = 0; j < OUTPUT_TILE / THREADS; ++j) {
        const int i = start + j * THREADS;
        if (i < vectors_per_batch) {
            float4 v = relu4(x[base + i]);
            v.x = fmaf(v.x, a.x, b.x);
            v.y = fmaf(v.y, a.y, b.y);
            v.z = fmaf(v.z, a.z, b.z);
            v.w = fmaf(v.w, a.w, b.w);
            x[base + i] = v;
        }
    }
}

} // namespace

torch::Tensor padded_channels_last(
    torch::Tensor x, int64_t pd, int64_t ph, int64_t pw)
{
    TORCH_CHECK(
        x.is_cuda() && x.scalar_type() == at::kFloat &&
        x.dim() == 5 && x.size(1) == 64 && x.is_contiguous(),
        "Expected contiguous CUDA float32 input with 64 channels");
    TORCH_CHECK(pd >= 0 && ph >= 0 && pw >= 0, "Invalid padding");

    const c10::cuda::CUDAGuard guard(x.device());

    const int64_t od = x.size(2) + 2 * pd;
    const int64_t oh = x.size(3) + 2 * ph;
    const int64_t ow = x.size(4) + 2 * pw;
    const int64_t input_spatial = x.size(2) * x.size(3) * x.size(4);
    const int64_t output_spatial = od * oh * ow;

    TORCH_CHECK(
        x.size(0) > 0 && x.size(0) <= 65535 &&
        output_spatial > 0 && output_spatial < 2147483647,
        "Unsupported input size");

    auto output = torch::empty(
        {x.size(0), 64, od, oh, ow},
        x.options().memory_format(at::MemoryFormat::ChannelsLast3d));

    const cudaStream_t stream = at::cuda::getCurrentCUDAStream();
    padded_channels_last_kernel<<<
        dim3((output_spatial + 31) / 32,
             static_cast<unsigned int>(x.size(0))),
        THREADS, 0, stream>>>(
        x.data_ptr<float>(), output.data_ptr<float>(),
        x.size(2), x.size(3), x.size(4), oh, ow,
        pd, ph, pw, input_spatial, output_spatial);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}

torch::Tensor reverse_weight(torch::Tensor weight)
{
    TORCH_CHECK(
        weight.is_cuda() && weight.scalar_type() == at::kFloat &&
        weight.dim() == 5 && weight.size(0) == 64 &&
        weight.size(1) == 128 &&
        weight.is_contiguous(at::MemoryFormat::ChannelsLast3d),
        "Expected channels-last CUDA float32 convolution weight");

    const c10::cuda::CUDAGuard guard(weight.device());

    auto output = torch::empty(
        {128, 64, weight.size(2), weight.size(3), weight.size(4)},
        weight.options().memory_format(at::MemoryFormat::ChannelsLast3d));

    const int volume = weight.size(2) * weight.size(3) * weight.size(4);
    const cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    reverse_weight_kernel<<<dim3(volume, 2, 4), THREADS, 0, stream>>>(
        weight.data_ptr<float>(), output.data_ptr<float>(), volume);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}

torch::Tensor relu_group_norm(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor bias,
    double eps)
{
    TORCH_CHECK(
        x.is_cuda() && x.scalar_type() == at::kFloat &&
        x.dim() == 5 && x.size(1) == CHANNELS &&
        x.is_contiguous(at::MemoryFormat::ChannelsLast3d),
        "Expected channels-last CUDA float32 input with 128 channels");
    TORCH_CHECK(
        weight.device() == x.device() && bias.device() == x.device() &&
        weight.scalar_type() == at::kFloat &&
        bias.scalar_type() == at::kFloat &&
        weight.numel() == CHANNELS && bias.numel() == CHANNELS &&
        weight.is_contiguous() && bias.is_contiguous(),
        "Invalid normalization parameters");

    const c10::cuda::CUDAGuard guard(x.device());

    const int batch = x.size(0);
    const int spatial = x.size(2) * x.size(3) * x.size(4);
    if (batch == 0 || spatial == 0)
        return x;

    const int chunks = (spatial + SPATIAL_TILE - 1) / SPATIAL_TILE;
    const int vectors_per_batch = spatial * VECTORS;

    auto partial = torch::empty({batch * GROUPS, chunks, 2}, x.options());
    auto affine = torch::empty({2, batch, CHANNELS}, x.options());

    float* scale = affine.data_ptr<float>();
    float* shift = scale + static_cast<int64_t>(batch) * CHANNELS;
    const cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    partial_stats_kernel<<<dim3(chunks, batch), THREADS, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        spatial, chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_stats_kernel<<<batch * GROUPS, THREADS, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        weight.data_ptr<float>(), bias.data_ptr<float>(),
        scale, shift, chunks, spatial, static_cast<float>(eps));
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    normalize_kernel<<<
        dim3((vectors_per_batch + OUTPUT_TILE - 1) / OUTPUT_TILE, batch),
        THREADS, 0, stream>>>(
        reinterpret_cast<float4*>(x.data_ptr<float>()),
        reinterpret_cast<const float4*>(scale),
        reinterpret_cast<const float4*>(shift),
        vectors_per_batch);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return x;
}
"""

_ops = load_inline(
    name="padded_conv3d_fused_relu_groupnorm_v6",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["padded_channels_last", "reverse_weight", "relu_group_norm"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(
        self, in_channels, out_channels, kernel_size, groups, bias=False
    ):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size, bias=bias
        ).to(memory_format=torch.channels_last_3d)
        self.relu = nn.ReLU()
        self.group_norm = nn.GroupNorm(groups, out_channels)

    def forward(self, x):
        conv = self.conv_transpose
        norm = self.group_norm

        fast = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.ndim == 5
            and x.shape[1] == 64
            and conv.in_channels == 64
            and conv.out_channels == 128
            and conv.groups == 1
            and conv.stride == (1, 1, 1)
            and conv.padding == (0, 0, 0)
            and conv.output_padding == (0, 0, 0)
            and conv.dilation == (1, 1, 1)
            and conv.weight.dtype == torch.float32
            and conv.weight.device == x.device
            and norm.num_groups == 8
            and norm.weight.dtype == torch.float32
            and norm.bias.dtype == torch.float32
            and norm.weight.device == x.device
            and norm.bias.device == x.device
            and 0 < x.shape[0] <= 65535
            and all(s > 0 for s in x.shape[2:])
        )

        if not fast:
            return norm(self.relu(conv(x)))

        padding = tuple(k - 1 for k in conv.kernel_size)
        padded_spatial = (
            (x.shape[2] + 2 * padding[0])
            * (x.shape[3] + 2 * padding[1])
            * (x.shape[4] + 2 * padding[2])
        )
        output_spatial = (
            (x.shape[2] + padding[0])
            * (x.shape[3] + padding[1])
            * (x.shape[4] + padding[2])
        )

        if padded_spatial >= 2**31 or output_spatial * 128 >= 2**31:
            return norm(self.relu(conv(x)))

        weight = _ops.reverse_weight(
            conv.weight.contiguous(memory_format=torch.channels_last_3d)
        )

        if x.is_contiguous():
            x = _ops.padded_channels_last(x, *padding)
            x = F.conv3d(x, weight, conv.bias, stride=1, padding=0)
        else:
            x = x.contiguous(memory_format=torch.channels_last_3d)
            x = F.conv3d(x, weight, conv.bias, stride=1, padding=padding)

        x = x.contiguous(memory_format=torch.channels_last_3d)

        if x.dtype != torch.float32:
            return norm(self.relu(x))

        return _ops.relu_group_norm(
            x, norm.weight.contiguous(), norm.bias.contiguous(), norm.eps
        )
```