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


_cpp_source = r"""
torch::Tensor to_channels_last16_cuda(torch::Tensor x);

torch::Tensor centered_bn32_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor conv_bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    double eps,
    double momentum,
    bool training,
    bool inplace);
"""

_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 GROUPS = 8;
constexpr int TILE = 2048;
constexpr int OUTPUT_VECTORS = 4096;

__device__ __forceinline__ float4 zero4() {
    return make_float4(0.f, 0.f, 0.f, 0.f);
}

__device__ __forceinline__ float4 add4(float4 a, float4 b) {
    return make_float4(a.x+b.x, a.y+b.y, a.z+b.z, a.w+b.w);
}

__device__ __forceinline__ float4 mul4(float4 a, float b) {
    return make_float4(a.x*b, a.y*b, a.z*b, a.w*b);
}

__device__ __forceinline__ float4 square_add4(float4 a, float4 v) {
    return make_float4(
        fmaf(v.x, v.x, a.x), fmaf(v.y, v.y, a.y),
        fmaf(v.z, v.z, a.z), fmaf(v.w, v.w, a.w));
}

__device__ __forceinline__ float4 shuffle4(
    float4 a, int offset, unsigned mask = 0xffffffffu, int width = 32)
{
    return make_float4(
        __shfl_down_sync(mask, a.x, offset, width),
        __shfl_down_sync(mask, a.y, offset, width),
        __shfl_down_sync(mask, a.z, offset, width),
        __shfl_down_sync(mask, a.w, offset, width));
}

__device__ __forceinline__ float4 streaming_load(const float4* p) {
    float4 v;
    asm volatile(
        "ld.global.cs.v4.f32 {%0, %1, %2, %3}, [%4];"
        : "=f"(v.x), "=f"(v.y), "=f"(v.z), "=f"(v.w)
        : "l"(p));
    return v;
}

__device__ __forceinline__ void streaming_store(float4* p, float4 v) {
    asm volatile(
        "st.global.cs.v4.f32 [%0], {%1, %2, %3, %4};"
        :
        : "l"(p), "f"(v.x), "f"(v.y), "f"(v.z), "f"(v.w)
        : "memory");
}

__global__ void transpose16(
    const float* __restrict__ x,
    float4* __restrict__ out,
    int64_t spatial)
{
    __shared__ float tile[16][65];

    const int t = threadIdx.x;
    const int c = t >> 4;
    const int s = (t & 15) * 4;
    const int64_t begin = int64_t(blockIdx.x) * 64;
    const int64_t sample = int64_t(blockIdx.y) * spatial * 16;

    float4 v = zero4();
    if (begin + s < spatial) {
        v = *reinterpret_cast<const float4*>(
            x + sample + int64_t(c) * spatial + begin + s);
    }

    tile[c][s] = v.x;
    tile[c][s+1] = v.y;
    tile[c][s+2] = v.z;
    tile[c][s+3] = v.w;
    __syncthreads();

    const int os = t >> 2;
    const int oc = (t & 3) * 4;
    if (begin + os < spatial) {
        out[sample / 4 + (begin + os) * 4 + (t & 3)] =
            make_float4(
                tile[oc][os], tile[oc+1][os],
                tile[oc+2][os], tile[oc+3][os]);
    }
}

torch::Tensor to_channels_last16_cuda(torch::Tensor x) {
    TORCH_CHECK(
        x.is_cuda() && x.scalar_type() == at::kFloat &&
        x.dim() == 5 && x.size(1) == 16 && x.is_contiguous(),
        "Expected contiguous CUDA float32 NCDHW with 16 channels");

    const int64_t spatial = x.size(2) * x.size(3) * x.size(4);
    TORCH_CHECK(spatial % 4 == 0 && x.size(0) <= 65535,
                "Unsupported dimensions");

    const c10::cuda::CUDAGuard guard(x.device());
    auto out = torch::empty(
        x.sizes(), x.options().memory_format(at::MemoryFormat::ChannelsLast3d));

    transpose16<<<dim3((spatial + 63) / 64, x.size(0)), 256, 0,
                  at::cuda::getCurrentCUDAStream()>>>(
        x.data_ptr<float>(),
        reinterpret_cast<float4*>(out.data_ptr<float>()), spatial);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return out;
}

