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


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

std::vector<torch::Tensor> centered_bn_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    bool training,
    double momentum,
    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>
#include <vector>

constexpr int THREADS = 256;
constexpr int ITEMS = 16;
constexpr int CHUNK = THREADS * ITEMS;

__device__ __forceinline__ float block_sum(float value) {
    __shared__ float scratch[THREADS / 32];

    __syncthreads();
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1)
        value += __shfl_down_sync(0xffffffff, value, offset);

    if ((threadIdx.x & 31) == 0)
        scratch[threadIdx.x >> 5] = value;
    __syncthreads();

    float total = threadIdx.x < THREADS / 32
        ? scratch[threadIdx.x] : 0.0f;

    if (threadIdx.x < 32) {
        #pragma unroll
        for (int offset = 16; offset > 0; offset >>= 1)
            total += __shfl_down_sync(0xffffffff, total, offset);
    }

    if (threadIdx.x == 0)
        scratch[0] = total;
    __syncthreads();
    return scratch[0];
}

__global__ void partial_statistics(
    const float* __restrict__ x,
    float* __restrict__ partial_mean,
    float* __restrict__ partial_m2,
    int spatial,
    int parts)
{
    const int nc = blockIdx.x / parts;
    const int part = blockIdx.x % parts;
    const int start = part * CHUNK;
    const int count = min(CHUNK, spatial - start);
    const int64_t base = static_cast<int64_t>(nc) * spatial + start;

    float values[ITEMS];
    float sum = 0.0f;

    #pragma unroll
    for (int j = 0; j < ITEMS; ++j) {
        const int offset = threadIdx.x + j * THREADS;
        values[j] = offset < count ? x[base + offset] : 0.0f;
        sum += values[j];
    }

    const float mean = block_sum(sum) / static_cast<float>(count);
    float m2 = 0.0f;

    #pragma unroll
    for (int j = 0; j < ITEMS; ++j) {
        const int offset = threadIdx.x + j * THREADS;
        if (offset < count) {
            const float difference = values[j] - mean;
            m2 += difference * difference;
        }
    }

    m2 = block_sum(m2);
    if (threadIdx.x == 0) {
        partial_mean[blockIdx.x] = mean;
        partial_m2[blockIdx.x] = m2;
    }
}

__global__ void finish_statistics(
    const float* __restrict__ partial_mean,
    const float* __restrict__ partial_m2,
    float* __restrict__ instance_mean,
    float* __restrict__ channel_mean,
    float* __restrict__ invstd,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    int batch,
    int channels,
    int spatial,
    int parts,
    bool training,
    float momentum,
    float eps)
{
    const int channel = blockIdx.x;
    const int entries = batch * parts;
    const float sample_count = static_cast<float>(batch) * spatial;

    if (training) {
        float sum = 0.0f;
        for (int i = threadIdx.x; i < entries; i += THREADS) {
            const int n = i / parts;
            const int part = i % parts;
            const int count = min(CHUNK, spatial - part * CHUNK);
            const int index = (n * channels + channel) * parts + part;
            sum += partial_mean[index] * static_cast<float>(count);
        }

        const float mean = block_sum(sum) / sample_count;
        float m2 = 0.0f;

        for (int i = threadIdx.x; i < entries; i += THREADS) {
            const int n = i / parts;
            const int part = i % parts;
            const int count = min(CHUNK, spatial - part * CHUNK);
            const int index = (n * channels + channel) * parts + part;
            const float difference = partial_mean[index] - mean;
            m2 += partial_m2[index]
                + static_cast<float>(count) * difference * difference;
        }

        m2 = block_sum(m2);
        if (threadIdx.x == 0) {
            channel_mean[channel] = mean;
            invstd[channel] = rsqrtf(m2 / sample_count + eps);

            running_mean[channel] =
                (1.0f - momentum) * running_mean[channel] + momentum * mean;
            running_var[channel] =
                (1.0f - momentum) * running_var[channel]
                + momentum * (m2 / (sample_count - 1.0f));
        }
    } else if (threadIdx.x == 0) {
        channel_mean[channel] = running_mean[channel];
        invstd[channel] = rsqrtf(running_var[channel] + eps);
    }

    for (int n = threadIdx.x; n < batch; n += THREADS) {
        float sum = 0.0f;
        const int base = (n * channels + channel) * parts;

        for (int part = 0; part < parts; ++part) {
            const int count = min(CHUNK, spatial - part * CHUNK);
            sum += partial_mean[base + part] * static_cast<float>(count);
        }

        instance_mean[n * channels + channel] =
            sum / static_cast<float>(spatial);
    }
}

