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


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

torch::Tensor centered_batch_norm_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>

constexpr int THREADS = 256;
constexpr int CHUNK = 2048;

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

__global__ void partial_stats_kernel(
    const float* __restrict__ x,
    float* __restrict__ partial_sum,
    float* __restrict__ partial_square,
    int spatial,
    int parts)
{
    int row = blockIdx.x / parts;
    int part = blockIdx.x % parts;
    int start = part * CHUNK;
    int end = min(start + CHUNK, spatial);

    float sum = 0.0f;
    float square = 0.0f;
    for (int i = start + threadIdx.x; i < end; i += THREADS) {
        float value = x[(int64_t)row * spatial + i];
        sum += value;
        square = fmaf(value, value, square);
    }

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

    __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;
    }
    __syncthreads();

    if (warp == 0) {
        sum = warp_sum(lane < 8 ? sums[lane] : 0.0f);
        square = warp_sum(lane < 8 ? squares[lane] : 0.0f);
        if (lane == 0) {
            partial_sum[blockIdx.x] = sum;
            partial_square[blockIdx.x] = square;
        }
    }
}

__global__ void finalize_stats_kernel(
    const float* __restrict__ partial_sum,
    const float* __restrict__ partial_square,
    const float* __restrict__ weight,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    float* __restrict__ instance_mean,
    float* __restrict__ scale,
    int batch,
    int channels,
    int spatial,
    int parts,
    bool training,
    float momentum,
    float eps)
{
    int channel = blockIdx.x;
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;

    float batch_sum = 0.0f;
    float batch_square = 0.0f;

    for (int n = warp; n < batch; n += 8) {
        int64_t offset = ((int64_t)n * channels + channel) * parts;
        float sum = 0.0f;
        float square = 0.0f;

        for (int p = lane; p < parts; p += 32) {
            sum += partial_sum[offset + p];
            if (training)
                square += partial_square[offset + p];
        }

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

        if (lane == 0) {
            instance_mean[n * channels + channel] = sum / spatial;
            batch_sum += sum;
            batch_square += square;
        }
    }

    __shared__ float sums[8];
    __shared__ float squares[8];
    if (lane == 0) {
        sums[warp] = batch_sum;
        squares[warp] = batch_square;
    }
    __syncthreads();

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

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

                running_mean[channel] =
                    (1.0f - momentum) * running_mean[channel]
                    + momentum * mean;
                running_var[channel] =
                    (1.0f - momentum) * running_var[channel]
                    + momentum * variance * (count / (count - 1.0f));
            } else {
                variance = running_var[channel];
            }

            scale[channel] = weight[channel] * rsqrtf(variance + eps);
        }
    }
}

__global__ void center_scale_kernel(
    const float* __restrict__ x,
    const float* __restrict__ instance_mean,
    const float* __restrict__ scale,
    float* __restrict__ output,
    int channels,
    int spatial,
    int parts)
{
    int row = blockIdx.x / parts;
    int part = blockIdx.x % parts;
    int channel = row % channels;
    int start = part * CHUNK;
    int end = min(start + CHUNK, spatial);

    float mean = instance_mean[row];
    float multiplier = scale[channel];

    for (int i = start + threadIdx.x; i < end; i += THREADS) {
        int64_t index = (int64_t)row * spatial + i;
        output[index] = (x[index] - mean) * multiplier;
    }
}

torch::Tensor centered_batch_norm_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    bool training,
    double momentum,
    double eps)
{
    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 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 parts = (spatial + CHUNK - 1) / CHUNK;

    auto output = torch::empty_like(x);
    if (x.numel() == 0)
        return output;

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

    auto partial_sum = torch::empty(
        {(int64_t)batch * channels * parts}, x.options());
    auto partial_square = torch::empty_like(partial_sum);
    auto instance_mean = torch::empty(
        {(int64_t)batch * channels}, x.options());
    auto scale = torch::empty({channels}, x.options());

    auto stream = at::cuda::getCurrentCUDAStream();
    int blocks = batch * channels * parts;

    partial_stats_kernel<<<blocks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        partial_sum.data_ptr<float>(),
        partial_square.data_ptr<float>(),
        spatial, parts);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    finalize_stats_kernel<<<channels, THREADS, 0, stream>>>(
        partial_sum.data_ptr<float>(),
        partial_square.data_ptr<float>(),
        weight.data_ptr<float>(),
        running_mean.data_ptr<float>(),
        running_var.data_ptr<float>(),
        instance_mean.data_ptr<float>(),
        scale.data_ptr<float>(),
        batch, channels, spatial, parts,
        training, (float)momentum, (float)eps);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    center_scale_kernel<<<blocks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        instance_mean.data_ptr<float>(),
        scale.data_ptr<float>(),
        output.data_ptr<float>(),
        channels, spatial, parts);
    C10_CUDA_KERNEL_LAUNCH_CHECK();

    return output;
}
"""

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


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 bn.weight.dtype != torch.float32
            or x.numel() == 0
        ):
            x = bn(x)
            return x - x.mean(dim=(2, 3, 4), keepdim=True)

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

        return _centered_bn.centered_batch_norm_cuda(
            x.contiguous(),
            bn.weight,
            bn.running_mean,
            bn.running_var,
            bn.training,
            0.0 if momentum is None else momentum,
            bn.eps,
        )
```