template<bool Training>
__global__ void partial_stats32(
    const float4* __restrict__ x,
    float4* __restrict__ sums,
    float4* __restrict__ squares,
    int spatial,
    int parts,
    int batch)
{
    const int group = threadIdx.x & 7;
    const int row = threadIdx.x >> 3;
    const int warp = threadIdx.x >> 5;
    const int n = blockIdx.y;
    const int p = blockIdx.x;
    const int begin = p * TILE;
    const int end = min(begin + TILE, spatial);
    const int64_t base = int64_t(n) * spatial * GROUPS;

    float4 sum0 = zero4();
    float4 sum1 = zero4();
    float4 square0 = zero4();
    float4 square1 = zero4();

    for (int start = begin + row; start < end; start += 128) {
        #pragma unroll
        for (int j = 0; j < 4; ++j) {
            const int s = start + j * 32;
            if (s < end) {
                const float4 v = streaming_load(
                    x + base + int64_t(s) * GROUPS + group);
                if (j & 1) {
                    sum1 = add4(sum1, v);
                    if (Training)
                        square1 = square_add4(square1, v);
                } else {
                    sum0 = add4(sum0, v);
                    if (Training)
                        square0 = square_add4(square0, v);
                }
            }
        }
    }

    float4 sum = add4(sum0, sum1);
    float4 square = add4(square0, square1);

    sum = add4(sum, shuffle4(sum, 16));
    sum = add4(sum, shuffle4(sum, 8));
    if (Training) {
        square = add4(square, shuffle4(square, 16));
        square = add4(square, shuffle4(square, 8));
    }

    __shared__ float4 shared_sum[64];
    __shared__ float4 shared_square[64];

    if ((threadIdx.x & 31) < 8) {
        shared_sum[warp * GROUPS + group] = sum;
        if (Training)
            shared_square[warp * GROUPS + group] = square;
    }
    __syncthreads();

    if (threadIdx.x < GROUPS) {
        sum = shared_sum[group];
        if (Training)
            square = shared_square[group];

        #pragma unroll
        for (int w = 1; w < 8; ++w) {
            sum = add4(sum, shared_sum[w * GROUPS + group]);
            if (Training)
                square = add4(square, shared_square[w * GROUPS + group]);
        }

        const int64_t dst = (int64_t(group) * batch + n) * parts + p;
        sums[dst] = sum;
        if (Training)
            squares[dst] = square;
    }
}

template<bool Training>
__device__ __forceinline__ float channel_scale(
    float mean,
    float second,
    float weight,
    float bias,
    float* running_mean,
    float* running_var,
    float eps,
    float momentum,
    float correction)
{
    float variance;
    if (Training) {
        variance = fmaxf(fmaf(-mean, mean, second), 0.f);
        *running_mean =
            (1.f - momentum) * (*running_mean) + momentum * (mean + bias);
        *running_var =
            (1.f - momentum) * (*running_var)
            + momentum * variance * correction;
    } else {
        variance = *running_var;
    }
    return weight * rsqrtf(variance + eps);
}

