```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,
    torch::Tensor counter,
    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 TILE = 2048;

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

template<bool Training>
__device__ __forceinline__ Moment merge(Moment a, Moment b) {
    if (a.count == 0) return b;
    if (b.count == 0) return a;
    float delta = b.mean - a.mean;
    int count = a.count + b.count;
    float ratio = float(b.count) / float(count);
    a.mean += delta * ratio;
    if (Training) {
        a.m2 += b.m2 + delta * delta * float(a.count) * ratio;
    }
    a.count = count;
    return a;
}

template<bool Training>
__device__ __forceinline__ Moment warp_reduce(Moment a) {
    #pragma unroll
    for (int offset = 16; offset > 0; offset >>= 1) {
        Moment 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);
        a = merge<Training>(a, b);
    }
    return a;
}

template<bool Training>
__global__ void partial_moments(
    const float* __restrict__ x,
    float2* __restrict__ partial,
    int spatial,
    int parts)
{
    int group = blockIdx.x / parts;
    int part = blockIdx.x % parts;
    int start = part * TILE;
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;

    Moment a{0.0f, 0.0f, 0};

    #pragma unroll
    for (int j = 0; j < TILE / THREADS; ++j) {
        int position = start + threadIdx.x + j * THREADS;
        if (position < spatial) {
            float value = x[(int64_t)group * spatial + position];
            float delta = value - a.mean;
            ++a.count;
            a.mean += delta / float(a.count);
            if (Training) {
                a.m2 = fmaf(delta, value - a.mean, a.m2);
            }
        }
    }

    a = warp_reduce<Training>(a);
    __shared__ Moment shared[THREADS / 32];
    if (lane == 0) shared[warp] = a;
    __syncthreads();

    if (warp == 0) {
        a = lane < THREADS / 32
            ? shared[lane] : Moment{0.0f, 0.0f, 0};
        a = warp_reduce<Training>(a);
        if (lane == 0) {
            partial[blockIdx.x] = make_float2(a.mean, a.m2);
        }
    }
}

template<bool Training>
__global__ void finish_moments(
    const float2* __restrict__ partial,
    float* __restrict__ instance_mean,
    float2* __restrict__ stats,
    float* __restrict__ running_mean,
    float* __restrict__ running_var,
    int64_t* __restrict__ counter,
    int batch,
    int channels,
    int spatial,
    int parts,
    float momentum,
    float eps)
{
    int channel = blockIdx.x;
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;
    Moment accumulated{0.0f, 0.0f, 0};

    for (int n = warp; n < batch; n += THREADS / 32) {
        Moment a{0.0f, 0.0f, 0};
        int group = n * channels + channel;

        for (int p = lane; p < parts; p += 32) {
            float2 value = partial[(int64_t)group * parts + p];
            int count = min(TILE, spatial - p * TILE);
            Moment b{value.x, value.y, count};
            a = merge<Training>(a, b);
        }

        a = warp_reduce<Training>(a);
        if (lane == 0) {
            instance_mean[group] = a.mean;
            accumulated = merge<Training>(accumulated, a);
        }
    }

    __shared__ Moment shared[THREADS / 32];
    if (lane == 0) shared[warp] = accumulated;
    __syncthreads();

    if (warp == 0) {
        Moment a = lane < THREADS / 32
            ? shared[lane] : Moment{0.0f, 0.0f, 0};
        a = warp_reduce<Training>(a);

        if (lane == 0) {
            if (Training) {
                float variance = fmaxf(a.m2 / float(a.count), 0.0f);
                stats[channel] = make_float2(a.mean, rsqrtf(variance + eps));

                running_mean[channel] =
                    (1.0f - momentum) * running_mean[channel]
                    + momentum * a.mean;
                running_var[channel] =
                    (1.0f - momentum) * running_var[channel]
                    + momentum * (a.m2 / float(a.count - 1));

                if (channel == 0) {
                    *counter += 1;
                }
            } else {
                stats[channel] = make_float2(
                    running_mean[channel],
                    rsqrtf(running_var[channel] + eps));
            }
        }
    }
}

__global__ void apply_centered_bn(
    const float* __restrict__ x,
    const float* __restrict__ weight,
    const float* __restrict__ instance_mean,
    const float2* __restrict__ stats,
    float* __restrict__ output,
    int channels,
    int spatial,
    int parts)
{
    int group = blockIdx.x / parts;
    int part = blockIdx.x % parts;
    int channel = group % channels;
    float mean = instance_mean[group];
    float scale = weight[channel] * stats[channel].y;
    int64_t base = (int64_t)group * spatial;

    #pragma unroll
    for (int j = 0; j < TILE / THREADS; ++j) {
        int position = part * TILE + threadIdx.x + j * THREADS;
        if (position < spatial) {
            output[base + position] = (x[base + position] - mean) * scale;
        }
    }
}

std::vector<torch::Tensor> centered_bn_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    torch::Tensor counter,
    bool training,
    double momentum,
    double eps)
{
    const 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);
    int parts = (spatial + TILE - 1) / TILE;
    int blocks = batch * channels * parts;

    auto output = torch::empty_like(x);
    auto partial = torch::empty({batch, channels, parts, 2}, x.options());
    auto means = torch::empty({batch, channels}, x.options());
    auto stats = torch::empty({channels, 2}, x.options());
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (training) {
        partial_moments<true><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            spatial, parts);

        finish_moments<true><<<channels, THREADS, 0, stream>>>(
            reinterpret_cast<const float2*>(partial.data_ptr<float>()),
            means.data_ptr<float>(),
            reinterpret_cast<float2*>(stats.data_ptr<float>()),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            counter.data_ptr<int64_t>(),
            batch, channels, spatial, parts, float(momentum), float(eps));
    } else {
        partial_moments<false><<<blocks, THREADS, 0, stream>>>(
            x.data_ptr<float>(),
            reinterpret_cast<float2*>(partial.data_ptr<float>()),
            spatial, parts);

        finish_moments<false><<<channels, THREADS, 0, stream>>>(
            reinterpret_cast<const float2*>(partial.data_ptr<float>()),
            means.data_ptr<float>(),
            reinterpret_cast<float2*>(stats.data_ptr<float>()),
            running_mean.data_ptr<float>(),
            running_var.data_ptr<float>(),
            counter.data_ptr<int64_t>(),
            batch, channels, spatial, parts, float(momentum), float(eps));
    }

    apply_centered_bn<<<blocks, THREADS, 0, stream>>>(
        x.data_ptr<float>(),
        weight.data_ptr<float>(),
        means.data_ptr<float>(),
        reinterpret_cast<const float2*>(stats.data_ptr<float>()),
        output.data_ptr<float>(),
        channels, spatial, parts);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return {output, means, stats};
}
"""

_centered_bn = load_inline(
    name="centered_bn3d_welford_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,
        counter, training, momentum, eps
    ):
        output, means, stats = _centered_bn.centered_bn_cuda(
            x, weight, running_mean, running_var,
            counter, training, momentum, eps
        )
        ctx.save_for_backward(x, weight, means, stats)
        ctx.training = training
        return output

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

        centered = x - means[:, :, None, None, None]
        invstd = stats[:, 1].reshape(channel_shape)

        grad_weight = (
            (grad_output * centered).sum(dim=batch_dims)
            * stats[:, 1]
        )
        grad_bias = torch.zeros_like(weight)

        grad_x = grad_output - grad_output.mean(
            dim=spatial_dims, keepdim=True
        )
        if ctx.training:
            batch_centered = x - stats[:, 0].reshape(channel_shape)
            covariance = (grad_output * centered).mean(
                dim=batch_dims, keepdim=True
            )
            grad_x = grad_x - batch_centered * invstd.square() * covariance

        grad_x = grad_x * weight.reshape(channel_shape) * invstd
        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

        use_custom = (
            x.is_cuda
            and x.dtype == torch.float32
            and x.numel() > 0
            and bn.affine
            and bn.track_running_stats
            and bn.momentum is not None
            and bn.weight.dtype == torch.float32
            and bn.running_mean.dtype == torch.float32
            and (not bn.training or x.numel() // x.shape[1] > 1)
        )

        if use_custom:
            return _CenteredBatchNorm.apply(
                x.contiguous(),
                bn.weight,
                bn.bias,
                bn.running_mean,
                bn.running_var,
                bn.num_batches_tracked,
                bn.training,
                bn.momentum,
                bn.eps,
            )

        x = bn(x)
        return x - x.mean(dim=(2, 3, 4), keepdim=True)
```