```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"""
torch::Tensor centered_bn32_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor conv_bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    double eps,
    double momentum,
    bool training);
"""

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

constexpr int CHANNELS = 32;
constexpr int TILE = 2048;
constexpr int OUTPUT_VECTORS = 4096;

__device__ __forceinline__ float4 add4(float4 a, float4 b) {
    return make_float4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w);
}

__device__ __forceinline__ float4 shuffle_down4(float4 a, int offset) {
    return make_float4(
        __shfl_down_sync(0xffffffff, a.x, offset),
        __shfl_down_sync(0xffffffff, a.y, offset),
        __shfl_down_sync(0xffffffff, a.z, offset),
        __shfl_down_sync(0xffffffff, a.w, offset));
}

__global__ void partial_stats32(
    const float4* __restrict__ x,
    float4* __restrict__ sums,
    float4* __restrict__ squares,
    int spatial,
    int parts)
{
    const int group = threadIdx.x & 7;
    const int row = threadIdx.x >> 3;
    const int warp = threadIdx.x >> 5;
    const int n = blockIdx.y;
    const int p = blockIdx.x;
    const int begin = p * TILE;
    const int end = min(begin + TILE, spatial);
    const int64_t base = int64_t(n) * spatial * 8;

    float4 sum = make_float4(0.f, 0.f, 0.f, 0.f);
    float4 square = make_float4(0.f, 0.f, 0.f, 0.f);

    for (int start = begin + row; start < end; start += 128) {
        #pragma unroll
        for (int j = 0; j < 4; ++j) {
            const int s = start + j * 32;
            if (s < end) {
                float4 v = x[base + int64_t(s) * 8 + group];
                sum = add4(sum, v);
                square.x = fmaf(v.x, v.x, square.x);
                square.y = fmaf(v.y, v.y, square.y);
                square.z = fmaf(v.z, v.z, square.z);
                square.w = fmaf(v.w, v.w, square.w);
            }
        }
    }

    sum = add4(sum, shuffle_down4(sum, 16));
    square = add4(square, shuffle_down4(square, 16));
    sum = add4(sum, shuffle_down4(sum, 8));
    square = add4(square, shuffle_down4(square, 8));

    __shared__ float4 shared_sum[64];
    __shared__ float4 shared_square[64];

    if ((threadIdx.x & 31) < 8) {
        shared_sum[warp * 8 + group] = sum;
        shared_square[warp * 8 + group] = square;
    }
    __syncthreads();

    if (threadIdx.x < 8) {
        sum = shared_sum[group];
        square = shared_square[group];
        #pragma unroll
        for (int w = 1; w < 8; ++w) {
            sum = add4(sum, shared_sum[w * 8 + group]);
            square = add4(square, shared_square[w * 8 + group]);
        }
        const int64_t dst = (int64_t(n) * parts + p) * 8 + group;
        sums[dst] = sum;
        squares[dst] = square;
    }
}

__global__ void sample_stats32(
    const float* __restrict__ sums,
    const float* __restrict__ squares,
    float* __restrict__ means,
    float* __restrict__ second_moments,
    int spatial,
    int parts)
{
    const int c = threadIdx.x & 31;
    const int row = threadIdx.x >> 5;
    const int n = blockIdx.x;
    const int64_t base = int64_t(n) * parts * CHANNELS + c;

    float sum = 0.f;
    float square = 0.f;

    for (int p = row; p < parts; p += 8) {
        sum += sums[base + p * CHANNELS];
        square += squares[base + p * CHANNELS];
    }

    __shared__ float shared_sum[8][32];
    __shared__ float shared_square[8][32];
    shared_sum[row][c] = sum;
    shared_square[row][c] = square;
    __syncthreads();

    if (row == 0) {
        sum = 0.f;
        square = 0.f;
        #pragma unroll
        for (int r = 0; r < 8; ++r) {
            sum += shared_sum[r][c];
            square += shared_square[r][c];
        }
        const float inverse_spatial = 1.f / float(spatial);
        means[n * CHANNELS + c] = sum * inverse_spatial;
        second_moments[n * CHANNELS + c] = square * inverse_spatial;
    }
}

__global__ void channel_stats32(
    const float* __restrict__ means,
    const float* __restrict__ second_moments,
    const float* __restrict__ weight,
    const float* __restrict__ conv_bias,
    float* __restrict__ scales,
    float* running_mean,
    float* running_var,
    int batch,
    int spatial,
    float eps,
    float momentum,
    bool training)
{
    const int c = threadIdx.x & 31;
    const int row = threadIdx.x >> 5;

    float sum = 0.f;
    float square = 0.f;

    if (training) {
        for (int n = row; n < batch; n += 8) {
            sum += means[n * CHANNELS + c];
            square += second_moments[n * CHANNELS + c];
        }
    }

    __shared__ float shared_sum[8][32];
    __shared__ float shared_square[8][32];
    shared_sum[row][c] = sum;
    shared_square[row][c] = square;
    __syncthreads();

    if (row == 0) {
        float variance;
        if (training) {
            sum = 0.f;
            square = 0.f;
            #pragma unroll
            for (int r = 0; r < 8; ++r) {
                sum += shared_sum[r][c];
                square += shared_square[r][c];
            }
            const float inverse_batch = 1.f / float(batch);
            const float mean = sum * inverse_batch;
            variance = fmaxf(
                fmaf(-mean, mean, square * inverse_batch), 0.f);

            const float count = float(batch) * float(spatial);
            const float unbiased = variance * (count / (count - 1.f));
            const float biased_mean =
                mean + (conv_bias != nullptr ? conv_bias[c] : 0.f);

            running_mean[c] =
                (1.f - momentum) * running_mean[c] + momentum * biased_mean;
            running_var[c] =
                (1.f - momentum) * running_var[c] + momentum * unbiased;
        } else {
            variance = running_var[c];
        }
        scales[c] = weight[c] * rsqrtf(variance + eps);
    }
}

__global__ void center32(
    const float4* __restrict__ x,
    const float4* __restrict__ means,
    const float4* __restrict__ scales,
    float4* __restrict__ output,
    int64_t sample_vectors)
{
    const int n = blockIdx.y;
    const int group = threadIdx.x & 7;
    const float4 mean = means[n * 8 + group];
    const float4 scale = scales[group];

    const int64_t sample_base = int64_t(n) * sample_vectors;
    const int64_t begin =
        int64_t(blockIdx.x) * OUTPUT_VECTORS + threadIdx.x;

    #pragma unroll 4
    for (int k = 0; k < OUTPUT_VECTORS / 256; ++k) {
        const int64_t i = begin + k * 256;
        if (i < sample_vectors) {
            const float4 v = x[sample_base + i];
            output[sample_base + i] = make_float4(
                (v.x - mean.x) * scale.x,
                (v.y - mean.y) * scale.y,
                (v.z - mean.z) * scale.z,
                (v.w - mean.w) * scale.w);
        }
    }
}

torch::Tensor centered_bn32_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor conv_bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    double eps,
    double momentum,
    bool training)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat,
                "Expected CUDA float32 input");
    TORCH_CHECK(x.dim() == 5 && x.size(1) == CHANNELS &&
                x.is_contiguous(at::MemoryFormat::ChannelsLast3d),
                "Expected channels-last 3D input with 32 channels");

    const c10::cuda::CUDAGuard guard(x.device());
    const int batch = x.size(0);
    const int64_t spatial64 = x.size(2) * x.size(3) * x.size(4);

    TORCH_CHECK(batch > 0 && batch <= 65535 &&
                spatial64 > 0 && spatial64 <= 2147481600LL,
                "Unsupported input dimensions");
    TORCH_CHECK(!training || int64_t(batch) * spatial64 > 1,
                "Expected more than one value per channel");

    const int spatial = static_cast<int>(spatial64);
    const int parts = (spatial + TILE - 1) / TILE;

    auto partial = torch::empty({2, batch, parts, CHANNELS}, x.options());
    auto sample_stats = torch::empty({2, batch, CHANNELS}, x.options());
    auto scales = torch::empty({CHANNELS}, x.options());
    auto output = torch::empty_like(x);

    float* sums = partial.data_ptr<float>();
    float* squares = sums + int64_t(batch) * parts * CHANNELS;
    float* means = sample_stats.data_ptr<float>();
    float* second_moments = means + int64_t(batch) * CHANNELS;
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    partial_stats32<<<dim3(parts, batch), 256, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        reinterpret_cast<float4*>(sums),
        reinterpret_cast<float4*>(squares),
        spatial, parts);

    sample_stats32<<<batch, 256, 0, stream>>>(
        sums, squares, means, second_moments, spatial, parts);

    channel_stats32<<<1, 256, 0, stream>>>(
        means, second_moments, weight.data_ptr<float>(),
        conv_bias.numel() ? conv_bias.data_ptr<float>() : nullptr,
        scales.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        batch, spatial, float(eps), float(momentum), training);

    const int64_t sample_vectors = spatial64 * 8;
    const dim3 grid(
        (sample_vectors + OUTPUT_VECTORS - 1) / OUTPUT_VECTORS, batch);

    center32<<<grid, 256, 0, stream>>>(
        reinterpret_cast<const float4*>(x.data_ptr<float>()),
        reinterpret_cast<const float4*>(means),
        reinterpret_cast<const float4*>(scales.data_ptr<float>()),
        reinterpret_cast<float4*>(output.data_ptr<float>()),
        sample_vectors);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_centered_bn = load_inline(
    name="conv_transpose3d_centered_bn_vector32_v3",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["centered_bn32_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class _CenteredBatchNorm(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, conv_bias,
                running_mean, running_var, eps, momentum, training):
        output = _centered_bn.centered_bn32_cuda(
            x, weight, conv_bias, running_mean, running_var,
            eps, momentum, training,
        )
        ctx.save_for_backward(
            x, weight,
            x.new_empty(0) if training else running_var.clone(),
        )
        ctx.eps = eps
        ctx.training = training
        ctx.has_conv_bias = conv_bias.numel() != 0
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, running_var = ctx.saved_tensors
        spatial_dims = (2, 3, 4)
        channel_dims = (0, 2, 3, 4)
        shape = (1, -1, 1, 1, 1)

        centered = x - x.mean(dim=spatial_dims, keepdim=True)
        grad_centered = (
            grad_output - grad_output.mean(dim=spatial_dims, keepdim=True)
        )

        if ctx.training:
            variance, mean = torch.var_mean(
                x, dim=channel_dims, unbiased=False, keepdim=True,
            )
            inv_std = torch.rsqrt(variance + ctx.eps)
        else:
            inv_std = torch.rsqrt(running_var.view(shape) + ctx.eps)

        correlation = (grad_output * centered).sum(
            dim=channel_dims, keepdim=True,
        )
        grad_x = grad_centered * weight.view(shape) * inv_std

        if ctx.training:
            count = x.numel() // x.size(1)
            grad_x = grad_x - (
                (x - mean) * weight.view(shape)
                * inv_std.pow(3) * correlation / count
            )

        grad_weight = (correlation * inv_std).reshape_as(weight)
        grad_bias = torch.zeros_like(weight)
        grad_conv_bias = (
            torch.zeros_like(weight) if ctx.has_conv_bias else None
        )

        return (
            grad_x, grad_weight, grad_bias, grad_conv_bias,
            None, None, None, None, None,
        )


class ModelNew(nn.Module):
    def __init__(self, in_channels, out_channels, kernel_size,
                 stride, padding, bias=True):
        super().__init__()
        self.conv_transpose = nn.ConvTranspose3d(
            in_channels, out_channels, kernel_size,
            stride=stride, padding=padding, bias=bias,
        )
        self.batch_norm = nn.BatchNorm3d(out_channels)
        self.conv_transpose.to(memory_format=torch.channels_last_3d)

    def forward(self, x):
        conv = self.conv_transpose
        bn = self.batch_norm

        if (
            not x.is_cuda
            or x.dim() != 5
            or x.dtype != torch.float32
            or conv.weight.dtype != torch.float32
            or conv.out_channels != 32
            or not bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
            or bn.running_mean.dtype != torch.float32
            or x.numel() == 0
            or x.size(0) > 65535
        ):
            y = bn(conv(x))
            return y - y.mean(dim=(2, 3, 4), keepdim=True)

        # Both constant channel shifts cancel under spatial centering.
        y = F.conv_transpose3d(
            x.contiguous(memory_format=torch.channels_last_3d),
            conv.weight,
            None,
            conv.stride,
            conv.padding,
            conv.output_padding,
            conv.groups,
            conv.dilation,
        )

        spatial = y.size(2) * y.size(3) * y.size(4)
        if y.dtype != torch.float32 or spatial > 2147481600:
            if conv.bias is not None:
                y = y + conv.bias.to(y.dtype).view(1, -1, 1, 1, 1)
            y = bn(y)
            return y - y.mean(dim=(2, 3, 4), keepdim=True)

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

        momentum = bn.momentum if bn.momentum is not None else 0.0
        if bn.training:
            bn.num_batches_tracked.add_(1)
            if bn.momentum is None:
                momentum = 1.0 / bn.num_batches_tracked.item()

        conv_bias = conv.bias if conv.bias is not None else y.new_empty(0)
        return _CenteredBatchNorm.apply(
            y, bn.weight, bn.bias, conv_bias,
            bn.running_mean, bn.running_var,
            bn.eps, momentum, bn.training,
        )
```