template<bool Training>
__global__ void finish_stats32(
    const float4* __restrict__ sums,
    const float4* __restrict__ squares,
    float4* __restrict__ means,
    float4* __restrict__ scales,
    const float* __restrict__ weight,
    const float* __restrict__ conv_bias,
    float* running_mean,
    float* running_var,
    int spatial,
    int parts,
    int batch,
    float eps,
    float momentum)
{
    const int group = blockIdx.x;
    const int lane = threadIdx.x & 15;
    const int row = threadIdx.x >> 4;
    const unsigned mask =
        (threadIdx.x & 16) ? 0xffff0000u : 0x0000ffffu;
    const float inverse_spatial = 1.f / float(spatial);

    float4 batch_sum = zero4();
    float4 batch_second = zero4();

    for (int n = row; n < batch; n += 16) {
        const int64_t base = (int64_t(group) * batch + n) * parts;
        float4 sum = zero4();
        float4 square = zero4();

        for (int p = lane; p < parts; p += 16) {
            sum = add4(sum, sums[base + p]);
            if (Training)
                square = add4(square, squares[base + p]);
        }

        #pragma unroll
        for (int offset = 8; offset > 0; offset >>= 1) {
            sum = add4(sum, shuffle4(sum, offset, mask, 16));
            if (Training)
                square = add4(square, shuffle4(square, offset, mask, 16));
        }

        if (lane == 0) {
            const float4 mean = mul4(sum, inverse_spatial);
            means[int64_t(n) * GROUPS + group] = mean;
            if (Training) {
                batch_sum = add4(batch_sum, mean);
                batch_second = add4(
                    batch_second, mul4(square, inverse_spatial));
            }
        }
    }

    __shared__ float4 shared_sum[16];
    __shared__ float4 shared_second[16];

    if (Training && lane == 0) {
        shared_sum[row] = batch_sum;
        shared_second[row] = batch_second;
    }
    __syncthreads();

    if (threadIdx.x < 16) {
        float4 mean = zero4();
        float4 second = zero4();

        if (Training) {
            mean = shared_sum[threadIdx.x];
            second = shared_second[threadIdx.x];

            #pragma unroll
            for (int offset = 8; offset > 0; offset >>= 1) {
                mean = add4(mean, shuffle4(mean, offset, 0xffffu, 16));
                second = add4(
                    second, shuffle4(second, offset, 0xffffu, 16));
            }
            mean = mul4(mean, 1.f / float(batch));
            second = mul4(second, 1.f / float(batch));
        }

        if (threadIdx.x == 0) {
            const int c = group * 4;
            const float count = float(batch) * float(spatial);
            const float correction =
                Training ? count / (count - 1.f) : 1.f;

            float4 scale;
            scale.x = channel_scale<Training>(
                mean.x, second.x, weight[c],
                conv_bias ? conv_bias[c] : 0.f,
                running_mean+c, running_var+c,
                eps, momentum, correction);
            scale.y = channel_scale<Training>(
                mean.y, second.y, weight[c+1],
                conv_bias ? conv_bias[c+1] : 0.f,
                running_mean+c+1, running_var+c+1,
                eps, momentum, correction);
            scale.z = channel_scale<Training>(
                mean.z, second.z, weight[c+2],
                conv_bias ? conv_bias[c+2] : 0.f,
                running_mean+c+2, running_var+c+2,
                eps, momentum, correction);
            scale.w = channel_scale<Training>(
                mean.w, second.w, weight[c+3],
                conv_bias ? conv_bias[c+3] : 0.f,
                running_mean+c+3, running_var+c+3,
                eps, momentum, correction);
            scales[group] = scale;
        }
    }
}

__device__ __forceinline__ float4 normalize4(
    float4 v, float4 mean, float4 scale)
{
    return make_float4(
        (v.x-mean.x)*scale.x,
        (v.y-mean.y)*scale.y,
        (v.z-mean.z)*scale.z,
        (v.w-mean.w)*scale.w);
}

