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

_cpp_source = r"""
torch::Tensor centered_bn_cuda(
    torch::Tensor x, torch::Tensor weight,
    torch::Tensor running_mean, torch::Tensor running_var,
    double eps, double momentum, bool batch_stats, bool update_running);
"""

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

constexpr int THREADS = 256;
constexpr int TILE = 4096;

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

__device__ __forceinline__ Moments combine(Moments a, Moments b) {
    if (!a.count) return b;
    if (!b.count) return a;
    int count = a.count + b.count;
    float delta = b.mean - a.mean;
    float fraction = float(b.count) / float(count);
    a.mean += delta * fraction;
    a.m2 += b.m2 + delta * delta * float(a.count) * fraction;
    a.count = count;
    return a;
}

__device__ __forceinline__ Moments warp_reduce(Moments a) {
    for (int offset = 16; offset > 0; offset >>= 1) {
        Moments b;
        b.mean = __shfl_down_sync(0xffffffff, a.mean, offset);
        b.m2 = __shfl_down_sync(0xffffffff, a.m2, offset);
        b.count = __shfl_down_sync(0xffffffff, a.count, offset);
        if ((threadIdx.x & 31) + offset < 32)
            a = combine(a, b);
    }
    return a;
}

__device__ __forceinline__ Moments block_reduce(Moments a) {
    __shared__ Moments shared[8];
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;
    a = warp_reduce(a);
    if (lane == 0) shared[warp] = a;
    __syncthreads();
    a = {0.0f, 0.0f, 0};
    if (warp == 0) {
        if (lane < 8) a = shared[lane];
        a = warp_reduce(a);
    }
    return a;
}

__global__ void partial_moments_kernel(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int spatial, int parts)
{
    int nc = blockIdx.y;
    int part = blockIdx.x;
    int begin = part * TILE;
    int end = min(begin + TILE, spatial);
    const float* row = x + int64_t(nc) * spatial;

    Moments a = {0.0f, 0.0f, 0};
    for (int i = begin + threadIdx.x; i < end; i += THREADS) {
        float value = row[i];
        ++a.count;
        float delta = value - a.mean;
        a.mean += delta / float(a.count);
        a.m2 += delta * (value - a.mean);
    }
    a = block_reduce(a);
    if (threadIdx.x == 0)
        partial[int64_t(nc) * parts + part] = make_float2(a.mean, a.m2);
}

__global__ void finish_moments_kernel(
    const float2* __restrict__ partial,
    const float* __restrict__ weight,
    float* __restrict__ sample_mean,
    float* __restrict__ scale,
    float* running_mean,
    float* running_var,
    int batch, int channels, int spatial, int parts,
    float eps, float momentum, bool batch_stats, bool update_running)
{
    int c = blockIdx.x;
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;
    Moments channel = {0.0f, 0.0f, 0};

    for (int n = warp; n < batch; n += 8) {
        Moments sample = {0.0f, 0.0f, 0};
        int nc = n * channels + c;
        for (int p = lane; p < parts; p += 32) {
            float2 value = partial[int64_t(nc) * parts + p];
            Moments other = {
                value.x, value.y, min(TILE, spatial - p * TILE)
            };
            sample = combine(sample, other);
        }
        sample = warp_reduce(sample);
        if (lane == 0) {
            sample_mean[nc] = sample.mean;
            channel = combine(channel, sample);
        }
    }

    channel = block_reduce(channel);
    if (threadIdx.x == 0) {
        float variance;
        if (batch_stats) {
            variance = fmaxf(channel.m2 / float(channel.count), 0.0f);
            if (update_running) {
                float unbiased = channel.m2 / float(channel.count - 1);
                running_mean[c] =
                    (1.0f - momentum) * running_mean[c]
                    + momentum * channel.mean;
                running_var[c] =
                    (1.0f - momentum) * running_var[c]
                    + momentum * unbiased;
            }
        } else {
            variance = running_var[c];
        }
        scale[c] = weight[c] * rsqrtf(variance + eps);
    }
}

__global__ void normalize_center_kernel(
    const float* __restrict__ x,
    const float* __restrict__ sample_mean,
    const float* __restrict__ scale,
    float* __restrict__ out,
    int channels, int spatial)
{
    int nc = blockIdx.y;
    int begin = blockIdx.x * TILE;
    int end = min(begin + TILE, spatial);
    int64_t base = int64_t(nc) * spatial;
    float mean = sample_mean[nc];
    float multiplier = scale[nc % channels];

    for (int i = begin + threadIdx.x; i < end; i += THREADS)
        out[base + i] = (x[base + i] - mean) * multiplier;
}

torch::Tensor centered_bn_cuda(
    torch::Tensor x, torch::Tensor weight,
    torch::Tensor running_mean, torch::Tensor running_var,
    double eps, double momentum, bool batch_stats, bool update_running)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat,
                "Expected a CUDA float32 tensor");
    TORCH_CHECK(x.dim() == 5 && x.is_contiguous(),
                "Expected contiguous NCDHW input");

    c10::cuda::CUDAGuard guard(x.device());
    int batch = x.size(0);
    int channels = x.size(1);
    int spatial = x.size(2) * x.size(3) * x.size(4);
    TORCH_CHECK(batch > 0 && spatial > 0, "Input must be nonempty");
    TORCH_CHECK(!batch_stats || int64_t(batch) * spatial > 1,
                "Batch normalization requires more than one value per channel");

    int parts = (spatial + TILE - 1) / TILE;
    auto partial = torch::empty({batch, channels, parts, 2}, x.options());
    auto means = torch::empty({batch, channels}, x.options());
    auto scale = torch::empty({channels}, x.options());
    auto out = torch::empty_like(x);
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    dim3 grid(parts, batch * channels);
    partial_moments_kernel<<<grid, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        spatial, parts);

    finish_moments_kernel<<<channels, THREADS, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        weight.data_ptr<float>(), means.data_ptr<float>(),
        scale.data_ptr<float>(), running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        batch, channels, spatial, parts,
        float(eps), float(momentum), batch_stats, update_running);

    normalize_center_kernel<<<grid, THREADS, 0, stream>>>(
        x.data_ptr<float>(), means.data_ptr<float>(),
        scale.data_ptr<float>(), out.data_ptr<float>(),
        channels, spatial);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return out;
}
"""

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


class _CenteredBatchNorm(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, running_mean, running_var,
                eps, momentum, batch_stats, update_running):
        out = _centered_bn.centered_bn_cuda(
            x, weight, running_mean, running_var,
            eps, momentum, batch_stats, update_running,
        )
        ctx.save_for_backward(x, weight, running_var)
        ctx.eps = eps
        ctx.batch_stats = batch_stats
        return out

    @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.batch_stats:
            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.batch_stats:
            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)
        return (
            grad_x, grad_weight, grad_bias,
            None, 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)

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

        if (
            not x.is_cuda
            or x.dtype != torch.float32
            or not x.is_contiguous()
            or x.numel() == 0
            or x.size(0) * x.size(1) > 65535
            or not bn.affine
            or not bn.track_running_stats
        ):
            x = bn(x)
            return x - x.mean(dim=(2, 3, 4), keepdim=True)

        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()

        return _CenteredBatchNorm.apply(
            x, bn.weight, bn.bias, bn.running_mean, bn.running_var,
            bn.eps, momentum, bn.training, bn.training,
        )
```