__global__ void normalize_and_center(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ instance_mean,
    const float* __restrict__ invstd,
    float* __restrict__ output,
    int64_t total,
    int spatial,
    int channels)
{
    const int64_t base =
        (static_cast<int64_t>(blockIdx.x) * blockDim.x + threadIdx.x) * 4;
    if (base >= total)
        return;

    const int64_t nc = base / spatial;
    const int channel = nc % channels;

    if (base + 3 < total && base % spatial + 3 < spatial) {
        const float mean = instance_mean[nc];
        const float scale = weight[channel] * invstd[channel];
        const float4 values = reinterpret_cast<const float4*>(x)[base / 4];

        float4 result;
        result.x = (values.x - mean) * scale;
        result.y = (values.y - mean) * scale;
        result.z = (values.z - mean) * scale;
        result.w = (values.w - mean) * scale;

        reinterpret_cast<float4*>(output)[base / 4] = result;
    } else {
        #pragma unroll
        for (int j = 0; j < 4; ++j) {
            const int64_t index = base + j;
            if (index < total) {
                const int64_t current_nc = index / spatial;
                const int current_channel = current_nc % channels;
                output[index] =
                    (x[index] - instance_mean[current_nc])
                    * (weight[current_channel] * invstd[current_channel]);
            }
        }
    }
}

std::vector<torch::Tensor> centered_bn_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    bool training,
    double momentum,
    double eps)
{
    const c10::cuda::CUDAGuard device_guard(x.device());
    const auto stream = at::cuda::getCurrentCUDAStream();

    const int batch = x.size(0);
    const int channels = x.size(1);
    const int spatial = x.size(2) * x.size(3) * x.size(4);
    const int parts = (spatial + CHUNK - 1) / CHUNK;

    TORCH_CHECK(!training || static_cast<int64_t>(batch) * spatial > 1,
                "Expected more than one value per channel when training");

    auto partial_mean = torch::empty({batch, channels, parts}, x.options());
    auto partial_m2 = torch::empty_like(partial_mean);
    auto instance_mean = torch::empty({batch, channels}, x.options());
    auto channel_mean = torch::empty({channels}, x.options());
    auto invstd = torch::empty({channels}, x.options());
    auto output = torch::empty_like(x);

    partial_statistics<<<batch * channels * parts, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        partial_mean.data_ptr<float>(),
        partial_m2.data_ptr<float>(),
        spatial,
        parts);

    finish_statistics<<<channels, THREADS, 0, stream>>>(
        partial_mean.data_ptr<float>(),
        partial_m2.data_ptr<float>(),
        instance_mean.data_ptr<float>(),
        channel_mean.data_ptr<float>(),
        invstd.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        batch,
        channels,
        spatial,
        parts,
        training,
        static_cast<float>(momentum),
        static_cast<float>(eps));

    const int64_t blocks = (x.numel() + THREADS * 4 - 1) / (THREADS * 4);
    normalize_and_center<<<blocks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        instance_mean.data_ptr<float>(),
        invstd.data_ptr<float>(),
        output.data_ptr<float>(),
        x.numel(),
        spatial,
        channels);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return {output, instance_mean, channel_mean, invstd};
}
"""

_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,
        training, momentum, eps
    ):
        output, instance_mean, channel_mean, invstd = (
            _centered_bn.centered_bn_cuda(
                x, weight, running_mean, running_var,
                training, momentum, eps
            )
        )
        ctx.training = training
        ctx.save_for_backward(
            x, weight, instance_mean, channel_mean, invstd
        )
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, instance_mean, channel_mean, invstd = ctx.saved_tensors
        channel_shape = (1, -1, 1, 1, 1)
        inverse = invstd.view(channel_shape)
        spatial_dims = (2, 3, 4)
        reduction_dims = (0, 2, 3, 4)

        grad_x = grad_weight = grad_bias = None

        if ctx.needs_input_grad[0]:
            centered_grad = grad_output - grad_output.mean(
                dim=spatial_dims, keepdim=True
            )
            if ctx.training:
                centered_x = x - channel_mean.view(channel_shape)
                correction = (centered_grad * centered_x).mean(
                    dim=reduction_dims, keepdim=True
                )
                centered_grad = (
                    centered_grad
                    - centered_grad.mean(dim=reduction_dims, keepdim=True)
                    - centered_x * inverse.square() * correction
                )
            grad_x = centered_grad * (weight.view(channel_shape) * inverse)

        if ctx.needs_input_grad[1]:
            centered_x = x - instance_mean[:, :, None, None, None]
            grad_weight = (grad_output * centered_x).sum(
                dim=reduction_dims
            ) * invstd

        if ctx.needs_input_grad[2]:
            grad_bias = torch.zeros_like(weight)

        return (
            grad_x, grad_weight, grad_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)

    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 not bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
            or bn.running_mean is None
            or bn.running_var is None
        ):
            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 / float(bn.num_batches_tracked.item())

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