__global__ void center32(
    const float4* x,
    const float4* __restrict__ means,
    const float4* __restrict__ scales,
    float4* output,
    int64_t sample_vectors)
{
    const int n = blockIdx.y;
    const int group = threadIdx.x & 7;
    const float4 mean = means[int64_t(n) * GROUPS + group];
    const float4 scale = scales[group];
    const int64_t base = int64_t(n) * sample_vectors;
    const int64_t begin =
        int64_t(blockIdx.x) * OUTPUT_VECTORS + threadIdx.x;

    if ((int64_t(blockIdx.x) + 1) * OUTPUT_VECTORS <= sample_vectors) {
        #pragma unroll
        for (int k = 0; k < OUTPUT_VECTORS; k += 1024) {
            const int64_t i = base + begin + k;
            const float4 a = streaming_load(x + i);
            const float4 b = streaming_load(x + i + 256);
            const float4 c = streaming_load(x + i + 512);
            const float4 d = streaming_load(x + i + 768);

            streaming_store(output + i, normalize4(a, mean, scale));
            streaming_store(output + i + 256, normalize4(b, mean, scale));
            streaming_store(output + i + 512, normalize4(c, mean, scale));
            streaming_store(output + i + 768, normalize4(d, mean, scale));
        }
    } else {
        for (int k = 0; k < OUTPUT_VECTORS; k += 256) {
            const int64_t i = begin + k;
            if (i < sample_vectors) {
                const float4 v = streaming_load(x + base + i);
                streaming_store(
                    output + base + i, normalize4(v, mean, scale));
            }
        }
    }
}

