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


_cpp_source = """
#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,
    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>
#include <vector>

constexpr int THREADS = 256;
constexpr int CHUNK = 4096;

__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;
}

__global__ void partial_moments(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int spatial,
    int chunks)
{
    int nc = blockIdx.x / chunks;
    int chunk = blockIdx.x % chunks;
    int start = chunk * CHUNK;

    float sum = 0.0f;
    float square_sum = 0.0f;

    for (int j = threadIdx.x; j < CHUNK; j += THREADS) {
        int s = start + j;
        if (s < spatial) {
            float v = x[(int64_t)nc * spatial + s];
            sum += v;
            square_sum += v * v;
        }
    }

    sum = warp_sum(sum);
    square_sum = warp_sum(square_sum);

    __shared__ float sums[8];
    __shared__ float squares[8];

    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;

    if (lane == 0) {
        sums[warp] = sum;
        squares[warp] = square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        sum = lane < 8 ? sums[lane] : 0.0f;
        square_sum = lane < 8 ? squares[lane] : 0.0f;
        sum = warp_sum(sum);
        square_sum = warp_sum(square_sum);
        if (lane == 0)
            partial[blockIdx.x] = make_float2(sum, square_sum);
    }
}

__global__ void finish_statistics(
    const float2* __restrict__ partial,
    float* __restrict__ instance_mean,
    float* __restrict__ inverse_std,
    float* __restrict__ global_mean,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    int batch,
    int channels,
    int spatial,
    int chunks,
    float eps,
    float momentum,
    bool training)
{
    int c = blockIdx.x;
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;

    float channel_sum = 0.0f;
    float channel_square_sum = 0.0f;

    for (int n = warp; n < batch; n += 8) {
        float sum = 0.0f;
        float square_sum = 0.0f;

        int64_t base = ((int64_t)n * channels + c) * chunks;
        for (int k = lane; k < chunks; k += 32) {
            float2 v = partial[base + k];
            sum += v.x;
            square_sum += v.y;
        }

        sum = warp_sum(sum);
        square_sum = warp_sum(square_sum);

        if (lane == 0) {
            instance_mean[n * channels + c] = sum / spatial;
            channel_sum += sum;
            channel_square_sum += square_sum;
        }
    }

    __shared__ float sums[8];
    __shared__ float squares[8];

    if (lane == 0) {
        sums[warp] = channel_sum;
        squares[warp] = channel_square_sum;
    }
    __syncthreads();

    if (warp == 0) {
        float sum = lane < 8 ? sums[lane] : 0.0f;
        float square_sum = lane < 8 ? squares[lane] : 0.0f;
        sum = warp_sum(sum);
        square_sum = warp_sum(square_sum);

        if (lane == 0) {
            if (training) {
                float count = (float)((int64_t)batch * spatial);
                float mean = sum / count;
                float variance =
                    fmaxf(square_sum / count - mean * mean, 0.0f);

                global_mean[c] = mean;
                inverse_std[c] = rsqrtf(variance + eps);

                running_mean[c] =
                    (1.0f - momentum) * running_mean[c] + momentum * mean;

                float unbiased = variance * count / (count - 1.0f);
                running_var[c] =
                    (1.0f - momentum) * running_var[c] + momentum * unbiased;
            } else {
                global_mean[c] = running_mean[c];
                inverse_std[c] = rsqrtf(running_var[c] + eps);
            }
        }
    }
}

__global__ void normalize_and_center(
    const float* __restrict__ x,
    const float* __restrict__ instance_mean,
    const float* __restrict__ inverse_std,
    const float* __restrict__ weight,
    float* __restrict__ output,
    int64_t total,
    int channels,
    int spatial)
{
    int64_t base = (int64_t)blockIdx.x * THREADS * 4 + threadIdx.x;

    #pragma unroll
    for (int j = 0; j < 4; ++j) {
        int64_t i = base + j * THREADS;
        if (i < total) {
            int64_t nc = i / spatial;
            int c = nc % channels;
            float scale = inverse_std[c] * weight[c];
            output[i] = (x[i] - instance_mean[nc]) * scale;
        }
    }
}

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

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

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

    auto partial = torch::empty({batch, channels, chunks, 2}, x.options());
    auto means = torch::empty({batch, channels}, x.options());
    auto inverse_std = torch::empty({channels}, x.options());
    auto global_mean = torch::empty({channels}, x.options());
    auto output = torch::empty_like(x);

    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    partial_moments<<<batch * channels * chunks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        reinterpret_cast<float2*>(partial.data_ptr<float>()),
        spatial,
        chunks);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finish_statistics<<<channels, THREADS, 0, stream>>>(
        reinterpret_cast<const float2*>(partial.data_ptr<float>()),
        means.data_ptr<float>(),
        inverse_std.data_ptr<float>(),
        global_mean.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        batch,
        channels,
        spatial,
        chunks,
        static_cast<float>(eps),
        static_cast<float>(momentum),
        training);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    int64_t total = x.numel();
    int blocks = (total + THREADS * 4 - 1) / (THREADS * 4);

    normalize_and_center<<<blocks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        means.data_ptr<float>(),
        inverse_std.data_ptr<float>(),
        weight.data_ptr<float>(),
        output.data_ptr<float>(),
        total,
        channels,
        spatial);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return {output, means, inverse_std, global_mean};
}
"""

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


class _CenteredBatchNorm(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, batch_norm, momentum):
        output, means, inverse_std, global_mean = (
            _centered_bn.centered_bn_cuda(
                x,
                weight,
                batch_norm.running_mean,
                batch_norm.running_var,
                batch_norm.eps,
                momentum,
                batch_norm.training,
            )
        )
        ctx.training = batch_norm.training
        ctx.save_for_backward(x, weight, means, inverse_std, global_mean)
        return output

    @staticmethod
    def backward(ctx, grad_output):
        x, weight, means, inverse_std, global_mean = ctx.saved_tensors
        spatial_dims = (2, 3, 4)
        channel_dims = (0, 2, 3, 4)
        channel_shape = (1, -1, 1, 1, 1)

        centered = x - means[:, :, None, None, None]
        correlation = (grad_output * centered).sum(channel_dims)
        inverse = inverse_std.view(channel_shape)

        grad_x = None
        if ctx.needs_input_grad[0]:
            grad_centered = grad_output - grad_output.mean(
                spatial_dims, keepdim=True
            )
            if ctx.training:
                count = x.numel() // x.shape[1]
                correction = (
                    correlation * inverse_std.square() / count
                ).view(channel_shape)
                grad_centered = grad_centered - (
                    x - global_mean.view(channel_shape)
                ) * correction
            grad_x = grad_centered * (
                weight * inverse_std
            ).view(channel_shape)

        grad_weight = (
            correlation * inverse_std if ctx.needs_input_grad[1] else None
        )
        grad_bias = (
            torch.zeros_like(weight) if ctx.needs_input_grad[2] else None
        )
        return grad_x, grad_weight, grad_bias, 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.dtype != torch.float32
        ):
            x = bn(x)
            return x - x.mean(dim=(2, 3, 4), keepdim=True)

        momentum = bn.momentum
        if bn.training:
            with torch.no_grad():
                bn.num_batches_tracked.add_(1)
            if momentum is None:
                momentum = 1.0 / bn.num_batches_tracked.item()
        elif momentum is None:
            momentum = 0.0

        return _CenteredBatchNorm.apply(
            x, bn.weight, bn.bias, bn, momentum
        )
```