torch::Tensor centered_bn32_cuda(
    torch::Tensor x,
    torch::Tensor weight,
    torch::Tensor conv_bias,
    torch::Tensor running_mean,
    torch::Tensor running_var,
    double eps,
    double momentum,
    bool training,
    bool inplace)
{
    TORCH_CHECK(x.is_cuda() && x.scalar_type() == at::kFloat,
                "Expected CUDA float32 input");
    TORCH_CHECK(x.dim() == 5 && x.size(1) == 32 &&
                x.is_contiguous(at::MemoryFormat::ChannelsLast3d),
                "Expected channels-last 3D input with 32 channels");

    const c10::cuda::CUDAGuard guard(x.device());
    const int batch = x.size(0);
    const int64_t spatial64 = x.size(2) * x.size(3) * x.size(4);

    TORCH_CHECK(batch > 0 && batch <= 65535 &&
                spatial64 > 0 && spatial64 <= 2147481600LL,
                "Unsupported input dimensions");
    TORCH_CHECK(!training || int64_t(batch) * spatial64 > 1,
                "Expected more than one value per channel");

    const int spatial = static_cast<int>(spatial64);
    const int parts = (spatial + TILE - 1) / TILE;
    const int64_t scratch_vectors = int64_t(GROUPS) * batch * parts;

    auto partial = torch::empty(
        {(training ? 2 : 1) * scratch_vectors * 4}, x.options());
    auto statistics = torch::empty(
        {int64_t(batch + 1) * 32}, x.options());
    auto output = inplace ? x : torch::empty_like(x);

    auto* sums = reinterpret_cast<float4*>(partial.data_ptr<float>());
    auto* squares = training ? sums + scratch_vectors : nullptr;
    auto* means = reinterpret_cast<float4*>(statistics.data_ptr<float>());
    auto* scales = means + int64_t(batch) * GROUPS;
    const float* bias_ptr =
        conv_bias.numel() ? conv_bias.data_ptr<float>() : nullptr;
    const auto* input = reinterpret_cast<const float4*>(x.data_ptr<float>());
    cudaStream_t stream = at::cuda::getCurrentCUDAStream();

    if (training) {
        partial_stats32<true><<<dim3(parts, batch), 256, 0, stream>>>(
            input, sums, squares, spatial, parts, batch);
        finish_stats32<true><<<GROUPS, 256, 0, stream>>>(
            sums, squares, means, scales,
            weight.data_ptr<float>(), bias_ptr,
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            spatial, parts, batch, float(eps), float(momentum));
    } else {
        partial_stats32<false><<<dim3(parts, batch), 256, 0, stream>>>(
            input, sums, nullptr, spatial, parts, batch);
        finish_stats32<false><<<GROUPS, 256, 0, stream>>>(
            sums, nullptr, means, scales,
            weight.data_ptr<float>(), bias_ptr,
            running_mean.data_ptr<float>(), running_var.data_ptr<float>(),
            spatial, parts, batch, float(eps), float(momentum));
    }

    const int64_t sample_vectors = spatial64 * GROUPS;
    center32<<<dim3(
        (sample_vectors + OUTPUT_VECTORS - 1) / OUTPUT_VECTORS,
        batch), 256, 0, stream>>>(
            input, means, scales,
            reinterpret_cast<float4*>(output.data_ptr<float>()),
            sample_vectors);

    C10_CUDA_KERNEL_LAUNCH_CHECK();
    return output;
}
"""

_ops = load_inline(
    name="conv_transpose3d_centered_bn_streaming_v5",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["to_channels_last16_cuda", "centered_bn32_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class _ChannelsLast16(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x):
        return _ops.to_channels_last16_cuda(x)

    @staticmethod
    def backward(ctx, grad_output):
        return grad_output


class _CenteredBatchNorm(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x, weight, bias, conv_bias,
                running_mean, running_var, eps, momentum, training):
        output = _ops.centered_bn32_cuda(
            x, weight, conv_bias, running_mean, running_var,
            eps, momentum, training, False,
        )
        ctx.save_for_backward(
            x, weight,
            x.new_empty(0) if training else running_var.clone(),
        )
        ctx.eps = eps
        ctx.training = training
        ctx.has_conv_bias = conv_bias.numel() != 0
        return output

    @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.training:
            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.training:
            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)
        grad_conv_bias = (
            torch.zeros_like(weight) if ctx.has_conv_bias else None
        )

        return (
            grad_x, grad_weight, grad_bias, grad_conv_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)
        self.conv_transpose.to(memory_format=torch.channels_last_3d)

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

        if (
            not x.is_cuda
            or x.dim() != 5
            or x.dtype != torch.float32
            or conv.weight.dtype != torch.float32
            or conv.out_channels != 32
            or not bn.affine
            or not bn.track_running_stats
            or bn.weight.dtype != torch.float32
            or bn.running_mean.dtype != torch.float32
            or bn.running_var.dtype != torch.float32
            or not bn.weight.is_contiguous()
            or not bn.running_mean.is_contiguous()
            or not bn.running_var.is_contiguous()
            or (conv.bias is not None and (
                conv.bias.dtype != torch.float32
                or not conv.bias.is_contiguous()
            ))
            or x.numel() == 0
            or x.size(0) > 65535
        ):
            y = bn(conv(x))
            return y - y.mean(dim=(2, 3, 4), keepdim=True)

        if not x.is_contiguous(memory_format=torch.channels_last_3d):
            input_spatial = x.size(2) * x.size(3) * x.size(4)
            if (
                x.size(1) == 16
                and x.is_contiguous()
                and input_spatial % 4 == 0
                and x.storage_offset() % 4 == 0
            ):
                if torch.is_grad_enabled() and x.requires_grad:
                    x = _ChannelsLast16.apply(x)
                else:
                    x = _ops.to_channels_last16_cuda(x)
            else:
                x = x.contiguous(memory_format=torch.channels_last_3d)

        # Spatial centering removes both constant channel shifts.
        y = F.conv_transpose3d(
            x, conv.weight, None, conv.stride, conv.padding,
            conv.output_padding, conv.groups, conv.dilation,
        )

        spatial = y.size(2) * y.size(3) * y.size(4)
        if y.dtype != torch.float32 or spatial > 2147481600:
            if conv.bias is not None:
                y = y + conv.bias.to(y.dtype).view(1, -1, 1, 1, 1)
            y = bn(y)
            return y - y.mean(dim=(2, 3, 4), keepdim=True)

        y = y.contiguous(memory_format=torch.channels_last_3d)

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

        conv_bias = conv.bias if conv.bias is not None else y.new_empty(0)
        needs_grad = torch.is_grad_enabled() and (
            y.requires_grad
            or bn.weight.requires_grad
            or bn.bias.requires_grad
            or conv_bias.requires_grad
        )

        if not needs_grad:
            return _ops.centered_bn32_cuda(
                y, bn.weight, conv_bias,
                bn.running_mean, bn.running_var,
                bn.eps, momentum, bn.training, True,
            )

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